ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2014, том 54, № 6, с. 901-904
УДК 519.61
АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ ОДНОГО КЛАССА ПОЛУТОРАЛИНЕЙНЫХ МАТРИЧНЫХ УРАВНЕНИЙ
© 2014 г. Ю. О. Воронцов, X. Д. Икрамов
(119992 Москва, Ленинские горы, МГУВМК) e-mail: vv@cs.msu.su; ikramov@cs.msu.ru Поступила в редакцию 26.12.2013 г.
Установлена связь между решениями полуторалинейного матричного уравнения X*DX + + AX + X*B + C = 0 с матричными n х «-коэффициентами, с одной стороны, и нейтральными
подпространствами 2« х 2«-матрицы M = C A
, с другой. Эта связь используется в пред-
в б:
лагаемом алгоритме решения матричных уравнений указанного вида. Приведены результаты численных экспериментов, выполненных с использованием данного алгоритма. Библ. 1. Фиг. 2.
Ключевые слова: полуторалинейное матричное уравнение, нейтральное подпространство, квазиопределенная матрица, приведение матрицы к диагональному виду.
БО1: 10.7868/80044466914060167
1. Рассмотрим полуторалинейное матричное уравнение
Х*БХ + АХ + Х*В + С = 0, (1)
в котором коэффициенты А, В, С, Б и искомое решение X суть п х п-матрицы. Сопоставим этому уравнению 2п х 2п-матрицу
M =
( \ C A
V В Б У
Предположим, что уравнение (1) совместно, и пусть Х0 — некоторое его решение:
X* DX0 + AX0 + X*B + C = 0.
Положим
Z =
v Xo;
Тогда (3) можно переписать в виде
Z*MZ0 = o.
(2)
(3)
(4)
(5)
Обозначим через подпространство в для которого является базисной матрицей. Из равенства (5) вытекает, что
(Мх, у) = 0 (6)
для любых векторов х, у е Подпространство с таким свойством принято называть нейтральным (или изотропным) подпространством матрицы М. Таким образом, разрешимость уравнения (1) влечет за собой наличие п-мерного нейтрального подпространства у ассоциированной с этим уравнением матрицы (2). Соответственно наличие нескольких решений означает суще-
ствование у M различных «-мерных нейтральных подпространств. Напротив, отсутствие таких подпространств указывает на неразрешимость уравнения (1).
Пусть, наоборот, — нейтральное подпространство размерности п матрицы М, а
г =
г, ^
22 )
(7)
есть его базисная матрица, разбитая на п х «-блоки 21 и 22. Если блок ^ невырожден, то матрицу 2 можно заменить другой базисной матрицей того же подпространства
Ж =
V X )
X = 2221 .
(8)
Нейтральность подпространства означает, что
W*MW = 0,
откуда следует, что блок X есть решение уравнения (1). Итак, всякое «-мерное нейтральное подпространство матрицы М дает решение исходного полуторалинейного уравнения при дополнительном условии невырожденности верхнего блока 21 в любой базисной матрице.
Если в матрице (7) невырожден блок 22, то таким же образом можно показать разрешимость уравнения
У*СУ + БУ + У*Л + Б = 0. (9)
Если невырожденны и 21, и 22, то порождаемые ими решения Xи Ууравнений (1) и (9) являются взаимно обратными матрицами.
2. В остальной части статьи мы рассматриваем только частный случай "эрмитова уравнения" (1). Под свойством эрмитовости понимаются следующие требования к матричным коэффициентам:
В = А *, С = С*, Б = Б*. (10)
Это в точности условия эрмитовости ассоциированной с (1) матрицы М.
Назовем эрмитову т х т-матрицу Б квазиопределенной матрицей типа (к, I), если в ее блочном представлении
Б =
( \ Н К
К* -в )
\
матрицы Ни О положительно определены и их порядки равны, соответственно, к и I (к + I = т). В [1] доказана следующая
Теорема 1. Уравнение (1) разрешимо, если ассоциирования с ним матрица Мявляется квазиопределенной матрицей типа («, «).
Замечание. Условие этой теоремы равносильно добавлению к соотношениям (10) левнеровских неравенств
С > 0, Б < 0 или Б > 0, С < 0.
Предполагая, что для уравнения (1) выполнено условие теоремы 1, укажем один из возможных алгоритмов вычисления соответствующего решения X.
1. Приведение матрицы М к диагональному виду. На этом этапе ищется унитарная матрица Р, трансформирующая Мв матрицу Р*МР = Л, где Л = ё1а§(^1, А,2, ..., А,2п).
2. Упорядочивание собственных значений матрицы Л. На данном этапе находится матрица-перестановка и, сортирующая собственные значения матрицы Л следующим образом:
и* Л и =
( \
Л+ 0 0 Л_ )
(11)
АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ ОДНОГО КЛАССА 903
где диагональные блоки
Л+ = ^(сть Ст2, ...,а„), Л- = diag(хь Т2, ...,т„)
составлены, соответственно, из положительных и отрицательных собственных значений матрицы М.
3. Обнуление блока (1, 1) матрицы (11). Для каждого к (1 < к< п) строится вращение Як с опорным квадратом (к, к), (к, п + к), (п + к, к), (п + к, п + к), трансформирующее главную 2 х 2-под-матрицу в матрицу с нулем в верхней диагональной позиции:
(
Ck —Sk
V Sk
(
tfk о
(
V 0 xk Л -sk ck
( \ 0 Vk
VVk % У
22
Элемент в позиции (1, 1) будет аннулирован, если выполнено равенство ск стк + sk тк = 0. Нужные
параметры вращения определяются из формул s2k =
и Ck =
-Tk
- Tk
- Tk
4. Матрица преобразует M в матрицу N вида
Q = PUR, R2... Rn
N = Q *MQ =
0 N,
12
V N12 N22
в которой блоки М12 и М22 являются диагональными п х п-матрицами. Первые п столбцов матрицы О образуют ортонормированный базис п-мерного нейтрального подпространства исходной матрицы М, и решение Xвычисляется в соответствии с формулами (7), (8).
2
3. Описанный в предыдущем разделе алгоритм был реализован в виде функции Sesquilinear языка Matlab версии R2011b. В обсуждаемых ниже экспериментах использовалась именно эта функция.
Поскольку все этапы алгоритма основаны на ортогональных процедурах, нет оснований сомневаться в его численной устойчивости. Наши эксперименты подтверждают эту априорную оценку.
Была проведена массовая серия расчетов для уравнений порядка n = 100. Коэффициент A генерировался как матрица с псевдослучайными элементами, равномерно распределенными в круге радиусом 10, после чего полагалось B = A*. Коэффициенты C и — D генерировались как эрмитовы положительно определенные матрицы применением унитарных подобий к диагональным матрицам с псевдослучайными положительными элементами. Упомянем и другой возможный способ добиться положительной определенности: в псевдослучайной эрмитовой матрице, сконструированной по аналогии с A (и с учетом эрмитовости), надлежащим образом увеличиваются диагональные элементы (например, чтобы обеспечить диагональное преобладание в каждой строке).
Точное решение полученного таким образом уравнения неизвестно, поэтому качество вычисленного решения X оценивалось посредством (евклидовой) нормы невязки
R (X) = X* dx + AX + X*B + C.
Типичное значение нормы невязки (полученное усреднением по 100 уравнений) составляло 2.2948 х 10-11. Время решения одного уравнения при том же усреднении составляло 0.1755 с. Эксперименты проводились на персональном компьютере с процессором Intel Core i7 2600 K и 16 гигабайт оперативной памяти.
Еще одна серия экспериментов была проведена с тем, чтобы проследить ухудшение качества вычисленного решения по мере возрастания порядка n матричных коэффициентов. Матричные коэффициенты генерировались тем же способом, но теперь их порядок n возрастал от 50 до 950 с шагом 50. Для каждого из этих значений n была вычислена норма невязки с усреднением по
Фиг. 1. Фиг. 2.
10 однотипным уравнениям. Даже для уравнений порядка п = 950 невязка остается очень малой (ее норма есть величина порядка 10-8). На фиг. 1 норма невязки представлена как функция от п. Зависимость времени (в секундах), затраченного на решение одного уравнения, от порядка п показана на фиг. 2.
СПИСОК ЛИТЕРАТУРЫ
1. Икрамов Х.Д. О разрешимости одного класса полуторалинейных матричных уравнений // Докл. АН. 2014. Т. 454. № 3. С. 265-267.
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.