[SciPy-user] PyDSTool and SciPy (integrators)

[Anne Archibald]

On 18/04/2008, Gabriel Gellner <ggellner at uoguelph.ca> wrote:
> On Fri, Apr 18, 2008 at 01:06:01PM -0400, Anne Archibald wrote:

>  > 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?

I think that the correct tool for this is a spline whose values and
(first?) derivatives are specified at each control point; higher
derivatives would be determined by smoothness and minimum bending
energy constraints as usual. I don't think scipy has tools for
constructing splines like this, but I would be willing to take a swing
at writing them using scipy's banded matrix solver. My idea would be
to produce a data structure that could be evaluated using scipy.splev.
This gives, for free, C implementations of evaluation, integration,
and derivatives.

>  > 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 . . .)

My complaints with it, apart from the mass of impenetrably-named
keyword arguments and "infodict" keys, have to do with control of the
t values. Since you can't stop at a place based on the y value, I
usually end up having to use it in a "running" mode, which it supports
rather badly. If I want to evaluate the solution at multiple points,
it is usually because I want an approximation to the true solution
function, but such an approximation should normally have points
distributed in a way that takes into account the behaviour of the
function - more where it behaves in a complicated fashion, and fewer
where it is simple.

On the other hand, for the lightweight uses I've had, it's quite
inconvenient to have to prepend the initial time to the list (often
only one element long) of times at which I want to evaluate the
result.

Anne