[SciPy-user] PyDSTool and SciPy (integrators)

[Anne Archibald]

On 21/04/2008, Rob Clewley <rob.clewley at gmail.com> wrote:

>  I am very interested in getting spline support to a point where it
>  could be used as a more accurate way of representing solutions to
>  ODEs. I will certainly use it in my Trajectory class if you give it a
>  similar interface to that already used by interp1d.

The interface I've used is similar to that used in interp1d, though a
bit more complicated as I rely on having derivative information (if
you want a continuous derivative). You can feed my code a list/array
of x values and a list the same length of y values and derivatives,
but the list of y values is of the form
[[y1, y1prime, y1doubleprime],
 [y2, y2prime],
 ...
]

You can construct a PiecewisePolynomial by giving it x and y values
(in the above format), or by appending new points onto the end of an
existing PiecewisePolynomial.

> In case this can help in your problem, our current interface to the
>  H&W codes allows you to specify time points at which the integrator is
>  guaranteed to step on to as it goes along.
>
>  If someone is willing to help improve the low level interfaces to the
>  H&W codes by exposing internal information like the polynomial
>  coefficients for the dense solution then I would include that and
>  other trajectory-related data in our ODE support.

When working out this interface design I looked at the FORTRAN
interface to LSODA/LSODAR. The way I envisioned the interface working
was that the inner loop would ask LSODAR to take a single internal
step then return; the inner loop would then read off the time that had
been reached, the order, and all the available derivatives from arrays
the FORTRAN code exposed. This information then allows you to add a
(single) new point to the end of the trajectory; the polynomial is
constructed to have the correct order and to match enough derivatives
at each end to completely determine it. (Thus it will not be *exactly*
the polynomial the solver uses, but should be close. The solver's
internal polynomial is probably not even continuous from step to
step.) Other information related to the solution and its convergence
can straightforwardly be added.

At the moment all this will be relatively slow, since the
PiecewisePolynomial object is implemented in pure python. But for
modest orders - 12 is the highest LSODAR uses - the operations are
fairly inexpensive, so conversion to cython should produce a
Trajectory object that is about as good an approximation to the
solution as the solver's internal polynomials, with exactly the same
knots as the solution steps.

I think it makes a certain amount of sense to use the
one-step-at-a-time approach if you want to keep convergence
information, since you need to store information per step in any case.

Time points that require special care are a valuable feature, I think
(certainly the FORTRAN code goes out of its way to support them), but
didn't solve my original problem, as I only knew where problems would
occur based on the y coordinates rather than the t coordinates.

Anne