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

조회 수: 4(최근 30일)
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 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)))
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표시 이전 댓글 수: 2숨기기 이전 댓글 수: 2
Udaya Kiran S 2023년 3월 21일
Ya its clear now. Thanks

댓글을 달려면 로그인하십시오.

### 범주

Find more on Event Functions in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by