[This article was first published on mind of a Markov chain » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

As per request, here is the code that I wrote to draw bifurcation plots in R.

Bifurcation diagrams for discrete maps can be done using this code by James Jones. It is a little easier since approximation is not needed.

In the following code, I used the deSolve library to draw bifurcation diagrams for a system of ODEs (continuous). You basically need to choose the parameter you’d like to perturb (param.name), for a given range (param.seq). When plotting, it plots the chosen parameter range on the x-axis and the variable of choice (plot.variable) on the y-axis. This is an example for the Lotka-Volterra.

library(deSolve)

LotVmod <- function (Time, State, Pars) {
with(as.list(c(State, Pars)), {
dx = x*(alpha - beta*y)
dy = -y*(gamma - delta*x)
return(list(c(dx, dy)))
})
}

n <- 100 # number of simulations
param.name <- "gamma" # choose parameter to perturb
param.seq <- seq(0,1,length = 50) # choose range of parameters

Pars <- c(alpha = 1, beta = .001, gamma = 1, delta = .001)
Time <- seq(0, 10, length = n)
State <- c(x = .5, y = .9)

param.index <- which(param.name == names(Pars))
out <- list()
for (i in 1:length(param.seq))
out[[i]] <- matrix(0, n, length(State))

for (i in 1:length(param.seq)) {
# set params
Pars.loop <- Pars
Pars.loop[param.index] <- param.seq[i]
# converge
init <- ode(State, Time, LotVmod, Pars.loop)
# get converged points
out[[i]] <- ode(init[n,-1], Time, LotVmod, Pars.loop)[,-1]
}

range.lim <- lapply(out, function(x) apply(x, 2, range))
range.lim <- apply(do.call("rbind", range.lim), 2, range)
plot.variable <- "x" # choose which variable to show
plot(0, 0, pch = "", xlab = param.name, ylab = plot.variable,
xlim = range(param.seq), ylim = range.lim[,plot.variable])
for (i in 1:length(param.seq)) {
points(rep(param.seq[i], n), out[[i]][,plot.variable])
} Filed under: deSolve, Food Web, R        