Scipy：加快二维复积分的计算

Question

我想使用 scipy.integrate 中的 dblquad 重复计算二维复积分。由于评估的数量会很高，我想提高我的代码的评估速度。

Dblquad 似乎无法处理复杂的被积函数。因此，我将复被积函数拆分为实部和虚部：

def integrand_real(x, y):
    R1=sqrt(x**2 + (y-y0)**2 + z**2)
    R2=sqrt(x**2 + y**2 + zxp**2)
    return real(exp(1j*k*(R1-R2)) * (-1j*z/lam/R2/R1**2) * (1+1j/k/R1))

def integrand_imag(x,y):
    R1=sqrt(x**2 + (y-y0)**2 + z**2)
    R2=sqrt(x**2 + y**2 + zxp**2)
    return imag(exp(1j*k*(R1-R2)) * (-1j*z/lam/R2/R1**2) * (1+1j/k/R1))

y0、z、zxp、k、lam是预先定义的变量。为了评估半径为 ra 的圆的面积上的积分，我使用以下命令：

from __future__ import division
from scipy.integrate import dblquad
from pylab import *

def ymax(x):
    return sqrt(ra**2-x**2)

lam = 0.000532
zxp = 5.
z = 4.94
k = 2*pi/lam
ra = 1.0

res_real = dblquad(integrand_real, -ra, ra, lambda x: -ymax(x), lambda x: ymax(x))
res_imag = dblquad(integrand_imag, -ra, ra, lambda x: -ymax(x), lambda x: ymax(x))
res = res_real[0]+ 1j*res_imag[0]

根据分析器，这两个被积函数被评估了大约 35000 次。总计算大约需要一秒钟，这对于我所考虑的应用程序来说太长了。

我是使用 Python 和 Scipy 进行科学计算的初学者，我很高兴 cmets 指出了提高评估速度的方法。是否有方法可以重写 integrand_real 和 integrand_complex 函数中的命令，从而显着提高速度？

使用 Cython 之类的工具编译这些函数是否有意义？如果是：哪种工具最适合此应用程序？

Answer 1

使用 Cython 可以将速度提高大约 10 倍，见下文：

In [87]: %timeit cythonmodule.doit(lam=lam, y0=y0, zxp=zxp, z=z, k=k, ra=ra)
1 loops, best of 3: 501 ms per loop
In [85]: %timeit doit()
1 loops, best of 3: 4.97 s per loop

这可能还不够，坏消息是这可能是 与 C/Fortran 中的一切速度非常接近（最多可能是 2 倍） --- 如果使用相同的算法进行自适应积分。 （scipy.integrate.quad 本身已经在 Fortran 中。）

为了更进一步，您需要考虑不同的方法来完成 一体化。这需要一些思考——不能提供太多 现在是我的头顶。

或者，您可以降低积分的公差 被评估。

# 在 Python 中执行 # # >>> 导入 pyximport; pyximport.install(reload_support=True) # >>> 导入 cython 模块 cimport numpy 作为 np cimport cython 来自“complex.h”的 cdef extern： 双复 csqrt(双复 z) nogil 双复数 cexp(双复数 z) nogil 双Creal（双复数z）诺吉尔 双 cimag（双复数 z）诺吉尔 从 libc.math cimport sqrt 从 scipy.integrate 导入 dblquad cdef 类参数： cdef public double lam, y0, k, zxp, z, ra def __init__(self, lam, y0, k, zxp, z, ra): self.lam = 拉姆 自我.y0 = y0 自我.k = k 自我.zxp = zxp 自我.z = z 自我.ra = ra @cython.cdivision（真） def integrand_real（双 x，双 y，参数 p）： R1 = sqrt(x**2 + (y-p.y0)**2 + p.z**2) R2 = sqrt(x**2 + y**2 + p.zxp**2) 返回 creal(cexp(1j*p.k*(R1-R2)) * (-1j*p.z/p.lam/R2/R1**2) * (1+1j/p.k/R1)) @cython.cdivision（真） def integrand_imag(double x, double y, 参数 p): R1 = sqrt(x**2 + (y-p.y0)**2 + p.z**2) R2 = sqrt(x**2 + y**2 + p.zxp**2) 返回 cimag(cexp(1j*p.k*(R1-R2)) * (-1j*p.z/p.lam/R2/R1**2) * (1+1j/p.k/R1)) def ymax（双 x，参数 p）： 返回 sqrt(p.ra**2 + x**2) def doit(lam, y0, k, zxp, z, ra): p = 参数(lam=lam, y0=y0, k=k, zxp=zxp, z=z, ra=ra) rr, err = dblquad(integrand_real, -ra, ra, lambda x: -ymax(x, p), lambda x: ymax(x, p), args=(p,)) ri, err = dblquad(integrand_imag, -ra, ra, lambda x: -ymax(x, p), lambda x: ymax(x, p), args=(p,)) 返回 rr + 1j*ri

Answer 2

您是否考虑过多处理（多线程）？似乎您不需要进行最终集成（在整个集合中），因此简单的并行处理可能就是答案。即使您确实必须集成，您也可以等待正在运行的线程完成计算，然后再进行最终集成。也就是说，您可以阻塞主线程，直到所有工作人员都完成为止。

http://docs.python.org/2/library/multiprocessing.html

Answer 3

quadpy（我的一个项目）支持许多磁盘功能的集成方案。它支持复值函数并且是完全向量化的。例如 83 号订单的Peirce's scheme：

from numpy import sqrt, pi, exp
import quadpy

lam = 0.000532
zxp = 5.0
z = 4.94
k = 2 * pi / lam
ra = 1.0
y0 = 0.0


def f(X):
    x, y = X
    R1 = sqrt(x ** 2 + (y - y0) ** 2 + z ** 2)
    R2 = sqrt(x ** 2 + y ** 2 + zxp ** 2)
    return exp(1j * k * (R1 - R2)) * (-1j * z / lam / R2 / R1 ** 2) * (1 + 1j / k / R1)


scheme = quadpy.disk.peirce_1957(20)
val = scheme.integrate(f, [0.0, 0.0], ra)

print(val)

(18.57485726096671+9.619636385589759j)