I read this post by Corey Chivers at bayesianbiologist on plotting equations with chaotic behavior. I'd already written some code to do this for an earlier project so I thought I'd throw it up as well. I'd guess his is a bit more effecient than mine, but here it goes.
Here are two examples of plotting Bifurcation plots. The first uses R base graphics and both methods use the Ricker equation:
\( N_t = rN_{t-1}e^{-cN_{t-1}} \)
This plotting method relies on a nested looping structure. First we call an empty plot. You have to be careful and set the ylim variable to a reasonable size. rmax is the maximum value that our rate parameter will take on, and a is the strength of density dependence. We can set the grain size of our r values with the “by” parameter in the last line below.The smaller the value, the more values of r we iterate over.
rmax <- 30
plot(-1, -1, xlim = c(0, rmax), ylim = c(0, 1000), xlab = "r", ylab = "N")
a <- 0.01
r <- seq(0, rmax, by = 0.1)
Next we simply iterate a time series function (the Ricker equation in this case) n times. So we simply start with the first value of r, run the population to n and then select all the unique values after we think it's come to a steady state (40 is a reasonable point but you can change it) and just plot them as points! That's it, you're done
n <- 100
for (z in 1:length(r)) {
xl <- vector()
xl[1] <- 10
for (i in 2:n) {
xl[i] <- xl[i - 1] * r[z] * exp(-a * xl[i - 1])
}
uval <- unique(xl[40:n])
points(rep(r[z], length(uval)), uval, cex = 0.1, pch = 19)
}
Another approach would be to save the data from the previous loops and use ggplot to create our plot. This is pretty easy, and we can generate a nice plot with this code. Happy bifurcating.
library(ggplot2)
rmax <- 30
out.df <- matrix(NA, ncol = 2, nrow = 0)
a <- 0.01
r <- seq(0, rmax, by = 0.01)
n <- 100
for (z in 1:length(r)) {
xl <- vector()
xl[1] <- 10
for (i in 2:n) {
xl[i] <- xl[i - 1] * r[z] * exp(-a * xl[i - 1])
}
uval <- unique(xl[40:n])
### Here is where we can save the output for ggplot
out.df <- rbind(out.df, cbind(rep(r[z], length(uval)), uval))
}
out.df <- as.data.frame(out.df)
colnames(out.df) <- c("r", "N")
ggplot(out.df, aes(x = r, y = N)) + geom_point(size = 0.5)
