# ode45 failing to produce a sensible solution

15 views (last 30 days)

Show older comments

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
on 28 Sep 2014

### Accepted Answer

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.

### More Answers (1)

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
on 30 Sep 2014

### See Also

### Community Treasure Hunt

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

Start Hunting!