Об обращении преобразования Лапласа, зависящего от некоторой степени параметра преобразования Текст научной статьи по специальности «Математика»

УДК 518:517.948

Вестник СПбГУ. Сер. 1, 2005, вып. 3

В. М. Рябов

ОБ ОБРАЩЕНИИ ПРЕОБРАЗОВАНИЯ ЛАПЛАСА, ЗАВИСЯЩЕГО ОТ НЕКОТОРОЙ СТЕПЕНИ ПАРАМЕТРА ПРЕОБРАЗОВАНИЯ*

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

/• ж

F(p) = e-px f (x) dx,

J 0

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

1 гс+гж

f(t) = — eptF(p) dp, О 0,

J с-ЪЖ

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

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

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

2. Построение квадратурных формул обращения. Предположим, что при некотором s > 0 функция уs(p) = psF(p) регулярна в полуплоскости Re (p) > 0. Запишем формулу обращения Римана—Меллина в виде

1 /■ с+гж fs-l f с+гж

f(t) = — / eP W»(P) dp = ~z г / ePp-\s(p/t) dp. (1)

J с-гж J с-гж

* Работа выполнена при финансовой поддержке РФФИ (грант №05-01-00984) и Совета по грантам Президента РФ для государственной поддержки молодых российских ученых и ведущих научных школ (грант НШ 2268.2003.1).

© В. М. Рябов, 2005

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

/•с+гж п

9т/ (2)

2пг ■) с—гж к=1

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

п 1

= ТйТ7)> з = о,1,...,п-1. (з)

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

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

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

P(s\x)= Е-1)"-^t) (n + S - 1)fc ^'

Ь = 0 v '

(a)k =

к=П

1, k = 0, a(a +1) ••• (a + k — 1), k > 1.

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

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

п

f (t) « ts-1Y] AkVs(pk/t), t> 0. (4)

k=i

По построению она точна для функций f (t) = ts-1+j ,j = 0,1,..., 2n — 1.

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

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

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

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

p—p—>0

КФНСТ (4) при t ^ 0 и t ^ж приводит к верным результатам лишь для s = 1. 3. Построение обобщенных квадратурных формул. В задачах линейной вяз-коупругости [2], описывающих напряженное состояние на основе определяющего соотношения Больцмана—Вольтерра (пространственные координаты для простоты опущены)

e(t) = i +ß jQK{t- r)a(r) dr^j , (5)

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

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

£=0 г((1 + а)(1 + к))

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

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

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

Г1(а,х)= Ь—аЭа(х), Г2(а,х)= г—а—1[ Эа (т) ¿т, х = га+1.

Jo

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

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

11

а = 1 + а.

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

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

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

/•с+гж п

9т/ (7)

¿п% ,] с—гж д,= 1

точные для функций <^(р) = р аз,з = 0,1,...,2п — 1, или для оригиналов вида tS — 1Q2n—l(ta).

Такие формулы были введены в работе [5] и подробно исследованы в статье [6]. В статье [6] доказана

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

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

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

n

"п(Х) = П (Ж " P-a)

k=1

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

/• о+гж

/ врр-3шп(р-а )р-ат ¿р = 0, т = 0,1,...,п - 1.

^ о-гж

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

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

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

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

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

Г±^) = Ф1^а ехр(±ш^)), г> 0.

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

Теорема 2. Пусть выполнены условия

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

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

F(rexp(±^)) р(±¥))|<а(г),

1 + г

где а(г) не зависит от у и a(r)exp(—Sr) G Li(R+) для любого S > 0. Тогда

1 Í

f{x) = C-l{F) = - e~xt\mF~ (t) dt. (8)

п J о

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

^ , . 1 ta sin па

lmF~ (t) = Im-

ta exp( —iпа) + 1 1 + 2ta cos па + t2a 70

так что выполнены все условия теоремы 2 и формула (8) дает

:in па Гœ ta dt

Эа (x)

n J0 1 + 2ta cos na + t2a

„ -, sin na zae-z dz . .

c - / ^-9

n J0 z2a + 2zaxa cos na + x2a

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

Положим x

9a(x)= í 3a(t) dt. (10)

Jo

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

i

Gaip)

p(pa + 1)

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

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

1 pa-1

G°(p) = - - 4-т (п)

p pa + 1

и положим

pa-1

Q°(p) = 4-т- (12)

pa + 1

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

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

„ . . ta-1 exp(—in(a — 1)) ta-1 sin па

ImF"(í) = Im- V V

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

ta exp(—¿na) + 1 í + 2ta cos na + t2a'

, ч sin na Гж _xt ta— dt

qa{x) = - / e -—

y ' n J0 1 + 2ta cos na + t2a

, sin -к a, f°° za-xe-zdz

'o

n ,l0 z2a + 2zaxa cos na + x2a

(13)

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

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

. sin na zae-z dz

Эа(х)=ха 1- / ---z-, (14)

w n J0 z a + 2zaxa cos na + x2a v 7

. . п Бт па [ж га 1е 2 ¿г . .

аа{х)=ха- / —-—. 15

п Уо г + 2гаха сов па + х2а У 7

Для их приближенного вычисления можно применить квадратурные формулы типа Гаусса [8] с весом Лагерра гае-2 для первого интеграла и с весом га- 1е-2 для второго интеграла. Однако, как отмечалось выше, с уменьшением а точность формул будет уменьшаться.

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

г 00 п

/ гве-2/(г) ¿г «V Л/(гк), ¡3 > —1, (16)

]о к=1

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

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

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

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

п

"п(г)=П(г — гка) (17)

к=1

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

/• ж

/ гв е- 2 ип(га)гат ¿г = 0, т = 0, 1,...,п — 1. (18)

о

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

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

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

/• ж

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

о

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

Шп(г) = гп + Ь1гп-1 + ••• + Ьп

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

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

п

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

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

n

Y^Ak (zak)j-1 =Г(в + (j - 1) а +1), j = 1, 2,...,n. к=1

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

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

Замечания.

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

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

Summary

V. M. Ryabov. On the Laplace transform inversion which is depending on some transform's parameter power.

Quadrature formulas for the Laplace transform inversion which are exact for inverse transform in the form of ts-1 P2n-1 (t), where s > 0, 0 < a < 1, P2n-1(t) —polynomial, n — number of nodes, are proposed.

Литература

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

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

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

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

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

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

7. Bobylev A. V., Cercignani G. 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.

8. Мысовских И. П. Лекции по методам вычислений. СПб., 1998. 472 с.

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