Научная статья на тему 'Численное решение матричных уравнений Сильвестра в самосопряженном случае'

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

CC BY
133
22
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
НЕПРЕРЫВНОЕ УРАВНЕНИЕ СИЛЬВЕСТРА / ДИСКРЕТНОЕ УРАВНЕНИЕ СИЛЬВЕСТРА / САМОСОПРЯЖЕННЫЙ ОПЕРАТОР / SELF-ADJOINT OPERATOR / ФУНКЦИИ MATLAB'А LYAP И DLYAP / MATLAB FUNCTIONS LYAP AND DLYAP / CONTINUOUS-TIME SYLVESTER EQUATION / DISCRETE-TIME SYLVESTER EQUATION

Аннотация научной статьи по математике, автор научной работы — Икрамов Х. Д., Воронцов Ю. О.

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

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «Численное решение матричных уравнений Сильвестра в самосопряженном случае»

УДК 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.

i Надоели баннеры? Вы всегда можете отключить рекламу.