Filter löschen
Filter löschen

Why ode45 returns NaN for a system of differential equations?

2 Ansichten (letzte 30 Tage)
Hou X.Y
Hou X.Y am 3 Apr. 2022
Kommentiert: Hou X.Y am 5 Apr. 2022
clc,clear;
R = 8.314;
g = 9.8;
N = 10;
name = {'n2';'o2';'ar';'h2o';'h2o';'co2';'co';'so2';'o3';'no2';'c6h6'};
Mi(1) = 14.01;
c0(1) = 0.7809;
Mi(2) = 16;
c0(2) = 0.2095;
Mi(3) = 39.95;
c0(3) = 0.0093
Mi(4) = 18.02
c0(4) = 0.01
Mi(5) = 44.01
c0(5) = 500*0.000001
Mi(6) = 28.01
c0(6) = 20*0.000001
Mi(7) = 64.07
c0(7) = 75*0.000001
Mi(8) = 48
c0(8) = 70*1E-9
Mi(9) = 46.01
c0(9) = 53*1E-9
Mi(10) = 78.11
c0(10) = 50*1E-9
dTdz = 0.027
T0 = 293
syms T z x
%%
%initial moler concentration
% sum_mirho0_1 = 0;
% for i = 1:N
% sum_mirho0 = sum_mirho0_1+mirho0(i);
% sum_mirho0_1 = sum_mirho0
% end
% for i = 1:N
% rho0av = 1/sum_mirho0
% c0(i) = rho0av*m(i)*1000/Mi(i);
% end
%%
% dTdz = input('temperature gradient(K/m):')
% T0 = input('entrance temperature(K):')
T = T0 + dTdz*1000*z;%温度随深度的变化
%thermal diffusion factor
%%
syms z
c = sym('c',[N,1])
sum_c = sum(c);
sum_c01 = 0;
for i = 1:N
sum_c0 = sum_c01+c0(i)
sum_c01 = sum_c0
end
%%
x = c./sum_c;
alpha_0 = 0;
i = 1;
for i = 1:10
j = 1;
while j <= 10
if j ~= i
alpha_psn(j,i) = x(j).*log((c(i)./c0(i))*(c0(j)./c(j))).*(log(T./T0)^-1)
% alpha(i) = alpha_0+x(j)*alpha_psn(j,i)
% alpha_0 = alpha(i)
end
j = j+1;
end
end
%%
%individual height
for i = 1:10
Hi(i) = 1./((Mi(i)*g)./(R*T));
alpha(i) = sum(alpha_psn(i,:))
%ode for every species
dc(i) = -c(i).*((dTdz*1000.*((1 + alpha(i))./T)-(1./Hi(i))))
end
%%
f = matlabFunction(dc');
F = @(z,c)f(c(1),c(2),c(3),c(4),c(5),c(6),c(7),c(8),c(9),c(10),z)
% F = @(z,c)f(c(1),c(2),c(3),z)
%%
[z,c] = ode23(F,[0:100],c0)
figure
sum_carray0 = 0;
for i = 1:N
sum_carray = sum_carray0+c(:,i);
sum_carray0 = sum_carray;
end
for i = 1:N
plot(z,(c(:,i)./sum_carray),'linewidth',1.5)
hold on
end
grid on
legend(name)
xlabel('Depth(km)','fontsize',15)
ylabel('molar fraction','fontsize',15)
It returns a lot of NaN,and when I switch tspan [0:100] to other value,it may have some numerical results,and the result changes with the value of tspan,especially the initial points t0,for example,it may return some numerical results when I set 1 as t0.
This is the problem in the literature, which says the numerical solution can be obtained by ode23.

Akzeptierte Antwort

Torsten
Torsten am 3 Apr. 2022
N = 10;
Mi = zeros(N,1);
c0 = zeros(N,1);
Mi(1) = 14.01;
c0(1) = 0.7809;
Mi(2) = 16;
c0(2) = 0.2095;
Mi(3) = 39.95;
c0(3) = 0.0093;
Mi(4) = 18.02;
c0(4) = 0.01;
Mi(5) = 44.01;
c0(5) = 500*0.000001;
Mi(6) = 28.01;
c0(6) = 20*0.000001;
Mi(7) = 64.07;
c0(7) = 75*0.000001;
Mi(8) = 48;
c0(8) = 70*1E-9;
Mi(9) = 46.01;
c0(9) = 53*1E-9;
Mi(10) = 78.11;
c0(10) = 50*1E-9;
%f = matlabFunction(dc');
%F = @(z,c)f(c(1),c(2),c(3),c(4),c(5),c(6),c(7),c(8),c(9),c(10),z)
% F = @(z,c)f(c(1),c(2),c(3),z)
[z,c] = ode23(@(z,c)F(z,c,Mi,c0),[0:100],c0)
figure
sum_carray0 = 0;
for i = 1:N
sum_carray = sum_carray0+c(:,i);
sum_carray0 = sum_carray;
end
for i = 1:N
plot(z,(c(:,i)./sum_carray),'linewidth',1.5)
hold on
end
grid on
legend(name)
xlabel('Depth(km)','fontsize',15)
ylabel('molar fraction','fontsize',15)
function dc = F(z,c,Mi,c0);
R = 8.314;
g = 9.8;
N = 10;
name = {'n2';'o2';'ar';'h2o';'h2o';'co2';'co';'so2';'o3';'no2';'c6h6'};
dTdz = 0.027;
T0 = 293;
%syms T z x
%%
%initial moler concentration
% sum_mirho0_1 = 0;
% for i = 1:N
% sum_mirho0 = sum_mirho0_1+mirho0(i);
% sum_mirho0_1 = sum_mirho0
% end
% for i = 1:N
% rho0av = 1/sum_mirho0
% c0(i) = rho0av*m(i)*1000/Mi(i);
% end
%%
% dTdz = input('temperature gradient(K/m):')
% T0 = input('entrance temperature(K):')
T = T0 + dTdz*1000*z;%温度随深度的变化
%thermal diffusion factor
%%
%syms z
%c = sym('c',[N,1])
sum_c = sum(c);
sum_c01 = 0;
for i = 1:N
sum_c0 = sum_c01+c0(i);
sum_c01 = sum_c0;
end
%%
x = c./sum_c;
alpha_0 = 0;
i = 1;
for i = 1:10
j = 1;
while j <= 10
if j ~= i
alpha_psn(j,i) = x(j).*log((c(i)./c0(i))*(c0(j)./c(j)));%.*(log(T./T0)^-1);
% alpha(i) = alpha_0+x(j)*alpha_psn(j,i)
% alpha_0 = alpha(i)
end
j = j+1;
end
end
%%
%individual height
for i = 1:10
Hi(i) = 1./((Mi(i)*g)./(R*T));
alpha(i) = sum(alpha_psn(i,:));
%ode for every species
dc(i) = -c(i).*((dTdz*1000.*((1 + alpha(i))./T)-(1./Hi(i))));
end
%%
dc = dc.';
end
  7 Kommentare
Torsten
Torsten am 4 Apr. 2022
I find that if I don't set tspan as [0,100] directly,and divide it into [0,10],[10,20],[20,30],...,[90,100],it doesn't return NaN,and I don't konw why.If you know why,could you please tell me?
The reason is not the splitting of tspan.
Your code is different now.
Last time, you multiplied by log(T./T0)^-1 which is infinity since T = T0 at z = 0.
Now you multiply by log(T./T0) which gives 0 for T = T0 at z=0.
Hou X.Y
Hou X.Y am 5 Apr. 2022
Sorry,I am too careless.
But if I switch ' log(T./T0)' to 'log(T./T0)^-1',it still has numerical results,but this result is not consistent with the literature.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Tags

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by