[SciPy-user] PyDSTool and SciPy (integrators)

Anne Archibald peridot.faceted@gmail....
Thu Apr 17 19:58:05 CDT 2008

On 17/04/2008, Gabriel Gellner <ggellner@uoguelph.ca> wrote:

>  > What I have needed is an integrator that offers more control. The
>  > particular problem I ran into was one where I needed to be able to set
>  > stopping conditions based on the dependent variables, not least
>  > because my model stopped making sense at a certain point. I ended up
>  > wrapping LSODAR from ODEPACK, and it made my life vastly easier. If
>  > you have integrators that can provide that kind of control, that would
>  > be very useful. Similarly something like a trajectory object that
>  > knows enough about the integrator to provide a spline with the right
>  > derivatives and knot spacing would be handy.
> The dopri5 and dop853 codes have a really good interface for this, that is
>  after each successful integration step they call a solution callback that
>  allows for dense interpolation of the solution between stepping points (with
>  error of the same order as the solution), and any kind of stopping or event
>  detection routines. The hard part is thinking of good ways of making default
>  callbacks that are extensible, but it is nice that the base codes were
>  designed with this functionality.

This sounds like it might have been a problem for me: what was
happening was I was integrating an ODE within a domain. Outside the
domain my RHS function raised exceptions, so it was not even safe to
evaluate it there. I tried hacking up a solution with scipy's
one-step-at-a-time integrator modes, but the best I could do was stop
"near" the boundary, which failed if the integrator took bigger steps
than my notion of "near".

>  Ultimately with my Cython wrapper you could just pass in an optional callback
>  as well (coded in C or Cython as otherwise it would be painfully slow), that
>  would give you any kind of control you might want. Though the way forward is
>  to have good default templates available, I have been looking at PyDSTool's
>  method as well as how matlab (not the code, rather the api) deals with this
>  control.
>  If you now of any other good examples of event detection or stopping rule
>  api's I would love to look at them!

I used LSODAR in ODEPACK and it did just what I wanted: before
evaluating the RHS at any point, it called a list of user-provided
functions. If any of them returned negative values, that point was
outside the domain and the integrator took a smaller step. It assumed
the user-provided stopping functions were somewhat smooth; in effect
it ended up doing some root-finding so that it could stop right at the
boundary. The API I came up with was simple and ugly: I just passed a
list of functions, each of which took the same kind of pair of values
as the RHS does, and each of which returned a float of appropriate

> What do people think about returning an object from the integrators? I have
>  been using R a lot lately and I really like how they use simple data objects
>  to represent the output of solutions. That way you could have an object that
>  would work like an array of the solution points, but would also have the full
>  output available as attributes (like the time vector, even locations, etc).
>  That way, just like in R, the complexity of the solution would be hidden, but
>  it gets rid of the awkward use of 'full_output' booleans, and dealing with the
>  resulting tuples (which seems like a suboptimal attempt to get matlab like
>  functionality with out the syntax that matlab gives for only taking the
>  information you need).

I think the right return format depends what people want from the
solution. I see two major kinds of problem one might use an ODE solver

* Integrate out to some stopping point and give the values of
dependent and independent variable there.
* Integrate out to some stopping point and give back a representation
of the solution.

(The first case can in principle be obtained by evaluating the second
at the endpoint, but at the cost of replacing O(1) space by O(n)

For the second it is clear that some kind of object to represent the
solution is a good idea. What exactly it should contain is open for
discussion; I envision keeping estimated function values and
derivatives at each step, along with convergence information, and
providing an interface to view it as either an array of values or a
function akin to our interpolating spline objects. Convergence
information and auxiliary function values may also be useful. Even for
the first, simple case, I think returning an object is a more
manageable way to provide convergence information.

Are there other styles of using an ODE solver that people want?
Perhaps a "running" mode, where the solver keeps track of whatever
state information it has and is asked to advance by a certain amount
(say, because more experimental data has become available)?


More information about the SciPy-user mailing list