[FIXED] Oscillatory integral in python

Issue

I wrote the following code to plot the intensity of light exiting an optical components, which is basically a spherical Fourier integral on the incident field, so it has a Bessel function. The argument of which depends on the integrating variable (x) and the plotting variable (r).

from sympy import *
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad
from scipy.special import jn

#constants
mm = 1
um = 1e-3 * mm
nm = 1e-6 * mm

wavelength = 532*nm
klaser = 2*np.pi / wavelength
waist = 3.2*mm
angle = 2 #degrees
focus = 125 * mm
ng = 1.5 # refractive index of axicon
upperintegration = 5

#integrals

def b(angle):
    radians = angle* np.pi/180
    return klaser * (ng-1) * np.tan(radians)

def delta(angle):
    return np.pi/(b(angle)*waist)    

def integrand(x, r): 
    return klaser/focus * waist**2 * np.exp(-x**2) * np.exp(-np.pi * 1j * x/delta(angle)) * jn(0, waist*klaser*r*x/focus) * x

def intensity1D(r):
    return np.sqrt(quad(lambda x: np.real(integrand(x, r)), 0, upperintegration)[0]**2 + quad(lambda x: np.imag(integrand(x, r)), 0, upperintegration)[0]**2)  

fig = plt.figure()
ax = fig.add_subplot(111)

t = np.linspace(-3.5, 3.5, 25)
plt.plot(t, np.vectorize(intensity1D)(t))

The issue is that the plot changes drastically as I change the number of points I am using in my linspace, when I plot it.

I suspect this may be because of the oscillatory nature of the integral, so the step-size taken can dramatically change the value of the exponent and hence of the integral.

How does quad deal with this? Are there better methods to integrate numerically for this particular application?

Solution

In the call to quad, set the limit argument to a large number. This increases the maximum number subintervals that quad is allowed to use to estimate the integral. When I use

def intensity1D(r):
    re = quad(lambda x: np.real(integrand(x, r)), 0, upperintegration, limit=8000)[0]
    im = quad(lambda x: np.imag(integrand(x, r)), 0, upperintegration, limit=8000)[0]
    return np.sqrt(re**2 + im**2)

and compute the function with the array t defined as

t = np.linspace(1.5, 3, 1000)

I get the following plot:

(I also removed the line from sympy import *. sympy does not appear to be used in your script.)

You should always check the error estimate that is the second return value of quad. For example:

In [14]: r = 3.0

In [15]: val, err = quad(lambda x: np.real(integrand(x, r)), 0, upperintegration, limit=8000)

In [16]: val
Out[16]: 2.975500141416676e-11

In [17]: err
Out[17]: 1.4590630152807049e-08

As you can see, the error estimate is much larger than the approximate integral. The estimates returned by quad might be conservative, but a result with such a large error estimate should still be treated with caution. Let's take a look at the corresponding imaginary part:

In [25]: val, err = quad(lambda x: np.imag(integrand(x, r)), 0, upperintegration, limit=8000)

In [26]: val
Out[26]: 0.0026492702707317257

In [27]: err
Out[27]: 1.4808416189183e-08

val is now orders of magnitude larger than the estimated error. So when the magnitude of the complex value is computed in intensity1D(), we end up with estimated relative error on the order of 1e-5. That may be sufficient for your calculation.

At the peak near r=2.1825, the magnitude of the error estimate is still small, and it is much smaller than the computed integral:

In [32]: r = 2.1825

In [33]: quad(lambda x: np.real(integrand(x, r)), 0, upperintegration, limit=8000)
Out[33]: (6.435730031424414, 8.801375195176556e-08)

In [34]: quad(lambda x: np.imag(integrand(x, r)), 0, upperintegration, limit=8000)
Out[34]: (-6.583055286038913, 9.211333259956749e-08)

Answered By - Warren Weckesser