为了练习,我使用了Python科学计算中提到的示例作为模板。因为odeint需要一个函数f作为参数,所以我编写了一个静态类型的Cython版本,并希望odeint的运行时间能够显著减少。
函数f包含在名为ode.pyx的文件中,如下所示:
import numpy as np
cimport numpy as np
from libc.math cimport sin, cos
def f(y, t, params):
cdef double theta = y[0], omega = y[1]
cdef double Q = params[0], d = params[1], Omega = params[2]
cdef double derivs[2]
derivs[0] = omega
derivs[1] = -omega/Q + np.sin(theta) + d*np.cos(Omega*t)
return derivs
def fCMath(y, double t, params):
cdef double theta = y[0], omega = y[1]
cdef double Q = params[0], d = params[1], Omega = params[2]
cdef double derivs[2]
derivs[0] = omega
derivs[1] = -omega/Q + sin(theta) + d*cos(Omega*t)
return derivs
我随后创建了一个名为setup.py
的文件来编译该函数:
from distutils.core import setup
from Cython.Build import cythonize
setup(ext_modules=cythonize('ode.pyx'))
解决微分方程的脚本(还包括Python版本的f
)称为solveODE.py
,其代码如下:
import ode
import numpy as np
from scipy.integrate import odeint
import time
def f(y, t, params):
theta, omega = y
Q, d, Omega = params
derivs = [omega,
-omega/Q + np.sin(theta) + d*np.cos(Omega*t)]
return derivs
params = np.array([2.0, 1.5, 0.65])
y0 = np.array([0.0, 0.0])
t = np.arange(0., 200., 0.05)
start_time = time.time()
odeint(f, y0, t, args=(params,))
print("The Python Code took: %.6s seconds" % (time.time() - start_time))
start_time = time.time()
odeint(ode.f, y0, t, args=(params,))
print("The Cython Code took: %.6s seconds ---" % (time.time() - start_time))
start_time = time.time()
odeint(ode.fCMath, y0, t, args=(params,))
print("The Cython Code incorpoarting two of DavidW_s suggestions took: %.6s seconds ---" % (time.time() - start_time))
我接下来运行:
python setup.py build_ext --inplace
python solveODE.py
在终端中。
Python版本的时间大约为0.055秒,而Cython版本需要大约0.04秒。
有人推荐使用Cython改进我的微分方程求解尝试吗?最好不要修改odeint例程本身。
编辑
我在两个文件ode.pyx
和solveODE.py
中采纳了DavidW的建议。使用这些建议只需大约0.015秒即可运行代码。
numba
而不是cython
,但任何差异可能很小。大部分计算时间可能是在odeint
调用您的函数时发生的上下文切换。您可能会从编写自己的数值积分函数(再次使用cython或numba)以避免上下文切换中获得最佳收益。 - Aaronf
和ode.f
是 Python 对象,每次调用至少需要一次上下文切换(0-200 步长为 0.05 的 4000 次调用),否则odeint
将无法使用任何旧的自定义用户函数。我已经通过 numba 获得了 4 倍的加速,但现在我正在努力获得更多... - Aaron