I am trying to solve an ODE in R using deSolve
. With the following code, I expected the parameter gamma0
takes the values 5 at time step 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 and 10, and 0 otherwise. However, the print(gamma0)
shows that gamma0
stays at 0.
Here is my ODE:
library(deSolve)
param <- c(a = 0.1, b = 1)
yini <- c(alpha0 = 6, beta0 = 2)
mod <- function(times, yini, param) {
with(as.list(c(yini, param)), {
gamma0 <- ifelse(times %in% seq(0,10,1), 5, 0)
## print(gamma0)
dalpha0 <- - a*alpha0 + gamma0
dbeta0 <- a*alpha0 - b*beta0
return(list(c(dalpha0, dbeta0)))
})}
times <- seq(from = 0, to = 10, by = 1/24)
out <- ode(func = mod, times = times, y = yini, parms = param)
plot(out, lwd = 2, xlab = "day")
What am I doing wrong?
This is a really simple modification of your function. If you're interested in knowing what you are doing wrong you can look below.
mod <- function(times, yini, param) {
dt = times[2] - times[1]
with(as.list(c(yini, param)), {
gamma0 <- ifelse(times <= 10*dt, 5, 0)
## print(gamma0)
dalpha0 <- - a*alpha0 + gamma0
dbeta0 <- a*alpha0 - b*beta0
return(list(c(dalpha0, dbeta0)))
})}
EDIT
Same as G5W's answer, the problem is with what you are comparing times
to.
When you are writing
times %in% seq(0,10,1)
you are not referring to time steps. You simply refer to the values of times
.
So, if you want to have it for the first 10 time steps you just need to go with my code or anything that considers the dt
.
But here's a question for you:
If you do not need gamma0
to change according to times
and want it to be 5 at the first 11 (10) time steps, why you are comparing it to times
? Why not simply set it to be 5 for those time steps?