[SciPy-user] PyDSTool and SciPy (integrators)

Erik Sherwood wesher@bu....
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  
so far:

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.

4. Compilation

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.

5. Misc

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  

Thanks again,


William Erik Sherwood, Ph. D.
Research Fellow
Center for BioDynamics, Boston University
111 Cummington Street
Boston, MA 02215
(617) 358-4684

More information about the SciPy-user mailing list