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.

`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]