我想使用来自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(双复数z)nogil
double cimag(double complex z)nogil
来自libc.math cimport sqrt
来自scipy.integrate import dblquad
cdef类参数:
cdef public double lam，y0，k，zxp，z，ra
def __init __(self，lam，y0，k，zxp，z，ra):
self.lam = lam
self.y0 = y0
self.k = k
self.zxp = zxp
self.z = z
self.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
关于python - Scipy:加快2D复数积分的计算，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/15793224/