rk4 {odesolve}  R Documentation 
Solving initial value problems for
systems of firstorder ordinary differential equations
(ODEs) using the classical RungeKutta 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 nonlinear models
in R such as optim
, nls
,
nlm
or nlme
.
rk4(y, times, func, parms, ...)
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 usersupplied function that computes the values of the
derivatives in the ode system (the model defininition) at time
t.
The usersupplied 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 
The method is implemented primarily for didactic purposes. Please use lsoda for your real work!
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.
Thomas Petzoldt thomas.petzoldt@tudresden.de
Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (1992) Numerical Recipes in C. Cambridge University Press.
## A simple resource limited LotkaVolterraModel 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")