有一个函数可以确定圆孔的 Fraunhofer 衍射图案的强度......(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 的衍射...我不知道..可能是缩放...可能是方法! :(
最佳答案
问题似乎出在接近零的函数行为中。如果绘制函数,它看起来很平滑:
但是,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 是总功率。
关于python - python中强度函数的积分,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/24470389/