О вычислении интеграла Римана-Меллина Текст научной статьи по специальности «Математика»

О ВЫЧИСЛЕНИИ ИНТЕГРАЛА РИМАНА—МЕЛЛИНА*

Н. И. Порошина1, В. М. Рябов2

1. С.-Петербургский государственный университет, аспирант, [email protected]

2. С.-Петербургский государственный университет, д-р физ.-мат. наук, профессор, [email protected]

1. Введение. Интегральное преобразование Лапласа

где функция Г (г) — изображение, / (£) — оригинал, представляет собой мощный инструмент для решения широкого класса прикладных задач математической физики. Одним из его главных достоинств является алгебраизация процедур математического анализа, с помощью которой удается свести интегральные и дифференциальные уравнения к более простым. Кроме того, изображение Лапласа является аналитической функцией в некоторой полуплоскости И,е г > А, что позволяет привлечь к исследованию решаемой задачи результаты теории функций комплексного переменного.

Как правило, при решении задач операционными методами наиболее трудным этапом является процесс обращения, т. е. определение оригинала по его изображению. Существуют таблицы соответствия функций-оригиналов и их изображений [1], «теоремы разложения», формула обращения Римана—Меллина, позволяющие точно или приближенно находить оригинал. Но решение практических задач приводит к изображениям, к которым не могут быть применены эти «классические» приемы обращения. Следовательно, возникает необходимость разработки и применения приближенных методов.

Наиболее полно возможные подходы к задаче обращения и их реализация описаны Крыловым и Скобля [2]. Обзор других способов обращения, не вошедших в [2], приведен в статье [3].

В настоящей работе мы исследуем метод вычисления интеграла Римана—Меллина

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

2. Деформирование контура интегрирования. Напомним, что интеграл (2) понимается в смысле главного значения, он не зависит от с, и в случае разрыва оригинала

* Работа выполнена при финансовой поддержке РФФИ (грант №08-01-00285).

© Н. И. Порошина, В. М. Рябов, 2009

(1)

с > 0,

(2)

в точке £ мы получаем полусумму предельных значений оригинала слева и справа от точки £.

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

С целью уменьшения осцилляций сомножителя ехр(гт£) и ограничения скорости роста сомножителя ехр(с£) заменим линию интегрирования в (2) эквивалентным контуром Ь, начинающемся и заканчивающемся в левой полуплоскости так, что Ие (г) ^ —то на обоих его концах. Такая замена возможна, если выполнены два условия:

1) внутри контура Ь содержатся все особенности изображения Г (г);

2) |Г(г)| ^ 0 равномерно в полуплоскости Ие (г) < 70 при |г| ^ то (70 — абсцисса сходимости интеграла (1)).

Далее всюду предполагается, что эти условия выполняются.

Телботом [4] был предложен контур

Он состоит из двух симметричных относительно вещественной оси плоскости г ветвей, исходящих из точки г(0) = 1 налево и при и ^ ±1 стремящихся к асимптотам г = ±пі.

Так как кривая Ь содержит внутри себя особые точки изображения и выполнены условия леммы Жордана [5],

Для приближенного вычисления этого интеграла в статье [4] применяется формула трапеций. Заметим, что при малых и с ростом £ сомножитель ехр(г(и)£) неограниченно растет. Однако это затруднение легко преодолевается предварительной заменой переменной вида г\ = в интеграле (2), в результате чего экспоненциальный сомножитель становится ограниченным и при и ^ ±1 быстро стремится к нулю.

В статье [6] предлагается в представлении

Этот контур состоит из двух симметричных относительно вещественной оси плоскости г ветвей, исходящих из точки г(0) = 1 налево и при и ^ ±то стремящихся к асимптотам г = ±*^ — параметр контура (5)).

В результате формула (4) принимает вид

(3)

где

Сг{и) = т-1—ехр (1(и)1'(и))Р(1(и)/г). (7)

2п*

Функция (7) не имеет особенностей как на линии интегрирования, так и в некоторой «полосе», содержащей внутри себя линию интегрирования. Как показано в работе [7], формула трапеций для вычислений интеграла (6) будет давать хорошую точность, при этом скорость сходимости зависит от ширины полосы и шага численного интегрирования [7].

Сформулируем полученный в работе [7] результат в удобном для дальнейших целей виде: пусть т = и + «V, и, V € Д, функция $(и>) аналитична в полосе — ! < V < с для некоторых с > 0, ! > 0, и при |и>| ^ то равномерно в этой полосе $(и>) ^ 0 с такой скоростью, что существуют интегралы

/Ж

|$(и + «г)| !и, г € [—!, с].

-Ж

Для вычисления интеграла

/Ж

^(и) !и

-Ж

применим формулу трапеций как с бесконечным числом узлов

Ж

1н = н ^2 #(&н),

к= — ж

так и с конечным числом 2^ + 1 узлов

N

=

1hJN = Н ^ #(*Н).

к= — N

Положим ДЕ = |1 — Д|, ТЕ = |/^, — |.

Теорема. Пусть для некоторых положительных чисел М+, М- справедливы нера-

венства

/ОО рЖ

|^(и + «г)| !и < М-, г € [—!, 0], / |$(и + «в)| !и < М+, в € [0, с]. (8)

-ОО </ —ОО

Тогда

где

|1 — Д| < ДЕ+ + ДЕ-,

м+ м_

ВЕ< =_____________-______ ВЕ- =_________________________. (9)

ехр(27г с/К) — 1 ’ ехр(2тгг1/к) — 1

Эта теорема без доказательства приведена в работе [8]. Для случая с = ! оно имеется в статье [7].

Константы М+,М-, входящие в оценки (8), не зависят от шага численного интегрирования Н и значения N, и потому вместо оценок (9) в дальнейшем ограничимся их качественным поведением с помощью соотношений

ДЕ+ = 0(ехр(—2пс/Н)), ДЕ- = 0(ехр(—2п!/Н)), Н ^ 0. (10)

Для погрешности ТЕ замены бесконечной суммы /^ метода трапеций суммой /ь,^, как и в работе [8], при постоянном Н полагаем

ТЕ = О(|0^)|), N ^то. (11)

Далее, для избранного контура интегрирования естественно исходить из равенства характеристик (10), (11) величин ДЕ+, ДЕ-, ТЕ, что приведет к некоторому способу выбора параметров контура и шага интегрирования в зависимости от N и возможности получения полной погрешности метода.

К сожалению, контуры (3) и (5) достаточно сложны для реализации полученных в [7] оценок погрешности формулы трапеций для вычисления интегралов вида (6).

3. Параболический контур. Рассмотрим кривую

г = ^(*ад + 1)2, w = и + «V, ^> 0. (12)

При V =0 получаем параболу

г(и) = ^(1 — и2) + 2*^и. (13)

Возьмем ее в качестве линии интегрирования, тогда формула обращения (интеграл Римана—Меллина) примет вид

1 (Ж

/(£) = -- ехр(г(и)Ь)Р(г(и))г (и) <1и. (14)

2п* ./-Ж

Прямые т = и — И, ! > 0, и т = и + «с, 0 < с < 1, при отображении (12) переходят, соответственно, в параболы

г (и) = м((1 + !)2 — и2) + 2*^и(1 + !), г(и) = ^((1 — с)2 — и2) + 2*^и(1 — с),

первая из которых расположена правее линии интегрирования (13), а вторая — левее (и при с =1 она вырождается в отрицательную полуось и < 0).

Предположим, что все особые точки изображения расположены в конечной части полуплоскости Ие г < 0 и не все они вещественны. Очевидно, при некотором ц все они попадут внутрь параболы (13) и, тем более, внутрь параболы г(и) = ^((1 + !)2 — и2) + 2*^и(1 + !) при любом ! > 0. С ростом ! знаменатель дроби в формуле (9) для ДЕ- растет, но и числитель М- тоже возрастает. В работе [8] показано, что существует оптимальное значение параметра !, при котором достигается наилучшая оценка

ДЕ- = 0(ехр(—п2/(^Н)2 + 2п/Н), Н ^ 0.

В работе [8] рассмотрен случай, когда все особые точки изображения расположены на полуоси и < 0. Следовательно, наилучшим значением параметра с, при котором полоса аналитичности подынтегральной функции в (14) будет самой широкой, равно единице. В нашем предположении ветви параболы не смыкаются, и наибольшее допустимое значение с удовлетворяет неравенству 0 < с < 1. Далее, критерий (11) дает

ТЕ = 0(| exp(z(uN)£)|) = 0(ехр(уи,£(1 — uN))) = 0(ехр(и£(1 — (НN)2))), N ^ то.

Приравнивание показателей величин ДЕ-, ДЕ+, ТЕ приводит к соотношению

2п п2 2п

----1~с = —-

Н ^Н2 Н

из которого находим

ь=^ЦТ4ф+Т)) * =--------/ " (15)

N ^ 2{с+1)Ш 2(с+ 1)^/1 + 4с(с+ 1) *

Напомним, что число с здесь не произвольно. Например, в случае /(£) = вт(^£) изображение Е(г) = ^/(г2 + ^2) имеет особые точки ±гш, которые должны лежать между точками пересечения параболы г(и) = м((1 — с)2 — и2) + 2г^и(1 — с) с мнимой осью, т. е. должно выполняться неравенство 2^,(1 — с)2 > |^| или равносильное ему неравенство

М> Н^+1Ь/1 + 4с(с+1)^ п(1 — с)2

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

° (“Р (~|С)) = ° (“Р (~ ^1 + 4^+1)")) ’ (1?)

Замечание 1. Из оценки (17) при с =1 получается результат статьи [8].

Замечание 2. Зависимость подынтегральной функции от £ в формуле (14) в соответствии с выражением для ц в (15) целиком сосредоточена в изображении Е(г) и, фактически, означает указанную выше замену переменной вида = г£, приводящую интеграл (14) к форме (4).

Замечание 3. Разумеется, в оценке (16) можно устремить с к единице в надежде получить более точный результат, но при этом резко возрастет допустимая нижняя граница числа N в соответствии с неравенством (16), что приведет как к возрастанию объема вычислений, так и к необходимости увеличения точности вычислений.

Замечание 4. В работе [9] для обращения преобразования Лапласа с особенностями на отрицательной полуоси (в частности, при наличии в начале координат точки ветвления) по существу предлагается замена вертикальной линии интегрирования на предельную параболу г(и) = м((1 — с)2 — и2) + 2г^и(1 — с) при с =1, что соответствует наихудшему выбору «ширины» полосы аналитичности интегрируемой функции. Правда, в работе [9] в результате исходный комплексный интеграл сводится к интегралу по вещественной полуоси от вещественной функции, о приближенном вычислении которого речи вовсе не идет.

4. Гиперболический контур. Рассмотрим кривую

г = ^,(1 + 8т(*'Ш — а)), т = и + «V, (18)

определяемую вещественными параметрами ^, а (^ > 0, 0 < а < п/2).

Прямые т = и + «с при отображении (18) переходят в гиперболы

-4^ЧУ-(^^У=М2, г = х + гу. (19)

вт(а + с) / \ео8(а + с) у

В качестве линии интегрирования возьмем левую ветвь гиперболы, получающейся из (19) при с = 0. При возрастании с ветви гиперболы сближаются, и при с = п/2 — а она переходит в отрицательную полуось. При убывании с гипербола расширяется и

при с = —а она переходит в вертикальную прямую. В работе [8], как и для параболического контура, рассмотрен случай принадлежности всех особых точек изображения отрицательной полуоси. В такой ситуации наибольшая ширина полосы аналитичности подынтегральной функции в (14) достигается, если границам этой полосы соответствуют крайние гиперболы при с = п/2—а (отрицательная полуось) и с = —а (вертикальная прямая) (см. [8]).

При наличии невещественных особых точек границам полосы аналитичности соответствуют две гиперболы: одна при некотором с € (0, п/2 — а) и другая при с = —а (вертикальная прямая). Повторяя рассуждения работы [8], в этом случае придем к следующим асимптотическим выражениям ошибок:

ДЕ+ = 0(ехр(—2пс/^)), ДЕ_ = 0(ехр(^£ — 2па/Н)), Н ^ 0,

ТЕ = 0(ехр(^£(1 — вт(а) ch(НN)))), N ^ то.

Приравняем порядки этих величин:

—2пс/Н = ^ — 2па/Н = ^(1 — вт(а) ch(НN)).

Отсюда получаем

А(а) N 2п(а — с) 1 ( а \

и=-------- —г—А(а) = сЬ ( ------------ --— . 20

N t А(а) \ (а — с) 81п(а) у

Для существования А(а) и положительности ц необходимо и достаточно выполнения неравенства с < а. Параметры а и с должны удовлетворять неравенству 0 < с < шш{а,п/2 — а}. В частности, при а = п/4 оно заведомо выполняется, так как при с = п/4 соответствующий контур превращается в отрицательную полуось.

Например, в случае f (£) = 8ш(^) изображение Е(г) = ш/(,г2 + ш2) имеет особые точки ±«ш, которые должны лежать между точками пересечения левой ветви гиперболы (19) с мнимой осью, т. е. должно выполняться неравенство ^сов2(а + с)/ вш(а + с) > |ш| или равносильное ему неравенство

N 2п(а — с) сов2 (а + с)

--------------:—-----г > М-

£ А(а) 8т(а + с)

Очевидно, для любых допустимых значений а и с оно выполняется для достаточно больших N.

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

O(exp(—2пс/Н)) = O(exp(—B(а)N)), В(а) = 2пс/А(а).

Пример. Дано изображение Е(г) = г/(г2 + ш2) оригинала /(£) = cos(wt). Пусть ш = 2, £ = 5. Возьмем с = 0.3. Для параболического контура в соответствии с неравенством (16) получаем N = 14, и приближенное значение отличается от точного на

2 • 10-7. Для гиперболического контура возьмем с = 0.3, а =1. Из неравенства (20) находим наименьшее допустимое значение N = 35, и отличие приближенного значения от точного составляет 4 • 10-9.

Литература

1. Диткин В. А., Прудников А. П. Интегральные преобразования и операционное исчисление. М., 1961. 524 с.

2. Крылов В. И., Скобля Н. С. Методы приближенного преобразования Фурье и обращения Лапласа. М., 1974. 224 с.

3. Порошина Н. И., Рябов В. М. Об обращении преобразования Лапласа некоторых специальных функций // Вестн. С.-Петерб. ун-та. Сер. 1. 2009. Вып. 3. С. 50-60.

4. Talbot A. The accurate numerical inversion of Laplace transform // J. Inst. Maths. Applies. 1979. Vol. 23. P. 97-120.

5. Лаврентьев М.А., Шабат Б. В. Методы теории функций комплексного переменного. М., 2002. 688 с.

6. Рябов В. М. Об одном способе обращения преобразования Лапласа // Кубатурные формулы и их приложения. Материалы X международного совещания-семинара. Улан-Удэ. 24-28 августа 2009. Улан-Удэ: Изд-во ВСГТУ, 2009.

7. Самокиш Б. А. Замечание о вычислении определенных интегралов // Методы вычислений. Вып. 2. Л., 1963. С. 45-49.

8. Weideman J.A. C., Trefethen L. N. Parabolic and hyperbolic contours for computing the Bromwich integral // Math. Comput. 2007. Vol. 76. N259. P. 1341-1356.

9. Bobylev A.V., Cercignani C. The inverse Laplace transform of some analytic functions with an application to the eternal solutions of the Boltzmann equation // Applied Mathematics Letters. 2002. Vol. 15. P. 807-813.

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