Strange behavior of dsolve({ },....,numeric,...

Posted on 2008-07-28 17:41 By lemelinm ( Maple Rating 3

227)

Categories:

How do I ...? with Maple

Hi all,

I have the following ODE

> deq := diff(x(t), t) = 1-t-x(t)/t;

                               d                 x(t)
                       deq := --- x(t) = 1 - t - ----
                               dt                 t  
> ci := x(1) = 0;

                               ci := x(1) = 0

If I solve for the exact solution, I have

> solex := rhs(dsolve({ci, deq}, x(t)));

                                    1  2   1      1 
                         solex := - - t  + - t - ---
                                    3      2     6 t

then I know that there will be a discontinuity at t=0

But if I do the following:

> p := dsolve({ci, deq}, x(t), type = numeric, range = -1 .. 4);

                        p := proc(x_rkf45)  ...  end;
> plots[odeplot](p, color = navy);


I have a plot that is not the right one if I compare with:

>plot(solex,t=-1..4)

What am I missing?

Mario

Strange bug

Comment on 2008-07-29 00:29 from Robert Israel ( Maple Rating 4

1460)

This is a strange and rather alarming bug in dsolve(..., numeric). The strangest part is that it hasn't been found and corrected long ago. Here's a slightly simpler example.

> dsolve({diff(x(t),t)=-x(t)/t, x(1)=2}, numeric)(0.9);

The answer should be essentially

[t = .9, x(t) = 2.22222222222222232]

and it is, up to Maple 7. But from Maple 8 up it is

[t = 0.9, x(t) = 1.90299400449141576]

I would guess that Maple is trying to make some sort of transformation, and ends up solving the wrong differential equation.

Not a transformation bug

Comment on 2008-08-05 11:23 from Allan Wittkopf (

75)

I have looked at this problem in detail, and the same issues can be observed for the simpler problem:

{diff(x(t),t)=x(t)/(1-t), x(0)=1}

which exhibits the same behavior for a singularity at t=1.

First of all, only the rkf45 method appears to have any difficulty with this problem, as choosing method=<anything else> appears to catch the singularity (or at least behave as though the solution is blowing up).

Tracing through the method, the problem is really that the step size chosen for the rkf45 method is dramatically increasing at each time step. The cause appears to be that the Runge-Kutta pair in use for this method dramatically underestimates the error - to such an extent that I believe that errors in equations of the form x'=x/(1-t) are in the effective nullspace of the pair in use.

Recall that the way that an rkf45 method works is to compute a pair of solutions of 4th and 5th order accuracy, then use their difference to estimate the error. The problem occurs when there is little or no difference between these solutions, but there is still significant error. I believe that this is one of those cases.

You can see that choosing a very small initial step size via 'initstep' causes the solution to be accurately computed over a (slightly) larger range.

This suggests the addition of a 'maxstep' option to restrict the maximum step size for a problem to avoid this difficulty.

- Allan Wittkopf

RK

Comment on 2008-08-05 14:07 from Robert Israel ( Maple Rating 4

1460)

Very interesting: the Runge-Kutta formulas of order 2, 3 and 4 (I didn't have the formula for order 5 handy) give the exact solution for this DE. Maybe all RK formulas (for order > 1) do.

Characterizing this behaviour

Comment on 2008-08-05 15:02 from Robert Israel ( Maple Rating 4

1460)

More and more interesting: I think I can characterize all first-order linear DE's that have this property. I'll try to write it up in a blog page.

Additional Information

Comment on 2008-08-05 14:02 from Allan Wittkopf (

75)

OK, this is a very interesting problem...

I have the exact coefficients for the RK pair, and for an equation as simple as x'=x/(1-t) I have obtained an exact form for the single step solution beginning at x(t0)=x0 for step size h as:

x0*(-1+t0)/(-1+t0+h)

and the error estimate is indeed exactly zero.

The pair in use is the well known Fehlberg pair from

Klassische Runge-Kutta-Formeln vierter und niedrigerer Ordnung mit Schrittweiten-Kontrolle und ihre Anwendung auf Waermeleitungsprobleme, Computing, 6 (1970) 61-71

In contrast, the form of the error estimate for x'=-x is given by:

1/780*h^4*x0+1/2080*h^5*x0
So here the error grows monotonically with the stepsize.

So this appears to be an issue with the method, and not the implementation.

- Allan Wittkopf

Final Comment

Comment on 2008-08-05 15:26 from Allan Wittkopf (

75)

OK, now I have all the information.

The real issue is that the solution is non-polynomial, but both the 5th and 4th order solutions compute the answer exactly (thanks for the clue Robert!)

The method, as it is currently programmed, computes the solution from x0 to x0+h, then obtains an interpolant that should approximate the solution over x0..x0+h to 5th order accuracy. The problem is that since we get an exact solution for a non-polynomial input the interpolant is quite wrong, and the solution becomes incorrect for a region.

The really strange thing is that if you integrate past the one bad step, the solution becomes correct again (e.g. run the revised problem with initstep=1e-8, over t=0..10, and you will see a correct solution for 0..0.4 and 3.7...) which shows that it is purely a case where the interpolant is bad, but the solution points are being computed correctly.

- Allan