Интегрирование быстро осциллирующих функций Текст научной статьи по специальности «Математика»

УДК 519.644; 517.584

Интегрирование быстро осциллирующих функций

К. П. Ловецкий, В. В. Петров

Кафедра систем телекоммуникаций Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия

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

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

Ключевые слова: интегрирование осциллирующих функций, метод Филона, метод Левина, чебышевская матрица дифференцирования, вырожденные матрицы, функции Бесселя.

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

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

Рассмотрим интеграл вида:

при этом область интегрирования такова, что константа осцилляций ш ^ 1 — «большая» величина; амплитуда f (х) и фаза д(х) — достаточно гладкие не осциллирующие функции. Даже в случае линейной функции д(х) = х функции И,е (/(х)ешхЛ) и 1т (/(х)ешх) имеют на [а,Ь] примерно ш(Ъ — а)/л нулей, поскольку фазовый множитель изменяется в области интегрирования на величину порядка ш(Ъ — а). Число периодов Р(х) ~ ш(Ъ — а)/(2л), и на каждом периоде имеется 2 корня. Таким образом в области интегрирования имеется порядка ш(Ъ — а)/л нулей.

Это означает, что величина, например равномерного, шага Ь = (Ь — а)/п ^ 1/ш при численном интегрировании должна быть мала или должно быть велико п ^ ш(Ъ — а) > ш(Ъ — а)/л. Таким образом, на каждом периоде функции Р(х)

Статья поступила в редакцию 30 декабря 2010 г.

Авторы выражают признательность за обсуждение данной работы и ценные замечания профессору Севастьянову Л.А.

1. Введение

(1)

необходимо брать много узлов сетки и высокой степени полином для аппроксимации, что невыгодно с любой точки зрения — это приводит к большому объёму вычислений [1].

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

ь

I хк е1шх&х

может быть проведено интегрированием по частям или с использованием соотношения

ъ

хкeinxdx =-L___ [Г(1 + к, —ша) — Г(1 + к, —шЬ)] ,

где Г — неполная гамма-функция.

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

ъ

J хк eiug(x)dx.

а

В этом случае все зависит от вида осциллятора д(х). В случае простой формы осциллятора можно применять метод Филона. Например, для полиномиальных функций вида д(х) = хг. Но даже для функций д(х) = х3 — х или д(х) = cos х метод Филона неприменим.

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

условию [р(х)ешд(х)] = f (х)ешд(х), позволяет вычислить осциллирующий интеграл

ь ь

I [f ] = J feiugdx = J^ [ре1шд] ёж = p(b)eiug(b) — p(a)eiug(a). (2)

a a

Можно переписать это условие в виде L [р] = f для оператора

L [р] = р' + гшд'р. (3)

Отметим, что метод не использует граничных условий, поскольку любое частное решение даёт решение задачи о значении определённого интеграла. Если удаётся хорошо аппроксимировать функцию р(х), то значение I [f ] легко вычисляется.

Рассмотрим оператор L(p), где р = 1 ск'Фк — разложение функции v по

векторам базиса {"ф\,..., фп]. При заданной системе узлов коллокации {xi,..., хп] коэффициенты Ск могут быть определены как решение системы уравнений метода коллокации:

L [p](xi) = f (xi), ..., L [p](xn) = f (xn). (4)

При определении приближенного значения интеграла I [/] в виде ь ь

Яь [/] = ! Ь(р)е1шдАх = I^ [ре1шд] ёж = р(Ьушд(Ъ) — р{а)е1шд(а) (5)

а а

справедлива следующая оценка ошибки аппроксимации [4]:

1) I [/] — 0Ь [/] = 0(ш-1) — в том случае, когда граничные точки не входят в число узлов сетки;

2) I [/] — 0Ь [/] = 0(ш-2) — в том случае, когда граничные точки входят в число узлов сетки.

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

2. Чебышевская дифференциальная матрица

Предположим, что нам известны значения полинома п-й степени в (п + 1)-й точках хо,... ,хп. Тогда эти значения определяют полином единственным образом и, следовательно, значения его производных р' (х) = ёр(х)/ёх в тех же точках. Значение производной в каждой точке может быть представлено в виде линейной комбинации значений полинома в этих точках. Эта зависимость может быть записана в матричной форме:

( р'(хо) \

/ ^0,0

^0,п \ ( Р(хо) \

\ &(Хп) ) V

й.

п,0

<1Г,

(6)

) \ Р(Хп) )

Матрица И = {dj,k} называется дифференциальной матрицей.

Если в качестве базисных функций выбираются полиномы Чебышева первого рода, а точки сетки являются узлами Гаусса-Лобатто, то элементы дифференциальной матрицы вычисляются по следующим формулам [5]:

(-1)

к — з

Xз Хк

й0к = 2

(-1)к

^кп о

1 - Хк 1 (-1)" — к

0 <3 = к <п,

2 1+ Хк '

0 < к <п, 0 <к <п,

^0п = 1 ( — 1)П,

, 1 Хк Лкк = — о

^0 = 6 (1 + 2п2), &

2 1 - хГ — 6 (1 + 2п2),

йк0=— 2т-Ь, (-1)™—к

1 + Хк

&пк — 2 ^п0 = — 1 ( — 1)П.

0 <к <п,

0 <к <п, 0 <к <п,

(7)

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

3. Универсальный метод квадратур

Если рассматривать интеграл на отрезке х € [а,Ь], то для перехода к стандартной области задания полиномов Чебышева [—1, 1] можно провести замену

Ь — а Ь + а г , переменных х = —-— ъ +---—, ъ Е [-1, 1]. Тогда

Р'(х) = ъ2~аР'

(8)

В соответствии с введённым линейным преобразованием узлы Гаусса-Лобатто

Ж]

' в исходных координатах имеют вид

Ъ = с™\Ы- 1

Ь — а (Ж]

■сой ' ^Т1

+ 3 = 0,1,...,М — 1.

Векторы значений функций и их производных в узлах Гаусса-Лобатто вычисляются по формулам

р = [р(хо),р(х1), ...,р(хп)] ,

т

$ = [р'(Хо),р'(хх), ...,р'(Хп)} , д1 = У (Х0),д'(Xl), ...,д'(Хп)]Т , / = [1 (х0),1 (Х\),...,1 (Хп)]Т .

(9)

Очевидно, что в соответствии с определением дифференциальной матрицы Чебышева мы можем записать р' в векторно-матричной форме из (8)

Р

2

а

Бр

(10)

В узлах сетки х^, уравнения (4) представим в виде

р' (хо) ' д'(хо)р(хо) " /(хо)

р' (жх) д' (Х1)р(Х1) ! (хл)

+ г =

р' (хм -1) д' (хм-1)р(хм-1) / (хм-1)

(11)

Подставим (9) и (10) в (11). Тогда система может быть записана в виде:

2

Бр + г ■ diag (д )р = /,

а

(12)

или

(Б + гЛ)р = X/, (13)

где Л = (Ь — а)/2, Л = diag (Хд'(х0), Хд'(х{),..., Хд'(хм-1)) является диагональной матрицей. Решение системы (13) содержит р(Ь) и р(а), а искомый интеграл вычисляется по формуле (5). Таким образом, вычисление интеграла (1) сводится к решению системы линейных уравнений (13).

X3 =

4. Метод решения системы линейных уравнений

Дифференциальная матрица О является сингулярной матрицей, но её число обусловленности улучшается при прибавлении Б к диагональной матрице гЛ. В этом случае ЦгЛ^2 показывает степень улучшения. Когда ||«Л^2 сравнительно велико, матрица О + становится хорошо обусловленной, в противном случае она остаётся плохо обусловленной, что ведёт к большой ошибке при решении (13). В

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

Пример. Рассмотрим в качестве примера вычисление функции Бесселя

-к

Jn(x) = je~i(nr~хsinтW

—-к

сотого порядка (п = 100) на интервале х £ [80,130] (рис. 1).

Ось Y 0.110.051-0.05-0.11 -

ВО 35 90 95 101 115 111 115 120 125

Ось X

Рис. 1. Значения функции Бесселя на интервале х G [80,130]

Вычисленные по описанному алгоритму значения функции сравнивались со значениями, определёнными по программе BesselJN известного пакета Cephes Math Library [6], автор Stephen L. Moshier. Максимальное отклонение по всем вычисленным значениям не превышает 1.0e-12 при вычислениях с двойной точностью.

5. Заключение

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

Литература

1. Приклонский В. И. Численные методы. Конспект лекций по курсу «Численные методы в физике». — М.: МГУ, 1999. — 146 с. [Priklonskiyj V. I. Chislennihe metodih. Konspekt lekciyj po kursu «Chislennihe metodih v fizike». — M.: MGU, 1999. — 146 s.]

2. Jeffrey G. B. Louis Napoleon, George Filon (1875-1937) // Obituary Notices of Fellows of the Royal Society. — 1939. — Vol. 2, No 7. — Pp. 500-509.

3. Levin D. Procedures for Computing One and Two-Dimensional Integrals of Functions with Rapid Irregular Oscillations // Math. Comp. — 1982. — Vol. 38, No 158. — Pp. 531-538.

4. Iserles A. On the Numerical Quadrature of Highly-Oscillatory Integrals I: Fourier Transforms // IMA J. Num. Anal. — 2004. — Vol. 24. — Pp. 1110-1123.

5. Mason J. C., Handscomb D. C. Chebyshev Polynomials. — Chapman and Hall/CRC, 2002. — 360 p.

6. www.netlib.org/cephes/.

UDC 519.644; 517.584

Integration of Highly Oscillatory Functions K. P. Lovetskiy, V. V. Petrov

Telecommunication Systems Department Peoples Friendship University of Russia Miklukho-Maklaya str., 6, Moscow, 117198, Russia

The paper demonstrates approximate methods of integral calculation of highly oscillatory functions. The paper describes a quadrature method which adopts Chebyshev differential matrix to solve the ordinary differential equation (ODE) and thus obtain integral values. This method make the system of linear equations well-conditioned for general oscillatory integrals. Furthermore, even if the system of linear equations is ill-conditioned, regularization method can be adopted to solve it properly and eventually obtain accurate integral results.

Key words and phrases: oscillatory integrals, Filon method, Levin method, Chebyshev differential matrix, ill-conditioned matrix, Bessel functions.