################################################### ### chunk number 1: ################################################### setwd("~/Documents/Projects/Book/draft/Chap08") rm(list=ls()) badfonts = FALSE .detach.stuff = function() { s1 = grep("package|Autoloads",search()) nattach = length(search()) xattach = search()[c(-1,-s1)] for (i in xattach) eval(substitute(detach(i),list(i=i))) } .detach.stuff() ## ps.options(family="NimbusSan") ps.options(colormodel="cmyk") options(width=75, SweaveHooks=list(fig=function() par(mar=c(5,4,.5,1.5))) ) ## palette(gray((0:8)/8)) library(primer) trellis.par.set(canonical.theme(color=FALSE)) #source("DMBpplane.R") require(deSolve); ################################################### ### chunk number 2: valleymba ################################################### z <- 2 * volcano x <- 10 * (1:nrow(z)) y <- 10 * (1:ncol(z)) par(mar=c(0,0,0,0)) persp(x, y, -z, theta = 120, phi = 30, col="green",shade=.75, scale = FALSE, axes = FALSE) ################################################### ### chunk number 3: contourmba ################################################### par(mar=c(0.5,0.5,0.5,0.5)) contour(x,y,z) ################################################### ### chunk number 4: ################################################### lvcompg <- function(t, n, parms) { r <- parms[[1]]; a <- parms[[2]] dns.dt <- r * n * (1 - (a%*%n)) return( list(c(dns.dt)) ) } ################################################### ### chunk number 5: ################################################### r <- c(r1=.6, r2=1, r3=2) a <- matrix(c(a11=.001, a12=.002, a13=.002, a21=.002, a22=.00101, a23=.002, a31=.002, a32=.002, a33=.00102), nrow=3, ncol=3) parms <- list(r,a) parms ################################################### ### chunk number 6: ################################################### t=seq(0,40, by=.1) ni <- 200; std=10 N0 <- sapply(1:30, function(i) rnorm(3, mean=ni,sd=std) ) ################################################### ### chunk number 7: ################################################### N0[,1] <- ni ################################################### ### chunk number 8: LVPriority ################################################### par(mar=c(2,2,1,1)) layout(matrix(1:30, ncol=5)) for(i in 1:30) {lvout <- ode(N0[,i], t, lvcompg, parms) matplot(t, lvout[,2:4], type="l", lwd=1.5, col=1) if(all(N0[,i]==200)) { text(3, 500, "Equal", srt=90) }else { text(3, 500, paste("Sp.", which.max(N0[,i])), srt=90)} lastN <- lvout[nrow(lvout),2:4] text(3, max(lastN), paste("Sp.", which.max(lastN)), adj=c(0,1)) } ################################################### ### chunk number 9: Hysteresis ################################################### par(mar=c(3,3,1,1), mgp=c(1,0,0), tck=0) x <- seq(-2, 5, by=0.01) y.ex <- expression(x^3 - 4*x^2) y <- eval(y.ex) plot(y, -x, type='n', axes=F, ylab="State Variable\n(e.g. Precipitation)", xlab="Environmental Variable\n(e.g. Temperature)") polygon(c(-10, 0, 0, -10), c(-4.5, -4.5, 2, 2), col="grey", border=NA) lines(y, -x) lines(y[x>=0 & x <= 8/3], -x[x>=0 & x <= 8/3], col=0, lty=2) box() xs <- x[c(1, 60, 701, 660)] x <- xs ys <- eval(y.ex) arrows(ys[1], -xs[1]-.5, ys[2], -xs[2]-.5, length=.1) arrows(ys[3], -xs[3]+.5, ys[4], -xs[4]+.5, length=.1) #segments(-9, -4.5, -1, -4.5, lwd=10, col="grey") ################################################### ### chunk number 10: ################################################### p <- c(N=1, as=0.01, af=0.01, b=0.02, qs=0.075, qf=0.005, hs=0, hf=0.2, ls=0.05, lf=0.05, rs=0.5, rf=0.5, W=0) t <- 1:200 Initial <- c(F=10, S=10) ################################################### ### chunk number 11: Scheffer1 ################################################### S.out1 <- ode(Initial, t, scheffer, p) matplot(t, S.out1[,-1], type='l') legend('right', c("F", "S"), lty=1:2, bty='n') ################################################### ### chunk number 12: Scheffer2 ################################################### p['N'] <- 4 S.out2 <- ode(Initial, t, scheffer, p) matplot(t, S.out2[,-1], type='l') ################################################### ### chunk number 13: ################################################### N.s <- seq(.5,4, by=0.1) t <- 1:1000 S.s <- t( sapply(N.s, function(x) {p['N'] <- x ode(Initial, t, scheffer, p)[length(t),2:3] }) ) ################################################### ### chunk number 14: Scheffer3 ################################################### matplot(N.s, S.s, type='l') legend('right', c('F','S'), lty=1:2, bty='n') ################################################### ### chunk number 15: ################################################### Initial.Eutrophic <- c(F=600, S=10) S.s.E <- t( sapply(N.s, function(x) {p['N'] <- x ode(Initial.Eutrophic, c(1,1000), scheffer, p)[2,2:3] }) ) ################################################### ### chunk number 16: Scheffer4 ################################################### matplot(N.s, S.s.E, type='l') ################################################### ### chunk number 17: ################################################### quartz(,4,4) ################################################### ### chunk number 18: ################################################### plot(N.s[1:23], S.s[1:23,1], type='l', lwd=2, xlim=c(0,4), ylim=c(0,900), main="Floating Plants", ylab=expression("Biomass (g m"^-2*")"), xlab="Nitrogen Supply Rate") lines(N.s[-(1:5)], S.s.E[-(1:5),1], lwd=2) ################################################### ### chunk number 19: ################################################### arrows(3, 10, 3, 620, length=.1); arrows(3, 820, 3, 720, length=.1) arrows(0.5, 620, 0.5, 50, length=.1) ################################################### ### chunk number 20: ################################################### arrows(2.5, -10, 2.5, 60, length=.1); arrows(2.5, 200, 2.5, 100, length=.1) text(2.5, 100, "Coexisting\nwith S", adj=c(1.1, 0)) ################################################### ### chunk number 21: ################################################### arrows(2, 480, 2, 580, length=.1); arrows(2, 750, 2, 650, length=.1) text(2, 700, "Monoculture", adj=c(1.1,0) ) ################################################### ### chunk number 22: ################################################### dev.print(pdf, "Scheffer5.pdf") ################################################### ### chunk number 23: ################################################### quartz(,4,4) ################################################### ### chunk number 24: ################################################### plot(N.s[1:23], S.s[1:23,2], type='l', lwd=2, xlim=c(0,4), ylim=c(0,900), main="Submerged Plants", ylab=expression("Biomass (g m"^-2*")"), xlab="Nitrogen Supply Rate") lines(N.s[-(1:5)], S.s.E[-(1:5),2], lwd=2) ################################################### ### chunk number 25: ################################################### arrows(.7, 30, .7, 830, length=.1) arrows(3.8, 830, 3.8, 30, length=.1) ################################################### ### chunk number 26: ################################################### arrows(2.3, 650, 2.3, 750, length=.1); arrows(2.3, 900, 2.3, 800, length=.1) text(2.4, 900, "Coexisting\nwith F", adj=c(0, 1) ) ################################################### ### chunk number 27: ################################################### arrows(2, 130, 2, 30, length=.1) text(2, 140, "Excluded\nDue to Light Comp.", adj=c(.5,-.3) ) ################################################### ### chunk number 28: ################################################### dev.print(pdf, "Scheffer6.pdf") ################################################### ### chunk number 29: ################################################### igp ################################################### ### chunk number 30: ################################################### params1 <- c(bpb= 0.032, abp=10^-8, bpn=10^-5, anp=10^-4, mp=1, bnb=0.04, abn=10^-8, mn=1, r=1, abb=10^-9.5) ################################################### ### chunk number 31: ################################################### t=seq(0, 60, by=.1) N.init <- cbind(B = rep(10^9, 4), N = 10^c(2, 5, 3, 4), P = 10^c(5,2,3, 4) ) ################################################### ### chunk number 32: ################################################### quartz(,4,4); layout(matrix(1:4, nr=2)); par(mar=c(4, 4, 1, 1)) for(i in 1:4) { igp.out <- ode(N.init[i,1:3], t, igp, params1) matplot(t, log10(igp.out[,2:4]+1), type="l", lwd=2, ylab="log(Abundance)") } ################################################### ### chunk number 33: ################################################### dev.print(pdf, file="IGPDyn.pdf") ################################################### ### chunk number 34: ################################################### logNP <- seq(2, 5, by=.1) N.inits <- cbind(B=rep(10^9, length(logNP)), N=10^logNP, P=10^rev(logNP) ) ################################################### ### chunk number 35: ################################################### t1 <- 1:500 MBAs <- t( sapply(1:nrow(N.inits), function(i) { tmp <- ode(N.inits[i,], t1, igp, params1, hmax=0.1) cbind(tmp[1,3:4],tmp[50,3:4],tmp[500,3:4]) } ) ) colnames(MBAs) <- c('N1','P1','N50', 'P50', 'N500','P500') ################################################### ### chunk number 36: IGPAltRA0 ################################################### matplot(log10(N.inits[,"N"]/N.inits[,"P"]), log10(MBAs[,1:2]+1), type="l", col=1, lty=2:3, lwd=2, ylab="log(Abundance+1)", xlab="log[N/P]") legend("right", c("N", "P"), lty=2:3, col=1, bty="n") ################################################### ### chunk number 37: IGPAltRA50 ################################################### matplot(log10(N.inits[,"N"]/N.inits[,"P"]), log10(MBAs[,3:4]+1), type="l", col=1, lty=2:3, lwd=2, ylab="log(Abundance+1)", xlab="log[N/P]") ################################################### ### chunk number 38: IGPAltRA500 ################################################### matplot(log10(N.inits[,"N"]/N.inits[,"P"]), log10(MBAs[,5:6]+1), type="l", col=1, lty=2:3, lwd=2, ylab="log(Abundance+1)", xlab="log[N/P]") ################################################### ### chunk number 39: ################################################### logAbs <- seq(2, 7, by=.2) N.abs.inits <- cbind(B=rep(10^9, length(logAbs)), N=10^logAbs, P=10^logAbs ) ################################################### ### chunk number 40: ################################################### t1 <- 1:500 MBA.abs <- t( sapply(1:nrow(N.abs.inits), function(i) { tmp <- ode(N.abs.inits[i,], t1, igp, params1, hmax=0.1) cbind(tmp[1,3:4],tmp[50,3:4],tmp[500,3:4]) } ) ) colnames(MBAs) <- c('N1','P1','N50', 'P50', 'N500','P500') ################################################### ### chunk number 41: ################################################### quartz(,7,3) ################################################### ### chunk number 42: ################################################### layout(matrix(1:3, nr=1) ) matplot(log10(N.abs.inits[,"N"]), log10(MBA.abs[,1:2]+1), type="l", main="Initial Abundances (t=1)", col=2:3, lty=2:3, lwd=2, ylab="log(Abundance+1)", xlab=expression(log[10]("N"))) legend("right", c("N", "P"), lty=2:3, col=2:3, bty="n") matplot(log10(N.abs.inits[,"N"]), log10(MBA.abs[,3:4]+1), type="l", main="At time = 50)", col=2:3, lty=2:3, lwd=2, ylab="log(Abundance+1)", xlab=expression(log[10]("N"))) matplot(log10(N.abs.inits[,"N"]), log10(MBA.abs[,5:6]+1), type="l", main="At time = 500", col=2:3, lty=2:3, lwd=2, ylab="log(Abundance+1)", xlab=expression(log[10]("N"))) ################################################### ### chunk number 43: ################################################### dev.print(pdf, file="IGPAltAbs.pdf") ################################################### ### chunk number 44: ################################################### params1 with(as.list(params1), { rbind(B=c(r-r*abb*10^9, -abn, -abp), N=c( bnb*abn, 0, -anp), P=c(bpb*abp, bpn*anp, 0)) }) with(as.list(params1), { rbind(B=c(r-r*abb*10^9, -abn, -abp), N=c( bnb, 0, -anp), P=c(bpb, bpn, 0)) }) ################################################### ### chunk number 45: ################################################### params2 <- params1 params2['anp'] <- 10^-7 ################################################### ### chunk number 46: IGPtrial1 ################################################### t <- 1:500 N.init.1 <- c(B=10^9, N=10^7, P=1) trial1 <- ode(N.init.1, t, igp, params2) matplot(t, log10(trial1[,-1]+1), type='l', col=1, ylab=quote(log[10]*("Density+1"))) legend('bottomright', c('B','N','P'), lty=1:3, bty='n') ################################################### ### chunk number 47: IGPtrial2 ################################################### params2['bpn'] <- .5 trial2 <- ode(N.init.1, t, igp, params2) matplot(t, log10(trial2[,-1]+1), type='l', col=1 , ylab=quote(log[10]*("Density+1")) )