rk4 {deSolve} | R Documentation |
Solving initial value problems for systems of first-order ordinary differential equations (ODEs) using the classical Runge-Kutta 4th order integration.
rk4(y, times, func, parms, hini=min(diff(times)), ...) euler(y, times, func, parms, hini=min(diff(times)), ...)
y |
the initial 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 |
an R-function that computes the values of the
derivatives in the ode system (the model defininition) at time
t.
The R-function func 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, and 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 second element is a vector of global values that are required at
each point in times . |
parms |
vector or list holding the parameters used in func . |
hini |
length of the (fixed) internal time step. |
... |
additional arguments, allowing this to be a generic function |
See help pages of rk
and rkMethod
for details.
rk4
and euler
are primarily provided for didactic reasons.
For most practical cases, solvers with flexible timestep (e.g. "rk45dp7")
and especially solvers of the Livermore family (e.g. lsoda
) are superior.
Thomas Petzoldt thomas.petzoldt@tu-dresden.de
rk
,
ode
,
rkMethod
,
lsoda
,lsode
, lsodes
,
lsodar
, daspk
## numeric solution of logistic growth 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) ## good solution out <- as.data.frame(rk4(x, time, logist, parms, hini=2)) points(out$time, out$N, pch=16, col="blue", cex=0.5) ## same time step, systematic under-estimation out <- as.data.frame(euler(x, time, logist, parms, hini=2)) points(out$time, out$N, pch=1) ## unstable result out <- as.data.frame(euler(x, time, logist, parms, hini=4)) 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")