C语言实现SIR疫情模型：核心函数设计与数值模拟实践164

作为一名资深的程序员，我非常理解在编程实践中，我们可能会遇到各种各样的函数命名，其中一些可能是标准库提供的，而另一些则是根据特定需求自定义的。关于您提及的“C语言sir函数”，首先需要明确指出，在标准C语言库中，并没有一个名为`sir`的内建函数。这通常意味着`sir`要么是一个用户自定义的函数，旨在完成特定任务；要么是一个在特定领域（如科学计算、嵌入式系统或某个第三方库）中约定俗成的函数名；又或者，它可能是一个拼写错误（例如，与`sin`函数、`sizeof`运算符等相似）。

鉴于“sir”这个缩写在科学计算领域，尤其是流行病学模型中有着一个非常著名的指代——即“Susceptible-Infected-Recovered”（易感者-感染者-康复者）模型，因此，本文将着重围绕“如何在C语言中设计并实现一个用于模拟SIR模型的核心函数，并将其命名为`sir`”这一主题展开深入探讨。这不仅能充分利用“sir”这个名称的深层含义，还能展示C语言在科学计算和数值模拟方面的强大能力。

我们将从SIR模型的基本原理出发，逐步介绍其在C语言中的数值离散化方法，核心`sir`函数的实现细节，以及一个完整的模拟程序框架。

SIR模型基础：疫情传播的数学描述

SIR模型是最经典的流行病学模型之一，它将人群分为三个相互关联的类别：

S (Susceptible)：易感者，指那些尚未感染疾病，但有可能被感染的个体。
I (Infected)：感染者，指已经感染疾病并能传播疾病的个体。
R (Recovered)：康复者（或移除者），指已经从疾病中康复并获得免疫力（或死亡）的个体，他们不再具有传染性，也不会再次被感染。

SIR模型通过一组常微分方程来描述这三类人群随时间的变化。假设总人口数N是恒定的（N = S + I + R），并且疾病是直接接触传播的。模型的数学表达式如下：

$\frac{dS}{dt} = -\beta \frac{S \cdot I}{N}$

$\frac{dI}{dt} = \beta \frac{S \cdot I}{N} - \gamma I$

$\frac{dR}{dt} = \gamma I$

其中：

$\beta$ (beta) 是感染率，表示每个感染者每天能够有效感染的易感者数量。
$\gamma$ (gamma) 是恢复率，表示每天康复的感染者所占的比例。其倒数 $1/\gamma$ 通常代表感染期的平均持续时间。
$\frac{S \cdot I}{N}$ 代表易感者与感染者接触并导致感染的概率，其中 $\frac{I}{N}$ 可以看作人群中感染者的密度。

这些微分方程描述了S、I、R的变化速度。由于这些方程通常没有解析解，我们需要使用数值方法进行模拟。

C语言实现SIR模型：核心函数 `sir` 的设计

在C语言中实现SIR模型，我们首先需要考虑如何表示模型的状态（S, I, R）和参数（$\beta$, $\gamma$, N）。结构体（struct）是理想的选择。

数据结构定义

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// 定义SIR模型在某一时刻的状态
typedef struct {
double S; // 易感者数量
double I; // 感染者数量
double R; // 康复者数量
} SIR_State;
// 定义SIR模型的参数
typedef struct {
double beta; // 感染率
double gamma; // 恢复率
double N; // 总人口数
} SIR_Params;

数值离散化：欧拉法（Euler Method）

由于SIR模型的微分方程是连续的，我们需要将其离散化才能在计算机上进行模拟。最简单的数值方法是欧拉法。它通过小步长来近似微分，即：

$X_{t+\Delta t} = X_t + \Delta t \cdot \frac{dX}{dt}$

其中 $\Delta t$ 是时间步长。将其应用于SIR模型，我们可以得到S、I、R在下一个时间步的状态：

$S_{t+\Delta t} = S_t + \Delta t \cdot (-\beta \frac{S_t \cdot I_t}{N})$

$I_{t+\Delta t} = I_t + \Delta t \cdot (\beta \frac{S_t \cdot I_t}{N} - \gamma I_t)$

$R_{t+\Delta t} = R_t + \Delta t \cdot (\gamma I_t)$

核心函数 `sir` 的实现

我们的核心 `sir` 函数将接收当前状态、模型参数和时间步长，然后计算并返回下一个时间步的状态。

/
* @brief 计算SIR模型在一个时间步长后的状态。
* 这个函数模拟了"SIR"模型的核心动力学。
* @param current_state 当前时刻的SIR状态。
* @param params SIR模型的参数（beta, gamma, N）。
* @param dt 时间步长。
* @return SIR_State 经过dt时间步长后的新状态。
*/
SIR_State sir_step(SIR_State current_state, SIR_Params params, double dt) {
SIR_State next_state;
// 计算当前时间步的S, I, R变化率
double dS_dt = - * current_state.S * current_state.I / params.N;
double dI_dt = * current_state.S * current_state.I / params.N - * current_state.I;
double dR_dt = * current_state.I;
// 使用欧拉法更新状态
next_state.S = current_state.S + dt * dS_dt;
next_state.I = current_state.I + dt * dI_dt;
next_state.R = current_state.R + dt * dR_dt;
// 确保S, I, R不出现负值（数值误差可能导致）
if (next_state.S < 0) next_state.S = 0;
if (next_state.I < 0) next_state.I = 0;
if (next_state.R < 0) next_state.R = 0;
// 确保总人口数保持近似常数（由于浮点数误差，可能会有微小偏差）
// 简单地重新归一化可以防止累积误差，但通常不严格必要，除非步长过大
double current_N_sum = next_state.S + next_state.I + next_state.R;
if (fabs(current_N_sum - params.N) > 1e-6) { // 如果总和偏离较大
// 可以选择进行归一化，例如：
// double scale = params.N / current_N_sum;
// next_state.S *= scale;
// next_state.I *= scale;
// next_state.R *= scale;
}
return next_state;
}

在这里，我们将函数命名为 `sir_step` 以明确其是计算一个时间步长的函数。在更广泛的语境中，一个管理整个模拟过程的函数也可以被称为 `sir_simulate` 或简称为 `sir`。

完整的SIR模拟程序框架

为了运行模拟并观察结果，我们需要一个 `main` 函数来初始化模型、循环调用 `sir_step` 函数，并输出结果。

初始化

首先，定义模型的初始状态和参数。

// 辅助函数：打印当前状态
void print_sir_state(int time_step, SIR_State state) {
printf("%4d | S: %10.2f | I: %10.2f | R: %10.2f", time_step, state.S, state.I, state.R);
}
// 辅助函数：将数据保存到CSV文件
void save_to_csv(FILE* fp, int time_step, SIR_State state) {
fprintf(fp, "%d,%.2f,%.2f,%.2f", time_step, state.S, state.I, state.R);
}
int main() {
// 1. 初始化模型参数
SIR_Params params = {
.beta = 0.2, // 感染率
.gamma = 0.05, // 恢复率 (平均感染期 = 1/0.05 = 20天)
.N = 1000.0 // 总人口
};
// 2. 初始化初始状态
SIR_State current_state = {
.S = params.N - 1.0, // 初始易感者：总人口减去少量感染者
.I = 1.0, // 初始感染者：1人
.R = 0.0 // 初始康复者：0人
};
// 3. 设置模拟参数
double total_time = 200.0; // 模拟总时长 (天)
double dt = 0.1; // 时间步长 (天)
int num_steps = (int)(total_time / dt); // 总步数
// 打印初始状态
printf("--- SIR Model Simulation ---");
printf("Initial State:");
print_sir_state(0, current_state);
printf("----------------------------");
printf("Step | S | I | R ");
printf("-----|------------|------------|-----------");
// 打开文件用于保存数据
FILE *fp = fopen("", "w");
if (fp == NULL) {
perror("Error opening file");
return EXIT_FAILURE;
}
fprintf(fp, "Time,S,I,R"); // CSV文件头
save_to_csv(fp, 0, current_state);
// 4. 模拟循环
for (int i = 1; i