Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. with Ode15s.

I'm writing a code about kinetic reaction in reactor solve 9 differiential equation, but now i'm stuck in Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

Here this code that I write

P0                                    =  [0.5]                                                ;%bar         (Intial pressure)
Ta0                                   =  [823]                                          ;%K           (Intial air heat temperature)              
repeat                                =  [1]                                                  ;%Repeat condition Ta constant and Ta non constant
sheet                                 =  0                                                      ;%Sheet in excel 
condition                             =  0                                                      ;%Number of condition
for a = 1:length(repeat)  
    for i = length(P0)
        for j = length(Ta0)
        T0                            = Ta0(j)                                                  ;%K           (Intial temperature before entering bed)
        yj0                           = [0.1 0.1 0 0 0 0 0.8]'	  
        pj0                           = yj0.*P0(i)                                              ;%bar         (Partial pressure)
        Fj0                           = (pj0.*v0_feed)./(R_feed.*T0)                            ;%mol/s       (Molar flowrate)
        Lspan                         = [0 0.03]                                                ;%m           (Bed length range)
        dep_var0                      = [Fj0' T0 Ta0(j)]                                        ;%[mol/s,K,K] (Intial dependent variable)
        condition                     = condition + 1                                           ;
        var                           = [P0(i) Ta0(j) condition T0]                             ;%[bar,K,-,K] (Another variable that nessery)          
        [L ,dep_var]                  = ode23s(@func_Ni,Lspan,dep_var0,[],var)                   ;
        FCH4                          = dep_var(:,1)                                            ;%mol/s
        FCO2                          = dep_var(:,2)                                            ;%mol/s
        FCO                           = dep_var(:,3)                                            ;%mol/s
        FH2                           = dep_var(:,4)                                            ;%mol/s
        FH2O                          = dep_var(:,5)                                            ;%mol/s
        FC                            = dep_var(:,6)                                            ;%mol/s
        FN2                           = dep_var(:,7)                                            ;%mol/s
        T                             = dep_var(:,8)                                            ;%K 
        Ta                            = dep_var(:,9)                                            ;%K     

For the functions

function dydx = func_Ni(y,x,op)
load Data_import_Ni
FCH4                    = x(1,1)                                                                                                                    ;
FCO2                    = x(2,1)                                                                                                                    ;
FCO                     = x(3,1)                                                                                                                    ;
FH2                     = x(4,1)                                                                                                                    ;
FH2O                    = x(5,1)                                                                                                                    ;
FC                      = x(6,1)                                                                                                                    ;
FN2                     = x(7,1)                                                                                                                    ;
T                       = x(8,1)                                                                                                                    ;
Ta                      = x(9,1)                                                                                                                    ;

P0                      = op(1)                                                                                                                     ;%bar         (Intial pressure)                   
Tac                     = op(2)                                                                                                                     ;%K           (Intial air heat temperature constant)
condition               = op(3)                                                                                                                     ; 
T0                      = op(4)                                                                                                                     ;%K           (Intial temperature)
FT                      = FCH4+FCO2+FCO+FH2+FH2O+FC+FN2                                                                                             ;%mol/s       (Total molar flowrate)                                                                                                            
%% Partial pressure (bar)
pCH4                    = P0*(FCH4/FT)*(T0/T)                                                                                                       ;
pCO2                    = P0*(FCO2/FT)*(T0/T)                                                                                                       ;
pCO                     = P0*(FCO/FT)*(T0/T)                                                                                                        ;
pH2                     = P0*(FH2/FT)*(T0/T)                                                                                                        ;
pH2O                    = P0*(FH2O/FT)*(T0/T)                                                                                                       ;
pC                      = P0*(FC/FT)*(T0/T)                                                                                                         ;
pN2                     = P0*(FN2/FT)*(T0/T)                                                                                                        ;
%% Heat capacity (J/(mol*K))                                                            
cpj                     = heatcp(:,1) + heatcp(:,2)*T + heatcp(:,3)*(T^2) + heatcp(:,4)*(T^-2)                                                      ;                                                      
%% Del_Heat of reaction
delH1                   = hrxi_Tr(1)+((2.*cpj(3)+2.*cpj(4)-cpj(1)-cpj(2)).*(T-Tr))                                                                  ;
delH2                   = hrxi_Tr(2)+((cpj(3)+cpj(5)-cpj(2)-cpj(4)).*(T-Tr))                                                                        ;
delH3                   = hrxi_Tr(3)+((cpj(6)+2.*cpj(4)-cpj(1)).*(T-Tr))                                                                            ;
delH4                   = hrxi_Tr(4)+((cpj(3)+cpj(4)-cpj(5)-cpj(6)).*(T-Tr))                                                                        ;
delH5                   = hrxi_Tr(5)+((2*cpj(3)-cpj(6)-cpj(2)).*(T-Tr))                                                                             ;
%% Reaction rate constant (mol/(kg*s))
ki                      = k0i.*exp(-Eai./T)                                                                                                         ;
%% Equilibrium constant
Kpi                     = Kp0i.*exp(-dHrxi_Tr./T)                                                                                                   ;%[bar^2,-,bar,bar,bar]
%% Adsorption equilibrium constant 
KCH4_1                  = K0j1(1)*exp(hdx1(1)/T)                                                                                                    ;%mol/(kg*s*bar)
KCO2_1                  = K0j1(2)*exp(hdx1(2)/T)                                                                                                    ;%mol/(kg*s*bar)

KCO2_2                  = K0j2(2)*exp(hdx2(2)/T)                                                                                                    ;%bar^-1
KH2_2                   = K0j2(4)*exp(hdx2(4)/T)                                                                                                    ;%bar^-1

KCH4_3                  = K0j3(1)*exp(hdx3(1)/T)                                                                                                    ;%bar^-1
KH2_3                   = K0j3(4)*exp(hdx3(4)/T)                                                                                                    ;%bar^1.5

KCH4_4                  = K0j4(1)*exp(hdx4(1)/T)                                                                                                    ;%bar^-1
KH2_4                   = K0j4(4)*exp(hdx4(4)/T)                                                                                                    ;%bar^1.5
KH2O_4                  = K0j4(5)*exp(hdx4(5)/T)                                                                                                    ;%-

KCO2_5                  = K0j5(2)*exp(hdx5(2)/T)                                                                                                    ;%bar
KCO_5                   = K0j5(3)*exp(hdx5(3)/T)                                                                                                    ;%bar^-1
%% Rate reaction (mol/(kg*s)
r1                      = ((ki(1)*KCH4_1*KCO2_1*pCH4*pCO2)/((1+KCO2_1*pCO2+KCH4_1*pCH4)^2))*(1-(((pCO*pH2)^2)/(Kpi(1)*pCH4*pCO2)))                  ; 
r2                      = ((ki(2)*KCO2_2*KH2_2*pCO2*pH2)/((1+KCO2_2*pCO2+KH2_2*pH2)^2))*(1-((pCO*pH2O)/(Kpi(2)*pCO2*pH2)))                         ;   
r3                      = ((ki(3)*KCH4_3)*(pCH4-((pH2^2)/Kpi(3))))/((1+KCH4_3*pCH4+(pH2^1.5/KH2_3))^2)                                              ;
r4                      = ((ki(4)/KH2O_4)*((pH2O/pH2)-(pCO/Kpi(4))))/((1+KCH4_4*pCH4+(pH2O/(KH2O_4*pH2))+(pH2^1.5/KH2_4))^2)                        ;
r5                      = ((ki(5)/(KCO_5*KCO2_5))*((pCO2/pCO)-(pCO/Kpi(5))))/((1+KCO_5*pCO+(pCO2/(KCO_5*KCO2_5*pCO)))^2)                            ;
%% rnet of species (mol/(kg*s))
rn_CH4                  =  -r1-r3                                                                                                                   ; 
rn_CO2                  =  -r1-r2-r5                                                                                                                ;
rn_CO                   =  2*r1+r2+r4+2*r5                                                                                                          ;
rn_H2                   =  2*r1-r2+2*r3+r4                                                                                                          ;
rn_H2O                  =  r2-r4                                                                                                                    ; 
rn_C                    =  r3-r4-r5                                                                                                                 ;
rn_N2                   =  0                                                                                                                        ;
%% 
sumFjcpj                =  FCH4*cpj(1)+ FCO2*cpj(2)+ FCO*cpj(3)+ FH2*cpj(4)+ FH2O*cpj(5)+ FC*cpj(6)+ FN2*cpj(7)                                     ;%Watt/K
sumrihi                 =  r1*delH1+ r2*delH2+ r3*delH3+ r4*delH4+ r5*delH5                                                                         ;%Watt/kg   
%% ODE equations 
dFCH4dL                 = Rho_bed*Ac*rn_CH4                                                                                                         ;%mol/(m*s)
dFCO2dL                 = Rho_bed*Ac*rn_CO2                                                                                                         ;%mol/(m*s)
dFCOdL                  = Rho_bed*Ac*rn_CO                                                                                                          ;%mol/(m*s)
dFH2dL                  = Rho_bed*Ac*rn_H2                                                                                                          ;%mol/(m*s)
dFH2OdL                 = Rho_bed*Ac*rn_H2O                                                                                                         ;%mol/(m*s)
dFCdL                   = Rho_bed*Ac*rn_C                                                                                                           ;%mol/(m*s)
dFN2dL                  = Rho_bed*Ac*rn_N2                                                                                                          ;%mol/(m*s)      
    if condition <= length(T0)*length(P0)
        dTdL            = ((Rho_bed*Ac*sumrihi)-(U*ain*Ac*(T-Tac)))/sumFjcpj                                                                        ;%K/m       
        dTadL           = 0                                                                                                                         ;%K/m                    
    else
        dTdL            = ((Rho_bed*Ac*sumrihi)-(U*ain*Ac*(T-Ta)))/sumFjcpj                                                                         ;%K/m  
        dTadL           = (U*aout*Ac_air*(T-Ta))/(mc_air*cpj(8))                                                                                    ;%K/m
    end
                                                               
%% Calculations
dydx(1,1)               = dFCH4dL                                                                                                                   ;
dydx(2,1)               = dFCO2dL                                                                                                                   ;
dydx(3,1)               = dFCOdL                                                                                                                    ;
dydx(4,1)               = dFH2dL                                                                                                                    ;
dydx(5,1)               = dFH2OdL                                                                                                                   ;
dydx(6,1)               = dFCdL                                                                                                                     ;
dydx(7,1)               = dFN2dL                                                                                                                    ;
dydx(8,1)               = dTdL                                                                                                                      ;
dydx(9,1)               = dTadL                                                                                                                     ;

end

[Walter Roberson]

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

Avoid load() inside your ode function.

[Jidapa Adam]

I don't understant. What do you mean?

[Jan]

load() inside the function to be integrated is extremely slow. Prefer to load the data once in the main code and provide them as inputs.

if condition <= length(T0)*length(P0)

This introduces a discontinuity in the function to be integrated. Remember, that Matlab's integrators are designed to handle smooth functions only. As a result the stepsize controller might drive crazy and maybe this does not even show a warning, but only the result is inaccurate.

Please post the complete error message, such that the readers do not have to guess, which line is failing.

[Torsten]

The if - statement does not depend on solution variables. So there is no discontinuity in the equations.

But testing your code is impossible since the calling part is incomplete and Data_import_Ni is missing. So it will be impossible for us to answer where the problem is.

[Jidapa Adam]

Okay. Thank you this is Data_import_Ni.

[Torsten]

The calling part of your program (the part above the function part) is still incomplete.

[Jidapa Adam]

I will add load Data_import_Ni part above the function part right?

load Data_import_Ni
P0                                    =  [0.5]                                                  %bar         (Intial pressure)
Ta0                                   =  [823]                                                  ;%K           (Intial air heat temperature)              
repeat                                =  [1]                                                    ;%Repeat condition Ta constant and Ta non constant
sheet                                 =  0                                                      ;%Sheet in excel 
condition                             =  0                                                      ;%Number of condition
for a = 1:length(repeat)  
    for i = length(P0)
        for j = length(Ta0)
        T0                            = Ta0(j)                                                  ;%K           (Intial temperature before entering bed)
        yj0                           = [0.1 0.1 0 0 0 0 0.8]'	  
        pj0                           = yj0.*P0(i)                                              ;%bar         (Partial pressure)
        Fj0                           = (pj0.*v0_feed)./(R_feed.*T0)                            ;%mol/s       (Molar flowrate)
        Lspan                         = [0 0.03]                                                ;%m           (Bed length range)
        dep_var0                      = [Fj0' T0 Ta0(j)]                                        ;%[mol/s,K,K] (Intial dependent variable)
        condition                     = condition + 1                                           ;
        var                           = [P0(i) Ta0(j) condition T0]                             ;%[bar,K,-,K] (Another variable that nessery)          
        [L ,dep_var]                  = ode23s(@func_Ni,Lspan,dep_var0,[],var)                   ;
        FCH4                          = dep_var(:,1)                                            ;%mol/s
        FCO2                          = dep_var(:,2)                                            ;%mol/s
        FCO                           = dep_var(:,3)                                            ;%mol/s
        FH2                           = dep_var(:,4)                                            ;%mol/s
        FH2O                          = dep_var(:,5)                                            ;%mol/s
        FC                            = dep_var(:,6)                                            ;%mol/s
        FN2                           = dep_var(:,7)                                            ;%mol/s
        T                             = dep_var(:,8)                                            ;%K 
        Ta                            = dep_var(:,9)                                            ;%K 

Answer 1

[Torsten]

P0                                    =  [0.5]                                                  %bar         (Intial pressure)
Ta0                                   =  [823]                                                  ;%K           (Intial air heat temperature)              
repeat                                =  [1]                                                    ;%Repeat condition Ta constant and Ta non constant
sheet                                 =  0                                                      ;%Sheet in excel 
condition                             =  0                                                      ;%Number of condition
for a = 1:length(repeat)  
    for i = length(P0)
        for j = length(Ta0)
            T0                            = Ta0(j)                                                  ;%K           (Intial temperature before entering bed)
            yj0                           = [0.1 0.1 1e-8 1e-8 1e-8 1e-8 0.8]'	  
            pj0                           = yj0.*P0(i)                                              ;%bar         (Partial pressure)
            Fj0                           = (pj0.*v0_feed)./(R_feed.*T0)                            ;%mol/s       (Molar flowrate)
            Lspan                         = [0 0.03]                                                ;%m           (Bed length range)
            dep_var0                      = [Fj0' T0 Ta0(j)]                                        ;%[mol/s,K,K] (Intial dependent variable)
            condition                     = condition + 1                                           ;
            var                           = [P0(i) Ta0(j) condition T0]                             ;%[bar,K,-,K] (Another variable that nessery)          
            [L ,dep_var]                  = ode15s(@(y,x)func_Ni(y,x,var),Lspan,dep_var0)                   ;
            FCH4                          = dep_var(:,1)                                            ;%mol/s
            FCO2                          = dep_var(:,2)                                            ;%mol/s
            FCO                           = dep_var(:,3)                                            ;%mol/s
            FH2                           = dep_var(:,4)                                            ;%mol/s
            FH2O                          = dep_var(:,5)                                            ;%mol/s
            FC                            = dep_var(:,6)                                            ;%mol/s
            FN2                           = dep_var(:,7)                                            ;%mol/s
            T                             = dep_var(:,8)                                            ;%K 
            Ta                            = dep_var(:,9)                                            ;%K 
        end
    end
end

By changing the initial conditions for the species and by modifying the call to the ODE solver, NaN values for the derivatives are no longer computed.

But you'll have to study the output from your function more carefully to make the solver really advance in time.

[Jidapa Adam]

Thank you very much. Where do you live? I will send the gift for you. I live in Thailand.

[Torsten]

I live in Germany. I'm happy if I could help - no need to send a gift.

Best

Torsten.