有一个函数可以确定圆形孔径的夫琅禾费衍射图案的强度......( more information )
距离 x= [-3.8317 , 3.8317] 函数的积分必须约为 83.8%(如果假设 I0 为 100),当您将距离增加到 [-13.33 , 13.33] 时,它应该约为 95%。
但是当我在 python 中使用积分时,答案是错误的..我不知道我的代码出了什么问题:(
from scipy.integrate import quad
from scipy import special as sp
I0=100.0
dist=3.8317
I= quad(lambda x:( I0*((2*sp.j1(x)/x)**2)) , -dist, dist)[0]
print I
积分的结果不能大于 100 (I0) 因为这是 I0 的衍射......我不知道......可能是缩放......可能是方法! :(
最佳答案
问题似乎在于函数的行为接近于零。如果绘制函数,它看起来很平滑:
但是,scipy.integrate.quad
会提示舍入错误,这对于这条漂亮的曲线来说很奇怪。但是,该函数未定义为 0(当然,您除以零!),因此集成不顺利。
您可以使用更简单的集成方法或对您的功能做一些事情。您也可以从两侧将其积分到非常接近于零。但是,对于这些数字,在查看您的结果时,积分看起来并不正确。
但是,我想我预感到您的问题是什么。据我所知,您所展示的积分实际上是夫琅禾费衍射的强度(功率/面积)与中心距离的函数关系。如果要对某个半径内的总功率进行积分,则必须在两个维度上进行。
通过简单的区域积分规则,您应该在积分之前将您的函数乘以 2 pi r(或在您的情况下使用 x 而不是 r)。然后就变成了:
f = lambda(r): r*(sp.j1(r)/r)**2
或者
f = lambda(r): sp.j1(r)**2/r
甚至更好:
f = lambda(r): r * (sp.j0(r) + sp.jn(2,r))
最后一种形式是最好的,因为它没有任何奇点。它基于 Jaime 对原始答案的评论(请参阅此答案下方的评论!)。
(请注意,我省略了几个常量。)现在您可以将它从零积分到无穷大(没有负半径):
fullpower = quad(f, 1e-9, np.inf)[0]
然后您可以从其他一些半径积分并按完整强度进行归一化:
pwr = quad(f, 1e-9, 3.8317)[0] / fullpower
你得到 0.839(非常接近 84%)。如果您尝试更远的半径(13.33):
pwr = quad(f, 1e-9, 13.33)
这给出了 0.954。
应该注意的是,我们通过从 1e-9 而不是 0 开始积分引入了一个小误差。可以通过尝试不同的起点值来估计误差的大小。积分结果在 1e-9 和 1e-12 之间变化很小,因此它们似乎是安全的。当然,您可以使用,例如 1e-30,但是除法中可能存在数值不稳定。 (在这种情况下没有,但通常奇点在数值上是邪恶的。)
让我们仍然做一件事:
import matplotlib.pyplot as plt
import numpy as np
x = linspace(0.01, 20, 1000)
intg = np.array([ quad(f, 1e-9, xx)[0] for xx in x])
plt.plot(x, intg/fullpower)
plt.grid('on')
plt.show()
这就是我们得到的:
至少这看起来是对的,艾里斑的暗条纹清晰可见。
问题的最后一部分是什么:I0 定义了最大强度(单位可能是,例如 W/m2),而积分给出了总功率(如果强度以 W/m2 为单位,则总功率以 W 为单位) .将最大强度设置为 100 并不能保证总功率。这就是为什么计算总功率很重要的原因。
对于辐射到圆形区域的总功率,实际上存在一个封闭形式的方程:
P(x) = P0 ( 1 - J0(x)^2 - J1(x)^2 ),
其中 P0 是总功率。
关于python - python中强度函数的积分,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/24470389/