问python中强度函数的积分

EN

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

​

​

然而，scipy.integrate.quad抱怨舍入错误，这是非常奇怪的这条美丽的曲线.但是，函数没有定义为0(当然，您要除以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()

这就是我们得到的

​

​

至少这看起来是正确的，黑色条纹的Airy磁盘是清楚可见的。

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

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

P(x) = P0 (1- J0(x)^2 - J1(x)^2 )，

其中P0是总功率。

请注意，您还可以使用西米获得用于集成的封闭表单解决方案。

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时执行jacobian时，它们可能是积分算法中的一个奇点。您可以从与“积分”的集成中排除这一点：

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

然后我得到以下结果(这是您想要的结果吗？)

331.4990321315221