SEIRD-gov微分方程组Python改写及时间步对应ind序列传入问题求解

我明白你的问题了——在Julia版本里，si(t)函数会根据当前时间t计算出对应的索引，从而获取self_isolation的对应值，但在Python的odeint调用中，你没办法直接让每个时间步自动匹配ind序列的对应元素。下面是完整的解决方案，完全对应原Julia逻辑，同时修复时间步索引匹配的问题：

核心修改点 #

• 重构c_sig函数，让它能根据当前时间t动态获取ind序列的对应元素（处理Python 0-based索引与Julia 1-based的差异）
• 调整函数参数传递，确保odeint求解时每个时间步都能拿到正确的self_isolation值
• 增加索引越界保护，避免时间t超过ind序列长度时出错

完整Python代码实现 #

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

# 假设你已经加载了df_region数据，比如：
# df_region = pd.read_csv("your_data.csv")

# 1. 生成ind序列（和原逻辑一致）
ind = df_region['self_isolation'].apply(lambda x: int(x)).values

# 2. 定义根据时间t获取对应ind值的辅助函数（对应Julia的si(t)）
def get_current_ind(t):
    # 对应Julia的convert(Int, round(t + 1))，转换为Python 0-based索引
    idx = int(round(t))
    # 防止索引越界：如果t超过ind序列长度，取最后一个元素
    if idx >= len(ind):
        idx = len(ind) - 1
    return ind[idx]

# 3. 修正c_sig函数，动态获取当前时间的ind值
def c_sig(t, c_1, c_2):
    current_ind = get_current_ind(t)
    sig = 1 / (1 + math.exp(c_1 * (current_ind - c_2)))
    return sig

# 4. gov函数保持原逻辑不变
def gov(t, tg, g_1, g_2):
    if t > tg:
        alpha = 1 - g_1
    else:
        alpha = 1 - g_2
    return alpha

# 5. beta_gov函数调整，不再传入整个ind序列
def beta_gov(t, beta_0, c_1, c_2, tg, g_1, g_2):
    beta_rez = beta_0 * gov(t, tg, g_1, g_2) * c_sig(t, c_1, c_2)
    return beta_rez

# 6. SEIRD_gov微分方程组函数
def SEIRD_gov(y, t, beta_0, c_1, c_2, sigma, gamma, dr, ro, tg, g_1, g_2, N):
    S, E, I, R, D = y
    current_beta = beta_gov(t, beta_0, c_1, c_2, tg, g_1, g_2)
    dSdt = -current_beta * S * I / N
    dEdt = current_beta * S * I / N - sigma * E
    dIdt = sigma * E - (1 - dr) * gamma * I - dr * ro * I
    dRdt = (1 - dr) * gamma * I
    dDdt = dr * ro * I
    return dSdt, dEdt, dIdt, dRdt, dDdt

# ---------------------- 求解与绘图部分 ----------------------
# 示例参数（你可以替换成自己的参数）
N = 100000  # 总人口
S0 = N - 1
E0 = 1
I0 = 0
R0 = 0
D0 = 0
beta_0 = 0.3
c_1 = 0.1
c_2 = 50
sigma = 1/5.2  # 潜伏期倒数
gamma = 1/2.9  # 感染期倒数
dr = 0.02      # 死亡率比例
ro = 1/7       # 死亡时间倒数
tg = 30        # 政策切换时间点
g_1 = 0.4      # 政策切换后的干预强度
g_2 = 0.1      # 政策切换前的干预强度

# 时间数组（假设和ind序列长度对应）
t = np.linspace(0, len(ind)-1, len(ind))

# 初始条件
y0 = (S0, E0, I0, R0, D0)

# 调用odeint求解（注意参数顺序，新增N参数）
ret = odeint(SEIRD_gov, y0, t, args=(beta_0, c_1, c_2, sigma, gamma, dr, ro, tg, g_1, g_2, N))
S, E, I, R, D = ret.T

# 绘制结果曲线
plt.figure(figsize=(12, 8))
plt.plot(t, S, label='Susceptible')
plt.plot(t, E, label='Exposed')
plt.plot(t, I, label='Infected')
plt.plot(t, R, label='Recovered')
plt.plot(t, D, label='Dead')
plt.xlabel('Time (days)')
plt.ylabel('Number of Individuals')
plt.title('SEIRD-gov Model Simulation')
plt.legend()
plt.grid(True)
plt.show()

关键说明 #

• 索引匹配：get_current_ind函数完全对应原Julia的si(t)逻辑，通过round(t)将时间t转换为0-based索引，解决了Python和Julia索引规则的差异
• 参数传递：不再向beta_gov和c_sig传递整个ind序列，而是通过闭包的方式让get_current_ind访问ind数组，避免了odeint无法逐时间步传参的问题
• 越界保护：增加了索引越界判断，确保当求解时间超过ind序列长度时，不会抛出索引错误，而是取最后一个元素的值

这样修改后，模型就能像原Julia版本一样，每个时间步都能正确获取对应self_isolation值，求解结果也会和原逻辑一致。

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