CC BY
17
УДК 519.61 Вестник СПбГУ. Математика. Механика. Астрономия. 2022. Т. 9 (67). Вып. 1
MSC 65F22
Метод моментов в задаче обращения преобразования Лапласа и его регуляризация*
А. В. Лебедева, В. М. Рябов
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Для цитирования: Лебедева А. В., Рябов В. М. Метод моментов в задаче обращения преобразования Лапласа и его регуляризация // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2022. Т. 9(67). Вып. 1. С. 46-52. https://doi.org/10.21638/spbu01.2022.105
Рассматриваются интегральные уравнения первого рода, относящиеся к классу некорректных задач. Сюда же относится задача обращения интегрального преобразования Лапласа, применяемого для решения широкого класса математических задач. Интегральные уравнения сводятся к плохо обусловленным системам линейных алгебраических уравнений, неизвестными в которых являются коэффициенты разложения в ряд по смещенным многочленам Лежандра некоторой функции, просто выражающейся через искомый оригинал. Эта функция находится как решение некоторой конечной проблемы моментов в гильбертовом пространстве. Для получения надежного решения системы используют методы регуляризации. Общей стратегией является использование стабилизатора Тихонова или его модификаций. Указан конкретный вид стабилизатора в методе регуляризации, ориентированный на априорно невысокую степень гладкости искомого оригинала. Приведены результаты численных экспериментов, подтверждающие эффективность предлагаемого алгоритма обращения. Ключевые слова: система линейных алгебраических уравнений, интегральные уравнения первого рода, некорректные задачи, плохо обусловленные задачи, число обусловленности, метод регуляризации.
1. Применение интегрального преобразования Лапласа приводит к более простому уравнению относительно изображения искомого оригинала. На следующем шаге возникает задача обращения, т. е. нахождения оригинала f (t) по его изображению F(p) из уравнения
/ e-pt f (t) dt = F(p). (1)
J 0
Теория преобразования Лапласа и аналитические методы его обращения содержатся в классических работах [1, 2].
Формула обращения задается интегралом Римана — Меллина:
1 f C+iTO
f(t) = — eptF(p) dp, с >7, (2)
JC-iTO
где y — абсцисса сходимости интеграла Лапласа (1).
* Статья подготовлена при поддержке гранта Санкт-Петербургского государственного университета (Мероприятие 3, Pure ID 75207094). Исследование выполнено с использованием оборудования ресурсного центра «Вычислительный центр СПбГУ» (http://cc.spbu.ru).
© Санкт-Петербургский государственный университет, 2022
К сожалению, формула (2) мало пригодна для вычислений, и потому возникает необходимость в разработке и использовании приближенных методов обращения.
Методам обращения посвящены книги [3, 4]. Выбор метода обращения определяется той априорной информацией об оригинале и его образе, которой мы располагаем. Известно, что образ ^(р), как функция комплексной переменной р, регулярна в полуплоскости Ие (р) > 7, так что мы можем при построении метода обращения использовать значения образа либо в окрестности начала координат, либо в окрестности бесконечно удаленной точки, либо на вещественной оси, либо во всей полуплоскости регулярности и т. д.
Частные случаи реализации возможных ситуаций и аналитический вид приближений подробно описаны в книге [4]. С целью указания особенностей задачи обращения остановимся на одном методе — квадратурных формулах наивысшей алгебраической степени точности, имеющих следующий вид [4]:
Функция фз(р>) предполагается регулярной в полуплоскости Ие (р) > 0. Узлы рк и коэффициенты Ак формулы (3) комплексны и таковы, что она точна для функций <ре(р) = р— ,3 = 0,1,..., 2п — 1. Узлы рк попарно различны, расположены в полуплоскости Ие (р) > 0 и являются корнями специальных ортогональных многочленов [4], обладающих теми же свойствами, что и классические ортогональные многочлены [5]. Скорость сходимости формул (3) такая же, как у классических формул Гаусса, и характеризуется величиной [4]
Существенное отличие формулы (3) от классических квадратур типа Гаусса состоит в характеристике ее устойчивости относительно погрешности в интегрируемой функции ^я(р): допущенная ошибка может возрасти в Мп = ^|Ак| раз. В работе [4] доказано, что Мп = 0(п1-я3.764п). Следовательно, вычислительный процесс по формуле (3) с ростом п становится неустойчивым.
Отметим, что задача обращения, т. е. решения интегрального уравнения первого рода (1), некорректна [6], что, в частности, проявляется в указанной неустойчивости вычислительного процесса, а также в том, что правая часть формулы (3) есть гладкая функция, а левая часть, как и функция-оригинал, даже может иметь разрывы.
В книге [4] предложен способ определения точек разрыва оригинала и величины скачка оригинала в них.
Заметим, что в названных книгах [1-4] проблема некорректности задачи обращения не рассматривалась.
Как правило, любой приближенный метод обращения приводит к системе линейных алгебраических уравнений (СЛАУ) либо относительно значений оригинала в некоторых точках, либо коэффициентов разложения оригинала в ряд по некоторой системе функций. Из-за некорректности задачи обращения получаемая СЛАУ плохо обусловлена. В случае применения квадратурных формул исследование числа обусловленности в зависимости от выбора узлов квадратуры было проведено в работе [7].
п
(3)
к=1
Общая теория решения некорректных и плохо обусловленных уравнений изложена в фундаментальных работах [6, 8], вопросы реализации общей теории применительно к конкретным прикладным задачам изложены в книге [9] и многочисленных статьях.
По сути, процедура регуляризации некорректно поставленной задачи для уравнения
Ая = и (4)
состоит из двух шагов: сначала переходим к уравнению
А* Аг = А*и,
которое всегда разрешимо, в отличие от исходного уравнения, а затем осуществляем сдвиг
(А*А + = А*и, а> 0, (5)
где Ь — некоторый симметричный положительно определенный оператор, содержащий в себе априорную информацию об искомом решении. Фактически эти переходы осуществляются путем минимизации функционала вида
Ма(А, и, г) = ||Аг — и||2 + (6)
где П — неотрицательный функционал, производная которого порождает оператор Ь.
Уравнение (5) относительно точки минимума функционала (6) стандартными приемами сводится к СЛАУ относительно каркасов приближенных решений [10].
Методы регуляризации СЛАУ и интегральных уравнений первого рода изучались нами в работах [11, 12].
2. Следуя работам [13, 14], положим в уравнении (1)
рк = a + гк, Ьк = F (рк), k = 1,...,т, (7)
сделаем замену x = ехр(— г£) и введем функцию
ад = (-'Щ . ,8,
гг
В результате из уравнения (1) получаем соотношения
I хк-1^(х) ¿х = Ьк, к = 1, .. ., т. (9)
ио
Напомним, что значения оригинала f (+0), f если они существуют, можно вы-
числить по формулам
f(+0) = lim pF(p), f(+ro) = lim pF(p).
Нахождение функции h(x), удовлетворяющей условиям (9), можно рассматривать как частный случай полной проблемы моментов [15].
Приближение к функции (8) будем разыскивать в виде разложения по смещенным многочленам Лежандра:
m
h(x) cjPj-i(2x - 1). (10)
j=i
Из формул (8) и (10) при x ^ 1 получаем равенство
/(+0) ™
j=i
Подстановка представления (10) в уравнения (9) приводит к СЛАУ вида (4) с матрицей
A = (akj j akj = J xk-íPj-í(2x - 1) dx,
правой частью u = (6i, 62,..., 6m)' и неизвестным вектором z = (ci, C2,..., cm)'.
Требовать от функции h(x) априорно какой-либо степени гладкости, накладывая некоторые условия на коэффициенты представления (10), не стоит, учитывая связь (8) функции h(x) с искомым оригиналом, необязательно гладким. В таком случае полагаем
m
= z2, L = E, j=i
и уравнение Эйлера (5) принимает вид
(A*A + aE)z = A*u, a> 0. (11)
Для решения этой системы можно применять описанные в работе [11] методы регуляризации. Разумный выбор значений входящих в (7) параметров a, r описан в работе [14].
В общем случае матрица A и правая часть и известны приближенно с погрешностью 5, и выбор параметра регуляризации a зависит от этой погрешности [6, 11]. В приводимых далее примерах матрица A и правая часть и известны практически точно, т.е. 5 — это погрешность представления чисел в компьютере, и параметр a имеет тот же порядок.
В приведенных ниже примерах вычисления проводились с помощью пакета Maple 15 с использованием оборудования Ресурсного центра «Вычислительный центр СПбГУ» (http://www.cc.spbu.ru/), где параметр Digits задает число десятичных знаков, с которым проводятся вычисления. На рис. 1, 2 представлены решения системы (11) с регуляризацией и без нее, а также значения всех использованных параметров.
Пример 1. Дано изображение
168 8640
F(p) =
(р + 2)5 (р + 3)7 ' Результаты вычислений см. на рис. 1.
Пример 2. Дано изображение
= ехр(—р) - ехр(-Зр) р2
Результаты вычислений см. на рис. 2.
Заметим, что при отсутствии регуляризации (а = 0) получаемые приближения совершенно не похожи на искомые оригиналы, в то время как регуляризация приводит ко вполне удовлетворительным результатам.
Литература
1. Лаврентьев М.А., Шабат Б. В. Методы теории функций комплексного переменного. Москва, Лань (2002).
2. Крылов В. И., Скобля Н. С. Методы приближенного преобразования Фурье и обращения преобразования Лапласа. Москва, Наука (1974).
3. Cohen A.M. Numerical methods for Laplace transform inversion. New York, Springer (2007).
4. Рябов В. М. Численное обращение преобразования Лапласа. Санкт-Петербург, Изд-во С.-Петерб. ун-та (2013).
5. Суетин П. К. Классические ортогональные многочлены. Москва, Наука (1976).
6. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. Москва, Наука (1979).
7. Gautschi W. On the condition of a matrix arising in the numerical inversion of the Laplace transform. Mathematics of computation 23 (105), 109—118 (1969).
8. Иванов В. К., Васин В. В., Танана В. П. Теория линейных некорректных задач и ее приложения. Москва, Наука (1978).
9. Кабанихин С. И. Обратные и некорректные задачи. Новосибирск, Сибирское научное изд-во (2009).
10. Даугавет И. К. Теория приближенных методов. Линейные уравнения. Санкт-Петербург, БХВ-Петербург (2006).
11. Лебедева А. В., Рябов В. М. О численном решении систем линейных алгебраических уравнений с плохо обусловленными матрицами. Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия 6 (64), вып. 4, 619—626 (2019). https://doi.org/10.21638/11701/spbu01.2019.407
12. Лебедева А. В., Рябов В.М. О регуляризации решения интегральных уравнений первого рода с помощью квадратурных формул. Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия 8 (66), вып. 4, 593—599 (2021). https://doi.org/10.21638/spbu01.2021.404
13. Papoulis A. A new method of inversion of the Laplace transform. Quarterly of applied mathematics 14 (4), 405-414 (1967).
14. Brianzi P., Frontini M. On the regularized inversion of Laplace transform. Inverse problems 7, 355-368 (1991).
15. Крейн М.Г., Нудельман А. А. Проблема моментов Маркова и экстремальные задачи. Москва, Наука (1973).
Статья поступила в редакцию 17 июля 2021 г.;
доработана 25 августа 2021 г.; рекомендована к печати 2 сентября 2021 г.
Контактная информация:
Лебедева Анастасия Владимировна — канд. физ.-мат. наук, доц.; [email protected] Рябов Виктор Михайлович — д-р физ.-мат. наук, проф.; [email protected]
Method of moments in the problem of inversion of the Laplace transform and its regularization*
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. Method of moments in the problem of inversion of the Laplace transform and its regularization. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2022, vol. 9(67), issue 1, pp. 46-52. https://doi.org/10.21638/spbu01.2022.105 (In Russian)
Integral equations of the first kind are considered, which belong to the class of ill-posed problems. This also includes the problem of inverting the integral Laplace transform, which is used to solve a wide class of mathematical problems. Integral equations are reduced to ill-conditioned systems of linear algebraic equations, in which the unknowns are the coefficients of the expansion in a series in the shifted Legendre polynomials of some function that simply expresses in terms of the sought original. This function is found as a solution to a certain finite moment problem in a Hilbert space. To obtain a reliable solution to the system, regularization methods are used. The general strategy is to use the Tikhonov stabilizer or its modifications. A specific type of stabilizer in the regularization method is indicated,
*This paper was prepared with the support by a grant from St Petersburg State University (Event 3, Pure ID 75207094). Research was carried out using computational resources provided by Resource Center "Computer Center of SPbU" (http://www.cc.spbu.ru/en).
focused on a priori low degree of smoothness of the desired original. The results of numerical
experiments are presented, confirming the effectiveness of the proposed inversion algorithm.
Keywords: system of linear algebraic equations, integral equations of the first kind, ill-posed
problems, ill-conditioned problems, condition number, regularization method.
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. Cohen A.M. Numerical methods for Laplace transform inversion. New York, Springer (2007).
4. Ryabov V. M. Numerical inversion of the Laplace transform. St Petersburg, St Petersburg Univ. Press (2013). (In Russian)
5. Suetin P. K. Classical orthogonal polynomials. Moscow, Nauka Publ. (1976). (In Russian)
6. Tikhonov A. N., Arsenin V. Ya. Metody resheniia nekorrektnykh zadach. Moscow, Nauka Publ. (1979). (In Russian) [Eng. transl.: Tikhonov A.N., Arsenin V. Ya. Solutions of Ill-Posed Problems. Winston (1977)].
7. Gautschi W. On the condition of a matrix arising in the numerical inversion of the Laplace transform. Mathematics of computation 23 (105), 109—118 (1969).
8. Ivanov V. K., Vasin V. V., Tanana V. P. Theory of linear ill-posed problems and its applications. Moscow, Nauka Publ. (1978). (In Russian)
9. Kabanikhin S. I. Inverse and ill-posed problems. Novosibirsk, Sibirskoe nauchnoe izdatel'stvo Publ. (2009). (In Russian)
10. Daugavet I. K. The theory of approximate methods. Linear equations. St Petersburg, BHV-Petersburg Publ. (2006). (In Russian)
11. Lebedeva A.V., Ryabov V. M. Numerical solution of systems of linear algebraic equations with ill-conditioned matrices. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy 6 (64), iss. 4, 619-626 (2019). https://doi.org/10.21638/11701/spbu01.2019.407 (In Russian) [Eng. transl.: Vestnik St Petersburg University, Mathematics 52 (4), 388-393 (2019). https://doi.org/10.1134/S1063454119040058].
12. Lebedeva A.V., Ryabov V. M. On the regularization of the solution of integral equations of the first kind using quadrature formulas. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy 8 (66), iss.4, 593-599 (2021). https://doi.org/10.21638/spbu01.2021.404 (In Russian) [Eng. transl.: Vestnik St Petersburg University, Mathematics 54 (4), 361-365 (2021). https://doi.org/10.1134/S1063454121040129].
13. Papoulis A. A new method of inversion of the Laplace transform. Quarterly of applied mathematics 14 (4), 405-414 (1967).
14. Brianzi P., Frontini M. On the regularized inversion of Laplace transform. Inverse problems 7, 355-368 (1991).
15. Krein M.G., Nudel'man A. A. Problema momentov Markova i ekstremal'nye zadachi. Moscow, Nauka Publ. (1973). (In Russian) [Eng. transl.: Krein M.G., Nudel'man A. A. The Markov moment problem and extremal problems. In Ser.: Translations of Mathematical Monographs, vol. 50. AMS (1977)].
Received: July 17, 2021 Revised: August 25, 2021 Accepted: September 2, 2021
A u t h o r s' i n fo r m a t i o n:
Anastasia V. Lebedeva — [email protected] Victor M. Ryabov — [email protected]
