[SciPy-Dev] ode and odeint callable parameters

[Benny Malengier]

2016-01-25 14:29 GMT+01:00 Marcel Oliver <m.oliver at jacobs-university.de>:

> Anne Archibald writes:
>  > <snip>
>  > to implement it myself. And anyway, some people need very different
>  > things from their ODE solvers (for example, solution objects
>  > evaluable anywhere).
>  >
>  > Anne
>
> <snip>
>
> 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...)
>

Yes rootfn does that.
A nicer interface with onroot functionality to that was added by a
contributor to CVODE, I just extended it to the IDA solver. Only available
through the solve method of the solvers (as STEP methods just stop),
example:
https://github.com/bmcage/odes/blob/master/docs/src/examples/ode/freefall.py
Same with fixed tstop added in a list of output times:
https://github.com/bmcage/odes/blob/master/docs/src/examples/ode/freefall_tstop.py


>
> 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.
>

agree. Also, just reinit the solver and solve again from closest start time
is an option too if for some reason interpolation error must be avoided.
CVODE can actually give you output backward if you have jumped too far
ahead, but that requires you to know where you want interpolation output
immediately after doing a STEP. All very problem depending.

>
> 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.
>

These things all require a lot of boilerplate code. In the end, my
experience is that the documentation of the original solvers is the best to
learn the details, so keeping variable names the same in options is more
important than something that is more pythonic. Eg, the rootfn above is a
bad name, but fits the original documentation at
https://computation.llnl.gov/casc/sundials/documentation/documentation.html
CVODE doc is 164 pages, in the end, for problems you want to optimize you
will grab that documentation to learn what options you want to tweak.
About transparant algebraic part. The backends need a lot of info on the
algebraic part. You also need to structure your variables specifically if
you eg want banded Jacobian. And the calling sequence of all functions
becomes different with an extra parameter in and out (ydot).
Quite a lot of boilerplate that slows your wrapper down will be needed I'm
afraid.


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

I'm not really following. You mean unwrap and wrap again like the
complex_ode solver does now in scipy?
http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.complex_ode.html

Writing your own wrapper for your problem does not seem like a big task to
me.

>
> 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?)
>

They are used for that too, but it is not the core. Large complex kinetic
systems with different time scales seems like the largest use case on the
mailing list, see also as first example in the example doc.
Sensitivities are not present in odes at the moment, but I think some
wrappers have that covered already. Let's say it is asked once or twice a
year to add.


> 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.)
>

Some extra tests seem in order. How does method=adams in ode succeed?
Should you be able to quickly make an ipython notebook ...


> 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.
>

Yes. We don't have our own backends, we only wrap backends. In the end,
those backends with continued development will win out. Wrapping solvers is
a tricky business, as it is very hard to avoid design of the solver to not
slip through


>
> Regards,
> Marcel
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> https://mail.scipy.org/mailman/listinfo/scipy-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20160125/e34dbc3b/attachment.html>