## =======================================================================
## Example from the book of Soetaert and Herman(2009)
## A practical guide to ecological modelling -
## using R as a simulation platform. Springer
## implemented by Karline Soetaert
## =======================================================================
library(rootSolve) #uniroot.all, gradient
r <- 0.05
K <- 10
bet <- 0.1
alf <- 1
# the model : 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)
# uniroot.all finds all roots within an interval, her [0,10]
Eq <- uniroot.all(rate, c(0, 10))
# jacobian evaluated at all roots:
# This is just one value - and therefore jacobian = eigenvalue
# the sign of eigenvalue: stability of the root: neg=stable, 0=saddle, pos=unstable
eig <- vector()
for (i in 1:length(Eq))
eig[i] <- sign (gradient(rate, Eq[i]))
curve(rate(x), ylab = "dx/dt", from = 0, to = 10,
main = "Budworm model, roots",
sub = "Example from Soetaert and Herman, 2009")
abline(h = 0)
points(x = Eq, y = rep(0, length(Eq)), pch = 21, cex = 2,
bg = c("grey", "black", "white")[eig+2] )
legend("topleft", pch = 21, pt.cex = 2,
c("stable", "saddle", "unstable"),
pt.bg = c("grey", "black", "white"))
