如何在R中使用deSolve求解含时变参数的常微分方程(ODE)方程组

实现逻辑 #

deSolve的ODE求解函数会自动将当前迭代的时间值传入你定义的微分方程函数的time参数，你可以直接基于这个参数判断时变参数c的取值，无需额外复杂配置。

修正后可运行代码（直接判断法，适合阶跃参数） #

# 先加载依赖包
library(deSolve)

time <- seq(0, 200, by=1)
# 修正参数a为需求的0.30，c为动态参数无需放入固定参数列表
parameters <- c(a=0.30, b=0.2)
state <- c(X = 0, Y = 0)

two_comp <- function(time, state, parameters){
  with(as.list(c(state, parameters)), {
    # 按时间区间给c赋值
    c <- ifelse(time >= 50 & time <= 70, 1, 0)
    dX = -X*a + (Y-X)*b + c
    dY = -Y*a + (X-Y)*b
    # 可选将c作为第二个返回值，方便后续校验参数取值是否正确
    return(list(c(dX, dY), c = c))
  })
}

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

通用时变参数方案（插值法，适配任意非规则时变参数） #

如果后续c的取值逻辑更复杂，可以预先生成时间-参数对应序列，用插值函数调用：

library(deSolve)

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

# 预定义c的时间节点和对应取值
c_time_nodes <- c(0, 49.999, 50, 70, 70.001, 200)
c_value_nodes <- c(0, 0, 1, 1, 0, 0)
# 生成常数插值函数，rule=2表示超出节点范围用最近值填充
c_interp <- approxfun(c_time_nodes, c_value_nodes, method = "constant", rule = 2)

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

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

内容的提问来源于stack exchange，提问作者14thTimeLord