如何在R中针对多组参数循环求解微分方程（SIR模型实例）

1. 确认单组参数SIR函数结构 #

假设你已有如下可运行的SIR求解函数（如果你的函数结构不同，只需保证输入参数能和参数数据框的列对应即可）：

library(deSolve)

# 单组参数SIR模型求解函数
solve_sir <- function(beta, gamma, N, I0, t_end) {
  S0 <- N - I0
  # 定义SIR微分方程
  sir_ode <- function(t, y, params) {
    with(as.list(c(y, params)), {
      dS <- -beta * S * I / N
      dI <- beta * S * I / N - gamma * I
      dR <- gamma * I
      return(list(c(dS, dI, dR)))
    })
  }
  # 初始条件与参数
  init <- c(S = S0, I = I0, R = 0)
  params <- c(beta = beta, gamma = gamma)
  # 求解ODE
  times <- seq(0, t_end, by = 1)
  out <- ode(y = init, times = times, func = sir_ode, parms = params)
  # 转为数据框返回
  return(as.data.frame(out))
}

2. 批量运行多场景参数 #

假设你的9个场景参数数据框名为scenario_params（列需包含模型所需的所有参数，如beta、gamma、N等），使用purrr的pmap函数按行传递参数，批量求解并合并结果：

library(purrr)
library(dplyr)
library(tidyr)

# 示例多场景参数数据框（替换为你的实际数据）
scenario_params <- data.frame(
  scenario = paste0("场景", 1:9),
  beta = c(0.3, 0.4, 0.5, 0.3, 0.4, 0.5, 0.3, 0.4, 0.5),
  gamma = c(0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3),
  N = rep(1000, 9),
  I0 = rep(10, 9),
  t_end = rep(100, 9)
)

# 批量求解所有场景
scenario_results <- scenario_params %>%
  mutate(sir_output = pmap(list(beta, gamma, N, I0, t_end), solve_sir)) %>%
  unnest(sir_output)

3. 保存结果 #

将合并后的结果保存为CSV文件，方便后续调用：

write.csv(scenario_results, "多场景SIR模型结果.csv", row.names = FALSE)

4. 用ggplot2和cowplot绘制对比图 #

分面展示每个场景的动力学曲线

library(ggplot2)
library(cowplot)

facet_plot <- ggplot(scenario_results, aes(x = time)) +
  geom_line(aes(y = S, color = "易感者(S)"), linewidth = 1) +
  geom_line(aes(y = I, color = "感染者(I)"), linewidth = 1) +
  geom_line(aes(y = R, color = "康复者(R)"), linewidth = 1) +
  facet_wrap(~scenario, ncol = 3) +
  scale_color_manual(values = c("blue", "red", "green")) +
  labs(x = "时间", y = "人数", color = "人群类别") +
  theme_bw() +
  theme(legend.position = "bottom")

同图对比所有场景的趋势

combined_plot <- ggplot(scenario_results, aes(x = time, color = scenario)) +
  geom_line(aes(y = S), linewidth = 0.8) +
  geom_line(aes(y = I), linewidth = 0.8) +
  geom_line(aes(y = R), linewidth = 0.8) +
  scale_color_viridis_d(name = "场景") +
  labs(x = "时间", y = "人数") +
  theme_bw() +
  theme(legend.position = "right")

组合两个图

plot_grid(facet_plot, combined_plot, labels = c("A", "B"), ncol = 1)

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