       Re: False divergence of the NDSolve solution: how to avoid

• To: mathgroup at smc.vnet.net
• Subject: [mg101842] Re: False divergence of the NDSolve solution: how to avoid
• From: Roland Franzius <roland.franzius at uos.de>
• Date: Sun, 19 Jul 2009 07:14:14 -0400 (EDT)
• References: <h3n5uh\$2gp\$1@smc.vnet.net> <h3pesm\$oib\$1@smc.vnet.net>

```Roland Franzius schrieb:
> Alexei Boulbitch schrieb:
>> 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.

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.

>> 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?
>
> 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.
>

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

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

```

• Prev by Date: Re: Player Pro Limitations
• Next by Date: Re: eigenvector centrality
• Previous by thread: Re: False divergence of the NDSolve solution: how to avoid
• Next by thread: eigenvector centrality