当传染病模型遇上季节轮回：用 SIR 模型解析流行病的周期性奥秘-CSDN博客

本文链接：https://blog.csdn.net/xyhxx99499/article/details/148228422

🌍 从病毒传播到数学建模：为什么我们需要 SIR？

2020 年新冠疫情暴发时，全球都在关注一个词 ——“基本再生数 R₀”。这个看似抽象的数学概念，背后其实是传染病模型的魔力。今天要聊的 SIR 模型（易感者 Susceptible - 感染者 Infectious - 康复者 Recovered），正是流行病学最经典的动力学框架之一。它像一台 “病毒计算器”，用微分方程刻画人群中三类个体的动态变化，而我们今天要探讨的，是加入季节性因素的 “强制版 SIR”，看看季节轮回如何让病毒传播变得更具戏剧性～

📐 SIR 模型的基础骨架：三幕式传播剧本

先抛开代码，用生活化的场景理解 SIR 的核心逻辑：

易感者（S）：像等待被点燃的干柴，对病毒毫无抵抗力的人群；
感染者（I）：正在释放病毒的 “传播源”，每天接触易感者就可能引发新感染；
康复者（R）：成功 “通关” 的人，获得免疫力后退出传播链（这里假设永久免疫哦～）。

三者的关系可以用一组微分方程表示： \(\frac{dS}{dt} = \mu N - \beta(t) \frac{SI}{N} - \mu S\) \(\frac{dI}{dt} = \beta(t) \frac{SI}{N} - \gamma I - \mu I\) \(\frac{dR}{dt} = \gamma I - \mu R\)

是不是有点晕？别急～这里的\(\beta(t)\)就是今天的主角 ——时变传播率，而代码里的theta.beta0和theta.beta1，正是在为它注入季节的灵魂✨

🍂 当传播率跳起 “季节圆舞曲”：强制 SIR 的魔法参数

回到代码的参数部分，每一个数字背后都藏着现实意义：

theta.beta0 = 10/7;       % 基线传播率：平均每人每周感染10/7人（超过1就意味着疫情会扩散哦）
theta.beta1 = 0.05;       % 季节性振幅：传播率随季节波动的幅度，5%的变化看似微小，却能让疫情曲线翻天覆地  
theta.omega = 2*pi/365;   % 季节频率：每年2π的周期性变化，对应365天的周期（和地球公转同步~）
theta.gamma = 1/7;        % 康复率：平均7天康复，和常见呼吸道疾病的病程吻合  
theta.mu    = 1/(70*365); % 出生率/死亡率：70年的人口更替周期，模拟人口的自然更新

💡 小思考：为什么要加入出生率和死亡率？因为真实世界中，人群不是静态的～新生儿会增加易感者，老年人离世会减少总人口，这些 “背景音” 虽缓慢，却能在长时间尺度上影响疫情的最终走向。

⏳ 从 100 年到 10 年：两步模拟里的疫情演化哲学

代码中分成了两个阶段模拟，这可不是随便写的哦～

第一阶段：100 年的 “生态平衡”

tspan1 = 0: 1: 100 * 365;
[~, Y1] = ode45(@(t,y) sir_forced_ode(t,y,theta), tspan1, y0);
y_eq = Y1(end, :)';

这一步是让模型 “跑” 到稳定状态，就像让一个生态系统经历足够长的时间，直到出生、死亡、感染、康复达到某种动态平衡。初始条件y0 = [99999; 1; 0]很有意思：10 万人中只有 1 个感染者，其余都是易感者 —— 这像极了传染病暴发的早期场景，一个 “零号病人” 闯入了纯净的人群🌰。

第二阶段：从平衡态开始的 10 年追踪

tspan2 = 0: 1: 10 * 365;
[time2, Y2] = ode45(@(t,y) sir_forced_ode(t,y,theta), tspan2, y_eq);

为什么不从 0 开始直接模拟 110 年？因为前 100 年是为了 “过滤” 掉初始条件的影响，让我们看到纯粹由季节驱动的疫情波动。就像观察钟摆时，要等它摆脱初始推动，才能看到真正的周期性摆动～

📊 揭开模拟结果的面纱：感染人数的季节性狂舞

最后一步是绘图，也是最直观的部分：

plot(time2, Y2(:,2), 'LineWidth', 2);
xlabel('Time (days)'); ylabel('Number of Infectious (I)');
title('Forced SIR: Seasonal Infectious Dynamics');

想象一下这幅图的样子：感染人数（Y2 (:,2)）随时间（time2）起伏，不是单调的增长或衰减，而是呈现出规律的周期性波动。当theta.beta1=0.05时，传播率会在基线值附近上下浮动 5%，但这种微小的变化却能让疫情出现 “旺季” 和 “淡季”：

冬天可能因为室内聚集、空气干燥，传播率beta(t)达到峰值，感染人数飙升📈；
夏天则可能因通风良好、人群户外活动多，传播率下降，疫情进入低谷🌞。

这和我们观察到的流感季节性非常吻合！比如北半球的流感季通常在冬春季节，而南半球则在秋冬，正是季节相位差的体现～

🔍 深入 ODE 函数：代码背后的微分方程灵魂

最后来拆解核心的 ODE 函数sir_forced_ode：

function dydt = sir_forced_ode(t, y, p)
    S = y(1); I = y(2); R = y(3); N = S + I + R;
    beta_t = p.beta0 * (1 + p.beta1 * sin(p.omega * t));
    dS = p.mu * N - beta_t * S * I / N - p.mu * S;
    dI = beta_t * S * I / N - p.gamma * I - p.mu * I;
    dR = p.gamma * I - p.mu * R;
    dydt = [dS; dI; dR];
end

这里的关键是beta_t的表达式：1 + p.beta1 * sin(...)让传播率随时间呈正弦曲线变化。当t增加，sin(p.omega * t)在 - 1 到 1 之间震荡，带动beta_t在beta0*(1-beta1)到beta0*(1+beta1)之间波动。而dS、dI、dR的计算，正是对前文微分方程的代码直译：

易感者减少的原因：被感染（beta_t * S * I / N）或自然死亡（mu S）；
感染者增加的来源：易感者被感染，减少的原因：康复（gamma I）或死亡（mu I）；
康复者增加的来源：感染者康复，减少的原因：自然死亡（现实中康复者也会老去呀～）。

🚀 从模型到现实：强制 SIR 的应用与局限

这个模型能帮我们做什么？

预测季节性传染病：如流感、诺如病毒，提前估算高峰时间以调配医疗资源；
评估防控措施：假设疫苗接种降低了易感者比例，或社交距离降低了beta0，看看季节波动是否被削弱；
研究气候影响：如果全球变暖改变了季节模式（比如omega或beta1变化），疫情会如何演变？

但它的局限也很明显：

假设康复者永久免疫，但很多传染病（如新冠）会二次感染；
季节性仅用正弦函数模拟，现实中可能受湿度、温度、人群行为等多因素影响；
没有考虑空间扩散，把整个群体当作 “混合均匀” 的一锅粥，而真实世界有城市、乡村的差异。

🌟 结语：当数学遇见病毒，我们在计算什么？

这段代码背后，是人类用理性对抗不确定性的努力。SIR 模型像一个简化的 “病毒实验室”，让我们在计算机中推演疫情的可能轨迹。而加入季节性因素后，它更像一面镜子，映出自然规律与传染病的微妙互动 —— 不是所有疫情都会 “大杀四方”，季节的轮回可能成为抑制病毒的天然盟友🌿。

当然，模型永远只是现实的近似，真正的防疫还需要科学措施与人文关怀的结合。下次当你看到新闻里的疫情曲线时，或许可以想起这个带着sin函数的 SIR 模型，想象一下：在那些起伏的线条里，既有病毒的狡猾，也有季节的韵律，还有数学模型默默计算的智慧✨。

完整代码（省流版）

% Forced SIR model with seasonal transmission rate

% Parameters
theta.beta0 = 10/7;       % baseline transmission rate
theta.beta1 = 0.05;       % amplitude of seasonal forcing
theta.omega = 2*pi/365;   % seasonal frequency
theta.gamma = 1/7;        % recovery rate
theta.mu    = 1/(70*365); % birth/death rate

% Initial conditions
y0 = [99999; 1; 0]; % [S; I; R]

% 1) Simulate until equilibrium (100 years)
tspan1 = 0: 1: 100 * 365;
[~, Y1] = ode45(@(t,y) sir_forced_ode(t,y,theta), tspan1, y0);
y_eq = Y1(end, :)';

% 2) Simulate 10 years from equilibrium
 tspan2 = 0: 1: 10 * 365;
[time2, Y2] = ode45(@(t,y) sir_forced_ode(t,y,theta), tspan2, y_eq);

% Plot Infectious individuals over time
figure;
plot(time2, Y2(:,2), 'LineWidth', 2);
xlabel('Time (days)');
ylabel('Number of Infectious (I)');
title('Forced SIR: Seasonal Infectious Dynamics');
grid on;


% --- ODE function ---
function dydt = sir_forced_ode(t, y, p)
    S = y(1);
    I = y(2);
    R = y(3);
    N = S + I + R;
    beta_t = p.beta0 * (1 + p.beta1 * sin(p.omega * t));
    dS = p.mu * N - beta_t * S * I / N - p.mu * S;
    dI = beta_t * S * I / N - p.gamma * I - p.mu * I;
    dR = p.gamma * I - p.mu * R;
    dydt = [dS; dI; dR];
end