使用deSolve求解疟疾延迟微分方程模型未达预期平衡求助

问题分析 #

你的代码存在两个关键问题，导致状态值未达到预期平衡反而骤降：

1. 参数未正确提取：模型函数中直接使用a、b等参数，但未从传入的parms列表中解析，可能依赖全局变量导致参数不匹配。
2. 延迟项处理错误：当t < taum时，错误地用当前时刻的状态值（y[1]/y[2]/y[3]）代替t-taum时刻的历史值，违背了延迟微分方程的逻辑，导致动力学行为异常。

修正后的代码 #

library(deSolve)

model1 <- function(t, y, parms) {
  # 从parms中提取参数，避免全局变量干扰
  with(as.list(parms), {
    # 直接用lagvalue获取延迟值，dede会自动处理t-taum早于初始时刻的情况
    lagIh <- lagvalue(t - taum, 1)
    lagEm <- lagvalue(t - taum, 2)
    lagIm <- lagvalue(t - taum, 3)
    
    dIh <- (a*b*m*y[3])*(1 - y[1]) - r*y[1]
    dEm <- (a*c*y[1])*(1 - y[2] - y[3]) - (a*c*lagIh)*(1 - lagEm - lagIm)*exp(-mu2*taum) - mu2*y[2]
    dIm <- (a*c*lagIh)*(1 - lagEm - lagIm)*exp(-mu2*taum) - mu2*y[3]
    
    list(c(dIh, dEm, dIm))
  })
}

times <- seq(-15, 10000, by = 1)

# 初始状态
y <- c(Ih = 0.1,
       Em = 0.1,
       Im = 0.1)

# 参数列表
parms <- c(a = 0.5,
           b = 0.5,
           c = 0.5,
           m = 20,
           r = 0.05,
           mu2 = 0.05,
           taum = 5)

# 求解延迟微分方程
output <- dede(y = y, times = times, func = model1, parms = parms)

# 绘图
plot(output, xlab = "time", ylab = "State value")

修正说明 #

• 参数解析：用with(as.list(parms), { ... })将参数列表转为局部环境变量，确保函数内部使用的是传入的参数，避免全局变量引发的参数混乱。
• 延迟项修正：移除错误的ifelse判断，直接调用lagvalue(t - taum, idx)获取历史状态值。dede会自动处理t-taum早于初始时刻（times[1] = -15）的情况，默认返回初始状态的常值，符合延迟微分方程的历史条件设定。

运行修正后的代码，状态值会逐渐收敛到平衡态，符合预期。

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