Solving initial value problems for nonstiff systems of firstorder ordinary differential equations (ODEs).
The R function rk
is a toplevel function that provides
interfaces to a collection of common explicit onestep solvers of the
RungeKutta family with fixed or variable time steps.
The system of ODE's is written as an R function (which may, of
course, use .C
, .Fortran
,
.Call
, etc., to call foreign code) or be defined in
compiled code that has been dynamically loaded. A vector of
parameters is passed to the ODEs, so the solver may be used as part of
a modeling package for ODEs, or for parameter estimation using any
appropriate modeling tool for nonlinear models in R such as
optim
, nls
, nlm
or
nlme
1 2 3 4 5 6 7  rk(y, times, func, parms, rtol = 1e6, atol = 1e6,
verbose = FALSE, tcrit = NULL, hmin = 0, hmax = NULL,
hini = hmax, ynames = TRUE, method = rkMethod("rk45dp7", ... ),
maxsteps = 5000, dllname = NULL, initfunc = dllname,
initpar = parms, rpar = NULL, ipar = NULL,
nout = 0, outnames = NULL, forcings = NULL,
initforc = NULL, fcontrol = NULL, events = NULL, ...)

y 
the initial (state) values for the ODE system. If 
times 
times at which explicit estimates for 
func 
either an Rfunction 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 The return value of If 
parms 
vector or list of parameters used in 
rtol 
relative error tolerance, either a scalar or an array as
long as 
atol 
absolute error tolerance, either a scalar or an array as
long as 
tcrit 
if not 
verbose 
a logical value that, when TRUE, triggers more verbose output from the ODE solver. 
hmin 
an optional minimum value of the integration stepsize. In
special situations this parameter may speed up computations with the
cost of precision. Don't use 
hmax 
an optional maximum value of the integration stepsize. If
not specified, 
hini 
initial step size to be attempted; if 0, the initial step
size is determined automatically by solvers with flexible time step.
For fixed step methods, setting 
ynames 
if 
method 
the integrator to use. This can either be a string
constant naming one of the predefined methods or a call to function

maxsteps 
average maximal number of steps per output interval
taken by the solver. This argument is defined such as to ensure
compatibility with the Livermoresolvers. 
dllname 
a string giving the name of the shared library
(without extension) that contains all the compiled function or
subroutine definitions refered to in 
initfunc 
if not 
initpar 
only when ‘dllname’ is specified and an
initialisation function 
rpar 
only when ‘dllname’ is specified: a vector with
double precision values passed to the dllfunctions whose names are
specified by 
ipar 
only when ‘dllname’ is specified: a vector with
integer values passed to the dllfunctions whose names are specified
by 
nout 
only used if 
outnames 
only used if ‘dllname’ is specified and

forcings 
only used if ‘dllname’ is specified: a list with
the forcing function data sets, each present as a twocolumned matrix,
with (time,value); interpolation outside the interval
[min( See forcings or package vignette 
initforc 
if not 
fcontrol 
A list of control parameters for the forcing functions.
See forcings or vignette 
events 
A list that specifies events, i.e. when the value of a
state variable is suddenly changed. See events for more information.
Not also that if events are specified, then polynomial interpolation
is switched off and integration takes place from one external time step
to the next, with an internal step size less than or equal the difference
of two adjacent points of 
... 
additional arguments passed to 
Function rk
is a generalized implementation that can be used to
evaluate different solvers of the RungeKutta family of explicit ODE
solvers. A predefined set of common method parameters is in function
rkMethod
which also allows to supply userdefined
Butcher tables.
The input parameters rtol
, and atol
determine the error
control performed by the solver. The solver will control the vector
of estimated local errors in y, according to an inequality of
the form maxnorm of ( e/ewt ) <= 1, where
ewt is a vector of positive error weights. The values of
rtol
and atol
should all be nonnegative. The form of
ewt is:
\bold{rtol} * abs(\bold{y}) + \bold{atol}
where multiplication of two vectors is elementbyelement.
Models can be defined in R as a usersupplied
Rfunction, that must be called 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.
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 second element contains output variables that are
required at each point in time. Examples are given below.
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.
Arguments rpar
and ipar
are provided for compatibility
with lsoda
.
Starting with version 1.8 implicit RungeKutta methods are also
supported by this general rk
interface, however their
implementation is still experimental. Instead of this you may
consider radau
for a specific full implementation of an
implicit RungeKutta method.
Thomas Petzoldt thomas.petzoldt@tudresden.de
Butcher, J. C. (1987) The numerical analysis of ordinary differential equations, RungeKutta and general linear methods, Wiley, Chichester and New York.
EngelnMuellges, G. and Reutter, F. (1996) Numerik Algorithmen: Entscheidungshilfe zur Auswahl und Nutzung. VDI Verlag, Duesseldorf.
Hindmarsh, Alan C. (1983) ODEPACK, A Systematized Collection of ODE Solvers; in p.55–64 of Stepleman, R.W. et al.[ed.] (1983) Scientific Computing, NorthHolland, Amsterdam.
Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (2007) Numerical Recipes in C. Cambridge University Press.
For most practical cases, solvers of the Livermore family (i.e. the ODEPACK solvers, see below) are superior. Some of them are also suitable for stiff ODEs, differential algebraic equations (DAEs), or partial differential equations (PDEs).
rkMethod
for a list of available RungeKutta
parameter sets,
rk4
and euler
for special
versions without interpolation (and less overhead),
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 1D models,
ode.2D
for integrating 2D models,
ode.3D
for integrating 3D models,
diagnostics
to print diagnostic messages.
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 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65  ## =======================================================================
## Example: Resourceproducerconsumer LotkaVolterra model
## =======================================================================
## Notes:
##  Parameters are a list, names accessible via "with" function
##  Function sigimp passed as an argument (input) to model
## (see also ode and lsoda examples)
SPCmod < function(t, x, parms, input) {
with(as.list(c(parms, x)), {
import < input(t)
dS < import  b*S*P + g*C # substrate
dP < c*S*P  d*C*P # producer
dC < e*P*C  f*C # consumer
res < c(dS, dP, dC)
list(res)
})
}
## The parameters
parms < c(b = 0.001, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)
## vector of timesteps
times < seq(0, 200, length = 101)
## external signal with rectangle impulse
signal < data.frame(times = times,
import = rep(0, length(times)))
signal$import[signal$times >= 10 & signal$times <= 11] < 0.2
sigimp < approxfun(signal$times, signal$import, rule = 2)
## Start values for steady state
xstart < c(S = 1, P = 1, C = 1)
## Euler method
out1 < rk(xstart, times, SPCmod, parms, hini = 0.1,
input = sigimp, method = "euler")
## classical RungeKutta 4th order
out2 < rk(xstart, times, SPCmod, parms, hini = 1,
input = sigimp, method = "rk4")
## DormandPrince method of order 5(4)
out3 < rk(xstart, times, SPCmod, parms, hmax = 1,
input = sigimp, method = "rk45dp7")
mf < par("mfrow")
## deSolve plot method for comparing scenarios
plot(out1, out2, out3, which = c("S", "P", "C"),
main = c ("Substrate", "Producer", "Consumer"),
col =c("black", "red", "green"),
lty = c("solid", "dotted", "dotted"), lwd = c(1, 2, 1))
## userspecified plot function
plot (out1[,"P"], out1[,"C"], type = "l", xlab = "Producer", ylab = "Consumer")
lines(out2[,"P"], out2[,"C"], col = "red", lty = "dotted", lwd = 2)
lines(out3[,"P"], out3[,"C"], col = "green", lty = "dotted")
legend("center", legend = c("euler", "rk4", "rk45dp7"),
lty = c(1, 3, 3), lwd = c(1, 2, 1),
col = c("black", "red", "green"))
par(mfrow = mf)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.