如何更好地使用 Cython 更快地求解微分方程？

Question

我想降低 Scipy 的 odeint 求解微分方程所需的时间。

为了练习，我使用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 例程本身？

编辑

我结合DavidW的建议在这两个文件ode.pyx，并solveODE.py只用了大约0.015秒运行这些建议的代码。

Answer 1

最简单的改变（这可能会给你带来很多好处）是使用 C 数学库sin并对cos单个数字而不是数字进行运算。调用numpy和确定它不是数组所花费的时间是相当昂贵的。

from libc.math cimport sin, cos

    # later
    -omega/Q + sin(theta) + d*cos(Omega*t)

我很想为输入分配一个类型d（在不更改界面的情况下，其他输入都无法轻松键入）：

def f(y, double t, params):

我想我也会像你在 Python 版本中那样返回一个列表。我认为使用 C 数组不会带来很多好处。