Scipy：加快2D复数积分的计算 [英] Scipy: Speeding up calculation of a 2D complex integral

问题描述

我想使用来自scipy.integrate的dblquad重复计算二维复数积分。由于评估的数量将非常多，我想提高代码的评估速度。

I want to repeatedly calculate a two-dimensional complex integral using dblquad from scipy.integrate. As the number of evaluations will be quite high I would like to increase the evaluation speed of my code.

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

Dblquad does not seem to be able to handle complex integrands. Thus, I have split the complex integrand into a real and an imaginary part:

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的圆的面积上的积分，请使用以下命令：

y0, z, zxp, k, and lam are variables defind in advance. To evaluate the integral over the area of a circle with radius ra I use the following commands:

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次。整个计算大约需要一秒钟，对于我所考虑的应用程序来说太长了。

According to the profiler the two integrands are evaluated about 35000 times. The total calculation takes about one second, which is too long for the application I have in mind.

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

I am a beginner to scientific computing with Python and Scipy and would be happy about comments that point out ways of improving the evaluation speed. Are there ways of rewriting the commands in the integrand_real and integrand_complex functions that could lead to siginficant speed improvements?

使用Cython之类的工具编译这些功能是否有意义？如果是，哪种工具最适合此应用？

Would it make sense to compile those functions using tools like Cython? If yes: Which tool would best fit this application?

推荐答案

使用Cython可以使速度提高约10倍，参见以下内容：

You can gain a factor of about 10 in speed by using Cython, see below:

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中。）

This is probably not enough, and the bad news is that this is probably quite close (maybe factor of 2 at most) to everything-in-C/Fortran speed --- if using the same algorithm for adaptive integration. (scipy.integrate.quad itself is already in Fortran.)

要想进一步了解，您需要考虑采用其他方法来制作
整合。这需要一些思考---现在不能从
拿出多少钱。

To get further, you'd need to consider different ways to do the integration. This requires some thinking --- can't offer much from the top of my head now.

或者，您可以降低容忍度，直到

Alternatively, you can reduce the tolerance up to which the integral is evaluated.

# Do in Python
#
# >>> import pyximport; pyximport.install(reload_support=True)
# >>> import cythonmodule

cimport numpy as np
cimport cython

cdef 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) nogil

from libc.math cimport sqrt

from scipy.integrate import dblquad

cdef 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

这篇关于Scipy：加快2D复数积分的计算的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！