deSolve (version 1.2-1)

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

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,...).

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.

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