# dopri853: Dormand-Prince Runge-Kutta of Order 8(5,3) In deTestSet: Test Set for Differential Equations

## Description

Solves the initial value problem for systems of ordinary differential equations (ODE) in the form:

dy/dt = f(t,y)

The R function `dopri853` provides an interface to the Fortran ODE solver DOP853, written by Hairer and Wanner.

It implements the explicit Runge-Kutta method of order 8(5,3) due to Dormand & Prince with stepsize contral and dense output

The system of ODE's is written as an R function or can be defined in compiled code that has been dynamically loaded.

## Usage

 ```1 2 3 4 5``` ```dopri853 (y, times, func, parms, rtol = 1e-6, atol = 1e-6, verbose = FALSE, hmax = NULL, hini = hmax, ynames = TRUE, maxsteps = 10000, dllname = NULL, initfunc = dllname, initpar = parms, rpar = NULL, ipar = NULL, nout = 0, outnames = NULL, forcings = NULL, initforc = NULL, fcontrol = NULL, ...) ```

## Arguments

 `y ` the initial (state) values for the ODE system. If `y` has a name attribute, the names will be used to label the output matrix. `times ` time sequence for which output is wanted; the first value of `times` must be the initial time; if only one step is to be taken; set `times` = `NULL`. `func ` either an R-function that computes the values of the derivatives in the ODE system (the model definition) at time t, or a character string giving the name of a compiled function in a dynamically loaded shared library. If `func` is an R-function, it must be defined as: `func <- function(t, y, parms,...)`. `t` is the current time point in the integration, `y` is the current estimate of the variables in the ODE system. If the initial values `y` has a `names` attribute, the names will be available inside `func`. `parms` is a vector or list of parameters; ... (optional) are any other arguments passed to the function. The return value of `func` should be a list, whose first element is a vector containing the derivatives of `y` with respect to `time`, and whose next elements are global values that are required at each point in `times`. The derivatives should be specified in the same order as the state variables `y`. If `func` is a string, then `dllname` must give the name of the shared library (without extension) which must be loaded before `lsode()` is called. See package vignette `"compiledCode"` for more details. `parms ` vector or list of parameters used in `func` or `jacfunc`. `rtol ` relative error tolerance, either a scalar or an array as long as `y`. See details. `atol ` absolute error tolerance, either a scalar or an array as long as `y`. See details. `verbose ` if `TRUE`: full output to the screen, e.g. will print the `diagnostiscs` of the integration - if the method becomes stiff it will rpint a message. `hmax ` an optional maximum value of the integration stepsize. If not specified, `hmax` is set to the largest difference in `times`. `hini ` initial step size to be attempted. `ynames ` logical, if `FALSE` names of state variables are not passed to function `func`; this may speed up the simulation especially for multi-D models. `maxsteps ` maximal number of steps taken by the solver, for the entire integration. This is different from the settings of this argument in the solvers from package deSolve! `dllname ` a string giving the name of the shared library (without extension) that contains all the compiled function or subroutine definitions refered to in `func` and `jacfunc`. See vignette `"compiledCode"` from package `deSolve`. `initfunc ` if not `NULL`, the name of the initialisation function (which initialises values of parameters), as provided in ‘dllname’. See vignette `"compiledCode"` from package `deSolve`. `initpar ` only when ‘dllname’ is specified and an initialisation function `initfunc` is in the dll: the parameters passed to the initialiser, to initialise the common blocks (FORTRAN) or global variables (C, C++). `rpar ` only when ‘dllname’ is specified: a vector with double precision values passed to the dll-functions whose names are specified by `func` and `jacfunc`. `ipar ` only when ‘dllname’ is specified: a vector with integer values passed to the dll-functions whose names are specified by `func` and `jacfunc`. `nout ` only used if `dllname` is specified and the model is defined in compiled code: the number of output variables calculated in the compiled function `func`, present in the shared library. Note: it is not automatically checked whether this is indeed the number of output variables calculed in the dll - you have to perform this check in the code - See vignette `"compiledCode"` from package `deSolve`. `outnames ` only used if ‘dllname’ is specified and `nout` > 0: the names of output variables calculated in the compiled function `func`, present in the shared library. These names will be used to label the output matrix. `forcings ` only used if ‘dllname’ is specified: a list with the forcing function data sets, each present as a two-columned matrix, with (time,value); interpolation outside the interval [min(`times`), max(`times`)] is done by taking the value at the closest data extreme. See forcings or package vignette `"compiledCode"`. `initforc ` if not `NULL`, the name of the forcing function initialisation function, as provided in ‘dllname’. It MUST be present if `forcings` has been given a value. See forcings or package vignette `"compiledCode"`. `fcontrol ` A list of control parameters for the forcing functions. See forcings or vignette `compiledCode`. `... ` additional arguments passed to `func` and `jacfunc` allowing this to be a generic function.

## Details

The work is done by the FORTRAN subroutine `dop853`, whose documentation should be consulted for details. The implementation is based on the Fortran 77 version fromOctober 11, 2009.

The input parameters `rtol`, and `atol` determine the error control performed by the solver, which roughly keeps the local error of y(i) below rtol(i)*abs(y(i))+atol(i).

The diagnostics of the integration can be printed to screen by calling `diagnostics`. If `verbose` = `TRUE`, the diagnostics will written to the screen at the end of the integration.

See vignette("deSolve") from the `deSolve` package for an explanation of each element in the vectors containing the diagnostic properties and how to directly access them.

Models may be defined in compiled C or FORTRAN code, as well as in an R-function. See package vignette `"compiledCode"` from package `deSolve` for details.

Information about linking forcing functions to compiled code is in forcings (from package `deSolve`).

## Value

A matrix of class `deSolve` with up to as many rows as elements in `times` and as many columns as elements in `y` plus the number of "global" values returned in the next elements of the return from `func`, plus and additional column for the time value. There will be a row for each element in `times` unless the FORTRAN routine ‘lsoda’ returns with an unrecoverable error. If `y` has a names attribute, it will be used to label the columns of the output value.

## Author(s)

Karline Soetaert <[email protected]>

## References

E. Hairer, S.P. Norsett AND G. Wanner, Solving Ordinary Differential Equations I. Nonstiff Problems. 2nd Edition. Springer Series In Computational Mathematics, SPRINGER-VERLAG (1993)

• `ode` for a general interface to most of the ODE solvers from package `deSolve`,

• `ode.1D` for integrating 1-D models,

• `ode.2D` for integrating 2-D models,

• `ode.3D` for integrating 3-D models,

• `mebdfi` for integrating DAE models,

• `gamd` for the generalised adams method

`diagnostics` to print diagnostic messages.

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45``` ```## ======================================================================= ## Example : ## The Arenstorff orbit model ## ======================================================================= Arenstorff <- function(t, y, parms) { D1 <- ((y[1]+mu)^2+y[2]^2)^(3/2) D2 <- ((y[1]-(1-mu))^2+y[2]^2)^(3/2) dy1 <- y[3] dy2 <- y[4] dy3 <- y[1] + 2*y[4]-(1-mu)*(y[1]+mu)/D1 -mu*(y[1]-(1-mu))/D2 dy4 <- y[2] - 2*y[3]-(1-mu)*y[2]/D1 - mu*y[2]/D2 list(c(dy1, dy2, dy3, dy4)) } #----------------------------- # parameters, initial values and times #----------------------------- mu <- 0.012277471 yini <- c(x = 0.994, y = 0, dx = 0, dy = -2.00158510637908252240537862224) times <- seq(0, 18, 0.01) #----------------------------- # solve the model #----------------------------- out <- dopri853 (times = times, y = yini, func = Arenstorff, parms = NULL, rtol = 1e-17, atol = 1e-17) plot(out[,c("x", "y")], type = "l", lwd = 2, main = "Arenstorff") #----------------------------- # First and last value should be the same #----------------------------- times <- c(0, 17.0652165601579625588917206249) Test <- dopri853 (times = times, y = yini, func = Arenstorff, parms = NULL) diagnostics(Test) ```

deTestSet documentation built on May 29, 2017, 2:56 p.m.