每日算法 — 使用java实现SIR模型：微分方程数值求解与疫情传播模拟

SIR模型是流行病学中最经典的仓室模型，将人群划分为易感者(Susceptible)、感染者(Infected)和康复者(Recovered)三个 compartment，通过常微分方程组描述疾病在人群中的传播动态。本文用Java实现SIR模型，并对比欧拉法与四阶龙格-库塔法(RK4)的数值精度。

模型建立

微分方程组

设总人口 $N = S + I + R$ 恒定，传播参数为：
– $\beta$：感染率，每个感染者每天有效接触并传染的人数
– $\gamma$：康复率，每天康复人数占感染者的比例（$1/\gamma$ 为平均感染周期）

$$
\begin{cases}
\frac{dS}{dt} = -\frac{\beta S I}{N} \
\frac{dI}{dt} = \frac{\beta S I}{N} – \gamma I \
\frac{dR}{dt} = \gamma I
\end{cases}
$$

基本再生数

$R_0 = \beta / \gamma$ 表示一个感染者在完全易感人群中平均能传染的人数。当 $R_0 > 1$ 时疫情扩散，$R_0 < 1$ 时疫情自然消退。

Java实现

/**
 * SIR流行病学模型
 * 支持欧拉法与四阶龙格-库塔法数值求解
 */
public class SIRModel {

    // 模型参数
    private final double beta;   // 感染率
    private final double gamma;  // 康复率
    private final double N;      // 总人口

    public SIRModel(double beta, double gamma, double N) {
        this.beta = beta;
        this.gamma = gamma;
        this.N = N;
    }

    /** 计算 dS/dt */
    public double dSdt(double S, double I) {
        return -beta * S * I / N;
    }

    /** 计算 dI/dt */
    public double dIdt(double S, double I) {
        return beta * S * I / N - gamma * I;
    }

    /** 计算 dR/dt */
    public double dRdt(double I) {
        return gamma * I;
    }

    /**
     * 欧拉法前向迭代
     * @param S0, I0, R0 初始状态
     * @param dt 时间步长
     * @param steps 迭代步数
     * @return 每一步的 [S, I, R] 数组
     */
    public double[][] euler(double S0, double I0, double R0, double dt, int steps) {
        double[][] result = new double[steps + 1][3];
        double S = S0, I = I0, R = R0;
        result[0] = new double[]{S, I, R};

        for (int i = 1; i <= steps; i++) {
            double sNew = S + dt * dSdt(S, I);
            double iNew = I + dt * dIdt(S, I);
            double rNew = R + dt * dRdt(I);
            S = sNew; I = iNew; R = rNew;
            result[i] = new double[]{S, I, R};
        }
        return result;
    }

    /**
     * 四阶龙格-库塔法 (RK4)
     * 精度为 O(dt⁴)，远优于欧拉法的 O(dt)
     */
    public double[][] rk4(double S0, double I0, double R0, double dt, int steps) {
        double[][] result = new double[steps + 1][3];
        double S = S0, I = I0, R = R0;
        result[0] = new double[]{S, I, R};

        for (int n = 1; n <= steps; n++) {
            // k1
            double k1S = dSdt(S, I);
            double k1I = dIdt(S, I);
            double k1R = dRdt(I);

            // k2
            double k2S = dSdt(S + 0.5 * dt * k1S, I + 0.5 * dt * k1I);
            double k2I = dIdt(S + 0.5 * dt * k1S, I + 0.5 * dt * k1I);
            double k2R = dRdt(I + 0.5 * dt * k1I);

            // k3
            double k3S = dSdt(S + 0.5 * dt * k2S, I + 0.5 * dt * k2I);
            double k3I = dIdt(S + 0.5 * dt * k2S, I + 0.5 * dt * k2I);
            double k3R = dRdt(I + 0.5 * dt * k2I);

            // k4
            double k4S = dSdt(S + dt * k3S, I + dt * k3I);
            double k4I = dIdt(S + dt * k3S, I + dt * k3I);
            double k4R = dRdt(I + dt * k3I);

            S += dt / 6.0 * (k1S + 2 * k2S + 2 * k3S + k4S);
            I += dt / 6.0 * (k1I + 2 * k2I + 2 * k3I + k4I);
            R += dt / 6.0 * (k1R + 2 * k2R + 2 * k3R + k4R);
            result[n] = new double[]{S, I, R};
        }
        return result;
    }

    /** 计算基本再生数 R₀ */
    public double getR0() {
        return beta / gamma;
    }

    /** 计算群体免疫阈值 */
    public double herdImmunityThreshold() {
        return 1.0 - 1.0 / getR0();
    }

    public static void main(String[] args) {
        // 模拟参数：总人口100万，初始1人感染
        // R0 = 2.5 (类似流感), 感染周期7天 -> gamma = 1/7
        double N = 1_000_000;
        double beta = 2.5 / 7.0;  // R0 = beta / gamma = 2.5
        double gamma = 1.0 / 7.0;
        double S0 = N - 1, I0 = 1, R0_init = 0;
        double dt = 0.1;  // 天
        int days = 200;
        int steps = (int) (days / dt);

        SIRModel model = new SIRModel(beta, gamma, N);
        System.out.printf("基本再生数 R₀ = %.2f%n", model.getR0());
        System.out.printf("群体免疫阈值 = %.1f%% %n%n", model.herdImmunityThreshold() * 100);

        double[][] rk4Result = model.rk4(S0, I0, R0_init, dt, steps);

        // 输出关键时间节点
        System.out.println("天数\t\t易感者\t\t感染者\t\t康复者");
        for (int d = 0; d <= days; d += 20) {
            int idx = (int) (d / dt);
            double[] state = rk4Result[idx];
            System.out.printf("第%3d天\t%.0f\t\t%.0f\t\t%.0f%n", d, state[0], state[1], state[2]);
        }

        // 找到感染峰值
        double maxI = 0;
        int peakDay = 0;
        for (int i = 0; i <= steps; i++) {
            if (rk4Result[i][1] > maxI) {
                maxI = rk4Result[i][1];
                peakDay = (int) (i * dt);
            }
        }
        System.out.printf("%n感染峰值出现在第 %d 天，约 %.0f 人同时感染%n", peakDay, maxI);
    }
}

数值方法对比

对于SIR模型，RK4仅需欧拉法约 $1/100$ 的步长即可达到相同精度，在200天模拟中可将计算量从 $20000$ 步降至 $2000$ 步。

扩展：SEIR模型

在SIR基础上增加潜伏期 compartment (E, Exposed)：

$$
\frac{dE}{dt} = \frac{\beta S I}{N} – \sigma E, \quad \frac{dI}{dt} = \sigma E – \gamma I
$$

其中 $\sigma$ 为潜伏期转阳率。SEIR更适用于具有明显潜伏期的新冠等传染病建模。

SIR模型虽然简化，但揭示了疫情传播的核心规律：控制 $R_0$ 是遏制疫情的唯一途径，而疫苗接种通过降低易感者比例实现群体免疫。