## MatLAB gives NaN result on simple loop

### 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 :)

### 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 :)

### per isakson (view profile)

on 20 Jun 2019

"where the f2 is always not equal to f1"
My steps
• set dbstop if naninf
• 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

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.