python的密度图使用贝塞尔积分制作衍射图案，但它不会停止运行 [英] Density plot with python making a Diffraction pattern with Bessel Integrals but it wont stop running

问题描述

我正在尝试制作圆形衍射图案，它的中心点被一系列环包围.它涉及到代码中定义的贝塞尔积分.

我的问题是它花费的时间太长，比如我等了 10 分钟让代码运行但没有显示任何内容.我知道这是因为我的贝塞尔积分每点有1000次迭代，任何人都可以帮忙吗?

我在正确的轨道上吗?

我正在尝试通过 Mark Newmans 的《计算物理学》自学自己的 Python 和计算物理学，练习是计算物理学的 5.4.这是该章节的链接.它在第9页上. 运行您的代码.PyPy擅长优化此类事情.您需要安装其numpy版本(请参见 http://pypy.readthedocs.org/en/latest/getting-started.html#).但是 PyPy 在我的机器上在 40 秒内完成了您的原始代码(删除了绘图内容).

编辑

我的系统上没有在 PyPy 中安装 plotlib，所以我在最后用

替换了你的 plt 调用

#plt.imshow(强度,vmax=0.01,cmap='hot')#plt.gray()#plt.show()np.savetxt('bessel.txt', 强度)

并创建了一个单独的程序，我使用普通 Python 执行该程序，其中包含:

将 numpy 导入为 np将 pylab 导入为 plt强度 = np.loadtxt('bessel.txt')plt.imshow(强度,vmax=0.01,cmap='hot')plt.gray()plt.show()

这会生成下面的图像，并对您的代码进行以下修改:

sepration = 0.008 # 减少分离强度 = np.empty([points,points],np.float)对于 i 在范围内(点):y =间隔*(i-points/2)#图像中心对于范围内的 j(点):x = 分离*(j - 点数/2)

I am trying to make a circular diffraction pattern, which has a central spot surrounded by a series of rings. It involves a Bessel integral to do it which is defined in the code.

My problems is that it takes too long like I waited 10 min for the code to run but didn't get anything to display. I understand it is because my Bessel integral has 1000 iterations per point can any one help with this ?

Am I on the right track?

I am trying to self teach myself python and computational physics via Mark Newmans book Computational Physics the exercise is 5.4 of Computational Physics.Here is a link to the chapter. It is on page 9. http://www-personal.umich.edu/~mejn/cp/chapters/int.pdf

Here is the Image I am trying to make.

.

My Code:

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

解决方案

Your code seems to run fine - if I reduce N down to 100 (from 1000) and the image size (points) down to 50 (from 500). I get the following image after around 4s of execution time:

Here's what we get when profiling your code with cProfile:

$ python -m cProfile -s time bessel.py | head -n 10
         361598 function calls (359660 primitive calls) in 3.470 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   252399    2.250    0.000    2.250    0.000 bessel.py:24(f)
     2499    0.821    0.000    3.079    0.001 bessel.py:23(J)
        1    0.027    0.027    3.472    3.472 bessel.py:15(<module>)
     2499    0.015    0.000    3.093    0.001 bessel.py:45(I)
        1    0.013    0.013    0.013    0.013 backend_macosx.py:1(<module>)

So it looks like most of your execution time is spent in f. You could possibly optimise this function, or alternatively try running your code using PyPy. PyPy is excellent at optimising this sort of thing. You need to install their version of numpy (see http://pypy.readthedocs.org/en/latest/getting-started.html#). But PyPy completes your original code in 40s on my machine (with the plotting stuff removed).

EDIT

I haven't got plotlib installed in PyPy on my system, so I replaced your plt calls at the end with

#plt.imshow(Intensity,vmax=0.01,cmap='hot')
#plt.gray()
#plt.show()
np.savetxt('bessel.txt', Intensity)

And created a separate program that I execute with normal Python containing:

import numpy as np
import pylab as plt

Intensity = np.loadtxt('bessel.txt')

plt.imshow(Intensity,vmax=0.01,cmap='hot')
plt.gray()
plt.show()

This produces the image below with the following modifications to your code:

sepration = 0.008 # decrease separation

Intensity = np.empty([points,points],np.float)

for i in range(points):
    y = sepration*(i - points/2) # centre on image
    for j in range(points):
        x = sepration*(j - points/2)

这篇关于python的密度图使用贝塞尔积分制作衍射图案，但它不会停止运行的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！