Change parameter values at time step in deSolve
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?
I get a slightly different result from you. If I uncomment the print(gamma0)
it prints out 5 twice then prints out 513 zeros. It is not hard to trace why in a superficial way, although you may want more than I will offer here.
Where you have the (commented out) statement print(gamma0)
instead, put the line:
cat("g:", gamma0, " t:", times, "\n")
and run the code. You will see that the first two times it displays are 0. Since those are on your list seq(0,10,1)
gamma0 is 5. After that, the times values displayed change. Notice that none of them that are printed are from your original list of times seq(from = 0, to = 10, by = 1/24)
and none of them are integers so none meet your condition to set gamma0 to 5. ode
is doing something with the times (interpolating?) but it is not simply using the values that you provided. In fact, it does not print out 241 values of gamma0 and times. It prints out 515 such values. I note that the result out
does have 241 values.
I think from your question that you assumed ode
would actually evaluate the function at your times
. It does not. It is treating times
like a continuous variable. But your condition
gamma0 <- ifelse(times %in% seq(0,10,1), 5, 0)
only tests for 11 specific values - not ranges of values. A continuous variable is quite unlikely to hit exactly those values.