MatLAB gives NaN result on simple loop

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]

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]

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

Answer 1

[Geoff Hayes]

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.

[Bertiningrum Devi]

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

Answer 2

[per isakson]

"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]

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]

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.