Budworm model, bifurcation

## =======================================================================
## use of function 'gradient' and uniroot.all
## Stability of the budworm model, as a function of its
## rate of increase.
##
## Example from the book of Soetaert and Herman(2009)
## A practical guide to ecological modelling,
## using R as a simulation platform. Springer
## code and theory are explained in this book
## implemented by Karline Soetaert
## =======================================================================
require(rootSolve)   # uniroot.all and gradient
# parameters
r   <- 0.05
K   <- 10
bet <- 0.1
alf <- 1

# density-dependent growth and sigmoid-type mortality rate
rate <- function(x, r = 0.05) 
  r*x*(1-x/K) - bet*x^2/(x^2+alf^2)

# Stability of the roots ~ sign of eigenvalue of Jacobian 
stability <- function (r)  {
  Eq <- uniroot.all(rate, c(0, 10), r = r) # finds all roots
  eig  <- vector()
  for (i in 1:length(Eq)) 
      eig[i] <- sign (gradient(f = rate, x = Eq[i], r = r))
  return(list(Eq = Eq, Eigen = eig))
}

# bifurcation diagram for a range of r-values
rseq <- seq(0.01, 0.07, by = 0.0001)

plot(0, xlim = range(rseq), ylim = c(0, 10), type = "n",
     xlab = "r", ylab = "B*", main = "Budworm model, bifurcation",
     sub = "Example from Soetaert and Herman, 2009")

for (r in rseq) {
  st <- stability(r) 
  points(rep(r, length(st$Eq)), st$Eq, pch = 15,
         col = c("darkblue", "black", "lightblue")[st$Eigen+2]) 
}

legend("topleft", pch = 15, pt.cex = 2, c("stable", "unstable"),
        col = c("darkblue", "lightblue"))

equi <- uniroot.all(rate, c(0, 10), r = 0.05) 
arrows(0.05, 10         , 0.05, equi[4]+0.2, length = 0.1)
arrows(0.05, equi[3]+0.2, 0.05, equi[4]-0.2, length = 0.1)
arrows(0.05, equi[3]-0.2, 0.05, equi[2]+0.2, length = 0.1)
arrows(0.05, equi[1]+0.2, 0.05, equi[2]-0.2, length = 0.1)  

equi <- uniroot.all(rate, c(0, 10), r = 0.038) 
arrows(0.038, 10         , 0.038, equi[2]+0.1, length = 0.1)
arrows(0.038, equi[1]+0.1, 0.038, equi[2]-0.1, length = 0.1)

equi <- uniroot.all(rate, c(0, 10), r = 0.07) 
arrows(0.07, 10         , 0.07, equi[2]+0.2, length = 0.1)
arrows(0.07, equi[1]+0.2, 0.07, equi[2]-0.2, length = 0.1)