False divergence of the NDSolve solution: how to avoid

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

  1. 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] == 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?

    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
    #1
    1. Advertisements

  2. 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] == 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`. >>

    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
    #2
    1. Advertisements

  3. Alexei Boulbitch

    sean Guest

    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]==0.669,y[0]==0.881,z[0]==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=
    r your co-operation.
     
    sean, Jul 18, 2009
    #3
  4. 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] == 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, Jul 19, 2009
    #4
    1. Advertisements

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.