Scipy odeint 在负时间值处评估

Question

我正在尝试从最后时间的边界条件中使用scipy.integrate.odeint 求解 ODE 系统，并一直工作到初始时间（如此处所述：Backward integration in time using scipy odeint）。

但是，odeint 迭代到负时间值 - 超出了我正在寻找的实际解决方案范围 - 这会导致错误，因为我的 ODE 取决于时间的平方根，而我的函数返回一个复数值一个实数。

这是一个重现问题的示例：

    import numpy as np 
    from scipy.integrate import odeint
    tmax = 4e4
    tmin = 1
    t = np.linspace(tmax,tmin,1e3)
    param0 = [1] #value of x at tmax

    def func(param,t):
        x = param[0]
        dxdt = 1e-10/np.sqrt(t)
        print(t) #show what values of t are tried by odeint
        return dxdt

    res = odeint(func,param0,t)

很快print(t) 显示负值，res 被nan 填充。

有没有办法防止odeint 变为负值？为什么它在我的输入数组 t 之外尝试值？

对于我的实际代码，我找到了一些避免nan 结果的方法，例如在函数中添加if t<0: t=0（这在此处不起作用）或通过施加非常小的最大时间步长（hmax<tmin）但这会使计算时间更长。

请注意，Scipy odeint Non-negative solution 是相关的，但我的问题略有不同：我不关心否定解决方案 x，而是关心否定论点 t。

Answer 1

我建议使用 solve_ivp 而不是旧的 odeint。请注意，solve_ivp 期望 f(t, y) 而不是 f(y, t) 作为 ode 的右手边：

import numpy as np 
from scipy.integrate import solve_ivp
tmax = int(4e4)
tmin = 1
t = np.linspace(tmax, tmin, int(1e3))
param0 = [1] #value of x at tmax

def func(t, param):
    dxdt = 1e-10/np.sqrt(t)
    return dxdt

res = solve_ivp(func, y0=param0, t_span=[tmax, tmin], t_eval=t)

其中t_eval=t 确保解决方案将存储在时间点t。否则求解器选择时间点。