了解scipy整合的内部行为 [英] Understanding scipy integrate's internal behavior

问题描述

我正在尝试了解 scipy.integrate 在内部进行的操作.也就是说，似乎正在发生一些奇怪且不一致的事情.

I am trying to understand what scipy.integrate is doing internally. Namely, it seems that something weird and inconsistent is happening.

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

How get it working properly? I need it to perform one integration step at a time, because I do some stuff with t inside the ODE and need it to be consistent

所以，这是我的MWE

So, here is my 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 增大而减小.

What I'm trying to do to perform one integration step at a time. Instead, integrate does something else. As you can see from the output below, it performs several steps at a time, and those steps are inconsistent: sometimes, t increases and the decreases again.

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

推荐答案

所有优于标准固定步长RK4方法的求解器都使用可变步长.将求解器视为黑匣子，不知道使用了什么内部步长.

All solvers that are better than the standard fixed-step RK4 method use a variable step size. Treating the solver as a black-box, one can not know what internal step sizes are used.

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

What is known, however, is that the explicit one-step methods have multiple stages, at least equal to their order, that each comprise a call to the ODE function at a point close to, but not necessarily on the solution trajectory. Implicit methods may have less stages than the order, but require an iterative approach to the solution of the implicit step equations.

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

The Dormand-Prince 45 method has 7 stages, where the last stage is also the first of the next step, so in the long-time average 6 evaluations per step. This is what you see in the ode(dopri) method.

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] .

There one can see that the minimal step of the scipy method consists of 2 DoPri steps. In the sequence of the first step, the first evaluation is just probing if the initial step size is appropriate, this is only done once. All the other step points are at the prescribed times t_n+c_i*dt where c=[0,1/5,3/10,4/5,8/9,1,1].

使用新类作为新界面 solve_ivp 的步进器，您可以获得正确的单个步骤.请注意，默认公差要比 ode(dopri)情况宽松得多，这可能是遵循Matlab产生足够好"信号的原则.只需很少的努力就可以完成.对于 RK45 ，它看起来像

You can get proper single steps with the new classes that are the steppers for the new interface solve_ivp. Take care that the default tolerances are here much looser than in the ode(dopri) case, probably following the Matlab philosophy of generating "good enough" plots with minimal effort. For RK45 this can look like

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}')

这使用了稍微不同的内部步骤，但是，如上所述，它具有真"的含义.单步输出.

This uses slightly different internal steps, but, as said, has a "true" single-step output.

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 对象作为整合边界.将最后一个值保存为下一步的初始值.也许也将最后一个步长作为下一步的初始步长.

If you have an input that is a step function, or a zero-order hold, the most expedient solution would be to loop over the steps and initialize one RK45 object per step with the step segment as integration boundaries. Save the last value as initial value for the next step. Perhaps also the last step size as initial step size in the next step.

直接在ODE函数内部使用阶跃函数效率低下，因为阶跃大小控制器期望最优步长序列具有非常平滑的ODE函数.在严重违反的跳跃动作上，可能会导致步幅急剧缩小，从而增加功能评估的次数.

Directly using a step function inside the ODE function is inefficient, as the step size controller expects a very smooth ODE function for an optimal step size sequence. At jumps that is grossly violated and can lead to stark local reductions in the step size, and accordingly an increased number of function evaluations.

这篇关于了解scipy整合的内部行为的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！