如何在SIR模型ODE系统中实现gamma随时间线性从0.4增至0.8？

我正在基于SIR模型编写ODE系统代码，期望参数gamma的值随时间t线性递增，从初始值0.4逐步以微小增量上升，直至时间序列结束（t=1000）时达到0.8。请问能否在当前的ODE系统框架内实现这一需求？

原代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint  # 原代码遗漏该导入，需补充

# 总人口数，N
N = 1
# 初始感染人数I0和康复人数R0
I0, R0 = 0.001, 0
# 其余个体初始均为易感人群U0
U0 = N - I0 - R0
J0 = I0
Lf0, Ls0 = 0, 0
# 接触率beta，平均康复率gamma（单位：1/天）
beta, gamma = 8, 0.4
int_gamma = np.linspace(0.4, 0.8, 1000+1)
mu, muTB, sigma, rho = 1/80, 1/6, 1/6, 0.03
u, v, w = 0.88, 0.083, 0.0006
t = np.linspace(0, 1000, 1000+1)

# SIR模型微分方程
def deriv(y, t, N, beta, gamma, mu, muTB, sigma, rho, u, v, w):
    U, Lf, Ls, I, R, cInc = y
    b = (mu * (U + Lf + Ls + R)) + (muTB * I)
    lamda = beta * I
    clamda = 0.2 * lamda
    dU = b - ((lamda + mu) * U)
    dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
    dLs = (u * Lf) - ((w + clamda + mu) * Ls)
    dI = w*Ls + v*Lf - ((gamma + muTB + sigma) * I) + (rho * R)
    dR = ((gamma + sigma) * I) - ((rho + clamda + mu) * R)
    cI = w*Ls + v*Lf + (rho * R)
    return dU, dLf, dLs, dI, dR, cI


# 在时间网格t上积分SIR方程
solve = odeint(deriv, (U0, Lf0, Ls0, I0, R0, J0), t, args=(N, beta, gamma, mu, muTB, sigma, rho, u, v, w))
U, Lf, Ls, I, R, cInc = solve.T

解决方案：可以实现 #

完全能在现有框架内实现gamma随时间线性递增的需求，核心是在微分方程函数deriv内部，根据当前传入的时间t实时计算对应gamma值，具体修改如下：

修改后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 总人口数，N
N = 1
# 初始感染人数I0和康复人数R0
I0, R0 = 0.001, 0
# 其余个体初始均为易感人群U0
U0 = N - I0 - R0
J0 = I0
Lf0, Ls0 = 0, 0
# 接触率beta，gamma的初始值、终值和时间序列终点
beta = 8
gamma_start = 0.4
gamma_end = 0.8
t_max = 1000
mu, muTB, sigma, rho = 1/80, 1/6, 1/6, 0.03
u, v, w = 0.88, 0.083, 0.0006
t = np.linspace(0, t_max, 1000+1)

# SIR模型微分方程
def deriv(y, t, N, beta, gamma_start, gamma_end, t_max, mu, muTB, sigma, rho, u, v, w):
    U, Lf, Ls, I, R, cInc = y
    # 实时计算当前时间t对应的gamma值，线性递增
    gamma = gamma_start + (gamma_end - gamma_start) * (t / t_max)
    b = (mu * (U + Lf + Ls + R)) + (muTB * I)
    lamda = beta * I
    clamda = 0.2 * lamda
    dU = b - ((lamda + mu) * U)
    dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
    dLs = (u * Lf) - ((w + clamda + mu) * Ls)
    dI = w*Ls + v*Lf - ((gamma + muTB + sigma) * I) + (rho * R)
    dR = ((gamma + sigma) * I) - ((rho + clamda + mu) * R)
    cI = w*Ls + v*Lf + (rho * R)
    return dU, dLf, dLs, dI, dR, cI


# 在时间网格t上积分SIR方程
solve = odeint(deriv, (U0, Lf0, Ls0, I0, R0, J0), t, args=(N, beta, gamma_start, gamma_end, t_max, mu, muTB, sigma, rho, u, v, w))
U, Lf, Ls, I, R, cInc = solve.T

关键修改说明 #

1. 移除预生成的gamma数组：无需提前创建int_gamma，改为在每个时间步实时计算gamma值
2. 动态计算gamma：通过线性公式gamma = gamma_start + (gamma_end - gamma_start) * (t / t_max)，确保t=0时gamma为0.4，t=1000时达到0.8
3. 更新参数传递：将原固定gamma参数替换为gamma_start、gamma_end和t_max，让函数能动态适配时间变化

修改后，ODE求解器每一步计算微分方程时，都会使用对应时间点的gamma值，完美实现线性递增需求。

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