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

1 次查看（过去 30 天）
Udaya Kiran S
Udaya Kiran S 2023-3-21
评论： Udaya Kiran S 2023-3-21
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),'.');

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

采纳的回答

Sam Chak
Sam Chak 2023-3-21
Ran in:
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更早的评论
Sam Chak
Sam Chak 2023-3-21
Ran in:
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! 🙏

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by