使用python进行复杂平面的区域集成dblquad

Question

我的问题与Quantum Optics有关。我正在试图确定一个Fock状态的Wigner函数，它涉及在复杂平面中的集成。我的代码很慢，所以我想问你能否帮我提高速度？

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() 

注：这其中大部分被复制和修改计算器，所以不要要求我拥有...

Answer 1

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