Сравнение методов вычисления интегралов от быстро осциллирующих функций

Ловецкий Константин Петрович; Мигаль Илья Александрович

Интернет-журнал «Науковедение» ISSN 2223-5167 http ://naukovedenie. ru/ Том 7, №3 (2015) http://naukovedenie.ru/index.php?p=vol7-3 URL статьи: http://naukovedenie.ru/PDF/70TVN315.pdf DOI: 10.15862/70TVN315 (http://dx.doi.org/10.15862/70TVN315)

УДК 519.644; 517.584

Ловецкий Константин Петрович

ФГБОУ ВПО «Российский университет дружбы народов»

Россия, Москва1 Доцент

Кандидат физико-математических наук E-mail: [email protected] РИНЦ: http://elibrary.ru/author profile.asp?id=584101

Мигаль Илья Александрович

ФГБОУ ВПО «Российский университет дружбы народов»

Россия, Москва Студент Бакалавр

E-mail: [email protected]

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

осциллирующих функций

1 115419, Москва, ул. Орджоникидзе, 3

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

Для отыскания первообразной осуществляется переход к решению системы обыкновенных дифференциальных уравнений (ОДУ) без граничных условий на интервале интегрирования. В работе исследованы известные методы, основанные на подходах, предложенных Филоном и Левиным. Предложены варианты построения матриц дифференцирования, приводящие к возможности устойчиво решать получающиеся системы линейных алгебраических уравнений с последующим вычислением интегралов от быстро осциллирующих функций. Преимущества предложенных методов продемонстрированы на ряде численных примеров.

Ключевые слова: осциллирующие функции; интегралы от быстро осциллирующих функций; метод коллокации Левина; ии разложение; сингулярное разложение; полиномы Чебышева; матрица дифференцирования; узлы Гаусса-Лобатто; устойчивые методы решения систем линейных алгебраических уравнений.

Ссылка для цитирования этой статьи:

Ловецкий К.П., Мигаль И.А. Сравнение методов вычисления интегралов от быстро осциллирующих функций // Интернет-журнал «НАУКОВЕДЕНИЕ» Том 7, №2 (2015) http://naukovedenie.ru/PDF/70TVN315.pdf (доступ свободный). Загл. с экрана. Яз. рус., англ. DOI: 10.15862/70TVN315

Интернет-журнал «НАУКОВЕДЕНИЕ» Том 7, №3 (май - июнь 2015)

http://naukovedenie.ru [email protected]

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

V

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

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

d [ p ( л ) eilBg (л) ] = f ( л ) eilBg (л) (0.1)

Если продифференцировать р(х)ешд(х\ то получим:

d[Р(л)ewg(л)] = [p (л) + ig (л)p(л)]e'mg(л) = f (x)ewg(л) (0.2)

В итоге для определения р (х) достаточно удовлетворить уравнению:

Р'(л) + ig'(л) p (л) = f(л) (0.3)

После нахождения функции р (х) можно вычислить и сам интеграл

I [ f] = ]f (л ) e^d = ]d [ p (л ) e^] dx = p (b ) ei(0g(b) - p(a)eiWg(a) (0.4)

dx

a a

Таким образом, главная цель применения метода коллокации Левина - нахождение функции р(х). Именно благодаря формуле (0.3) задачу можно свести к решению ОДУ.

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

Суть метода - отыскания первообразной р(х)удовлетворяющей формуле (0.3). Нахождение неизвестной функции заключается в разложении ее

ад

p (л )=&% (л) (а5)

k=0

по базису {фк(х))1° в некотором гильбертовом пространстве. Использование достаточно большого количества членов ряда позволяет достичь приемлемой точности. Пусть существует оператор Ь(р): Ь[р](х) = /(х). Требуется, чтобы при некоторых коэффициентах ак,к = 1, ...,п выполнялись следующие равенства

L

Srn (x)

k=0

= f (x) J = 0,---n.

Метод коллокаций относится к спектральным методам. Он заключается в нахождении коэффициентов разложения решения дифференциального уравнения (0.3) по значениям L[p](x) = f(x) на заданной системе узлов {х0, ...,хп}. Эти коэффициенты а^ могут быть представлены как решение системы уравнений:

L[p]{x0) = f(x0) L[p]{x„) = f(xn)

и задача о вычислении интеграла от быстро осциллирующих функций может быть сведена к решению системы линейных алгебраических уравнений (СЛАУ).

Матрица дифференцирования. Если в качестве базисных функций выбираются полиномы Чебышева первого рода, а точки сетки являются узлами Гаусса-Лобатто, тогда производную искомой р(х) можно представить в виде:

p'(x) = Dp( x)

Матрица D называется матрицей дифференцирования Чебышева. Элементы матрицы дифференцирования вычисляются по следующим формулам [1]:

d}.,k = (-l)k-j ,0 < j ф k < n, dkk = -,0 < k < n

' 21 - Ук

d-0,0 = - (l + 2n2), dn,n =--(l + 2n2),

(-l)k 1 (-l)k d0k = 2^ —,0 < k < n, dk0 =--^-,0 < k < n, (0.6)

l - yk 2l - yk

l (-l)n-k (-l)n-k dk,n = <k <n, dnk = -2^^,0<k <n,

2 l + yk l + yk

d0,n = l (-l)n, dn, =--(-1)п.

Метод квадратур. Если рассматривать интеграл на отрезке x е [a, b], то для перехода к стандартной области задания полиномов Чебышева [-l, l] необходимо провести замену переменных. Это возможно благодаря следующему преобразованию:

b - a b + a

x =-1 +-

2 2

2 2 (0.7)

P' (x ) = т-P '(t)

b - a

Подставляя узлы Гаусса-Лобатто получаем:

Ь - а х. =-соб

/ = соб

г

N

ж у I Ь + а . .т

— 1 +-, У = 0,1,..., N

V N ) 2

(0.8)

Векторы значений функций и их производных в точках Гаусса-Лобатто можно записать следующим образом [2], [1]:

(0.9)

Р = [Р ( хо) Р ( х1 )--Р ( ^ )]Т Р' = [р'( хо ) Р' ( Х1 )--Р'( ^ )]Т

g' = [g'( Хо ) g'( Х )... g'( XN )]' / = [/ ( Хо ) / ( Х1).../( XN )]Т

Следовательно, исходя из формулы (0.9), формулу (0.3) можно представить в виде:

Р (Ху) + V (х у) Р (Ху ) = / (Ху) ,у = 0,...я. (°Л°)

2

Вектор Р можно представить в виде р' =-Бр , и формула (0.10) принимает вид:

Ь - а

(011)

Р ' (Х0) g' (Х0 ) Р(Х0) " / (Х0)

Р'( Х1) +1 g'(Х1) р( х1) = / (Х1)

_ Р' ( XN )] _ g' (^) р( ^) ] / (XN )

Подставив формулы (0.7), (0.9) в (0.11) получаем:

( Б + Л ) Р = 1/,

Я =

Ь - а 2

(012)

Матрица Л = й1ад{Хд'{х0),Хд'{х1), ..,Ад'(хы), является диагональной. Решая систему линейных уравнений (0.12) относительно вектора р, значения искомого интеграла вычисляются в соответствии с формулой (0.4). Для решения СЛАУ были использованы методы ЬИ и БУО. Причиной выбора ЬИ метода стала его скорость, однако он не всегда стабилен. SVD не такой быстрый, тем не менее, более надежен, чем ЬИ метод.

Альтернативный метод квадратур. Для вычисления производной искомого решения будем считать, что в качестве базисных функций выбраны полиномы Чебышева первого рода

Т(х) . Выражение для производной р'(х) имеет вид:

( ( [ " I "

(р(Х))=~т Е акТк(х) I=Е ЬкТк(Х)=Р ' (Х).

(Х V к=0 ) к=0

(х

(0.13)

Известно, что полиномы Чебышева первого рода могут быть вычислены с использованием рекуррентного соотношения Тк ( х) = 2ХТк_ х ( х) - Тк_2 ( х), к = 2,3,..., с начальными условиями т0 (х) = 1, Т1 (х) = х . Данное соотношение позволяет вывести

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

рекуррентное соотношение между коэффициентами разложения самой функции и коэффициентами ее производной [3]:

10 -1/2

1/4 0 -1/4

1/6 0

1/8

-1/6 0

-1/8

1

2п - 2

0

1

2п

Используя (0.14), можно вывести формулу, связывающую эти коэффициенты в «обратном» порядке с помощью матрицы дифференцирования А в пространстве коэффициентов:

" Ь0 " а

Ь а2

Ъ2 а

¿3

К-2 ап-1

_ Ьп-Х _ _ ап _

(0.14)

где А =

Аа = Ь

(1/ с)х 2у если / > у и сумма / + у - нечетна

(0.15)

0

иначе,

Здесь 0 < /, у < п, с ='

2 / = 0,

1 / > 0.

Тогда формулу (0.3) можно представить в виде:

ТАа + iЛdiag (g ')Та = Л /

(0.16)

где Т - матрица размерности (п +1) х (п +1), элементами которой являются значения полиномов Чебышева первого рода в узлах сетки. Подробнее формулу (0.16) можно записать в виде:

Т Т Т

10,0

1,0

2,0

т„

и,0

Т Т Т

Т„

Т Т Т

-'-по -1! о о

0,2

1,2

2,2

Т.

п,2

Т Т Т

Т 0,п Т1,п Т2,п

т

0 1 о

0 4 0

3 0 6

0 а,

ап

а

а0

"g0 " Гт Т0,0 Т Т1,0 т 2,0 т 1 ао

g; Т Т Т1,1 т 2,1 т и,1 а\ /1

/Л g 0 Т Т 1,2 т 2,2 т а2 = Л ./2

_ ^п _ т 10,п Г Т1,п т 2,и ■ т п,п ап_ _/п _

(0.17)

Здесь использованы обозначения Т , = Т (Х) и g' = g'(х ) .

Решение этой системы линейных алгебраических уравнений относительно коэффициентов а = (а0,а,•••,а„) разложения решения по базисным функциям позволяет определить приближенное значение интеграла по формуле (0.4).

Преобразование интеграла в нерегулярном случае к каноническому интегралу.

Пусть дан интеграл вида:

ъ

|/ (х) ег(Ве (х)ёх (0.18)

а

В случае, когда производная фазовой функции не обращается в ноль в точках используемой сетки Гаусса-Лобатто на интервале интегрирования, интеграл можно преобразовать к (более устойчивому с точки зрения решения соответствующей СЛАУ) виду:

е(ъ) (е4 ( V))

1 [? 1= 1 (е -1 (У)) ^¿У (0.19)

е (е (у))

( а )

где введено обозначение у = д (х). Интеграл можно упростить, применяя

гг д(Ь)-д(а) д(Ъ)+д(а)

преобразование координат у =---х +----:

1

I [ /1 = 1^ ( х ) е*хах, (0.20)

-1

где F(x) = Л^ и * = ^^Ъ^.

у у 2 д'(д~1(у)) ъ 2

Для решения поставленной задачи необходимо найти обратную функцию от д(х) -функцию д-1(х). Численное значение д-1(х) можно найти [4] применив метод Ньютона. Часто можно легко найти и аналитический вид функции. Например, для интеграла

1

о

обратная функция на отрезке [0,1] имеет вид

у2 -106 3 =х = Т^ЮЗ'

Однако отыскание аналитического вида функции порой бывает проблематичным. Для нахождения значений обратной функции был адаптирован метод дихотомии. Предполагается, что функция д(х) является гладкой, непрерывно дифференцируемой на отрезке [а,Ь] и на этом отрезке нет стационарных точек. Ниже представлен псевдокод алгоритма:

§'(аггау[] х, аггау[] у, У_ТО_БШБ) {

//массив значений х //массив значений у

Х_ТО_ЕШБ; х1 = 0; х2 = 0; тёех1 = 0;

index2 = 0;

//если значение Y_TO_FIND было подсчитано - возвращаем соответствующее

значение х

if (y.Contains(Y_TO_FIND))

{

for (int i = 0; i < y. Length; i++) {

if (Y_TO_FIND == y[i]) return x[i];

}

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

}

//если y[0] - самое большое число в массиве

//ищем такие значения y[i] и y[i-1], чтобы удовлетворяло условию // y [i]>Y_TO_FIND>y [i-1] //начинаем с начала массива

else if (y[0] > y[y.Length - 1])

{

for (int i = 1; i < y. Length; i++) {

if (Y_TO_FIND > y[i] && Y_TO_FIND < y[i - 1])

{

index2 = i-1; index1 = i; break;

}

//если y[last] - самое большое число в массиве

//ищем такие значения y[i] и y[i-1], чтобы удовлетворяло условию // y [i]>Y_TO_FIND>y [i-1]

//начинаем с конца массива

else if (y[0] < y[y.Length - 1])

{

for (int i = y.Length-1; i > 0; i--)

{

if (Y_TO_FIND < y[i] && Y_TO_FIND > y[i - 1])

{

index2 = i; index1 = i - 1; break;

}

//когда были найдены значения y[i] и y[i-1] удовлетворяющие условию // можно выбрать такие x[i] > X_TO_FIND > x[i-1] // индексы реализованы так. что g(x[i]) = y[i] x1 = x[index1];

x2 = x[index2];

// находим точку посередине

X_TO_FIND = (x1 + x2) * 0.5;

//подсчитваем значени y = g(X_TO_FIND) в данной точке

double y_test = gFunction(X_TO_FIND);

//error

double err = 10e-16;

//если ( |Y_TO_FIND-Y_TEST| )A2 > err

while (Math.Pow(Math.Abs(Y_TO_FIND - y_test),2) > err) {

//if calculates y_test if (Y_TO_FIND >= y_test)

x1 = X_TO_FIND;

else if (Y_TO_FIND <= y_test)

x2 = X_TO_FIND;

else { throw new Exception("Problem finding inverse of g(x)"); } X_TO_FIND = (x2 + x1) * 0.5; y_test = gFunction(X_TO_FIND);

}

return X_TO_FIND;

}

Программная реализация. В рамках исследования разработана программа на языке программирования С#, использующая математические библиотеки ALGLIB и Math.NET Numerics. Ниже представлена краткая логика работы программы:

1. Задать границы отрезка на оси х (значения а и b).

2. Указать количество узлов на сетке (должно быть >= 2).

3. Указать метод решения СЛАУ (LU, SVD).

4. Выбрать вид интеграла (формула (0.18) или (0.20)).

Вычисление интеграла на интервале без стационарных точек. Пример для оценки работы программы был взят из статьи [2].

i

I = j sin(x)e500(x2+x)dx .

о

Метод решения СЛАУ - SVD. Применяется SVD метод так, как он более устойчив при решении систем с плохо обусловленными матрицами.

Точное численное значение интеграла:

/ = (4.59859397840143 - i X 3.15443542737400) X 10-4

Ниже представлена оценка точности вычисления интеграла на основе данных полученных при помощи программы.

Интернет-журнал «НАУКОВЕДЕНИЕ» Том 7, №3 (май - июнь 2015)

http://naukovedenie.ru [email protected]

Таблица 1

Зависимость погрешности от количества узлов при вычислении интеграла на сетке с граничными значениями по формуле (0.20)

Узлы Яеа1-Егг-Метод Левина с преобразованием Img-Err-Метод Левина с преобразованием

15 4,17731299176523Е-06 2,55537162377136Е-06

30 2,28380628197562Е-10 3,57509279549190Е-10

45 7,65363661147884Е-11 2,98424293106377Е-10

60 5,80199577692515Е-11 1,68334038918193Е-11

75 5,40202093913642Е-10 1,89875544053685Е-10

90 1,66935783836325Е-10 1,33843193119391Е-10

Таблица 2

Зависимость погрешности от количества узлов при вычислении интеграла на сетке с граничными значениями по формуле (0.18)

Узлы ЯеаЬЕгг-Метод Левина 1т^-Егг-Метод Левина

15 1,92143293374891Е-09 2,79268047206601Е-09

30 2,28341490678251Е-13 5,32574308444421Е-13

45 2,23979779188785Е-12 7,66862322223794Е-12

60 5,50813430662947Е-12 1,63578849787939Е-11

75 1,73900258244059Е-11 1,83804300049579Е-11

90 7,48057097262249Е-12 8,06474577175856Е-12

1,00Е-01

1,00Е-03

*

Ч 1,00Е-05

3 о

о 1,00Е-07

4 к о.

£ 1,00Е-09

1,00Е-11

1,00Е-13

15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100

Количество Узлов

1^еа!-Егг-Метод Левина с преобразованием ■ 1^еа!-Егг-Метод Левина !т§-Егг-Метод Левина с преобразованием ^^ !т§-Егг-Метод Левина

Рис. 1. Зависимость погрешности от количества узлов при вычислении интеграла на

интервале с граничными значениями.

Оценка погрешности вычисляется по формулам:

Er -

jnum т real real

I

real

Er -

img

num

img img

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

I

(0.21)

/кит *■» тпи^т

геа1 - действительная часть полученного численного значения интеграла, ¡тд -мнимая часть полученного численного значения интеграла, 1геа1 - действительная часть точного значения интеграла, 1тд - мнимая часть точного значения интеграла.

В данном примере отсутствуют стационарные точки. Из Таблица 1 и Таблица 2, а так же из Рис. 1 можно заметить, что результаты, полученные при вычислении интеграла по формуле (0.18) обладают большей точностью, нежели результаты при вычислении интеграла по формуле (0.20).

Алгоритм вычисления интеграла по формуле (0.18) достигает наименьшей погрешности при количестве узлов равном 35.

Вычисление интеграла на интервале со стационарными точками.

Пример был взят из статьи [4]. Ниже представлен сам интеграл:

К - + -) ^

х+1)

(0.22)

Метод решения СЛАУ - БУО. Применяется БУО метод поскольку он более устойчив и прекрасно работает с плохо обусловленными матрицами.

Точное численное значение интеграла:

! = -0.393011626656505 + 0.601601971947752 *1

Ниже представлена оценка точности вычисления интеграла на основе данных, полученных при помощи программы на интервале [-1,1].

Таблица 3

Зависимость погрешностей от количества узлов при вычислении интеграла на интервале с граничными значениями по формуле (0.18)

Узлы Real-Err-Метод Левина Img-Err-Метод Левина

20 3,94146E-12 9,5129E-12

30 6,28741E-12 7,54602E-13

40 9,8208E-13 1,526E-12

50 1,58011E-12 1,16448E-13

60 5,72468E-13 8,3045E-15

70 1,20087E-12 1,96171E-13

80 6,69222E-13 3,84037E-13

90 2,22645E-12 3,80826E-12

100 4,83484E-13 1,87866E-13

5

*

Ю 5

3

о *

о

4

к о. о

1,00E-01 1,00E-03 1,00E-05 1,00E-07 1,00E-1,00E-11 1,00E-13 1,00E-

20 30 40 50 60 70 Количество Узлов

80

90

100

Real-Err

Img-Err

Рис. 2. Зависимость погрешностей от количества узлов

Оценка погрешности была подсчитана при помощи формул (0.21).

В данном примере имеется стационарная точка х = -1. Поэтому вычисление обратной функции в ней невозможно, поскольку значение производной от фазовой функции g'(х) в данной точке равно 0.

Эту проблему можно решить, если для определения первообразной р(х) использовать сетку без граничных значений, узлы которой можно вычислить по формуле

2п

Xj = cos , j = 1,..., n (0.23)

Используя альтернативный метод квадратур (0.16) для подсчета коэффициентов а = ( а0, а,..., а ) на сетке, узлы которой вычисляются по формуле (0.23), и формулу (0.4) для вычисления интеграла, были получены следующие результаты:

Таблица 4

Зависимость погрешности от количества узлов при вычислении интеграла на сетке

без граничных значений по формуле (0.20)

Узлы Яеа1-Егг-Метод Левина с Img-Err-Метод Левина с

преобразованием преобразованием

10 0,000135723 0,000148638

20 0,000326519 0,000782133

30 0,000119972 0,000406786

40 3,2008Е-05 3,78245Е-05

50 3,78409Е-06 2,18732Е-05

60 0,000113298 5,32805Е-05

70 8,23595Е-05 7,38884Е-05

80 3,16898Е-05 3,95543Е-05

По результатам из Таблица 4 видно, что погрешность данного метода вычисления интеграла достигает минимального стабильного значения 105 при количестве узлов равном 40. Тем не менее, точность результатов на сетке без граничных значений хуже, чем на сетке с граничными значениями.

В то же время прямой алгоритм без применения преобразований дает более точные результаты на сетке (со стационарной точкой х = -1), включающей граничные точки. Уже на сетке из 20 узлов погрешность составляет 10-12. Это можно увидеть из Таблица 3 и Рис. 2. Зависимость погрешностей от количества узлов

Алгоритм SVD решения системы (0.12) оказался в этом случае весьма устойчивым. Это показывает, что используя метод без преобразования интеграла на сетке, включающей стационарные точки фазовой функции, можно получить результаты с малой погрешностью.

Заключение. В работе представлен сравнительный анализ методов вычисления интегралов от быстро осциллирующих функций и разработана программа на языке C# с применением библиотек ALGLIB и Math.NET Numerics.

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

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

ЛИТЕРАТУРА

1. Mason J.C., Handscomb D.C. Chebyshev Polynomials. — New York: A CRC Press Company, 2003.

2. Wang XueSong, Li JianBing, Wang Tao A universal solution to one-dimensional oscillatory integrals // Science in China Series F: Information Sciences. — 2008. — V. 51, 10. — P. 1614-1622.

3. Fornberg Bengt A Practical Guide to Pseudospectral Methods. — New York: Cambridge University Press, 1996.

4. Liu Ying Fast Evaluation of Canonical Oscillatory Integrals // Applied Mathematics & Information Sciences. — Natural Sciences Publishing Cor., 2012. — V. 6. — P. 245251.

5. Canuto C., Hussaini M.Y., Quarteroni A Spectral Methods in Fluid Dynamics // Springer. — New York, 1988.

6. Olver Sheehan Fast, numerically stable computation of oscillatory integrals with stationary points. — Oxford, UK: Oxford University Computing Laboratory, 2009.

7. Olver Sheehan Numerical Approximation of Highly Oscillatory Integrals. — Cambridge: University of Cambridge, June 14, 2008.

8. Xiang Shuhuang Efficient Filon-type methods for f^ f(x)elwg(x)dx // Numerische Mathematik. — Hunan, China, 2007. — Vol. 105, 4. — P. 633-658.

9. Shen Jie, Tang Tao, Wang Li-Lian Spectral Methods. Algorithms, Analysis and Applications. — Berlin Heidelberg : Springer-Verlag, 2011.

10. Boyd John P. Chebyshev and Fourier Spectral Methods. — Mineola, New York 11501: DOVER Publications, Inc., 2000.

11. Goncharsky A.V. Goncharsky A.A. Computer optics & computer holography. — M.: : Moscow University Press, 2004.

12. Гончарский А.А., Дурлевич С.Р. Нанооптические элементы для формирования 2D изображений // вычислительные методы и программирование. — 2010. — Т. 11. — C. 246-249.

13. Гончарский А.А., Дурлевич С.Р. Об одной задаче синтезананооптических элементов для формирования динамических изображений // вычислительные методы и программирование. — 2013. — Т. 14. — C. 343-347.

14. Lovetskiy K.P., Sevastyanov L.A., Sevastyanov A.L., Mekeko N.M. Integration of highly oscillatory functions // Mathematical Modelling and Geometry. — 2014. — V 3. — No. 3. — P. 11-24.

Рецензент: Севастьянов Леонид Антонович, профессор, доктор физико-математических наук.

Lovetskiy Konstantin Petrovich

The Federal State Budgetary Educational Higher Professional Institution of Education

"People's Friendship University of Russia" (PFUR)

Moscow, Russian Federation E-mail: [email protected]

Migal Ilya Aleksandrovich

The Federal State Budgetory Educatinal Higher Professional Insitution of Education

"People's Friendship University of Russia" (PFUR)

Moscow, Russian Federation E-mail: [email protected]

Comparison of methods for calculation of oscillatory integrals

Annotation. In this article numerical methods of calculating oscillatory integrals are presented. The tasks, which are reduced to integration of such functions, arise in many appendices. In recent years new methods appeared, allowing changing the task from integration of oscillating functions to restoration of primitive function and calculation of required integral with their help. It appeared that the problem could be solved very effectively, and, the stronger are the oscillation, the more accurate is the result.

To obtain the antiderivative, a transition to the system of the ordinary differential equations (ODE) without boundary conditions on an integration interval is carried out. In this work the known methods based on the approaches offered by Fillon and Levin are investigated. The work overviews proposed algorithms for creating differentiation matrices, leading to opportunity of steadily solving the linear algebraic equations with the subsequent calculation of oscillatory integrals. Advantages of the offered methods are shown on a number of numerical examples.

Keywords: oscillatory functions; oscillatory integrals; collocation method; LU decomposition; Singular value decomposition; Chebyshev Polynomials; differential matrix; Gauss-Lobatto nodes; stable methods of solving algebraic linear system of equations.

REFERENCES