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))
}
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'))
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)
rmarkdown::render('sheep-phaseR.Rmd')