Search code examples
rparametersodedifferential-equationsdesolve

Change parameter values at time step in deSolve


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?


Solution

  • 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?