[SciPy-Dev] ODE solvers

Xavier Barthelemy xabart at gmail.com
Thu Jun 4 20:32:48 EDT 2015


Dear Scipy contributors,

I just would like to make a side step of this discussion.
I have worked and badly implemented multiple classical method to solve ode.
my point is that I can provide (with tests) the coefficients and the method
for
RK31, RK32, RK41,RK42, RK5, RK61,RK62, RK7, RK8 from the book and
discussions with John Butcher (http://jcbutcher.com/d/)

and
I derived predictor-corrector coefficients for Adams-Bashforth and
Adams-Bashforth-Moulton for non uniform (time) steps. it also comes with
tests.
I have predictor AB and corrector ABM for 3, 4, 5 and 6 non uniform time
steps.

The main difficulty is not the implementation, but finding the info and
testing it.
I can provide that if needed and/or interested

all the best
Xavier




2015-06-05 1:25 GMT+10:00 Sturla Molden <sturla.molden at gmail.com>:

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



-- 
 « Quand le gouvernement viole les droits du peuple, l'insurrection est,
pour le peuple et pour chaque portion du peuple, le plus sacré des droits
et le plus indispensable des devoirs »

Déclaration des droits de l'homme et du citoyen, article 35, 1793
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20150605/8c0856ce/attachment.html>


More information about the SciPy-Dev mailing list