Getting NaN with ODE45/15s

i have a system of five coupled odes which i need to solve and i'm trying to use ode45/15s to solve it.

i put the odes in a function and i'm using a script to solve the given function

i'm getting NaN as the answer

double checked the equations, nothing's wrong with them (same equations in several papers)

double checked if i have a zero in denominator of any term, didn't have any so im all out of options and here i am asking you guys to kindly review my code and see if you can find what's wrong with it.

thank you all in advance

. . . FUNCTION . . .

function xdot=thermal_mass(t,x)

global pa f R0 T0 P0 Nwa

s=0.0725; %Surface Tension

dl=998; %Water Density

c=1481; %Speed of Sound

u=1e-3; %Viscosity

pv=2.33e3; %Vapor Pressure

ps=1.01e5; %Static Pressure

K=1.38e-23; %Boltzman Constant

b=5.1e-29; %Water/Argon Covolume

e1=2295;

e2=5255;

e3=5400;

P0=ps-pv+2*s/R0;

Ntot=4*pi*P0*(R0^3)/(3*K*T0+3*P0*b);

Nwa=0.95*Ntot;

Nar=0.05*Ntot;

Cv=1.5*Nar*K+(((e1/x(3))*exp(e1/x(3)))/(((exp(e1/x(3))-1)^2))+((e2/x(3))*exp(e2/x(3)))/(((exp(e2/x(3))-1)^2))+((e3/x(3))*exp(e3/x(3)))/(((exp(e3/x(3))-1)^2)))*K*x(4);

v=(4/3)*pi*x(1)^3;

M=(x(4)*18e-3+Nar*40e-3)/(6.022e23);

den=M/v;

nAr=Nar/v;

nH2O=5.9e23;

n0=2.446e25;

D0=23.55e-6;

D=D0*(n0/(nH2O+nAr));

ni0=5.9e23;

ni=x(4)/v;

ld=min(sqrt((x(1)*D)/(abs(x(2)))),x(1)/pi);

Nd=4*pi*(x(1)^2)*D*(ni0-ni)/ld;

L=17.9e-3;

Cp=2.5*Nar*K+4*x(4)*K;

a=L/(den*Cp);

Lth=min(sqrt(x(1)*a/abs(x(2))),x(1)/pi);

Td=((4*pi*(x(1)^2))/Cv)*((L/Lth)*(T0-x(3))-x(5)*x(2))+(4*T0-3*x(3)-x(3)*(((e1/x(3))/(exp(e1/x(3))-1))+((e2/x(3))/(exp(e2/x(3))-1))+((e3/x(3))/(exp(e3/x(3))-1))))*K*Nd/Cv;

Pgd=((Nd*K*x(3)+x(4)*K*Td)/((4/3)*pi*x(1)^3-(x(4)+Nar)*b))+(((x(4)+Nar)*K*x(3))*(Nd*b-4*pi*x(2)*(x(1)^2))/(((4/3)*pi*x(1)^3-(x(4)+Nar)*b)^2));

xdot(1,1)=x(2);

xdot(2,1)=(-0.5*(x(2)^2)*(3-x(2)/c)+(1+x(2)/c)*(x(5)/dl)+(x(1)/(dl*c))*Pgd-(2*s/(dl*x(1)))-(4*u*x(2)/(dl*x(1)))-(1+x(2)/c)*((ps-pv+pa*sin(2*pi*f*t))/dl)-(2*x(1)*pi*f*cos(2*pi*f*t)/(dl*c)))/(x(1)*(1-x(2)/c)+4*u/(dl*c));

xdot(3,1)=((4*pi*(x(1)^2))/Cv)*((L/Lth)*(T0-x(3))-x(5)*x(2))+(4*T0-3*x(3)-x(3)*(((e1/x(3))/(exp(e1/x(3))-1))+((e2/x(3))/(exp(e2/x(3))-1))+((e3/x(3))/(exp(e3/x(3))-1))))*K*Nd/Cv;

xdot(4,1)=4*pi*(x(1)^2)*D*(ni0-ni)/ld;

xdot(5,1)=((Nd*K*x(3)+x(4)*K*Td)/((4/3)*pi*x(1)^3-(x(4)+Nar)*b))+(((x(4)+Nar)*K*x(3))*(Nd*b-4*pi*x(2)*(x(1)^2))/(((4/3)*pi*x(1)^3-(x(4)+Nar)*b)^2));

. . . SCRIPT . . .

clear;

global pa f R0 T0 P0 Nwa

pa=1.2e5;

f=26.5e3;

R0=4e-6;

T0=300;

options= odeset('RelTol',1e-12,'AbsTol',1e-13);

[t,r]=ode45(@thermal_mass,(0:0.001/f:10/f),[R0 0 T0 Nwa P0],options);

Answer 1

[Walter Roberson]

Your script does not initialize Nwa before it uses it, so Nwa will have its default value of [] (because global variables default to [] ). This results in only 4 values being passed as your initial conditions to thermal_mass, which expects 5 values for x. Your thermal_mass function assigns to Nwa but by then it is too late.

You should avoid using global variables. If you have fixed variables then parameterize your function. If you need the updated Nwa to be available afterwards (there is no evidence of that at the moment) then you should construct nested functions with shared variables.

Your x(2) starts out as 0, so the line

ld=min(sqrt((x(1)*D)/(abs(x(2)))),x(1)/pi);

contains a division by 0, resulting in a temporary nan. The min() ignores that NaN so the overall statement produces a non-nan value. You have the same difficulty at

    Lth=min(sqrt(x(1)*a/abs(x(2))),x(1)/pi);

After a few iterations, x(3) gets to about +15, and at that time, the expression (((exp(e2/x(3))-1)^2)) in the definition of Cv reaches about 10E+300; with the other items being multiplied, you overflow to infinity there in a temporary result, but Cv itself does not go infinity or nan. Not yet. But by the time x(3) gets down to about 7.47 then Cv does go completely to NaN because of multiple inf in the expression. After that, the rest of your calculations are toast.

[The.Hossein]

thanks for the quick reply

i checked the global variable issue and also the the lth\ld issue you mentioned. but as it turns out they were not whats wrong with my code. about the last part, as far as i am getting from your explanation, you think something may be is wrong with the definition of Cv. if you look up the formula for Heat Capacity at Constant volume for Vapor, you are gonna get the equation i put there in Cv. so i think its not the issue. i also you used the same Cv in another Code but with a little difference, in that code x(4) was constant (the initial value was assigned for all the steps), and worked like a charm in there, adding the varying x(4) somehow messes up the calculations in matlab. and im not sure why. its just another equation coupled with the rest, shouldn't be an issue i think :/ but something's clear ! i messed something up, and i can't find it

Answer 2

[Torsten]

Check xdot at the beginning of the integration. My guess is that they result in NaN, Inf or something similar.

Note that you divide by x(2) which has zero initial condition.

Best wishes

Torsten.