[SciPy-user] PyDSTool and SciPy (integrators)

[Gabriel Gellner]

On Fri, Apr 18, 2008 at 01:06:01PM -0400, Anne Archibald wrote:
> On 18/04/2008, Erik Sherwood <wesher at bu.edu> wrote:
> 
> >  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.
> 
> Just to check: radau and dopri853 require you to know whether your
> problem is stiff or not? An extremely convenient feature of the
> ODEPACK integrator is that it can switch between stiff and non-stiff
> solvers as needed.
> 
True, using the radau and dopri codes is much more like using the matlab codes,
you must explicitly choose the method.

> >  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.
> 
> Most of the extra information I had in mind was intended to reduce
> memory usage and improve smoothness of the (interpolated) trajectory.
> Since the trajectory is the result of solving an ODE, it will normally
> be at least differentiable, but linear interpolation will produce a
> non-differentiable trajectory. I can envision many situations where
> the derivative of the objective function is wanted. This can at worst
> be obtained by evaluating the RHS, but if the integrator is internally
> working with something that is locally a polynomial, ideally the
> trajectory would follow that. Using polynomials, at least to take the
> derivative into account, should reduce the density with which one must
> sample the trajecotry.
> 
Does scipy have the required interpolation routines? Could we do just as you
say, save the information and then have the option to interpolate as needed?

> The other information, on convergence, I suggested mostly so that
> users could evaluate the trajectory and look for places where the
> solution risked being inaccurate. The current ODEINT provides a fair
> amount of such instrumentation, in its messy FORTRAN way.
> 
I think outputting this information is necessary, if the size is an issue we
can just have a flag that turns it off, but in general I think having
convergence codes available to be essential.

> The "running" mode I suggested is based on the way ODEPACK is designed
> to be called from FORTRAN: one call initializes the integrator, then
> it maintains internal state. It can be told to advance either one step
> (of its choice) or a fixed amount (obtained by interpolation based on
> the solution). The result can then be "walked" forward, keeping the
> integrator's internal state (derivatives, appropriate step size,
> knowledge about stiffness) at each step. This kind of calling
> convention is useful for (say) particle systems, where for every frame
> of an animation one wants to advance the system by a thirtieth of a
> second, continuing indefinitely and forgetting the past.
> 
The problem with this of course is that it would be really slow in python. But
I don't see any reason why we can't do this if the user wants. I imagine this
might be a nice use of a generator.

> >  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.
> 
> The odeint interface is, frankly, a disaster. I have had to abuse it
> every time I needed to solve an ODE. The underlying ODEPACK interface
> is fairly horrible, but I think it's valuable to look at it because
> it's seen a lot of use and has been adjusted to be able to do all the
> various things that users want from it.
> 
Could you give some examples why it is a disaster? As I would like avoid
repeating said mistakes. (Above and beyond the example of stepping into
regions that might be undefined . . .)

Gabriel