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

2 visualizzazioni (ultimi 30 giorni)
Udaya Kiran S il 21 Mar 2023
Commentato: Udaya Kiran S il 21 Mar 2023
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 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Risposta accettata

Sam Chak il 21 Mar 2023
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 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Sam Chak il 21 Mar 2023
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 il 21 Mar 2023
Ya its clear now. Thanks

Accedi per commentare.

### Categorie

Scopri di più su Programming in Help Center e File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by