[SciPy-user] ODE solvers and Scipy

A. M. Archibald peridot.faceted at gmail.com
Sat Aug 12 11:51:05 CDT 2006


On 12/08/06, JJ <josh8912 at yahoo.com> wrote:

> I am just about to start building a somewhat complex
> dynamic model that will incorporate a large  number
> (~50-100) of differential equations and parameters to
> be estimated.  This is my first project of this type,
> so I am on a steep learning curve.  I could build it
> using Matlab's Simulink, but I would rather build it
> in Python.  I have looked for good graphical
> interfaces to ODE solvers in Python but have not found
> any that seem suitable.  I would not mind using the
> Scipy ODE solvers, but I have a couple of quick,
> simple questions before I invest a lot of time in that
> direction.  Any thoughts would be appreciated.
>
> -- The (few) examples I have seen in using Scipy's ODE
> functions are only small, simple ones.  Is there
> anything to prevent use of the solvers for larger,
> more complex systems of ODEs?

Scipy's ODE solver is based on ODEPACK's lsoda.f, which is robust,
well-tested, and well-documented (see, for example,
http://www.netlib.org/alliant/ode/prog/lsoda.f ). All it does is
integrate a system of first-order ODEs over a user-specified range; if
problems occur it carries the integration as far as it can and then
stops. The Python wrapper, however, has a few serious flaws. In an
attempt to provide Numeric-style calling, the ability to stop partway
through an interval appears to have been lost, and any error or
exception (including KeyboardInterrupt) is trapped and an error is
printed to the standard output, rendering debugging extremely
cumbersome.

More seriously, if there are parts of the solution space that cause
problems (the solution approaches a singularity, or the edge of a
coordinate patch), there is no really effective way to bring the
solution up to the edge and then stop. This is because lsoda doesn't
do that - that feature is in lsodar.f (
http://www.netlib.org/alliant/ode/prog/lsodar.f ), which has not
(yet?) been wrapped in python.

I'm actually hoping to get a wrapper for lsodar written, but as I've
no experience with FORTRAN, we'll see.

A. M. Archibald


More information about the SciPy-user mailing list