ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2014, том 54, № 11, с. 1707-1710
УДК 519.61
АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ ОДНОГО КЛАССА КВАДРАТИЧНЫХ МАТРИЧНЫХ УРАВНЕНИЙ
© 2014 г. Ю. О. Воронцов, X. Д. Икрамов
(119992 Москва, Ленинские горы МГУ, ВМ и К) e-mail: vv@cs.msu.su; ikramov@msu.su Поступила в редакцию 05.03.2014 г.
Установлена связь между решениями квадратичного матричного уравнения XTDX + AX + + X^B + C = 0 с матричными n х «-коэффициентами, с одной стороны, и нейтральными под-
пространствами 2n х 2«-матрицы M =
C A B D
, с другой. Эта связь используется в предлагае-
мом алгоритме решения матричных уравнений указанного вида. Приведены результаты численных экспериментов, выполненных с использованием данного алгоритма. Библ. 3. Фиг. 2.
Ключевые слова: квадратичное матричное уравнение, нейтральное подпространство, разложение Такаги, вычислительный алгоритм.
DOI: 10.7868/S004446691411012X
1. Рассмотрим квадратичное матричное уравнение
XTDX + AX + ХтВ + C = 0,
(1)
в котором коэффициенты А, В, С, Б и искомое решение X суть п х п-матрицы. Сопоставим этому уравнению 2п х 2п-матрицу
M =
C A B D
Предположим, что уравнение (1) совместно, и пусть Х0 — некоторое его решение:
X DX0 + AX0 + X B + C = 0.
Положим
Z =
v x;
(2)
(3)
(4)
Тогда (3) можно переписать в виде
z0 mz0 = 0.
(5)
Будем рассматривать С2п как комплексное евклидово пространство, т.е. определим скалярное произведение векторов х = (х1, х2, ..., х2п) и у = (у1, у2, ..., у2п) как
(х, у) = Х1У1 + Ху + ... + Х2пУ2п.
Обозначим через подпространство в С2п, для которого является базисной матрицей. Из равенства (5) вытекает, что
(Mx, y) = 0
(6)
1707
2*
1708
ВОРОНЦОВ, ИКРАМОВ
для любых векторов х, у е Будем называть подпространство с таким свойством нейтральным подпространством матрицы М. Таким образом, разрешимость уравнения (1) влечет за собой наличие п -мерного нейтрального подпространства у ассоциированной с этим уравнением матрицы (2). Соответственно, наличие нескольких решений означает существование у M различных ^мерных нейтральных подпространств. Напротив, отсутствие таких подпространств указывает на неразрешимость уравнения (1).
Пусть, наоборот, — нейтральное подпространство размерности п матрицы М, а
г =
г, ^
)
(7)
есть его базисная матрица, разбитая на п х п-6локи ^ и Z2. Если блок Z1 невырожден, то матрицу Z можно заменить другой базисной матрицей того же подпространства
Ж =
V X )
X = .
(8)
Нейтральность подпространства означает, что
WтMW = 0,
откуда следует: блок X есть решение уравнения (1). Итак, всякое п-мерное нейтральное подпространство матрицы Мдает решение исходного квадратичного уравнения при дополнительном условии невырожденности верхнего блока Z1 в любой базисной матрице.
Если в матрице (7) невырожден блок Z2, то таким же образом можно показать разрешимость уравнения
У* СУ + BY + УтА + D = 0. (9)
Если невырожденны и Z1, и Z2, то порождаемые ими решения X и Гуравнений (1) и (9) являются взаимно обратными матрицами.
2. В остальной части статьи мы рассматриваем только симметричный случай уравнения (1). Матричные коэффициенты такого уравнения подчинены соотношениям
В = Ат, С = Ст, D = Dт.
(10)
Матрица М, ассоциированная с таким уравнением, симметрична.
В [1] предложено строить нейтральное подпространство путем приведения матрицы М к матрице N
(
N = 0ТИ0 =
N,1 N,2
V N,2 N22 )
(11)
в которой блок порядка п нулевой. Укажем один из возможных алгоритмов вычисления решения X уравнения (1), использующий такое приведение.
1. Приведение матрицы М к диагональному виду. На этом этапе ищется унитарная матрица Р, трансформирующая Мв матрицу РТМР = Л, где Л = ё1а§(Х, Х2, ..., Х2п) — диагональная матрица с неотрицательными диагональными элементами. Средством вычисления матрицы Р служит разложение Такаги (см. [2, гл. 4, следствие 4.4.4]).
2. Обнуление блока матрицы N. Для каждого к, 1 < к < п, такого, что элемент в позиции (к, к) ненулевой, строится унитарная матрица Як, отличающаяся от единичной лишь в узлах опорного квадрата (к, к), (к, п + к), (п + к, к), (п + к, п + к) и трансформирующая главную 2 х 2-подматрицу матрицы N в матрицу с нулем в верхней диагональной позиции:
(
ск
V Якск )
'к о и
о х
\
ск
V Якск )
( \ о Ук
V Ук Пк )
АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ
1709
Элемент в позиции (1, 1) будет аннулирован, если выполнено равенство е2к Хк + $2к ХК + к = 0. Полагая
Ck =
получаем искомую матрицу Rk. 3. Матрица
X
n + к
Xk + Xn
и sk = I
Xk
Xk + Xn
Q = PRR2..R
(12)
преобразует М в матрицу N с нулевым блоком Первые п столбцов матрицы Q образуют орто-нормированный базис п-мерного нейтрального подпространства исходной матрицы М, и решение X вычисляется в соответствии с формулами (7), (8) (после проверки блока Z1 матрицы Z на невырожденность).
3. Описанный в предыдущем разделе алгоритм был реализован в виде функции Quad языка Matlab версии R2011b. Обсуждаемые ниже эксперименты использовали именно эту функцию. Ввиду отсутствия в системе Matlab стандартной функции, реализующей разложение Такаги, было решено использовать для экспериментов матричные коэффициенты уравнения (1), не только удовлетворяющие условиям (10), но и делающие матрицу Мунитарной. В этом случае для вычисления разложения Такаги унитарной симметричной матрицы М можно использовать конечный алгоритм, предложенный в [3].
Фиг. 1.
Фиг. 2.
1710
ВОРОНЦОВ, ИКРАМОВ
Поскольку все этапы алгоритма основаны на ортогональных процедурах, нет оснований сомневаться в его численной устойчивости. Наши эксперименты подтверждают эту априорную оценку.
Была проведена массовая серия расчетов для уравнений порядка n = 100. Матрица М генерировалась как псевдослучайная унитарная матрица U, которая затем превращалась в симметричную (и по-прежнему унитарную) по формуле М = UUT. Из полученной матрицы выделялись матричные коэффициенты А, В, С и D уравнения (1).
Точное решение полученного таким образом уравнения неизвестно, поэтому качество вычисленного решения X оценивалось посредством (евклидовой) нормы невязки
R (X) = X dx + AX + X B + C.
Типичное значение нормы невязки (полученное усреднением по 100 уравнениям) составляло 5.7378 х 10-10. Время решения одного уравнения при том же усреднении составляло 0.7241 с. Эксперименты проводились на персональном компьютере с процессором Intel Core i7 2600K и 16 гигабайтами оперативной памяти.
Еще одна серия экспериментов была проведена с тем, чтобы проследить ухудшение качества вычисленного решения по мере возрастания порядка n матричных коэффициентов. Матричные коэффициенты уравнения (1) генерировались тем же способом, но теперь их порядок n возрастал от 50 до 950 с шагом 50. Для каждого из этих значений n была вычислена норма невязки с усреднением по 10 однотипным уравнениям. Даже для уравнений порядка n = 950 невязка остается очень малой (ее норма есть величина порядка 10-7). На фиг. 1 норма невязки представлена как функция от n. Зависимость времени (в секундах), затраченного на решение одного уравнения, от порядка n показана на фиг. 2.
СПИСОК ЛИТЕРАТУРЫ
1. Икрамов X. Д. О разрешимости одного класса квадратичных матричных уравнений // Докл. АН. 2014. Т. 455. № 2. С. 135-137.
2. Хорн Р., Джонсон Ч. Матричный анализ. М.: Мир, 1989.
3. Икрамов X. Д. Разложение Такаги симметричной унитарной матрицы как конечный алгоритм // Ж. вычисл. матем. и матем. физики. 2012. Т. 52. № 1. С. 4-7.
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.