0

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.

Laura
  • 35
  • 5
  • You can use `approxfun` to define a time-dependent input, see https://tpetzoldt.github.io/deSolve-forcing/deSolve-forcing.html or search StackOverflow for a related question (there are several). For more specific help, please post a [minimal reproducible example](https://stackoverflow.com/a/5963610/3677576). – tpetzoldt Jan 26 '22 at 21:44

1 Answers1

1

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.

forcing

tpetzoldt
  • 5,338
  • 2
  • 12
  • 29
  • Thanks this is really helpful and it works for my first influx! I am wondering now how I can incorporate multiple different influxes in my model. Do I make a flux2 and flux3 etc, incorporate them in my differential equations, and then add them as additional arguments to ode (or dede) as flux=c(flux1, flux2, flux3)? – Laura Jan 31 '22 at 10:01
  • There are several ways how to do this: use several fluxes (flux1, flux2, flux3, but put it in a list, not a vector). Another possible way is [approxTime](https://www.rdocumentation.org/packages/simecol/versions/0.8-13/topics/approxTime) from package simecol, and you can also indirect vector indexing to make it faster. In case of interest, write a new question and don't forget to add a minimal reproducible example. – tpetzoldt Jan 31 '22 at 21:34