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] == 0.669, y[0] == 0.881,
>> z[0] == 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] == 0.669, x'[0] == 0.881,
z[0] == 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] == 0.669, x'[0] == 0.881,
z[0] == 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