在 SciPy 中动态结束 ODE 集成

计算科学 Python scipy 网格
2021-12-23 02:15:47

我有一条穿过时空的光线,即一条曲线R4,由某个变量 λ 参数化。确切的轨迹,即曲线的坐标函数由一些 ODExμ(λ)dxμdλ=f(λ)

现在,我的代码要满足两个要求:

  1. 结果必须在规则网格上给出,即步长 dλ 必须在任何地方都相同。xμ(λ)

  2. 我不知道区间 λ 预先运行,即积分必须在满足某个条件后立即动态结束(基本上我正在测试坐标是否在我不知道时空属性的某个坐标区间之外)。g(λ,xμ(λ))xμ(λ)

要求 2 显然不包括scipy.odeint(). 这给我留下了一个类,如果我选择算法scipy.ode,它允许我通过solout回调将中断条件传递给积分器。dopri5dop853

不过,这有点复杂:由于这些(和所有其他)积分器动态调整它们的内部步长(这是不手动进行积分的一种点),如果我愿意,我不能简单地存储传递给回调的值满足要求1。(更新:事实上,这样做是没有意义的。请参阅下面的我的注释(*)。)相反,我仍然需要循环λ并scipy.ode.integrate()为每个步骤调用方法然而,这意味着我必须以某种方式停止循环:λ+dλ

from scipy.integrate import ode

def deriv(λ, x, params):
    # …
    return dx_dλ

end_reached = False

def solout(λ, x):
    if g(λ, x):
        # Stop integration
        # Returning -1 will NOT make solver.successful() return False, so
        # we need to inform the integration loop (see below) separately.
        nonlocal end_reached
        end_reached = True
        return -1

    # Everything's ok. Move on.
    return 0

λ0 = 0
x0 = …
params = …
dλ = 0.1

solver = ode(deriv)
solver.set_integrator('dopri5', nsteps=500)
solver.set_solout(solout)
solver.set_initial_value(x0, λ0)
solver.set_f_params(params)

λ = [λ0]
x = [x0]

while True:
    solver.integrate(solver.t + dλ)

    if not solver.successful() or end_reached:
        break

    λ.append(solver.t)
    x.append(solver.y)

总而言之,这是相当丑陋的。不仅逻辑基本上无处不在(例如,我无法测试它所属的循环中的条件)。但是由于solout回调(我想),它也至少比没有回调慢一个数量级。更新:事实上,即使没有回调,dopri5积分器也比vode(至少对于我想使用这种方法的另一个 ODE)慢几个数量级。

我是否遗漏了什么,或者这真的是满足上述要求的唯一方法吗?我还尝试引发异常deriv()(并随后在循环中捕获它)以停止积分器。虽然这适用于所有集成商,但它仍然会导致一堆错误被打印到标准输出。

(*) 据我了解,传递给的值deriv()并不是一个可靠的近似值,而只是求解器返回的最终值。例如,考虑到(快速)变化的斜率,求解器可能会进行过冲以选择合理的步长,从而得出最佳估计。

4个回答

您需要的功能在 Matlab ODE 求解器包中称为事件位置,或在SUNDIALS求解器套件术语中称为寻根。本质上,此功能允许在自由变量和因变量的某些向量函数具有根的点处停止积分。即,对于系统 具有事件定位能力的 ODE 求解器将是能够找到某个函数何时过零——这个过零过程称为事件

dydt=f(t,y),yRn,
g(t,y)

不幸的是,scipy.integrate缺乏支持事件位置的良好 ODE 求解器。您可以使用求解器solout的功能,但在这种情况下您会过冲,因为只有在 step 完成后才会调用。相反,当 ODE 求解器能够进行事件定位时,它将努力寻找的根的确切位置,以便更精确地停止积分。dopri5soloutg(t,y)

我的计算也需要这个功能。我使用以下包https://github.com/bmcage/odes这个包很难编译,因为它依赖于 BLAS、LAPACK 和 SUNDIALS。但是,接口(API)非常好,它比使用fromsolout的功能要好得多dopri5scipy.integrate

解决方案是 https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

来自文档:用于非刚性问题的“RK45”或“RK23”方法和用于刚性问题的“Radau”或“BDF”方法

来自 scipy 的文档:

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, **options)[来源] 解决系统的初始值问题常微分方程。

参数:. . . .

events:可调用的,或可调用的列表,可选的要跟踪的事件。如果 None(默认),则不会跟踪任何事件。每个事件都发生在时间和状态的连续函数的零点处。每个函数都必须具有签名 event(t, y) 并返回一个浮点数。求解器将使用寻根算法在 event(t, y(t)) = 0 处找到 t 的准确值。默认情况下,将找到所有零。求解器会在每个步骤中寻找符号变化,因此如果在一个步骤中发生多个过零,则可能会错过事件。此外,每个事件函数可能具有以下属性:

终端:布尔,可选如果发生此事件是否终止集成。如果未指定,则为 False。

方向:浮动,可选的过零方向。如果方向为正,则事件只会在从负变为正时触发,反之亦然,如果方向为负。如果为 0,则任一方向都会触发事件。如果未分配,则隐式为 0。

您可以将 event.terminal = True 之类的属性分配给 Python 中的任何函数。

我最终在deriv()函数中引发了一个异常(并在循环中捕获它),因为事实证明它是在每个案例和每个 ODE 中结束集成的唯一可靠方法。

(虽然我在最初的问题中已经考虑过这个选项,但似乎没有其他办法。如果有人提出更好的选择,我很乐意将她/他的回答标记为已接受。更新:我接受了 Dmitry 的回答。)

恕我直言,一个更简单的答案是使用PyDSTool解决您的系统。它具有对准确事件查找的内置支持事实上,如果您使用基于 C 的求解器,它比 Matlab 的方法更容易设置和更高效(包括更快)。我是这个包的作者。

这里给出了使用事件的教程