Open GoogleCodeExporter opened 9 years ago
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
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:
Original comment by fredrik....@gmail.com
on 13 Nov 2008 at 7:10
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
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
Original issue reported on code.google.com by
ondrej.c...@gmail.com
on 24 Mar 2008 at 11:23