使用deSolve包ode()求解SIR模型并绘制相平面报错求助

错误原因分析 #

• 缺失SIR模型导数函数：ode()要求传入定义微分方程的函数，你未定义SIR_dydt，且该函数必须返回包含各变量导数的列表
• 初始条件违背种群守恒：N=100，但初始R设为1000，不满足S+I+R=N的约束，正确初始R应为0
• 绘图语法错误：output[, "S", "I"]的索引方式无效，且题目要求绘制(S,I)相平面（S为横轴、I为纵轴），而非时间序列图
• 字符串引号嵌套错误：main="SIR solutions"中双引号嵌套会触发语法报错，需改用单引号或调整引号层级

修正后的完整代码 #

library(deSolve)

# 定义SIR模型微分方程函数，必须返回列表格式的导数
SIR_dydt <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    dS <- -beta * S * I / N
    dI <- beta * S * I / N - gamma * I
    dR <- gamma * I
    # 返回由导数组成的列表
    return(list(c(dS, dI, dR)))
  })
}

# 参数与初始条件设置
beta <- 0.02
gamma <- 1
N <- 100
# 满足S+I+R=N，初始R设为0
initial <- c(S = 99.9, I = 0.1, R = 0)
parms <- c(beta = beta, gamma = gamma, N = N)

# 生成时间点序列
times <- seq(0, 20, by = 0.1)

# 求解ODE
output <- ode(y = initial, times = times, func = SIR_dydt, parms = parms)

# 绘制(S,I)相平面
plot(output[, "S"], output[, "I"], type = "l", 
     xlab = "易感者数量(S)", ylab = "感染者数量(I)", 
     main = "(S,I)相平面：SIR模型解", col = "darkred", lwd = 2)

代码说明 #

1. SIR导数函数：SIR_dydt函数接收时间t、当前状态y（S/I/R）和参数parms，通过with()简化变量调用，计算各变量的导数后以列表形式返回，这是ode()要求的标准格式
2. 初始条件修正：将R初始值设为0，保证S+I+R=100，符合种群总数N的约束
3. 相平面绘制：直接提取output中的S列和I列，绘制S随I变化的曲线，完全匹配题目要求的(S,I)相平面展示
4. 语法修正：修复了plot函数中的引号嵌套问题，调整坐标轴标签和标题以贴合相平面的含义

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