Why is there noise for higher bandwidth intensity distributions?

Question

I have been trying to computationally evaluate the following E-field computationally:

It is basically the sum of waves of a given wavelength, weighted by the square root of a gaussian distribution for the wavelengths.

I compute it via Python by performing gauss quadrature integrals for each value of $x$, through the package scipy.integrate.quad function. The code is stated below:

# Imports
import numpy as np
import scipy as sp
from scipy import integrate

# Parameters
mu = 0.635 # mean wavelength
sigma = 0.01 # std dev of wavelength distribution
# wl is wavelength
x_ara = np.arange(0, 1.4, 0.01)

# Limits of Integration
lower, upper = mu - 4*sigma, mu+4*sigma
if lower < 0 :
    print('lower limit met')
    lower = 1e-15 # cannot evaluate sigma = 0 due to singularity for gaussian function

# Functions
def Iprofile_func(wl, mu, sigma):
    profile = np.exp(-( ((wl-mu) / (np.sqrt(2)*sigma))**2))
    return profile

def E_func(x_ara, wl, mu, sigma):
    return np.sqrt(Iprofile_func(wl, mu, sigma)) * np.cos(2*np.pi/wl * (x_ara))

# Computation
field_ara = np.array([])
for x in x_ara:
    def E(wl):
        return E_func(x, wl, mu, sigma)

    field = sp.integrate.quad(E, lower, upper)[0]
    field_ara = np.append(field_ara, field)

I fixed the value of $\mu$ = 0.635, and performed the same computation for two values of $\sigma$, $\sigma$ = 0.01 and $\sigma$ = 0.2. The arrays that I get, I have plotted below, the plot above is the wavelength distribution, while the plot below is the computed field array:

Why does noise appear in the computed field when the value of sigma increases?

Answer 1

For large x even small changes in lambda do already lead to quickly fluctuating integrands. At some point the numeric integration routine will either take very, very long to converge or will not consider sufficiently many integration points so that the contributions from each integration point will not cancel out completely and will show exactly the noise that you see. When I run the code I actually get a warning from scipy about reaching a limit ("IntegrationWarning: The maximum number of subdivisions (50) has been achieved.").

The good thing: you know that for sufficiently large x the integration must go to zero. There is no need to compute it outside a reasonable range.

Example:

x = 10, mu = 0.635, sigma = 0.01

The integration bounds are mu+/-4sigma = [0.595, 0.675]

2Pi/0.595*10=105.6, 2Pi/0.675*10=93.08

That means roughly two oscillations of the integrand over the wavelength range at x=10.

x = 100, everything else the same

That means 20 oscillations of the integrand over the wavelength range.

x = 10, mu = 0.635, sigma = 0.1

The integration bounds are mu+/-4sigma=[0.235, 1.035]

2Pi/0.235*10=267.37, 2Pi/1.035*10=60.71

That means 33 oscillations of the the integrand over the wavelength range already at x=10.

x = 100, everything else the same

That means 329 oscillations of the the integrand over the wavelength range.

More and more integration points wold be needed if x or sigma get large. Therefore there is no alternative to increase the limits in scipy.integrate for larger x.