MatLAB gives NaN result on simple loop
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
the cyclist
am 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.
Antworten (2)
Geoff Hayes
am 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.
per isakson
am 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 Kommentare
Steven Lord
am 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.
Siehe auch
Kategorien
Mehr zu Loops and Conditional Statements finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!