如何在 R 中求解具有时间相关参数的 ODE 系统？

Question

我试图通过 deSolve 求解这个常微分方程组，dX/dt = -X*a + (YX) b + c 和 dY/dt = -Y a + (XY)*b 时间为 [0,200]，a=0.30 ，b=0.2，但时间 [50,70] 时 c 为 1，否则为 0。我一直使用的代码是，

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)

我省略了 c 参数的随时间变化的部分，因为我无法找到包含它并顺利运行它的方法。我尝试将其包含在函数定义中，但无济于事。

Answer 1

标准方法是使用approxfun，即创建一个时间相关信号，我们也将其称为强制变量：

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)

另请注意 deSolve 特定plot函数以及时间因变量cc用作附加输出变量。

有关此内容的更多信息，请参阅：

• 在?forcings帮助页面和
• 在 Github 上的一个简短教程中。