--- output: prettydoc::html_pretty: theme: architect highlight: vignette --- # Discrete Logistic Model: Bifurcation Diagram ## 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. ```{r} 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]) } par(mfrow = c(1, 1)) ``` 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. ``` {r} 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]) } par(mfrow = c(1, 1)) ``` 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$. ```{r} 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) } par(mfrow = c(1, 1)) ``` ## Task 4 Make a bifurcation diagram. ```{r} 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 = ".") } ```