implement precise ODE solvers

[GoogleCodeExporter]

I'd like to have something to solve any ODE, where I'd tell it:

I want 15 digits (or 100 digits) and I want them all right.

And the algorithm needs to automatically determine the steps size and the
precision, so that the result is correct to 15 (or 100) digits, whatever I
specify.

Original issue reported on code.google.com by ondrej.c...@gmail.com on 24 Mar 2008 at 11:23

[GoogleCodeExporter]

I think the only way to do this efficiently is to use compute a high order 
Taylor
polynomial with which a relatively large step size can be used.

Googling turned up something called the Parker-Sochacki method; possibly that is
something that can be used.

Original comment by fredrik....@gmail.com on 24 Mar 2008 at 12:26

[GoogleCodeExporter]

I've started implementing a custom Taylor series based high precision ODE 
solver.

It uses Euler's method with extremely small step size and extremely high 
precision to
generate a Taylor polynomial valid over a large interval. Euler's method is 
then used
again at the endpoint of the interval to generate a second Taylor polynomial, 
etc.
The Taylor polynomials are computed as needed and cached on the fly. The 
approach
makes it really easy to generate a solution function that can be evaluated 
anywhere;
you don't have to provide all function values at the start.

Here is a test with the harmonic oscillator equation:

from time import clock
mp.dps = 30
f = odefun(lambda x, y: (y[1], -y[0]), 0, [0, 1], degree=30)
t1=clock(); q=f(20); t2=clock(); t2-t1  # Fill in cache
t1=clock(); q=f(20); t2=clock(); t2-t1  # Cached
print q[0], q[1]
print sin(20), cos(20)

Output:
1.633488613776332
0.0021075304263531436
0.912945250727627654376099983846 0.408082061813391986062267860928
0.912945250727627654376099983846 0.408082061813391986062267860928

So computing the data required for 30-digit accuracy on the interval [0, 20] 
takes
1.6 seconds, and after that, function evaluation for any 0 < x < 20 takes only 
a few
milliseconds. I think both results are excellent. And if you want less than 
15-digit
accuracy, it could be made *really* fast by converting the Taylor coefficients 
to
floats once they have been computed.

There is a lot of work to do still:

* The automatic choice of degree and step size needs to be tuned
* The determination of accuracy radius for the Taylor series needs to be made 
robust
* The degree of the Taylor polynomial could be made adaptive
* It should be fixed to handle negative values of x
* It should also be fixed to handle complex values of x (this can be done by 
taking
pure real steps followed by pure imaginary steps, so you get caching over a 
grid in
the complex plane)
* The algorithm and the degree/step/radius settings need to be tested for a 
large set
of differential equations

I'm not an expert on numerical ODEs, so if anyone is interested in looking at 
any of
the above, feel encouraged to do so.

Original comment by fredrik....@gmail.com on 13 Nov 2008 at 7:03

Attachments:

• taylorode.py

[GoogleCodeExporter]

Original comment by fredrik....@gmail.com on 13 Nov 2008 at 7:10

[GoogleCodeExporter]

This is excellent!

I currently don't need this, but let's get this in mpmath. The next time I need 
to
solve some ODE, I'll look at that.

Original comment by ondrej.c...@gmail.com on 13 Nov 2008 at 7:19

[GoogleCodeExporter]

Committed in r783, docs at
http://mpmath.googlecode.com/svn/trunk/doc/build/calculus/odes.html

Now this just needs some more stress testing, and the RK4 solver should to be
integrated in the same framework somehow.

Original comment by fredrik....@gmail.com on 26 Nov 2008 at 11:44