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：

^{pr2}$

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

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

问题2：解算器.积分（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

由于您已经将t0与第一个调用solver.integrate(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

解算器在一次调用中执行的最大（内部）步骤数不会引发错误。这可以用nsteps参数set_integrator设置。如果一次集成大量时间，即使没有任何错误，也将超过nsteps，并引发以下错误消息：

^{pr2}$

积分器会在任何时候停止。在

问题2.2

如果设置nsteps=10**10，则集成运行时不会出现问题。不过，它还是相当快的（在我的机器上大约是1 s）。原因如下：

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

• 积分器中的向量和矩阵运算。在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)，这是一个巨大的数字。在

另一方面，整合速度之快也不足为奇。大多数提供的方法使用步长自适应，并且在标准误差公差下，这可能只导致几十个内部步骤。减小误差容限将增加内部步骤的数量，并可能显著改变数值结果。在