Сравнение методов численного интегрирования в задаче дифракции плоской электромагнитной волны на прямоугольном отверстии Текст научной статьи по специальности «Физика»

Сравнение методов численного интегрирования в задаче дифракции плоской электромагнитной волны на прямоугольном отверстии

А.С. Мокеев '■2, В.М. Ямщиков1,2

1 Федеральное государственное унитарное предприятие

«Российский федеральный ядерный центр Всероссийский научно-исследовательский институт экспериментальной физики»

(ФГУП «РФЯЦ - ВНИИЭФ»);

2 Институт лазерно-физических исследований (ИЛФИ), 607188, Нижегородская область, г. Саров, пр. Мира, д. 37

Аннотация

На примере классической задачи дифракции плоской электромагнитной волны на прямоугольном отверстии рассмотрены особенности вычисления интеграла Гюйгенса-Френеля в дальней зоне стандартными квадратурными методами численного интегрирования и специализированным методом коллокаций Левина. Для квадратурных методов численного интегрирования получен критерий оценки шага интегрирования в зависимости от размеров области наблюдения на экране и требуемой точности вычислений. Показаны преимущества использования специализированного метода коллокаций Левина над стандартными методами численного интегрирования.

Ключевые слова: дифракционный интеграл, интегрирование осциллирующих функций, метод прямоугольников, метод трапеций, метод Левина, дифракция в дальней зоне.

Цитирование: Мокеев, А.С. Сравнение методов численного интегрирования в задаче дифракции плоской электромагнитной волны на прямоугольном отверстии / А.С. Мокеев, В.М. Ямщиков // Компьютерная оптика. - 2021. - Т. 45, № 5. - С. 773-778. - DOI: 10.18287/2412-6179-CO-877.

Citation: Mokeev AS, Yamschikov VM. Comparison of numerical integration methods for calculating diffraction of a plane electromagnetic wave diffraction on rectangular aperture. Computer Optics 2021; 45(5): 773-778. DOI: 10.18287/2412-6179-CO-877.

Введение

Во многих задачах о распространении лазерного излучения в среде возникает необходимость учета волновой природы света, в частности, это задачи дифракции, фокусировки лазерного излучения, оценки расходимости и др. Их решение основано на вычислении интеграла Гюйгенса-Френеля [1 - 4]:

ер=± г Ем ехннм &, (1)

2к1 р

где Ер - амплитуда дифрагированной волны в точке наблюдения Р, Ем - амплитуда падающей волны в произвольной точке М волновой поверхности, к - модуль волнового вектора падающей волны, ёБ - площадь участка волновой поверхности, а р - абсолютная величина радиус-вектора от точки м до точки р.

Вычисление дифракционного интеграла - трудоемкая задача. Во-первых, подынтегральное выражение в нем представляет собой быстро осциллирующую функцию, что накладывает определенные требования как на выбор метода численного интегрирования, так и на шаг разбиения области интегрирования. Во-вторых, интеграл вычисляется по поверхности, площадь которой пропорциональна квадрату характерного размера пучка, а время, затрачиваемое на расчет, также пропорционально квадрату этого раз-

мера и обратно пропорционально шагу разбиения. Большой шаг разбиения неизбежно приводит к большим погрешностям вычислений, а малый - к большим временным затратам. Поэтому выбор оптимального шага и метода численного интегрирования является весьма актуальной задачей.

Вопросам разработки, исследования, тестирования и применения численных методов приближенного вычисления определенных интегралов посвящено большое количество работ [5 - 12].

Использование традиционных алгоритмов для вычисления интеграла от быстро осциллирующих функций подразумевает расчет по квадратурным формулам (метод прямоугольников, трапеций), что требует мелкого шага разбиения и приводит к значительным затратам вычислительных ресурсов [5 - 7].

Существуют и другие, более сложные и специализированные методы. Для интегралов от функций с фазой линейной формы часто используют метод Филона [6]. Он основан на построении составных квадратурных формул, в которых на каждом частичном интервале используется интерполяционный многочлен для амплитуды, и дальнейшее интегрирование выполняется точно.

Метод коллокаций Левина [9] предназначен для вычисления осциллирующих интегралов с более сложной фазовой функцией и может быть применим

для интеграла вида (1). Он заключается в переходе к решению обыкновенного дифференциального уравнения (ОДУ). Использование матрицы дифференцирования Чебышева позволяет далее свести задачу к решению системы линейных алгебраических уравнений (СЛАУ), чаще всего невырожденной. Основным недостатком данного метода является трудоемкость решения СЛАУ, поскольку получающаяся в результате матрица является плохо обусловленной, что может приводить к большой ошибке при получении решения. Тем не менее, во многих случаях использование метода кол-локаций Левина позволяет получить более точный результат по сравнению с другими методами [13].

В данной работе приведены результаты сравнения стандартных квадратурных методов интегрирования дифракционного интеграла (1) в приближении Фраун-гофера со специализированным методом коллокаций Левина и аналитическим решением на примере классической задачи дифракции плоской электромагнитной волны на прямоугольном отверстии в дальней зоне.

Аналитическое решение задачи дифракции плоской электромагнитной волны

на прямоугольном отверстии в дальней зоне

Для сравнения методов численного интегрирования с точным аналитическим решением рассмотрим классическую задачу дифракции плоской электромагнитной волны на прямоугольном отверстии в дальней зоне. Ее решение подробно изложено, например, в [1 - 2]. При нормальном падении плоской электромагнитной волны на прямоугольное отверстие шириной 2а0 и длиной 2Ь0, лежащее в плоскости ху, а также разложении величины р по малому параметру теории дифракции и пренебрежении квадратичными членами по сравнению с линейными р = 2 - (Хх+Уу) / 2 амплитуда дифрагированной волны Ер в точке наблюдения Р равна:

EP =

] ехрф dy J exp(^)dx, (2)

ikXxч

где (Х, У, 2) - координаты точки наблюдения Р, (х, у, 0) - координаты произвольной точки М волновой поверхности.

Интегралы (2) легко вычисляются аналитически, и тогда амплитуда электромагнитной волны:

EP =

iE04a0b0 7Л . ( 2a0X Л . ( 2b0Y

1Z

-exp(-ikZ) sinc

1Z

Sinc

1Z

. (3)

Здесь X - длина электромагнитной волны. Физической и измеряемой величиной является интенсивность света:

Ip = 10 Х| iO0b0Г Sinc2( 2a°X

2b„Y

, . _____ . „ . вшс21 —^ I, (4)

Ш ) У Ш ) у \2 )

где 10 = се0 Е02/2 - интенсивность падающей волны, с -скорость света в среде, е0 - электрическая постоянная.

Выражение (4) представляет собой аналитическое решение рассматриваемой задачи.

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

Рассмотрим построение решения интеграла (2) численным методом левых прямоугольников с постоянным шагом [5 - 6]:

| /(х^ « ^ /(х1) (X;+1 - х1) = /(х1),

а ;=1 ;=1

где х; = а + (( -1) И, N = (Ь - а) / И.

Применяя (5) к выражению (2), получим:

(5)

EP = Ahxhy

^ i kXhxmx Л Nr-' • (kXhxmx

Xcosl — J+Ъ™ {—

mx =1

x cosl |+i > Sin

i X sii

kYhymy

(6)

где

. ikE0 „

A =-exp(-ik

2nZ

v X (ao + hx) Y (bo + hy) Z +—--- + —-—

Z Z

\

Nx = 2a0/hx +1, Ny = 2b0¡hy +1. Как известно [14]:

sin (N a/2) • cos ([(N +1)/2] a)

Xcos (ja)=

j=1 N

Xsin (ja)=

j=1

sin (a/2) sin (N a/2) • sin ([(N + 1)/2]a)

(7)

sin (a/2)

При подстановке (7) в (6) получим выражение для поля электромагнитной волны:

EP = Ahx

81П

(Nx - 1)kXhx

2Z

NX7v1 cos| —kXhx 2Z

Б1П

kXhx

xhy

Б1П

(Ny - 1)kYhy

2Z

cos| NLkYh. 2Z y

81П

(Nx - 1)kXhx

2z

- +1-

NX1V1 sin| —kXhx 2Z

Б1П

kXhx ~TZ

81П

(Ny-1)kYhy 2Z

81П

kYhy

- +1-

sin| —^kYh, 2Z y

81П

kYhy ~2Z

(8)

Тогда аналитическое выражение для интенсивности, характеризующее решение рассматриваемой задачи численным методом левых прямоугольников:

I = I ( Кк>

12

2 БШ"

2ка0 X 12

2лЬ0У 12

лХИх 12

Б1П

2, лУИу

12

(9)

Как видно из решения (9), в окрестности точек:

пХ„Нх пУтИу

-= лп, -- = лт, п,т = 1,2,...,

12 12

(10)

возникает неопределенность вида 0 / 0. Устраняя эту неопределенность, в пределе е ^ 0 получим:

I=I • (£1 -

б1п2 {(Ых - 1)(лп + е)}

б1п2 (лп + е)

б1п2 {((у - 1)(лт + е)}

б1п2 (лт + е)

(11)

= 10

4арЬр 12

Формула (11) показывает, что значение интенсивности электромагнитной волны в точках (Хп, Ут) существенно отличается от аналитического решения (4) и совпадает со значением для центрального максимума. В пределе Их ^ 0, Иу ^ 0 формула (9) становится в точности равной аналитическому решению (4), которое не подразумевает наличие каких-либо максимумов вдали от центра экрана. То есть (Хп, Ут) -координаты центров артефактных максимумов. Их возникновение в распределении интенсивности при численном интегрировании дифракционного интеграла методом левых прямоугольников неизбежно и связано с шагом интегрирования, а период определяется выражением:

АХ = Хп+1 - Хп =12 / Их, АУ = Ут+1 - Ут =12 / Иу. (12)

Если в заданных пределах области наблюдения на экране Хтах и Утах выполняются соотношения АХ/Хтах < 1, АУ / Утах < 1, то метод левых прямоугольников неприменим для решения рассматриваемой задачи.

В связи с этим возникает необходимость введения ограничения на шаг интегрирования и количество узлов расчетной сетки:

2а0 ККтах

Иг = -

1

кХ,

тах ^ х

Иу =

кУ а

тах ^у

Nх = 1 + а х

Nу = 1 + а у

2

2Ь0 кУтах

(13)

где ах = f (Ах), Пу =/(Ау) - характерные параметры точности, зависящие от заданной абсолютной погрешности квадратурного метода Ах, Ау, в частности, для метода левых прямоугольников и трапеций

ах = а0/Ах, ау = Ь0/Ау и ах =а0/6Ах, ау = ^Ь0/6Ау соответственно.

Формулы (13) представляют собой критерий для оценки шага интегрирования или количества узлов расчетной сетки в зависимости от размеров области наблюдения на экране и требуемой точности вычислений.

В табл. 1 приведена зависимость относительной погрешности е,- вычисления центрального и двух последующих максимумов в распределении интенсивности на одномерном экране от характерного параметра точности ах. Расчет приведен для следующих параметров задачи: интенсивность падающей волны 10 = 250 Вт / см2, длина волны X = 1 мкм, размеры пучка по оси х и у а0 = Ь0 = 1 см, расстояние от отверстия до экрана 2 = 1 км, пределы наблюдения на экране

Хтах = 20 см, Утах = 0 см.

Табл. 1. Зависимость относительной погрешности е, от характерного параметра точности ах

ах, [ - ] е0, % (Х = 0 см) е1, % (Х = 7,15 см) е2, % (Х = 12,3 см)

0,3 4,66Е-14 11,21 37,90

0,5 1,28Е-13 4,81 15,05

0,8 1,16Е-14 1,69 5,13

1,6 8,15Е-14 0,42 1,25

8,0 6,75Е-13 0,017 0,05

Как видно из табл. 1, метод левых прямоугольников описывает центральный максимум в распределении интенсивности с точностью ~ 10 - 13 %, которая практически не зависит от параметра От. В то же время с увеличением ах погрешность вычисления интенсивности во втором и третьем максимумах уменьшается.

Таким образом, в заданных пределах области наблюдения на экране Хтах = 20 см в точке У = 0 см для описания центрального и двух последующих максимумов в распределении интенсивности I (Х, 0) с погрешностью не более 2 % характерный параметр точности должен составлять ах = 1,6, что соответствует ( = 41 узлу расчетной сетки.

Построение численного решения методом коллокаций Левина

Рассмотрим решение поставленной задачи численного интегрирования выражения для электрического поля (2) специализированным методом коллокаций Левина. Как известно [9], данный метод приближенного вычисления интегралов от быстро осциллирующих функций основан на переходе к численному решению ОДУ с последующим сведением к СЛАУ для определения первообразной от подынтегральной функции р (х).

Функция р (х), удовлетворяющая условию:

ё [р( х)ехр(,щ (х))] = / (х)ехр(,щ (х)), (14)

ах

позволяет вычислить осциллирующий интеграл:

b

J f (x) exp(iff)g (x)) dx =

a

= p(b)exp(irng (b)) - p(a)exp(io>g (a)).

(15)

Для определения функции p (x) достаточно удовлетворить уравнению:

dp(x) . dg (x)

'- + irn^-^p(x) = f (x).

dx

dx

(16)

Производную искомой функции йр / йх в дискретном наборе точек [х0,..., х^] можно представить в виде:

dp dx

= DP,

(17)

(-1)* - j

dj ,k - , 0 < j ф k < N, dKk = -—

tj — tk

1

21-tf

где р = [р (х0),...,р (х^]г Д - матрица дифференцирования Чебышева. Элементы матрицы дифференцирования вычисляются по следующим формулам [15]:

1

0 < к < N, d00 = б(1 + 2N2),

:<N, dk0 = 0 <k < N,

2 1 - U

dNtN = --1 (1 + 2N2), d0,k = 21-1-, 0 < k ■ 6 1 - tk 1 (- 1)N-k (-1)N-k 1 N 1 N

v ' 0 < k < N, dNk = -2^—, 0 < k < N, d0,N = -(-1), dNfi =— (-1)

1 + tk 2 2

(18)

dk,N = "

2 1 + tk

где

= ^(л; / N), 0 < ; < N -

система узлов Гаусса-Лобатто.

Так как интеграл рассматривается на отрезке [а, Ь], то для перехода к стандартной области задания полиномов Чебышева [-1, 1] необходимо провести соответствующие преобразования переменных:

b - a b + a

0 < j < N.

(19)

Таким образом, при подстановке (17), (18) и (19) в уравнение (16) получим СЛАУ следующего вида:

(D + /иЛ) р = f,

(20)

где

Р = [Р(xq),...,p(xN)f , f = [ f (x0),...,f (xn)f Л = digb -adg(xq) b - a dg(xN)

dx

dx

b - a г d

Дифференциальная матрица Д является сингулярной матрицей, но ее число обусловленности СоМ(Д) = ||Д|| ||Д-1|| улучшается при прибавлении к Д диагональной матрицы /юЛ. В этом случае ||/юЛ||2 показывает степень улучшения. Когда ||/юЛ||2 сравнительно велико, матрица (Д + /юЛ) становится хорошо обусловленной СоМ(Д + /юЛ) ~ 1, в противном случае она остается плохо обусловленной

СоМ(Д + /юЛ) >> 1, что ведет к большой ошибке при решении СЛАУ (20).

Применительно к рассматриваемой задаче с учетом (19) выражение (15) окончательно примет вид (21):

| / (х) ехрО'щд (х)) йх = | йх [ х) ехР('®Я (х))]йх = Ь--а [р( х„) ехр(/щд (х0)) - р( XN )ехр(/юя (XN))]. (21)

а -1

Тогда выражение (2) для амплитуды электромагнитной волны будет иметь вид:

EP =

ikE0a0b0 exp(-ikZ)

2л

Р1( xo)exp(ik

xq X

) - Р1( xn )exp(ik

y Y y Y

Р2 (yo) exp(ik^-j-) - p2 (yN) exp(ik

. (22)

Нетрудно заметить, что в окрестности точки X = Y = 0 матрица (D + /юЛ) ~ D, а значит, решение, полученное методом коллокаций Левина, может приводить к большой ошибке в распределении интенсивности вблизи центра экрана, поэтому возникает необходимость в выборе подходящего метода решения плохо обусловленной СЛАУ (20). В работе [13] отмечено, что надежным для решения данного типа задач является метод SVD (singular value decomposition), который и будет использоваться в дальнейшем.

Сравнение методов численного интегрирования с точным решением

Для сравнения двух квадратурных методов численного интегрирования (методы левых прямоуголь-

ников и трапеций), а также специализированного метода коллокаций Левина с точным аналитическим решением (4) были выбраны следующие параметры задачи: интенсивность падающей волны 10 = 250 Вт / см2, длина волны X = 1 мкм, размеры пучка по оси х и у а0 = Ь0 = 1 см, расстояние от отверстия до одномерного экрана 2 = 1 км, пределы области наблюдения на экране Хт£1х = 20 см, Утах = 0 см.

На рис. 1 представлены графики зависимости интенсивности электромагнитной волны I от координаты Х на экране при У = 0 см, полученные в результате численного расчета методами левых прямоугольников, трапеций и коллокаций Левина, а также точное

xj =

b

аналитическое решение. Численные расчеты квадратурными методами проводились на равномерной сетке с характерным параметром точности ах = 0,16, а специализированным методом - на сетке узлов Гаус-са-Лобатто. Количество узлов расчетной сетки во всех случаях составляло Ых = 5. Для решения СЛАУ (20) применялся метод 8УБ.

/, Вт/см2_

-о- Метод левых

прямоугольников о Метод трапеций А Метод коллокаций

-Аналитическое

решение

-20 -15 -10 -5 0 5 10 15 Х,см Рис. 1. Сравнение приближенного численного и точного аналитического решения: центральный максимум в увеличенном масштабе (а); второй максимум в увеличенном масштабе (б)

Из рис. 1а, б видно хорошее качественное совпадение рассмотренных методов приближенного численного интегрирования с точным аналитическим решением при описании центрального максимума. Метод коллокаций Левина дает более точный результат при описании второго и третьего максимума в распределении интенсивности. Метод левых прямоугольников при этом завышает результат, а метод трапеций - занижает. Также из данных, представленных на рис. 1, виден артефактный максимум в распределении интенсивности, координата которого в соответствии с выбранными параметрами расчета совпадает с периодом АХ = 20 см.

В табл. 2 приведены результаты сравнения методов численного интегрирования по относительной погрешности е, для центрального и двух последующих максимумов в распределении интенсивности. Количество узлов расчетной сетки во всех случаях составляло Ых = 41.

Табл. 2. Сравнение методов численного интегрирования по относительной погрешности е,

Метод численного интегрирования е0, % (Х = 0 см) е1, % (Х = 7,15 см) е2, % (Х = 12,3 см)

Левых прямоугольников 8,15Е - 14 0,42 1,25

Трапеций 8,15Е - 14 0,84 2,48

Коллокаций Левина 2,68Е - 13 4,32Е - 13 3,75Е - 13

Из табл. 2 видно, что на расчетной сетке, состоящей из Ых = 41 узла, все рассмотренные методы количественно описывают центральный максимум с точ-

ностью ~ 10 -13 %. Квадратурные методы оказываются в ~ 3 раза точнее метода коллокаций Левина. Второй и третий максимумы при численном интегрировании методом коллокаций Левина описываются с точностью ~ 10 -13 %. В то же время максимальная погрешность вычислений методами левых прямоугольников и трапеций не превышает 3 %.

На рис. 2 показан график зависимости интенсивности I во втором максимуме с координатой Х = 7,15 см при У = 0 см от количества узлов расчетной сетки Ых для трех методов численного интегрирования.

I,Вт/см2

3,0 2,5 2,0 1,5 1,0 0,5

,-'- —□ - Метод левых прямоугольников ......о- Метод трапеций ▲ Метод коллокаций Левина ---Точное значение

\

•

о

10° 101 102 103 Рис. 2. Исследование методов численного интегрирования на сходимость

Из рис. 2 видно, что специализированный метод коллокаций Левина обладает наилучшей сходимостью и позволяет получать решение поставленной задачи уже на сетке, состоящей из Ых ~ 5 узлов. В то же время для описания второго максимума в распределении интенсивности стандартными квадратурными методами численного интегрирования с погрешностью не более 1 % потребуется сетка из Ых ~ 40 узлов.

На рис. 3 представлена зависимость отношения времени t, затрачиваемого на расчет различными методами, ко времени затрачиваемому на расчет методом левых прямоугольников от количества узлов расчетной сетки Ых.

Шла

-— —□- Метод левых прямоугольников ......о- Метод трапеций —а- Метод коллокаций Левина -4

-— /

—Г

(

ж

/

Л.......<5!................ о.........о-а-а

К I -и-и

"А—'

0

100 10' 102 103

Рис. 3. Сравнение методов численного интегрирования по времени расчёта

Из данных, представленных на рис. 3, видно, что по временным затратам метод коллокаций Левина не уступает квадратурным методам вплоть до Nx ~ 40 узлов расчетной сетки. При Nx > 40 точек преимуществом по времени расчета обладают стандартные квадратурные методы численного интегрирования.

Заключение

На примере классической задачи дифракции плоской электромагнитной волны на прямоугольном отверстии в дальней зоне проведено сравнение трех алгоритмов численного интегрирования: метод левых прямоугольников, метод трапеций и специализированный метод коллокаций Левина.

Для квадратурных методов численного интегрирования получен критерий, определяющий выбор шага интегрирования в зависимости от размеров области наблюдения на экране и требуемой точности вычислений.

Показано, что для описания центрального и двух последующих максимумов в распределении интенсивности излучения на экране с погрешностью не более 2% метод коллокаций Левина требует ~ в 8 раз меньшее количество узлов расчетной сетки по сравнению с методами левых прямоугольников и трапеций и не уступает им по временным затратам вплоть до Nx ~ 40 узлов. При Nx > 40 точек преимуществом по времени расчета обладают стандартные квадратурные методы численного интегрирования.

Таким образом, на примере решения задачи дифракции плоской электромагнитной волны на прямоугольном отверстии в дальней зоне показаны преимущества использования специализированного метода коллокаций Левина над стандартными квадратурными методами численного интегрирования.

Литература

1. Ахманов, С.А. Физическая оптика [Электронный ресурс] / С.А. Ахманов. - М.: Московский государственный университет, 2004. - 656 с. - URL: https://ibooks.ru/bookshelf/27350/reading (дата обращения: 06.11.2020).

2. Борн, М. Основы оптики / М. Борн, Э. Вольф; под ред. Г.П. Мотулевич; пер. с англ. С.Н. Бреус, А.И. Головашкина, А. А. Шубина. - Изд. 2-е, испр. - М.: Наука, 1973. - 720 с.

3. Устинов, А. В. Быстрый способ вычисления интеграла Релея-Зоммерфельда первого типа / А.В. Устинов // Компьютерная оптика. - 2009. - Т. 33, № 4. - С. 412-419.

4. Veerman, J.A.C. Calculation of the Rayleigh-Sommerfeld diffraction integral by exact integration of the fast oscillating factor / J.A.C. Veerman, J.J. Rusch, H.P. Urbach // JJournal of the Optical Society of America. - 2005. -Vol. 22, Issue 4. - P. 636-646.

5. Калиткин, Н.Н. Численные методы / Н.Н. Калиткин. -М.: Наука, 1978. - 512 с.

6. Самарский, А.А. Численные методы / А.А. Самарский, А.В. Гулин. - М.: Наука, 1989. - 432 с.

7. Соболь, И.М. Метод Монте-Карло / И.М. Соболь. - М.: Наука, 1978.

8. Jeffrey, G.B. Louis Napoleon George Filon, 1875-1937 / G.B. Jeffrey // Obituary Notices of Fellows of the Royal Society. - 1939. - Vol. 2, Issue 7. - P. 500-509. - DOI: 10.1098/rsbm.1939.0010.

9. Levin, D. Procedures for computing one and two-dimensional integrals of functions with rapid irregular oscillations / D. Levin // Mathematics of Computation. - 1982. -Vol. 38, No. 158. - P. 531-538.

10. Li, J. A universal solution to one-dimensional oscillatory integrals / J. Li, X. Wang, T. Wang // Science in China Series F: Information Sciences. - 2008. - Vol. 51, Issue 10. -P. 1614-1622. - DOI: 10.1007/s11432-008-0121-2.

11. Lovetskiy, K.P. Integration of highly oscillatory functions / K.P. Lovetskiy, L.A. Sevastyanov, A.L. Sevastyanov, N.M. Mekeko // Mathematical Modelling and Geometry. -2014. - Vol. 2, Issue 3. - P. 11-27.

12. Liu, Y. Fast evaluation of canonical oscillatory integrals Y. Liu // Applied Mathematics & Information Sciences. -2012. - Vol. 6, No 2. - P. 245-251.

13. Ловецкий, К.П. Сравнение методов вычисления интегралов от быстро осциллирующих функций [Электронный песурс] / К.П. Ловецкий, И.А. Мигаль // Науковедение. - 2015. - Т. 7, № 2. - 16 с. - URL: http://naukovedenie.ru/PDF/70TVN315.pdf (дата обращения 15.04.2021). - DOI: 10.15862/70TVN315.

14. Фихтенгольц Г.М. Курс дифференциального и интегрального исчисления. Том 2 : Учеб. пособие / Г.М. Фихтенгольц. - 7-е изд. - М.: Наука, 1969. - 800 с.

15. Mason, J.C. Chebyshev polynomials / J.C. Mason, D.C. Handscomb. - Chapman and Hall/CRC, 2002. - 360 p.

Сведения об авторах

Мокеев Александр Сергеевич, 1993 года рождения, в 2016 году окончил Национальный исследовательский ядерный университет «МИФИ» по специальности 140305 «Ядерные реакторы и энергетические установки», работает научным сотрудником в институте лазерно-физических исследований ИЛФИ РФЯЦ - ВНИИЭФ. Область научных интересов: компьютерная оптика, численные методы, программирование. E-mail: mokeev-a@mail.ru .

Ямщиков Виталий Михайлович, 1994 года рождения, учился в Национальном исследовательском ядерном университете «МИФИ». Работает младшим научным сотрудником в институте лазерно-физических исследований ИЛФИ РФЯЦ - ВНИИЭФ. Область научных интересов: перенос излучения, физика плазмы, лазерная физика. E-mail: _ feynman94@yandex.ru.

ГРНТИ: 27.41.41(алгоритмырешения задач вычислительной и дискретной математики) Поступила в редакцию 9 февраля 2021 г. Окончательный вариант - 21 апреля 2021 г.

Comparison of numerical integration methods for calculating diffraction of a plane electromagnetic wave by a rectangular aperture

A.S. Mokeev1, V.M. Yamshchikov1 1 Federal State Unitary Enterprise RUSSIAN FEDERAL NUCLEAR CENTER All-Russian Research Institute of Experimental Physics (FSUE "RFNC - VNIIEF"),

The Institute of Laser Physics (ILFI), 607188, Sarov, Nizhny Novgorod region, Mira Ave, 37

Abstract

We discuss features of the calculation of a Fraunhofer integral by traditional quadrature numerical integration methods and a special collocation Levin method when calculating the diffraction of a plane electromagnetic wave by a rectangular aperture. For the quadrature numerical integration methods, a criterion for the assessment of the integration step is derived depending on the screen size and required calculation accuracy. Advantages of the use of the special collocation Levin method in comparison with the traditional quadrature numerical integration methods are shown.

Keywords: diffraction integral, integration of oscillatory functions, method of rectangles, trapezium method, Levin method, Fraunhofer diffraction.

Citation: Mokeev AS, Yamschikov VM. Comparison of numerical integration methods for calculating diffraction of a plane electromagnetic wave diffraction on rectangular aperture. Computer Optics 2021; 45(5): 773-778. DOI: 10.18287/2412-6179-CO-877.

References

[1] Akhmanov SA, Nikitin SYu. Physical optics. Oxford: Clarendon Press; 1997.

[2] Born M, Wolf E. Principles of optics. Electromagnetic theory of propagation, interference and diffraction of light. 7th ed. Cambridge: Cambridge University Press, 1999.

[3] Ustinov AV. The fast way for calculation of first class Rayleigh-Sommerfeld integral [In Russian]. Computer Optics 2009; 33(4): 412-419.

[4] Veerman JAC. Calculation of the Rayleigh-Sommerfeld diffraction integral by exact integration of the fast oscillating factor. J Opt Soc Am 2005; 22(4): 636-646.

[5] Kalitkin NN. Numerical methods [In Russian]. Moscow: "Nauka" Publisher; 1986.

[6] Samarskii AA, Gulin AV. Numerical methods [In Russian]. Moscow: "Nauka" Publisher; 1989.

[7] Sobol IM. A primer for the Monte Carlo method. Boca Raton, FL: CRC Press; 1994.

[8] Jeffery GB. Louis Napoleon George Filon, 1875-1937. Obit Not Fell R Soc 1939; 2(7): 501-509. DOI: 10.1098/rsbm.1939.0010.

[9] Levin D. Procedures for computing one and two-dimensional integrals of functions with rapid irregular oscillations. Math Comp 1982, 38(158): 531-538.

[10] Li J, Wang X, Wang T. A universal solution to one-dimensional oscillatory integrals. Sci China Ser F-Inf Sci 2008; 51(10): 1614-1622. DOI: 10.1007/s11432-008-0121-2.

[11] Lovetskiy KP, Sevastyanov LA, Sevastyanov AL, Mekeko NM. Integration of highly oscillatory functions. Mathematical Modelling and Geometry 2014; 2(3): 11-27. DOI: 10.26456/mmg/2015-312.

[12] Liu Y. Fast evaluation of canonical oscillatory integrals. Appl Math Inf Sci 2012; 6(2): 245-251.

[13] Lovetskiy KP, Migal IA. Comparison of methods for calculation of oscillatory integrals {In Russian]. Nau-kovedenie 2015; 7(2). Source: <http://naukovedenie.ru/PDF/70TVN315.pdf).

[14] Fikhtengol'ts GM. Differential- und Integralrechnung. II [In German]. Berlin: VEB Deutscher Verlag der Wissenschaften; 1978.

[15] Mason JC, Handscomb DC. Chebyshev polynomials. Chapman & Hall /CRC Press; 2002.

Authors' information

Aleksandr Sergeevich Mokeev (b. 1993) graduated from National Research Nuclear University "MEPHI" in 2016, majoring in Nuclear Reactors and Power Plants. Currently he works as the junior researcher at the Institute of Laser Physics (ILFI) FSUE "RFNC - VNIIEF". Research interests are computer optics, numerical methods, and programming. E-mail: mokeev-a@mail.ru .

Vitaliy Mikhailovich Yamshchikov (b. 1994), studied at the National Research Nuclear University "MEPHI". Works as research engineer. Research interests: transfer of radiation, physics of plasma, laser physics. E-mail: _ feynman94@yandex.ru .

Received February 9, 2021. The final version - April 21, 2021.