python - python中强度函数的积分

Question

有一个函数可以确定圆孔的 Fraunhofer 衍射图案的强度......(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 的衍射...我不知道..可能是缩放...可能是方法! :(

Answer 1

问题似乎出在接近零的函数行为中。如果绘制函数，它看起来很平滑:

但是，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 是总功率。