# R Tools for Dynamical Systems ~ bifurcation plot in R for system of ODEs

[This article was first published on

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

**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

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 about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.