您好, 欢迎来到 !    登录 | 注册 | | 设为首页 | 收藏本站

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

5b51 2022/1/14 8:22:54 python 字数 4130 阅读 630 来源 www.jb51.cc/python

我想从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

概述

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 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,lambda x: ymax(x))
res = res_real[0]+ 1j*res_imag[0]

根据剖析器,两个被积函数被评估约35000次.总计算大约需要一秒钟,这对我应用程序来说太长了.

我是使用Python和Scipy进行科学计算的初学者,对于提出评估速度的评论意见很高兴.有没有办法重写integrand_real和integrand_complex函数中的命令,这些功能可能会导致速度提升?

使用Cython等工具编译这些函数是否有意义?如果是的:哪个工具最适合这个应用程序?

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 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,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复数积分的计算所遇到的程序开发问题。


如果您也喜欢它,动动您的小指点个赞吧

除非注明,文章均由 laddyq.com 整理发布,欢迎转载。

转载请注明:
链接:http://laddyq.com
来源:laddyq.com
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。


联系我
置顶