如何在SIR模型中获取t=50时的易感、感染、康复人群数量？

一、提取t=50对应数值的方法 #

你通过np.linspace生成了均匀分布的时间向量time，可以用以下两种方式找到t=50对应的S、I、R值：

方法1：直接计算索引 #

因为时间范围是0到100，共1000个点，t=50正好是中间位置，索引为data//2（即500），直接提取：

t50_idx = data // 2
S50 = S[t50_idx]
I50 = I[t50_idx]
R50 = R[t50_idx]
print(f"S(50) = {S50:.4f}, I(50) = {I50:.4f}, R(50) = {R50:.4f}")

方法2：通用查找（适配非整数/非均匀时间点） #

如果时间点不是恰好落在生成的数组中，用numpy的函数找到最接近目标值的索引：

t50_idx = np.argmin(np.abs(time - 50))
S50 = S[t50_idx]
I50 = I[t50_idx]
R50 = R[t50_idx]
print(f"S(50) = {S50:.4f}, I(50) = {I50:.4f}, R(50) = {R50:.4f}")

二、代码中的关键错误修正 #

你的四阶龙格-库塔（RK4）公式存在两处错误，会导致计算结果完全偏离正确值：

1. 状态更新公式写错：比如计算S_k2时，正确的中间状态是S[i] + (h/2)*S_k1，你写成了S[i] + h + (1/2)*S_k1；
2. 中间时刻的状态未同步：计算S的中间导数时，未使用I的中间状态值，仍用初始的I[i]，不符合RK4的精度要求。

修正后的RK循环代码如下：

for i in range(data-1):
    # 计算S的RK4系数
    S_k1 = derivada_S(time[i], I[i], S[i])
    # 中间时刻的I状态需同步更新
    I_mid1 = I[i] + (h/2)*I_k1
    S_k2 = derivada_S(time[i] + h/2, I_mid1, S[i] + (h/2)*S_k1)
    
    I_mid2 = I[i] + (h/2)*I_k2
    S_k3 = derivada_S(time[i] + h/2, I_mid2, S[i] + (h/2)*S_k2)
    
    I_mid3 = I[i] + h*I_k3
    S_k4 = derivada_S(time[i] + h, I_mid3, S[i] + h*S_k3)
    
    S[i+1] = S[i] + (h/6)*(S_k1 + 2*S_k2 + 2*S_k3 + S_k4)

    # 计算I的RK4系数
    I_k1 = derivada_I(time[i], I[i], S[i])
    S_mid1 = S[i] + (h/2)*S_k1
    I_k2 = derivada_I(time[i] + h/2, I[i] + (h/2)*I_k1, S_mid1)
    
    S_mid2 = S[i] + (h/2)*S_k2
    I_k3 = derivada_I(time[i] + h/2, I[i] + (h/2)*I_k2, S_mid2)
    
    S_mid3 = S[i] + h*S_k3
    I_k4 = derivada_I(time[i] + h, I[i] + h*I_k3, S_mid3)
    
    I[i+1] = I[i] + (h/6)*(I_k1 + 2*I_k2 + 2*I_k3 + I_k4)
    
    # 计算R的RK4系数（同步使用I的中间状态）
    R_k1 = derivada_R(time[i], I[i])
    R_k2 = derivada_R(time[i] + h/2, I_mid1)
    R_k3 = derivada_R(time[i] + h/2, I_mid2)
    R_k4 = derivada_R(time[i] + h, I_mid3)
    
    R[i+1] = R[i] + (h/6)*(R_k1 + 2*R_k2 + 2*R_k3 + R_k4)

三、使用流程 #

将提取数值的代码放在绘图代码之前，运行后即可得到t=50时的S、I、R数值。

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