rk4 {deSolve} R Documentation

## Solve System of ODE (Ordinary Differential Equation)s by Euler's Method or Classical Runge-Kutta 4th Order Integration.

### Description

Solving initial value problems for systems of first-order ordinary differential equations (ODEs) using Euler's method or the classical Runge-Kutta 4th order integration.

### Usage

```euler(y, times, func, parms, verbose = FALSE, ynames = TRUE,
dllname = NULL, initfunc = dllname, initpar = parms,
rpar = NULL, ipar = NULL, nout = 0, outnames = NULL,
forcings = NULL, initforc = NULL, fcontrol = NULL, ...)

rk4(y, times, func, parms, verbose = FALSE, ynames = TRUE,
dllname = NULL, initfunc = dllname, initpar = parms,
rpar = NULL, ipar = NULL, nout = 0, outnames = NULL,
forcings = NULL, initforc = NULL, fcontrol = NULL, ...)

euler.1D(y, times, func, parms, nspec = NULL, dimens = NULL,
names = NULL, verbose = FALSE, ynames = TRUE,
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 ` times at which explicit estimates for `y` are desired. The first value in `times` must be the initial time. `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 must 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 `rk4` is called. See package vignette `"compiledCode"` for more details. `parms ` vector or list of parameters used in `func`. `nspec ` for 1D models only: the number of species (components) in the model. If `NULL`, then `dimens` should be specified. `dimens` for 1D models only: the number of boxes in the model. If `NULL`, then `nspec` should be specified. `names ` for 1D models only: the names of the components; used for plotting. `verbose ` a logical value that, when `TRUE`, triggers more verbose output from the ODE solver. `ynames ` if `FALSE`: names of state variables are not passed to function `func` ; this may speed up the simulation especially for large models. `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`. See package vignette `"compiledCode"`. `initfunc ` if not `NULL`, the name of the initialisation function (which initialises values of parameters), as provided in ‘dllname’. See package vignette `"compiledCode"`, `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 calculated in the DLL - you have to perform this check in the code. See package vignette `"compiledCode"`. `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. `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` allowing this to be a generic function.

### Details

`rk4` and `euler` are special versions of the two fixed step solvers with less overhead and less functionality (e.g. no interpolation and no events) compared to the generic Runge-Kutta codes called by `ode` resp. `rk`.

If you need different internal and external time steps or want to use events, please use: `rk(y, times, func, parms, method = "rk4")` or `rk(y, times, func, parms, method = "euler")`.

See help pages of `rk` and `rkMethod` for details.

Function `euler.1D` essentially calls function `euler` but contains additional code to support plotting of 1D models, see `ode.1D` and `plot.1D` for details.

### 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 integration routine returns with an unrecoverable error. If `y` has a names attribute, it will be used to label the columns of the output value.

### Note

For most practical cases, solvers with flexible timestep (e.g. `rk(method = "ode45")` and especially solvers of the Livermore family (ODEPACK, e.g. `lsoda`) are superior.

### Author(s)

Thomas Petzoldt thomas.petzoldt@tu-dresden.de

• `rkMethod` for a list of available Runge-Kutta parameter sets,

• `rk` for the more general Runge-Code,

• `lsoda`, `lsode`, `lsodes`, `lsodar`, `vode`, `daspk` for solvers of the Livermore family,

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

• `ode.band` for solving models with a banded Jacobian,

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

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

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

• `dede` for integrating models with delay differential equations,

`diagnostics` to print diagnostic messages.

### Examples

```## =======================================================================
## Example: Analytical and numerical solutions of logistic growth
## =======================================================================

## the derivative of the logistic
logist <- function(t, x, parms) {
with(as.list(parms), {
dx <- r * x * (1 - x/K)
list(dx)
})
}

time  <- 0:100
N0    <- 0.1; r <- 0.5; K <- 100
parms <- c(r = r, K = K)
x <- c(N = N0)

## analytical solution
plot(time, K/(1 + (K/N0-1) * exp(-r*time)), ylim = c(0, 120),
type = "l", col = "red", lwd = 2)

## reasonable numerical solution with rk4
time <- seq(0, 100, 2)
out <- as.data.frame(rk4(x, time, logist, parms))
points(out\$time, out\$N, pch = 16, col = "blue", cex = 0.5)

## same time step with euler, systematic under-estimation
time <- seq(0, 100, 2)
out <- as.data.frame(euler(x, time, logist, parms))
points(out\$time, out\$N, pch = 1)

## unstable result
time <- seq(0, 100, 4)
out <- as.data.frame(euler(x, time, logist, parms))
points(out\$time, out\$N, pch = 8, cex = 0.5)

## method with automatic time step
out <- as.data.frame(lsoda(x, time, logist, parms))
points(out\$time, out\$N, pch = 1, col = "green")

legend("bottomright",
c("analytical","rk4, h=2", "euler, h=2",
"euler, h=4", "lsoda"),
lty = c(1, NA, NA, NA, NA), lwd = c(2, 1, 1, 1, 1),
pch = c(NA, 16, 1, 8, 1),
col = c("red", "blue", "black", "black", "green"))
```

[Package deSolve version 1.30 Index]