[SciPy-Dev] Improvements (scipy.integrate._ode) -- accessing the dense output of DOPRI5 and DOP853.

[Jim Martin]

** Summary **

The wrappers around the differential equation solvers DOPRI5 and DOP853 in
"scipy/integrate/_ode.py" do not currently allow access to their "dense
output" option (iout=2).  This option provides the coefficients of an
interpolating polynomial for the solution between the last and current
integration step.  Dense output interpolation is commonly used for event
location and improved plotting.  The "_ode.py" code can be modified to
optionally allow access to the dense output of DOPRI5/DOP853.  The proposed
changes are here:

https://github.com/scipy/scipy/compare/master...jddmartin:dense_output_from_dopri5_and_dop853
with example usage here:
  https://github.com/jddmartin/dense_output_example_usage
Existing code using "_ode.py" would not be impacted by the changes.

** Background **

The efficient numerical solution of differential equations requires
adaptive stepsizes.  Modern high order Runge-Kutta code --- such as DOPRI5
and DOP853 --- can take large variable stepsizes, so that plots made simply
by connecting lines between steps are of poor visual quality.  As discussed
by [Shampine, 1997, see Section 5 and Figure 1, and Shampine et al., 2003,
pages 54, 55] this issue is addressed in the Matlab solvers by
interpolation between integration steps, with a default of 4 output
solutions provided per step (but can be modified using the "Refine"
option). Interpolation of the solution output between steps is also
desirable for event location (see discussion in Shampine et al., 2003,
Section 2.3.1).

The Fortran codes DOPRI5/DOP853 have a dense output option (iout=2), as
discussed in Hairer et al.'s book [Hairer, 1993].  This dense output option
is illustrated in the example usage code for these routines:  dr_dopri5.f
and dr_dop853.f (see http://www.unige.ch/~hairer/software.html ). However,
without implementation of the iout=2 option, these examples cannot be
replicated using scipy.  Some popular sources such as Press et al.'s
"Numerical Recipes" and Boost's odeint have solvers based on DOPRI5/DOP853
and provide access to dense output.

** Details **

The original DOPRI5/DOP853 codes provide the option for a user-specified
callback function to be called after each internal integration step.
Currently this callback function is only passed the solution at the current
step.  With the proposed change, the callback function will be provided
(optionally) with the coefficients of an interpolating polynomial for the
solution between the last and current step by using the "dense output"
option (iout=2).

As far as the user of scipy.integrate is concerned, the main change is the
addition of an optional "dense_components" keyword argument to the
".set_solout" method of classes "ode" and "complex_ode".  Behavior should
be unchanged in the absence of this keyword.  In addition, I have put in
some additional error handling code for the bug reported here:
  https://github.com/scipy/scipy/issues/4118
and made some minor changes (see git commit messages) that should not
impact existing code.

The new "dense_components" keyword argument is used to specify the
components for which dense output is desired.  If this keyword option is
present, the callback function will now be passed the additional
information required for interpolation.  A convenience function "dense_dop"
has been provided in _ode.py for interpolation using this information.

The introduction of the "dense_dop" function seems a bit awkward ---
suggestions for alternative approaches are welcome.  Most users will simply
want to evaluate the interpolation function to determine the solution at
points between the last and current step position.   My motivation for
development of access to the dense output is an application in which I will
use the coefficients of the interpolating polynomial directly, rather than
the interpolated solution.  Access to the interpolating polynomial may also
be useful for event location.

I've updated the relevant docstrings and added test code for this new
option (I have only tested this code on my Linux system with Python 2.7).

Comments are welcome.  Thanks -- Jim.

References:
[Shampine, 1997]: Lawrence F. Shampine and Mark W. Reichelt, The MATLAB ODE
Suite, SIAM J. Sci. Comput., 18(1), 1–22. 1997 (22 pages)
http://dx.doi.org/10.1137/S1064827594276424

[Shampine, 2003], L. F. Shampine, I. Gladwell, S. Thompson, "Solving ODEs
with Matlab", 2003, Cambridge University Press.

[Hairer, 1993] Hairer et al., "Solving Ordinary Differential Equations,
Nonstiff Problems",  Second Revised Edition, Springer, 1993.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20150418/43f49294/attachment.html>