python – Scipy：加快2D复数积分的计算

概述我想从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 我想从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是事先确定的变量.要评估半径为圆的区域的积分,我使用以下命令：

from __future__ import divisionfrom scipy.integrate import dblquadfrom pylab import *def ymax(x):    return sqrt(ra**2-x**2)lam = 0.000532zxp = 5.z = 4.94k = 2*pi/lamra = 1.0res_real = dblquad(integrand_real,-ra,ra,lambda x: -ymax(x),lambda x: ymax(x))res_imag = dblquad(integrand_imag,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 loopIn [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,ra    def __init__(self,lam,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,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,ra):    p = Params(lam=lam,ra=ra)    rr,err = dblquad(integrand_real,lambda x: -ymax(x,p),lambda x: ymax(x,args=(p,))    ri,err = dblquad(integrand_imag,))    return rr + 1j*ri

总结

以上是内存溢出为你收集整理的python – Scipy：加快2D复数积分的计算全部内容，希望文章能够帮你解决python – Scipy：加快2D复数积分的计算所遇到的程序开发问题。

如果觉得内存溢出网站内容还不错，欢迎将内存溢出网站推荐给程序员好友。