Why is the ODE45 producing NaN values for all the set of solutions?
조회 수: 2 (최근 30일)
이전 댓글 표시
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
채택된 답변
Sam Chak
2023년 3월 21일
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)))
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))
댓글 수: 3
Sam Chak
2023년 3월 21일
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
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개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!