python强度公式计算_python中强度函数的积分

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

然而，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是总功率。在