Почему существует шум для более высоких распределений интенсивности полосы пропускания?

Question

Я пытался вычислительно оценить следующее E-поле в вычислительном отношении:

Это в основном сумма волн данной длины волны,взвешенный по квадратному корню из гауссовского распределения для длин волн.

Я вычисляю его через Python, выполняя квадратурные интегралы Гаусса для каждого значения $ x $, через пакетную функцию scipy.integrate.quad.Код приведен ниже:

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

Я установил значение $ \ mu $ = 0.635 и выполнил одно и то же вычисление для двух значений $ \ sigma $, $ \ sigma $ = 0.01 и $ \сигма $ = 0,2.Массивы, которые я получаю, я изобразил ниже, график выше - это распределение длины волны, а график ниже - это вычисленный массив полей:

Почему в вычисляемом поле появляется шум, когда значение сигмы увеличивается?

Answer 1

Для больших x даже небольшие изменения в лямбде уже приводят к быстрым колебаниям подынтегральных выражений. В какой-то момент процедура числовой интеграции либо сойдет очень, очень долго, либо не будет учитывать достаточно много точек интеграции, так что вклады от каждой точки интеграции не будут полностью исключены и покажут именно тот шум, который вы видите. Когда я запускаю код, я на самом деле получаю предупреждение от scipy о достижении предела («IntegrationWarning: максимальное количество подразделений (50) было достигнуто».).

Хорошая вещь: вы знаете, что для достаточно большого x интеграция должна идти до нуля. Нет необходимости вычислять его за пределами разумного диапазона.

Пример:

х = 10, мю = 0,635, сигма = 0,01

Границы интегрирования: mu +/- 4sigma = [0,595, 0,675]

2Pi / 0,595 * 10 = 105,6, 2Pi / 0,675 * 10 = 93,08

Это означает примерно два колебания подынтегрального выражения в диапазоне длин волн при x = 10.

х = 100, все остальное тоже самое

Это означает 20 колебаний подынтегральной функции в диапазоне длин волн.

х = 10, мю = 0,635, сигма = 0,1

Границы интегрирования: mu +/- 4sigma = [0,235, 1,035]

2Pi / 0,235 * 10 = 267,37, 2Pi / 1,035 * 10 = 60,71

Это означает 33 колебания подынтегральной функции в диапазоне длин волн уже при x = 10.

х = 100, все остальное тоже самое

Это означает 329 колебаний подынтегральной функции в диапазоне длин волн.

Требуется все больше и больше точек интеграции, если x или sigma становятся большими. Поэтому нет альтернативы увеличению пределов в scipy.integrate для больших x.