11

我想使用 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 之类的工具编译这些函数是否有意义？如果是：哪种工具最适合此应用程序？

4

3 回答 3

14

通过使用 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 + (yp.y0)**2 + pz**2)
    R2 = sqrt(x**2 + y**2 + p.zxp**2)
    返回 creal(cexp(1j*pk*(R1-R2)) * (-1j*pz/p.lam/R2/R1**2) * (1+1j/pk/R1))

@cython.cdivision（真）
def integrand_imag(double x, double y, 参数 p):
    R1 = sqrt(x**2 + (yp.y0)**2 + pz**2)
    R2 = sqrt(x**2 + y**2 + p.zxp**2)
    返回 cimag(cexp(1j*pk*(R1-R2)) * (-1j*pz/p.lam/R2/R1**2) * (1+1j/pk/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
于 2013-04-03T17:37:48.530 回答
4

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

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

于 2013-04-03T16:47:40.783 回答
0

quadpy（我的一个项目）支持许多磁盘功能的集成方案。它支持复值函数并且是完全向量化的。例如，Peirce 的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)
于 2019-10-24T09:20:55.447 回答