MATLAB has a nice program called pplane that draws phase planes of differential equations models. pplane on MATLAB is an elaborate program with an interactive GUI where you can just type the model to draw the phase planes. The rest you fidget by clicking (to grab the initial conditions) and it draws the dynamics automatically.

As far as I know, R doesn’t have a program of equal stature. R’s GUI itself is non-interactive (maybe because creating a good GUI require money), and you can’t fiddle around with the axes graphically, for example. The closest I could find was code from Prof. Kaplan from Macalester College in his program, pplane.r.

Below is a slight modification of his program that uses the deSolve package for a more robust approximation of the trajectory, and I made it so you can draw the trajectories by clicking, using the locator() function.

The pplane.r program takes in a 2D differential equation model, initial values and parameter value specifications to draw the dynamics on a plane. It draws arrows at evenly spaced out points at a certain resolution to see the general shape of the dynamics. This is done by using a crude method to create the Jacobian matrix. The next step is to give in initial values to draw the trajectory.

The only changes that made were to change the phasetraj() function, which draws the trajectories after you’ve made the arrow plot. Instead of using a self-made Runge Kutta method, I replaced it with a more robust ode() from the deSolve package. I also made it possible to point click multiple points (initial values) to draw the trajectories from. The code is shown below and it’s a little bit redundant because of different model specifications between pplane.r and deSolve, but it works. Also, the nullclines() function that draws the nullclines seems to not be working for whatever reason.

I could make the code more coherent, but I am lazy. Point is, the code can reproduce the pplane package in MATLAB to the best of my knowledge.

When I run draw.traj(), it asks for the number of initial points you’d like to give (denoted by loc.num). If that number is 5, I click the graph 5 times, and it automatically runs the model. I run a predator-prey model from a previous post with parameters: .

I could specify the color of the trajectory, and its time range. As analyzed by linearization, the predator-prey dynamics of this model is a center (you could see the dynamics go in a circle with the initial value up top). One can imagine running all kinds of 2D models. It’s not as interactive as the MATLAB version, but I think it works well enough as a first step.

Code is follows (please source in pplane.r and deSolve package):

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)))
	})
}

# please source() from http://www.macalester.edu/~kaplan/math135/pplane.r
nullclines(predatorprey(alpha, beta, gamma, delta),c(-10,100),c(-10,100),40)
phasearrows(predatorprey(alpha, beta, gamma, delta),c(-10,100),c(-10,100),20);

# modification of phasetraj() in pplane.r
draw.traj <- function(func, Pars, tStart=0, tEnd=1, tCut=10, loc.num=1, color = "red") {
	traj <- list()
	print(paste("Click", loc.num, "initial values"))
	x0 <- locator(loc.num, "p")
	for (i in 1:loc.num) {
		out <- as.data.frame(ode(func=func, y=c(x=x0$x[i], y=x0$y[i]), parms=Pars, times = seq(tStart, tEnd, length = tCut)))
		lines(out$x, out$y, col = color)
		traj[[i]] <- out
	}
	return(traj)
}

alpha = 1; beta = .001; gamma = 1; delta = .001

nullclines(predatorprey(alpha, beta, gamma, delta),c(-10,100),c(-10,100),40)
phasearrows(predatorprey(alpha, beta, gamma, delta),c(-10,100),c(-10,100),20, col = "grey")
draw.traj(func=LotVmod, Pars=c(alpha = alpha, beta = beta, gamma = gamma, delta = delta), tEnd=10, tCut=100, loc.num=5)