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

June 12, 2010
By

(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...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , ,

Recent popular posts

Contact us if you wish to help support R-bloggers, and place your banner here.

Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)