> It seems to me it would be useful for scipy to have an "odetools"
> package, providing convenience features on top of our existing ODE
> solvers (as well as a better fortran ODE solver, namely lsodar, but
> that's beside the point).
Sure, and I've said before that the SciPy powers-that-be are welcome to
use any parts of our PyDSTool package as the basis of such a package.
e.g., our SWIG-based interfaces to two great integrators (Dopri and
Radau). Our interface includes events at the C-level and external inputs
interpolated from arrays. Our ODEGenerator classes provide a lot of
convenience features on top of the solvers (including the existing Scipy
VODE solver).
>> Features could include:
> * Implementation of the shooting method
Yes, that would be useful.
> * class ODESolution, which provides an interpolated function object
That's exactly what our Variable and Trajectory classes do, built over
scipy's interp1d.
> * An implementation (necessarily slow and limited) of the finite
> difference method
Our diff function provides that with various options for
functions from R^N -> R^M. Here's the signature and docstring
from PyDSTool/common.py:
def diff(func, x0, vars=None, axes=None, dir=1, dx=None):
"""Numerical 1st derivative of R^N -> R^M function about x0 by finite
differences.
dir=1 uses finite forward difference.
dir=-1 uses finite backward difference.
List valued dx rescales finite differencing in each axis separately.
vars argument specifies which elements of x0 are to be treated as
variables for the purposes of taking the Jacobian.
If axes argument is unused or set to be all axes, the Jacobian of the
function evaluated at x0 with respect to the variables is returned,
otherwise a sub-matrix of it is returned.
If dx is not given an appropriate step size is chosen (proportional to
sqrt(machine precision)).
"""
> * Tools for working with linear ODEs (generate a complete set of solutions, say)
> * Tools for solving eigenvalue problems
>> Any other suggestions?
A set of phase plane tools, like finding nullclines, fixed points,
stability.... I have a whole bunch of code working for this
that will end up in the next release of PyDSTool (for Scipy 0.5.1,
finally) and I'll be adding a few additional features like calculating
stability of fixed points.
The port to Numpy and new Scipy is will be done in the next few weeks so
code will be available for people to clean up to their liking for use in
SciPy.
-Rob