## =======================================================================
## Example1: Predator-Prey Lotka-Volterra model (with logistic prey)
## =======================================================================
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 <- ode(yini, times, LVmod, pars)
summary(out)
## Default plot method
plot(out)
## User specified plotting
matplot(out[ , 1], 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: Substrate-Producer-Consumer Lotka-Volterra model
## =======================================================================
## Note:
## Function sigimp passed as an argument (input) to model
## (see also lsoda and rk examples)
SPCmod <- function(t, x, parms, input) {
with(as.list(c(parms, x)), {
import <- input(t)
dS <- import - b*S*P + g*C # substrate
dP <- c*S*P - d*C*P # producer
dC <- e*P*C - f*C # consumer
res <- c(dS, dP, dC)
list(res)
})
}
## The parameters
parms <- c(b = 0.001, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)
## vector of timesteps
times <- seq(0, 200, length = 101)
## external signal with rectangle impulse
signal <- data.frame(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, C = 1)
## Solve model
out <- ode(y = xstart, times = times,
func = SPCmod, parms = parms, input = sigimp)
## Default plot method
plot(out)
## User specified plotting
mf <- 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", "C"))
plot(out[,"P"], out[,"C"], type = "l", lwd = 2, xlab = "producer",
ylab = "consumer")
par(mfrow = mf)
## =======================================================================
## Example3: Discrete time model - using method = "iteration"
## The host-parasitoid model from Soetaert and Herman, 2009,
## Springer - p. 284.
## =======================================================================
Parasite <- function(t, y, ks) {
P <- y[1]
H <- y[2]
f <- A * P / (ks + H)
Pnew <- H * (1 - exp(-f))
Hnew <- H * exp(rH * (1 - H) - f)
list (c(Pnew, Hnew))
}
rH <- 2.82 # rate of increase
A <- 100 # attack rate
ks <- 15 # half-saturation density
out <- ode(func = Parasite, y = c(P = 0.5, H = 0.5), times = 0:50, parms = ks,
method = "iteration")
out2<- ode(func = Parasite, y = c(P = 0.5, H = 0.5), times = 0:50, parms = 25,
method = "iteration")
out3<- ode(func = Parasite, y = c(P = 0.5, H = 0.5), times = 0:50, parms = 35,
method = "iteration")
## Plot all 3 scenarios in one figure
plot(out, out2, out3, lty = 1, lwd = 2)
## Same like "out", but *output* every two steps
## hini = 1 ensures that the same *internal* timestep of 1 is used
outb <- ode(func = Parasite, y = c(P = 0.5, H = 0.5),
times = seq(0, 50, 2), hini = 1, parms = ks,
method = "iteration")
plot(out, outb, type = c("l", "p"))
## =======================================================================
## Example4: Playing with the Jacobian options - see e.g. lsoda help page
##
## IMPORTANT: The following example is temporarily broken because of
## incompatibility with R 3.0 on some systems.
## A fix is on the way.
## =======================================================================
## a stiff equation, exponential decay, run 500 times
stiff <- function(t, y, p) { # y and r are a 500-valued vector
list(- r * y)
}
N <- 500
r <- runif(N, 15, 20)
yini <- runif(N, 1, 40)
times <- 0:10
## Using the default
print(system.time(
out <- ode(y = yini, parms = NULL, times = times, func = stiff)
))
# diagnostics(out) shows that the method used = bdf (2), so it it stiff
## Specify that the Jacobian is banded, with nonzero values on the
## diagonal, i.e. the bandwidth up and down = 0
print(system.time(
out2 <- ode(y = yini, parms = NULL, times = times, func = stiff,
jactype = "bandint", bandup = 0, banddown = 0)
))
## Now we also specify the Jacobian function
jacob <- function(t, y, p) -r
print(system.time(
out3 <- ode(y = yini, parms = NULL, times = times, func = stiff,
jacfunc = jacob, jactype = "bandusr",
bandup = 0, banddown = 0)
))
## The larger the value of N, the larger the time gain...Run the code above in your browser using DataLab