矢量化 SciPy ode 求解器答案

【问题标题】：Vectorized SciPy ode solver矢量化 SciPy ode 求解器
【发布时间】：2014-02-11 14:08:50
【问题描述】：

我的问题是关于当前的 scipy ode 求解器。从scipy doc page来看，它们的用法是：

# A problem to integrate and the corresponding jacobian:

from scipy.integrate import ode
y0, t0 = [1.0j, 2.0], 0
def f(t, y, arg1):
    return [1j*arg1*y[0] + y[1], -arg1*y[1]**2]
def jac(t, y, arg1):
    return [[1j*arg1, 1], [0, -arg1*2*y[1]]]

# The integration:
r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True)
r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
t1 = 10
dt = 1
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    print("%g %g" % (r.t, r.y))

我的问题是：它使用了大量的 python 循环（while 循环），这实际上会减慢程序的速度。我可以尝试编写 C 代码并使用 ctypes 使其更快，但我将无法在 scipy 中访问像 dopri5 这样的好算法（除非我自己实现它）。

是否有任何矢量化的编码方式可以加快速度？

谢谢！

【问题讨论】：

标签： python scipy ode

【解决方案1】：

这有点晚了，但我认为你需要的是scipy.integrate.odeint。它需要一组时间点并计算每个时间点的解决方案。这使得循环在 Python 之外（在 FORTRAN 中）进行计算。

【讨论】：

永远不会太晚，伙计！这就是我最终使用的，其他人在未来查看这篇文章可能会发现它很有帮助。非常感谢！

【解决方案2】：

“向量化”意味着一次并行地进行大量计算。是的，详细的实现将涉及迭代，但它在C 中，并且顺序对你来说并不重要，Python 程序员。

但是像这样的 ode 解决方案本质上是串行操作。您必须在时间t 解决问题，然后才能在时间t+dt 解决问题。您无法通过时间矢量化解决方案。您可以做的最好的事情是选择一个ode 求解器，它可以对时间步长 (dt) 做出明智的选择，尽可能大步长，在需要捕捉快速变化时小步长。

一个好的 ode 求解器可以让您对空间维度进行矢量化 - 即并行求解 10 个 ode。您还可以向量化计算雅可比矩阵，同时返回 y+dy 和 y-dy。基本上，您希望尽可能快地计算f 和jac。

【讨论】：

非常感谢！这对我来说很有意义。