Специальные квадратурные формулы обращения преобразования Лапласа Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 11, Специальный выпуск, 2006

СПЕЦИАЛЬНЫЕ КВАДРАТУРНЫЕ ФОРМУЛЫ ОБРАЩЕНИЯ ПРЕОБРАЗОВАНИЯ ЛАПЛАСА*

В. М. Рябов

Санкт-Петербургский государственный университет, Россия e-mail: Riabov@VR1871. spb. edu

Methods for the Laplace transform inversion by means of quadrature formulas are described. Some quadrature formulas specially adjusted for recovering of a slowly changing inverse transform with approximation in the form ts-1P2n-i (ta), where s > 0, 0 < a < 1, P2n-1 (t) is polynomial and n is the number of formula's nodes are proposed.

Введение

Задача обращения интегрального преобразования Лапласа состоит в нахождении решения уравнения

х

F (p) = / е-/М

о

в котором F(p) — известное изображение; /(x) — искомый оригинал. Для простоты будем считать, что функция F(p) регулярна в полуплоскости Re(p) > 0, чего всегда можно добиться домножением оригинала на соответствующую экспоненту. Как правило, точное обращение осуществить не удается, и потому возникает необходимость разработки и применения приближенных методов. Наиболее полно возможные подходы к задаче обращения и их реализация описаны в [1]. В первую очередь следует назвать построение квадратурных формул для приближенного вычисления интеграла Римана — Меллина

e+ix

f(t) = (L~1F)(t) = -^ J e^F{p) dp, с> 0,

е—ix

задающего обращение преобразования Лапласа.

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00984).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

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

1. Построение квадратурных формул обращения

Предположим, что при некотором 5 > 0 функция ^Др) = (р) регулярна в полуплоскости Ее(р) > 0. Запишем формулу обращения Римана — Меллина в виде

е+гте е+гте

1С I3-1 г

■П1)=2ш ] у (1)

е-гте е-гте

Выберем произвольные попарно различные числа р1,..., рп в иолуплоскости Ее(р) > 0 и построим интерполяционную квадратурную формулу (ИКФ) вида

е+гте

1 7 п

— / ерр~3ф)йр^^Акфк), (2)

к—1

е-гте к—1

точную для функций ^(р) = р-, ] =0,1,..., п — 1, что равносильно выполнению условий

п1

= ¿ = (з)

к—1

Очевидно, система (3) однозначно разрешима.

Теперь потребуем, чтобы ИКФ вида (2) имела наивысшую степень точности 2п — 1, т. е. чтобы равенства (3) выполнялись для ] = 0,1,..., 2п — 1.

Этими условиями квадратурная формула наивысшей степени точности (КФНСТ) вида (2) определяется однозначно [1]: ее узлы рк суть корни уравнения рПз)(1/рк) = 0, где

РП3)(х) = ¿( — 1)п-к( П) (п + 5 — 1)к Хк к—0 ^ '

1, к = 0,

ак =

а(а +1) ••• (а + к — 1), к > 1.

Коэффициенты КФНСТ определяются из условия интерполяционности (3).

Применяя КФНСТ вида (2) для вычисления последнего интеграла в формуле (1), получаем КФНСТ обращения преобразования Лапласа

п

/(*) « Ак^з(ркА), *> 0. (4)

к—1

По построению она точна для функций /(¿) = , ] = 0,1,..., 2п — 1.

Итак, сначала находятся узлы КФНСТ из уравнения Pis)(1/pk) = 0, а затем определяются ее коэффициенты из системы (3).

Выбор в качестве s максимально возможного положительного числа, равного скорости убывания изображения на бесконечности, позволяет точно описать поведение оригинала при t ^ +0 то при t ^ ж приближенное решение никаким априорным условиям удовлетворять не будет.

Значения искомого оригинала f (+0) f (+ж) в предположении их существования можно вычислить по формулам

f (+0) = lim pF(p), f (+ж) = lim pF(p).

Квадратурная формула наибольшей степени точности (4) при t ^ 0 и t ^ ж приводит к верным результатам лишь для s = 1.

2. Построение обобщенных квадратурных формул

Описанные выше квадратурные формулы обращения обладают высокой степенью точности, если искомый оригинал хорошо приближается функциями вида ts-1Q2n-1(t). Однако такое предположение не всегда справедливо.

Так, в задачах линейной вязкоупругости [2], описывающих напряженное состояние на основе определяющего соотношения Больцмана — Вольтерра (пространственные координаты для простоты опущены)

e(t) = | + Р J K(t- т)а(т) dr j , (5)

сравнительно просто находятся изображения по Лапласу решений (деформаций e(t) и напряжений a(t)) из соответствующих интегральных уравнений вида (5). Сами решения медленно изменяются во времени и допускают хорошие приближения вида ts-1Q(ta) при 0 < a < 1, где Q(t) — некоторый многочлен. Такая ситуация характерна для длительных процессов деформирования [2, 3]. Изображения таких функций имеют вид p-sR(p-a), где R(x) — многочлен той же степе ни, что и Q(x).

Формально для обращения изображений такого вида можно использовать вышеописанные формулы обращения наивысшей степени точности, полагая, что a ~ 1.

Ядро K должно иметь интегрируемую особенность в точке t = 0 [2]. Чаще всего в качестве такового берут дробно-экспоненциальную функцию Работнова [2]

Способ определения параметров дробно-экспоненциальной функции по измеренной функции ползучести описан в работе [3].

Интеграл от этого ядра по полуоси t > 0 должен быть конечным, для чего необходимо в < 0. Не умаляя общности, далее считаем, что в = -1, и пусть символ 3a(t) означает

Эа(-1,t). "

В наследственной механике твердого тела наряду с функцией (6) широко используется и интеграл от пее с переменным верхним пределом [2, !|. Для облегчения использования этих величин составлены таблицы функций [2, 4]

^(а,х) = ГаЭа(г), (а,х) = Г"-1/ Эа(т) ¿т, х = га+1

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

Изображения по Лапласу функции Эа(£) и интеграла от нее равны соответственно

1 1

а = 1 + а.

ра + Г р(ра + 1)'

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

Для устранения этого недостатка следует использовать обобщенные квадратурные формулы наивысшей степени точности (ОКФНСТ) вида

с+гм

1 7 п

— / ерр~3ф)йр^^Акфк), (7)

к—1

с-гм к—1

точные для функций ^(р) = р-а\ ] = 0,1,..., 2п— 1, пли для оригиналов вида ¿5-1^2п-1 (¿а). Такие формулы введены в работе [5] и подробно исследованы в статье [6]. В статье [6] доказана

Теорема 1. Для того чтобы, формула (7) была точна для функций р(р) = р-а\ ] = 0,1,..., 2п — 1, необходимо и достаточно выполнение двух условий:

— формула (7) — интерполяционная;

— построенный по узлам, формулы (7) многочлен

п

ип(х) = П (х — Р-а

к—1

удовлетворяет условиям

с+гм

врр-3ип(р-а)р-ат ¿р = 0, т = 0,1,..., п — 1.

г

с- гм

Доказано [6], что такой многочлен существует и определяется однозначно, а его корни, т.е. узлы ОКФНСТ, удовлетворяют неравенствам 11е (р а) > 0 к = 1, 2,..., п. Отметим, что узлы и коэффициенты формулы (7) суть комплексные числа. Эти формулы оказались весьма эффективными при решении задач линейной вязкоупругости [3].

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

Приведенные выше преобразования Лапласа функции Эа(£) и интеграла от нее не имеют особенностей на комплексной плоскости С \ Я- с разрезом вдоль полупрямой

Я- = {р € С : 1т(р) = 0, 11е(р) < 0}.

Наши изображения F(р) фактически зависят от ра, т. е. F(р) = ФДр"). Введем в рассмотрение функции

F±(г) = ФДГехр(±гап)), г > 0.

Очевидно, F + (t) = F- (t) в силу вещественности функции-оригинала. Воспользуемся получеппым в работе [7] следующим результатом. Теорема 2. Пусть выполнены условия:

(A) F(p) = o(l) при |p| ^ F(p) = o(|p|-1) при |p| ^ 0 равномерно в любом секторе | arg p| < п — п п > п > 0;

(Б) существует £ > 0 такое, что для любого у, удовлетворяющего неравенству п — £ < у < п, справедливы, соотношения

F(rexр(±»у)) g |ir(rexp(±¿y>))| < а(г)>

gde a(r) не зависит от у и a(r) exp(—ir) G L1 (R+) для, любо го ö > 0. Тогда

те

/(ж) = (L~1F){x) = i J e~xtlmF-(t)dt. (8)

0

Пусть F (p) = 1/(ptt + 1), тогда

1 ttt sin па

ImF"(í) = Im-

ta exp(-ina) + 1 1 + 2ta cos na + t2a так что выполнены все условия теоремы 2 и формула (8) дает

те

^ . . sin na [ ta dt

Э„(х) =- / e

= xa-1

п J 1 + 2ta cos па + t2a 0 те

sin па f zae-z dz

(9)

п J z2a + 2zaxa cos па + x2a 0

При x ^ 0 последний интеграл в (9) стремится к величине Г(1 —а), и с учетом формулы Г(а)Г(1 — а) = п/ sin па при х ^ 0 го представления (9) получаем Эа(х) ~ ха-1/Г(а), что совпадает с первым членом ряда (6). Положим

x

ga(x)= / 3a(t) di. (10)

Изображение этой функции равно

Ga(p) =

p(pa + 1)

Для нее не выполняется условие (А) теоремы 2 (при p ^ 0 величина |Ga(p)| слишком быстро возрастает).

Представим Ga(p) в виде

1 pa-1

= - - 4-7 (n)

p pa + 1

и положим

pa- 1

Q°(P) = ^TT- (!2)

Обозначим через qa(x) функцию-оригинал с изображением (12). Формула (11) означает, ЧТО ga(x) = 1 - qa (x).

Изображение (12) удовлетворяет условиям теоремы 2, для него находим

л ,, т ta-1 exp(-in(a — 1)) ta-1 sin na

Im F~(t) = Im- v

ta exp(-ina) + 1 1 + 2ta cos na + t2a

Следовательно,

с

, , sin na f -xt ta-1 dt

qa(x) =- / e

xa

n J 1 + 2ta cos na + t2a 0

сс

sin na [' za-1 e-z dz

(13)

n j z2a + 2zaxa cos na + x2a 0

Из определения (10) находим ga(0) = 0, поэтому необходимо qa(0) = 1. Подставив x = 0 в первый интеграл в представлении (13) и сделав замену ta = z, придем к табличному легко вычисляемому интегралу и таким образом убедимся в справедливости равенства qa(0) = 1 при всех a > 0.

Итак, паши задачи обращения преобразования Лапласа свелись к вычислению интегралов

со

эа(х) = х^ ™ [ --^^-г; (14)

w n J z +2zaxa cos na + x2a v ;

0

со

, , a sin na f za-1 e-z dz , , qa{x) = x - / -—. (15)

n J z2a + 2zaxa cos na + x2a 0

Для их приближенного вычисления можно применить квадратурные формулы типа Гаусса [8] с весом Лагерра zae-z для первого интеграла и с весом za-1 e-z для второ-

a

умепьшаться. Поэтому для приближенного вычисления интегралов (14) и (15) построим обобщенные квадратурные формулы вида

zee-z f (z) dz Ak f (zk), в> -1, (16)

k=1

с

точные для функций / (г) = гат, т = 0,1..., 2п — 1.

Теорема 3. Для, того чтобы, формула (16) была, точна для, функций /(г) = гат, т = 0,1..., 2п — 1; необходимо и достаточно выполнение двух условий:

— формула (16) — интерполяционная;

— построенный, по узлам, формулы (16) многочлен

п

и«(*) = П(* — га) (17)

к=1

удовлетворяет условиям,

со

/ Л-^а)гат * = °, т = М,....П — 1. (18)

0

Доказательство этой теоремы проводится точно так же, как в случае классических формул типа Гаусса [8], и мы не будем его здесь повторять.

Покажем, что многочлен (17), удовлетворяющий условиям (18), существует и определяется однозначно.

После замены переменной га = х условия (18) принимают вид

(в+1)/а-1 ехр(—х1/а) ^п(х) хт йх = 0, т = 0,1,..., п — 1.

х

Функция ю(х) = х(в+1)/а-1 ехр(—х1/а) обладает свойствами веса на полуоси (0, то), поскольку (в + 1)/а > 0, следовательно, искомый многочлен существует и единствен, а его корни, т.е. гак, к = 1, 2,...,п, попарно различны и положительны. Все коэффициенты формулы положительны. Итак, квадратурная формула типа Гаусса вида (16) существует. Опишем способ вычисления узлов и коэффициентов формулы (16). Будем искать многочлен (17) в виде

^п(г) = гп + Ь^-1 +----+ Ьп

с неизвестными коэффициентами Ьк.

Условия (18) приводят к системе линейных алгебраических уравнений

п

Е

3 = 1

Г(в + (к + п — з)а + 1) Ьз = —Г(в + (к + п)а +1), к = 0,1,... ,п — 1.

Ее решение существует и единственно, как показано выше. Далее находим корни уравнения шп(г) = 0, т. е. числа к = 1, 2,..., п. Коэффициенты формулы (16) определяем из системы уравнений

£ А (гаГ1 = Г(в + (з — 1) а +1), з = 1, 2,..., п. к=1

Заметим, что в отличие от ОКФНСТ (2) узлы и коэффициенты формулы (16) вещественны.

с

Были проведены численные эксперименты по вычислению функций (14) и (15), для чего использовался математический пакет MAPLE, позволяющий проводить вычисления с указанным в параметре Digits пакета количеством десятичных знаков. Найденные с помощью формулы (16) с пятнадцатью узлами при Digits=50 значения этих функций практически совпали с приведенными в книге [2] результатами.

Замечание 1. Описанная схема построения специальных квадратурных формул обращения преобразования Лапласа пригодна и для более общих изображений, изученных в работе [3].

Замечание 2. Параметр а может принимать любые положительные нецелые значения. Малые а характерны для изображений, описывающих медленно протекающие процессы.

Список литературы

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

[2] Работнов Ю.Н. Элементы наследственной механики твердых тел. М.: Наука, 1977. 384 с.

[3] Екельчик B.C., Рябов В.М. Об использовании одного класса наследственных ядер в линейных уравнениях вязкоупругости // Механика композитных материалов. 1981. № 3. С. 393-404.

[4] Работнов Ю.Н., Паперник Л.Х., Звонов E.H. Таблицы дробно-экспоненциальной функции отрицательных параметров и интеграла от нее. М.: Наука, 1969. 132 с.

[5] Рябов В.М. О многочленах, возникающих при численном обращении преобразования Лапласа // Методы вычислений. Вып. 12. Л.: Изд-во Ленингр. гос. ун-та, 1981. С. 46-53.

[6] Рябов В.М. Свойства квадратурных формул наивысшей степени точности, применяемых для обращения преобразования Лапласа // Журн. вычисл. математатики и мат. физики. 1989. Т. 29, № 7. С. 1083-1087.

[7] Bobylev А. V., Cercignani G. The inverse Laplace transform of some analytic functions with an application to the eternal solutions of the Boltzmann equation // Appl. Math. Lett. 2002. Vol. 15. P. 807-813.

[8] Мысовских И.П. Лекции по методам вычислений. СПб.: Изд-во С.-Петербург, гос. ун-та, 1998. 472 с.

Поступила в редакцию 15 сентября 2006 г.