getting NaN in ODE45

8 views (last 30 days)
Hossein
Hossein on 18 Oct 2014
Edited: Star Strider on 18 Oct 2014
i have a set of nonlinear odes to solve but im getting NaN in output matrix r
i dont know what is wrong with my code, it should work,
ill be more than happy if you share your knowledge on such issue if you've faced before
cheers
here is my function and solver
function dr=Thermal_Bubble2(t,r)
dr=zeros(3,1);
global pa f r0 T0; a=0.150e-6;
l=0.0184;
kb=1.38e-23;
ps=1e5;
c=1500;
d=998;
s=0.0725;
cv=1460;
u=0.000001;
pv=2.33e3;
w=2*pi*f;
b=5.1e-29;
%Equations n=4*pi*(r0^3)*(ps+2*(s/r0)+pv)/(3*kb*T0+(ps+2*(s/r0)+pv)*b) ;
lth=min(r(1)/pi,sqrt(a*r(1)/abs(r(2))));
A=4*pi*(r(1)^3)-3*n*b;
pg=(3*n*kb*r(3))/A;
dr(1)=r(2);
dr(2)=(((1+r(1)/c)/d)*(pg-(2*s/r(1))-(4*u*r(2)/r(1))-ps+pv-pa*sin(w*t))+(r(1)/(d*c))*(((3*n*kb*((4*pi*((r(1))^2)/cv)*((l*(T0-r(3))/lth)-pg*r(2))))/A)-(36*pi*(r(1)^2)*r(2)*n*kb*(r(3))/(A^2))-(2*s*r(2)/(r(1)^2))-(4*u*(((r(2))/(r(2)))^2))-w*pa*cos(w*t))-(1.5*((r(2))^2)*(1-(r(2))/(3*c))))/((r(1))*(1-(r(2))/c));
dr(3)=(4*pi*((r(1))^2)/cv)*((l*(T0-r(3))/lth)-pg*r(2));
solver
clear;
global f pa r0 T0
f=25e3;
pa=50e3;
r0=10e-6;
T0=300;
[t,r]=ode45(@Thermal_Bubble2,(0:0.01/f:20/f),[r0 0 T0]);
plotmatrix(t,r);

Answers (1)

Star Strider
Star Strider on 18 Oct 2014
Edited: Star Strider on 18 Oct 2014
I believe the problem is here:
dr(2)=(((1+r(1)/c)/d)*(pg-(2*s/r(1))-(4*u*r(2)/r(1))-ps+pv-pa*sin(w*t))+(r(1)/(d*c))*(((3*n*kb*((4*pi*((r(1))^2)/cv)*((l*(T0-r(3))/lth)-pg*r(2))))/A)-(36*pi*(r(1)^2)*r(2)*n*kb*(r(3))/(A^2))-(2*s*r(2)/(r(1)^2))-(4*u*(((r(2))/(r(2)))^2))-w*pa*cos(w*t))-(1.5*((r(2))^2)*(1-(r(2))/(3*c))))/((r(1))*(1-(r(2))/c));
NOTE: r(2)/r(2) in the expression ‘(4*u*(((r(2))/(r(2)))^2))-w*pa*cos(w*t))’
Since ‘r(2)’ is initially zero, this will cause the entire expression to be NaN, which will then propagate throughout your calculations.

Tags

Community Treasure Hunt

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

Start Hunting!