python - 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 进行科学计算的初学者，如果评论指出提高评估速度的方法，我会很高兴。是否有重写 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

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

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

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

# Do in Python## >>> import pyximport; pyximport.install(reload_support=True)# >>> import cythonmodulecimport numpy as npcimport cythoncdef extern from "complex.h":    double complex csqrt(double complex z) nogil    double complex cexp(double complex z) nogil    double creal(double complex z) nogil    double cimag(double complex z) nogilfrom libc.math cimport sqrtfrom scipy.integrate import dblquadcdef class Params:    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(True)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)    return creal(cexp(1j*p.k*(R1-R2)) * (-1j*p.z/p.lam/R2/R1**2) * (1+1j/p.k/R1))@cython.cdivision(True)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)    return 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):    return sqrt(p.ra**2 + x**2)def doit(lam, y0, k, zxp, z, ra):    p = Params(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,))    return rr + 1j*ri