Численное решение матричных уравнений Сильвестра с нормальными коэффициентами Текст научной статьи по специальности «Математика»

УДК 519.6

Х. Д. Икрамов1, Ю. О. Воронцов2

ЧИСЛЕННОЕ РЕШЕНИЕ МАТРИЧНЫХ УРАВНЕНИЙ СИЛЬВЕСТРА С НОРМАЛЬНЫМИ КОЭФФИЦИЕНТАМИ

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

Ключевые слова: непрерывное уравнение Сильвестра, дискретное уравнение Сильвестра, нормальный оператор, функции Matlab'a lyap и dlyap.

1. Введение. Матричные уравнения

АХ + ХВ = С, (1)

X - АХ В = С (2)

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

Существуют эффективные и численно устойчивые методы для решения матричных уравнений Сильвестра умеренного порядка, например, алгоритмы Бартелса-Стьюарта и Голуба Ibinn Ван Лоана. Функции Matlab'a lyap и dlyap являются хорошо продуманными реализациями этих алгоритмов. В п. 2 мы даем краткое описание структуры алгоритма Бартелса-Стьюарта.

Сопоставим уравнениям (1) и (2) линейные операторы

1-л.и :Х^АХ + ХВ, С л. и :Х ^ Х^ АХВ, (3)

действующие в пространстве МТО)П(С) комплексных т xn-матриц. Будем рассматривать МТО)П(С) как унитарное пространство со скалярным произведением, заданным правилом

(R,S)=tr(S*R). (4)

Зададимся вопросом: при каких условиях на матрицы А и В операторы (3) нормальны? Ответ на этот вопрос дан в п. 3. В п. 4 будет показано, как следует изменить алгоритм Бартелса-Стьюарта, чтобы использовать нормальность операторов (3). Численные результаты, представленные в п. 5, иллюстрируют превосходство модифицированных алгоритмов над исходными.

2. Алгоритм Бартелса^Стьюарта. Удобно разделить описание данного алгоритма на четыре этапа. Для определенности будем рассматривать уравнение (1). С незначительными и очевидными изменениями это же описание пригодно для алгоритма, применяемого к уравнению (2).

1. Приведение матриц А и В к форме Шура. На этом этапе вычисляются унитарные матрицы U и V, такие, что матрицы R = U* AU, S = V*BV являются (верхними или нижними) треугольными. Выбор конкретной комбинации треугольных форм (из четырех возможных вариантов) определяет организацию этапа 3.

2. Преобразование правой части. На этом этапе вычисляется матрица D = U*CV.

Результатом первых двух этапов является замена исходного уравнения (1) новым уравнением Сильвестра

RY + YS = D. (5)

1 Факультет ВМК МГУ, проф., д.ф-м.н., e-mail: ikramovQcs.msu.su

2 ООО "Глобус Медиа", инженер-разработчик, к.ф-м.н., e-mail: vvQcs.msu.su

Неизвестная матрица Y связана с X соотношением

Y = U* XV. (6)

3. Решение уравнения (5). Это уравнение можно рассматривать как блочно-треугольную систему, состоящую из тп линейных уравнений относительно тп коэффициентов у у матрицы Y. Если бы эта система не имела никакой специфики, то ее решение потребовало бы по меньшей мере 0((тп)2) арифметических операций. Однако специфика здесь имеется, и ее остроумное использование, предложенное Бартелсом и Стьюартом, снижает сложность этого этапа до 0(тп(т + п)) операций. Подробности реализации этого этапа можно найти, например, в [1, § 7.6].

4. Вычисление искомой матрицы X. В соответствии с (6) имеем X = UYV*.

Алгоритм Голуба Ibinn Ван Лоана имеет сходную структуру. Наиболее существенное его отличие от алгоритма Бартелса-Стьюарта касается этапа 1, где теперь одна из матриц А и В приводится к форме Хессенберга, а не к форме Шура. Как следствие, изменяется процесс решения системы линейных уравнений на этапе 3. Однако этот этап по-прежнему имеет сложность 0(тп(т + п)) операций.

3. Нормальные матричные коэффициенты. Так как Mm,n(С) — унитарное пространство со скалярным произведением (4), то операторы, сопряженные к операторам (3), задаются формулами

fa,B х А*х + А'/Г. G*a b : X X - А*ХВ*.

Теорема 1. Оператор Fa,b тогда и только тогда является нормальным, когда нормальны матрицы А и В.

Теорема 2. Предположим, что обе матрицы А и В ненулевые. Оператор Ga,b тогда и только тогда является нормальным, когда нормальны матрицы А и В.

Доказательства этих теорем приведены в [2, 3].

4. Модифицированный алгоритм Бартелса^Стьюарта. Нормальные матрицы отличаются от всех остальных следующим замечательным свойством.

Теорема 3. Форма Шура нормальной матрицы А есть диагональная матрица ее собственных значений.

Приведение к форме Шура осуществляется в системе Matlab функцией schür. QR-алгоритм, реализованный этой функцией, не способен извлечь пользы из информации о нормальности коэффициентов А и В. В частности, форма Шура, вычисляемая для нормальной матрицы А функцией schür, будет треугольной, а не диагональной матрицей. Однако теоретический факт, формулируемый теоремой 3, проявляется в том, что внедиагональные элементы этой формы, возникающие вследствие округлений, очень малы. Пусть, например, D — диагональная матрица порядка 2000, а N — нормальная матрица, полученная из D унитарным подобием с трансформирующей матрицей U. Если D' — треугольная матрица, полученная в арифметике двойной точности применением к N функции schür, то для типичной матрицы D с элементами, не превосходящими по модулю 10, справедливо неравенство ||D — D'\\p < Ю-11, где ||-||_р — норма Фробениуса.

Для уравнения (1) с нормальными коэффициентами А и В предлагается следующая модификация алгоритма Бартелса-Стьюарта: не изменяя первых двух этапов, выделить из матриц R и S их диагональные части Dr = diag(pi,... ,рт) и D$ = diag(ci,... ,от). Как следствие, уравнение (5) на третьем этапе заменится уравнением

DrY + YDS = D. Это уравнение решается тривиально, а именно:

Угз = ———> 1 ^ г ^ m, 1 < j < п.

Рг + Vj

Как видно, стоимость этого (в общем случае наиболее трудоемкого) этапа снижается до 0(тп) арифметических операций. Четвертый этап остается неизменным.

Аналогичным образом модифицируется метод Бартелса-Стьюарта для уравнения (2) с нормальными коэффициентами Л и IL

5. Численные результаты. Алгоритмы решения нормальных уравнений (1) и (2), описанные в предыдущем разделе, были реализованы в виде функций языка МаМаЬ (соответственно В8погта1 и 81>етпогта1). Было проведено сравнение времени работы этих функций с временем работы стандартных процедур МаМаЬ'а 1уар и сНуар.

Для этого сравнения был взят случай т = п. Матричные уравнения для тестов конструировались следующим образом: матричные коэффициенты А и В генерировались как нормальные п х п-матрицы путем произвольного унитарного преобразования произвольных диагональных матриц, а матрица С как п х п-матрица с псевдослучайными элементами, равномерно распределенными в круге радиуса 10.

Для каждого значения п решались десять нормальных уравнений типа (1), используя вначале 1уар, а затем В8погта1. Время решения, усредненное по этим десяти уравнениям, обозначим через ^(п) в случае 1уар и через ¿2(п) Для В8погта1. Начальное значение порядка п равнялось 50. В дальнейшем п возрастало с шагом 100 до значения 2950. На рис. 1 отношение tl{n)/ 1-2{п) показано как функция от порядка п.

Рис. 1

Точность решений оценивалась с помощью (евклидовой) нормы матричной невязки Я(Х) = = АХ + ХВ — С. С ростом п значение ||Д(Х)|| возрастает, при п = 3000 не превосходя 10~9 для 1уар и Ю-' для В8погта1. Этот уровень невязки считаем вполне приемлемым.

Аналогичные тесты были проведены для нормальных уравнений типа (2). При этом использовались функции сНуар и 81>етпогта1. Усредненное время решения обозначается через (п) в случае сНуар и через ¿4(п) для 81>етпогта1.

На рис. 2 показано отношение ¿3(п)/£4(п) как функция от п. Это отношение несколько меньше отношения ¿1 (п)/£2(п). что можно объяснить высоким быстродействием функции сНуар. По крайней мере, в тестах эта функция всегда работала быстрее, чем 1уар. Невязка в этом случае вычислялась по формуле ЩХ) = X — АХВ — С, и при п = 3000 ее норма не превосходила 10~8 для сНуар и 10~6 для 81>етпогта1.

Рис. 2

СПИСОК ЛИТЕРАТУРЫ

1. Голуб Дж.. Ван Лоун Ч. Матричные: вычисления. М.: Мир. 1999.

2. И к рам о в Х.Д. Условия нормальности матричных уравнений типа Сильвестра // Докл. АН. 2014. 459. № 4. С. 403 405.

3. И к рам о в Х.Д. Условия нормальности линейных матричных операторов типа Стейна // Докл. АН. 2015. 460. № 3. С. 209 271.

Поступила в редакцию 14.10.1G