library(Revolve)
## Figures from Huisman 2001:
plot_huisman <- function(obj, ylim=NULL, lty=1:n, col=1:n) {
np <- length(obj)
n <- ncol(obj[[1]]$y)
if (is.null(ylim)) {
ylim <- range(sapply(obj, function(x) range(x$y)))
}
op <- par(mfrow=c(np, 1), mar=rep(.5, 4), oma=c(3, 4, 0, 0), cex=1)
on.exit(par(op))
for (i in seq_len(np)) {
matplot(obj[[i]]$t, obj[[i]]$y, type="l", lty=lty, col=col,
ylim=ylim, las=1, xaxt="n", xlab="", ylab="")
axis(1, labels=i == np)
}
mtext("Time (d)", 1, line=2.1, xpd=NA)
mtext("Population abundance", 2, line=2, outer=TRUE)
}
# First, results from the paper. These are the matrices given in
# Appendix C.
# Figure2:
K.fig2 <- rbind(c(1.0, 0.6, 0.3),
c(0.3, 1.0, 0.6),
c(0.6, 0.3, 1.0))
C.fig2a <- rbind(c(0.07, 0.04, 0.04),
c(0.08, 0.10, 0.08),
c(0.10, 0.10, 0.14))
C.fig2b <- rbind(c(0.04, 0.07, 0.04),
c(0.08, 0.08, 0.10),
c(0.14, 0.10, 0.10))
C.fig2c <- rbind(c(0.04, 0.04, 0.07),
c(0.10, 0.08, 0.08),
c(0.10, 0.14, 0.10))
mat.fig2a <- rstar_matrices_fixed(K.fig2, C.fig2a)
mat.fig2b <- rstar_matrices_fixed(K.fig2, C.fig2b)
mat.fig2c <- rstar_matrices_fixed(K.fig2, C.fig2c)
sys.fig2 <- sys(NULL, rep(1, 3))
S.fig2 <- c(6, 10, 14)
t.fig2 <- seq(0, 100, length=201)
m.fig2 <- list(rstar(mat.fig2a, S=S.fig2),
rstar(mat.fig2b, S=S.fig2),
rstar(mat.fig2c, S=S.fig2))
res.fig2 <- lapply(m.fig2, function(m) m$run(sys.fig2, t.fig2))
plot_huisman(res.fig2, ylim=c(0, 100), lty=c(1, 2, 4))
# Figure 3
K.fig3a <- rbind(c(1.0, 0.9, 0.3),
c(0.3, 1.0, 0.9),
c(0.9, 0.3, 1.0))
K.fig3b <- K.fig3a
K.fig3b[rbind(c(1,2), c(2,3), c(3,1))] <- 0.5
K.fig3c <- K.fig3a
K.fig3c[rbind(c(1,2), c(2,3), c(3,1))] <- 0.4
C.fig3 <- rbind(c(0.04, 0.07, 0.04),
c(0.08, 0.08, 0.10),
c(0.14, 0.10, 0.10))
sys.fig3 <- sys(NULL, rep(1, 3))
S.fig3 <- c(6, 10, 14)
t.fig3 <- seq(0, 400, length=201)
mat.fig3a <- rstar_matrices_fixed(K.fig3a, C.fig3)
mat.fig3b <- rstar_matrices_fixed(K.fig3b, C.fig3)
mat.fig3c <- rstar_matrices_fixed(K.fig3c, C.fig3)
m.fig3 <- list(rstar(mat.fig3a, S=S.fig3),
rstar(mat.fig3b, S=S.fig3),
rstar(mat.fig3c, S=S.fig3))
res.fig3 <- lapply(m.fig3, function(m) m$run(sys.fig3, t.fig3))
# Note that this is different to the paper, because they run the last
# case out to to 2000, not 400. But close enough for now.
plot_huisman(res.fig3, ylim=c(0, 150), lty=c(1, 2, 4))
# Figure 4:
K.fig4a <- rbind(c(1.0, 0.8, 0.4, 0.3, 0.2),
c(0.2, 1.0, 0.8, 0.4, 0.3),
c(0.3, 0.2, 1.0, 0.8, 0.4),
c(0.4, 0.3, 0.2, 1.0, 0.8),
c(0.8, 0.4, 0.3, 0.2, 1.0))
C.fig4 <- rbind(c(0.04, 0.07, 0.04, 0.04, 0.04),
c(0.08, 0.08, 0.10, 0.08, 0.08),
c(0.10, 0.10, 0.10, 0.14, 0.10),
c(0.03, 0.03, 0.03, 0.03, 0.05),
c(0.09, 0.07, 0.07, 0.07, 0.07))
K.fig4b <- K.fig4a
K.fig4b[rbind(c(1,2), c(2,3), c(3,4), c(4,5), c(5,1))] <- 0.6
sys.fig4 <- sys(NULL, rep(1, 5))
S.fig4 <- c(6, 10, 14, 4, 9)
t.fig4a <- seq(0, 500, length=201)
t.fig4b <- seq(0, 3000, length=201)
mat.fig4a <- rstar_matrices_fixed(K.fig4a, C.fig4)
mat.fig4b <- rstar_matrices_fixed(K.fig4b, C.fig4)
m.fig4 <- list(rstar(mat.fig4a, S=S.fig4),
rstar(mat.fig4b, S=S.fig4))
res.fig4 <- list(m.fig4[[1]]$run(sys.fig4, t.fig4a),
m.fig4[[2]]$run(sys.fig4, t.fig4b))
# Looks OK except that we start overstepping at some point and drive
# numbers negative, but apart from that we're all good really.
plot_huisman(res.fig4, ylim=c(0, 100))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.