[SciPy-Dev] ode and odeint callable parameters

[Marcel Oliver]

Anne Archibald writes:
 > There is periodically discussion of the mess of interfaces and
 > solvers for ODEs in scipy (see the archives about six months ago,
 > for example). One of the concerns is that people want to do very
 > different things, so settling on a good interface is not at all
 > easy. I don't just mean the underlying algorithms, but also the
 > calling interface. It's easy to support someone who drops in an RHS
 > and asks for the solution at a hundred predefined points. It's also
 > clear that a general toolkit shouldn't support my use case: a
 > 22-dimensional solver with a compiled RHS generated by a symbolic
 > algebra package, with internal precision switchable between 80 and
 > 128 bits, with a root-finder attached to the output to define
 > stopping places, that needs to all run with compiled speed. So
 > where along that scale do you aim scipy's interface? I have my own
 > ideas (see aforementioned archive thread) but don't have the energy
 > to implement it myself. And anyway, some people need very different
 > things from their ODE solvers (for example, solution objects
 > evaluable anywhere).
 > 
 > Anne

Thanks for all the comments.  I had not been aware of the existence of
scikits odes - the additional capabilities are very nice, but API-wise
it just adds to the mess...

Being unsure what has previously been discussed, I'd just like to
share some random thoughts, as I don't claim to understand the problem
domain completely.  Basically, the API should make simple thing simple
and natural, and complicated things possible.  So I think your
requirements stated above should not be completely out of reach - in
principle - even though high precision arithmetic would require
completely new backends and whether one can get close to compiled
speed is probably very problem dependent.

So, here an unordered list of thoughts:

1. I don't think there is a big difference in practice between the
   ode, odes, and odeint interfaces.  Personally, think that

     odeint (f, y0, t)

   is just fine.  But the object oriented approach of ode/odes is
   good, too; however, the various parameters should be controlled by
   (keyword) arguments as a statement like

     r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)

   is just not pretty.

   So I like the odes approach best, with a basic object oriented
   structure, but the parameters set as arguments, not via methods.

   Not sure if odes allows to set an initial time different from zero,
   that is essential.  I also dislike that specifying the solver is
   mandatory and first in the argument list; I think this is
   specialist territory as many of the existing solvers perform
   sufficiently well over a large range of problems so that changing
   the solver is more of an exception.

   If the interface is object oriented (as in ode or odes), then at
   least the default solver should be reentrant, otherwise the
   API raises false expectations.

   I would also suggest that, independent of the rest of the API, the
   final time t can be an array or a scalar; the result is either
   returned as a vector of shape y0 when t is scalar and as an array
   with one additional axis for time when t is a 1-d-array.

2. It would be very nice to control output time intervals or exit
   times.  Mathematica has a fairly rich "WhenEvent" API, although it
   seems to me it may also be no more than a fancy wrapper around a
   standard solver.

   The Mathematica API feels almost a bit too general, but there are
   certainly two independent ways triggers could be implemented, both
   seem useful to me:

   (a) Simple Boolean trigger for "output event" and "termination
       event" based on the evaluation of some function of (y,t,p)
       which is passed to the solver and which is evaluated after each
       time step.

   (b) Trigger for "output event" and/or "termination" based on
       finding the root of a scalar function of the (y,t,p).  This
       needs either a separate root finder or a modification to the
       implicit solver (the latter may be hard as that would require
       nontrivial modification of the embedded solver code).

       (odes "rootfn" may be doing just that, but there seems to be no
       detailed documentation...)

3. Returning interpolating function objects as, e.g., Mathematica
   does, seems a bit of overkill to me.  It would be easier to set up
   an independent general interpolating function library (as possibly
   even exists already) and pass a discrete time series to it.  It is
   likely that by deriving the interpolation function directly from
   the interpolation on which the solver is based one could gain some
   efficiency and give better control of the interpolation accuracy,
   but to me that seems to require a lot of effort for relatively
   little improvement.

4. odes includes DAE solvers.  I wonder if it would not be more
   natural to pass the algebraic part of the equation as a separate
   keyword argument and let the solver wrapper decide, based on the
   presence or absence of this argument, whether to call a DAE
   backend.

5. As I wrote before, the API should deal transparently with
   d-dimensional vector fields.

Just a question about solver backends: it seems that the SUNDIALS
solvers are particularly strong for large stiff systems which come
from parabolic PDEs?  (And of course the ability to propagate
sensitivities, is that ability made available through odes, too?)
On smaller system of ODEs, I have found that ode "dop853" is extremely
accurate and seems capable of dealing with rather stiff problems even
though I have no clue how it does it.  (E.g., for the Van der Pol
oscillator with large stiffness parameter where other explicit solvers
either blow up or grind to a halt, dop853 and also dopri5 are doing
just fine while Matlab's dopri5 solver fails as expected.  It's a
mystery to me.)

In any case, it would be very nice to develop a clear plan for
cleaning up ODE solver backend in Scipy.  Once there is a clear
target, it might be easier to see what is easy to do and what is
missing on the backend side.

Regards,
Marcel