[SciPy-user] Solving Two boundary value problems with Python/Scipy

[Robert Clewley]

> 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