Rcode_8.r

#########################
# Rosenzweig - McArthur #
#########################
library(simecol)

## Loading required package: deSolve

rm_mod <- new("odeModel", 
              main= function(time, init, parms) { 
                with(as.list(c(init, parms)), {
                  dn1dt <- b * n1 * (1-n1/K) - w * n1 / (D + n1) * n2
                  dn2dt <- -d * n2 + w * e * n1 / (D + n1) * n2
                  list(c(dn1dt, dn2dt))
                })
              },
              parms = c(b = 0.7,K = 1000, w = 5, D = 400, d = 0.4, e = 0.15),
              times  = c(from = 0, to = 200, by = 0.1),
              init = c(n1 = 20, n2 = 20),
              solver= "lsoda"
)   

res_rm <- sim(rm_mod)

par(mfrow=c(1,1))
matplot(x = res_rm@out[,"time"], y = res_rm@out[,c("n1","n2")], type = "l", 
        lty = 1, lwd = 3, ylab = "abundance", xlab = "time", ylim = c(0,1000))
legend("top", ncol = 2, col = 1:2, lwd = 3, legend = c("prey", "predator"), bty = "n")

b <- 0.7
K <- 1000
w <- 5
D <- 400
d <- 0.4
e <- 0.15

# derivatives for N1 and N2
dN1dt <- expression(b * N1 * (1 - N1/K) - w * N2 * N1/(D + N1))
dN2dt <- expression(e * w * N2 * N1/(D + N1) - d * N2)

# partial derivation in R
(ddN1dN1 <-D(dN1dt, "N1"))

## b * (1 - N1/K) - b * N1 * (1/K) - (w * N2/(D + N1) - w * N2 * 
##     N1/(D + N1)^2)

(ddN1dN2 <- D(dN1dt, "N2"))

## -(w * N1/(D + N1))

(ddN2dN1 <- D(dN2dt, "N1"))

## e * w * N2/(D + N1) - e * w * N2 * N1/(D + N1)^2

(ddN2dN2 <- D(dN2dt, "N2"))

## e * w * N1/(D + N1) - d

# Jacobian matrix
J <- expression(
  matrix(
    c(eval(ddN1dN1), eval(ddN1dN2),
      eval(ddN2dN1), eval(ddN2dN2)),
    nrow=2, byrow=TRUE
  ) 
)

#trivial equilibrium solutions
N1 <- 0
N2 <- 0

J_01 <- eval(J)
eigen(J_01)$values

## [1]  0.7 -0.4

# unstable


# non-trivial equilibrium solutions
N1 <- d * D/(e * w -d)
N2 <- b/w * (D + N1*(1-D/K)-N1^2/K)

J_02 <- eval(J)
eigen(J_02)$values

## [1] -0.05866667+0.2597914i -0.05866667-0.2597914i

# real parts all negative, imaginary parts present => damped oscillation


# plot trajectories
plot(x = res_rm@out[,"n1"], y = res_rm@out[,"n2"], type = "l", xlab = "Prey", 
     ylab = "Predator", lwd = 3, col = "darkblue")
points(20,20, col = "darkblue", pch = 19, cex = 1.5)

# different starting values
init(rm_mod)["n1"] <- 60
init(rm_mod)["n2"] <- 30
res_rm <- sim(rm_mod)
lines(x = res_rm@out[,"n1"], y = res_rm@out[,"n2"], col = "red3", lty = 2, lwd=2)
points(60,30, col = "red3", pch = 19, cex = 1.5)

# matplot(x = res_rm@out[,"time"], y = res_rm@out[,c("n1","n2")], type = "l", 
#         lty = 2, lwd = 3, add = T)