ode45 failing to produce a sensible solution

15 views (last 30 days)
Hossein
Hossein on 25 Sep 2014
Commented: Mike Hosea on 30 Sep 2014
i have a code to solve a set of nonlinear ordinary differential equations for me, it has 3 different control parameters to simulate the behavior of system, so far everything is fine, when i set one of the control parameters over 50*10^3 system behavior makes no sense and Matlab fails to solve the equations correctly somehow, results coming out of matlab's computations makes no sense comparing to lab results i checked the code and equation derivations dozens of times and couldn't find anything wrong with them,and i also used other ode functions, non worked
ill be more than happy if you share your knowledge on such issue if you've faced before
cheers
here is my code
function dy=thermal_bubble3(t,y)
dy=zeros(4,1);
d=998;
dg=23;
c=1500;
global r0
r0=10^-5;
global s
s=0.0725;
k=0.016;
u=0.000001;
p0=10^5;
global pa;
global f;
cv=1.5*10^3;
a=2.338*10^-5;
Tinf=300;
g=1.31;
global pv
pv=2.33*10^3;
dy(1)=y(2);
dy(2)=((1+y(2)/c)*((y(4)-(2*s/y(1))-(4*u*y(2)/y(1))-p0+pv+-pa*sin(2*pi*f*t))/d)+(g*y(4)*(-3*y(2)/y(1))+(g-1)*(3*k/y(1))*(y(3)-Tinf)*sqrt(3*(g-1)*abs(y(2))/(a*y(1)))+2*s*y(2)/((y(1))^2)-4*u*((y(2)/y(1))^2)-2*pi*f*pa*cos(2*pi*f*t))*(y(1)/(d*c))-1.5*((y(2))^2)*(1-y(2)/(3*c)))*((y(1)*(1-y(2)/c))+4*u/(d*c))^(-1);
dy(3)=((3*(y(1))^2)/(dg*cv*(r0^3)))*((k*(y(3)-Tinf)/(min(y(1)/pi,sqrt((a*y(1))/abs(y(2))))))-y(4)*y(2));
dy(4)=g*y(4)*(-3*y(2)/y(1))+(g-1)*(3*k/y(1))*(y(3)-Tinf)*sqrt(3*(g-1)*abs(y(2))/(a*y(1)));
and here is my solver
clear;
global pa f pv r0 s
pa=50*1000;
f=250000;
r0=10*10^-6;
[t,y]=ode45(@thermal_bubble3,(0:0.001/f:20/f),[r0 0 300 1.01*10^5+2*(s/r0)+pv]);
plotmatrix(t,y);
  1 Comment
John D'Errico
John D'Errico on 28 Sep 2014
I think it is time for you to start learning about other tools in matlab than ODE45. Odds are your problem is stiff when you set that parameter too large. So learn about stiff solvers.

Sign in to comment.

Accepted Answer

Jan
Jan on 26 Sep 2014
It is hard to guess if there is a problem at all, because the question is vague only. It is not clear which parameter is meant and of course we cannot debug the formula without any further knowledge.
But you use the function min and abs, which are not smooth, but step-size controlled ODE solvers need smooth functions. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047.
  3 Comments
Hossein
Hossein on 29 Sep 2014
thanks for the advices, they helped the code to run much smoother about ODE45 is there any alternative in Matlab to solve nonlinear unsteady coupled ODEs ? a parallel processing would be a bless !

Sign in to comment.

More Answers (1)

Mike Hosea
Mike Hosea on 29 Sep 2014
I somewhat disagree with Jan's conception of the limitations of the code. They should be able to handle non-smoothness. Jump discontinuities and functions like min and abs will slow them down. It will be faster if the cross-over points were known in advance and the integrator called to integrate from one point to the next. Now they can't integrate over singularities, obviously. I reformulated the problem like so
function [t,y] = odeprob
r0 = 10^-5;
s = 0.0725;
pv = 2.33*10^3;
pa = 50*1000;
f = 250000;
ic = [r0,0,300,1.01e5+2*(s/r0)+pv];
dydt = @(t,y)thermal_bubble3(t,y,r0,s,pa,f,pv);
tinterval = 0:0.001/f:20/f;
[t,y] = ode45(dydt,tinterval,ic);
plotmatrix(t,y);
function dy=thermal_bubble3(t,y,r0,s,pa,f,pv)
% Constants
d = 998;
dg = 23;
c = 1500;
k = 0.016;
u = 0.000001;
p0 = 10^5;
cv = 1.5e3;
a = 2.338e-5;
Tinf = 300;
g = 1.31;
% dydt
dy=zeros(4,1);
dy(1) = y(2);
dy(2) = ((1+y(2)/c)*((y(4)-(2*s/y(1))-(4*u*y(2)/y(1))-p0+pv+-pa*sin(2*pi*f*t))/d)+(g*y(4)*(-3*y(2)/y(1))+(g-1)*(3*k/y(1))*(y(3)-Tinf)*sqrt(3*(g-1)*abs(y(2))/(a*y(1)))+2*s*y(2)/((y(1))^2)-4*u*((y(2)/y(1))^2)-2*pi*f*pa*cos(2*pi*f*t))*(y(1)/(d*c))-1.5*((y(2))^2)*(1-y(2)/(3*c)))*((y(1)*(1-y(2)/c))+4*u/(d*c))^(-1);
dy(3) = ((3*(y(1))^2)/(dg*cv*(r0^3)))*((k*(y(3)-Tinf)/(min(y(1)/pi,sqrt((a*y(1))/abs(y(2))))))-y(4)*y(2));
dy(4) = g*y(4)*(-3*y(2)/y(1))+(g-1)*(3*k/y(1))*(y(3)-Tinf)*sqrt(3*(g-1)*abs(y(2))/(a*y(1)));
I solved with ODE23, ODE45, and ODE113 and compared results. There are some shapes in the generated results that remind me of the cotangent and secant functions, which leads me to wonder if the problem doesn't have singularities in the interval, or at least that the limiting case of some parameter values is singular where those rapid changes in position or first derivative occur. Be that as it may, what I also see is an warning in every case like
Warning: Failure at t=6.823945e-05. Unable to meet integration tolerances without reducing the step size below the smallest value
allowed (2.168404e-19) at time t.
> In ode45 at 308
In odeprob at 10
One can't just ignore a warning like that. It means that the integration was NOT successful. In that case, results that make no sense past that point are to be expected. You have been warned...literally.
I also see signs that all the integrators are trying to cope with the instability of the underlying problem. If it were instability of the method then we might find another method, but instability of the problem is a fundamental issue that cannot really be side-stepped by a superior numerical method. If the problem itself is unstable, then the codes can be used over short distances (despite mild non-smoothness), but integration for very long is hopeless because of the instability. This would be true whether the problem were smooth or not. Such a problem must be approached in another mathematical way.
  2 Comments
Mike Hosea
Mike Hosea on 30 Sep 2014
Some methods can integrate farther than others before failing, but all numerical initial value solvers will stray from the correct, mathematical solution curve eventually. The problem is that even if the numerical method is theoretically exact, relative errors on the order of the machine epsilon are inevitable. If your ODE takes a perturbation on the order of 1e-16*y in the solution component y and magnifies it considerably, then the problem cannot be solved accurately by any numerical initial value solver for ODEs. This is a general principle.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!