rk4 {deSolve} | R Documentation |
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.
euler(y, times, func, parms, verbose = FALSE, ynames = TRUE, dllname = NULL, initfunc = dllname, initpar = parms, rpar = NULL, ipar = NULL, nout = 0, outnames = NULL, ...) rk4(y, times, func, parms, verbose = FALSE, ynames = TRUE, dllname = NULL, initfunc = dllname, initpar = parms, rpar = NULL, ipar = NULL, nout = 0, outnames = NULL, ...)
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:
yprime = func(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 lsoda() is called. See package vignette
"compiledCode" for more
details.
|
parms |
vector or list of parameters used in func |
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 and
jacfunc . 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 calculed 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.
|
... |
additional arguments passed to func allowing this
to be a generic function.
|
rk4
and euler
are special versions of the two fixed step
solvers with less overhead and less functionality (e.g. no interpolation)
compared to the generic Runge-Kutta codes called by rk
.
If you need different internal and external time steps, you may 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.
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.
Thomas Petzoldt thomas.petzoldt@tu-dresden.de
rkMethod
for a list of available Runge-Kutta
parameter sets,
rk
and euler
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,
diagnostics
to print diagnostic messages.
## =============================================================== ## 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] * (1 - x[1]/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 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, 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"))