rk4 {odesolve}R Documentation

Solve System of ODE (ordinary differential equation)s by classical Runge-Kutta 4th order integration.

Description

Solving initial value problems for systems of first-order ordinary differential equations (ODEs) using the classical Runge-Kutta 4th order integration. The system of ODE's may be written as an R function (which may, of course, use .C, .Fortran, .Call, etc., to call foreign code). 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 non-linear models in R such as optim, nls, nlm or nlme.

Usage

  rk4(y, times, func, parms, ...)

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 a user-supplied function that computes the values of the derivatives in the ode system (the model defininition) at time t. The user-supplied function func 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, and parms is a vector of parameters (which may have a names attribute, desirable in a large 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 is a vector (possibly with a names attribute) of global values that are required at each point in times.
parms vector or list holding the parameters used in func that should be modifiable without rewriting the function.
... additional arguments, allowing this to be a generic function

Details

The method is implemented primarily for didactic purposes. Please use lsoda for your real work!

Value

A matrix 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 second element of the return from func, plus and additional column for the time value. There will be a row for each element in times If y has a names attribute, it will be used to label the columns of the output value.

Author(s)

Thomas Petzoldt thomas.petzoldt@tu-dresden.de

References

Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (1992) Numerical Recipes in C. Cambridge University Press.

See Also

lsoda

Examples

## A simple resource limited Lotka-Volterra-Model
lvmodel <- function(t, x, parms) {
  s <- x[1] # substrate
  p <- x[2] # producer
  k <- x[3] # consumer
  with(as.list(parms),{
    import <- approx(signal$times, signal$import, t)$y
    ds <- import - b*s*p + g*k
    dp <- c*s*p  - d*k*p
    dk <- e*p*k  - f*k
    res<-c(ds, dp, dk)
    list(res)
  })
}

## vector of timesteps
times  <- seq(0, 100, length=101)

## external signal with rectangle impulse
signal <- as.data.frame(list(times = times,
                            import = rep(0,length(times))))

signal$import[signal$times >= 10 & signal$times <=11] <- 0.2

## Parameters for steady state conditions
parms  <- c(a=0.0, b=0.0, c=0.1, d=0.1, e=0.1, f=0.1, g=0.0)

## Start values for steady state
y<-xstart <- c(s=1, p=1, k=1)

## Classical RK4 with fixed time step
out1  <- as.data.frame(rk4(xstart, times, lvmodel, parms))

## LSODA (default step size)
out2 <- as.data.frame(lsoda(xstart, times, lvmodel, parms))

## LSODA: with fixed maximum time step
out3 <- as.data.frame(lsoda(xstart, times, lvmodel, parms, hmax=1))

par(mfrow=c(2,2))
plot (out1$time, out1$s,  type="l",   ylim=c(0,3))
lines(out2$time, out2$s, col="red",   lty="dotted")
lines(out3$time, out3$s, col="green", lty="dotted")

plot (out1$time, out1$p, type="l",    ylim=c(0,3))
lines(out2$time, out2$p, col="red",   lty="dotted")
lines(out3$time, out3$p, col="green", lty="dotted")

plot (out1$time, out1$k, type="l",    ylim=c(0,3))
lines(out2$time, out2$k, col="red",   lty="dotted")
lines(out3$time, out3$k, col="green", lty="dotted")

plot (out1$p, out1$k, type="l")
lines(out2$p, out2$k, col="red",   lty="dotted")
lines(out3$p, out3$k, col="green", lty="dotted")

[Package odesolve version 0.5-20 Index]