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

УДК 519.61 Вестник СПбГУ. Математика. Механика. Астрономия. 2024. Т. 11 (69). Вып. 1 МБС 65F22

Характеристики сходимости и устойчивости некоторых методов обращения преобразования Лапласа*

А. В. Лебедева, В. М. Рябов

Санкт-Петербургский государственный университет,

Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9

Для цитирования: Лебедева А. В., Рябов В. М. Характеристики сходимости и устойчивости некоторых методов обращения преобразования Лапласа // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2024. Т. 11 (69). Вып. 1. С. 115— 130. https://doi.org/10.21638/spbu01.2024.107

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

Ключевые слова: преобразование Лапласа, обращение преобразования Лапласа, интегральные уравнения первого рода, квадратурные формулы, некорректные задачи, плохо обусловленные задачи.

Введение. Интегральное преобразование Лапласа Г(р) функции /(Ь),

представляет собой мощный инструмент для решения широкого класса прикладных задач математической физики. Одним из его главных достоинств является алгебра-изация процедур математического анализа, с помощью которой удается свести интегральные и дифференциальные уравнения к более простым. Кроме того, изображение Лапласа является регулярной функцией в некоторой полуплоскости Re p > y, что позволяет привлечь к исследованию решаемой задачи результаты теории функций комплексного переменного.

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

* Статья подготовлена при поддержке гранта Санкт-Петербургского государственного университета Мероприятие 3 (Pure ID 75207094).

© Санкт-Петербургский государственный университет, 2024

(1)

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

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

1) известны значения изображения Г(р) и его производных в некоторой фиксированной точке, отличной от бесконечности;

2) известны значения изображения Г(р) и его производных в некоторой окрестности бесконечно удаленной точки;

3) известны значения изображения Г(р) на вещественной полуоси р ^ 0;

4) известны значения изображения Г(р) во всей полуплоскости Иер >

Теория преобразования Лапласа и аналитические методы его обращения содержатся в классических работах [1, 2].

Выбор подходящих методов обращения для указанных ситуаций, их описание либо отсылка к соответствующей литературе рассмотрены в работе [3].

Задача обращения преобразования Лапласа относится к классу некорректных задач. Нас интересует характер ее неустойчивости, т. е. чувствительность искомого решения к ошибкам задания изображения Г(р), а также скорость сходимости конкретных методов обращения.

Формула обращения задается интегралом Римана — Меллина:

где 7 — абсцисса сходимости интеграла Лапласа (1). Интеграл (2) понимается в смысле главного значения и на самом деле он доставляет значение (I + 0) + ](I —

Ниже рассмотрены некоторые методы обращения и указаны их характеристики сходимости и устойчивости.

1. Обращение преобразования Лапласа с помощью разложения в ряды Лагерра. Дано уравнение

где х(Ь) — искомая вещественная функция-оригинал, а Г(р) — ее изображение по

(2)

0))/2.

(3)

Лапласу.

Рассмотрим систему функций Лагерра:

{ШУ^О, Ш = ехР(—г/2)ьк(г),

ш = -— {Ьке г) — многочлен Лагерра. к! ¿г

Эта система полна и ортонормирована в пространстве ¿2(0, то).

Пусть а — произвольное положительное число. Система функций

ы*)}£=0, = (4)

также будет полна и ортонормирована в ¿2(0, то).

Предположим, что решение х(г) уравнения (3) при некотором а > 0 представимо в виде ряда по системе функций (4):

х(г) = скт(£).

к=0

Подставляя его в уравнение (1), получаем

]ГСк (Апк )(р) = F (р). (5)

к=0

Следовательно, задача нахождения х(г) сводится к определению коэффициентов Ск из уравнения (5).

Как известно, изображение F(р) является регулярной функцией комплексной переменной р в полуплоскости Иер ^ ро при некотором вещественном ро.

Заметим, что в случае конечного р о переход к новому оригиналу /1^) = /(г) ехр(— р0t) приводит нас к случаю р0 = 0. Пусть выполнены следующие условия:

1) F(р) регулярна в полуплоскости Иер ^ 0,

2) существует конечный предел pF(р),

3) все особые точки F(р) можно заключить в замкнутый круг, целиком лежащий в полуплоскости Ие р < 0.

По изображениям многочленов Лагерра

/• то

(АЬк )(р)= ехр(—рг)Ьк(г) ¿г =

0

вычисляем изображения функций системы (2):

Их подстановка в (5) приводит к уравнению

то / /п \ к

Положим

Ък = (-1)кск^Б, к = 0,1,...

и сделаем замену переменной

г = (6)

а/2+р К '

Вестник СПбГУ. Математика. Механика. Астрономия. 2024■ Т. 11 (69). Вып. 1

Преобразование (6) конформно отображает правую полуплоскость Rep ^ 0 на круг \z\ ^ 1. Функция

ад=(£+р>(р) = "Г~—F ■ 7—— ^ (7)

V2 / p=a(1-z )/(2(i+z)) 1 + z \2 1 + z)

в силу сделанных относительно F (p) предположений будет регулярна на круге

\z\ < 1.

Итак, функция (7) имеет представление

то

G(z) = YJ buzk, (8)

k=0

причем радиус сходимости ряда (8) не менее единицы.

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

Замечание 1. Если множество особых всех точек F (p) лежит в полуплоскости Rep < 0, но неограничено, то при любом выборе параметра а радиус сходимости ряда (8) не превосходит единицы.

Замечание 2. Функции системы (4) стремятся к нулю при t ^ +ж. Естественно считать, что и оригинал f (t) при t ^ ж также стремится к нулю. Если значение f (+ж) существует и не равно нулю, то можно ввести новый оригинал f1(t) = f (t) — f (+ж) либо f2(t) = f (t) — (1 — exp(—t))f (+ж). Значение f (+ж) можно вычислить по формуле f (+ж) = limp^o pF(p).

Пусть радиус сходимости ряда (8) больше единицы. Тогда на единичном круге \z\ ^ 1 он сходится абсолютно и равномерно, и его там можно почленно дифференцировать любое число раз. Очевидно, что в этом случае выполнено неравенство Sfc=0 \ak \ < ж, которое гарантирует равномерную сходимость ряда (8), поскольку функции Лагерра (4) равномерно ограничены.

Предположим, что мы выбрали узлы интерполирования на круге \z\ ^ 1 так, что интерполяционный процесс по ним сходится равномерно к G(z). В таком случае наряду со сходимостью интерполяционных многочленов к G(z) будут сходиться и производные интерполяционных многочленов к соответствующим производным функции G(z), в частности и их значения при z = 0, т.е. коэффициенты интерполяционных многочленов будут стремиться к коэффициентам ряда (8), которые мы и разыскиваем.

Приведем необходимые для дальнейшего сведения из теории сходимости интерполяционных процессов на комплексной плоскости. Для детального изучения вопроса рекомендуем книги [5, 6].

Пусть K — компактное множество в C, дополнение Kc = C\K которого является односвязной областью расширенной комплексной плоскости C. Тогда существует конформное отображение

z = ip{w) = cw + со Н---Ь • • •, с > 0, (9)

w

внешности единичного круга \w\ > 1 на Kc. 118

Выберем узлы интерполирования гкп)(к = 0,1,...п) для всех п = 0, 1,..., принадлежащие К, и положим

"п(г)=П (* - **) . (10)

к=0

Рассмотрим числа

Мп = шах{ |шп(г)| : г € К}, п = 0,1,...

Они удовлетворяют неравенству Мп ^ еп+1.

Определение. Узлы г^^ называются равномерно распределенными (или рав-нораспределенными) на К, если

-> с. (И)

Теорема 1. Для того чтобы для каждой регулярной на компакте K функции f интерполяционный процесс равномерно сходился, т. е.

Ln(z) ^ f (z) (n z е K),

необходимо и достаточно, чтобы узлы интерполяции z^ были равномерно распределены на K.

Введем обратное отображение:

w = f(z)

и линии уровня

Cp = {z \\ф)\ = р,р> 1} .

Теорема 2. Пусть Cp = {z : \^(z) \ = р} — линии уровня, Bp = Int Cp и f регулярна на Bp; тогда

рП

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

Итак, для обеспечения равномерной сходимости интерполяционных многочленов узлы интерполирования, т.е. корни многочленов (10), следует выбирать так, чтобы выполнялось условие (11).

Известны три принципа выбора узлов на компакте, обеспечивающие равномерную сходимость: Вандермонда (Фекете), Чебышёва, Фейера.

Принцип Вандермонда. Для заданных точек е K (k = 0,1,...,n) вычисляем определитель Вандермонда:

' 1 1 ... 1

zo zi ... zn 2 2 2 z0 z1 ... zn

W (zo,zi,...,zn) =

nn o zi

n

Узлами Вандермонда называется такой набор узлов, на котором достигается максимум модуля определителя. Такие узлы попарно различны и расположены на границе дК = дКс компакта К.

Если К — единичный круг плоскости г, то узлы Вандермонда суть корни многочлена (с точностью до поворота на любой угол) шп(г) = гп+1 — 1. Очевидно, в этом случае Мп = 2, конформное отображение (9) имеет вид г = л, и условие (11) выполняется.

Принцип Чебышёва. Как и выше, выберем узлы интерполирования

)

к

zkn\k = 0,1,...n) для всех n = 0,1,..., принадлежащие K, и положим

>(z) = П (z -

(n)

■4 )

к=0

Функция |wn(z)| ограничена и непрерывна на K, и потому существует набор точек *0n), z*[n)..., *Пп) из K, на котором она принимает минимальное значение. Этот набор называют n-й системой узлов Чебышёва.

Если K — единичный круг плоскости z, то узлы Чебышёва суть корни многочлена шп(г) = zn+1. Следовательно, Mn = 1, конформное отображение (9) имеет вид z = w, и условие (11) выполняется. Очевидно, в этом случае речь идет об интерполировании Эрмита с одним кратным узлом z = 0.

Если K = [-1,1], то многочленами вида (10), наименее уклоняющимися от нуля на K, являются многочлены Чебышёва первого рода un(z) = Tn+i(z)/2n с корнями zkn) = cos((2k + 1)n/(2n + 2)).

Принцип Фейера. Пусть компакт K таков, что функция ф, отображающая множество {w : |w| > 1} на Kc, непрерывно продолжается на множество {w : |w| > 1}. Пусть набор узлов {w^^Yn^ равнораспределен на |w| = 1. Тогда набор точек

z(n) = Ф^) (к = 0, 1,...,n)

называется n-й системой узлов Фейера.

Пусть K = [-1,1]. В этом случае отображение (9) задается функцией Жуковского z = (w + 1/w)/2, отображающей окружность |w| = 1 на K = [-1,1].

Пусть wkn) = exp((4k + 1)ni/(2n + 2)), тогда z(n) = cos((4k + 1)n/(2n + 2)). Их удобнее перенумеровать, положив

(n) 2k + 1 , n 1 n 1

z, = cos-7Г, к = 0, 1,. .., n, n = 0,1, .. .

k 2n + 2 ' ' ' '

Следовательно, un(z) = Tn+i(z)/2n, Mn = 2-n, и условие (11) выполнено.

Замечание 3. Далее в качестве компакта K плоскости z рассматривается единичный круг или некоторое его подмножество, на которых порожденная изображением функция (7) регулярна. Скорость сходимости равномерно сходящихся интерполяционных процессов зависит от «размера» наибольшей области регулярности функции G(z) и определяется теоремой 2 в случае максимальной сходимости. В случае K = [-1,1] такой областью будет эллипс с фокусами в точках ±1, для которого сумма полуосей р максимальна.

Подробнее на этих вопросах мы не будем останавливаться и отсылаем читателей к работам [5, 6]. Там же можно найти доказательства теорем 1, 2. Методы фактического построения интерполяционных многочленов подробно рассмотрены в [4].

Перейдем к изучению влияния погрешности задания Г (р) на приближенное решение в виде отрезка ряда

г-1

Х" И = Е (г). (12)

к=0

Теорема 3 ([4]). Пусть в узлах Фейера гк = сов((2к + 1)п/(2п)), к = 0,1,..., п— 1, значения функции С(гк) известны с погрешностью ек, причем \ек\ ^ е. Тогда возникающая при этом ошибка в правой части (12) по модулю не превосходит величины е ■ ю0-384" для любого г ^ 0.

Узлами Вандермонда для круга \г\ ^ 1 являются корни из единицы, т.е. точки гк = тк, т = ещ)(—2п1/п), к = 0,1,...,п — 1. На плоскости р им соответствуют точки мнимой оси

а кп Рк = г-tg—.

2 п

В случае четного п среди них находится и бесконечно удаленная точка.

Теорема 4 ([4]). Пусть в узлах Вандермонда гк = тк значения функции С(гк) известны с погрешностью ек, причем \ек\ ^ е. Тогда возникающая при этом ошибка в правой части (12) по модулю не превосходит величины пе для любого г ^ 0.

Замечание 4. Погрешность вычисления функции О(г) складывается из ошибок в функциях Г (р) и рГ(р). Из представления

1п х

Р

и вытекающего из него неравенства

рр(р) = ^ / ^ р>0

!тъп ^ рГ (р) ^ Iтах

становится понятным смысл теорем 3, 4.

2. Обращение преобразования Лапласа с помощью квадратурных формул. Пусть Г(р) регулярна в правой полуплоскости Иер > 0 наряду с функцией ^в(р) = р'эГ(р) при некотором в > 0 и пусть существует значение Запишем формулу обращения (2) в виде

= с> 0. (13)

2п% с—г<х>

Построим квадратурную формулу вида

п

I(г) «Лк^в(рк/г). (14)

к=1

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

1) формула (14) интерполяционная, т. е. она точна для функций (р) = р—3, 3 = 0,— 1,

2) выполняются условия «ортогональности»

,-с+гж / 1 \

/ ер*р-*ш (-\р-к (1р = 0, к = 0,1,...,п- 1,

о С — 1(Х> \р/

Р

где

(15)

\р/ ^ Vр Рк)

Доказывается эта теорема так же, как и для классических квадратурных формул типа Гаусса, и оно приведено в книге [2]. Итак, существование квадратурной формулы наивысшей степени точности (КФНСТ) обращения преобразования Лапласа сводится к существованию многочлена (16), удовлетворяющего условию (15).

Оказывается, что многочлен (16) с точностью до постоянного сомножителя совпадает с многочленом

корни которого рк попарно различны и лежат в полуплоскости Ие г > 0. В записи многочлена РЩ использован символ Похгаммера:

, , ¡1, к = 0; (а)к = <

V 'к |^а(а +1)---(а + к — 1), к > 1.

Коэффициенты этой квадратурной формулы можно найти как решение системы уравнений

п1

У2АкРкт = 7Гг—:-ч, то = 0,1,.. ,,п - 1,

к Г(в + т)

эквивалентной условию интерполяционности квадратурной формулы. Для погрешности еп(Ь) КФНСТ, записанной в виде

I (V= ?—1

была получена оценка [4]

Ак ф8(рк/г) + £п(ь)

. я

,к=1

2 - 2т

1/2

+ т) к \р/ ' (17)

т=2п\ ' к=1 ' К> ' 1 У '

0 <г < р< г■

Здесь Мря — величина, зависящая от изображения и не зависящая от квадратурной формулы.

Устойчивость КФНСТ по отношению к ошибкам в функции определяется величиной

n

Mn = \Ak\ = O(n1-s3.764n),

k=i

а скорость сходимости КФНСТ — величиной [4]

Bon = Oln

1-s (3.764\

2„ = (J | П | -j

n2

В случае натурального в узлы КФНСТ совпадают с полюсами аппроксимаций Паде экспоненциальной функции ехр(р) [4].

К сожалению, КФНСТ плохо приспособлены для обращения изображений, соответствующих медленно протекающим длительным процессам. В таком случае целесообразно вместо КФНСТ построить и использовать обобщенные квадратурные формулы наивысшей степени точности (ОКФНСТ), точные для функций р) = = 0,2п — 1, или для оригиналов вида ^^(Згп-!^") (<32п-1 — произвольный многочлен), где а — любое положительное число (наибольший интерес представляет случай а € (0,1]).

Такие формулы были описаны в работе [7] и исследованы в статье [8]. Результаты их применения к решению практических задач содержатся в [9].

3. Дельта-методы обращения преобразования Лапласа и метод Видде-

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

п т

/(*) ^ЕАу(*)), (18)

к=0 з = \

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

/• ж

^ (р)= ехр(-рг)/(г) ¿г,

0

перепишем формулу (18) иначе:

г n m

f (t) ^ / E Y,(-1)mAkjeM-xPkj(t))

j 0 _ï„_П _1

l-fc = Q j=1

0

Положим

nm

m

ikj

f (x) dx. (19)

Snm(x,t) = YJY,(~1)mAkj (t)xk exp(-xpkj (t)). (20)

к=0 з = 1

Приближенное равенство (19) означает, что параметры метода обращения т, п, Акз, Ркз (г) желательно подобрать так, чтобы функция (20) была близка к дельта-функции в точке х = г. Такие ядра назовем дельтообразными.

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

К дельта-методам относятся операторы Виддера

pn+1 t ^ , wn(f, t) = p=n/t ,n = 1,2,...,

n! 1

и операторы Поста

Sn(f,t) = (-1 T^-rF{n\p) |p=(n+i)/t , П = 0,1,...

Нас интересует не только факт сходимости, но и скорость сходимости, т. е. поведение при n ^œ величин

(f,t) = Wn(f,t) - f (t), £Sn (f,t)= Sn(f,t) - f (t). (21)

Лемма 1 ([4]). Пусть f (x) непрерывно дифференцируема на полуоси x ^ 0 и \f '(x)\ < Mi. Тогда

I ff„(/,i)-/(i)|aitl- + 4 о.

V n n2

Лемма 2 ([4]). Пусть f e C2(0, œ) и \f "(x)\ < M2. Тогда

t2

Перейдем к более детальному изучению погрешностей формул (21). Как правило, выводы будем делать лишь для одной из них, ибо для второй все делается аналогично, и к тому же получающиеся результаты качественно одинаковы. Пусть f € C(0, ж), тогда погрешность (f,t) представима в виде

w,* ^ 1 ^ (1 2 A t2f"(t) ( 5 6 \ t4f"'(t)

t) = -tf'{t) + - + -П -Ш + "2 + "а -^Г1 + ■■■

п ' п \п п2) 2! \п2 п3/ 3!

В частности, отсюда при больших п получаем главный член погрешности

п2

Как видим, метод Виддера быстро насыщаем, а именно наличие производных оригинала выше второго порядка не улучшает сходимости.

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

Выберем к различных натуральных чисел и\ < < ■■■ < и положим ¿3 = пз/п1, 3 = 1, 2,...,к (впрочем, в знаменателе можно написать любое натуральное число), после чего образуем числа

= з = 1>2>->к>

г=3

и составим линейную комбинацию

к 3=1

Числа Скз выбираются из условий

= I

0, т = 1, 2,. ..,к — 1.

Ej-m т =

ckj d- m

3 = 1

Пусть a, b, ai,bi — произвольные фиксированные числа такие, что 0 < a < ai < bi < b < ж.

Теорема 6 ([4]). Пусть f (x) — непрерывный ограниченный на полуоси x > 0 оригинал, имеющий на отрезке [a, Ъ] непрерывную производную порядка 2k. Тогда при n ^ж равномерно относительно t G [ai, bi] справедливо равенство

2k

nk [Wn(k, f, t) - f (t)] = £ Mmtmf (m)(t) + o(1),

m=k

где Mm — не зависящие от t и f (t) постоянные. В частном случае k = 1 получаем

lim nk[Wn(f,t)-f(t)] = l(t2f'(t))',

п^ж 2

а при k = 2

2

lim пк [Wn(2J,t)-f(t)} = - (t3/»(t))' + - (^/'"(^'Е^Л72-

3=i

Замечание 5. Значения величин Wn(f,t) следует вычислять с возможно большей точностью, ибо в них содержится информация о характере зависимости ошибки (f,t) от 1/n, которая используется существенным образом на шаге ускорения сходимости. Отметим, что Wn(f,t) заведомо существуют при достаточно больших n (как только точка n/t попадает в область регулярности изображения).

Замечание 6. Выше поведение операторов Виддера исследовалось в предположении достаточной гладкости оригинала. Отметим следующий известный результат: если функция f (t) имеет ограниченную вариацию на отрезке [0, T] при любом T > 0 и интеграл Лапласа (1) сходится абсолютно при некотором p, то

lim L2Lxn+iF(n){Xn) = 1 (/(t + 0) + f(t - 0)),

п^ж n! 2

где хп = (п + вп)/г, 0 ^ 9п ^ 1. (см. [10], с. 289). Очевидно, в левой части под знаком предела стоят операторы Виддера Бп при вп = 0,1 соответственно.

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

Пусть г — некоторое число, 0 < г < 1, и т — натуральное число, т > п. Введем функции

ет (х)

! — 6Хр I % X

т

2п

и рассмотрим величину

Wnm(f,t)

1 т ( г \

тп 3=1 х 7

(22)

Теорема 7 ([4]). Если исходный оригинал ограничен на полуоси х ^ 0 и его изображение регулярно в полуплоскости Ие р > 0, то справедливо неравенство

т

1- гт

при всех г > 0, где М = пжх \f (г)\.

Итак, значение Wn(f,t) можно с любой точностью заменить величиной №пт(^,г), если взять достаточно большое значение т, либо достаточно малое г. Расчетную формулу (22) можно записать иначе: с вычислительной точки зрения правую часть формулы (22) удобнее выразить через значения функции у(р) = рГ(р) (как правило, именно она находится из уравнения в изображениях, и к тому же она, в отличие от Г(р), не стремится к нулю с ростом \р\):

1

И7пт — У. [гет{])\

т < ^

1

3 = 1

1 - гет(])

п

^(1 - гет

(а))

(23)

Теорема 8 ([4]). Если значения функции у(р) = рГ(р) вычисляются с погрешностью, по модулю не превосходящей числа е, и оригинал ограничен, то справед-

ливы неравенства

где

^птШ) - Wn(f,t)\ < А,

д

Мгт

1 - г"

+

(1 - г) г

(24)

е

п

В общем случае, когда некоторые из входящих в формулу (24) величин М, г, т, п, £ известны, остальные следует выбирать из условия минимума величины Д. Практически алгоритм выбора параметров таков: по известным характеристикам: п — номеру искомого '№п(/,1) и £ — точности задания (вычисления) изображения — выбираем параметр г так, чтобы второе слагаемое в сумме (24) имело желаемый порядок малости, после чего находим наименьшее значение параметра т, при котором первое слагаемое в (28) примерно такое же, как и второе. Дальнейшее увеличение т ведет к росту объема вычислений для нахождения '№п(/,1), не увеличивая его точности.

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

Устойчивость вычислений по формуле (23) определяется величиной ^3=1 \ск] \. В работе [4] показано, что при использовании комбинации приближений по Виддеру с номерами от П1 до пк возможен такой выбор к штук из них, при котором величина 5^к=1 \ск^ \ будет минимальна и равна

к

Т.\ск! \ = (-1)к-1Тк-1(0),

3 = 1

где Тк-1 (х) — смещенные на [1/пк, 1/п1] многочлены Чебышёва первого рода.

5. Деформирование линии интегрирования в формуле обращения.

В этом разделе исследуется метод вычисления интеграла Римана — Меллина:

1 о+гж

/(¿) = — / (25)

2пг ■/ о-гж

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

Если все особые точки изображения расположены в левой полуплоскости, то в формуле обращения (25) можно положить с = 0. В противном случае вместо Е(г) будем работать с функцией Е(г + а), а > 0 (соответствующей функции-оригиналу ](1)ехр(—а;1)), с особыми точками в левой полуплоскости.

В качестве контура интегрирования в формуле (25) возьмем кривую, состоящую из прямолинейных участков ¿1,Ь2,Ьз. На участке Ь параметр г изменяется по закону г = хМ, х € [-1,1], Ь > 0, на участке ¿1 имеем г = х — Ьг, х € (-то, 0], а на участке Ьз г = —х + Ьг, х € [0, то). Формула обращения записывается в виде

т = ¿//'ад* ^ ¿? (I \L\fJ «■*<«>*■

После замен переменной 2 получим

[ вг*Е(г) 3г =в-т[ в-х*Е(—х — Ьг) 3х, ■ 'Ьг ио

/ ех1Е(г) 3г =Ьг / вШхЕ(гЬх) 3х, (26)

пЖ

I е7ЛЕ(г) 3г = — вш в-хЛЕ(—х + Ьг) 3х.

и Ь3 ио

Ьо •>-1

Эффективные методы вычисления интегралов (26) предложены в [11]. Для вычисления интегралов от сильно осциллирующих функций вида

J ¿шхд(х) 3х

следует применять предложенные в работе [12] специальные квадратурные формулы.

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

/' /'Ж

/ в2 ^(г) 3г = в-ы е-х1¥(-х - Ы) 3х JL1 ■! 0

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

Литература

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

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

3. Порошина Н. И., Рябов В. М. О методах обращения преобразования Лапласа. Вестник С.-Петерб. ун-та. Сер. 1. Математика 3, 55—64 (2011).

4. Рябов В.М. Численное обращение преобразования Лапласа. Санкт-Петербург, Изд-во С.-Петерб. ун-та (2013).

5. Смирнов В. И., Лебедев Н. А. Конструктивная теория функций комплексного переменного. Москва, Наука (1964).

6. Гайер Д. Лекции по аппроксимации в комплексной плоскости, пер. с англ. Москва, Мир (1986).

7. Рябов В.М. О многочленах, возникающих при численном обращении преобразования Лапласа. Методы вычислений 12, 46—53 (1981).

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

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

10. Widder D. V. The Laplace transform. Princeton (1946).

11. Lebedeva A. V., Ryabov V. M. On Integration Contour Deformation in a Laplace Transform inversion Formula. Computational Mathematics and Mathematical Physics 55 (7), 1103-1109 (2015)

12. Бахвалов Н. С., Васильева Л. Г. Вычисление интегралов от осциллирующих функций при помощи интерполяции по узлам формул Гаусса. Журнал вычислительной математики и математической физики 8 (1), 175-181 (1968).

Статья поступила в редакцию 3 февраля 2023 г.;

доработана 2 мая 2023 г.; рекомендована к печати 31 августа 2023 г.

Контактная информация:

Лебедева Анастасия Владимировна — канд. физ.-мат. наук, доц.; a.v.lebedeva@spbu.ru Рябов Виктор Михайлович — д-р физ.-мат. наук, проф.; v.ryabov@spbu.ru

Characteristics of convergence and stability of some methods for inverting the Laplace transform*

A. V. Lebedeva, V. M. Ryabov

St. Petersburg State University, 7—9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation

For citation: Lebedeva A. V., Ryabov V. M. Characteristics of convergence and stability of some methods for inverting the Laplace transform. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2024, vol. 11 (69), issue 1, pp. 115-130. https://doi.org/10.21638/spbu01.2024.107 (In Russian)

The problem of inversion of the integral Laplace transform, which belongs to the class of ill-posed problems, is considered. Integral equations are reduced to ill-conditioned systems of linear algebraic equations, in which the unknowns are either the coefficients of the series expansion in terms of special functions, or the approximate values of the desired original at a number of points. Various inversion methods are considered and their characteristics of accuracy and stability are indicated, which must be known when choosing an inversion method for solving applied problems. Quadrature inversion formulas are constructed, which are adapted for the inversion of long-term and slow processes of linear viscoelasticity. A method of deformation of the integration contour in the Riemann-Mellin inversion formula is proposed, which leads the problem to the calculation of certain integrals and allows obtaining error estimates.

Keywords: Laplace transform, Laplace transform inversion, integral equations of the first kind, quadrature formulas, ill-posed problems, ill-conditioned problems.

References

1. Lavrent'ev M. A., Shabat B. V. Methods of the theory of functions of a complex variable. Moscow Lan' Publ. (2002). (In Russian)

2. Krylov V. I., Skoblya N. S. Methods of the approximate Fourier transform and the inversion of the Laplace transform. Moscow, Nauka Publ. (1974). (In Russian)

3. Poroshina N. I., Ryabov V. M. On methods for inverting the Laplace transform. Vestnik St Petersburg University. Ser.1. Mathematics 3, 55—64 (2011). (In Russian)

4. Ryabov V. M. Numerical inversion of the Laplace transform. St. Petersburg, St. Petersburg University Press (2013). (In Russian)

5. Smirnov V. I., Lebedev N.A. Constructive theory of functions of a complex variable. Moscow, Nauka Publ (1964). (In Russian)

6. Gaier D. Vorlesungen uber Approximation im Komplexen. Birkhauser Verlag, Basel (1980) [Rus. ed.: Gaier D. Lektsii po approksimatsii v kompleksnoi ploskosti. Moscow, Mir Publ. (1986)].

7. Ryabov V. M. On polynomials arising from the numerical inversion of the Laplace transform. Numerical methods 12, 46—53 (1981). (In Russian)

8. Ryabov V. M. Properties of Quadrature Formulas of the Highest Degree of Accuracy Used to Invert the Laplace Transform. Computational Mathematics and Mathematical Physics 29 (7), 1083—1087 (1989). (In Russian)

9. Ekelchik V. S., Ryabov V. M. On the use of one class of hereditary kernels in linear equations of viscoelasticity. Mechanics of Composite Materials 3, 393—404 (1981). (In Russian)

10. Widder D. V. The Laplace transform. Princeton (1946).

11. Lebedeva A. V., Ryabov V. M. On Integration Contour Deformation in a Laplace Transform inversion Formula. Computational Mathematics and Mathematical Physics 55 (7), 1103—1109 (2015).

*This paper was prepared with the support by a grant from St. Petersburg State University Event 3 (Pure ID 75207094).

12. Bakhvalov N.S., Vasil'eva L. G. Calculation of integrals of oscillating functions using interpolation over the nodes of Gauss formulas. Computational Mathematics and Mathematical Physics 8 (1), 175-181 (1968). (In Russian)

Received: February 3, 2023 Revised: May 2, 2023 Accepted: August 31, 2023

Authors' information:

Anastasia V. Lebedeva — a.v.lebedeva@spbu.ru Victor M. Ryabov — v.ryabov@spbu.ru