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

Question

Udaya Kiran S 2023-3-21

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1932375-why-is-the-ode45-producing-nan-values-for-all-the-set-of-solutions

评论： Udaya Kiran S 2023-3-21

采纳的回答： Sam Chak

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),'.');

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Sam Chak 2023-3-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1932375-why-is-the-ode45-producing-nan-values-for-all-the-set-of-solutions#answer_1197435

在 MATLAB Online 中打开

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Sam Chak 2023-3-21

在 MATLAB Online 中打开

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 2023-3-21