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

Question

显然，我已经阅读了文档，但我没有能够找到一个更详细的描述，正在发生的情况下的封面。具体来说，有几种行为让我很困惑：

一般设置

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

由于您已经将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，并引发以下错误消息：

/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，则集成运行时不会出现问题。不过，它还是相当快的(在我的机器上大约有一个 )。其原因如下：

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

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

另一方面，集成速度快也就不足为奇了。大多数提供的方法都使用步长调整，而对于标准的误差公差，这可能只会导致几十个内部步骤。减少误差公差将增加内部步骤的数目，并可能大大改变数值结果。