[SciPy-user] better control over odeint

A. M. Archibald peridot.faceted at gmail.com
Sat Aug 12 01:55:24 CDT 2006


Hi,

I'm writing a program that requires extensive use of
scipy.integrate.odeint (tracing paths around a black hole). It mostly
works well but has certain frustrating features.

(1) It is apparently impossible to prevent it from printing to the
standard output. This is wrong; if it must output error messages, they
should go on standard error, so as not to contaminate the output of my
program. It would also be nice to be able to turn them off! Or, more
specifically, to signal the error conditions in a way I can reasonably
test for. String matching on info["messages"] is not at all
satisfactory.

(2) It is very difficult to control what happens when an error occurs.
When a function called from inside odeint throws an exception
(undefined variable, overflow error, what have you), odeint catches
it, prints a backtrace (to standard error, at least, thankfully), and
continues trying more points. Since what has usually happened is that
odeint has run into a singularity, it grinds away endlessly, spewing
page after page of the same error message, before finally giving up.
An error is then signaled by leaving info["mused"][n]==0.

Ideal for my purposes would be for odeint simply not to catch
exceptions, or to treat exceptions as an instruction to return early.
Failure to integrate over the whole range of t could be signaled by
raising an exception. This would allow easy catching of integration
failures, and it would also (for example) allow me to signal that the
ray had proceeded into a forbidden region by throwing an exception.

I realize that nobody is very keen on rooting around inside the
FORTRAN code for odeint, and I can't see directly how to make it stop
immediately. It should be fine for the user function to simply return
NaNs.

Another option is, if the integration falls apart at some point, to
simply stop and provide information about how far the integration
successfully progressed. lsoda.f is definitely capable of this without
modification.

Thanks.
A. M. Archibald


More information about the SciPy-user mailing list