Implement a chemostat model with one resource and Monod kinetics.
library("deSolve")
chemostat <- function(time, y, parms) {
with(as.list(c(y, parms)), {
mu <- mumax * P/(kp + P)
dAlg <- mu * Alg - D * Alg
dP <- D *(P0 - P) - 1/Y * mu * Alg
list(c(dAlg, dP), mu=mu)
})
}
parms <- c(
mumax = 0.5, # 1/d
kp = 0.01, # half saturation constant (mg/L)
Y = 41, # Yield coefficient
D = 0.1, # Dilution rate (1/d)
P0 = 0.05 # Phosphorus in the medium (mg/L)
)
times <- seq(0, 40, 0.1) # (d)
y0 <- c(Alg=2, P=0.05) # Phytoplankton C and Phosphorus P (mg/L)
out <- ode(y0, times, chemostat, parms)
plot(out)
Run simulations with different dilution rates D ans plot \(P\) and \(Alg\) as a function of \(D\).
As an alternative, you may also consider to write a loop in R with more than 10 intermediate steps.
out[nrow(out),]
library(rootSolve)
times <- c(0, Inf)
out <- stode(y0, times, chemostat, parms)
out
stode
(that solves the system with a Newton-Raphson method by setting the derivatives to zero) and runsteady
that runs the model numerically until an equilibrium is reached.For the chemostat model, an analytical solution would also be possible, i.e. mathematics with pencil and paper instead of the computer. Though this is in gereral the preferred method, we use simulation here as a didactical reasons, because an analytical solution is ofteh imposible for complex ecological models.