Nothing
inst/doc/deSolve.R
In deSolve: Solvers for Initial Value Problems of Differential Equations ('ODE', 'DAE', 'DDE')
### R code from vignette source 'deSolve.Rnw'
###################################################
### code chunk number 1: preliminaries
###################################################
library("deSolve")
options(prompt = "> ")
options(width=70)
###################################################
### code chunk number 2: deSolve.Rnw:180-183
###################################################
parameters <- c(a = -8/3,
b = -10,
c = 28)
###################################################
### code chunk number 3: deSolve.Rnw:191-194
###################################################
state <- c(X = 1,
Y = 1,
Z = 1)
###################################################
### code chunk number 4: deSolve.Rnw:221-232
###################################################
Lorenz<-function(t, state, parameters) {
with(as.list(c(state, parameters)),{
# rate of change
dX <- a*X + Y*Z
dY <- b * (Y-Z)
dZ <- -X*Y + c*Y - Z
# return the rate of change
list(c(dX, dY, dZ))
}) # end with(as.list ...
}
###################################################
### code chunk number 5: deSolve.Rnw:242-243
###################################################
times <- seq(0, 100, by = 0.01)
###################################################
### code chunk number 6: deSolve.Rnw:258-261
###################################################
library(deSolve)
out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
head(out)
###################################################
### code chunk number 7: ode
###################################################
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")
plot(out[, "X"], out[, "Z"], pch = ".")
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)
###################################################
### code chunk number 8: figode
###################################################
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")
plot(out[, "X"], out[, "Z"], pch = ".")
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)
###################################################
### code chunk number 9: deSolve.Rnw:315-318
###################################################
outb <- radau(state, times, Lorenz, parameters, atol = 1e-4, rtol = 1e-4)
outc <- ode(state, times, Lorenz, parameters, method = "radau",
atol = 1e-4, rtol = 1e-4)
###################################################
### code chunk number 10: deSolve.Rnw:334-340
###################################################
print(system.time(out1 <- rk4 (state, times, Lorenz, parameters)))
print(system.time(out2 <- lsode (state, times, Lorenz, parameters)))
print(system.time(out <- lsoda (state, times, Lorenz, parameters)))
print(system.time(out <- lsodes(state, times, Lorenz, parameters)))
print(system.time(out <- daspk (state, times, Lorenz, parameters)))
print(system.time(out <- vode (state, times, Lorenz, parameters)))
###################################################
### code chunk number 11: deSolve.Rnw:358-359
###################################################
rkMethod()
###################################################
### code chunk number 12: deSolve.Rnw:368-369
###################################################
rkMethod("rk23")
###################################################
### code chunk number 13: deSolve.Rnw:382-403
###################################################
func <- function(t, x, parms) {
with(as.list(c(parms, x)),{
dP <- a * P - b * C * P
dC <- b * P * C - c * C
res <- c(dP, dC)
list(res)
})
}
rKnew <- rkMethod(ID = "midpoint",
varstep = FALSE,
A = c(0, 1/2),
b1 = c(0, 1),
c = c(0, 1/2),
stage = 2,
Qerr = 1
)
out <- ode(y = c(P = 2, C = 1), times = 0:100, func,
parms = c(a = 0.1, b = 0.1, c = 0.1), method = rKnew)
head(out)
###################################################
### code chunk number 14: deSolve.Rnw:437-439
###################################################
diagnostics(out1)
diagnostics(out2)
###################################################
### code chunk number 15: deSolve.Rnw:443-444
###################################################
summary(out1)
###################################################
### code chunk number 16: deSolve.Rnw:518-526
###################################################
Aphid <- function(t, APHIDS, parameters) {
deltax <- c (0.5, rep(1, numboxes - 1), 0.5)
Flux <- -D * diff(c(0, APHIDS, 0)) / deltax
dAPHIDS <- -diff(Flux) / delx + APHIDS * r
# the return value
list(dAPHIDS )
} # end
###################################################
### code chunk number 17: deSolve.Rnw:531-538
###################################################
D <- 0.3 # m2/day diffusion rate
r <- 0.01 # /day net growth rate
delx <- 1 # m thickness of boxes
numboxes <- 60
# distance of boxes on plant, m, 1 m intervals
Distance <- seq(from = 0.5, by = delx, length.out = numboxes)
###################################################
### code chunk number 18: deSolve.Rnw:543-547
###################################################
# Initial conditions: # ind/m2
APHIDS <- rep(0, times = numboxes)
APHIDS[30:31] <- 1
state <- c(APHIDS = APHIDS) # initialise state variables
###################################################
### code chunk number 19: deSolve.Rnw:554-558
###################################################
times <-seq(0, 200, by = 1)
print(system.time(
out <- ode.1D(state, times, Aphid, parms = 0, nspec = 1, names = "Aphid")
))
###################################################
### code chunk number 20: deSolve.Rnw:564-565
###################################################
head(out[,1:5])
###################################################
### code chunk number 21: deSolve.Rnw:569-570
###################################################
summary(out)
###################################################
### code chunk number 22: deSolve.Rnw:603-605
###################################################
data <- cbind(dist = c(0,10, 20, 30, 40, 50, 60),
Aphid = c(0,0.1,0.25,0.5,0.25,0.1,0))
###################################################
### code chunk number 23: matplot1d
###################################################
par (mfrow = c(1,2))
matplot.1D(out, grid = Distance, type = "l", mfrow = NULL,
subset = time %in% seq(0, 200, by = 10),
obs = data, obspar = list(pch = 18, cex = 2, col="red"))
plot.1D(out, grid = Distance, type = "l", mfrow = NULL,
subset = time == 100,
obs = data, obspar = list(pch = 18, cex = 2, col="red"))
###################################################
### code chunk number 24: matplot1d
###################################################
par (mfrow = c(1,2))
matplot.1D(out, grid = Distance, type = "l", mfrow = NULL,
subset = time %in% seq(0, 200, by = 10),
obs = data, obspar = list(pch = 18, cex = 2, col="red"))
plot.1D(out, grid = Distance, type = "l", mfrow = NULL,
subset = time == 100,
obs = data, obspar = list(pch = 18, cex = 2, col="red"))
###################################################
### code chunk number 25: deSolve.Rnw:671-686
###################################################
daefun <- function(t, y, dy, parameters) {
res1 <- dy[1] + y[1] - y[2]
res2 <- y[2] * y[1] - t
list(c(res1, res2))
}
library(deSolve)
yini <- c(1, 0)
dyini <- c(1, 0)
times <- seq(0, 10, 0.1)
## solver
system.time(out <- daspk(y = yini, dy = dyini,
times = times, res = daefun, parms = 0))
###################################################
### code chunk number 26: dae
###################################################
matplot(out[,1], out[,2:3], type = "l", lwd = 2,
main = "dae", xlab = "time", ylab = "y")
###################################################
### code chunk number 27: figdae
###################################################
matplot(out[,1], out[,2:3], type = "l", lwd = 2,
main = "dae", xlab = "time", ylab = "y")
###################################################
### code chunk number 28: deSolve.Rnw:719-729
###################################################
pendulum <- function (t, Y, parms) {
with (as.list(Y),
list(c(u,
v,
-lam * x,
-lam * y - 9.8,
x^2 + y^2 -1
))
)
}
###################################################
### code chunk number 29: deSolve.Rnw:732-733
###################################################
yini <- c(x = 1, y = 0, u = 0, v = 1, lam = 1)
###################################################
### code chunk number 30: deSolve.Rnw:736-739
###################################################
M <- diag(nrow = 5)
M[5, 5] <- 0
M
###################################################
### code chunk number 31: deSolve.Rnw:743-747
###################################################
index <- c(2, 2, 1)
times <- seq(from = 0, to = 10, by = 0.01)
out <- radau (y = yini, func = pendulum, parms = NULL,
times = times, mass = M, nind = index)
###################################################
### code chunk number 32: pendulum
###################################################
plot(out, type = "l", lwd = 2)
plot(out[, c("x", "y")], type = "l", lwd = 2)
###################################################
### code chunk number 33: pendulum
###################################################
plot(out, type = "l", lwd = 2)
plot(out[, c("x", "y")], type = "l", lwd = 2)
###################################################
### code chunk number 34: deSolve.Rnw:781-794
###################################################
ZODE2 <- function(Time, State, Pars) {
with(as.list(State), {
df <- 1i * f
dg <- -1i * g * g * f
return(list(c(df, dg)))
})
}
yini <- c(f = 1+0i, g = 1/2.1+0i)
times <- seq(0, 2 * pi, length = 100)
out <- zvode(func = ZODE2, y = yini, parms = NULL, times = times,
atol = 1e-10, rtol = 1e-10)
###################################################
### code chunk number 35: deSolve.Rnw:806-808
###################################################
analytical <- cbind(f = exp(1i*times), g = 1/(exp(1i*times)+1.1))
tail(cbind(out[,2], analytical[,1]))
###################################################
### code chunk number 36: deSolve.Rnw:821-832
###################################################
f1 <- function (t, y, parms) {
ydot <- vector(len = 5)
ydot[1] <- 0.1*y[1] -0.2*y[2]
ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3]
ydot[3] <- -0.3*y[2] +0.1*y[3] -0.2*y[4]
ydot[4] <- -0.3*y[3] +0.1*y[4] -0.2*y[5]
ydot[5] <- -0.3*y[4] +0.1*y[5]
return(list(ydot))
}
###################################################
### code chunk number 37: deSolve.Rnw:837-839
###################################################
yini <- 1:5
times <- 1:20
###################################################
### code chunk number 38: deSolve.Rnw:846-847
###################################################
out <- lsode(yini, times, f1, parms = 0, jactype = "fullint")
###################################################
### code chunk number 39: deSolve.Rnw:854-863
###################################################
fulljac <- function (t, y, parms) {
jac <- matrix(nrow = 5, ncol = 5, byrow = TRUE,
data = c(0.1, -0.2, 0 , 0 , 0 ,
-0.3, 0.1, -0.2, 0 , 0 ,
0 , -0.3, 0.1, -0.2, 0 ,
0 , 0 , -0.3, 0.1, -0.2,
0 , 0 , 0 , -0.3, 0.1))
return(jac)
}
###################################################
### code chunk number 40: deSolve.Rnw:868-870
###################################################
out2 <- lsode(yini, times, f1, parms = 0, jactype = "fullusr",
jacfunc = fulljac)
###################################################
### code chunk number 41: deSolve.Rnw:877-879
###################################################
out3 <- lsode(yini, times, f1, parms = 0, jactype = "bandint",
bandup = 1, banddown = 1)
###################################################
### code chunk number 42: deSolve.Rnw:884-891
###################################################
bandjac <- function (t, y, parms) {
jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE,
data = c( 0 , -0.2, -0.2, -0.2, -0.2,
0.1, 0.1, 0.1, 0.1, 0.1,
-0.3, -0.3, -0.3, -0.3, 0))
return(jac)
}
###################################################
### code chunk number 43: deSolve.Rnw:896-898
###################################################
out4 <- lsode(yini, times, f1, parms = 0, jactype = "bandusr",
jacfunc = bandjac, bandup = 1, banddown = 1)
###################################################
### code chunk number 44: deSolve.Rnw:904-905
###################################################
out5 <- lsode(yini, times, f1, parms = 0, mf = 10)
###################################################
### code chunk number 45: deSolve.Rnw:936-942
###################################################
eventmod <- function(t, var, parms) {
list(dvar = -0.1*var)
}
yini <- c(v1 = 1, v2 = 2)
times <- seq(0, 10, by = 0.1)
###################################################
### code chunk number 46: deSolve.Rnw:949-953
###################################################
eventdat <- data.frame(var = c("v1", "v2", "v2", "v1"), time = c(1, 1, 5, 9),
value = c(1, 2, 3, 4), method = c("add", "mult", "rep", "add"))
eventdat
###################################################
### code chunk number 47: deSolve.Rnw:958-960
###################################################
out <- ode(func = eventmod, y = yini, times = times, parms = NULL,
events = list(data = eventdat))
###################################################
### code chunk number 48: event1
###################################################
plot(out, type = "l", lwd = 2)
###################################################
### code chunk number 49: figevent1
###################################################
plot(out, type = "l", lwd = 2)
###################################################
### code chunk number 50: deSolve.Rnw:982-987
###################################################
ballode<- function(t, y, parms) {
dy1 <- y[2]
dy2 <- -9.8
list(c(dy1, dy2))
}
###################################################
### code chunk number 51: deSolve.Rnw:994-995
###################################################
root <- function(t, y, parms) y[1]
###################################################
### code chunk number 52: deSolve.Rnw:1000-1005
###################################################
event <- function(t, y, parms) {
y[1]<- 0
y[2]<- -0.9 * y[2]
return(y)
}
###################################################
### code chunk number 53: deSolve.Rnw:1011-1016
###################################################
yini <- c(height = 0, v = 20)
times <- seq(from = 0, to = 20, by = 0.01)
out <- lsode(times = times, y = yini, func = ballode, parms = NULL,
events = list(func = event, root = TRUE), rootfun = root)
###################################################
### code chunk number 54: event2
###################################################
plot(out, which = "height", type = "l",lwd = 2,
main = "bouncing ball", ylab = "height")
###################################################
### code chunk number 55: figevent2
###################################################
plot(out, which = "height", type = "l",lwd = 2,
main = "bouncing ball", ylab = "height")
###################################################
### code chunk number 56: deSolve.Rnw:1065-1067
###################################################
times <- seq(0, 1, 0.1)
eventtimes <- c(0.7, 0.9)
###################################################
### code chunk number 57: deSolve.Rnw:1072-1073
###################################################
eventtimes %in% times
###################################################
### code chunk number 58: deSolve.Rnw:1080-1082
###################################################
times2 <- round(times, 1)
times - times2
###################################################
### code chunk number 59: deSolve.Rnw:1093-1094
###################################################
eventtimes %in% times2
###################################################
### code chunk number 60: deSolve.Rnw:1099-1100
###################################################
all(eventtimes %in% times2)
###################################################
### code chunk number 61: deSolve.Rnw:1110-1113
###################################################
times <- 1:10
eventtimes <- c(1.3, 3.4, 4, 7.9, 8.5)
newtimes <- sort(unique(c(times, eventtimes)))
###################################################
### code chunk number 62: deSolve.Rnw:1119-1122
###################################################
times <- 1:10
eventtimes <- c(1.3, 3.4, 4, 7.9999999999999999, 8.5)
newtimes <- sort(c(eventtimes, cleanEventTimes(times, eventtimes)))
###################################################
### code chunk number 63: deSolve.Rnw:1151-1185
###################################################
library(deSolve)
#-----------------------------
# the derivative function
#-----------------------------
derivs <- function(t, y, parms) {
if (t < 0)
lag <- 19
else
lag <- lagvalue(t - 0.74)
dy <- r * y * (1 - lag/m)
list(dy, dy = dy)
}
#-----------------------------
# parameters
#-----------------------------
r <- 3.5; m <- 19
#-----------------------------
# initial values and times
#-----------------------------
yinit <- c(y = 19.001)
times <- seq(-0.74, 40, by = 0.01)
#-----------------------------
# solve the model
#-----------------------------
yout <- dede(y = yinit, times = times, func = derivs,
parms = NULL, atol = 1e-10)
###################################################
### code chunk number 64: dde
###################################################
plot(yout, which = 1, type = "l", lwd = 2,
main = "Lemming model", mfrow = c(1,2))
plot(yout[,2], yout[,3], xlab = "y", ylab = "dy", type = "l", lwd = 2)
###################################################
### code chunk number 65: figdde
###################################################
plot(yout, which = 1, type = "l", lwd = 2,
main = "Lemming model", mfrow = c(1,2))
plot(yout[,2], yout[,3], xlab = "y", ylab = "dy", type = "l", lwd = 2)
###################################################
### code chunk number 66: deSolve.Rnw:1222-1234
###################################################
Stages <- c("DS 1yr", "DS 2yr", "R small", "R medium", "R large", "F")
NumStages <- length(Stages)
# Population matrix
A <- matrix(nrow = NumStages, ncol = NumStages, byrow = TRUE, data = c(
0, 0, 0, 0, 0, 322.38,
0.966, 0, 0, 0, 0, 0 ,
0.013, 0.01, 0.125, 0, 0, 3.448 ,
0.007, 0, 0.125, 0.238, 0, 30.170,
0.008, 0, 0.038, 0.245, 0.167, 0.862 ,
0, 0, 0, 0.023, 0.75, 0 ) )
###################################################
### code chunk number 67: deSolve.Rnw:1239-1243
###################################################
Teasel <- function (t, y, p) {
yNew <- A %*% y
list (yNew / sum(yNew))
}
###################################################
### code chunk number 68: deSolve.Rnw:1246-1248
###################################################
out <- ode(func = Teasel, y = c(1, rep(0, 5) ), times = 0:50,
parms = 0, method = "iteration")
###################################################
### code chunk number 69: difference
###################################################
matplot(out[,1], out[,-1], main = "Teasel stage distribution", type = "l")
legend("topright", legend = Stages, lty = 1:6, col = 1:6)
###################################################
### code chunk number 70: difference
###################################################
matplot(out[,1], out[,-1], main = "Teasel stage distribution", type = "l")
legend("topright", legend = Stages, lty = 1:6, col = 1:6)
###################################################
### code chunk number 71: deSolve.Rnw:1291-1295
###################################################
library(deSolve)
combustion <- function (t, y, parms)
list(y^2 * (1-y) )
###################################################
### code chunk number 72: deSolve.Rnw:1297-1299
###################################################
yini <- 0.01
times <- 0 : 200
###################################################
### code chunk number 73: deSolve.Rnw:1301-1305
###################################################
out <- ode(times = times, y = yini, parms = 0, func = combustion)
out2 <- ode(times = times, y = yini*2, parms = 0, func = combustion)
out3 <- ode(times = times, y = yini*3, parms = 0, func = combustion)
out4 <- ode(times = times, y = yini*4, parms = 0, func = combustion)
###################################################
### code chunk number 74: plotdeSolve
###################################################
plot(out, out2, out3, out4, main = "combustion")
legend("bottomright", lty = 1:4, col = 1:4, legend = 1:4, title = "yini*i")
###################################################
### code chunk number 75: plotdeSolve
###################################################
plot(out, out2, out3, out4, main = "combustion")
legend("bottomright", lty = 1:4, col = 1:4, legend = 1:4, title = "yini*i")
###################################################
### code chunk number 76: deSolve.Rnw:1335-1336
###################################################
head(ccl4data)
###################################################
### code chunk number 77: deSolve.Rnw:1339-1342
###################################################
obs <- subset (ccl4data, animal == "A", c(time, ChamberConc))
names(obs) <- c("time", "CP")
head(obs)
###################################################
### code chunk number 78: deSolve.Rnw:1348-1362
###################################################
parms <- c(0.182, 4.0, 4.0, 0.08, 0.04, 0.74, 0.05, 0.15, 0.32, 16.17,
281.48, 13.3, 16.17, 5.487, 153.8, 0.04321671,
0.40272550, 951.46, 0.02, 1.0, 3.80000000)
yini <- c(AI = 21, AAM = 0, AT = 0, AF = 0, AL = 0, CLT = 0, AM = 0)
out <- ccl4model(times = seq(0, 6, by = 0.05), y = yini, parms = parms)
par2 <- parms
par2[1] <- 0.1
out2 <- ccl4model(times = seq(0, 6, by = 0.05), y = yini, parms = par2)
par3 <- parms
par3[1] <- 0.05
out3 <- ccl4model(times = seq(0, 6, by = 0.05), y = yini, parms = par3)
###################################################
### code chunk number 79: plotobs
###################################################
plot(out, out2, out3, which = c("AI", "MASS", "CP"),
col = c("black", "red", "green"), lwd = 2,
obs = obs, obspar = list(pch = 18, col = "blue", cex = 1.2))
legend("topright", lty = c(1,2,3,NA), pch = c(NA, NA, NA, 18),
col = c("black", "red", "green", "blue"), lwd = 2,
legend = c("par1", "par2", "par3", "obs"))
###################################################
### code chunk number 80: plotobs
###################################################
plot(out, out2, out3, which = c("AI", "MASS", "CP"),
col = c("black", "red", "green"), lwd = 2,
obs = obs, obspar = list(pch = 18, col = "blue", cex = 1.2))
legend("topright", lty = c(1,2,3,NA), pch = c(NA, NA, NA, 18),
col = c("black", "red", "green", "blue"), lwd = 2,
legend = c("par1", "par2", "par3", "obs"))
###################################################
### code chunk number 81: deSolve.Rnw:1388-1390
###################################################
obs2 <- data.frame(time = 6, MASS = 12)
obs2
###################################################
### code chunk number 82: obs2
###################################################
plot(out, out2, out3, lwd = 2,
obs = list(obs, obs2),
obspar = list(pch = c(16, 18), col = c("blue", "black"),
cex = c(1.2 , 2))
)
###################################################
### code chunk number 83: plotobs2
###################################################
plot(out, out2, out3, lwd = 2,
obs = list(obs, obs2),
obspar = list(pch = c(16, 18), col = c("blue", "black"),
cex = c(1.2 , 2))
)
###################################################
### code chunk number 84: hist
###################################################
hist(out, col = grey(seq(0, 1, by = 0.1)), mfrow = c(3, 4))
###################################################
### code chunk number 85: plothist
###################################################
hist(out, col = grey(seq(0, 1, by = 0.1)), mfrow = c(3, 4))
###################################################
### code chunk number 86: deSolve.Rnw:1449-1451
###################################################
options(prompt = " ")
options(continue = " ")
###################################################
### code chunk number 87: deSolve.Rnw:1454-1478
###################################################
lvmod <- function (time, state, parms, N, rr, ri, dr, dri) {
with (as.list(parms), {
PREY <- state[1:N]
PRED <- state[(N+1):(2*N)]
## Fluxes due to diffusion
## at internal and external boundaries: zero gradient
FluxPrey <- -Da * diff(c(PREY[1], PREY, PREY[N]))/dri
FluxPred <- -Da * diff(c(PRED[1], PRED, PRED[N]))/dri
## Biology: Lotka-Volterra model
Ingestion <- rIng * PREY * PRED
GrowthPrey <- rGrow * PREY * (1-PREY/cap)
MortPredator <- rMort * PRED
## Rate of change = Flux gradient + Biology
dPREY <- -diff(ri * FluxPrey)/rr/dr +
GrowthPrey - Ingestion
dPRED <- -diff(ri * FluxPred)/rr/dr +
Ingestion * assEff - MortPredator
return (list(c(dPREY, dPRED)))
})
}
###################################################
### code chunk number 88: deSolve.Rnw:1480-1482
###################################################
options(prompt = " ")
options(continue = " ")
###################################################
### code chunk number 89: deSolve.Rnw:1485-1499
###################################################
R <- 20 # total radius of surface, m
N <- 100 # 100 concentric circles
dr <- R/N # thickness of each layer
r <- seq(dr/2,by = dr,len = N) # distance of center to mid-layer
ri <- seq(0,by = dr,len = N+1) # distance to layer interface
dri <- dr # dispersion distances
parms <- c(Da = 0.05, # m2/d, dispersion coefficient
rIng = 0.2, # /day, rate of ingestion
rGrow = 1.0, # /day, growth rate of prey
rMort = 0.2 , # /day, mortality rate of pred
assEff = 0.5, # -, assimilation efficiency
cap = 10) # density, carrying capacity
###################################################
### code chunk number 90: deSolve.Rnw:1502-1512
###################################################
state <- rep(0, 2 * N)
state[1] <- state[N + 1] <- 10
times <- seq(0, 200, by = 1) # output wanted at these time intervals
print(system.time(
out <- ode.1D(y = state, times = times, func = lvmod, parms = parms,
nspec = 2, names = c("PREY", "PRED"),
N = N, rr = r, ri = ri, dr = dr, dri = dri)
))
###################################################
### code chunk number 91: deSolve.Rnw:1515-1516
###################################################
summary(out)
###################################################
### code chunk number 92: deSolve.Rnw:1520-1522
###################################################
p10 <- subset(out, select = "PREY", subset = time == 10)
head(p10, n = 5)
###################################################
### code chunk number 93: deSolve.Rnw:1568-1573
###################################################
Simple2D <- function(t, Y, par) {
y <- matrix(nrow = nx, ncol = ny, data = Y) # vector to 2-D matrix
dY <- - r_x2y2 * y # consumption
return(list(dY))
}
###################################################
### code chunk number 94: deSolve.Rnw:1577-1584
###################################################
dy <- dx <- 1 # grid size
nx <- ny <- 100
x <- seq (dx/2, by = dx, len = nx)
y <- seq (dy/2, by = dy, len = ny)
# in each grid cell: consumption depending on position
r_x2y2 <- outer(x, y, FUN=function(x,y) ((x-50)^2 + (y-50)^2)*1e-4)
###################################################
### code chunk number 95: deSolve.Rnw:1588-1591
###################################################
C <- matrix(nrow = nx, ncol = ny, 1)
ODE3 <- ode.2D(y = C, times = 1:100, func = Simple2D, parms = NULL,
dimens = c(nx, ny), names = "C", method = "ode45")
###################################################
### code chunk number 96: deSolve.Rnw:1594-1597
###################################################
summary(ODE3)
t50 <- matrix(nrow = nx, ncol = ny,
data = subset(ODE3, select = "C", subset = (time == 50)))
###################################################
### code chunk number 97: twoD
###################################################
par(mfrow = c(1, 2))
contour(x, y, r_x2y2, main = "consumption")
contour(x, y, t50, main = "Y(t = 50)")
###################################################
### code chunk number 98: twoD
###################################################
par(mfrow = c(1, 2))
contour(x, y, r_x2y2, main = "consumption")
contour(x, y, t50, main = "Y(t = 50)")
###################################################
### code chunk number 99: deSolve.Rnw:1627-1635
###################################################
PCmod <- function(t, x, parms) {
with(as.list(c(parms, x)), {
dP <- c*P - d*C*P # producer
dC <- e*P*C - f*C # consumer
res <- c(dP, dC)
list(res)
})
}
###################################################
### code chunk number 100: deSolve.Rnw:1642-1643
###################################################
parms <- c(c = 10, d = 0.1, e = 0.1, f = 0.1)
###################################################
### code chunk number 101: deSolve.Rnw:1649-1655
###################################################
xstart <- c(P = 0.5, C = 1)
times <- seq(0, 200, 0.1)
out <- ode(y = xstart, times = times,
func = PCmod, parms = parms)
tail(out)
###################################################
### code chunk number 102: deSolve.Rnw:1676-1680
###################################################
out <- ode(y = xstart,times = times, func = PCmod,
parms = parms, atol = 0)
matplot(out[,1], out[,2:3], type = "l",
xlab = "time", ylab = "Producer, Consumer")
###################################################
### code chunk number 103: deSolve.Rnw:1736-1760
###################################################
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(func = LVmod, y = yini,
parms = pars, times = times)
###################################################
### code chunk number 104: deSolve.Rnw:1772-1775
###################################################
pars["rIng"] <- 100
out2 <- ode(func = LVmod, y = yini,
parms = pars, times = times)
###################################################
### code chunk number 105: err
###################################################
plot(out2, type = "l", lwd = 2, main = "corrupt Lotka-Volterra model")
###################################################
### code chunk number 106: deSolve.Rnw:1824-1827
###################################################
pars["rIng"] <- 100
out3 <- ode(func = LVmod, y = yini, parms = pars,
times = times, method = "ode45", atol = 1e-14, rtol = 1e-14)
###################################################
### code chunk number 107: err
###################################################
plot(out2, type = "l", lwd = 2, main = "corrupt Lotka-Volterra model")
Try the deSolve package in your browser
Any scripts or data that you put into this service are public.
deSolve documentation built on Nov. 28, 2023, 1:11 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
### R code from vignette source 'deSolve.Rnw'
###################################################
### code chunk number 1: preliminaries
###################################################
library("deSolve")
options(prompt = "> ")
options(width=70)
###################################################
### code chunk number 2: deSolve.Rnw:180-183
###################################################
parameters <- c(a = -8/3,
b = -10,
c = 28)
###################################################
### code chunk number 3: deSolve.Rnw:191-194
###################################################
state <- c(X = 1,
Y = 1,
Z = 1)
###################################################
### code chunk number 4: deSolve.Rnw:221-232
###################################################
Lorenz<-function(t, state, parameters) {
with(as.list(c(state, parameters)),{
# rate of change
dX <- a*X + Y*Z
dY <- b * (Y-Z)
dZ <- -X*Y + c*Y - Z
# return the rate of change
list(c(dX, dY, dZ))
}) # end with(as.list ...
}
###################################################
### code chunk number 5: deSolve.Rnw:242-243
###################################################
times <- seq(0, 100, by = 0.01)
###################################################
### code chunk number 6: deSolve.Rnw:258-261
###################################################
library(deSolve)
out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
head(out)
###################################################
### code chunk number 7: ode
###################################################
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")
plot(out[, "X"], out[, "Z"], pch = ".")
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)
###################################################
### code chunk number 8: figode
###################################################
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")
plot(out[, "X"], out[, "Z"], pch = ".")
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)
###################################################
### code chunk number 9: deSolve.Rnw:315-318
###################################################
outb <- radau(state, times, Lorenz, parameters, atol = 1e-4, rtol = 1e-4)
outc <- ode(state, times, Lorenz, parameters, method = "radau",
atol = 1e-4, rtol = 1e-4)
###################################################
### code chunk number 10: deSolve.Rnw:334-340
###################################################
print(system.time(out1 <- rk4 (state, times, Lorenz, parameters)))
print(system.time(out2 <- lsode (state, times, Lorenz, parameters)))
print(system.time(out <- lsoda (state, times, Lorenz, parameters)))
print(system.time(out <- lsodes(state, times, Lorenz, parameters)))
print(system.time(out <- daspk (state, times, Lorenz, parameters)))
print(system.time(out <- vode (state, times, Lorenz, parameters)))
###################################################
### code chunk number 11: deSolve.Rnw:358-359
###################################################
rkMethod()
###################################################
### code chunk number 12: deSolve.Rnw:368-369
###################################################
rkMethod("rk23")
###################################################
### code chunk number 13: deSolve.Rnw:382-403
###################################################
func <- function(t, x, parms) {
with(as.list(c(parms, x)),{
dP <- a * P - b * C * P
dC <- b * P * C - c * C
res <- c(dP, dC)
list(res)
})
}
rKnew <- rkMethod(ID = "midpoint",
varstep = FALSE,
A = c(0, 1/2),
b1 = c(0, 1),
c = c(0, 1/2),
stage = 2,
Qerr = 1
)
out <- ode(y = c(P = 2, C = 1), times = 0:100, func,
parms = c(a = 0.1, b = 0.1, c = 0.1), method = rKnew)
head(out)
###################################################
### code chunk number 14: deSolve.Rnw:437-439
###################################################
diagnostics(out1)
diagnostics(out2)
###################################################
### code chunk number 15: deSolve.Rnw:443-444
###################################################
summary(out1)
###################################################
### code chunk number 16: deSolve.Rnw:518-526
###################################################
Aphid <- function(t, APHIDS, parameters) {
deltax <- c (0.5, rep(1, numboxes - 1), 0.5)
Flux <- -D * diff(c(0, APHIDS, 0)) / deltax
dAPHIDS <- -diff(Flux) / delx + APHIDS * r
# the return value
list(dAPHIDS )
} # end
###################################################
### code chunk number 17: deSolve.Rnw:531-538
###################################################
D <- 0.3 # m2/day diffusion rate
r <- 0.01 # /day net growth rate
delx <- 1 # m thickness of boxes
numboxes <- 60
# distance of boxes on plant, m, 1 m intervals
Distance <- seq(from = 0.5, by = delx, length.out = numboxes)
###################################################
### code chunk number 18: deSolve.Rnw:543-547
###################################################
# Initial conditions: # ind/m2
APHIDS <- rep(0, times = numboxes)
APHIDS[30:31] <- 1
state <- c(APHIDS = APHIDS) # initialise state variables
###################################################
### code chunk number 19: deSolve.Rnw:554-558
###################################################
times <-seq(0, 200, by = 1)
print(system.time(
out <- ode.1D(state, times, Aphid, parms = 0, nspec = 1, names = "Aphid")
))
###################################################
### code chunk number 20: deSolve.Rnw:564-565
###################################################
head(out[,1:5])
###################################################
### code chunk number 21: deSolve.Rnw:569-570
###################################################
summary(out)
###################################################
### code chunk number 22: deSolve.Rnw:603-605
###################################################
data <- cbind(dist = c(0,10, 20, 30, 40, 50, 60),
Aphid = c(0,0.1,0.25,0.5,0.25,0.1,0))
###################################################
### code chunk number 23: matplot1d
###################################################
par (mfrow = c(1,2))
matplot.1D(out, grid = Distance, type = "l", mfrow = NULL,
subset = time %in% seq(0, 200, by = 10),
obs = data, obspar = list(pch = 18, cex = 2, col="red"))
plot.1D(out, grid = Distance, type = "l", mfrow = NULL,
subset = time == 100,
obs = data, obspar = list(pch = 18, cex = 2, col="red"))
###################################################
### code chunk number 24: matplot1d
###################################################
par (mfrow = c(1,2))
matplot.1D(out, grid = Distance, type = "l", mfrow = NULL,
subset = time %in% seq(0, 200, by = 10),
obs = data, obspar = list(pch = 18, cex = 2, col="red"))
plot.1D(out, grid = Distance, type = "l", mfrow = NULL,
subset = time == 100,
obs = data, obspar = list(pch = 18, cex = 2, col="red"))
###################################################
### code chunk number 25: deSolve.Rnw:671-686
###################################################
daefun <- function(t, y, dy, parameters) {
res1 <- dy[1] + y[1] - y[2]
res2 <- y[2] * y[1] - t
list(c(res1, res2))
}
library(deSolve)
yini <- c(1, 0)
dyini <- c(1, 0)
times <- seq(0, 10, 0.1)
## solver
system.time(out <- daspk(y = yini, dy = dyini,
times = times, res = daefun, parms = 0))
###################################################
### code chunk number 26: dae
###################################################
matplot(out[,1], out[,2:3], type = "l", lwd = 2,
main = "dae", xlab = "time", ylab = "y")
###################################################
### code chunk number 27: figdae
###################################################
matplot(out[,1], out[,2:3], type = "l", lwd = 2,
main = "dae", xlab = "time", ylab = "y")
###################################################
### code chunk number 28: deSolve.Rnw:719-729
###################################################
pendulum <- function (t, Y, parms) {
with (as.list(Y),
list(c(u,
v,
-lam * x,
-lam * y - 9.8,
x^2 + y^2 -1
))
)
}
###################################################
### code chunk number 29: deSolve.Rnw:732-733
###################################################
yini <- c(x = 1, y = 0, u = 0, v = 1, lam = 1)
###################################################
### code chunk number 30: deSolve.Rnw:736-739
###################################################
M <- diag(nrow = 5)
M[5, 5] <- 0
M
###################################################
### code chunk number 31: deSolve.Rnw:743-747
###################################################
index <- c(2, 2, 1)
times <- seq(from = 0, to = 10, by = 0.01)
out <- radau (y = yini, func = pendulum, parms = NULL,
times = times, mass = M, nind = index)
###################################################
### code chunk number 32: pendulum
###################################################
plot(out, type = "l", lwd = 2)
plot(out[, c("x", "y")], type = "l", lwd = 2)
###################################################
### code chunk number 33: pendulum
###################################################
plot(out, type = "l", lwd = 2)
plot(out[, c("x", "y")], type = "l", lwd = 2)
###################################################
### code chunk number 34: deSolve.Rnw:781-794
###################################################
ZODE2 <- function(Time, State, Pars) {
with(as.list(State), {
df <- 1i * f
dg <- -1i * g * g * f
return(list(c(df, dg)))
})
}
yini <- c(f = 1+0i, g = 1/2.1+0i)
times <- seq(0, 2 * pi, length = 100)
out <- zvode(func = ZODE2, y = yini, parms = NULL, times = times,
atol = 1e-10, rtol = 1e-10)
###################################################
### code chunk number 35: deSolve.Rnw:806-808
###################################################
analytical <- cbind(f = exp(1i*times), g = 1/(exp(1i*times)+1.1))
tail(cbind(out[,2], analytical[,1]))
###################################################
### code chunk number 36: deSolve.Rnw:821-832
###################################################
f1 <- function (t, y, parms) {
ydot <- vector(len = 5)
ydot[1] <- 0.1*y[1] -0.2*y[2]
ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3]
ydot[3] <- -0.3*y[2] +0.1*y[3] -0.2*y[4]
ydot[4] <- -0.3*y[3] +0.1*y[4] -0.2*y[5]
ydot[5] <- -0.3*y[4] +0.1*y[5]
return(list(ydot))
}
###################################################
### code chunk number 37: deSolve.Rnw:837-839
###################################################
yini <- 1:5
times <- 1:20
###################################################
### code chunk number 38: deSolve.Rnw:846-847
###################################################
out <- lsode(yini, times, f1, parms = 0, jactype = "fullint")
###################################################
### code chunk number 39: deSolve.Rnw:854-863
###################################################
fulljac <- function (t, y, parms) {
jac <- matrix(nrow = 5, ncol = 5, byrow = TRUE,
data = c(0.1, -0.2, 0 , 0 , 0 ,
-0.3, 0.1, -0.2, 0 , 0 ,
0 , -0.3, 0.1, -0.2, 0 ,
0 , 0 , -0.3, 0.1, -0.2,
0 , 0 , 0 , -0.3, 0.1))
return(jac)
}
###################################################
### code chunk number 40: deSolve.Rnw:868-870
###################################################
out2 <- lsode(yini, times, f1, parms = 0, jactype = "fullusr",
jacfunc = fulljac)
###################################################
### code chunk number 41: deSolve.Rnw:877-879
###################################################
out3 <- lsode(yini, times, f1, parms = 0, jactype = "bandint",
bandup = 1, banddown = 1)
###################################################
### code chunk number 42: deSolve.Rnw:884-891
###################################################
bandjac <- function (t, y, parms) {
jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE,
data = c( 0 , -0.2, -0.2, -0.2, -0.2,
0.1, 0.1, 0.1, 0.1, 0.1,
-0.3, -0.3, -0.3, -0.3, 0))
return(jac)
}
###################################################
### code chunk number 43: deSolve.Rnw:896-898
###################################################
out4 <- lsode(yini, times, f1, parms = 0, jactype = "bandusr",
jacfunc = bandjac, bandup = 1, banddown = 1)
###################################################
### code chunk number 44: deSolve.Rnw:904-905
###################################################
out5 <- lsode(yini, times, f1, parms = 0, mf = 10)
###################################################
### code chunk number 45: deSolve.Rnw:936-942
###################################################
eventmod <- function(t, var, parms) {
list(dvar = -0.1*var)
}
yini <- c(v1 = 1, v2 = 2)
times <- seq(0, 10, by = 0.1)
###################################################
### code chunk number 46: deSolve.Rnw:949-953
###################################################
eventdat <- data.frame(var = c("v1", "v2", "v2", "v1"), time = c(1, 1, 5, 9),
value = c(1, 2, 3, 4), method = c("add", "mult", "rep", "add"))
eventdat
###################################################
### code chunk number 47: deSolve.Rnw:958-960
###################################################
out <- ode(func = eventmod, y = yini, times = times, parms = NULL,
events = list(data = eventdat))
###################################################
### code chunk number 48: event1
###################################################
plot(out, type = "l", lwd = 2)
###################################################
### code chunk number 49: figevent1
###################################################
plot(out, type = "l", lwd = 2)
###################################################
### code chunk number 50: deSolve.Rnw:982-987
###################################################
ballode<- function(t, y, parms) {
dy1 <- y[2]
dy2 <- -9.8
list(c(dy1, dy2))
}
###################################################
### code chunk number 51: deSolve.Rnw:994-995
###################################################
root <- function(t, y, parms) y[1]
###################################################
### code chunk number 52: deSolve.Rnw:1000-1005
###################################################
event <- function(t, y, parms) {
y[1]<- 0
y[2]<- -0.9 * y[2]
return(y)
}
###################################################
### code chunk number 53: deSolve.Rnw:1011-1016
###################################################
yini <- c(height = 0, v = 20)
times <- seq(from = 0, to = 20, by = 0.01)
out <- lsode(times = times, y = yini, func = ballode, parms = NULL,
events = list(func = event, root = TRUE), rootfun = root)
###################################################
### code chunk number 54: event2
###################################################
plot(out, which = "height", type = "l",lwd = 2,
main = "bouncing ball", ylab = "height")
###################################################
### code chunk number 55: figevent2
###################################################
plot(out, which = "height", type = "l",lwd = 2,
main = "bouncing ball", ylab = "height")
###################################################
### code chunk number 56: deSolve.Rnw:1065-1067
###################################################
times <- seq(0, 1, 0.1)
eventtimes <- c(0.7, 0.9)
###################################################
### code chunk number 57: deSolve.Rnw:1072-1073
###################################################
eventtimes %in% times
###################################################
### code chunk number 58: deSolve.Rnw:1080-1082
###################################################
times2 <- round(times, 1)
times - times2
###################################################
### code chunk number 59: deSolve.Rnw:1093-1094
###################################################
eventtimes %in% times2
###################################################
### code chunk number 60: deSolve.Rnw:1099-1100
###################################################
all(eventtimes %in% times2)
###################################################
### code chunk number 61: deSolve.Rnw:1110-1113
###################################################
times <- 1:10
eventtimes <- c(1.3, 3.4, 4, 7.9, 8.5)
newtimes <- sort(unique(c(times, eventtimes)))
###################################################
### code chunk number 62: deSolve.Rnw:1119-1122
###################################################
times <- 1:10
eventtimes <- c(1.3, 3.4, 4, 7.9999999999999999, 8.5)
newtimes <- sort(c(eventtimes, cleanEventTimes(times, eventtimes)))
###################################################
### code chunk number 63: deSolve.Rnw:1151-1185
###################################################
library(deSolve)
#-----------------------------
# the derivative function
#-----------------------------
derivs <- function(t, y, parms) {
if (t < 0)
lag <- 19
else
lag <- lagvalue(t - 0.74)
dy <- r * y * (1 - lag/m)
list(dy, dy = dy)
}
#-----------------------------
# parameters
#-----------------------------
r <- 3.5; m <- 19
#-----------------------------
# initial values and times
#-----------------------------
yinit <- c(y = 19.001)
times <- seq(-0.74, 40, by = 0.01)
#-----------------------------
# solve the model
#-----------------------------
yout <- dede(y = yinit, times = times, func = derivs,
parms = NULL, atol = 1e-10)
###################################################
### code chunk number 64: dde
###################################################
plot(yout, which = 1, type = "l", lwd = 2,
main = "Lemming model", mfrow = c(1,2))
plot(yout[,2], yout[,3], xlab = "y", ylab = "dy", type = "l", lwd = 2)
###################################################
### code chunk number 65: figdde
###################################################
plot(yout, which = 1, type = "l", lwd = 2,
main = "Lemming model", mfrow = c(1,2))
plot(yout[,2], yout[,3], xlab = "y", ylab = "dy", type = "l", lwd = 2)
###################################################
### code chunk number 66: deSolve.Rnw:1222-1234
###################################################
Stages <- c("DS 1yr", "DS 2yr", "R small", "R medium", "R large", "F")
NumStages <- length(Stages)
# Population matrix
A <- matrix(nrow = NumStages, ncol = NumStages, byrow = TRUE, data = c(
0, 0, 0, 0, 0, 322.38,
0.966, 0, 0, 0, 0, 0 ,
0.013, 0.01, 0.125, 0, 0, 3.448 ,
0.007, 0, 0.125, 0.238, 0, 30.170,
0.008, 0, 0.038, 0.245, 0.167, 0.862 ,
0, 0, 0, 0.023, 0.75, 0 ) )
###################################################
### code chunk number 67: deSolve.Rnw:1239-1243
###################################################
Teasel <- function (t, y, p) {
yNew <- A %*% y
list (yNew / sum(yNew))
}
###################################################
### code chunk number 68: deSolve.Rnw:1246-1248
###################################################
out <- ode(func = Teasel, y = c(1, rep(0, 5) ), times = 0:50,
parms = 0, method = "iteration")
###################################################
### code chunk number 69: difference
###################################################
matplot(out[,1], out[,-1], main = "Teasel stage distribution", type = "l")
legend("topright", legend = Stages, lty = 1:6, col = 1:6)
###################################################
### code chunk number 70: difference
###################################################
matplot(out[,1], out[,-1], main = "Teasel stage distribution", type = "l")
legend("topright", legend = Stages, lty = 1:6, col = 1:6)
###################################################
### code chunk number 71: deSolve.Rnw:1291-1295
###################################################
library(deSolve)
combustion <- function (t, y, parms)
list(y^2 * (1-y) )
###################################################
### code chunk number 72: deSolve.Rnw:1297-1299
###################################################
yini <- 0.01
times <- 0 : 200
###################################################
### code chunk number 73: deSolve.Rnw:1301-1305
###################################################
out <- ode(times = times, y = yini, parms = 0, func = combustion)
out2 <- ode(times = times, y = yini*2, parms = 0, func = combustion)
out3 <- ode(times = times, y = yini*3, parms = 0, func = combustion)
out4 <- ode(times = times, y = yini*4, parms = 0, func = combustion)
###################################################
### code chunk number 74: plotdeSolve
###################################################
plot(out, out2, out3, out4, main = "combustion")
legend("bottomright", lty = 1:4, col = 1:4, legend = 1:4, title = "yini*i")
###################################################
### code chunk number 75: plotdeSolve
###################################################
plot(out, out2, out3, out4, main = "combustion")
legend("bottomright", lty = 1:4, col = 1:4, legend = 1:4, title = "yini*i")
###################################################
### code chunk number 76: deSolve.Rnw:1335-1336
###################################################
head(ccl4data)
###################################################
### code chunk number 77: deSolve.Rnw:1339-1342
###################################################
obs <- subset (ccl4data, animal == "A", c(time, ChamberConc))
names(obs) <- c("time", "CP")
head(obs)
###################################################
### code chunk number 78: deSolve.Rnw:1348-1362
###################################################
parms <- c(0.182, 4.0, 4.0, 0.08, 0.04, 0.74, 0.05, 0.15, 0.32, 16.17,
281.48, 13.3, 16.17, 5.487, 153.8, 0.04321671,
0.40272550, 951.46, 0.02, 1.0, 3.80000000)
yini <- c(AI = 21, AAM = 0, AT = 0, AF = 0, AL = 0, CLT = 0, AM = 0)
out <- ccl4model(times = seq(0, 6, by = 0.05), y = yini, parms = parms)
par2 <- parms
par2[1] <- 0.1
out2 <- ccl4model(times = seq(0, 6, by = 0.05), y = yini, parms = par2)
par3 <- parms
par3[1] <- 0.05
out3 <- ccl4model(times = seq(0, 6, by = 0.05), y = yini, parms = par3)
###################################################
### code chunk number 79: plotobs
###################################################
plot(out, out2, out3, which = c("AI", "MASS", "CP"),
col = c("black", "red", "green"), lwd = 2,
obs = obs, obspar = list(pch = 18, col = "blue", cex = 1.2))
legend("topright", lty = c(1,2,3,NA), pch = c(NA, NA, NA, 18),
col = c("black", "red", "green", "blue"), lwd = 2,
legend = c("par1", "par2", "par3", "obs"))
###################################################
### code chunk number 80: plotobs
###################################################
plot(out, out2, out3, which = c("AI", "MASS", "CP"),
col = c("black", "red", "green"), lwd = 2,
obs = obs, obspar = list(pch = 18, col = "blue", cex = 1.2))
legend("topright", lty = c(1,2,3,NA), pch = c(NA, NA, NA, 18),
col = c("black", "red", "green", "blue"), lwd = 2,
legend = c("par1", "par2", "par3", "obs"))
###################################################
### code chunk number 81: deSolve.Rnw:1388-1390
###################################################
obs2 <- data.frame(time = 6, MASS = 12)
obs2
###################################################
### code chunk number 82: obs2
###################################################
plot(out, out2, out3, lwd = 2,
obs = list(obs, obs2),
obspar = list(pch = c(16, 18), col = c("blue", "black"),
cex = c(1.2 , 2))
)
###################################################
### code chunk number 83: plotobs2
###################################################
plot(out, out2, out3, lwd = 2,
obs = list(obs, obs2),
obspar = list(pch = c(16, 18), col = c("blue", "black"),
cex = c(1.2 , 2))
)
###################################################
### code chunk number 84: hist
###################################################
hist(out, col = grey(seq(0, 1, by = 0.1)), mfrow = c(3, 4))
###################################################
### code chunk number 85: plothist
###################################################
hist(out, col = grey(seq(0, 1, by = 0.1)), mfrow = c(3, 4))
###################################################
### code chunk number 86: deSolve.Rnw:1449-1451
###################################################
options(prompt = " ")
options(continue = " ")
###################################################
### code chunk number 87: deSolve.Rnw:1454-1478
###################################################
lvmod <- function (time, state, parms, N, rr, ri, dr, dri) {
with (as.list(parms), {
PREY <- state[1:N]
PRED <- state[(N+1):(2*N)]
## Fluxes due to diffusion
## at internal and external boundaries: zero gradient
FluxPrey <- -Da * diff(c(PREY[1], PREY, PREY[N]))/dri
FluxPred <- -Da * diff(c(PRED[1], PRED, PRED[N]))/dri
## Biology: Lotka-Volterra model
Ingestion <- rIng * PREY * PRED
GrowthPrey <- rGrow * PREY * (1-PREY/cap)
MortPredator <- rMort * PRED
## Rate of change = Flux gradient + Biology
dPREY <- -diff(ri * FluxPrey)/rr/dr +
GrowthPrey - Ingestion
dPRED <- -diff(ri * FluxPred)/rr/dr +
Ingestion * assEff - MortPredator
return (list(c(dPREY, dPRED)))
})
}
###################################################
### code chunk number 88: deSolve.Rnw:1480-1482
###################################################
options(prompt = " ")
options(continue = " ")
###################################################
### code chunk number 89: deSolve.Rnw:1485-1499
###################################################
R <- 20 # total radius of surface, m
N <- 100 # 100 concentric circles
dr <- R/N # thickness of each layer
r <- seq(dr/2,by = dr,len = N) # distance of center to mid-layer
ri <- seq(0,by = dr,len = N+1) # distance to layer interface
dri <- dr # dispersion distances
parms <- c(Da = 0.05, # m2/d, dispersion coefficient
rIng = 0.2, # /day, rate of ingestion
rGrow = 1.0, # /day, growth rate of prey
rMort = 0.2 , # /day, mortality rate of pred
assEff = 0.5, # -, assimilation efficiency
cap = 10) # density, carrying capacity
###################################################
### code chunk number 90: deSolve.Rnw:1502-1512
###################################################
state <- rep(0, 2 * N)
state[1] <- state[N + 1] <- 10
times <- seq(0, 200, by = 1) # output wanted at these time intervals
print(system.time(
out <- ode.1D(y = state, times = times, func = lvmod, parms = parms,
nspec = 2, names = c("PREY", "PRED"),
N = N, rr = r, ri = ri, dr = dr, dri = dri)
))
###################################################
### code chunk number 91: deSolve.Rnw:1515-1516
###################################################
summary(out)
###################################################
### code chunk number 92: deSolve.Rnw:1520-1522
###################################################
p10 <- subset(out, select = "PREY", subset = time == 10)
head(p10, n = 5)
###################################################
### code chunk number 93: deSolve.Rnw:1568-1573
###################################################
Simple2D <- function(t, Y, par) {
y <- matrix(nrow = nx, ncol = ny, data = Y) # vector to 2-D matrix
dY <- - r_x2y2 * y # consumption
return(list(dY))
}
###################################################
### code chunk number 94: deSolve.Rnw:1577-1584
###################################################
dy <- dx <- 1 # grid size
nx <- ny <- 100
x <- seq (dx/2, by = dx, len = nx)
y <- seq (dy/2, by = dy, len = ny)
# in each grid cell: consumption depending on position
r_x2y2 <- outer(x, y, FUN=function(x,y) ((x-50)^2 + (y-50)^2)*1e-4)
###################################################
### code chunk number 95: deSolve.Rnw:1588-1591
###################################################
C <- matrix(nrow = nx, ncol = ny, 1)
ODE3 <- ode.2D(y = C, times = 1:100, func = Simple2D, parms = NULL,
dimens = c(nx, ny), names = "C", method = "ode45")
###################################################
### code chunk number 96: deSolve.Rnw:1594-1597
###################################################
summary(ODE3)
t50 <- matrix(nrow = nx, ncol = ny,
data = subset(ODE3, select = "C", subset = (time == 50)))
###################################################
### code chunk number 97: twoD
###################################################
par(mfrow = c(1, 2))
contour(x, y, r_x2y2, main = "consumption")
contour(x, y, t50, main = "Y(t = 50)")
###################################################
### code chunk number 98: twoD
###################################################
par(mfrow = c(1, 2))
contour(x, y, r_x2y2, main = "consumption")
contour(x, y, t50, main = "Y(t = 50)")
###################################################
### code chunk number 99: deSolve.Rnw:1627-1635
###################################################
PCmod <- function(t, x, parms) {
with(as.list(c(parms, x)), {
dP <- c*P - d*C*P # producer
dC <- e*P*C - f*C # consumer
res <- c(dP, dC)
list(res)
})
}
###################################################
### code chunk number 100: deSolve.Rnw:1642-1643
###################################################
parms <- c(c = 10, d = 0.1, e = 0.1, f = 0.1)
###################################################
### code chunk number 101: deSolve.Rnw:1649-1655
###################################################
xstart <- c(P = 0.5, C = 1)
times <- seq(0, 200, 0.1)
out <- ode(y = xstart, times = times,
func = PCmod, parms = parms)
tail(out)
###################################################
### code chunk number 102: deSolve.Rnw:1676-1680
###################################################
out <- ode(y = xstart,times = times, func = PCmod,
parms = parms, atol = 0)
matplot(out[,1], out[,2:3], type = "l",
xlab = "time", ylab = "Producer, Consumer")
###################################################
### code chunk number 103: deSolve.Rnw:1736-1760
###################################################
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(func = LVmod, y = yini,
parms = pars, times = times)
###################################################
### code chunk number 104: deSolve.Rnw:1772-1775
###################################################
pars["rIng"] <- 100
out2 <- ode(func = LVmod, y = yini,
parms = pars, times = times)
###################################################
### code chunk number 105: err
###################################################
plot(out2, type = "l", lwd = 2, main = "corrupt Lotka-Volterra model")
###################################################
### code chunk number 106: deSolve.Rnw:1824-1827
###################################################
pars["rIng"] <- 100
out3 <- ode(func = LVmod, y = yini, parms = pars,
times = times, method = "ode45", atol = 1e-14, rtol = 1e-14)
###################################################
### code chunk number 107: err
###################################################
plot(out2, type = "l", lwd = 2, main = "corrupt Lotka-Volterra model")
Try the deSolve package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.