Fig_10

#---- Nicholson & Bailey + binomial negativa
# fig_10

# variables
N <- numeric()
P <- numeric()

# coeficientes
lambda <- 2
c <- 1.3
a <- 0.01
k <- 0.95

# condiciones simulación
ngen <- 500    # generaciones
N[1] <- 195    # N inicial
P[1] <- 200    # P inicial

# simulación
for (t in 1:ngen) {
  f <- (1 + a * P[t] / k)^(-k)
  N[t + 1] <- lambda * N[t] * f
  P[t + 1] <- lambda * c * N[t] * (1 - f)
}

# equilibrio
Peq <- k * (lambda^(1/k) - 1) / a
Neq <- Peq / (c * lambda * (1 - (1 + a * Peq / k)^(-k)))

# gráfico
pdf("fig_10.pdf", width = 7.5, height = 6)
par(pty = "s",
    bty = "n",
    mar = c(5, 2, 1, 1),
    xaxs = "i",
    yaxs = "i",
    las = 1)
plot(N, P,
    type = "l",
    lty = 2,
    col = "lightblue",
    xlim = c(0, 200),
    ylim = c(0, 400))
points(N, P,
    cex = 0.8,
    pch = 21,
    col = "royalblue",
    bg = "cyan")
abline(h = Peq,
    col = "orange",
    lwd = 0.5)
abline(v = Neq,
    col = "orange",
    lwd = 0.5)
dev.off()