Matlab ODEs: unable to meet integration tolerances

Red Star

unread,

Oct 26, 2008, 9:20:02 AM10/26/08

to

My ODEs system is:

function dy=odesystem(s,y)
global U_inf rho_inf rho T_inf
dy = zeros(4,1); % a column vector

beta=0.4;
Cd_for=1.2;
g=9.81; %g [m/s^2]
k=0.021;
mu=1e-5;
Cp=1.7646e3;
d=sqrt(4*y(3)/pi);
Re=rho*y(2)*d/mu;
Pr=Cp*mu/k;
Nu=0.023*Re^(4/5)*Pr^0.4;

h=d;
E_star=0.2*(s/d)^1.37*(y(2)/U_inf)^(-0.6);
E=y(3)/(pi*d)*rho_inf*E_star*(y(2)-U_inf);
q_inf=0.5*rho_inf*U_inf^2;

% ODEs
dy(1)=(Cd_for*q_inf*h*(sin(y(1)))^2+...
g*y(3)*(rho-rho_inf)*cos(y(1))...
+E*U_inf*sin(y(3)))/...
(-rho*y(3)*y(2)^2);
d_alfa=(Cd_for*q_inf*h*(sin(y(1)))^2+...
g*y(3)*(rho-rho_inf)*cos(y(1))...
+E*U_inf*sin(y(3)))/...
(-rho*y(3)*y(2)^2);
dy(2)=(q_inf*y(3)*sin(y(1))*cos(y(1))*d_alfa...
-g*y(3)*(rho-rho_inf)*sin(y(1))+...
E*U_inf*cos(y(1))-...
pi*h*rho*beta*(y(2)-U_inf*cos(y(1)))^2-...
y(2)*E)/(rho*y(2)*y(3));
d_v= (q_inf*y(3)*sin(y(1))*cos(y(1))*d_alfa...
-g*y(3)*(rho-rho_inf)*sin(y(1))+...
E*U_inf*cos(y(1))-...
pi*h*rho*beta*(y(2)-U_inf*cos(y(1)))^2-...
y(2)*E)/(rho*y(2)*y(3));
dy(3)=(E-y(3)*rho*d_v)/(rho*y(2));
d_A=(E-y(3)*rho*d_v)/(rho*y(2));
dy(4)=(E*Cp*y(4)+Nu*d/k*(T_inf-y(4))...
-rho*y(2)*Cp*y(4)*d_A-rho*y(3)*Cp*y(4)*d_v)...
/(rho*y(3)*Cp*y(2));

I call this function from my main script:

global U_inf rho_inf rho T_inf;
U_inf=9.67;
rho_inf=5.84;
rho=52.40;
[S,Y]=ode45(@odesystem,[0 12e-3],[pi/2 17.99 1.96e-7 1000]);

Matlab output is:
Warning: Failure at t=2.174547e-003. Unable to meet integration
tolerances without reducing the step size below the smallest
value allowed (6.938894e-018) at time t.

How can solve my system?
An approximate solution is sufficient and I can wait Matlab
processing...

Steven Lord

unread,

Oct 26, 2008, 11:00:03 PM10/26/08

to

"Red Star" <matteo.d...@gmail.com> wrote in message
news:a736b687-c5ee-40ed...@b1g2000hsg.googlegroups.com...
> My ODEs system is:

*snip code*

> I call this function from my main script:
>
> global U_inf rho_inf rho T_inf;
> U_inf=9.67;
> rho_inf=5.84;
> rho=52.40;
> [S,Y]=ode45(@odesystem,[0 12e-3],[pi/2 17.99 1.96e-7 1000]);
>
> Matlab output is:
> Warning: Failure at t=2.174547e-003. Unable to meet integration
> tolerances without reducing the step size below the smallest
> value allowed (6.938894e-018) at time t.
>
> How can solve my system?
> An approximate solution is sufficient and I can wait Matlab
> processing...

Often that warning means that your system of ODEs is stiff. See section 7.9
of the ODE chapter of Cleve's "Numerical Computations with MATLAB" for more
information about stiffness:

http://www.mathworks.com/moler/chapters.html

Try a stiffer solver, like ODE15S or ODE23S.

http://www.mathworks.com/access/helpdesk/help/techdoc/ref/ode113.html

--
Steve Lord
sl...@mathworks.com

Allamarein

unread,

Oct 29, 2008, 9:20:09 AM10/29/08

to

> Often that warning means that your system of ODEs is stiff. See section 7.9
> of the ODE chapter of Cleve's "Numerical Computations with MATLAB" for more
> information about stiffness:
>
> http://www.mathworks.com/moler/chapters.html
>
> Try a stiffer solver, like ODE15S or ODE23S.
>
> http://www.mathworks.com/access/helpdesk/help/techdoc/ref/ode113.html
>
> --
> Steve Lord
> sl...@mathworks.com

Thank you Steve.
ODE15s and ODE23s but I get a warning similar to the previous one.
(i.e.: Using ODE23s:"Warning: Failure at t=3.528573e-003. Unable to

meet integration tolerances
without reducing the step size below the smallest value allowed

(1.253601e-017)at time t.")

I'm going to read Cleve's chapter....it seems interesting.
Writing my ODE system in dimensionless variables could help Matlab to
find a solution?

Torsten Hennig

unread,

Oct 29, 2008, 10:21:16 AM10/29/08

to

I suspect that the exponentiation of solution variables
which may become negative during integration cause the difficulties you encounter
e.g. d=sqrt(...), Re^(4/5) etc.

Best wishes
Torsten.

Allamarein

unread,

Oct 29, 2008, 12:06:06 PM10/29/08

to

> Torsten.

You should be right,Torsten.

I can define an average Nu and Re, but I need to define
E_star=0.2*(s/d)^1.37*(y(2)/U_inf)^(-0.6);

How can I bypass this trouble?
I'm sure someone solved this system by a Range-Kutta, but I don't know
which solver he used and which ICs and costants he had.

Torsten Hennig

unread,

Oct 30, 2008, 3:39:38 AM10/30/08

to

The dirty way ?
Write the yi's in local variables.
Before entering the equations with the local variables,
check if the relevant ones are not negative.
If they are negative, set them to a small positive value.
Check also if the time-derivative for the negative
variables becomes positive again. Otherwise, they will
become even more negative after the next time step.

Best wishes
Torsten.