使用solve_ivp求解僵尸爆发时间时遇索引越界错误求助

我需要计算僵尸数量≥人类数量1/4的爆发时间随 zombification率γ的变化关系，使用的核心代码如下：

arrOutbreakT = []
arrGamma = np.linspace(10**(-4), 1, 2001)

for i in range (0, len(arrGamma)):
    gamma = arrGamma[i]
    NstepC = (gamma*10000).astype(int)
    tlistC = np.linspace(iniTime, finalTime, NstepC + 1)
    solC = solve_ivp(gradF, t_span2, ini2, t_eval = tlistC)
    
    k = 0
    klist = np.zeros(NstepC+1)
    while solC.y[1, k] < (solC.y[0, k] / 4):
        k = k + 1
        
    arrOutbreakT.append(tlistC[k])

相关参数与函数定义：

alpha = 0.0 # 人类自然出生率
beta = 1.0*(10**(-4)) # 人类自然死亡率
gamma = 9.5*(10**(-3)) # zombification率（全局初始值）
delta = 1.0*(10**(-4)) # 僵尸销毁率
epsilon = 1.0*(10**(-4))
iniTime = 0.0
h = 0.01 # 步长
Nstep = 500
finalTime = iniTime + Nstep*h # = 5.0天

H0 = 1000 # 初始人类数量
Z0 = 5 # 初始僵尸数量
D0 = 0 # 初始死亡人数

ini2 = [H0, Z0, D0]
t_span2 = [time[0], finalTime] # 此处time变量未定义，存在问题

def gradF (time, HZD):
    return [(alpha - beta - gamma*HZD[1])*HZD[0], 
            gamma*HZD[0]*HZD[1] - delta*HZD[0]*HZD[1] + epsilon/HZD[1], 
            beta*HZD[0] + delta*HZD[0]*HZD[1] - epsilon*HZD[2]]

运行后出现如下错误：

IndexError                                Traceback (most recent call last)
    321 k = 0
    322 klist = np.zeros(NstepC+1)
---> 323 while solC.y[1, k] < (solC.y[0, k] / 4):
    324 k = k + 1
    326 arrOutbreakT.append(tlistC[k])

IndexError: index 2 is out of bounds for axis 1 with size 2

问题分析 #

1. 索引越界直接原因：当γ取值较小时，NstepC = (gamma*10000).astype(int)的结果会非常小（比如γ=1e-4时，NstepC=1），此时tlistC仅包含2个时间点（索引0和1）。如果循环条件solC.y[1, k] < solC.y[0, k]/4在k=1时仍然成立，k会自增到2，而solC.y的列数只有2，索引超出范围导致报错。
2. 潜在参数问题：t_span2中使用了未定义的time[0]，应替换为iniTime；gradF函数依赖全局变量gamma，循环中修改gamma可能导致函数读取旧值（写法不规范）。

解决方案 #

1. 给循环添加边界限制 #

修改while循环条件，确保k不会超出solC.y的有效索引范围：

k = 0
max_k = len(solC.y[0]) - 1  # 获取最大有效索引
while k < max_k and solC.y[1, k] < (solC.y[0, k] / 4):
    k += 1
arrOutbreakT.append(tlistC[k])

2. 避免时间步长过少 #

设置NstepC的最小值，防止gamma过小时tlistC的点数不足：

NstepC = max(int(gamma * 10000), 50)  # 确保至少50个时间步

3. 修正未定义变量问题 #

把t_span2的定义改为：

t_span2 = [iniTime, finalTime]

4. 优化函数参数传递 #

将gamma作为参数传入gradF，避免依赖全局变量，修改循环内的求解代码：

def gradF(time, HZD, gamma):
    return [(alpha - beta - gamma*HZD[1])*HZD[0], 
            gamma*HZD[0]*HZD[1] - delta*HZD[0]*HZD[1] + epsilon/HZD[1], 
            beta*HZD[0] + delta*HZD[0]*HZD[1] - epsilon*HZD[2]]

# 循环内调用solve_ivp时，用lambda包装参数
solC = solve_ivp(lambda t, hzd: gradF(t, hzd, gamma), t_span2, ini2, t_eval=tlistC)

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