我想使用来自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进行科学计算的初学者，并且对指出提高评估速度方法的评论感到高兴。是否有方法可以重写integrand_real和integrand_complex函数中的命令，从而可以显着提高速度？
使用Cython之类的工具编译这些功能是否有意义？如果是：哪种工具最适合此应用程序？
通过使用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中执行＃＃>>> import pyximport; pyximport.install（reload_support = True）＃>>>导入cythonmodule将numpy导入为np导入cython来自“ complex.h”的cdef extern：双复数csqrt（双复数z）nogil双复数cexp（双复数z）nogil双creal（double complex z）nogildouble cimag（double complex z）nogil来自libc.math cimport sqrt来自scipy.integrate import dblquadcdef类参数：cdef public double lam，y0，k，zxp，z，radef __init __（self，lam，y0，k，zxp，z，ra）：self.lam = lamself.y0 = y0self.k = kself.zxp = zxpself.z = zself.ra = ra@ cython.cdivision（真实）def integrand_real（double x，double y，Params 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，Params 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（double x，Params 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
您是否考虑过多处理（多线程）？似乎您不需要进行最终集成（在整个集合中），因此简单的并行处理可能是答案。即使必须集成，也可以等待运行线程完成计算之后再进行最终集成。也就是说，您可以阻塞主线程，直到所有工作人员都完成为止。
quadpy（我的一个项目）支持磁盘上功能的许多集成方案。它支持复数值函数，并且已完全向量化。例如，以Peirce's scheme的顺序为83：
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)