Python中强度函数的积分

Question

有一个函数可以用来确定圆形孔径的弗朗霍夫衍射图案的强度... (更多信息)

在距离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

在计算雅可比矩阵的时候，当x=0时，可能会出现一个特殊情况，这会影响到积分算法。你可以通过使用"points"来排除这些点，不让它们参与积分：

f = lambda x:( I0*((2*sp.j1(x)/x)**2))
I = quad(f, -dist, dist, points = [0])

这样我得到的结果是（这就是你想要的结果吗？）

331.4990321315221

Answer 2

注意，你也可以使用 Sympy 来得到你积分的一个明确解：

import sympy as sy

sy.init_printing()  # LaTeX like pretty printing in IPython

x,d = sy.symbols("x,d", real=True)

I0=100
dist=3.8317
f = I0*((2*sy.besselj(1,x)/x)**2)  # the integrand
F = f.integrate((x, -d, d))  # symbolic integration
print(F.evalf(subs={d:dist}))  # numeric evalution

F 的计算结果是：

1600*d*besselj(0, Abs(d))**2/3 + 1600*d*besselj(1, Abs(d))**2/3 - 800*besselj(1, Abs(d))**2/(3*d)

其中 besselj(0,r) 对应于 sp.j0(r)。

Answer 3

这个问题似乎出在函数在零附近的表现。如果把这个函数画出来，它看起来是平滑的：

但是，scipy.integrate.quad却抱怨出现了舍入误差，这让人很奇怪，因为这个曲线看起来很美。不过，函数在0点是没有定义的（当然，因为你在除以零！），所以积分就不太顺利。

你可以使用更简单的积分方法，或者对你的函数做点什么。你也可以尝试从两边接近零进行积分。不过，使用这些数字时，结果看起来不太对劲。

不过，我觉得我大概知道你的问题是什么。根据我的记忆，你展示的积分实际上是弗朗霍夫衍射的强度（功率/面积）与距离中心的关系。如果你想在某个半径内计算总功率，就需要在二维中进行积分。

根据简单的面积积分规则，你应该在积分之前把你的函数乘以2π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/m²），而积分给出了总功率（如果强度是W/m²，总功率就是W）。将最大强度设置为100并不能保证总功率有什么特定的值。这就是为什么计算总功率很重要。

实际上，存在一个封闭形式的方程，用于计算辐射到圆形区域的总功率：

P(x) = P0 ( 1 - J0(x)² - J1(x)² ),

其中P0是总功率。