MatLAB gives NaN result on simple loop
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
the cyclist
le 20 Juin 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.
Réponses (2)
Geoff Hayes
le 20 Juin 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
le 20 Juin 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 commentaires
Steven Lord
le 20 Juin 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.
Voir également
Catégories
En savoir plus sur Loops and Conditional Statements dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!