Phase plane analysis of Sheep–Rabbit system

Here we will look in detail at the rabbit-sheep system (Strogatz, p155). We first load the packages that we need:

require(phaseR)
require(deSolve)

The state of the system is a two-dimensional vector, v=(r,s), where r is the rabbit population, and s is the sheep population. There are no extra parameters to pass to this system.

Sheep=function(t,y,parameters) {
  r=y[1]; s=y[2]
  drdt = r*(3 - r - (2*s))
  dsdt = s * (2 - r - s)
  list(c(drdt, dsdt))
}

Evolution from one initial condition

First we simulate the evolution of the system from one initial condition, and find that in this case, the Rabbits win, and the sheep lose.

init.cond = c(1.5, 1.3)
times=seq(from=0, to=20, length=100)
out = lsoda(init.cond, times, Sheep)
matplot(out[,1], out[,-1], type='l',xlab='time',
         ylab='y1,y2',lwd=2)
legend('topleft', c('Rabbits', 'Sheep'), lty=c(1,2), col=c('black', 'red'))

Phase plane analysis

Now we do the phase plane analysis. In the following, y0 is a matrix of different initial conditions to try out.

flow <- flowField(Sheep, x.lim=c(0,3), y.lim=c(0,3), points=19,
                  xlab='Rabbits', ylab='Sheep', add=F)
grid()
nullclines <- nullclines(Sheep, x.lim=c(-1, 3.5), y.lim=c(-1,3.5), points=500)

y0 <- matrix( c(0.5, 0.5,  2.5, 1.0, 0.5, 1.5,  2.5, 2.5), 4, 2, TRUE)
traj <- trajectory(Sheep, y0=y0, t.end=10)

## [1] "Note: colour has been reset as required"

From the nullclines, we can read off the four steady-state points which we hard code into ss matrix (one steady-state per row). We add an extra function mystablility to examine the classification of the steady sate, to return a suitable point style.

mystability = function(...) {
  st = stability(...)$classification
  res = switch(st,
               "Unstable node"=1,
               "Stable node"=19,
               "Saddle"=13)
  res
}

ss = matrix( c(0,0,  0,2,   3,0,  1,1),ncol=2, byrow=T)
stab = apply(ss, 1, function(s) {mystability(Sheep, y.star=s)})
## 
## T: 5   Delta: 6   Discriminant: 1   Classification: Unstable node
## T: -3   Delta: 2   Discriminant: 1   Classification: Stable node
## T: -4   Delta: 3   Discriminant: 4   Classification: Stable node
## T: -2   Delta: -0.9999998   Discriminant: 8   Classification: Saddle
flow <- flowField(Sheep, x.lim=c(0,3), y.lim=c(0,3), points=19,
                  xlab='Rabbits', ylab='Sheep', add=F,
                  asp=1,
                  las=1)
nullclines <- nullclines(Sheep, x.lim=c(-1, 3.5), y.lim=c(-1,3.5), points=500)
points(ss[,1], ss[,2], pch=stab,cex=3)

Compilation notes

rmarkdown::render('sheep-phaseR.Rmd')