有一个函数可以确定圆形孔径的夫琅和费衍射图案的强度... (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的衍射……我不知道……可能是缩放……可能是方法！ :(
【问题讨论】：
标签： python scipy integral numerical-integration
问题似乎在于函数的行为接近于零。如果函数被绘制出来，它看起来很平滑：
然而，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是总功率。
【讨论】：
2 * j1(x) / x = j0(x) + j2(x)，所以你也可以通过使你的函数 I0 * (sp.j0(x) + sp.jn(2, x))**2 来解决所有的数值问题。我相信这应该让你从 0 开始整合。
请注意，您还可以使用 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)。
【讨论】：
在 x = 0 处进行雅可比时，它们可能是积分算法中的奇点。您可以使用“points”从积分中排除这些点：
f = lambda x:( I0*((2*sp.j1(x)/x)**2))
I = quad(f, -dist, dist, points = [0])
然后我得到以下结果（这是您想要的结果吗？）
331.4990321315221
【讨论】：