python - 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 时返回一个成功的积分。重复此操作会返回正确的数字，但 solver.successful() 方法返回 false:

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) 立即返回几乎所有 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 中如何处理上述两个方面，但如果上述任一挑战是使用解释循环而不是编译循环处理的，这将解释您观察到的运行时差异。