scipy.integrate.ode.integrate() 是如何工作的？

Question

我显然已经通读了 documentation，但我无法找到关于幕后发生的事情的更详细描述。具体来说，有几个行为我很困惑：

一般设置

import numpy as np
from scipy.integrate import ode

#Constants in ODE
N = 30                                 
K = 0.5                               
w = np.random.normal(np.pi, 0.1, N)

#Integration parameters
y0 = np.linspace(0, 2*np.pi, N, endpoint=False)   
t0 = 0                                                                         

#Set up the solver
solver = ode(lambda t,y: w + K/N*np.sum( np.sin( y - y.reshape(N,1) ), axis=1))
solver.set_integrator('vode', method='bdf')
solver.set_initial_value(y0, t0)

问题 1：`solver.integrate(t0)` 失败

正在设置积分器，并在 t0 第一次询问值 returns 成功积分。重复此 returns 正确的数字，但 solver.successful() 方法 returns 错误：

solver.integrate(t0)
>>> array([ 0.        ,  0.20943951,  0.41887902, ...,  5.65486678,
            5.86430629,  6.0737458 ])

solver.successful()
>>> True

solver.integrate(t0)
>>> array([ 0.        ,  0.20943951,  0.41887902, ...,  5.65486678,
            5.86430629,  6.0737458 ])

solver.successful()
>>> False

我的问题是，solver.integrate(t) 方法中发生了什么导致它第一次成功，然后失败，“不成功”的集成是什么意思？此外，为什么积分器默默地失败，并继续产生看起来有用的输出，直到我明确询问它是否成功？

相关，有没有办法重置失败的集成，或者我需要从头开始重新实例化求解器？

问题 2：solver.integrate(t) 立即 returns 几乎任何 `t`

值的答案

即使我的 y0 初始值在 t0=0 给出，我也可以请求 t=10000 的值并立即得到答案。我希望在如此大的时间跨度内进行数值积分至少需要几秒钟（例如，在 Matlab 中，要求积分超过 10000 个时间步长将需要几分钟）。

例如，重新运行上面的设置并执行：

solver.integrate(10000)
>>> array([ 2153.90803383,  2153.63023706,  2153.60964064, ...,  2160.00982959,
            2159.90446056,  2159.82900895])

Python真的有那么快吗，还是这个输出完全是胡说八道？

Answer 1

问题 0

不要忽略错误信息。是的，ode 的错误消息有时可能很含糊，但您仍然希望避免它们。

问题 1

由于您已经在 solver.integrate(t0) 的第一次调用中集成了 t0，因此您在第二次调用中集成了 0 的时间步长。这会引发神秘错误：

 DVODE--  ISTATE (=I1) .gt. 1 but DVODE not initialized      
      In above message,  I1 =         2
/usr/lib/python3/dist-packages/scipy/integrate/_ode.py:869: UserWarning: vode: Illegal input detected. (See printed message.)
  'Unexpected istate=%s' % istate))

问题2.1

解算器在一次调用中不抛出错误的情况下要执行的（内部）步数有上限。这可以用 set_integrator 的 nsteps 参数设置。如果一次积分很大，就算没问题也会超出nsteps，抛出如下错误信息：

/usr/lib/python3/dist-packages/scipy/integrate/_ode.py:869: UserWarning: vode: Excess work done on this call. (Perhaps wrong MF.)
  'Unexpected istate=%s' % istate))

每当发生这种情况时，积分器就会停止。

问题2.2

如果您设置 nsteps=10**10，集成运行没有问题。不过它仍然非常快（在我的机器上大约为 1 秒）。原因如下：

对于像您这样的多维系统，集成时有两个主要的运行时接收器：

积分器内的矢量和矩阵运算。在scipy.ode中，这些都是用NumPy运算或移植的Fortran或C代码实现的。无论如何，它们是通过编译代码实现的，没有 Python 开销，因此非常高效。
评估导数（lambda t,y: w + K/N*np.sum( np.sin( y - y.reshape(N,1) ), axis=1) 在你的例子中）。您通过 NumPy 操作实现了这一点，这些操作同样通过编译代码实现并且非常高效。您可以使用纯编译函数对此进行一点改进，但这最多会给您一个小因素。如果您改用 Python 列表和循环，速度会非常慢。

因此，对于您的问题，所有相关的事情都是在幕后用编译代码处理的，并且集成处理的效率与纯 C 程序相当。我不知道上述两个方面在 Matlab 中是如何处理的，但是如果上述任何一个挑战是使用解释循环而不是编译循环来处理的，这就可以解释您观察到的运行时差异。

Answer 2

对于第二个问题，是的，输出可能是废话。局部误差，无论是来自离散化还是浮点运算，都会随着复合因子累积，该复合因子大约是 ODE 函数的 Lipschitz 常数。初步估计，这里的 Lipschitz 常数是 K=0.5。早期误差的放大率，即它们作为全局误差的一部分的系数，因此可以大到exp(0.5*10000)，这是一个巨大的数字。

另一方面，集成速度很快也就不足为奇了。大多数提供的方法都使用步长自适应，并且使用标准误差容限，这可能只会导致几十个内部步骤。降低误差容限会增加内部步骤的数量，并且可能会显着改变数值结果。

scipy.integrate.ode.integrate() 是如何工作的？

How does scipy.integrate.ode.integrate() work?

python

scipy

ode