我试着做一个圆形衍射图, 有一个中心点被一系列的环包围着。 它涉及到一个贝塞尔积分来做,这是在代码中定义的。在
我的问题是,我花了太长时间等待代码运行,但没有得到任何显示。我理解这是因为我的贝塞尔积分每点有1000次迭代,谁能帮上忙?在
我走对了吗?在
我试图通过马克·纽曼的书《计算物理》自学python和计算物理学练习是计算的5.4物理。这里是本章的链接。在第9页。 http://www-personal.umich.edu/~mejn/cp/chapters/int.pdf
这是我要做的形象。在
。在
我的代码:
import numpy as np
import pylab as plt
import math as mathy
#N = number of slicices
#h = b-a/N
def J(m,x): #Bessel Integral
def f(theta):
return (1/mathy.pi)*mathy.cos(m*theta - x*mathy.sin(theta)) #I replaced np. with mathy. for this line
N = 1000
A = 0
a=0
b=mathy.pi
h = ((b-a)/N)
for k in range(1,N,2):
A += 4*f(a + h*k)
for k in range(2,N,2):
A += 2*f(a + h*k)
Idx = (h/3)*(A + f(a)+f(b))
return Idx
def I(lmda,r): #Intensity
k = (mathy.pi/lmda)
return ((J(1,k*r))/(k*r))**2
wavelength = .5 # microm meters
I0 = 1
points = 500
sepration = 0.2
Intensity = np.empty([points,points],np.float)
for i in range(points):
y = sepration*i
for j in range(points):
x = sepration*j
r = np.sqrt((x)**2+(y)**2)
if r < 0.000000000001:
Intensity[i,j]= 0.5 #this is the lim as r -> 0, I -> 0.5
else:
Intensity[i,j] = I0*I(wavelength,r)
plt.imshow(Intensity,vmax=0.01,cmap='hot')
plt.gray()
plt.show()
如果我将
N
减少到100(从1000)并将图像大小(points
)减少到50(从500)的话,您的代码似乎运行得很好。在执行大约4s之后,我得到了以下图像:以下是使用cProfile分析代码时得到的结果:
因此,您的大部分执行时间似乎都花在
f
中。您可以优化此函数,或者尝试使用PyPy运行代码。PyPy在优化这类事情上很在行。您需要安装他们的numpy版本(请参见http://pypy.readthedocs.org/en/latest/getting-started.html#)。但是PyPy在我的机器上用40秒就完成了你的原始代码(去掉了绘图的东西)。在编辑
我的系统上没有在PyPy中安装plotlib,所以最后我用
^{pr2}$并创建了一个我用普通Python执行的单独程序,其中包含:
这将生成下面的图像,并对代码进行以下修改:
相关问题 更多 >
编程相关推荐