python中强度函数的积分 [英] Integral of Intensity function in python

问题描述

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

There is a function which determine the intensity of the Fraunhofer diffraction pattern of a circular aperture... (more information)

距离 x= [-3.8317 , 3.8317] 函数的积分必须约为 83.8%(如果假设 I0 为 100)，当您将距离增加到 [-13.33 , 13.33] 时，它应该约为 95%.但是当我在python中使用积分时，答案是错误的..我不知道我的代码出了什么问题:(

Integral of the function in distance x= [-3.8317 , 3.8317] must be about 83.8% ( If assume that I0 is 100) and when you increase the distance to [-13.33 , 13.33] it should be about 95%. But when I use integral in python, the answer is wrong.. I don't know what's going wrong in my code :(

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

Result of the integral can't be bigger than 100 (I0) because this is the diffraction of I0 ... I don't know.. may be scaling... may be the method! :(

推荐答案

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

The problem seems to be in the function's behaviour near zero. If the function is plotted, it looks smooth:

然而，scipy.integrate.quad 抱怨四舍五入错误，这对于这条美丽的曲线来说很奇怪.但是，该函数未定义为 0(当然，您是在除以零！)，因此集成不顺利.

However, scipy.integrate.quad complains about round-off errors, which is very strange with this beautiful curve. However, the function is not defined at 0 (of course, you are dividing by zero!), hence the integration does not go well.

您可以使用更简单的集成方法或对您的函数做一些事情.您也可以从两侧将其积分到非常接近于零.但是，对于这些数字，在查看您的结果时，积分看起来并不正确.

You may use a simpler integration method or do something about your function. You may also be able to integrate it to very close to zero from both sides. However, with these numbers the integral does not look right when looking at your results.

但是，我想我预感到您的问题是什么.据我所知，您所展示的积分实际上是夫琅禾费衍射的强度(功率/面积)与中心距离的函数关系.如果要对某个半径内的总功率进行积分，则必须在两个维度上进行.

However, I think I have a hunch of what your problem is. As far as I remember, the integral you have shown is actually the intensity (power/area) of Fraunhofer diffraction as a function of distance from the center. If you want to integrate the total power within some radius, you will have to do it in two dimensions.

根据简单的面积积分规则，您应该在积分之前将您的函数乘以 2 pi r(或在您的情况下使用 x 而不是 r).然后变成:

By simple area integration rules you should multiply your function by 2 pi r before integrating (or x instead of r in your case). Then it becomes:

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 对原始答案的评论(请参阅此答案下方的评论！).

The last form is best as it does not suffer from any singularities. It is based on Jaime's comment to the original answer (see the comment below this answer!).

(请注意，我省略了几个常量.)现在您可以将其从零积分到无穷大(无负半径):

(Note that I omitted a couple of constants.) Now you can integrate it from zero to infinity (no negative radii):

fullpower = quad(f, 1e-9, np.inf)[0]

然后您可以从其他半径积分并按完整强度归一化:

Then you can integrate from some other radius and normalize by the full intensity:

pwr = quad(f, 1e-9, 3.8317)[0] / fullpower

你得到 0.839(非常接近 84%).如果您尝试更远的半径(13.33):

And you get 0.839 (which is quite close to 84 %). If you try the farther radius (13.33):

pwr = quad(f, 1e-9, 13.33)

得到 0.954.

需要注意的是，我们通过从 1e-9 而不是 0 开始积分引入了一个小误差.可以通过尝试不同的起点值来估计误差的大小.积分结果在 1e-9 和 1e-12 之间变化很小，因此它们似乎是安全的.当然，您可以使用，例如 1e-30，但是除法中可能存在数值不稳定.(在这种情况下没有，但通常奇点在数值上是邪恶的.)

It should be noted that we introduce a small error by starting the integration from 1e-9 instead of 0. The magnitude of the error can be estimated by trying different values for the starting point. The integration result changes very little between 1e-9 and 1e-12, so they seem to be safe. Of course, you could use, e.g., 1e-30, but then there may be numerical instability in the division. (In this case there isn't, but in general singularities are numerically evil.)

让我们再做一件事:

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()

这就是我们得到的:

至少这看起来是对的，艾里斑的暗条纹清晰可见.

At least this looks right, the dark fringes of the Airy disk are clearly visible.

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

What comes to the last part of the question: I0 defines the maximum intensity (the units may be, e.g. W/m2), whereas the integral gives total power (if the intensity is in W/m2, the total power is in W). Setting the maximum intensity to 100 does not guarantee anything about the total power. That is why it is important to calculate the total power.

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

There actually exists a closed form equation for the total power radiated onto a circular area:

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

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

其中 P0 是总功率.

这篇关于python中强度函数的积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！