使用scipy solve_ivp的Radau方法时出现数组索引过多错误

错误根源 #

刚性求解器（如Radau、BDF）在计算数值雅可比矩阵时，会向你的dHdt函数传入二维数组格式的H参数（用于批量计算多个扰动点的导数，以近似雅可比矩阵）。但你的dHdt函数始终返回标量值，而非与输入H同维度的数组，导致后续数组运算时维度不匹配，触发IndexError。

而非刚性求解器（如RK45）通常单次仅传入一维数组形式的H，此时返回标量可以被隐式转换为兼容的数组格式，因此不会报错。

修复方法 #

修改dHdt函数，确保返回值与输入H的形状一致。可以通过以下两种方式实现：

方式1：直接返回数组形式结果 #

import numpy as np

def dHdt(t, H):
    if H>H_p(t) and data2_sleep(t)>0:
        return np.array([-0.323*24])
    elif H>H_p(t) and data2_sleep(t)==0:
        return np.array([0.116*24])
    elif H<=H_m(t) and data2_sleep(t)>0:
        return np.array([-0.278*24])
    elif H<=H_m(t) and data2_sleep(t)==0:
        return np.array([0.150*24])
    elif H<=H_p(t) and  H>H_m(t) and data2_sleep(t) > 0:
        return np.array([-0.274*24])
    elif H<=H_p(t) and  H>H_m(t) and data2_sleep(t) == 0:
        return np.array([0.096*24])

方式2：统一处理返回值维度 #

import numpy as np

def dHdt(t, H):
    if H>H_p(t) and data2_sleep(t)>0:
        res = -0.323*24
    elif H>H_p(t) and data2_sleep(t)==0:
        res = 0.116*24
    elif H<=H_m(t) and data2_sleep(t)>0:
        res = -0.278*24
    elif H<=H_m(t) and data2_sleep(t)==0:
        res = 0.150*24
    elif H<=H_p(t) and  H>H_m(t) and data2_sleep(t) > 0:
        res = -0.274*24
    elif H<=H_p(t) and  H>H_m(t) and data2_sleep(t) == 0:
        res = 0.096*24
    return np.asarray(res).reshape(H.shape)

修改后，无论输入H是一维还是二维数组，函数都会返回匹配维度的导数数组，满足刚性求解器的计算要求。

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