ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 7, с. 1118-1124
УДК 519.65
O ДЕФОРМИРОВАНИИ КОНТУРА ИНТЕГРИРОВАНИЯ В ФОРМУЛЕ ОБРАЩЕНИЯ ПРЕОБРАЗОВАНИЯ ЛАПЛАСА
© 2015 г. А. В. Лебедева, В. М. Рябов
(198504 Санкт-Петербург, Ст. Петергоф, Университетский пр-т, 28, СПбГУ) e-mail: natulich@bk.ru, victor.ryabov@mail.ru Поступила в редакцию 03.09.2014 г.
Предложены формулы обращения интегрального преобразования Лапласа с помощью деформирования контура интегрирования в формуле обращения Римана—Меллина, последующего применения квадратурных формул и получения оценок погрешности. Библ. 18. Фиг. 2.
Ключевые слова: преобразование Лапласа, обращение преобразования Лапласа, квадратурные формулы, сходимость квадратурных формул.
DOI: 10.7868/S0044466915050142
1. ВВЕДЕНИЕ
Интегральное преобразование Лапласа F(z) функции-оригиналаf(t),
да
F(z) = Je-tf( t) dt,
0
представляет собой мощный инструмент для решения широкого класса прикладных задач математической физики. Одним из его главных достоинств является алгебраизация процедур математического анализа, с помощью которой удается свести интегральные и дифференциальные уравнения к более простым.
Как правило, при решении задач операционными методами наиболее трудным этапом является процесс обращения, т.е. определение оригинала по его изображению. Существуют таблицы соответствия функций-оригиналов и их изображений (см., например, [1]), теоремы разложения, формула обращения Римана—Меллина, позволяющие теоретически точно находить оригинал. Но решение практических задач часто приводит к изображениям, к которым не могут быть применены эти классические приемы обращения. Следовательно, возникает необходимость разработки и применения приближенных методов.
Возможные подходы к задаче обращения и их реализация описаны в [2]. Обзор других способов обращения, не вошедших в [2], приведен в [3] и [4], [5]. Теоретические основы операционного исчисления содержатся в классических работах [1], [6], [7], [8]. Вопросам приложения операционного исчисления к решению прикладных задач, среди прочих, посвящены фундаментальные труды [9], [10].
Не существует универсального метода обращения, дающего удовлетворительные результаты для произвольного изображения F(z). Любой конкретный метод обращения должен учитывать специфику поведения изображения (или функции-оригинала), что прежде всего находит отражение в выборе подходящих систем функций в пространствах оригиналов и изображений, с которыми легко работать и с помощью которых могут быть хорошо приближены заданные образы и оригиналы. Выбору методов обращения в зависимости от способа задания информации об изображении искомого оригинала и их описанию посвящена работа [3].
Цель настоящей работы состоит в рассмотрении методов обращения с помощью деформирования контура интегрирования в формуле обращения Римана—Меллина, последующего применения квадратурных формул и получения оценок погрешности.
2. ДЕФОРМИРОВАНИЕ ЛИНИИ ИНТЕГРИРОВАНИЯ В ФОРМУЛЕ ОБРАЩЕНИЯ
В этом пункте исследуется метод вычисления интеграла Римана—Меллина
c + i да
f(t) = -П Г eztF(z)dz, (1)
2 я i J
c - i да
задающего обращение преобразования Лапласа, с помощью замены линии интегрирования и последующего применения специальных квадратурных формул.
Напомним, что интеграл (1) понимается в смысле главного значения, он не зависит от c, и в случае разрыва оригинала в точке t мы получаем полусумму предельных значений оригинала слева и справа в точке t.
Положим в формуле обращения (1) г = c + ix, тогда exp(zt) = ехр(сОехр(г'хО. При фиксированном t первый сомножитель постоянен, а второй пробегает единичную окружность на комплексной плоскости бесконечное число раз. К тому же в случае с > 0 с ростом t первый сомножитель неограниченно возрастает, так что попытка приблизить интеграл в (1) римановыми суммами вряд ли приведет к цели.
С целью уменьшения осцилляции сомножителя ехр(/тО и ограничения скорости роста сомножителя exp(ct) заменим линию интегрирования в (1) эквивалентным контуром L, начинающимся и заканчивающимся в левой полуплоскости так, что Re(z) —► —да на обоих его концах, внутри которого содержатся все особенности изображения F(z).
Предположим, что изображение F(z) удовлетворяет следующей лемме.
Лемма Жордана (см. [2]). Если на некоторой последовательности дуг окружностей CR :
|z| = Rn, Rez < a, Rn—► да, а фиксировано,
функция F(z) стремится к нулю равномерно относительно arg z, то для любого положительного t справедливо предельное соотношение
lim Г F(z) eztdz = 0.
n ^ да J
Ч
В этих предположениях интеграл (1) представим в виде
c + i да
f(t) = -П Г ez'F(z)dz = -П. \ezt F(z)dz. (2)
2ni J 2яiJ
c - i да L
В [11] был предложен контур
L = J z|z = 2nui u e [ i i] l 1 1 - exp(-2nui) ' J
Он состоит из двух симметричных относительно вещественной оси плоскости z ветвей, исходящих из точки z(0) = 1 налево и при u —»- ±1 стремящихся к асимптотам z = ±ni. После замены переменной в интеграле (2) получаем
1
f(t) = -1- I"exp(z(u)t)F(z(u))z'(u)du.
2 П i
2 п/
-1
Для приближенного вычисления этого интеграла в [11] применяется формула трапеций.
Заметим, что при малых и вещественная часть функции z(u) положительна и с ростом I сомножитель ехр(г(и)0 неограниченно растет по модулю, в отличие от остальных значений и. Эту особенность поведения интегрируемой функции надо учитывать при численном интегрировании. В [12] предлагается в представлении
яО = тП-^Дг!ойг, (> 0, (3)
2 П 1И
L
использовать контур
Ь = {z\z = I(и), и е(-да, +да)}, I(и) = 2 - еИ(и) - ф&(и).
Этот контур состоит из двух симметричных относительно вещественной оси плоскости г ветвей, исходящих из точки г(0) = 1 налево и при и —»- ±да стремящихся к асимптотам г = (ц — параметр контура).
В результате формула (3) принимает вид
да
Я О = 1О,( и) йи, (4)
—да
где
О, ( и) = -П-ехр (1( и)) Г( и) Д 1( и)/0. (5)
2я 11
Функция (5) не имеет особенностей как на линии интегрирования, так и в некоторой "полосе", содержащей внутри себя линию интегрирования. Как показано в [13], формула трапеций для вычислений интеграла (4) может давать хорошую точность, при этом скорость сходимости зависит от ширины полосы и шага численного интегрирования. К сожалению, приведенные выше две линии интегрирования оказываются неудачными — для них полоса имеет переменную ширину регулярности интегрируемой функции и даже стремится к нулю при неограниченном возрастании модуля переменной интегрирования.
Этого недостатка лишены предложенные в [14] параболические и гиперболические контуры — в плоскости переменной и им соответствуют полосы постоянной ширины. Для них погрешность применения формулы трапеций для вычисления интеграла (4) убывает при увеличении ширины полосы регулярности функции (5), порожденной изображением F(z) и выбранным контуром интегрирования 1(и).
3. НОВЫЕ КОНТУРЫ ИНТЕГРИРОВАНИЯ В ФОРМУЛЕ ОБРАЩЕНИЯ
Если все особые точки изображения расположены в левой полуплоскости, то в формуле обращения (1) можно положить с = 0. В противном случае вместо F(z) будем работать с функцией F(z + а), а > 0 (соответствующей функции-оригиналу/(?)ехр(—Ш)) с особыми точками в левой полуплоскости.
Итак, считаем, что в формуле обращения (1) имеем с = 0.
В качестве контуров Ь в формуле (2) возьмем приведенные на фиг. 1, 2 кривые, состоящие из прямолинейных участков Ь1, Ь2, Ь3.
На участке Ь2 для обоих контуров параметр г изменяется по закону г = хЫ, х е [—1, 1], Ь > 0.
Очевидно, что первый контур является частным случаем второго (если углы наклона участков Ь1, Ь3 второго контура равны нулю), тем не менее, целесообразно его рассмотреть отдельно.
Ь3
Ь2
Ь1
-Ы
Фиг. 1. Первый контур.
Ы
0
Фиг. 2. Второй контур.
1. Случай первого контура. На участке L1 имеем z = х — bi, x е (—да, 0], а на участке L3 z = —х + bi, х е [0, да). Равенство (2) записывается в виде
f{ 0 = 2Л- \ezF{z) dz = ¿1 1+ 1+ i J ez'F(z) dz. (6)
После замен переменной z получаем
да
ieztF(z)dz = e~'bt Je~xtF(- x - bi)dx,
L, 0
1
JeztF(z)dz = bi JeibtxF(ibx)dx, (7)
-i
да
г) йг = -е1" - х + Ы) йх.
о
Опишем эффективные методы вычисления интегралов (7). Начнем с интегралов вида
1
jвlaxg(x) йх. (8)
-1
Число ю = Ы может быть большим, так что подынтегральная функция в (8) сильно колеблющаяся, что затрудняет вычисление интеграла с помощью стандартных квадратурных формул.
Воспользуемся предложенным в [15] способом приближенного вычисления интегралов вида (8): дана квадратурная формула Гаусса (см. [16]) с п узлами алгебраической степени точности 2п — 1:
1
п
|^{х)йх - ^АV(Xj). (9)
1 j = 1
Ее узлы — корни многочлена Лежандра Рп(х).
Построим интерполяционный многочлен Ьп-1(х) для функции g(х) по значениям g(х1), ..., g(хn) и представим его в виде разложения по многочленам Лежандра:
п - 1
£п - 1(х ) = ^ СкРк {х ),
L
коэффициенты которого точно вычисляются с помощью формулы Гаусса (9):
п
2 к + 1-
ск =
£ (X) Рк (X) .
2
1 = 1
Погрешность интерполяции представима в виде
g(x) - Ьп 1 (х) = ¿-ш Рп (X), ^ 6 [-1, 1 ]. (10)
п! (2п)!
Полагая в формуле (8) g(х) ~ Ь„ _ х(х) и используя известное равенство
1
¡е*хРк(х)йх = I ^Jk +1/2(га),
1 П г— П-1
]У^(х)йх * £Б^га^), Я/га) = А^ £ 2к+-1 ^к +1/2(га)Рк(х,), (11)
М*) = ¡е"**^(х) - Ln- 1(х))dx
-1
где ,Д,(га) — функция Бесселя, придем к интерполяционной квадратурной формуле
1 п I- п -1
, !а х ч , ( ^ ^ ч ^ Ч ^Ч .12 П
-1 1 = 1 'к = 0 погрешность которой равна
1 -1
и в силу преставления (10) оценивается неравенством
К< (х)|. (12)
(2 п -1 )!!Т2пЛ [-1, и
В роли g(х) выступает функция Ш(1Ьх). Для вычисления коэффициентов Dj(гa) нужны значения функций Бесселя. Один из способов их вычисления указан в [15]. Современные математические пакеты (например Мар1е) позволяют находить их без всяких затруднений. Эффективный алгоритм нахождения корней многочленов Лежандра и построения формулы Гаусса (9) приведен в классической работе [17].
В случае недостаточной точности оценки (12) для избранного п можно поступить двумя способами: 1) построить формулу Гаусса с большим числом узлов; 2) применить составную квадратурн
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.