### 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")