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

Question

我正在尝试通过 deSolve 求解这个 ODE 系统，dX/dt = -X*a + (YX) b + c 和 dY/dt = -Y a + (XY)*b 时间 [0,200], a=0.30 , b=0.2 但 c 在时间 [50,70] 时为 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 上的简短教程中。

Answer 2

c等于 1的区间限制可以作为 传递parameters。然后，在微分函数内部，使用它们创建一个逻辑值

time >= lower & time <= upper

由于FALSE/TRUE被编码为整数0/1，因此每次该条件为假时，c都将乘以零，然后技巧就完成了。

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")