MATLAB Answers

MatLAB gives NaN result on simple loop

Asked by Bertiningrum Devi

Bertiningrum Devi (view profile)

on 20 Jun 2019
Latest activity Commented on by Steven Lord

Steven Lord (view profile)

on 20 Jun 2019
Hi everyone!
I'm trying to solve non linear equation using a kind of simple loop, but then one of the variables (a2) gives NaN result which impact the other variables to be NaN.
as far as I know matLAb will give NaN if there's 0/0 or number/0. how could a2 be NaN since a2 is constant... and a3 = a2-(f2*(a2-a1))/(f2-f1) where the f2 is always not equal to f1.
Any help would be appreciate. Wish you good day, thank you!
Here is my code:
clc;
clear;
tol = 1e-5;
T = 94.56 + 273.15; %TL in K
Tref = 298.15;
P = 2; %colum pressure atm
R = 0.0821; %gas constant in L.atm/mol.K
k_CO2 = exp(231.465-(12092.1/T)-(36.7816*log(T)));
k_MDEA = exp(-46.09-(4756.9/T)+(6.4268*log(T)));
k_PZ = 4e10*exp(-4059.4/T);
k_H2S = exp(214.58-(12995.4/T)-(33.5471*log(T)));
K1 = k_CO2/k_MDEA;
K2 = k_H2S/k_MDEA;
K7 = 10^(285.52-(11988/T)+(0.061656*T)-(48.737*log(T)));
K4 = exp(220.067-(35.4189*log(T))-(12431.7/T))/K7;
K5 = exp(6.018-(4794.1/T)-(1.6744*log(T)));
K6 = exp(11.91-(4351/T));
wt_MDEA = 0.23039;
rho_CO2 =(P*44)/(R*T);
rho_H2S = (P*34)/(R*T);
rho_MDEA = (((-7.6/10^4)*(T-273.15))+1.057)*1000;
rho_H2O =(999.83952+(16.945176*(T-273.15))-((7.9870401*10^-3)*(T-273.15)^2)-((46.170461*10^-6)*(T-273.15)^3)+((105.56302*10^-9)*(T-273.15)^4)-((280.54253*10^-12)*(T-273.15)^5))/(1+(16.89785*10^-3*(T-273.15)));
V_H2O = (18/rho_H2O)*10^3;
V_CO2 = 22400;
V_H2S = V_CO2;
V_MDEA = (119.163/rho_MDEA)*10^3;
miu_H2O = 0.00002414*10^(247.8/(T-140))*1000;
miu_CO2 = 14.918*exp((0.61942*(log(T/Tref)))-(0.483713/(T/Tref))+(0.0652807/((T/Tref)^2))+0.418465);
miu_H2S = (5.448325+(0.0022148*exp((0.0028+(0.0072/T)+(72.734/(T^2)))*rho_H2S)));
miu_MDEA = -4.7986+(1476.9/(T-136.3343));
miu_MDEAX = exp(-26.12779397+5379.925131/T+0.029796484/T);
D_MDEAW = ((7*10^-8*((2.26*18)^0.5)*T)/(miu_H2O*V_MDEA^0.6))*10^-4;
D_MDEAL = (D_MDEAW*(miu_H2O^0.6)/(miu_MDEAX^0.6))*10^-4;
D_H2SMDEA = ((7*10^-8)*((1*119.163)^0.5)*T/(miu_MDEA*V_H2S^0.6))*10^-2;
D_CO2MDEA = ((7*10^-8*((1*119.163)^0.5)*T)/(miu_MDEA*V_CO2^0.6))*10^-4;
kL_CO2MDEA = 8.5*0.42*(((9.8*miu_CO2*0.001)/rho_CO2)^(1/3))*((D_CO2MDEA*rho_CO2)/(miu_CO2*0.001))^0.5;
kL_H2SMDEA = (8.5*0.42*(((9.8*miu_H2S*0.001)/rho_H2S)^(1/3))*((D_H2SMDEA*rho_H2S)/(miu_H2S*0.001))^0.5)/1000;
He_CO2 = exp(18.194+(-2808.5/T)+(2.5629*log(T))+(-0.0187*T));
He_H2S = exp(110.04+(-6789/T)+(-11.452*log(T))+(-0.0105*T));
%kmol/s
n_MDEAH_in = 0.07796;
n_PZH_in = 4.17e-9;
n_HCO3_in = 0.07429;
n_HS_in = 0.00366;
n_MDEA_in = 0.09576;
n_H2O_in = 1.26199;
n_PZ_in = 0.01602;
n_CO3_in = 1.77e-6;
n_H_in = 1.39e-7;
n_OH_in = 2.12e-10;
n_gasH2O_in = 0.0458;
v_Lin = 0.04414; %m3/s
v_Lout = 0.04458; %m3/s
x1 = 1;
y1 = 1;
z1 = 1;
fx1 = 0.0037;
fy1 = 0.0743;
fz1 = -7.7395e-05;
x2 = 0.3289;
y2 = 0.0790;
z2 = 0.4795;
fx2 = 0.0011;
fy2 = 0.0056;
fz2 = 0.0238;
x3 = (x1+x2)/2;
y3 = (y1+y2)/2;
z3 = (z1+z2)/2;
error=1;
while error > tol
x3s = x3;
y3s = y3;
z3s = z3;
x3 = x2-(x2-x1)*fx2/(fx2-fx1);
y3 = y2-(y2-y1)*fy2/(fy2-fy1);
z3 = z2-(z2-z1)*fz2/(fz2-fz1);
n_H2S_out = x3*n_HS_in;
n_CO2_out = y3*n_HCO3_in;
n_gasH2O_out = z3*n_gasH2O_in;
n_HS_reacted = n_H2S_out;
n_HCO3_reacted = n_CO2_out;
%~~~~~~~START EQUILIBRIUM CALCULATION USING SECANT~~~~~~~~~
a1=0.363;
a2=0.01;
M_PZH = a1 - (n_PZ_in/v_Lin + n_PZH_in/v_Lin);
M_MDEAH = (n_MDEAH_in - (n_HCO3_reacted - n_PZH_in) - n_HS_reacted)/v_Lout - M_PZH;
M_MDEA = (n_MDEA_in + (n_HCO3_reacted - n_PZH_in) + n_HS_reacted)/v_Lout + M_PZH;
M_OH = K5*M_MDEA/M_MDEAH;
M_H = M_PZH/(K6*a1);
f1 = K7-M_H*M_OH;
M_PZH = a2 - (n_PZ_in/v_Lin + n_PZH_in/v_Lin);
M_MDEAH = (n_MDEAH_in - (n_HCO3_reacted - n_PZH_in) - n_HS_reacted)/v_Lout - M_PZH;
M_MDEA = (n_MDEA_in + (n_HCO3_reacted - n_PZH_in) + n_HS_reacted)/v_Lout + M_PZH;
M_OH = K5*M_MDEA/M_MDEAH;
M_H = M_PZH/(K6*a2);
f2 = K7-M_H*M_OH;
e=1;
while e>=tol
a3 = a2-(f2*(a2-a1))/(f2-f1);
M_PZH = a3 - (n_PZ_in/v_Lin + n_PZH_in/v_Lin);
M_MDEAH = (n_MDEAH_in - (n_HCO3_reacted - n_PZH_in) - n_HS_reacted)/v_Lout - M_PZH;
M_MDEA = (n_MDEA_in + (n_HCO3_reacted - n_PZH_in) + n_HS_reacted)/v_Lout + M_PZH;
M_OH = K5*M_MDEA/M_MDEAH;
M_H = M_PZH/(K6*a3);
f3 = K7-M_H*M_OH;
e=abs(f3);
if abs(f1)>abs(f2)
a1=a2;
f1=f2;
a2=a3;
f2=f3;
else
a2=a3;
f2=f3;
a1=a1;
f1=f1;
end
end
%~~~~~~~END EQUILIBRIUM CALCULATION USING SECANT~~~~~~~~~
M_PZ = a3;
n_HS_out = n_HS_in - n_HS_reacted;
M_HS_out = n_HS_out/v_Lout;
n_HCO3_out = n_HCO3_in - n_HCO3_reacted;
M_HCO3_out = n_HCO3_out/v_Lout;
M_CO3 = K4*M_HCO3_out*M_OH;
M_H2S = M_MDEAH*M_HS_out/(M_MDEA*K2);
M_CO2 = M_MDEAH*M_HCO3_out/(M_MDEA*K1);
n_gas_out = n_H2S_out + n_CO2_out + n_gasH2O_in + 2.1161;
y_H2S_out = n_H2S_out/n_gas_out;
y_CO2_out = n_CO2_out/n_gas_out;
M_out = M_MDEA+M_MDEAH+M_PZ+M_PZH+M_H+M_OH+M_HS_out+M_HCO3_out+M_CO3+(n_gasH2O_in-n_gasH2O_out)/v_Lin+n_H2O_in/v_Lin;
x_H2S_out = y_H2S_out*P*101325/He_H2S;
x_CO2_out = y_CO2_out*P*101325/He_CO2;
M_H2S_intf = x_H2S_out*M_out;
M_CO2_intf = x_CO2_out*M_out;
k1 = k_CO2*M_MDEA+k_PZ*M_PZ;
Hatta_CO2 = D_CO2MDEA*k1/(kL_CO2MDEA^2);
E_CO2 = sqrt(Hatta_CO2);
N_CO2 = kL_CO2MDEA*E_CO2*(M_CO2-M_CO2_intf);
E_H2S = 1 + (M_MDEA*D_MDEAL/(D_H2SMDEA*M_H2S_intf));
N_H2S = kL_H2SMDEA*E_H2S*(M_H2S - M_H2S_intf);
Tg = 96+273.15;
N_H2O = abs(1000*(Tg - T)/(22569*1000*n_gasH2O_in*18));
intfarea = 25.71877213;
Vliquid = 0.72598;
fx3 = (n_HS_in - n_HS_out) - N_H2S;
fy3 = (n_HCO3_in - n_HCO3_out) - N_CO2;
fz3 = (n_gasH2O_in - n_gasH2O_out) - N_H2O;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
e1 = abs((x3-x3s)/x3s);
e2 = abs((y3-y3s)/y3s);
e3 = abs((z3-z3s)/z3s);
if e1>e2
error = e1;
else
error = e2;
end
if e3 > error
error = e3;
end
x1=x2;
fx1=fx2;
x2=x3;
fx2=fx3;
y1=y2;
fy1=fy2;
y2=y3;
fy2=fy3;
z1=z2;
fz1=fz2;
z2=z3;
fz2=fz3;
end
dbstop if naninf;

the cyclist

the cyclist (view profile)

on 20 Jun 2019
I find it a little confusing that in your question, you refer to
a3 = a2-(f2*(a2-a1))/(f2-f1)
but your line of code that triggers the NaN is
z3 = z2-(z2-z1)*fz2/(fz2-fz1);
However, setting my confusion aside for a moment ...
You claim that "f2 is always not equal to f1", but isn't true. When I run your code, that's exactly what happens. On the 5th iteration of your while loop,
fz2 == fz1
I haven't tried to diagnose it further than that yet.
Bertiningrum Devi

Bertiningrum Devi (view profile)

on 20 Jun 2019
hello, you're right it was the z3, I forgot to add some command there that's why it went NaN. it's solved. thank you for your help :)

Sign in to comment.

2 Answers

Geoff Hayes (view profile)

Answer by Geoff Hayes

Geoff Hayes (view profile)

on 20 Jun 2019

Bertiningrum - but a2 isn't constant
if abs(f1)>abs(f2)
a1=a2;
f1=f2;
a2=a3;
f2=f3;
else
a2=a3;
f2=f3;
a1=a1;
f1=f1;
end
and so in this case, a2 becomes NaN since a3 is NaN. If you trace this back far enough, I think that the problem is with
z3 = z2-(z2-z1)*fz2/(fz2-fz1);
where fz1 and fz2 are (near) identical and so their difference is zero. You probably want to add some sort of guard/check to ensure that there are no division by zero errors in your code.

1 Comment

Bertiningrum Devi

Bertiningrum Devi (view profile)

on 20 Jun 2019
hello Geoff, thank you for your help. you're right it was the z3, I forgot to add some command there that's why it went NaN. it's already solved :)

Sign in to comment.

per isakson (view profile)

Answer by per isakson

per isakson (view profile)

on 20 Jun 2019

"where the f2 is always not equal to f1"
My steps
• set dbstop if naninf
• started your script
• execution halted at line 77, z3 = z2-(z2-z1)*fz2/(fz2-fz1);
• steps past line 77
• inspected variable values of line 77
K>> z3
z3 =
NaN
K>> fz1
fz1 =
4.065758146820642e-19
K>> fz2
fz2 =
4.065758146820642e-19
K>>
Eventually the NaN spread to M_OH, f2 and a3

2 Comments

Bertiningrum Devi

Bertiningrum Devi (view profile)

on 20 Jun 2019
hello per isakson, thank you for your help. you're right it was the z3, I forgot to add some command there that's why it went NaN. it's already solved :)
can you teach me how to track what went NaN like you just did? just in case I experience another error someday.
Steven Lord

Steven Lord (view profile)

on 20 Jun 2019
Set a breakpoint. Specifically, set an "error breakpoint" to stop on NaN or Inf as described in the "Error Breakpoints" section on that documentation page.

Sign in to comment.