Problema 6.2 del libro Small System Dynamics Models for Big Issues: Triple Jump towards Real-World Complexity de Erik Pruyt.
Suppose there is a muskrat plague in a particular area. At first, there were 100 muskrats. The autonomous increase in the number of muskrats per muskrat per year amounts to an average of 20 muskrats per muskrat per year. Suppose that each year, 10 licences are granted to set muskrat traps. The licences are only valid for one year and each person holding a licence may set 10 traps. Assume the number of muskrats caught per trap is proportional to the number of muskrats and a catch rate per trap which is close to 0.2, say, on average 0.2, minimally 0.195, and maximally 0.205.
La variable de estado (la población de muskrats) puede crecer sin limite, estabilizarse o desaparecer dependiendo de la tasa de captura. Existe un punto en donde la población se mantiene pero con la minima variación el sistema cambia y por lo tanto cambia la población de muskrats.
library("deSolve")
# Variables exógenas o parámetros
parameters <- c(autonomous.increase = 0.2, # muskrats/muskrat/year
licences.granted = 10, # licences/year
traps.per.licence = 10, # traps/licence
catch.rate.per.trap = 0.2) # muskrats/trap/year
# Condiciones iniciales de las variables de estado (aquí se escriben sin "d")
InitialConditions <- c(muskrats = 100)
# Tiempo de modelación
times <- seq(0, # tiempo inicial
10, # fin de la simulación
0.1) # intervalo de tiempo
# Método de integración
intg.method <- c("rk4")
# Función del caso
muskrat.plague <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# Variables endógenas
traps.set <- licences.granted * traps.per.licence # traps/year
catch.rate <- catch.rate.per.trap * muskrats # muskrats/trap/year
# Variables de flujo
increase.rate <- autonomous.increase * muskrats # muskrats/year
catch.rate.total <- catch.rate * traps.set # muskrats/year
# Variables de estado (siempre inician con un "d")
dmuskrats <- increase.rate - catch.rate.total
list(c(dmuskrats)) # En la lista deben ir las variables que quieres graficar, dentro de la c() las variables de estado (las que inician con "d") y fuera de la c(), pero dentro del paréntesis de la lista, las demás variables que se deseen plotear.
})
}
# Simulación
muskrat.simulation <- ode(y = InitialConditions,
times = times,
func = muskrat.plague,
parms = parameters,
method = intg.method)
# Gráfica
plot(muskrat.simulation, col = c("blue"))