ode {deSolve} | R Documentation |
Solves a system of ordinary differential equations.
ode(y, times, func, parms, method = c("lsoda", "lsode", "lsodes", "lsodar", "vode", "daspk", "euler", "rk4", "ode23", "ode45"), ...)
y |
the initial (state) values for the ODE system, a vector. If
y has a name attribute, the names will be used to label the
output matrix.
|
times |
time sequence for which output is wanted; the first
value of times must be the initial time.
|
func |
either an R-function that computes the values of the
derivatives in the ODE system (the model definition) at time t, or a
character string giving the name of a compiled function in a
dynamically loaded shared library.
If func is an R-function, it must be defined 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. If the initial values y has a
names attribute, the names will be available inside func .
parms is a vector or list of parameters; ... (optional) are
any other arguments passed to the function.
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 next elements are global values
that are required at each point in times .The derivatives
should be specified in the same order as the state variables y .
If func is a string, then dllname must give the name
of the shared library (without extension) which must be loaded
before lsoda() is called. See package vignette
"compiledCode" for more details.
|
parms |
parameters passed to func . |
method |
the integrator to use, either a string ("lsoda" ,
"lsode" , "lsodes" , "lsodar" , "vode" ,
"daspk" , "euler" , "rk4" , "ode23" or
"ode45" ) or a function that performs integration, or a list
of class rkMethod . |
... |
additional arguments passed to the integrator. |
This is simply a wrapper around the various ode solvers.
See package vignette for information about specifying the model in compiled code.
See the selected integrator for the additional options.
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 an
additional column (the first) for the time value. There will be one
row for each element in times
unless the integrator returns
with an unrecoverable error. If y
has a names attribute, it
will be used to label the columns of the output value.
The output will have the attributes istate
, and rstate
,
two vectors with several useful elements. The first element of istate
returns the conditions under which the last call to the integrator
returned. Normal is istate = 2
. If verbose = TRUE
, the
settings of istate and rstate will be written to the screen. See the
help for the selected integrator for details.
Karline Soetaert <k.soetaert@nioo.knaw.nl>
ode.band
for solving models with a banded
Jacobian,
ode.1D
for integrating 1-D models,
ode.2D
for integrating 2-D models,
ode.3D
for integrating 3-D models,
aquaphy
, ccl4model
, where
ode
is used,
lsoda
, lsode
,
lsodes
, lsodar
, vode
,
daspk
,
rk
, rkMethod
diagnostics
to print diagnostic messages.
## ============================================= ## Example1: Predator-Prey Lotka-Volterra model ## ============================================= LVmod <- function(Time, State, Pars) { with(as.list(c(State, Pars)), { Ingestion <- rIng * Prey*Predator GrowthPrey <- rGrow * Prey*(1-Prey/K) MortPredator <- rMort * Predator dPrey <- GrowthPrey - Ingestion dPredator <- Ingestion*assEff -MortPredator return(list(c(dPrey, dPredator))) }) } pars <- c(rIng = 0.2, # /day, rate of ingestion rGrow = 1.0, # /day, growth rate of prey rMort = 0.2 , # /day, mortality rate of predator assEff = 0.5, # -, assimilation efficiency K = 10 ) # mmol/m3, carrying capacity yini <- c(Prey = 1, Predator = 2) times <- seq(0, 200, by = 1) out <- as.data.frame(ode(func = LVmod, y = yini, parms = pars, times = times)) matplot(out$time,out[,2:3],type = "l",xlab = "time",ylab = "Conc", main = "Lotka-Volterra",lwd = 2) legend("topright", c("prey", "predator"), col = 1:2, lty = 1:2) ## ========================================================== ## Example2: Resource-producer-consumer Lotka-Volterra model ## ========================================================== ## Note: ## 1. parameter and state variable names made ## accessible via "with" statement ## 2. function sigimp passed as an argument (input) to model ## (see also lsoda and rk examples) lvmodel <- function(t, x, parms, input) { with(as.list(c(parms, x)), { import <- input(t) dS <- import - b*S*P + g*K # substrate dP <- c*S*P - d*K*P # producer dK <- e*P*K - f*K # consumer res <- c(dS, dP, dK) list(res) }) } ## The parameters parms <- c(b = 0.0, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0) ## 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 sigimp <- approxfun(signal$times, signal$import, rule = 2) ## Start values for steady state xstart <- c(S = 1, P = 1, K = 1) ## Solve model out <- as.data.frame(ode(y = xstart,times = times, func = lvmodel, parms, input = sigimp)) ## Display results par(mfrow=c(1,2)) matplot(out[,1], out[,2:4], type="l", xlab="time", ylab="state") legend("topright", col = 1:3, lty = 1:3, legend = c("S", "P", "K")) plot(out$P, out$K, type = "l", lwd = 2, xlab = "producer", ylab = "consumer")