ode45 failing to produce a sensible solution

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);

[John D'Errico]

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.

Answer 1

[Jan]

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.

[Hossein]

this system of odes describe a system with HIGH nonlinear and unsteady properties (a micro bubble in a water solution under the effect of ultrasound waves, Pa is driving pressure properties of the ultrasound wave) the nature of the problem and system is highly nonlinear and highly unsteady i'm guessing i need to use extremely small steps which i'm correctly running on a code, to calculate a 4ms response of system i have 10^16 steps, on one system i'm running the above code with the option i just mentioned and on another system im running the same code and the forementioned option with more error control properties, the problem im facing with this code is, it is obviously gonna take a lot of time to calculate the results. do you have any advice on making my code more efficient ? about abs and min syntax you mentioned above i'm aware they are not smooth functions but with my current calculations and derivings unfortunately i dont have any other option to replace them with, and another question which is bugging my mind a lot ! is it possible to run ode codes on multiple cores to enhance their performance? speaking of hardware i only found out that i can turn intel hyper threading technology off in the bios to get the best out of one core, but is it possible to use more cores ?

some of my question may seem dumb ! but im kind of new to matlab thanks a ton for your advice !

[Jan]

Some general hints:

• There is no reason to use a time consuming global for r0 . Avoid the other globals also. If you need to define parameters externally, see Answers: Anonymous functions for parameters
• Avoid 10^-5, which is an expensive calculation, while 1e-5 is a constant. r0^3 wastes time also, because this is a constant.
• x^-1 is much slower than 1/x for the same reason, at least in older Matlab versions.
• Avoid repeated calculations of e.g. 2*pi*f*t and y(2)/y(1). Evaluate such terms once and store them in a variable.

But the main problem remains: ODE45 is simply not designed to handle non-smooth functions. The results cannot be trustworthy, although many scientists simply ignore this details, because sometimes the results look fine. But the result of an integration is not reliable without an analysis of the sensitivity at all - the difference of the results, when the inital values and parameters are varied.

Matlab's ODE integrators will not produce reliable results for non-smooth function. After 10e16 steps, the accumulation of round-off errors will dominate the solution, such that the result is more or less random.

[Hossein]

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 !

Answer 2

[Mike Hosea]

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.

[Hossein]

yeah mike that warning is my problem also tried much much shorter steps and error control options but i couldn't make it go away i reached the same conclusion that ODE codes in Matlab wont be able to solve my equations, thank for your advice anyway, but do you know another way of solving such equations in my case nonlinear unstable coupled ODEs ?

[Mike Hosea]

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.