基于NumPy与Matplotlib的SIR模型RK4实现结果异常求助

你代码共有3处核心逻辑错误，修正后即可得到正确的SIR曲线：

• 1、RungeKutta4函数的return语句缩进错误：当前return被放在for循环内部，循环仅执行1次就直接返回结果，除第一个时间步外所有x值都为初始的0，是出现竖线的直接原因。需要将return语句调整到for循环的同级缩进位置。
• 2、SIR微分方程定义错误：感染人群I（对应x[1]）的变化率应为β*S*I - α*I，你当前代码误写为减去α乘以恢复人群R（x[2]），完全不符合SIR动力学逻辑。
• 3、四阶龙格库塔k4项参数错误：k4对应的计算时间点应为t[k]+dt，你写为t[k]+dt/2，会导致数值计算精度出现偏差。

修正后完整代码 #

import numpy as np
from matplotlib import pyplot as plt

# 模型定义
def SIR(x, params):
    alpha = params["alpha"]
    beta = params["beta"]
    S, I, R = x
    # 修正微分方程
    xdot = np.array([-beta*S*I, beta*S*I - alpha*I, alpha*I])
    return xdot

def RungeKutta4(f, x0, t0, tf, dt):
    t = np.arange(t0, tf, dt)
    nt = t.size
    nx = x0.size
    x = np.zeros((nx, nt))
    x[:, 0] = x0
    for k in range(nt-1):
        k1 = dt*f(t[k], x[:, k])
        k2 = dt*f(t[k]+dt/2, x[:, k]+k1/2)
        k3 = dt*f(t[k]+dt/2, x[:, k]+k2/2)
        # 修正k4的时间参数
        k4 = dt*f(t[k]+dt, x[:, k]+k3)
        dx = (k1+2*k2+2*k3+k4)/6
        x[:, k+1] = x[:, k]+dx
    # 修正return缩进位置
    return x, t

params = {"alpha": 0.001 , "beta": 46.48}
def f(t, x): return SIR(x, params)

if __name__ == "__main__":
    x0 = np.array([900, 1, 0])
    t0 = 0
    tf = 10
    dt = 0.0001
    x, t = RungeKutta4(f, x0, t0, tf, dt)
    plt.plot(t, x[0, :], "r", label="Susceptible")
    plt.plot(t, x[1, :], "g", label="Infected")
    plt.plot(t, x[2, :], "b", label="Recovered")
    plt.xlabel("Time(t)")
    plt.grid()
    plt.legend()
    plt.show()

注：原代码中冗余的plt.subplot(1,2,1)已移除，易感人群标签拼写错误Suceptible已修正为Susceptible

内容的提问来源于stack exchange，提问作者Baturay Çoksaygılı