Change parameter values at time step in deSolve

Question

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?

Answer 1

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?

Answer 2

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.

Answer 3

M--'s answer doesn't works for me, what if you just try this?

mod <- function(times, yini, param) {
  with(as.list(c(yini, param)), {

    if (times < 10) {
      gamma0 = 5
    } else if (times >= 10) {
      gamma0 = 0
    }

    dalpha0 <- -a * alpha0 + gamma0
    dbeta0 <- a * alpha0 - b * beta0
    return(list(c(dalpha0, dbeta0)))

  })
}

This works as a Headvise type function