[SciPy-Dev] ODE solvers

[Sturla Molden]

On 04/06/15 14:55, Anne Archibald wrote:

> 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.

We can get a function pointer as a void* from ctypes too. An ODE solver 
could accept this too. I would also suggest that we allow an f2py object 
as callback as it has a _cpointer attribute (PyCObject) which stores the 
function pointer in a void*.

So now we have at least four different things that can be allowed as 
callback:

1. Any Python callable object
2. Subclass of a certain Cython cdef class
3. ctypes function object
4. f2py Fortran object

The latter three we can optimize for maximum efficiency.

With ctypes and f2py we need to handle the GIL. It should not be 
released if the ctypes callback is from a PyDLL, but it can be released 
if from a C DLL. And if it is an f2py object we must check to see if the 
callable is threadsafe. (Is this attribute exposed somehow?) With Cython 
it is easier because we can just use a different class with nogil callback.


>     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?

Sometimes I have wondered if I got a weird solution because my ODE did 
not converge correctly or because it should be be this weird.

We do have a dense RK3/8 solver in SciPy, though. It was actually 
published in the same paper as RK4 and does the same thing. RK4 is the 
classical one that is in every textbook. It might not be that those who 
look for RK4 knows that RK3/8 is almost equivalent.


> Still, a clean pure-python RK45 with adaptive step-size
> control and dense output wouldn't hurt.

Yes, similar to ode45 is Matlab. Fehlberg method or whatever it is called.


> 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.

B-S is not difficult to implement.

There is also a C++ implementation in last edition of the "book that 
must not be named". Thus far I am untainted by not having looked at it.


Personally I am mostly interested in solvers for sparse systems, e.g. 
for compartmental models of neurons. But I think all solvers in SciPy 
are intended for dense sets of equations.


Sturla