python中强度函数的积分

Question

有一个函数可以确定圆孔夫琅和费衍射图样的强度。。。（more information）

函数在距离x=[-3.8317，3.8317]中的积分必须是83.8%（如果假设I0是100），当你把距离增加到[-13.33，13.33]时，应该是95%。 但是当我在python中使用integral时，答案是错误的。。我不知道我的代码出了什么问题：（

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（当然，您要除以0！）所以整合不好。在

您可以使用更简单的积分方法，或者对您的函数做些什么。你也可以从两边把它积分到非常接近于零的位置。然而，在这些数字下，当你看到你的结果时，积分看起来并不正确。在

不过，我想我对你的问题有预感。据我所知，你所示的积分实际上是夫琅和费衍射强度（功率/面积）与中心距离的函数。如果你想在某个半径范围内整合总能量，你必须在两个维度上进行。在

根据简单的面积积分规则，在积分之前，你应该用2πr乘以你的函数（或者用x代替r）。然后变成：

f = lambda(r): r*(sp.j1(r)/r)**2

或者

^{pr2}$

甚至更好：

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

我们得到的是：

{2美元^

至少这看起来是对的，艾里圆盘的黑色条纹清晰可见。在

---

问题的最后一部分是什么：I0定义了最大强度（单位可能是，例如W/m2），而积分给出了总功率（如果强度以W/m2为单位，则总功率以W为单位）。将最大强度设置为100并不保证总功率。这就是为什么计算总功率很重要。在

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

p（x）=P0（1-J0（x）^2-J1（x）^2）

式中，P0是总功率。在

Answer 2

当在x=0时，它们可能是积分算法中的一个奇点。您可以从与“points”的集成中排除这些点：

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

我得到以下结果（这是你想要的结果吗？）在

^{pr2}$

Answer 3

请注意，您还可以使用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计算结果为：

^{pr2}$

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