1種の個体群動態をRで可視化する

個体群増殖モデル
以下の「呪文」は長いので、コピペして利用することをお勧めする。
ただし、コピペすると改行箇所などにRでエラーの原因となるスペースができることがある。そのスペースは削除してから実行してほしい。

環境収容力を想定しないモデルと想定したモデル
#図6.1(a)
R <- 0.2
N0 <- 2000
t <- 0:30
popsize <- N0*(1+R)^t
plot(t, popsize, xlim=c(0, 30), ylim=c(0, 10000), type="l", lwd=1.5, xlab="Time (generation)", ylab="Population size", cex.lab=1.5)
points(t, popsize, pch=15)
R <- 0
popsize <- N0*(1+R)^t
lines(t, popsize, lwd=1.5, col="red")
points(t, popsize, pch=16, col="red")
R <- -0.2
popsize <- N0*(1+R)^t
lines(t, popsize, lwd=1.5, col="blue")
points(t, popsize, pch=17, col="blue")

#図6.1(b)
r <- 0.2
N0 <- 2000
t <- seq(0, 30, 0.01)
popsize <- N0*exp(r*t)
plot(t, popsize, xlim=c(0, 30), ylim=c(0, 10000), type="l", lwd=4, xlab="Time (generation)", ylab="Population size", cex.lab=1.5)
r <- 0
popsize <- N0*exp(r*t)
lines(t, popsize, lty=3, lwd=3, col="red")
r <- -0.2
popsize <- N0*exp(r*t)
lines(t, popsize, lty=1, lwd=2, col="blue")

#図6.1(c)
R <- 0.3
N <- 500
K <- 6000
plot(x=0,y=0, xlim=c(0,30), ylim=c(0,10000), type="n", xlab="Time (generation)", ylab="Population size", cex.lab=1.5)
lines(c(0, 1), c(N, N+R*N*(K-N)/K), lwd=1.5)
points(0, N, pch=15)
for (t in 1:30){
        N <- N+R*N*(K-N)/K
    lines(c(t, t+1), c(N, N+R*N*(K-N)/K), lwd=2)
    points(t, N, pch=15)
    }
N <- 6000
lines(c(0, 1), c(N, N+R*N*(K-N)/K), lwd=2, lty=2, col="red")
points(0, N, pch=16, col="red")
for (t in 1:30){
        N <- N+R*N*(K-N)/K
    lines(c(t, t+1), c(N, N+R*N*(K-N)/K), lwd=2, lty=2, col="red")
    points(t, N, pch=16, col="red")
    }
N <- 10000
lines(c(0, 1), c(N, N+R*N*(K-N)/K), lwd=2, lty=3, col="blue")
points(0, N, pch=17, col="blue")
for (t in 1:30){
        N <- N+R*N*(K-N)/K
    lines(c(t, t+1), c(N, N+R*N*(K-N)/K), lwd=2, lty=3, col="blue")
    points(t, N, pch=17, col="blue")
    }   

#図6.1(d)
r <- 0.3
No <- 500
t <- seq(0, 30, 0.01)
K <- 6000
popsize <- K/(1+(K-No)*exp(-r*t)/No)
plot(t, popsize, type="l", lwd=4, las=1, ylim=c(0, 10000), xlim=c(0,30), xlab="Time (generation)", ylab="Population size", cex.lab=1.5)
No <- 10000
popsize <- K/(1+(K-No)*exp(-r*t)/No)
lines(t, popsize, lty=1, lwd=2, col="blue")
No <- 6000
popsize <- K/(1+(K-No)*exp(-r*t)/No)
lines(t, popsize, lty=3, lwd=3, col="red")