数值积分 Python 与 Matlab

Question

我的 python 代码运行大约需要 6.2 秒。 Matlab 代码在 0.05 秒内运行。为什么会这样，我能做些什么来加快 Python 代码的速度？ Cython 是解决方案吗？

Matlab：

function X=Test

nIter=1000000;
Step=.001;
X0=1;

X=zeros(1,nIter+1); X(1)=X0;

tic
for i=1:nIter
    X(i+1)=X(i)+Step*(X(i)^2*cos(i*Step+X(i)));
end
toc

figure(1) plot(0:nIter,X)

Python：

nIter = 1000000
Step = .001
x = np.zeros(1+nIter)
x[0] = 1
start = time.time()
for i in range(1,1+nIter):
      x[i] = x[i-1] + Step*x[i-1]**2*np.cos(Step*(i-1)+x[i-1])
end = time.time()
print(end - start)

Answer 1

如何加速你的 Python 代码

您最大的时间接收器是np.cos，它对输入格式执行多项检查。 这些对于高维输入是相关的并且通常可以忽略不计，但对于您的一维输入，这将成为瓶颈。 解决方案是使用math.cos，它只接受一维数字作为输入，因此速度更快（虽然不太灵活）。

另一个时间接收器正在多次索引x。 您可以通过更新一个状态变量并在每次迭代中只写入一次 x 来加快这一速度。

通过所有这些，您可以将速度提高大约 10 倍：

import numpy as np
from math import cos

nIter = 1000000
Step = .001
x = np.zeros(1+nIter)
state = x[0] = 1
for i in range(nIter):
    state += Step*state**2*cos(Step*i+state)
    x[i+1] = state

现在，您的主要问题是您真正最内层的循环完全发生在 Python 中，也就是说，您有很多包装操作会占用时间。 您可以通过使用 uFuncs（例如，使用 SymPy 的 ufuncify 创建）和使用 NumPy 的 accumulate 来避免这种情况：

import numpy as np
from sympy.utilities.autowrap import ufuncify
from sympy.abc import t,y
from sympy import cos

nIter = 1000000
Step = 0.001
state = x[0] = 1
f = ufuncify([y,t],y+Step*y**2*cos(t+y))

times = np.arange(0,nIter*Step,Step)
times[0] = 1
x = f.accumulate(times)

这几乎在瞬间运行。

……为什么这不是你应该担心的

如果您关心的是您的确切代码（并且仅此代码），那么无论如何您都不应该担心运行时间，因为无论哪种方式它都很短。 另一方面，如果您使用它来衡量运行时间相当长的问题的效率，那么您的示例将失败，因为它只考虑一个初始条件并且是一个非常简单的动态。

此外，您正在使用 Euler 方法，该方法不是非常有效或稳健，具体取决于您的步长。 后者（Step）在您的情况下非常低，产生的数据比您可能需要的多得多： 步长为 1 时，您可以看到正在发生的事情。

如果您想要在这种情况下进行稳健的集成，几乎总是最好使用现代自适应积分器，它可以自行调整步长，例如，这是使用原生 Python 积分器解决您的问题的方法：

from math import cos
import numpy as np
from scipy.integrate import solve_ivp

T = 1000
dt = 0.001

x = solve_ivp(
        lambda t,state: state**2*cos(t+state),
        t_span = (0,T),
        t_eval = np.arange(0,T,dt),
        y0 = [1],
        rtol = 1e-5
    ).y

这会根据容错rtol 自动将步长调整为更高的值。 它仍然返回相同数量的输出数据，但这是通过解决方案的插值。 对我来说它在 0.3 秒内运行。

如何以可扩展的方式加快速度

如果您仍然需要加快这样的速度，很可能您的导数 (f) 比您的示例复杂得多，因此它是瓶颈。 根据您的问题，您可以对其计算进行矢量化（使用 NumPy 或类似方法）。

如果你不能向量化，我写了一个module，专门通过在底层硬编码你的导数来关注这一点。 这是您的示例，采样步骤为 1。

import numpy as np
from jitcode import jitcode,y,t
from symengine import cos

T = 1000
dt = 1

ODE = jitcode([y(0)**2*cos(t+y(0))])
ODE.set_initial_value([1])
ODE.set_integrator("dop853")
x = np.hstack([ODE.integrate(t) for t in np.arange(0,T,dt)])

这会在瞬间再次运行。虽然这可能不是相关的速度提升，但它可以扩展到大型系统。

Answer 2

区别在于 jit 编译，Matlab 默认使用。让我们用Numba（一个Python jit-compiler）试试你的例子

代码

import numba as nb
import numpy as np
import time

nIter = 1000000
Step = .001

@nb.njit()
def integrate(nIter,Step):
  x = np.zeros(1+nIter)
  x[0] = 1
  for i in range(1,1+nIter):
    x[i] = x[i-1] + Step*x[i-1]**2*np.cos(Step*(i-1)+x[i-1])
  return x

#Avoid measuring the compilation time,
#this would be also recommendable for Matlab to have a fair comparison
res=integrate(nIter,Step)

start = time.time()
for i in range(100):
  res=integrate(nIter,Step)

end=time.time()
print((end - start)/100)

这导致每次调用的运行时间为 0.022 秒。