使用 python dblquad 进行复平面中的区域积分

Question

我的问题与量子光学有关。我正在尝试确定福克状态的维格纳函数，这涉及复平面中的积分。我的代码很慢，所以我想问你是否可以帮助我提高速度？

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import pylab
import numpy as np
import mpmath
import sys
mpmath.dps = 5
from sympy import laguerre_l
import scipy
from scipy.integrate import quad,dblquad

def complex_quadrature(func, a, b, **kwargs):
    def real_func(x,y):
    return scipy.real(func(complex(x,y)))
def imag_func(x,y):
    return scipy.imag(func(complex(x,y)))
real_integral = dblquad(real_func, a, b, lambda x: a, lambda x: b, **kwargs)
imag_integral = dblquad(imag_func, a, b, lambda x: a, lambda x: b, **kwargs)
return complex(real_integral[0], imag_integral[0])

f = lambda z: 1/(np.pi**2)*complex_quadrature(lambda x:scipy.exp(x.conjugate()*z-        x*z.conjugate())*complex(laguerre_l(1,0,abs(x)**2),0)*scipy.exp(-abs(x)**2), -5, 5)

fig = pylab.figure()
ax = Axes3D(fig)
X = np.arange(-4, 4, 0.1)
Y = np.arange(-4, 4, 0.1)
X, Y = np.meshgrid(X, Y)
xn, yn = X.shape
W = X*0

for xk in range(xn):
for yk in range(yn):
    try:
        z = complex(X[xk,yk],Y[xk,yk])
        print (f(z))
        w = abs(z)
                if w != w:
            raise ValueError
            W[xk,yk] = w
        except (ValueError, TypeError, ZeroDivisionError):
                pass

ax.plot_surface(X, Y, W, rstride=1, cstride=1, cmap=cm.jet,linewidth=0,alpha=.9)
pylab.show()

注意：其中大部分内容是从 stackoverflow 复制和修改的，所以我不声明所有权...

Answer 1

z = complex(X[xk,yk],Y[xk,yk])
##  print (f(z))    ## <---- Commenting this line will solve your problem
w = abs(z)