[SciPy-Dev] ODE solvers

[Anne Archibald]

Hi folks,

Over on PR 4904 there has been much discussion about the right interface
for ODE solvers. I have some experience with ODE solvers under various
conditions, including one fairly heavy-duty problem, but it would be good
to have a wider discussion about how ODE solvers ought to work from within
scipy.
https://github.com/scipy/scipy/pull/4904

Please let me know if you think this is a good design, or if you think I've
left out anything important. I can probably get a pure-python version
driving one of the existing solvers this weekend.

My suggestion has two layers:

Each solver gets wrapped in an object that does the following:
* Setup and initialization of the solver: RHS, initial conditions, error
tolerance, direction of t, solver parameters
* Manage reentrancy constraints (some FORTRAN solvers store state in global
variables)
* Advance one step (size set by solver, or shorter if requested)
* Know t, y, yprime at beginning and end of current step
* Report an estimate of the accumulated error
* Evaluate y (maybe yprime and higher derivatives) anywhere within the
current step
Most existing scipy solvers, and most modern implementations, support all
of these basic operations.

Then there is one user-level driver object that uses an underlying solver
object to provide a friendlier interface that can:
* Advance in a single step to any requested time.
* Take an ordered array of times and return y values (and/or derivatives of
various orders, perhaps) for each time.
* Return the same for a list of times chosen by the solver to adequately
represent the solution.
* Use additional stopping criteria: for example, stop when any of a list of
"event functions" changes sign or passes one of a list of values of
interest. Some of these would actually stop the solver even if the solver
would otherwise have tried to return more values, others would simply
indicate places to report the solution state.

Concerns for implementation:
* Users should be able to provide compiled callables of some kind and have
them called without passing through the python interpreter. This presumably
needs to apply to stopping conditions/event functions too.
* For fast solutions, a python driver loop may be too slow.

Interesting but out-of-scope ideas:
* Symplectic solvers
* ODE solutions as function objects evaluatable anywhere
* Solvers using long doubles or quad precision internally


For concreteness about the "stopping time"/"event functions", let me
mention two situations I have encountered. The first was I needed to detect
when the ray I was integrating got too close to a black hole and switch
coordinate systems. The second is where I need to integrate a solution not
to a particular list of times, but until a function of the time (time plus
position-dependent propagation delay) hits a particular list of values.


Anne
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20150604/35045b8e/attachment.html>