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

Posted 2023-02-18

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

【英文标题】：How to solve a system of ODE with time dependent parameters in R?

【发布时间】：2021-12-18 21:47:09

【问题描述】：

我正在尝试通过 deSolve 求解这个 ODE 系统，dX/dt = -X*a + (YX)b + c 和 dY/dt = -Ya + (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 参数的时变部分，因为我想不出一种方法来包含它并顺利运行它。我尝试将它包含在函数定义中，但无济于事。请帮忙。

【参考方案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 上的 short tutorial 中。

【参考方案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")