I am making a population distribution model with ODEs and DDEs in DeSolve with the method lsoda. In this model I want to have an influx of the population at a specific time (a specific day). A very simple example:
dn1dt=influx - mortality
The influx (x) needs to occur at time (t) = y (in days). If it is not day y, I don't want an influx.
At the moment I have written the influx as influx=function(t,y,x){ifelse((t==y), x, 0))
, but I am running into the problem of a changing time step due to the method I am using (lsoda). Due to the changing time step, I will not hit the specific time (y) that will trigger the influx.
I am just starting out to work with a changing time step, so I am not sure how to deal with this problem. Please let me know if anything is unclear.
The influx can be implemented as time dependent input (also called "forcing", that lasts a certain time, e.g. a day) or as event, where a state variable is changed immediately. Here an example using a rectangular signal as forcing. Interpolation of time is implemented with approxfun
. The flux function can simply be added as additional argument to ode
and to the model
:
library("deSolve")
model <- function(t, y, p, flux) {
influx <- flux(t)
dN <- influx - p["d"] * y["N"]
list(dN, influx = influx)
}
y0 <- c(N = 1)
p <- c(d = 0.2)
times <- seq(0, 10, length.out = 100)
flux <- approxfun(x = c(0, 5, 6, 10),
y = c(0, 1, 0, 0), rule = 2, method = "constant")
out <- ode(y = y0, times = times, model, parms = p,
flux = flux, method = "lsoda")
plot(out, las=1)
More can be found on the ?forcings
help page or at this page.