# False divergence of the NDSolve solution: how to avoid

Discussion in 'Mathematica' started by Alexei Boulbitch, Jul 16, 2009.

1. ### Alexei BoulbitchGuest

Dear Community,

I am simulating a system of ODE using v6. Here are the equations:

eq1 = x'[t] == y[t];
eq2 = y'[t] == 1/x[t] - 1.4 - (4.5 + y[t])*(1 + z[t]^2);
eq3 = z'[t] == 18*z[t] - 0.75*(4.5 + y[t])^2*z[t] - z[t]^3;

It is simulated at x>0. This system at x>0 seems to be globally stable.
To understand it observe that at large x, y, and z one finds
y' ~ - y*z^2 and z' ~ - z^3. In other words, there is a kind of a
non-linear "returning force" for y and z, while x follows the dynamics
of y.

However, when solving it on Mathematica I sometimes find trajectories
that counterintuitively diverge.
Check this for example:

NDSolve[{eq1, eq2, eq3, x == 0.669, y == 0.881,
z == 0.988}, {x, y, z}, {t, 0, 40}];

Plot[{Evaluate[x[t] /. s], Evaluate[y[t] /. s],
Evaluate[z[t] /. s]}, {t, 0, 45}, PlotRange -> All,
PlotStyle -> {Red, Green, Blue},
AxesLabel -> {Style["t", 16], Style["x, y, z", 16]}]

My guess is that this is due to some peculiarity in the numeric method
used, and the method should be probably changed, or its parameters
specified. I am however, not experienced in numeric approaches for
solving ODEs.

Now comes the question:
Can you give me a hint, of
(i) what may be the reason of such a behavior?
and
(ii) What should I do to avoid such a false divergence?

Thank you, Alexei

--
Alexei Boulbitch, Dr., habil.
Senior Scientist

IEE S.A.
ZAE Weiergewan
11, rue Edmond Reuter
L-5326 Contern
Luxembourg

Phone: +352 2454 2566
Fax: +352 2454 3566

Website: www.iee.lu

This e-mail may contain trade secrets or privileged, undisclosed or otherwise confidential information. If you are not the intended recipient and have received this e-mail in error, you are hereby notified that any review, copying or distribution of it is strictly prohibited. Please inform us immediately and destroy the original transmittal from your system. Thank you for your co-operation.

Alexei Boulbitch, Jul 16, 2009

2. ### Roland FranziusGuest

y is a dummy. You have to consider eigenvalues of the matrix for
|x,y,z|->oo. At |(x,z)| -> oo easily you will find by elimination of y

x'' == - x' z^2
z' == 18 z - x'^2 z - z^3

If, by chance of rounding errors z < 0 z will run away to -oo.
s=NDSolve[{eq1, eq2, eq3, x == 0.669, y == 0.881,
z == 0.988}, {x, y, z}, {t, 0, 40}];

NDSolve::mxst: Maximum number of 10000 steps reached at the point t == \
27.928120676960916`. >>

The Plot explodes at t=27 exponentially
Its always extremely difficult to control solutions of osciallating
nonlinear systems for a long time:-(

At least in this case use Abs[z[t]] on the right to cross the line
immediately if an error of z<0 is occurs or choose a t-step adaptive
method if z is approaching 0.

Roland Franzius, Jul 17, 2009

3. ### seanGuest

I think this arises from NDSolve not having the solutions beyond the
point where it reaches the max number of steps. You can increase the
number of steps NDSolve takes to a higher number of you can just set
it to infinity.(obviously is not a good idea if your system is
weird...)

eq1=x'[t]==y[t];
eq2=y'[t]==1/x[t]-1.4-(4.5+y[t])*(1+z[t]^2);
eq3=z'[t]==18*z[t]-0.75*(4.5+y[t])^2*z[t]-z[t]^3;

DynamicModule[{},
Manipulate[
s=NDSolve[{eq1,eq2,eq3,x==0.669,y==0.881,z==0.988}, =
{x,y,z},
{t,0,tf}, MaxSteps->Infinity];
Plot[{Evaluate[x[t]/.s], Evaluate[y[t]/.s], Evaluate[z[t]/.s]}, {t,
0,tf}, PlotRange-> All, PlotStyle->{Red,Green,Blue}, AxesLabel->{Style
["t",16], Style["x, y, z",16]}], {tf, 10, 100}]
]

wise confidential information. If you are not the intended recipient and ha=
ve received this e-mail in error, you are hereby notified that any review, =
copying or distribution of it is strictly prohibited. Please inform us imme=
diately and destroy the original transmittal from your system. Thank you fo=

sean, Jul 18, 2009
4. ### Roland FranziusGuest

Since x occurs in eq2 this is not correct. Elimination of dummy y gives

x'' ~ x' z^2 and z'= -z^3

but the problem is at x=0

x'' = 1/x
z' = z( 18 -0.75(4.5+x') -z^2)

If by error you have x<0 and x' <4.5 you have a problem.

#############

After a hint of Alois Steindl I reconsidered the full second order
system once again.

It is possible to integrate beyond t=28 if x is restricted to x>0 somehow.

The second order system

eq1 = x''[t] == 1/x[t] - 1.4 - (4.5 + x'[t])*(1 + z[t]^2)
eq2 = z'[t] == 18*z[t] - 0.75*(4.5 + x'[t])^2 *z[t] - z[t]^3

explodes numerically somewhere at t=28

s = NDSolve[{eq1, eq2, x == 0.669, x' == 0.881,
z == 0.988}, {x, z}, {t, 0, 28.5}];

To see the point of explosion at a logarithmic scale look at

Plot[{Evaluate[ArcSinh[x[t] /. s]], Evaluate[ArcSinh[z[t] /. s]]}, {t,
27., 28.1}, PlotRange -> All,
PlotStyle -> {Red, Green},
AxesLabel -> {Style["t", 16], Style["x, y, z", 16]}]

The reason is x becoming negative and a moment later z follows.

This is analytically impossible because of the hard reflection at x=0 by
1/x. Numerically of course a big z-term with x'<-4.5 may override it.

A try with absolute values

eq1 = x''[t] == 1/x[t] - Abs[1.4 - (4.5 + x'[t])*(1 + z[t]^2)]
eq2 = z'[t] == 18*z[t] - 0.75*(4.5 + x'[t])^2 *z[t] - z[t]^3

s = NDSolve[{eq1, eq2, x == 0.669, x' == 0.881,
z == 0.988}, {x, z}, {t, 0, 40}];

Plot[{Evaluate[x[t] /. s], Evaluate[z[t] /. s]}, {t, 0, 3},
PlotRange -> All,
PlotStyle -> {Red, Green, Blue},
AxesLabel -> {Style["t", 16], Style["x, y, z", 16]}]

The apparent singularities at x=0 is caused by the hard reflection of x
at x=0 by the non-Lipshitz x'' = 1/x force term.

The Abs-approach seamingly eliminates the problem for the given starting
values. If it is working for all starting values, who knows.

Roland Franzius, Jul 19, 2009