Add isoclines and/or direction field to plot

Question

I was solving differential equations with R when I came across a problem: I needed to add Isoclines and a direction field to my plot, but I didn't know how to. What package should I install/what function should I call/should I do it all manually?

Answer 1

Here is not an answer, but just a summary of the answers given in comments:

• the LVCompGames in primer package give this plot:

1. This is the solution given by @Ben Bolker. The solution is based on the graphics package and filled.contour function.

• Then @Josh O'Brien suggestion the rasterVis package with vectorplot function

I think all answers can be customized if you give us more information of what you wish as output.

Answer 2

There is a great package called phaseR. It uses the ODE solver from the package deSolve and let's you easily add nullclines, trajectories and so forth. Just be sure to name the parameters of your ODE function "t", "y" and "parameters". Otherwise the phaseR functions will give an error:

Error in deriv(t = 0, y = c(x[i], y[j]), parameters = parameters) :
unused arguments (y = c(x[i], y[j]), parameters = parameters)

Here is a working code example for the Lotka-Volterra predator-prey system:

require(deSolve) 
require(phaseR)

model.LV <- function(t, y, parameters){
  with(as.list(parameters),{

    N<-y[1] 
    P<-y[2]
    dN <- a*N - b*N*P  
    dP <- c*N*P -d*P

    list(c(dN,dP))
  })
}

params.LV<-c(a=0.4, b=0.3, c=0.1, d=0.2)

data.LV<-as.data.frame(lsoda(c(N=1,P=1),seq(1,250,by=0.5), model.LV, params.LV))
# plot the time series of both populations
plot(data.LV$time,data.LV$N, main="Time series of L-V equations", xlab="time", 
     ylab="Population densities N, P",
     type="l", col="green", 
     ylim=c(0,max(data.LV$N,data.LV$C)))
lines(data.LV$time,data.LV$P,col="red")

# plot the trajectories of the system
plot(data.LV$N, data.LV$P, type="l", col="blue", main="Trajectory of L-V equations",
     xlab="Prey density N", ylab="Predator density P", xlim=c(0,5), ylim=c(0,3))
#add Nullclines
nullclines(model.LV, x.lim=c(0.1,5),y.lim=c(0.1,3), parameters=params.LV, system="two.dim", colour=c("green","red"), add=TRUE)