Журнал вычислительной математики и математической физики, 2021, T. 61, № 4, стр. 531-538

ВВЕДЕНИЕ

Среди инвариантных по отношению к оператору $A$ пространств особо важную роль играют корневые пространства. Напомним кратко в связи с этим понятием некоторые используемые ниже сведения.

В классическом определении (см., например, [1]) корневым вектором для оператора $A$, соответствующим собственному значению $\lambda $ кратности $m$, называется вектор $x$ такой, что ${{(A - \lambda I)}^{m}}x = 0$; $I$ – единичная матрица. Совокупность корневых векторов образует корневое пространство $\ker {{(A - \lambda I)}^{m}}$ размерности $m$ для этого значения $\lambda $. Все рассматриваемое пространство разбивается в прямую сумму корневых пространств, соответствующих различным собственным значениям оператора $A$. В результате любой вектор $x$ может быть разложен в сумму векторов, взятых из каждого корневого пространства и называемых проекциями на эти пространства. Матрица оператора $A$ в базисе, составленном из базисов корневых пространств, является блочно-диагональной.

В работе рассматривается вопрос о выделении корневых пространств для спектральной задачи в обобщенной постановке. Для такой задачи вводится понятие корневого вектора, которое совпадает с приведенным выше определением в случае классической спектральной задачи.

Под выделением конечномерного подпространства в линейном пространстве обычно понимается построение какого-либо его базиса. Для корневых пространств часто рассматривается задача вычисления жорданового базиса, в котором матрица оператора $A$ имеет наиболее простой (канонический) вид. Однако при решении ряда задач бывает нужно определение корневого пространства в целом, без исследования структуры жордановых цепочек. Для этого случая в работе рассматриваются способы выделения корневого пространства, соответствующего собственному значению $\lambda $, с помощью нахождения произвольного базиса. Эти способы не требуют вычисления собственных и присоединенных векторов. Такой подход в данном случае более оправдан, поскольку построение цепочек корневых векторов (особенно для немалых значений их длин) требует весьма кропотливых вычислений. Кроме того, для приближенно заданных матриц это построение может быть неустойчивой задачей, если элементы матрицы возмущаются в пределах допустимой погрешности (см., например, в [1]).

В работе показано (см. разд. 4), что рассматриваемые методы определения корневого пространства путем описанного построения произвольного (не жорданового) базиса являются в случае малых возмущений исходной задачи численно устойчивыми.

1. ПОСТАНОВКА ЗАДАЧИ

Пусть $A$ и $B$ – заданные комплексные $n \times n$-матрицы. Далее рассматривается обобщенная задача на собственные значения (см., например, [2]) вида

где $x$ есть $n$-столбец. Здесь предполагается, что матрицы $A$ и $B$ могут быть вырожденными.

Введем в рассмотрение функцию от $\lambda $:

очевидно, $\omega (\lambda )$ – многочлен степени не выше $n.$

Пусть ${{\lambda }_{0}}$ есть $m$-кратный корень уравнения $\omega (\lambda ) = 0,\quad m \geqslant 1.$ Возьмем ${{x}_{0}}$ – решение задачи (1.1) при $\lambda = {{\lambda }_{0}}$, т.е. ее собственный вектор, отвечающий этому значению ${{\lambda }_{0}}$. Соответствующими этому вектору присоединенными векторами будем называть $n$-столбцы ${{x}_{1}},\;{{x}_{2}},\;...,\;{{x}_{p}}$ такие, что для $\lambda $, близких к значению ${{\lambda }_{0}}$, имеет место представление

Совокупность векторов ${{x}_{0}},\;{{x}_{1}},\;{{x}_{2}},\;...,\;{{x}_{p}}$_, удовлетворяющих соотношению (1.2), называется жордановой цепочкой. В [3] для более общего, чем рассматриваемый здесь, случая, а именно, для нелинейной по спектральному параметру задачи доказано, что длина определенной таким образом цепочки не превосходит кратности собственного значения, т.е.

Если $B$ – единичная матрица, то данные определения совпадают с принятыми (см., например, [1]). Действительно, в этом случае соотношения между векторами жордановой цепочки задаются следующим образом:

Для краткости выкладок положим ${{\lambda }_{0}} = 0{\kern 1pt} $. Тогда, преобразуя левую часть (1.2) с учетом (1.4), получаем что и требовалось.

Методы вычисления собственных и присоединенных векторов в таком случае см., например, также в [1].

Пространство, порожденное собственными векторами задачи (1.1) при $\lambda = {{\lambda }_{0}}$ и всеми отвечающими им присоединенными векторами, называется корневым пространством этой задачи, соответствующим ${{\lambda }_{0}}$.

Далее, как сказано во введении, рассматриваются численные методы выделения корневого пространства, не требующие нахождения соответствующих собственных и присоединенных векторов.

2. ПРЕОБРАЗОВАНИЕ ИСХОДНОЙ ЗАДАЧИ

Зафиксируем нужные для дальнейшего свойства рассматриваемой задачи. Будем преобразовывать уравнение (1.1), умножая матрицы $A$ и $B$ слева и справа на невырожденные $n \times n - $-матрицы (одни и те же для $A$ и $B$). Нетрудно видеть, что умножение матриц $A$ и $B$ на такую матрицу $C$ слева, как следует из соотношений (1.1), (1.2), не меняет задачу. Умножение этих матриц на матрицу $C$ справа соответствует замене переменных

где ${{X}_{i}}$ – новые искомые $n$-столбцы. Функция $\omega (\lambda )$ при этом умножается на ненулевое число.

В дальнейшем изложении, преобразуя матрицы, мы для простоты не будем менять обозначений.

Пусть $B$ – вырожденная матрица: ${\text{rank}}\,B = r$, $r < n$. Тогда, используя указанные преобразования, приведем матрицу $B$ к виду

где $I$ – единичная $r \times r$-матрица. Разобьем получившуюся в результате преобразований матрицу $A$, а также новые искомые векторы $X$ на соответствующие блоки, В новых переменных задача (1.1) принимает вид Отметим, что ранг матрицы $\left\| {{{a}_{{21}}}\;{{a}_{{22}}}} \right\|$ равен $n - r,$ в противном случае $\omega (\lambda ) \equiv 0$.

Далее еще раз перейдем к новым искомым величинам, умножив матрицы $A$ и $B$ справа на невырожденную $n \times n$-матрицу так, чтобы получить ${{a}_{{21}}} = 0$, ${{a}_{{22}}} = I$; здесь, как и раньше, $I$ обозначает единичную матрицу соответствующего размера. Тогда вместо системы (2.1) мы получим задачу вида

где через ${{b}_{1}}$ и ${{b}_{2}}$ обозначены блоки преобразующей матрицы размеров $(r \times r)$ и $(r,n - r)$ соответственно. Таким образом, задача (2.1) перейдет в задачу

Задача (2.2) имеет ту же структуру, что и исходная задача (1.1), но рассматривается в пространстве меньшей размерности.

Если $\det {{b}_{1}} = 0$, то будем продолжать аналогичную обработку задачи, т.е. применять к задаче (2.2) преобразования, подобные тем, которые были применены к задаче (1.1). И так далее, пока не получим матрицу $b{}_{1},$ такую, что $\det {{b}_{1}} \ne 0$. Число таких этапов преобразований и переходов к обработке последующей задачи конечно, так как размерность пространства, в котором рассматривается задача, на каждом этапе уменьшается. В результате мы придем к системе, которая имеет вид

Эта задача имеет вид исходной постановки (1.1) с матрицами $A$ и $B$, которые выражаются следующим образом:

Для преобразованной описанным образом задачи (2.3) рассмотрим разложение матричной функции ${{(a + \lambda I)}^{{ - 1}}}$ в окрестности изолированной особой точки ${{\lambda }_{0}}$ по степеням $\lambda - {{\lambda }_{0}}$ (положительным и отрицательным), т.е. в ряд Лорана в кольце с центром ${{\lambda }_{0}}$. Тогда коэффициент ${{a}_{{ - 1}}}$ при степени ${{(\lambda - {{\lambda }_{0}})}^{{ - 1}}}$ является идемпотентной или проекционной матрицей и, таким образом, в некотором базисе представляет собой матрицу оператора проектирования на подпространство.

Напомним, что квадратная матрица $P$ называется проекционной, если ${{P}^{2}} = P$. В любом базисе проективному оператору соответствует проекционная матрица.

Сформулированное выше утверждение непосредственно следует из вида спектрального разложения для матричной функции ${{(a + \lambda I)}^{{ - 1}}}$, если привести матрицу $a$ к жордановой форме (см., например, [4, §§ 5.4–5.5]). Там же доказано, что образ идемпотентной компоненты спектрального разложения, в нашем случае это $\operatorname{Im} \;{{a}_{{ - 1}}}$, совпадает с $\ker {{(a + \lambda I)}^{m}}$, т.е. с соответствующим корневым пространством (см. введение). Таким образом, ${{a}_{{ - 1}}}$ является матрицей проектора на корневое пространство задачи (2.3), соответствующее $\lambda = {{\lambda }_{0}}.$

Из сказанного, в частности, также следует, что размерность этого корневого пространства равна $m$.

Далее для задачи (1.1) составим матрицу

и рассмотрим, как меняется эта матрица в результате описанных преобразований исходной задачи.

При умножении матриц $A$ и $B$ на квадратную невырожденную $n \times n\, - $-матрицу $C$ слева матрица $Q(\lambda )$ не меняется:

При умножении матриц $A$ и $B$ на такие квадратные $n \times n\, - $-матрицы справа матрица $Q(\lambda )$ меняется так, как должна меняться матрица линейного преобразования при переходе к новому базису:

т.е. новая матрица подобна исходной. Для окончательно полученных в преобразованной задаче (2.3) величин матрица $Q(\lambda )\quad $ имеет вид Отсюда следует, что при разложении определяемой равенством (2.5) матрицы $Q(\lambda )\,$ по степеням $\lambda - {{\lambda }_{0}}$ матрица ${{Q}_{{ - 1}}}$, т.е. коэффициент разложения при ${{(\lambda - {{\lambda }_{0}})}^{{ - 1}}}$, является матрицей проектора на корневое пространство для полученной после преобразования задачи (2.3). Но поскольку исходная матрица $Q(\lambda )\quad $ в результате описанной трансформации меняется как матрица линейного преобразования, то ${{Q}_{{ - 1}}}$ – коэффициент разложения при ${{(\lambda - {{\lambda }_{0}})}^{{ - 1}}}$ для исходной матрицы (2.4), – это проектор на соответствующее корневое пространство в исходных величинах, т.е. для задачи (1.1). В результате доказана следующая

Теорема. Составим для задачи (1.1) матричную функцию $Q(\lambda ) = {{(A + \lambda B)}^{{ - 1}}}B\quad $ и разложим ее в ряд Лорана в окрестности $\lambda = {{\lambda }_{0}}$. Тогда коэффициент ${{Q}_{{ - 1}}}$ при ${{(\lambda - {{\lambda }_{0}})}^{{ - 1}}}$ – это матрица проектора на корневое пространство рассматриваемой задачи, соответствующее $\lambda = {{\lambda }_{0}}$.

Далее для сокращения выкладок положим ${{\lambda }_{0}} = 0$ так, что указанное разложение для $Q(\lambda )$ будет иметь вид

Очевидно, степеней ниже, чем $ - m$ в разложении $Q(\lambda )$ быть не может.

3. ИЗЛОЖЕНИЕ МЕТОДОВ

Из доказанной теоремы вытекают следующие методы нахождения корневого пространства задачи (1.1), не требующие построения жордановых цепочек собственных и присоединенных векторов. Эти методы сводятся к построению произвольного базиса с помощью проекционной матрицы ${{Q}_{{ - 1}}}$.

3.1. Для нахождения матрицы ${{Q}_{{ - 1}}}$ возьмем соотношение

Из этого соотношения выведем совокупность формул, которые будем использовать для вычисления последовательности матричных коэффициентов указанного разложения

Для этого, подставляя в соотношение (3.1) приведенное выше разложение для $Q(\lambda )$ и приравнивая коэффициенты при одинаковых степенях $\lambda $, получаем следующую цепочку уравнений:

При решении системы (3.2), в отличие от случая, когда $A$ – невырожденная матрица, каждая величина ${{Q}_{s}}$ определяется не из одного очередного соотношения, а требует использования некоторой совокупности этих соотношений.

Остановимся на практической реализации этих вычислений.

Матрица_{${{Q}_{{ - m}}}$, удовлетворяющая первому уравнению, – это произвольная$n \times n$-матрица, столбцы которой порождают $\ker A$. Нa второе уравнение накладывается ограничение: ${{Q}_{{ - m}}}$ должна быть такой, чтобы это уравнение было разрешимо относительно матрицы ${{Q}_{{ - m + 1}}}$. Матрица ${{Q}_{{ - m + 1}}}$ определяется из второго уравнения с точностью до слагаемого, удовлетворяющего первому уравнению, но она должна быть такой, чтобы третье уравнение разрешалось относительно матрицы ${{Q}_{{ - m + 2}}}$ и т.д. При использовании полной совокупности уравнений все матрицы ${{Q}_{s}}$ определяются однозначно.}

Поскольку систему (3.2) нельзя решать последовательно (прогонкой), рассмотрим вопрос о том, где можно оборвать эту цепочку уравнений. Другими словами, сколько уравнений нужно использовать при практическом решении задачи этим методом, чтобы однозначно определить окончательную величину ${{Q}_{{ - 1}}}$.

Пусть матрицы ${{Q}_{{ - m}}},\; \ldots ,\;{{Q}_{N}}$ и ${{\tilde {Q}}_{{ - m}}},\; \ldots ,\;{{\tilde {Q}}_{N}}$ удовлетворяют первым $N + m + 1$ уравнениям системы (3.2). Пусть ${{L}_{s}} = {{Q}_{s}} - {{\tilde {Q}}_{s}}$. Пусть матрица ${{L}_{{ - 1}}} \ne 0$ и пусть $k$ – наименьший индекс, для которого все матрицы ${{L}_{s}} \ne 0;$ как было отмечено, $k \leqslant - 1$. Возьмем какой-либо $n$-столбец $d$, для которого ${{L}_{k}}d \ne 0$ и образуем $n$-столбцы ${{x}_{q}} = {{L}_{{q + k}}}d,\;q = 0,\; \ldots ,\;N - k$. Нетривиальные матрицы ${{L}_{s}}$, $k \leqslant s \leqslant N$, удовлетворяют однородной линейной системе, соответствующей указанной совокупности уравнений из (3.2):

Из этой системы для векторов ${{x}_{q}}$, $0 \leqslant q \leqslant N - k$, следуют очевидные соотношения При сравнении этих соотношений с формулами (1.4) для ${{\lambda }_{0}} = 0$ видно, что векторы ${{x}_{0}},\;{{x}_{1}},\; \ldots ,\;{{x}_{{N - k}}}$ образуют жорданову цепочку и, следовательно, для малых значений $\lambda $ имеет место представление Но в этом случае, как указано в разд. 1, выполняется соотношение $N - k + 1 \leqslant m$ (см. формулу (1.3)). Поэтому справедлива оценка $N \leqslant m + k - 1 \leqslant m - 2$. Следовательно, если взять с запасом $N = m - 1$, то матрица ${{Q}_{{ - 1}}}$ определяется однозначно.

Так что для вычисления окончательно определенного значения ${{Q}_{{ - 1}}}$ достаточно использовать первые $2m$ уравнений.

В результате для определения ${{Q}_{{ - 1}}}$ мы получили линейную систему $2m$ матричных уравнений следующего специфического вида:

Методы решения таких систем известны, на них не останавливаемся.

Получив матрицу ${{Q}_{{ - 1}}}$, находим максимальную совокупность ее существенно линейно независимых столбцов, которая и составляет базис искомого корневого пространства. Действительно, проекция $x' = {{Q}_{{ - 1}}}x$ любого вектора $x = {{({{x}_{1}},...,{{x}_{n}})}^{{\text{т}}}}$ на это корневое пространство представляет собой линейную комбинацию

столбцов ${{q}_{1}},\; \ldots ,\;{{q}_{n}}$ матрицы ${{Q}_{{ - 1}}}$. Эта формула, кстати, служит иллюстрацией того, что образ оператора называют также пространством столбцов его матрицы в некотором базисе.

3.2. Рассмотрим подобный предыдущему способ построения корневого пространства, использующий проекционную матрицу, вычисление которой основано на применении интегральных формул от матрицы. Известно, что интегральная формула Коши для аналитических функций скалярного аргумента, определенных на спектре матрицы, может быть распространена на функции от матрицы (см. [5]). В результате выражения коэффициентов ряда Лорана, в частности, интересующая нас формула вычисления вычета аналитической функции в изолированной особой точке, переносятся на матричные ряды. Таким образом, для преобразованной задачи с матрицей $Q(\lambda )$ из (2.5) справедлива следующая формула (см., например, [6], [7]).

Если взять замкнутый контур $\Gamma $, отделяющий ${{\lambda }_{0}}$ от остальных собственных значений и пробегаемый против часовой стрелки, то имеет место равенство

где матрица ${{a}_{{ - 1}}}$ – это коэффициент при степени ${{(\lambda - {{\lambda }_{0}})}^{{ - 1}}}$ для разложения в ряд Лорана матричной функции ${{(a + \lambda I)}^{{ - 1}}}\,$в окрестности ${{\lambda }_{0}}$. А поскольку, как показано в работе (см. разд. 2), матрица (2.5) получается подобным преобразованием исходной матрицы $Q(\lambda ) = {{(A + \lambda B)}^{{ - 1}}}B$, то соответствующий коэффициент ряда Лорана для матрицы $Q(\lambda )$ из (2.4) выражается следующим образом:

Согласно доказанной в работе теореме, ${{Q}_{{ - 1}}}$ в (3.4) является проектором на корневое пространство, соответствующее ${{\lambda }_{0}}$ для исходной задачи. Этой формулой можно воспользоваться для численного нахождения матрицы ${{Q}_{{ - 1}}}$. На методах вычисления такого интеграла (его определение см., например, в [6]) мы не останавливаемся. Отметим лишь, что формулы рекуррентного вычисления подынтегральной матрицы (“обобщенной резольвенты”), приведенные также в [6], могут быть полезными.

Как и раньше, в качестве базиса искомого корневого пространства выбирается совокупность существенно линейно независимых столбцов этой матрицы.

В [6] также показано, что если внутри контура $\Gamma $ расположено несколько собственных значений, то интеграл в формуле (3.4) представляет собой проектор на сумму корневых пространств, соответствующих этим собственным значениям. Это, аналогично сказанному ранее, непосредственно следует из формулы для матрицы ${{(a + \lambda I)}^{{ - 1}}}\,$ после приведения матрицы $a$ к жордановой форме и подобия матриц $Q(\lambda )$, определяемых формулами (2.4) и (2.5).

Относительно рассмотренных в разд. 3 методов сделаем следующее

Замечание. На практике можно использовать упрощение вычислений, если проводить предварительное преобразование задачи, реализуя путь, которым была доказана приведенная в работе теорема. Напомним, что в этом доказательстве мы сохраняли исходный размер $n \times n$ фигурирующих в задаче матриц, в то время как фактическая размерность рассматриваемого пространства уменьшалась. При практической реализации алгоритма переход от системы (2.1) к системе вида (2.2) можно осуществить следующим образом. Как было отмечено выше, ранг матрицы $(a{}_{{21}}\,,\;{{a}_{{22}}})$ во втором уравнении системы (2.1) равен $n - r$. Выбирая из этой матрицы невырожденный $(n - r) \times (n - r)$-минор, выразим соответствующие $n - r$ переменных через остальные. Результат запишем в виде

где $\varsigma $ – искомый $r$-столбец, а матрицы ${{P}_{1}}$ и ${{P}_{2}}$, соответствующих размеров $r \times r$ и $(n - r) \times r$, вычисляются. А тогда, подставив полученные выражения для $\xi $ и $\eta $ в первое уравнение (2.1), мы получим задачу имеющую вид исходного уравнения (1.1), принятого за основу, но с матрицами $A$ и $B$ меньшего размера $r \times r\;$.

4. ОБ УСТОЙЧИВОСТИ МЕТОДОВ

В заключение отметим еще раз целесообразность использования предложенных в работе методов в случае, когда для выделения корневого пространства достаточно определить какой-либо (не обязательно жорданов) базис. Как было отмечено, такой подход предпочтителен, поскольку построение жорданового базиса – более трудоемкая задача и к тому же может быть неустойчивой, если в процессе вычислений в матрицы $A$ и $B$ вносятся малые погрешности.

Известно, что при возмущении элементов матрицы в результате ошибок округления четкая грань между случаями простых и кратных собственных значений может стираться. При этом кратное собственное значение может “распадаться”, т.е. превращаться в очень близкие, но различные значения, и, наоборот, близкие собственные значения могут “слипаться”, т.е. становиться кратным собственным значением в пределах точности. Подобные ситуации нередко встречаются в вычислительной практике. В этих случаях канонический вид матрицы претерпевает качественное изменение. Так, для слипшихся собственных значений чисто диагональная матрица переходит к общей канонической форме, так что даже число собственных векторов может меняться скачкообразно. Понятно, что в этом случае нахождение жорданового базиса теряет смысл. В условиях, близких к описанной ситуации, построение такого базиса, очевидно, является неустойчивой задачей.

Рассмотренные в работе способы численного построения корневого пространства лишены описанного выше недостатка. Действительно, так как искомая матрица ${{Q}_{{ - 1}}}$ лишь приближенно представляет собой искомый проектор, а в результате выделения совокупности линейно независимых столбцов этой матрицы мы также приближенно получаем корневой базис, мы имеем дело с возмущенной задачей. Однако если, например, ${{\lambda }_{0}}$ – слипшееся собственное значение, то как следует из приведенного в разд. 3 описания методов, полученная как решением системы (3.3), так и вычислением интеграла (3.4) матрица ${{Q}_{{ - 1}}}$ определяет для некоторой близкой к исходной (возмущенной) пары матриц $A$ и $B$ проекцию на суммарное корневое пространство для соответствующих собственных значений, близких к слипшемуся значению ${{\lambda }_{0}}$. Это на практике и означает устойчивость.