[SciPy-user] PyDSTool and SciPy (integrators)
Fri Apr 18 10:34:57 CDT 2008
Thanks to everyone who has offered their thoughts about PyDSTool.
To address some of the points that have been raised in the discussion
1. Best integration method
There is no single best integrator or class of integrator to use;
optimality is determined by the problem and the computing resources
available. Every 10 years or so, people have surveyed the integrator
landscape, tested the performance of the available methods on an
evolving battery of test problems, and come up with heuristics to
guide users who want a black box solution to apply to their ODE/DAE
problems. The Adams-type method used by scipy was considered the best
all-around solution in the 1970s (Hull, et al., SIAM J. Num. Anal.
1972) if function evaluations were relatively expensive; otherwise
Bullirsch and Stoer's extrapolation method was better. Hairer and
Wanner's (H&W) Radau methods and Petzold's extensions to Gear's
version of a variable-order Adams method basically seem to be the best
general purpose methods for stiff systems right now (Cash, Proc. Royal
Soc. 2003 & references). The advantage of H&W's Radau method is that
it provides significantly higher accuracy for the same choices of
tolerances than the DASSL (Petzold, Brennan) codes. The trade off is
that it requires more function calls per step, though it takes fewer
steps. If the function calls are expensive, as they would be if done
through callbacks to Python functions, this can be a significant
For the kinds of problems Rob and I deal with, high accuracy with
reasonable speed is essential, and our function evaluations are not
excessively expensive, hence our preference for Radau. We have found
the HW dopri853 implementation to be fast and robust for non-stiff
(and some moderately stiff) systems and the Radau5 implementation to
be able to handle almost all stiff problems that stymie dopri.
2. Integrator control
I would say the codes provided by H&W provide essentially the same
level of control as the codes in ODEPACK, based on my experience with
the H&W codes and an inspection of the ODEPACK documentation.
3. Stopping events, auxiliary functions
PyDSTool's event functions, which can be set as 'terminal' (halting
integration) or 'non-terminal' (not stopping integration, but saved
and returned to the caller), may be any C function which returns a
double. The event is considered to have happened when the event
function returns 0.0. There is a layer of code between the interface
to Python and the interface to the H&W codes where the events are
checked. After each step of the integration, the event handling
routines check if any event functions have changed sign; if so, the
exact point at which the sign change (and hence event) occurred is
found using bisection and the dense output function provided by the
H&W integrator, to user specified accuracy, which can be set
independently for each event function.
Anne: This feature of PyDSTool would seem to be what you wanted, and
does not require restarting the integration with a smaller step size.
However, it is possible that your problem involves other features that
I didn't take into consideration when designing this part of PyDSTool.
If you could provide us the example, we would find it useful.
At the end of an integration, the intermediate layer of the PyDSTool
integration code will compute any function of the phase space
variables, parameters, etc. you want on the computed trajectory and
return the result. Performing this internally rather than in python
makes it much faster.
Using compiled RHSs is essential for fast integration. I prefer to
keep this feature in; for many applications, the using callbacks to
Python is simply too slow.
That said, I can see a way to modify what I have written to support
(1) an all python callback interface to the H&W integrators,
eventually with event handling, etc.
(2) a version of what PyDSTool has now which allows for a mix of C
functions and callbacks to python functions for events, etc.
Note that (1) will nullify the speed increases we currently get,
particularly for Radau. Also, implementing (1) would likely be much
easier than (2), though both would take me some time.
I think it is better to return an entire trajectory rather than the
endpoint of a
Returning an object is what we do in PyDSTool, but it goes against the
lightweight advantages of the odeint interface as it stands now. I
think our Trajectory class is overkill for many people.
Anne: Providing so much information at every trajectory point would
require much revision of the current code, for which I see little
additional benefit. Could you explain how you would use all of the
additional information? Also, could you elaborate on your "running"
mode? I do not understand what you are suggesting.
I would like to provide an interface for all of the H&W codes that
have dense output. This includes other methods for stiff integration,
a delay differential equation integrator, and possibly a symplectic
integrator. However, these integrators require the systems to be input
with particular forms which do not necessarily conform so nicely to
the odeint interface, for example. I would like to settle on an API/
calling format before moving ahead with this project.
I have explored possibilities for adding some automatic
differentiation and Taylor series integration functionality to
PyDSTool. There is an AD package for python, which I have not used.
The applications I had in mind would probably require RHS compilation,
along the lines of the Jorba & Zhou package, in order for it to run in
reasonable time. At this point, AD is fairly low priority, though it
would certainly give you (Anne) all of the information you are looking
for along the trajectory. Taylor series integrators do not avoid
stiffness problems, of course.
Please keep up the discussion, as it is very useful for us to get lots
of input as we make design decisions. It would be nice to hear from
people doing dynamical systems work, as well as those who use ODEs
from a different perspective. The more examples we have of how people
want to use the integrators, the better we can plan interfaces and
William Erik Sherwood, Ph. D.
Center for BioDynamics, Boston University
111 Cummington Street
Boston, MA 02215
More information about the SciPy-user