(This article was first published on

**mind of a Markov chain » R**, and kindly contributed to R-bloggers)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

To

**leave a comment**for the author, please follow the link and comment on their blog:**mind of a Markov chain » R**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...