CC BY
5
УДК 519.61 Вестник СПбГУ. Математика. Механика. Астрономия. 2022. Т. 9 (67). Вып. 4
MSC 65F22
Регуляризация процедуры обращения преобразования Лапласа с помощью квадратурных формул*
А. В. Лебедева, В. М. Рябов
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Для цитирования: Лебедева А. В., Рябов В. М. Регуляризация процедуры обращения преобразования Лапласа с помощью квадратурных формул // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2022. Т. 9 (67). Вып. 4. С. 636-643. https://doi.org/10.21638/spbu01.2022.406
Рассматривается задача обращения интегрального преобразования Лапласа, относящаяся к классу некорректных задач. Интегральные уравнения сводятся к плохо обусловленным системам линейных алгебраических уравнений (СЛАУ), неизвестными в которых являются либо коэффициенты разложения в ряд по смещенным многочленам Лежандра, либо приближенные значения искомого оригинала в ряде точек. Первый шаг сведения к СЛАУ состоит в применении квадратурных формул, доставляющих минимальные значения числа обусловленности СЛАУ. Для получения надежного решения системы используют методы регуляризации. Общей стратегией является использование стабилизатора Тихонова или его модификаций. Приведен вариант метода регуляризации систем с матрицами осцилляционного типа, существенно уменьшающий обусловленность задачи по сравнению с классической схемой Тихонова. Приведен способ фактического построения специальных квадратур, приводящих к задачам с ос-цилляционными матрицами.
Ключевые слова: система линейных алгебраических уравнений, интегральные уравнения первого рода, некорректные задачи, плохо обусловленные задачи, число обусловленности, осцилляционные матрицы, метод регуляризации.
1. Применение интегрального преобразования Лапласа приводит к более простому уравнению относительно изображения искомого оригинала. После нахождения изображения возникает задача обращения, т.е. нахождения оригинала ](£) по его изображению Е(р) из уравнения
Теория преобразования Лапласа и аналитические методы его обращения содержатся в классических работах [1, 2].
Формула обращения задается интегралом Римана — Меллина:
* Статья подготовлена при поддержке гранта Санкт-Петербургского государственного университета «Мероприятие 3» (Pure ID 75207094).
© Санкт-Петербургский государственный университет, 2022
(1)
(2)
где 7 — абсцисса сходимости интеграла Лапласа (1). Интеграл (2) понимается в смысле главного значения. Если в точке £ оригинал ](£) имеет разрыв первого рода, то интеграл (2) доставляет значение (/(£ + 0) + /(£ — 0))/2. Способы нахождения точек разрыва и величины скачка оригинала ](£ + 0) — f (£ — 0) в них описаны в книге [3].
К сожалению, формула (2) мало пригодна для вычислений, и потому возникает необходимость в разработке и использовании приближенных методов обращения.
Методам обращения посвящены книги [3, 4]. Выбор метода обращения определяется той априорной информацией об оригинале и его образе, которой мы располагаем. Известно, что образ Е(р) как функция комплексной переменной р регулярна в полуплоскости Ив (р) > так что мы можем при построении метода обращения использовать значения образа либо в окрестности начала координат, либо в окрестности бесконечно удаленной точки, либо на вещественной оси, либо во всей полуплоскости регулярности и т. д.
2. Частные случаи реализации возможных ситуаций и аналитический вид приближений подробно описаны в книге [3]. С целью указания особенностей задачи обращения остановимся на одном методе — квадратурных формулах, в том числе наивысшей алгебраической степени точности (КФНСТ), имеющих следующий вид [3]:
Здесь в — скорость убывания изображения и функция (р) предполагается регулярной в полуплоскости Ив (р) > 0. Эти формулы получаются в результате построения квадратурной формулы для вычисления интеграла в формуле Римана — Меллина, переписанной в виде
Узлы рк и коэффициенты Лк формулы (3) комплексны и таковы, что она точна для функций ^я(р) = р— ,3 =0,1,..., 2п — 1. Узлы рк попарно различны, расположены в полуплоскости Ив (р) > 0 и являются корнями специальных ортогональных многочленов [3], обладающих теми же свойствами, что и классические ортогональные многочлены [5]. Скорость сходимости формул (3) такая же, как у классических формул Гаусса, и характеризуется величиной [3]
Существенное отличие формулы (3) от классических квадратур типа Гаусса состоит в характеристике ее устойчивости относительно погрешности в интегрируемой функции ^я(р) : допущенная ошибка может возрасти в Мп = ^п=1 \Лк \ раз. В работе [3] доказано, что Мп = 0(п1-я3.764п). Следовательно, вычислительный процесс по формуле (3) с ростом п становится неустойчивым.
Задача обращения преобразования Лапласа как решение интегрального уравнения первого рода (1) некорректна [6], что, в частности, проявляется в указанной неустойчивости вычислительного процесса, а также в том, что правая часть в формуле (3) есть гладкая функция, а левая часть даже может иметь разрывы.
п
(3)
к=1
Заметим, что в названных книгах [1-4] проблема некорректности задачи обращения не рассматривалась.
Как правило, любой приближенный метод обращения приводит к системе линейных алгебраических уравнений (СЛАУ) либо относительно значений оригинала в некоторых точках (т.е. каркаса приближенного решения), либо коэффициентов разложения оригинала в ряд по некоторой системе функций. Из-за некорректности задачи обращения получаемая СЛАУ плохо обусловлена.
Общая теория решения некорректных и плохо обусловленных уравнений изложена в фундаментальных работах [6, 7], вопросы реализации общей теории применительно к конкретным прикладным задачам изложены в книге [8] и многочисленных статьях.
3. Пусть СЛАУ имеет вид
Az = u, (4)
где A — вещественная квадратная матрица; u — заданный вещественный вектор; z — искомый вектор размерности п. Матрица и правая часть в уравнении (4) могут быть заданы приближениями к ним Ag и ug такими, что
||ДА|| = IIA - AgII < S, ||Ди|| = ||u - ugII < S.
Эти погрешности повлекут ошибку Дz определения вектора z. В качестве векторной нормы используем евклидову норму, а в качестве матричной — подчиненную ей норму _
Введем относительные погрешности:
_ ЦДА|| Л -UM л -ИМ
М\ ' INI ' INI •
Заметим, что эти погрешности присутствуют практически всегда, поскольку вычисления, как правило, проводятся с конечным числом значащих цифр.
Если матрица A не вырождена, то существует единственное решение системы (4) и можно ввести в рассмотрение число обусловленности ¡(A) = cond(A) = ||A|| •
IIA-11|.
Важность этой характеристики задачи видна из следующего результата.
Теорема 1 (см. [9, с. 111]). Для произвольных Su и достаточно малых Sa справедлива оценка погрешности
Следовательно, независимо от метода решения СЛАУ погрешность результата может оказаться примерно в ¡ (A) раз больше погрешности исходной информации.
4. Далее описаны некоторые алгоритмы, которые формально укладываются в общую схему регуляризации А. Н. Тихонова [6], но имеют меньшую обусловленность.
В методе Тихонова вводится функционал
Ma(A,z,u) = II Az - uI + aIIzII2, а > 0,
и ставится задача нахождения элемента га, минимизирующего этот функционал. Такой элемент однозначно определяется как решение уравнения Эйлера:
(Л*Л + аЕ )га = Л* и.
Теорема 2 (см. [6, с. 119]). Пусть г° есть нормальное решение системы Лг = и, и вместо вектора и мы имеем вектор и$ такой, что \\и — и$ \\ ^ 5. Пусть далее /31(5) и /32(5) — какие-либо непрерывные и положительные на [0, 52] функции,
монотонно стремящиеся к нулю при 5 ^ 0 и такие, что
5
</32(<5), /?2 (0) = 0.
ßi(5)
Тогда для любой положительной на [0,52] функции а = а(5), удовлетворяющей условию
с
(5)
векторы za(s), являющиеся решением уравнения Эйлера
(A*A + a(5)E)za(s) = A*u, сходятся к нормальному решению системы Az = u при 5 ^ 0, т. е.
lim - z°|| =0.
Замечание. Если система (4) имеет единственное решение, то неравенство (5) можно заменить на такое:
щг^м^т, (ß)
откуда, в частности, при ß\{5) = ß2(S) = S из (6) находим a(S) = S. Аналогично из (5) при ß\ (6) = /?2(S) = \/5 находим a(S) = \/5.
5. Сделаем в уравнении (1) замену переменной x = exp(-t) и введем функцию
h(x) = f (-ln(x)), (7)
и в результате из уравнения (1) получаем уравнение
i xp-1 h(x) dx = F(p). (8)
0
Напомним, что значения оригинала f (+0), f (+то), если они существуют, можно вычислить по формулам
f(+0) = lim pF(p), f(+то) = lim pF(p). Выберем некоторую квадратурную формулу вида
1 п
/ g(x) dx Akg(xk) (9)
Jo k=i
Вестник СПбГУ. Математика. Механика. Астрономия. 2022. Т. 9(67). Вып. 4
и применим ее для приближенного вычисления интеграла в формуле (8):
J2Ak*Pk-1 h(xk ) = F (p). (10)
k=i
Обозначим zk = Akh(xk) и положим в уравнении (9) p = pj, j = 1, 2,...,n, в результате чего придем к СЛАУ вида CZ = U, где
C = Х-Л n , Z =(zi,...,zn)', U =(F (pi),..., F(pn))'.
V / k,j=i
Исследование числа обусловленности матрицы C в зависимости от выбора узлов квадратуры было проведено в работе [10].
Если 0 < xi < x2 < ■■■ < xn, pi <p2 < ■■■ < pn, то C есть обобщенная матрица Вандермонда осцилляционного типа [11], обладающая следующими свойствами.
Теорема 3 (см. [11, с. 400]). 1. Осцилляционная матрица всегда имеет n различных положительных собственных чисел Ai > A2 > ■■■ > An > 0.2. У первого собственного вектора матрицы, отвечающего наибольшему собственному числу, все координаты отличны от нуля и одного знака. У собственного вектора, отвечающего второму по величине собственному числу, в ряду координат имеется точно одна перемена знака, и вообще в ряду координат k-го собственного вектора, соответствующего k-му собственному числу, имеется точно k — 1 перемена знака (k =1, 2,...,n).
В случае pj = j матрица C превращается в матрицу Вандермонда. Как отмечено в работе [10], плохая обусловленность системы вполне может потребовать высокоточных расчетов, особенно при достаточно больших n. В той же работе изучался вопрос о том, можно ли устранить или смягчить пагубные последствия плохой обусловленности, используя различные варианты выбора n и квадратурной формулы (9). Показано, что при выборе в качестве узлов квадратурной формулы (9) корней многочленов Якоби Pna'e)( х) число cond(C) ведет себя как величина 0((3+л/8)п), а в случае равноотстоящих узлов как O(8n). Следовательно, необходимо использование методов регуляризации при решении подобных СЛАУ и интегральных уравнений первого рода (см. работы [12, 13]).
Отметим, что классическая квадратурная формула Гаусса для отрезка [—1, 1] при любом n легко строится следующим образом: берем начальные приближения к узлам
0 (k — 1/4 \
хк = - cos I + тгj , к = 1,2,... ,п,
уточняем их стандартным методом Ньютона, а затем вычисляем коэффициенты
л к-12 п
Пусть A — осцилляционная матрица, а V есть матрица ее собственных векторов, Л = [Ai,..., An] — диагональная матрица собственных чисел матрицы A. Тогда AV = VA или А = VAV-1. Пусть Ai = [л/Ai, • • •, VK}- Положим В = VkiV~l. Понятно, что B2 = VЛ1V-i = Vti{V-i = A.
Следовательно, B есть корень из A, и собственные векторы этих матриц совпадают.
Умножая исходное уравнение Az = u слева на матрицу B-1, получим
Bz = B-1u.
Запишем функционал Тихонова для этого уравнения:
Ma(B, z, B-1u) = \\Bz - B-1u\\2 + а^Ц2, a > 0, уравнение Эйлера для него примет вид
(B*B + aE )za = B*(B-1u). (11)
Число обусловленности задачи (11) существенно меньше, чем у стандартной формы метода регуляризации. Отметим, что матрица B*B-1 в правой части уравнения (11) невелика (в случае симметрии A она равна единичной матрице) и не сильно увеличивает погрешность вектора u. Методы вычисления необходимого нам квадратного корня матрицы описаны в работе [14].
Литература
1. Лаврентьев М.А., Шабат Б. В. Методы теории функций комплексного переменного. Москва, Лань (2002).
2. Крылов В. И., Скобля Н.С. Методы приближенного преобразования Фурье и обращения преобразования Лапласа. Москва, Наука (1974).
3. Рябов В.М. Численное обращение преобразования Лапласа. Санкт-Петербург, Изд-во С.-Петерб. ун-та (2013).
4. Cohen A.M. Numerical methods for Laplace transform inversion. New York, Springer (2007).
5. Суетин П. К. Классические ортогональные многочлены.. Москва, Наука (1976).
6. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. Москва, Наука (1979).
7. Иванов В. К., Васин В. В., Танана В. П. Теория линейных некорректных задач и ее приложения. Москва, Наука (1978).
8. Кабанихин С. И. Обратные и некорректные задачи. Новосибирск, Сибирское научное изд-во (2009).
9. Воеводин В. В., Кузнецов Ю.А. Матрицы и вычисления. Москва, Наука (1984).
10. 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).
11. Гантмахер Ф.Р. Теория матриц. Москва, Наука (1967).
12. Лебедева А. В., Рябов В.М. О численном решении систем линейных алгебраических уравнений с плохо обусловленными матрицами. Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия 6 (64), вып. 4, 619—626 (2019). https://doi.org/10.21638/11701/spbu01.2019.407
13. Лебедева А. В., Рябов В. М. О регуляризации решения интегральных уравнений первого рода с помощью квадратурных формул. Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия 8 (66), вып.4, 593—599 (2021). https://doi.org/10.21638 /spbu01.2021.404
14. Higham N. J. Functions of matrices: Theory and computation. Philadelphia, Society for Industrial and Applied Mathematics (2008).
Статья поступила в редакцию 5 мая 2022 г.;
доработана 3 июня 2022 г.; рекомендована к печати 9 июня 2022 г.
Контактная информация:
Лебедева Анастасия Владимировна — канд. физ.-мат. наук, доц.; a.v.lebedeva@spbu.ru
Рябов Виктор Михайлович — д-р физ.-мат. наук, проф.; v.ryabov@spbu.ru
Regularization of the procedure for inverting the Laplace transform using quadrature formulas*
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. Regularization of the procedure for inverting the Laplace transform using quadrature formulas. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2022, vol. 9 (67), issue 4, pp. 636-643. https://doi.org/10.21638/spbu01.2022.406 (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 (SLAE), in which the unknowns are either the expansion coefficients in a series in terms of shifted Legendre polynomials, or the approximate values of the desired original at a number of points. The first step of reduction to SLAE is to apply quadrature formulas that provide the minimum values of the condition number of SLAE. Regularization methods are used to obtain a reliable solution of the system. A common strategy is to use the Tikhonov stabilizer or its modifications. A variant of the regular-ization method for systems with oscillatory-type matrices is presented, which significantly reduces the conditionality of the problem in comparison with the classical Tikhonov scheme. A method is given for actually constructing special quadratures leading to problems with oscillation matrices.
Keywords: system of linear algebraic equations, integral equations of the first kind, ill-posed problems, ill-conditioned problems, condition number, oscillation matrices, 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. Ryabov V. M. Numerical inversion of the Laplace transform. St Petersburg, St Petersburg University Press (2013). (In Russian)
4. Cohen A.M. Numerical methods for Laplace transform inversion. New York, Springer (2007).
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. Ivanov V. K., Vasin V. V., Tanana V. P. Theory of linear ill-posed problems and its applications, Moscow, Science Publ. (1978). (In Russian)
8. Kabanikhin S.I. Inverse and ill-posed problems. Novosibirsk, Sib. Science Publ. (2009). (In Russian)
9. Voevodin V. V., Kuznetsov Yu. A. Matrices and computations. Moscow, Science Publ. (1984). (In Russian)
10. 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).
11. Gantmakher F.R. Teoriv matriz, Moscow, Nauka Publ. (1967) (In Russian) [Eng. transl.: Gant-makher F. R. The Theory of Matrices, Chelsea Publishing Company (1989)].
12. Lebedeva A.V., Ryabov V. M. On the numerical solution of system of linear algebraic equations with ill-conditioned matrices Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy 6 (64), iss.4, 619-629 (2019). https://doi.org/10.21638/11701/spbu01.2019.407 (In
*This paper was prepared with the support by a grant from St Petersburg State University "Event 3" (Pure ID 75207094).
Russian) [Eng. transl.: Vestnik St Petersburg University, Mathematics 52, iss.4, 388—393 (2019). https://doi.org/10.1134/S1063454119040058].
13. 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, Mathematics, 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, iss.4, 361-365 (2021) https://doi.org/10.1134/S1063454121040129].
14. Higham N. J. Functions of matrices: Theory and computation. Philadelphia, Society for Industrial and Applied Mathematics (2008).
Received: May 5, 2022 Revised: June 3, 2022 Accepted: June 9, 2022
Authors' information:
Anastasia V. Lebedeva — a.v.lebedeva@spbu.ru Victor M. Ryabov — v.ryabov@spbu.ru
ХРОНИКА
27 апреля 2022 г. на заседании секции теоретической механики им. проф. Н. Н. По-ляхова в санкт-петербургском Доме ученых им. М. Горького выступили д-р физ.-мат. наук, проф. М. П. Юшкова и магистрант С. О. Бондаренко (математико-механический факультет СПбГУ) с докладом на тему «Применение принципа максимума Понтря-гина и обобщенного принципа Гаусса для гашения колебаний тележки с тройным математическим маятником».
Краткое содержание доклада:
Рассматривалась задача о гашении колебаний горизонтально перемещающегося носителя, несущего тройной математический маятник. Требовалось найти оптимальную горизонтальную силу, приложенную к носителю, переводящую рассматриваемую механическую систему за указанное время из состояния покоя в новое состояние покоя с условием перемещения носителя на заданное расстояние. Решения осуществлялись с применением принципа максимума Понтрягина и с использованием обобщенного принципа Гаусса. Показаны преимущества применения второго метода по отношению к первому. В частности, ставя и решая расширенную краевую задачу, можно построить сколь угодно плавное движение системы в начале и в конце ее движения. В отличие от этого при первом методе всегда существуют скачки управления в эти моменты времени.
