Why is the ODE45 producing NaN values for all the set of solutions?

FUNCTION FILE:

function f = as1(z,F)

if z>=0 && z < 9.5

qz = 96;

elseif z >= 9.5 && z < 19.0

qz = 84;

elseif z >= 19.0 && z < 28.5

qz = 80;

elseif z >= 28.5 && z < 38.0

qz = 71;

elseif z>= 38.0 && z < 47.5

qz = 63;

else

qz = 59;

end

dt = 0.108;

rb = 0.178;

mu = (0.4/1.4)*(0.035*(10^(-3))) + (1/1.4)*(2.43*(10^(-5))) ;

G = 1.4*68.68;

omega = (3.14159*(dt^2)/4);

alpha = -1; %taking it as sum of all stoichiometric coefficients

R = 8.314;

Re = dt*G/mu ;

Fric = 0.046*(Re^(-0.2));

zeta = (0.7 + (180/90)*0.35)*(0.051 + 0.19*(dt/rb));

a = [0 0 1 -1 0 0 0 1 ; 0 0 -1 1 0 0 0 -1 ; 1 0 0 -2 0 1 0 0 ; 1 1 0 0 -1 0 0 0 ; -1 -1 0 0 1 0 0 0 ; 0 -1 -1 0 0 0 1 0 ; 1 0 -1 -1 1 0 0 0 ];

cp = [35.8,68.25,42.9,3.846,1.48254,1.68,3.1031,14.307];

Hf0 = [-74.87, 227.4, 52.3, -84.68, 20.4, -103.3, 145.3, 0 ];

H = [Hf0(1,1) + cp(1,1)*(F(9) - 298),Hf0(1,2) + cp(1,2)*(F(9) - 298),Hf0(1,3) + cp(1,3)*(F(9) - 298),Hf0(1,4) + cp(1,4)*(F(9) - 298),...

Hf0(1,5) + cp(1,5)*(F(9) - 298),Hf0(1,6) + cp(1,6)*(F(9) - 298),Hf0(1,7) + cp(1,7)*(F(9) - 298),Hf0(1,8) + cp(1,8)*(F(9) - 298)];

k1 = (4.65*10^13)*exp(-273020/(R*F(9)));

k2 = (8.75*10^8)*exp(-136870/(R*F(9)));

k3 = (3.85*10^11)*exp(-273190/(R*F(9)));

k4 = (9.81*10^8)*exp(-154680/(R*F(9)));

k5 = (5.87*10^4)*exp(-29480/(R*F(9)));

k6 = (1.03*10^12)*exp(-172750/(R*F(9)));

k7 = (7.08*10^13)*exp(-253010/(R*F(9)));

Ft = (F(1) + F(2) + F(3) + F(4) + F(5) + F(6) + F(7) + F(8));

c1 = (F(1)/Ft)*(F(10)/(R*F(9))) ;

c2 = (F(2)/Ft)*(F(10)/(R*F(9))) ;

c3 = (F(3)/Ft)*(F(10)/(R*F(9))) ;

c4 = (F(4)/Ft)*(F(10)/(R*F(9))) ;

c5 = (F(5)/Ft)*(F(10)/(R*F(9))) ;

c6 = (F(6)/Ft)*(F(10)/(R*F(9))) ;

c7 = (F(7)/Ft)*(F(10)/(R*F(9))) ;

c8 = (F(8)/Ft)*(F(10)/(R*F(9))) ;

r1 = k1*(c3^a(1,3))*(c4^a(1,4))*(c8^a(1,8));

r2 = k2*(c3^a(2,3))*(c4^a(2,4))*(c8^a(2,8));

r3 = k3*(c1^a(3,1))*(c4*a(3,4))*(c6*a(3,6));

r4 = k4*(c1^a(4,1))*(c2^a(4,2))*(c5^a(4,5));

r5 = k5*(c1^a(5,1))*(c2^a(5,2))*(c5^a(5,5));

r6 = k6*(c2^a(6,2))*(c3^a(6,3))*(c7^a(6,7));

r7 = k7*(c1^a(7,1))*(c3^a(7,3))*(c4^a(7,4))*(c5^a(7,5));

f1 = (r3*a(3,1)+r4*a(4,1)+r5*a(5,1)+r7*a(7,1))*(3.14159*(dt^2)/4) ;

f2 = (r4*a(4,2)+r5*a(5,2)+r6*a(6,2))*(3.14159*(dt^2)/4) ;

f3 = (r1*a(1,3)+r2*a(2,3)+r6*a(6,3)+r7*a(7,3))*(3.14159*(dt^2)/4) ;

f4 = (r1*a(1,4)+r2*a(2,4)+r3*a(3,4)+r7*a(7,4))*(3.14159*(dt^2)/4) ;

f5 = (r4*a(4,5)+r5*a(5,5)+r7*a(7,5))*(3.14159*(dt^2)/4) ;

f6 = (r3*a(3,6))*(3.14159*(dt^2)/4) ;

f7 = (r6*a(6,7))*(3.14159*(dt^2)/4) ;

f8 = (r1*a(1,8)+r2*a(2,8))*(3.14159*(dt^2)/4) ;

FCpt = F(1)*cp(1,1) + F(2)*cp(1,2) + F(3)*cp(1,3) + F(4)*cp(1,4) + F(5)*cp(1,5) + F(6)*cp(1,6) + F(7)*cp(1,7) + F(8)*cp(1,8) ;

delH = [-sum(H.*a(1,:)),-sum(H.*a(2,:)),-sum(H.*a(3,:)),-sum(H.*a(4,:)),-sum(H.*a(5,:)),-sum(H.*a(6,:)),-sum(H.*a(7,:))];

delHR = delH(1,1)*r1 + delH(1,2)*r2 + delH(1,3)*r3 + delH(1,4)*r4 + delH(1,5)*r5 + delH(1,6)*r6 + delH(1,7)*r7;

f9 = (1/FCpt)*( (qz*3.14159*dt) + ((3.14159*(dt^2)/4)*(delHR)) ) ;

f11 = (f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8)/(G*omega) ;

f10 = (f11 + (F(11))*((1/F(9))*f9 + ((2*Fric/dt)+ (zeta/(3.14159*rb)))))/((F(11)/F(10))-(F(10)/(alpha*(G^2)*R*F(9))));

f = [f1;f2;f3;f4;f5;f6;f7;f8;f9;f10;f11];

end

MAIN FILE:

clc;

clear;

close all;

[z,X] = ode45(@as1,[0 95],[0;0;0.209;20.5238;0.1672;0;0;0;953.13;3.03*(10^5);33.249]);

plot(z,(X(1,4)-X(:,4))/X(1,4),'.');

Answer 1

[Sam Chak]

Hi @Udaya Kiran S

This is how I checked. It is because some parameters in produce Inf and NaN, where your ODEs depend on.

You need to check whether the values and formulas are entered correctly or not.

F = [0; 0; 0.209; 20.5238; 0.1672; 0; 0; 0; 953.13; 3.03*(10^5); 33.249];
   
dt    = 0.108;
rb    = 0.178;
mu    = (0.4/1.4)*(0.035*(10^(-3))) + (1/1.4)*(2.43*(10^(-5)));
G     = 1.4*68.68;
omega = (3.14159*(dt^2)/4);
alpha = -1; %taking it as sum of all stoichiometric coefficients
R     = 8.314;
Re    = dt*G/mu;
Fric  = 0.046*(Re^(-0.2));
zeta  = (0.7 + (180/90)*0.35)*(0.051 + 0.19*(dt/rb));
a     = [0 0 1 -1 0 0 0 1 ; 0 0 -1 1 0 0 0 -1 ; 1 0 0 -2 0 1 0 0 ; 1 1 0 0 -1 0 0 0 ; -1 -1 0 0 1 0 0 0 ; 0 -1 -1 0 0 0 1 0 ; 1 0 -1 -1 1 0 0 0 ];
cp    = [35.8,68.25,42.9,3.846,1.48254,1.68,3.1031,14.307];
Hf0   = [-74.87, 227.4, 52.3, -84.68, 20.4, -103.3, 145.3, 0 ];
H     = [Hf0(1,1) + cp(1,1)*(F(9) - 298),Hf0(1,2) + cp(1,2)*(F(9) - 298),Hf0(1,3) + cp(1,3)*(F(9) - 298),Hf0(1,4) + cp(1,4)*(F(9) - 298),...
    Hf0(1,5) + cp(1,5)*(F(9) - 298),Hf0(1,6) + cp(1,6)*(F(9) - 298),Hf0(1,7) + cp(1,7)*(F(9) - 298),Hf0(1,8) + cp(1,8)*(F(9) - 298)];

k1 = (4.65*10^13)*exp(-273020/(R*F(9)));
k2 = (8.75*10^8)*exp(-136870/(R*F(9)));
k3 = (3.85*10^11)*exp(-273190/(R*F(9)));
k4 = (9.81*10^8)*exp(-154680/(R*F(9)));
k5 = (5.87*10^4)*exp(-29480/(R*F(9)));
k6 = (1.03*10^12)*exp(-172750/(R*F(9)));
k7 = (7.08*10^13)*exp(-253010/(R*F(9)));
Ft = (F(1) + F(2) + F(3) + F(4) + F(5) + F(6) + F(7) + F(8));
c1 = (F(1)/Ft)*(F(10)/(R*F(9)))

c1 = 0

c2 = (F(2)/Ft)*(F(10)/(R*F(9)))

c2 = 0

c3 = (F(3)/Ft)*(F(10)/(R*F(9))) 

c3 = 0.3824

c4 = (F(4)/Ft)*(F(10)/(R*F(9))) 

c4 = 37.5484

c5 = (F(5)/Ft)*(F(10)/(R*F(9))) 

c5 = 0.3059

c6 = (F(6)/Ft)*(F(10)/(R*F(9))) 

c6 = 0

c7 = (F(7)/Ft)*(F(10)/(R*F(9))) 

c7 = 0

c8 = (F(8)/Ft)*(F(10)/(R*F(9))) 

c8 = 0

r1 = k1*(c3^a(1,3))*(c4^a(1,4))*(c8^a(1,8))

r1 = 0

r2 = k2*(c3^a(2,3))*(c4^a(2,4))*(c8^a(2,8))

r2 = Inf

r3 = k3*(c1^a(3,1))*(c4*a(3,4))*(c6*a(3,6))

r3 = 0

r4 = k4*(c1^a(4,1))*(c2^a(4,2))*(c5^a(4,5))

r4 = 0

r5 = k5*(c1^a(5,1))*(c2^a(5,2))*(c5^a(5,5))

r5 = Inf

r6 = k6*(c2^a(6,2))*(c3^a(6,3))*(c7^a(6,7))

r6 = NaN

r7 = k7*(c1^a(7,1))*(c3^a(7,3))*(c4^a(7,4))*(c5^a(7,5))

r7 = 0

[Udaya Kiran S]

oh thanks man

btw how did you check it like that (just wanna know)

[Sam Chak]

You are welcome, @Udaya Kiran S.

To check, you need know that Inf is generally caused by division-by-zero. So, I searched for that. Since the user-supplied parameter values are given and the initial values are known, then making substitutions to check if the ODEs at contain Inf.

Look up, has a term . Since and , it is interpreted as that returns Inf. Another one, has a product term that returns NaN.

Inf*0

ans = NaN

In conclusion, the problem come from elements in matrix a and the parameters (that depends on the initial values). Change the values so that they don't result in Inf and NaN. If you are happy with the problem diagnosis, please consider accepting ✔ and voting 👍 the Answer. Thanks a bunch! 🙏

[Udaya Kiran S]

Ya its clear now. Thanks