9.1 科学计算案例：疫情传播仿真

9.1 科学计算案例：SIR 疫情传播仿真

用一个 SIR 模型走完"建模—求解—参数扫描—出图"全流程。涉及第 2 章的 struct 建模、第 6 章的 DifferentialEquations 与 Plots，代码不足五十行。

问题与模型

SIR 把人群分三箱：易感者 S、感染者 I、康复者 R。变化率由两个参数驱动：感染率 β 与康复率 γ。用第 2 章的方式给参数建模：

using DifferentialEquations, Plots

struct SIRParams
    β::Float64      # 感染率
    γ::Float64      # 康复率
    N::Float64      # 总人口
end

function sir!(du, u, p, t)
    S, I, R = u
    du[1] = -p.β * S * I / p.N
    du[2] =  p.β * S * I / p.N - p.γ * I
    du[3] =  p.γ * I
end
​

求解与可视化

p = SIRParams(0.4, 0.1, 1e6)
u0 = [p.N - 10, 10.0, 0.0]
prob = ODEProblem(sir!, u0, (0.0, 120.0), p)
sol = solve(prob)

plot(sol, labels = ["易感" "感染" "康复"],
     title = "SIR 传播曲线", xlabel = "天", ylabel = "人数")
​

进阶：参数扫描

问"感染率翻倍会怎样"，答案用扫描给出而不是手算：

curves = map([0.2, 0.4, 0.8]) do β
    prob = ODEProblem(sir!, u0, (0.0, 120.0), SIRParams(β, 0.1, 1e6))
    maximum(solve(prob)[2, :])          # 峰值感染人数
end
# β 越大峰值越高、到峰越早 —— 压平曲线的数学解释
​

案例结构一览

💡 关键直觉：把参数收进 struct 而不是散落成全局变量，扫描代码才可能干净——这是第 2 章"类型稳定"在第 5 章性能上的直接回报。

⚠️ 常见坑：总人口写成整数 1_000_000 而 u0 是浮点数组，相减会触发类型提升报错或隐式转换。数值仿真里统一用 Float64 是省心的约定。

结果解读与模型校验

算出曲线只是第一步，科学计算的一半功夫在"解读与校验"。三个校验动作：其一，总量守恒检验——三箱之和应恒等于总人口，数值解里逐时刻检查 sum(sol.u, dims=1) 的漂移量，漂移超过求解容差说明容差设松了；其二，解析极限对照——疫情结束后感染者应趋于零、康复者趋于一个稳定值，仿真终点值应符合这个定性结论，不符合多半是方程符号写错；其三，量纲自查——β 的单位是"每天"、γ 也是"每天"，两者之比 R0 = β/γ 是无量纲的基本再生数，本例 R0 = 4，大于 1 必然出现流行，这与曲线形状一致。三步都过了，模型才算"数值正确"而不只是"代码能跑"：

# 校验一：总量漂移
drift = maximum(abs.(sum.(sol.u) .- p.N))     # 应接近 0（如 1e-8 量级）

# 校验二：终点状态
sol[end]     # 第三箱（康复）应占绝大多数

# 校验三：基本再生数
p.β / p.γ    # 4.0 > 1，与曲线出现感染峰一致
​

变式扩展：从 SIR 到 SIRd 与干预仿真

真实建模很少停在基础版，两个方向值得动手做。方向一，加入死亡箱或潜伏箱（SEIR）：状态从 3 维升到 4 维，方程多两行，代码结构完全不变——这正体现了"方程函数与求解器分离"的设计红利，扩展模型不需要动求解框架。方向二，干预仿真：假设第 30 天起管控使 β 下降一半，把 β 写成时间的函数：

function sir_intervene!(du, u, p, t)
    S, I, R = u
    β = t < 30 ? p.β : p.β / 2        # 第 30 天起干预
    du[1] = -β * S * I / p.N
    du[2] =  β * S * I / p.N - p.γ * I
    du[3] =  p.γ * I
end
​

对比有无干预两条曲线的峰值与到峰时间，"压平曲线"四个字就有了定量版本。把干预时点也纳入参数扫描（第 20、30、40 天各扫一遍），还能看到"晚一周干预，峰值多多少"这类对决策直接有用的结论——科学计算的价值正是在这种"问一句、算一片"的循环里兑现的。

常见错误与排错

本案例现场翻车点提前列出。其一，初值量纲错：u0 里塞了 10 个感染者却忘了 N 写成 1e5，S 项算出负数，解曲线出现"负人口"——数值求解器不检查物理意义，守恒校验（三条曲线和恒为 N）能当场抓住。其二，β、γ 顺序传反：struct 字段顺序与构造调用对不上，模型行为完全不同但不报错；防法是构造时用具名关键字式写法或在 struct 上加一个"校验构造函数"，β、γ 均为正、N 大于初值人口。其三，容差不适配：默认容差在尖锐峰值附近积分误差偏大，峰值高度对政策解读敏感时显式收紧 reltol，并对比两种容差下的峰值差异来估计数值误差量级。这三条对应"物理校验、构造校验、数值校验"三道闸，闸闸都过，仿真结果才敢往报告里放。

数值细节：刚性系统与求解器选择

ODEProblem 默认求解器 Tsit5 是显式 Runge-Kutta，对本文这类平滑问题又快又准；但把模型改成 SIRd 加死亡箱、或引入快慢时间尺度后，系统会变"刚性"——显式方法被迫把步长压到极小才能稳住，求解时间暴涨。此时换隐式求解器，solve(prob, Rosenbrock23()) 或 Rodas5() 两步改动即可，DifferentialEquations 的自动选型也会在检测到刚性时切换。判断刚性的土办法：同一问题分别用默认与 Rosenbrock23() 求解，隐式求解器用时优势超过一个数量级就是刚性信号。另外每次求解都应检查 sol.retcode：返回 Success 才算收敛，出现 MaxIters 或 Unstable 时，先查方程里是否有除数为 0。

本节要点回顾

• struct 装参数、带 ! 函数写方程、solve 出解，三段式适用于任何 ODE 建模；
• 参数扫描回答"如果"类问题，比单次求解信息量大得多；
• 峰值指标（maximum(sol[2, :])）是传播模型的标配观测量；
• 守恒校验、解析极限对照、量纲自查三道闸过了，数值结果才可信；
• 扩展模型（SEIR、时变参数干预）不动求解框架，这是方程与求解器分离的红利。

这个案例还有个值得回味的对比：全部代码不到六十行，却覆盖了建模、求解、扫描、校验、扩展五个环节——Julia 科学计算的 productivity 不在于代码量少，而在于每个环节之间没有翻译成本，struct 直接进方程、方程直接进求解器、解直接进绘图，中间不需要任何胶水层。