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)

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

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! :(

回答1:

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

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.

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

or

f = lambda(r): sp.j1(r)**2/r

or even better:

f = lambda(r): r * (sp.j0(r) + sp.jn(2,r))

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

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)

which gives 0.954.

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

Let us do one thing still:

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

And this is what we get:

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

---

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

where P0 is the total power.

回答2:

Note that you also can get a closed form solution for your integration using Sympy:

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 evaluates to:

1600*d*besselj(0, Abs(d))**2/3 + 1600*d*besselj(1, Abs(d))**2/3 - 800*besselj(1, Abs(d))**2/(3*d)

with besselj(0,r) corresponding to sp.j0(r).

回答3:

They might be a singularity in the integration algorithm when doing the jacobian at x = 0. You can exclude this points from the integration with "points":

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

I get then the following result (is this your desired result?)

331.4990321315221