MatLAB gives NaN result on simple loop
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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;
2 comentarios
the cyclist
el 20 de Jun. de 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.
Respuestas (2)
Geoff Hayes
el 20 de Jun. de 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.
per isakson
el 20 de Jun. de 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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/225582/image.png)
2 comentarios
Steven Lord
el 20 de Jun. de 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.
Ver también
Categorías
Más información sobre Loops and Conditional Statements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!