deSolve (version 1.2-3)

rk4: 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, hini = min(diff(times)), verbose = FALSE, ...)
rk4(y, times, func, parms, hini = min(diff(times)), verbose = FALSE, ...)

Arguments

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

parms

vector or list holding the parameters used in func.

hini

length of the (fixed) internal time step. Setting hini = 0, forces setting of internal time steps identically to external time steps provided by times.

verbose

a logical value that, when TRUE, triggers more verbose output from the ODE solver.

...

additional arguments, passed to rk

Details

See help pages of rk and rkMethod for details.

See Also

rk, ode, rkMethod, lsoda,lsode, lsodes, lsodar, daspk

Examples

## 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)

## reasonable numerical 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")
legend("bottomright",
  c("analytical","rk4, hini=2", "euler, hini=2",
    "euler, hini=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"))