Why ode45 returns NaN for a system of differential equations?

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.

Answer 1

[Torsten]

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

[Hou X.Y]

thanks a lot for your answer,but the definition of alpha is

alpha_psn(j,i) = x(j).*log((c(i)./c0(i))*(c0(j)./c(j)))/(log(T./T0));

anyway，thanks for your help.

[Torsten]

Yes, the log(T/T0) term is the reason for the NaN's.

That's what I alluded to when I removed this term from the computation.

[Hou X.Y]

Thank you.So,how can I modify the code to get numerical results without removing log(T/T0)?

[Torsten]

By not setting T=T0 and thus dividing by log(1) = 0.

[Hou X.Y]

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?

clc,clear;
R = 8.314;
g = 9.8;
N = 10;
name = {'n2','o2','ar','h2o','co2','co','so2','o3','no2','c6h6'};
Mi(1) = 28.01;
c0(1) = 0.7809;
Mi(2) = 32;
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 = 294 + 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
%%
for i = 1:N
    alpha(i) = -log((c(i)/c0(i))*((sum_c0-c0(i))/(sum_c-c(i))))*((log(T./T0)));
    %individual height
    Hi(i) = 1/((Mi(i)*g)/(R*T));
    %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)
%%
% opts = odeset('RelTol',1e-1,'AbsTol',1e-1);
% for i = 10:10:100
%     intv(1) = c0;
%     while j <= 10
intv = c0;
D = [];
depth = [];
for i = 1:10
    tspan = [10*(i-1):0.1:10*i];
    [z,c] = ode23(F,tspan,intv)
    intv = c(101,:);
    D = [D;c];%求出depth 0：0.1：100
    depth = [depth,tspan]
end
figure
% sum_carray0 = 0;
sum_carray0 = 0;
for i = 1:size(D,1)
    sum_carray(i) = sum(D(i,:));%求D中每行的和，即每个深度上不同组分的浓度和
end
for i = 1:size(D,1)
    for j = 1:N
        concentn(i,j) = D(i,j)./sum_carray(i)
        hold on
    end
end
%%
grid on
for i = 1:N
    plot(depth,(concentn(:,i)),'color',[rand(),rand(),rand()],'linewidth',1.5)
    hold on
end
% for i = 1:N
%     semilogy(depth,(concentn(:,i)),'color',[rand(),rand(),rand()],'linewidth',1.5)
%     hold on
% end
legend(name)
xlabel('Depth(km)','fontsize',15)
ylabel('molar fraction','fontsize',15)

[Torsten]

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]

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.