了解 scipy 集成的内部行为答案

【问题标题】：Understanding scipy integrate's internal behavior了解 scipy 集成的内部行为
【发布时间】：2020-10-02 05:58:48
【问题描述】：

我想了解scipy.integrate 在内部做什么。也就是说，似乎正在发生一些奇怪和不一致的事情。

如何让它正常工作？我需要它一次执行一个集成步骤，因为我在 ODE 中使用 t 做了一些事情，并且需要它保持一致

所以，这是我的 MWE

import numpy as np
from scipy.integrate import ode

t0 = 0
t1 = 1

def myODE(t, x):
    print('INTERNAL t = {time:2.3f}'.format(time=t))

    Dx = np.zeros([2, 1])
    Dx[0] = -x[0]**2
    Dx[1] = -x[1]**2

    return Dx

simulator = ode(myODE).set_integrator('dopri5')
simulator.set_initial_value(np.ones([2,1]), t0)

t = simulator.t

while t < t1:
    t = simulator.t
    print('Outside integrate t = {time:2.3f}'.format(time=t))
    x = simulator.integrate(2, step=True) 

    print('x1 = {x1:2.3f}'.format(x1=x[0,0]))

我试图一次执行一个集成步骤。相反，integrate 做了其他事情。从下面的输出可以看出，它一次执行多个步骤，并且这些步骤不一致：有时，t 增加然后又减少。

Outside integrate t = 0.000
INTERNAL t = 0.000
INTERNAL t = 0.010
INTERNAL t = 0.004
INTERNAL t = 0.006
INTERNAL t = 0.016
...
INTERNAL t = 1.969
INTERNAL t = 1.983
INTERNAL t = 2.000
INTERNAL t = 2.000
x1 = 0.333
Outside integrate t = 2.000
INTERNAL t = 2.000
...
INTERNAL t = 2.000
x1 = 0.333

【问题讨论】：

ode 求解器使用可变时间步长是正常的。根据特定的求解器，它可能会计算多次接近以获得数字jacobian（导数）。当函数变化快时，它也可以使用精细的步骤，而当它变化较慢时，它可能会使用更长的步骤。目标是在保持准确性的同时最小化评估的总数。文档提供了期刊参考。
看odeint。它允许您将数组指定为t 点。
您可以为dopri5 求解器提供nsteps=1 参数。然而，这仍然导致采取 2 个步骤。请注意，dopri5 每一步有 6+1 次 ODE 函数评估，可以重用一步的最后一次评估，因为它与下一步的第一步相同。
无论你想用 ODE 做什么，它几乎肯定不应该依赖于所用方法的内部及其步长控制器。使用事件和密集输出选项，如果可能的话，积极开发的接口solve_ivp。
@LutzLehmann 1solve_ivp` 不允许您一次执行一个步骤并存储中间数据。说，RK45 是一个在我看来要优越得多的类。问题是 events 不适用于我的情况：我需要在 RHS 中缓冲 t 变量

标签： scipy integration ode integrate

【解决方案1】：

所有优于标准固定步长 RK4 方法的求解器都使用可变步长。将求解器视为黑盒，无法知道使用了哪些内部步长。

然而，众所周知的是，显式一步法有多个阶段，至少等于它们的顺序，每个阶段都包含在接近但不一定在解轨迹上的点调用 ODE 函数.隐式方法的阶数可能少于阶数，但需要迭代方法来求解隐式阶跃方程。

Dormand-Prince 45 方法有 7 个阶段，其中最后一个阶段也是下一步的第一个阶段，因此长期平均每步 6 次评估。这就是您在ode(dopri) 方法中看到的内容。

INTERNAL t =   0.00000000
INTERNAL t =   0.01000000
INTERNAL t =   0.00408467
INTERNAL t =   0.00612700
INTERNAL t =   0.01633866
INTERNAL t =   0.01815407
INTERNAL t =   0.02042333
INTERNAL t =   0.02042333
INTERNAL t =   0.03516563
INTERNAL t =   0.04253677
INTERNAL t =   0.07939252
INTERNAL t =   0.08594465
INTERNAL t =   0.09413482
INTERNAL t =   0.09413482
Outside integrate t =   0.09413482
...

可以看到，scipy 方法的最小步骤由 2 个 DoPri 步骤组成。在第一步的序列中，第一次评估只是探测初始步长是否合适，这只进行一次。所有其他步骤点都在规定的时间t_n+c_i*dt 其中c=[0,1/5,3/10,4/5,8/9,1,1]。

您可以使用作为新接口solve_ivp 的步进器的新类获得正确的单步。注意这里的默认容差比ode(dopri) 的情况要宽松得多，可能遵循Matlab 的哲学，即以最小的努力生成“足够好”的图。对于RK45，这可能看起来像

simulator = RK45(myODE, t0, [1,1], t1, atol=6.8e-7, rtol=2.5e-8)

t = simulator.t

while t < t1:
    simulator.step() 
    t = simulator.t
    x = simulator.y
    print(f'Outside integrate t = {t:12.8f}')
    print(f'x1 = {x[0]:12.10f}, err = {x[0]-1/(1+t):8.6g}')

这使用了稍微不同的内部步骤，但如前所述，具有“真正的”单步输出。

INTERNAL t =   0.00000000
INTERNAL t =   0.01000000
INTERNAL t =   0.00408223
INTERNAL t =   0.00612334
INTERNAL t =   0.01632891
INTERNAL t =   0.01814323
INTERNAL t =   0.02041114
INTERNAL t =   0.02041114
Outside integrate t =   0.02041114
x1 = 0.9799971436, err = 5.2347e-13
INTERNAL t =   0.04750541
INTERNAL t =   0.06105254
INTERNAL t =   0.12878821
INTERNAL t =   0.14083011
INTERNAL t =   0.15588248
INTERNAL t =   0.15588248
Outside integrate t =   0.15588248
x1 = 0.8651399668, err = 1.13971e-07
...

如果您的输入是阶跃函数或零阶保持，最方便的解决方案是循环遍历这些阶跃并在每个阶跃初始化一个 RK45 对象，并将阶跃段作为积分边界。将最后一个值保存为下一步的初始值。或许也可以将最后一个步长作为下一步的初始步长。

在 ODE 函数中直接使用阶跃函数效率低下，因为步长控制器期望一个非常平滑的 ODE 函数来实现最优步长序列。在被严重违反的跳跃处，可能导致步长的局部明显减小，从而增加函数评估的次数。

【讨论】：