### Improvements (scipy.integrate._ode) -- accessing the dense output of DOPRI5 and DOP853.

2015-04-18 17:43:09 GMT

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.

_______________________________________________ SciPy-Dev mailing list SciPy-Dev <at> scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev