###################################################
### chunk number 1: eval=FALSE
###################################################
## help.start()
###################################################
### chunk number 2:
###################################################
library(deSolve)
library(primer)
## Loading required package: lattice
library(lattice)
###################################################
### chunk number 3:
###################################################
b <- .5 # What is this and what are the units?
a <- .01 # What is this and what are the units?
e <- .05 # What is this and what are the units?
s <- .2 # What is this and what are the units?
###################################################
### chunk number 4:
###################################################
### Prey isocline
### 0 = bH -aHP; P=b/a
Hiso <- b/a
Hiso
## [1] 50
### Predator isocline
### 0 = eaHP - sP; H=s/(ea)
Piso <- s/(e*a)
Piso
## [1] 400
### draw a graph with axes from zero to 3 x the magnitude relevant isocline.
plot(x=NULL, ylab="Predators", xlab="Prey", xlim=c(0,3*Piso), ylim=c(0, 3*Hiso), type='n')
abline(h=Hiso, lty=2)
abline(v=Piso, lty=3)
###################################################
### chunk number 5:
###################################################
### Write the growth equations as R "expressions"
dHdt <- expression(b*H - a*P*H)
dPdt <- expression(e*a*P*H - s*P)
### PDEs, the effects of an individual of species i on the growth rate of species j
### Determine the partial derivatives of each population with
### respect to each other.
ddHdH <- D(dHdt, "H")
ddHdP <- D(dHdt, "P")
ddPdP <- D(dPdt, "P")
ddPdH <- D(dPdt, "H")
### Create an expression of expressions.
pdes <- expression(ddHdH, ddHdP, ddPdH, ddPdP)
### Tell R what H* and P* are
P <- b/a
H <- s/(e*a)
### Evaluate (at the equilibrium) the PDEs,
### and stick them into a matrix.
### This is the Jacobian matrix.
JacLV <- matrix(sapply(pdes, function(x) eval(eval(x))),
nrow=2, byrow=TRUE)
### Note on code - the nested expressions require nested "evaluation" - eval(eval(x))
JacLV
## [,1] [,2]
## [1,] 0.000 -4
## [2,] 0.025 0
###################################################
### chunk number 6:
###################################################
ev <- eigen(JacLV)$values
ev
## [1] 0+0.3162278i 0-0.3162278i
### Use just the imaginary parts to calculate this quantity
2*pi/Im(ev)
## [1] 19.86918 -19.86918
###################################################
### chunk number 7:
###################################################
time <- seq(from=0, to=20, by=0.5)
parameters <- list(b=b, a=a, e=e, s=s)
initial.N <- c(H=800, P=50)
### Now run the ODE solver and plot it.
out <- ode(y=initial.N, times=time, func=predpreyLV, parms=parameters)
plot(out[,"H"], out[,"P"], type='b',xlab="Prey",ylab="Predators")
abline(h=Hiso, lty=2)
abline(v=Piso, lty=3)
###################################################
### chunk number 8:
###################################################
initial.N <- c(H=650, P=50)
out <- ode(y=initial.N, times=time, func=predpreyLV, parms=parameters)
lines(out[,"H"], out[,"P"], type='b', pch=2)
###################################################
### chunk number 9:
###################################################
time <- 1:100
out <- ode(y=initial.N, times=time, func=predpreyLV, parms=parameters)
matplot(out[,"time"], out[,2:3], type="l", ylim=c(0, max(out)),xlab="time",ylab="Predators and Prey")
###################################################
### chunk number 10:
###################################################
predpreyRM
## function (t, y, p)
## {
## H <- y[1]
## P <- y[2]
## with(as.list(p), {
## dH.dt <- b * H * (1 - alpha * H) - w * P * H/(D + H)
## dP.dt <- e * w * P * H/(D + H) - s * P
## return(list(c(dH.dt, dP.dt)))
## })
## }
## <environment: namespace:primer>
###################################################
### chunk number 11: rm1
###################################################
w <- 5
D <- 150
alpha <- 0.001
###################################################
### chunk number 12:
###################################################
curve(b/w*(D+(1-alpha*D)*x - alpha*x^2),
from=0, to=1/alpha, lty=2, ylim=c(0,100),xlab="Prey",ylab="Predators")
abline(v=s*D/(e*w-s), lty=3)
###################################################
### chunk number 13:
###################################################
h <- s*D/(e*w-s)
h
## [1] 600
###################################################
### chunk number 14:
###################################################
p <- b/w*(D+(1-alpha*D)*h - alpha*h^2)
p
## [1] 30
###################################################
### chunk number 15:
###################################################
dhdt <- expression(b*h*(1-alpha*h) - w*p*h/(D+h))
dpdt <- expression(e*w*p*h/(D+h) - s*p)
### PDEs, the effects of an individual of species i
### on the growth rate of species j
### Determine the partial derivatives of each population with
### respect to each other.
ddhdh <- D(dhdt, "h")
ddhdp <- D(dhdt, "p")
ddpdp <- D(dpdt, "p")
ddpdh <- D(dpdt, "h")
### Create an expression of expressions.
pdes <- expression(ddhdh, ddhdp, ddpdh, ddpdp)
### Evaluate (at the equilibrium) the PDEs,
### and stick them into a matrix.
### This is the Jacobian matrix.
JacRM <- matrix(sapply(pdes, function(x) eval(eval(x))),
nrow=2, byrow=TRUE)
### Note on code - the nested expressions require nested "evaluation" - eval(eval(x))
JacRM
## [,1] [,2]
## [1,] -0.140 -4
## [2,] 0.002 0
###################################################
### chunk number 16:
###################################################
ev2 <- eigen(JacRM)$values
ev2
## [1] -0.07+0.05567764i -0.07-0.05567764i
2*pi/Im(ev2) # using just the imaginary parts
## [1] 112.8493 -112.8493
###################################################
### chunk number 17:
###################################################
time <- seq(from=0, to=113, by=1)
parameters <- list(b=b, w=w, alpha=alpha, D=D, e=e, s=s)
initial.N <- c(h=700, p=40)
out <- ode(y=initial.N, times=time, func=predpreyRM, parms=parameters)
plot(out[,"h"], out[,"p"], type='b',xlab="Prey",ylab="Predators")
curve(b/w*(D+(1-alpha*D)*x - alpha*x^2),
from=0, to=1/alpha, lty=2, add=T)
abline(v=s*D/(e*w-s), lty=3)
