## Budworm model, roots

``````## =======================================================================
## 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
## =======================================================================

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))

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,
pt.bg = c("grey", "black", "white"))``````

## 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
## =======================================================================
# 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+0.2, length = 0.1)
arrows(0.05, equi+0.2, 0.05, equi-0.2, length = 0.1)
arrows(0.05, equi-0.2, 0.05, equi+0.2, length = 0.1)
arrows(0.05, equi+0.2, 0.05, equi-0.2, length = 0.1)

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

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