2012-04-24 5 views
1

моя проблема связана с квантовой оптикой. Я пытаюсь определить функцию Вигнера для состояния Фока, которая включает интеграцию в комплексную плоскость. Мой код очень медленный, поэтому я хотел бы спросить, можете ли вы помочь мне улучшить скорость?Интеграция области в сложной плоскости с использованием 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() 

Примечание: Большая часть этого была скопирована и изменена с StackOverflow, так что я не претендует на право собственности ...

+0

Я предполагаю, что это будет намного проще и быстрее, по крайней мере, чтобы отобразить комплексная плоскость к двумерному вещественному пространству. – pedda

+0

'if w! = W'? шутки в сторону? Также вставляемый код нуждается в улучшенном отступе – Tobia

ответ

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

 Смежные вопросы

  • Нет связанных вопросов^_^