1
我的問題與Quantum Optics有關。我正在試圖確定一個Fock狀態的Wigner函數,它涉及在複雜平面中的集成。我的代碼很慢,所以我想問你能否幫我提高速度?使用python進行復雜平面的區域集成dblquad
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()
注:這其中大部分被複制和修改計算器,所以不要要求我擁有...
我想這將是一個更容易或至少更快的映射覆雜的平面到二維真實空間。 – pedda 2012-05-01 15:18:05
'如果w!= w'?當真?此外,粘貼的代碼需要更好的縮進 – Tobia 2013-03-17 15:30:10