MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • 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