Discrete Logistic Model: Bifurcation Diagram
The Rmd assignment solution file is availabe here.
Task 1
Simulate \(N_{t+1} = N_t + rN_t(1 - \frac{N_t}{K})\) using \(K = 500\), \(N_0 = 100\), \(r = 1.8, 2.3, 2.45, 2.56, 2.8\). Compare the results.
K <- 500
r <- c(1.8, 2.3, 2.45, 2.56, 2.8)
t <- 1:100
N <- matrix(nrow = length(t), ncol = length(r))
N[1, ] <- 100
for(i in 1:length(r)) {
for(j in 1:(length(t)-1)) {
N[j+1, i] <- N[j, i] + r[i] * N[j, i] * (1 - (N[j, i] / K))
}
}
colnames(N) <- paste("r =", r)
par(mfrow = c(2, 2))
for(i in 1:ncol(N)) {
plot(t, N[, i], type = "l", ylab = "N", main = colnames(N)[i])
}
Series with r = 1.8, converges to 500, while other series fluctuate around some number and the size of those fluctuations increases as r increases.
Task 2
Try \(N_0 = 101\) and compare with Task 1.
M <- matrix(nrow = length(t), ncol = length(r))
M[1, ] <- 101
for(i in 1:length(r)) {
for(j in 1:(length(t)-1)) {
M[j+1, i] <- M[j, i] + r[i] * M[j, i] * (1 - (M[j, i] / K))
}
}
colnames(M) <- paste("r =", r)
par(mfrow = c(2, 2))
for(i in 1:ncol(M)) {
plot(t, M[, i], type = "l", ylab = "M", main = colnames(M)[i])
}
Results are very similar compared to Task 1.
Task 3
Plot time step 1 return map for \(r = 1.8, 2.3, 2.45, 2.56, 2.8\).
N <- matrix(nrow = length(t), ncol = length(r))
N[1, ] <- 0.0000001
for(i in 1:length(r)) {
for(j in 1:(length(t)-1)) {
N[j+1, i] <- N[j, i] + r[i] * N[j, i] * (1 - (N[j, i] / K))
}
}
colnames(N) <- paste("r =", r)
par(mfrow = c(2, 2))
for(i in 1:ncol(N)) {
plot(N[1:99, i], N[2:100, i], type = "s",
xlab = expression(N[t]), ylab = expression(N[t+1]), main = colnames(N)[i],
xlim = c(0, 1000), ylim = c(0, 1000))
x <- 0:1000
y <- x + r[i] * x * (1 - (x / K))
lines(x, y)
}
Task 4
Make a bifurcation diagram.
r <- seq(1.8, 2.8, by = 0.001)
t <- 1:200
data <- matrix(nrow = length(t), ncol = length(r))
for(i in 1:length(r)) {
# random initial condition
N <- runif(1, 50, 150)
# first converge to attractor
for(j in t){
N <- N + r[i] * N * (1 - (N / K))
}
# collect points on attractor
for(j in t){
N <- N + r[i] * N * (1 - (N / K))
data[j, i] <- N
}
}
# plot the bifurcation diagram
plot(r, data[1, ], pch = ".", xlab = "r", ylab = "N")
for(i in 2:length(t)) {
points(r, data[i, ], pch = ".")
}