CC BY
25
УДК 519.6
Х. Д. Икрамов1, Ю. О. Воронцов2
ЧИСЛЕННОЕ РЕШЕНИЕ МАТРИЧНЫХ УРАВНЕНИЙ СИЛЬВЕСТРА В САМОСОПРЯЖЕННОМ СЛУЧАЕ
Алгоритмы типа Бартелса-Стьюарта для численного решения матричных уравнений Сильвестра умеренного порядка модифицируются для случаев, когда линейные операторы, ассоциированные с этими уравнениями, самосопряжены. Представлены численные результаты, иллюстрирующие превосходство модифицированных алгоритмов над исходными.
Ключевые слова: непрерывное уравнение Сильвестра, дискретное уравнение Сильвестра, самосопряженный оператор, функции Matlab'a lyap и dlyap.
1. Введение. Матричные уравнения
АХ + ХВ = С (1)
и
X - АХ В = С (2)
называются соответственно непрерывным и дискретным уравнениями Сильвестра. В обоих уравнениях А и В — квадратные матрицы, возможно, различных порядков тип, тогда как С и X суть матрицы размера тхп. Будем считать, что для этих уравнений выполнены хорошо известные условия однозначной разрешимости.
Существуют эффективные и численно устойчивые методы для решения матричных уравнений Сильвестра умеренного порядка, например алгоритмы Бартелса-Стьюарта и Голуба Ibinn Ван Лоа-на. Функции Matlab'a lyap и dlyap являются хорошо продуманными реализациями этих алгоритмов. В разделе 2 мы даем краткое описание структуры алгоритма Бартелса-Стьюарта. Сопоставим уравнениям (1) и (2) линейные операторы
1-л.и : X ^ АХ + ХВ (3)
и
Ga,b ■ X X - AXIL (4)
действующие в пространстве МТО)П(С) комплексных т х n-матриц. Будем рассматривать МТО)П(С) как унитарное пространство со скалярным произведением, заданным правилом
(R,S) = tr(S*R). (5)
Зададимся таким вопросом: при каких условиях на матрицы А и В операторы (3) и (4) самосопряжены? Ответ на этот вопрос дан в разделе 3.
В разделе 4 мы показываем, как следует изменить алгоритм Бартелса-Стьюарта, чтобы использовать самосопряженность уравнений (1) и (2). Численные результаты, представленные в разделе 5, иллюстрируют превосходство модифицированных алгоритмов над исходными.
2. Алгоритм Барт еле а Стьюарта. Удобно разделить описание данного алгоритма на четыре этапа. Для определенности будем рассматривать уравнение (1). С незначительными и очевидными изменениями это же описание пригодно для алгоритма, применяемого к уравнению (2).
1. Приведение матриц А и В к форме Шура. На этом этапе вычисляются унитарные матрицы U и V, такие, что матрицы
R = U*AU и S = V*BV (6)
являются (верхними или нижними) треугольными. Выбор конкретной комбинации треугольных форм (из четырех возможных вариантов) определяет организацию этапа 3.
2. Преобразование правой части. На этом этапе вычисляется матрица
D = U*CV. (7)
1 Факультет ВМК МГУ, проф., д.ф-м.н., e-mail: ikramovQcs.msu.su
2 Факультет ВМК МГУ, асп., e-mail: vvQcs.msu.su
Результатом первых двух этапов является замена исходного уравнения (1) новым уравнением Сильвестра
КУ + УБ = Б. (8)
Неизвестная матрица У связана с X соотношением
у = и* XV. (9)
3. Решение уравнения (8). Это матричное уравнение можно рассматривать как блочно-треуголь-ную систему, состоящую из тп линейных уравнений относительно тп коэффициентов у^ матрицы У. Если бы эта система не имела никакой специфики, то ее решение потребовало бы по меньшей мере 0((тп)2) арифметических операций. Однако специфика здесь имеется, и ее остроумное использование, предложенное Бартелсом и Стьюартом, снижает сложность этого этапа до 0(тп(т + п)) операций. (О подробностях реализации этого этапа можно прочесть, например, в [1, § 7.6].)
4. Возврат к исходной матрице X. В соответствии с (9) искомое решение X связано с уже вычисленной матрицей У соотношением
.V ГП". (10)
Алгоритм Голуба 11:>шп Ван Лоана имеет сходную структуру. Наиболее существенное его отличие от алгоритма Бартелса-Стьюарта касается этапа 1, где теперь одна из матриц А и В приводится к форме Хессенберга, а не к форме Шура. Как следствие изменяется процесс решения системы линейных уравнений на этапе 3. Однако этот этап по-прежнему имеет сложность 0(тп(т + п)) операций.
3. Самосопряженный случай. Напомним, что МТО)П(С) рассматривается как унитарное пространство со скалярным произведением (5). Нетрудно показать, что операторы, сопряженные к операторам (3) и (4), задаются формулами
рА,в -Х^А*Х + ХВ* (И)
и
& а,в ■ х х ~ А*ХВ*. (12)
Отсюда мы выводим следующие заключения.
Теорема 1. Оператор Ра,в тогда и только тогда является самосопряженным, когда А и В суть эрмитовы матрицы с точностью до диагональных сдвигов соответственно на а и —а. Здесь а — произвольное чисто мнимое число.
Теорема 2. Предположим, что обе матрицы А и В ненулевые. Оператор О а,в тогда и только тогда является самосопряженным, когда А и В суть эрмитовы матрицы с точностью до скалярных множителей с модулем 1.
Доказательства этих теорем даны в [2, 3].
4. Модифицированный алгоритм Бартелса^Стьюарта. Рассмотрим вначале уравнение (1). В самосопряженном случае измененный алгоритм Бартелса-Стьюарта состоит из следующих операций.
1. Приведение матриц А и В к диагональному виду. На этом этапе вычисляются унитарные матрицы 17 и V, такие, что матрицы
Я = II* А11 и 5 = \'*ВУ (13)
являются диагональными.
Замечание. Если программа, используемая для приведения к диагональному виду, требует, чтобы задаваемая матрица была эрмитовой, то требуется некоторая подготовка к этапу 1, а именно: нужно удалить чисто мнимую компоненту диагональных элементов с тем, чтобы главная диагональ матрицы стала вещественной. Позднее этот диагональный сдвиг компенсируется на этапе 3.
Второй этап и соотношение (9) остаются неизменными.
3. Решение уравнения (8). Это уравнение представляет собой теперь диагональную систему линейных уравнений относительно тп коэффициентов у у матрицы У. Процесс решения такой системы очевиден и требует лишь 0(тп) арифметических операций вместо прежнего числа 0(тп(т + п)) операций.
Четвертый этап остается неизменным.
Аналогичным образом изменяется метод Бартелса-Стьюарта для уравнения (2). От приведенного выше описания эта модификация отличается лишь двумя моментами:
• предварительная подготовка к этапу 1 (если она необходима) состоит в удалении комплексных скалярных множителей из .4 и В с тем, чтобы обе матрицы стали эрмитовыми. Позднее эти множители восстанавливаются на этапе 3;
• вместо уравнения (8) на этапе 3 решается уравнение
У - ДУ5 = £>.
Оно снова представляет собой диагональную систему линейных уравнений относительно тп коэффициентов //¿у. Процесс ее решения очевиден и требует 0(тп) арифметических операций.
5. Численные результаты. Алгоритмы решения самосопряженных уравнений (1) и (2), описанные в предыдущем разделе, были реализованы в виде функций языка МаМаЬ (соответственно ВБкк и Б^т««). Мы провели сравнение времени работы этих функций с временем работы стандартных процедур МаМаЬ'а 1уар и сНуар.
Для этого сравнения был взят случай т = п. Матричные уравнения для наших тестов конструировались следующим образом: матричные коэффициенты А, В и С генерировались как п х п-матрицы с псевдослучайными элементами, равномерно распределенными в круге радиуса 10. Кроме того, при построении матриц .4 и В учитывались условия теорем 1 и 2.
Для каждого значения п решались десять самосопряженных уравнений типа (1), вначале использовалась 1уар, а затем ВБкк. Время решения, усредненное по этим десяти уравнениям, обозначим через ¿1 (п) в случае 1уар и через ¿2(п) для ВБкк. Начальное значение порядка п равнялось 50. В дальнейшем п возрастало с шагом 100 до значения 2950. На рис. 1 отношение tl(n)/1-2(п) показано как функция от порядка п.
Аналогичные тесты были проведены для самосопряженных уравнений типа (2). При этом использовались функции сНуар и 31>ет88. Усредненное время решения обозначается через Ц(п) в случае сНуар и через ¿4(п) для 31>ет88. На рис. 2 показано отношение t■з(n)/¿4(п) как функция от п. Это отношение несколько меньше отношения t\(n)/t2(n), что можно объяснить высоким быстродействием функции сНуар. По крайней мере, в наших тестах эта функция всегда работала быстрее, чем 1уар.
СПИСОК ЛИТЕРАТУРЫ
1. Голуб Дж., Ван Лоун Ч. Матричные: вычисления. М.: Мир. 1999.
2. И к рам о в X. Д. Линейные: матричные уравнения типа Сильвестра в самосопряженном случае // Докл. РАН. 2013. 450. № 2. С. 147 149.
3. И к рам о в X. Д. Условия самосопряженности матричных уравнений типа Стейна // Докл. РАН. 2013. 451. № 5. С. 498 500.
Поступила в редакцию 11.11.13
60
ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2014. № 2
SOLVING SYLVESTER MATRIX EQUATIONS IN THE SELF-ADJOINT CASE
Ikramov Kh. D., Vorontsov Yu. O.
Algorithms of the Bartels-Stewart type for solving medium-size Sylvester matrix equations are modified for the case where the linear operators associated with these equations are self-adjoint. We present numerical results illustrating the superior performance of the modified algorithms compared to the original ones.
Keywords: continuous-time Sylvester equation, discrete-time Sylvester equation, self-adjoint operator, Matlab functions lyap and dlyap.
