Asked by Bertiningrum Devi
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;

Answer by Geoff Hayes
on 20 Jun 2019

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.

Bertiningrum Devi
on 20 Jun 2019

Sign in to comment.

Answer by per isakson
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

Bertiningrum Devi
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
on 20 Jun 2019

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## the cyclist (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/467981-matlab-gives-nan-result-on-simple-loop#comment_716469

## Bertiningrum Devi (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/467981-matlab-gives-nan-result-on-simple-loop#comment_716538

Sign in to comment.