How to solve a system of ODE with time dependent parameters in R?

Question

I am trying to solve this system of ODEs through deSolve, dX/dt = -X*a + (Y-X)b + c and dY/dt = -Ya + (X-Y)*b for time [0,200], a=0.30, b=0.2 but c is 1 for time [50,70] and 0 otherwise. The code I have been using is,

time <- seq(0, 200, by=1)
parameters <- c(a=0.33, b=0.2, c=1)
state <- c(X = 0, Y = 0)

    two_comp <- function(time, state, parameters){
      with(as.list(c(state, parameters)), {
        dX = -X*a + (Y-X)*b + c
        dY = -Y*a + (X-Y)*b
        return(list(c(dX, dY)))
      })
    }

out <- ode(y = state, times = time, func = two_comp, parms = parameters)
out.df = as.data.frame(out)

I have left out the time varying part of the c parameter since I can't figure out a way to include it and run it smoothly. I tried including it in the function definitions, but to no avail.

Answer 1

The standard way is to use approxfun, i.e. create a time dependent signal, that we also call forcing variable:

library("deSolve")
time <- seq(0, 200, by=1)
parameters <- c(a=0.33, b=0.2, c=1)
state <- c(X = 0, Y = 0)

two_comp <- function(time, state, parameters, signal){
  cc <- signal(time)
  with(as.list(c(state, parameters)), {
    dX <- -X * a + (Y - X) * b + cc
    dY <- -Y * a + (X - Y) * b
    return(list(c(dX, dY), c = cc))
  })
}

signal <- approxfun(x = c(0, 50, 70, 200), 
                    y = c(0, 1,  0,  0), 
                    method = "constant", rule = 2)

out <- ode(y = state, times = time, func = two_comp, 
           parms = parameters, signal = signal)
plot(out)

Note also the deSolve specific plot function and that the time dependent variable cc is used as an additional output variable.

More about this can be found:

• in the ?forcings help page and
• in a short tutorial on Github.

Answer 2

The interval limits where c is equal to 1 can be passed as parameters. Then, inside the differential function, use them to create a logical value

time >= lower & time <= upper

Since FALSE/TRUE are coded as the integers 0/1, every time this condition is false, c is multiplied by zero and the trick is done.

library(deSolve)

two_comp <- function(time, state, parameters){
  with(as.list(c(state, parameters)), {
    dX = -X*a + (Y-X)*b + c*(time >= lower & time <= upper)
    dY = -Y*a + (X-Y)*b
    return(list(c(dX, dY)))
  })
}

time <- seq(0, 200, by=1)
parameters <- c(a=0.33, b=0.2, c=1, lower = 50, upper = 70)
state <- c(X = 0, Y = 0)

out <- ode(
  y = state, 
  times = time, 
  func = two_comp, 
  parms = parameters
)

out.df <- as.data.frame(out)
head(out.df)

matplot(out.df$time, out.df[-1], type = "l", lty = "solid", ylim = c(0, 3))
legend("topright", legend = names(out.df)[-1], col = 1:2, lty = "solid")