[SciPy-Dev] ODE solvers

[Yuxiang Wang]

+1 for the RK4 solver. I have manually written some RK4 ode solvers where:

1) In class, RK4 let the students understand more what is happening
while still a maintainable size;
2) I tried to reproduce some results in 60s paper where they used RK4.

Shawn

On Thu, Jun 4, 2015 at 11:25 AM, Sturla Molden <sturla.molden at gmail.com> wrote:
> 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
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev



-- 
Yuxiang "Shawn" Wang
Gerling Research Lab
University of Virginia
yw5aj at virginia.edu
+1 (434) 284-0836
https://sites.google.com/a/virginia.edu/yw5aj/