问题似乎出在函数接近零的行为上。如果绘制函数,它看起来很平滑:
然而,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是总功率。在