python中Intensity函数的积分

Question

有一个函数可以确定圆形孔径的夫琅和费衍射图的强度......(更多信息)

距离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的衍射......我不知道......可能是缩放...可能是方法!:(

Answer 1

问题似乎是函数的行为接近于零.如果绘制函数,它看起来很平滑:

然而,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是总功率.