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: [mg101802] Re: False divergence of the NDSolve solution: how to avoid
  • From: Roland Franzius <roland.franzius at uos.de>
  • Date: Fri, 17 Jul 2009 05:05:44 -0400 (EDT)
  • References: <h3n5uh$2gp$1@smc.vnet.net>

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.

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.

> 
> 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}];

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

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


> 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]}]

The Plot explodes at t=27 exponentially

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

-- 

Roland Franzius


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