[SciPy-Dev] ODE solvers

[Anne Archibald]

On Thu, Jun 4, 2015 at 1:52 PM Sturla Molden <sturla.molden at gmail.com>
wrote:

> Anne Archibald <archibald at astron.nl> wrote:
>
> > * Manage reentrancy constraints (some FORTRAN solvers store state in
> global
> > variables)
>
> Use a threading.Lock instead of the GIL.
>

Right. And also copy all the relevant information to some kind of local
storage between calls.


> > 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.
>
> A suggestion:
>
> Provide a Cython cdef class as permitted RHS callback. Users can then
> cimport and subclass this Cython class in their own Cython module. There
> will only be some tiny overhead (a vtab lookup). The superclass' callback
> method should just be "pure virtual" and raise NotImplementedError. We will
> have to provide a .pxd similarly to what we now do for sharing BLAS and
> LAPACK. It is also possible to have another Cython callback class which
> assumes nogil during the ODE solver loop.
>

This seems like a fairly clean and flexible solution. integrate.quad also
seems to have some kind of arrangement involving ctypes? I guess you could
just have a particular subclass that called via ctypes.


> Another thing:
>
> We should have more solvers, even simple ones like Euler and 4th order
> Runge-Kutta for reference. We also lack Bulirsch-Stoer. Ot seems we just
> have those that are in odepack.
>

Why the simple ones? Is the idea to provide some baseline solver for
comparisons? The design I was proposing only makes sense for solvers with
dense output, at least, and preferably adaptive step-size control. It's not
that hard to add basic dense output to simple solvers using
KroghInterpolator, but if the goal is educational this might be unwanted
complexity. Still, a clean pure-python RK45 with adaptive step-size control
and dense output wouldn't hurt.

That said, a good B-S integrator would be valuable, since there doesn't
seem to be anything too suitable for high-accuracy solution of smooth
problems. It also doesn't seem too complicated to implement, even in dense
form. That plus the Dormand-Prince ones plus ODEPACK should cover most
problems (apart from symplectic).

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