DOPRI5 is an explicit one-step integration code that we are currently working on connecting with ASCEND. It was originally written in Fortran by Ernst Hairer & Gerhard Wanner at the Université de Genève, later translated to C by J Colinge. The C code version describes itself as follows:
This code computes the numerical solution of a system of first order ordinary
differential equations y'=f(x,y). It uses an explicit Runge-Kutta method oforder (4)5 due to Dormand & Prince with step size control and dense output.
We are implementing a wrapper for the essentially unmodified DOPRI5 code as an additional choice for integration algorithms for a couple of reasons:
- a one-step explicit integrator is the simplest kind and in many ways should be the starting point for solving ODE problems; if difficulties are detected, then using a one-step explicit integrator might help to indicate the nature of those difficulties.
- DOPRI5 can detect likely stiff problems and alert the user in that case
- it is stable and well-tested
- it provides predictor-corrector adjustment of step length
- it allows us to test the design of our Integrator API and identify shortcomings in it
- the developer has a problem that he is concerned may not be as stiff as he thought it was; this will aid in understanding whether or not that is true.
As part of this work, the ability to load external Integrators at runtime will be provided. At present, this facility has not yet been implemented for Solvers, but that solver code is older and may require more re-working before that is possible. So this will provide some useful experience before that is attempted.
Changes made to the DOPRI5 code were simply those required to allow a 'user pointer' to be passed through DOPRI5 to the FcnEqnDiff and SolTrait hook functions, allowing DOPRI5 to interface with ASCEND without the need for global variables. These changes can be browsed online at solvers/dopri5.
Implementation of DOPRI5 is basically complete, but this integrator needs careful testing. Also there is some discrepancy between our integration results and the results of test problem that hasn't yet been resolved (see plot below). DOPRI5 will be available in standard ASCEND distributions but will need to be IMPORTed before it wil be accessible. In other words, if writing an ODE model, add the following at the start:
The option 'DOPRI5' should then appear in the Integrator dialog, at least in the PyGTK GUI.
The sample Arenstorf Orbit problem is available in the Model Library as models/test/dopri5/aren.a4c. The expected solution of this model is that at x=xend, the values of y and y are equal to their starting values. The FORTRAN implementation gives this result accurate to 1e-5 absolute.
File:Diverging.png Arenstorf orbit integrated with the experimental DOPRI5 integrator. We're still working on why this solution is diverging as currently shown.
Parameters for the DOPRI5 integrator will be documented here as they are implemented. The full set of possible parameters can be seen in the source code for solvers/dopri5/dopri5.h
- tolvect: whether to use a single tolerance for all variables (0), or whether to use a per-variable tolerances (1). Default is 1.
- rtol: relative tolerance for all variables if tolvect is 0. Otherwise ignored.
- atol: absolute tolerance for all variables if tolvect is 0. Otherwise ignored.
- nstiff: number of steps considered stiff. If too many (>nstiff) small steps are being required by the predictor-corrector pair, integration will be aborted with a warning about a likely stiff problem. Default is 1000.
The initial and maximum step size can also be configured using the standard Integrator API functions.
Parameters can be set using the Integrator window, or using Python scripts.
Need to fix asc_dopri5.c to make use of the contd5 function provided by DOPRI5 for polynomial interpolation within timesteps. Currently that code is disabled.