Журнал вычислительной математики и математической физики, 2022, T. 62, № 5, стр. 742-756

1.1. Роль псевдоспектра при решении прикладных задач

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

Другое общее свойство прикладных задач связано с тем, что модуль непрерывности собственных значений несимметричных операторов ведет себя непредсказуемым образом. Для его исследования вводится понятие ε-спектра или псевдоспектра ${{\Lambda }_{\varepsilon }}(A)$ (см. [1], [3]). По определению к ε-спектру $n \times n$ матрицы A относятся те комплексные числа λ, для которых матрица $A - \lambda {{I}_{n}}$ “почти вырождена”

где ${{I}_{n}}$ – единичная матрица, ${{\sigma }_{{\min }}}$ – минимальное сингулярное число.

Визуализировать структуру $\varepsilon $-спектра можно, изображая график функции

например, с помощью линий уровня (см. примеры в разд. 5).

Размеры пятен ε-спектра определяют, насколько могут отличаться собственные значения $\left| {{{\lambda }_{j}}(A) - {{\lambda }_{j}}(A + \Delta A)} \right|$ исходной и возмущенной матриц, $\left\| A \right\| \leqslant \varepsilon $. Расположением пятен объясняются на первый взгляд парадоксальные примеры, когда решения системы $y{\text{'}} = Ay$ устойчивы, а возмущенной системы $y{\text{'}} = (A + \Delta A)y$ – неустойчивы. Этим же объясняется тот факт, что решения устойчивых систем могут на начальных временных интервалах показывать большой рост (“практическая” неустойчивость), что на практике может приводить к разрушениям инженерных конструкций, развитию турбулентности течений и т.д. [4].

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

1.2. Исследования псевдоспектра и критерии расположения спектра в заданной области

Движение в сторону альтернативных постановок задач можно видеть во многих работах. Так, целый ряд статей посвящен методам изображения псевдоспектра. С точки зрения вычислений это весьма затратная процедура. Поэтому предпринимаются попытки создать более быстрые алгоритмы для вычисления некоторого приближения к функции (1.1) (см., например, [5]). Авторы других работ концентрируются на конкретных операторах и численно исследуют, какие из собственных значений наименее устойчивы к возмущению матрицы [6]. Заметим, что изображение псевдоспектра матрицы не гарантирует, что на нем будут видны все собственные значения, лежащие в данной области. Некоторые пятна ε-спектра могут быть настолько небольшими, что их легко пропустить при вычислении значений функции (1.1) даже на сетке с довольно мелкой ячейкой.

Другой подход может быть связан с подбором некоторого численного критерия (или критериев), который определенным образом характеризует часть спектра, находящуюся в фиксированной области в целом. В качестве известного примера можно привести теорему Раусса-Гурвица о том, что все корни полинома с вещественными коэффициентами лежат строго в левой полуплоскости тогда и только тогда, когда положительны все главные миноры матрицы Гурвица [7]. Этот факт является одним из многих результатов, связывающих расположение корней полиномов и свойства определенных матриц и знакопеременных квадратичных форм, собранных в подробном обзоре [8]. В частности, там упоминается работа [9] о расположении всех корней полинома в угловой области $ - \theta < {\text{arg}}\;x < \theta $.

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

является проектором на инвариантное подпространство матрицы A, соответствующее собственным значениям ${{\lambda }_{j}}(A)$, лежащим внутри контура γ. Причем след проектора равен числу этих собственных значений с учетом кратности [1]. Можно сказать, что этот факт является матричным вариантом принципа аргумента.

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

1.3. Задача дихотомии матричного спектра

К конструктивному решению проблемы новой постановки спектральной задачи подошла научная группа под руководством С.К. Годунова. Ими была сформулирована задача о дихотомии матричного спектра. Суть ее состоит в том, чтобы определить

1) есть ли на заданной кривой точки спектра матрицы,

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

Заметим, что сходство названия задачи с ε-дихотомией пространства решений дифференциальных уравнений, описанной в [10], неслучайно, так как и в том и в другом случае производится разделение пространства на прямую сумму подпространств.

При этом критерием дихотомии и ответом на первый вопрос задачи является (возможно, с дополнительной нормировкой) интеграл

который сходится только при отсутствии собственных значений матрицы на γ, а норма которого зависит от расположения пятен псевдоспектра вблизи кривой γ. Можно видеть, что в отличие от матричных критериев [7], [8], матрица H является самосопряженной и положительно-определенной. Для ответа на второй вопрос задачи дихотомии используется интеграл (1.2). Заметим, что в случае разделения спектра относительно единичной окружности или мнимой оси матрицы H и P являются решениями уравнений Ляпунова и их обобщений [1], [11].

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

В данный момент разработаны алгоритмы дихотомии относительно окружности и мнимой оси, которые можно назвать базовыми, так как они являются основной составляющей частью всех других алгоритмов дихотомии [12]–[17]. Кроме того, существуют алгоритмы дихотомии относительно кривых второго порядка [18]–[20] и алгоритм дихотомии корней полинома относительно окружности [21]. Более подробный обзор алгоритмов дихотомии и их приложений представлен в статье [22].

2.1. Математическая основа

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

(здесь и далее обозначение нормы $\left\| {\; \cdot \;} \right\|$, примененное к вектору, означает евклидову норму, примененное к матрице – операторную норму). Основную роль при этом играет тот факт, что существование и единственность решения задачи (2.1) обусловлены расположением спектра матрицы A относительно единичной окружности. При этом разрешимость такой задачи тесно связана с разрешимостью задачи для матричной разностной функции Грина

Заметим, что утверждения, которые будут приведены ниже, в известной степени аналогичны утверждениям о разрешимости краевой задачи для системы линейных обыкновенных дифференциальных уравнений, заданной на бесконечном интервале (см. [14], [23]–[25]), правда, в последнем случае спектр разделяется относительно мнимой оси.

В первую очередь поставим вопрос об условиях единственности решений (2.1) и (2.2). Для ответа на него используется следующая лемма о разрешимости однородной задачи:

Лемма 1. Однородная краевая задача

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

Доказательство следует из представления общего решения однородных разностных уравнений в виде $U_{j}^{{[{\text{од}}]}} = {{A}^{j}}U_{0}^{{[{\text{од}}]}}.$ Решение уравнения из (2.1) будет нетривиальном в случае $\left\| {U_{0}^{{[{\text{од}}]}}} \right\| \ne 0$ и ограниченным, если вектор $U_{0}^{{[{\text{од}}]}}$ является линейной комбинацией собственных векторов, соответствующих собственным значениям $\left| {\lambda (A)} \right| = 1$. В противном случае нетривиальных решений однородной задачи нет.

Следующее утверждение определяет условия единственности решений задач (2.1) и (2.2).

Лемма 2. Решения краевых задач (2.1) и (2.2) единственны тогда и только тогда, когда собственные значения матрицы $A$ не лежат на единичной окружности

Теперь рассмотрим вопрос о существовании решений. Заметим, что если спектр матрицы A делится единичной окружностью на два непересекающихся множества, то с помощью подобного преобразования она может быть приведена к блочно-диагональной форме (см., например, [1]):

Лемма 3. Если выполнено условие (2.4), то решение задачи (2.2) существует и может быть представлено в явном виде

То, что заданная в (2.6) функция ${{G}_{k}}$ является решением задачи (2.2), проверяется непосредственной подстановкой.

Замечание. Если у матрицы A есть собственные значения, лежащие на единичной окружности, то представления, аналогичные (2.5) и (2.6), отличающиеся только нестрогими неравенствами для спектров подматриц ${{\Lambda }_{0}}$ или ${{\Lambda }_{\infty }}$, также имеют место. Однако условия ограниченности решения (2.2) и/или его единственности при этом нарушаются.

Лемма 4. Если выполнено условие (2.4), то решение задачи (2.2) существует и может быть представлено в виде

Сходимость ряда в формуле (2.7) для любой ограниченной разностной функции правой части ${{f}_{k}}$ является следствием леммы 3. Справедливость представления (2.7) для решения задачи (2.1) может быть проверена непосредственной подстановкой.

Приведенные выше леммы 1–4 представляют собой отдельные этапы доказательства следующей теоремы.

Теорема 1. Для существования и единственности решения задачи (2.1) необходимо и достаточно выполнение условия (2.4).

Рассмотрим величину

и заметим, что этот ряд также как и (2.7) сходится только при условии (2.4). А значит, можно сформулировать следующий критерий разделения спектра матрицы $A$ единичной окружностью как следствие теоремы (1).

Утверждение 1. Условие (2.4) выполняется тогда и только тогда, когда $\left\| H \right\| < \infty $.

Обратим внимание, что матрица ${{P}_{{{\text{in}}}}} = {{G}_{1}}$ является проектором на инвариантное подпространство матрицы A, соответствующее собственным значениям, лежащим внутри единичной окружности. Проектор на подпространство, соответствующее собственным значениям, лежащим вне окружности, имеет вид ${{P}_{{{\text{out}}}}} = {{I}_{n}} - {{P}_{{{\text{in}}}}}$. Размерность инвариантных подпространств вычисляется по формулам ${{n}_{{{\text{in}}}}} = {\text{tr}}\;{{P}_{{{\text{in}}}}}$, ${{n}_{{{\text{out}}}}} = {\text{tr}}\;{{P}_{{{\text{out}}}}}$. C помощью сингулярного разложения (см., в частности, [1])

вычисляются ортогональные базисы ${{U}_{1}}$, ${{U}_{2}}$ инвариантных подпространств и матрица перехода $T = [{{U}_{1}},{{U}_{2}}]$, которая приводит матрицу A к клеточно-диагональному виду (2.5).

Обобщая вышесказанное на случай матричных пучков $A - \lambda B$, можно сформулировать следующее

Утверждение 2. На единичной окружности отсутствуют собственные значения пучка $A - \lambda B$ и возможно представление

где ${{I}_{\Lambda }}$, ${{I}_{M}}$ – единичные матрицы, тогда и только тогда, когда существует единственное решение краевой задачи

Причем численным критерием качества дихотомии спектра $A - \lambda B$ является величина (2.8).

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

2.2. Дихотомия единичной окружностью и мнимой осью

Приведенные выше факты лежат в основе алгоритма дихотомии единичной окружностью, который, в свою очередь, является ключевой частью алгоритмов дихотомии относительно остальных кривых. Существует несколько вариантов этого алгоритма. Однако основной цикл каждого из них представляет собой метод удвоений, т.е. переход от уравнений ${{A}_{k}}{{U}_{j}} - {{B}_{k}}{{U}_{{j + {{2}^{k}}}}} = 0$ к уравнениям ${{A}_{{k + 1}}}{{U}_{j}} - {{B}_{{k + 1}}}{{U}_{{j + {{2}^{{k + 1}}}}}} = 0$ путем исключения промежуточного слагаемого при помощи QR-разложения соответствующей матрицы.

Заметим также, что условия сходимости основного цикла могут быть сформулированы по-разному. Можно следить за абсолютной сходимостью, как это сделано ниже, или за относительной. В условие можно также включать сходимость проектора ${{P}_{k}}$. Кроме того, существуют оценки, позволяющие априорно оценить число итераций, необходимых для сходимости алгоритма при условии $\left\| H \right\| \leqslant {{\omega }_{{{\text{max}}}}}$. Эту оценку также можно использовать в качестве верхней границы числа итераций. При этом нужно принимать во внимание, что она, как правило, очень завышена, так как в обычной ситуации для сходимости с высокой точностью достаточно нескольких итераций.

3.1. Используемые обозначения

Основными математическими объектами, которые будут рассматриваться далее, являются лучи (или полупрямые) на комплексной плоскости $\mathbb{C}$, начало которых будет находиться в точке 0. Обозначать лучи будем малыми буквами $a,\;b,\;c,\; \ldots $ . Угол, который образуют два луча a и b, будем обозначать $ab$ или $\angle ab$ (см. фиг. 1). Положительный отсчет меры угла будет направлен против часовой стрелки. Тогда, очевидно, верно равенство для сопряженных углов $\angle ab = 2\pi - \angle ba$.

Прямые, проходящие через точку $(0,\;0)$, будем обозначать двумя буквами $aa{\text{'}}$, где a и $a{\text{'}}$ – лучи с началом в точке (0, 0), составляющие вместе данную прямую. Такое двойное обозначение позволяет отождествлять прямую $aa{\kern 1pt} '$ и развернутый угол $aa{\text{'}}$. Кроме того, по аналогии с вещественной и мнимой осью для каждой прямой будем задавать направление. То есть один из лучей будем называть положительным, другой – отрицательным. Эти направления можно задавать или полностью произвольно, или исходя из удобства выкладок. При обозначениях мы будем учитывать направление прямой следующим образом: у прямой $aa{\text{'}}$ полупрямая a будет положительной, а полупрямая $a{\text{'}}$ – отрицательной. Так как основным инструментом для решения основной задачи будет являться дихотомия матричного спектра мнимой осью, которая делит плоскость на правую и левую полуплоскости, то по аналогии мы определим правые и левые полуплоскости для всех прямых, у которых мы определим направление. Правой назовем полуплоскость, в которой находится угол $a{\text{'}}a$ (от отрицательного направления $a{\text{'}}$ прямой $aa{\text{'}}$ к положительному a), и обозначим через $\mathbb{C}_{{aa'}}^{r}$. Соответственно вторую полуплоскость, в которой находится угол $aa{\text{'}}$ (от положительного направления к отрицательному), будем называть левой и обозначать через $\mathbb{C}_{{aa'}}^{l}$. В этих обозначениях часть плоскости внутри угла $ab$ является пересечением левой полуплоскости относительно прямой $aa{\text{'}}$ и правой полуплоскости относительно прямой $bb'$: $\mathbb{C}_{{aa'}}^{l} \cap \mathbb{C}_{{bb'}}^{r}$.

3.2. Дихотомия спектра матрицы произвольной прямой

Для того, чтобы исследовать расположение спектра матрицы $A$ относительно произвольной прямой $aa{\text{'}}$, проходящей через начало координат, необходимо совершить поворот комплексной плоскости на угол φ между прямой $aa{\text{'}}$ и мнимой осью $\xi = {{e}^{{i\varphi }}}\lambda $. При таком преобразовании прямая $aa{\text{'}}$ перейдет в мнимую ось, а исходная спектральная задача $Av = \lambda v$ в задачу для нового спектрального параметра ${{A}_{\varphi }}v = \xi v$, где ${{A}_{\varphi }} = {{e}^{{i\varphi }}}A$. Из этих равенств видно, что собственные векторы, а значит, и инвариантные подпространства у матриц A и ${{A}_{\varphi }}$ совпадают, а собственные значения отличаются на угол φ. А значит, дихотомия спектра ${{A}_{\varphi }}$ относительно мнимой оси позволит определить проекторы ${{P}^{r}}$ и ${{P}^{l}}$ на инвариантные подпространства матрицы A, соответствующие собственным значениям, лежащим в правой $\mathbb{C}_{{aa'}}^{r}$ и левой $\mathbb{C}_{{aa{\text{'}}}}^{l}$ полуплоскостях. С помощью сингулярного разложения получаем ортогональные базисы ${{U}^{r}}$, ${{U}^{l}}$в инвариантных подпространствах, после чего можем представить матрицу в клеточно-диагональном виде:

3.3. Проверка наличия собственных значений матрицы на луче

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

Обозначим через a луч с началом в точке (0, 0), а угол между ним и вещественной положительной полуосью $\alpha $. Причем отсчет угла производим от луча a к полуоси против часовой стрелки. Тогда задача о локализации спектра $n \times n$ матрицы A на заданном луче a при помощи поворота сводится к вопросу о расположении спектра матрицы ${{A}_{\alpha }} = {{e}^{{i\alpha }}}A$ на вещественной положительной полуоси.

Далее заметим, что среди собственных значений матрицы ${{A}_{\alpha }}$ отсутствуют положительные вещественные числа тогда и только тогда, когда среди собственных значений матричного пучка ${{A}_{\alpha }} - {{\xi }^{2}}{{I}_{n}}$, квадратичного относительно спектрального параметра ξ, нет вещественных чисел. Из равенства определителей

где следует, что спектр квадратичного матричного пучка ${{A}_{\alpha }} - {{\xi }^{2}}{{I}_{n}}$ совпадает со спектром матрицы A удвоенного размера.

И, наконец, очевидно, что наличие собственных значений матрицы A на вещественной оси эквивалентно наличию собственных значений матрицы $i{\mathbf{A}}$ на мнимой оси.

Итак, для того чтобы определить наличие или отсутствие точек спектра матрицы A на луче a, нужно вычислить критерий дихотомии относительно мнимой оси для матрицы $i{\mathbf{A}}$. Заметим, что вычисления проекторов в данном случае не происходит, так как луч не разделяет плоскость на две непересекающиеся области.

4.1. Формулировка задачи

На комплексной плоскости задан угол с вершиной в начале координат, образованный лучами a, b (см. фиг. 1). Задача заключается в том, чтобы разделить спектр заданной матрицы A на две части в зависимости от расположения относительно угла $ab$. Для этого нужно будет определить, лежат ли на сторонах угла $ab$ собственные значения матрицы A. В случае, если стороны угла свободны от точек спектра, нужно вычислить проекторы на инвариантные подпространства матрицы, соответствующие собственным значениям, лежащим, соответственно, внутри и вне угла $ab$, определить базисы этих инвариантных подпространств и привести матрицу A к клеточно-диагональному виду.

Так как угол состоит из двух лучей, то для того, чтобы определить критерий дихотомии спектра углом, нужно определить численный критерий отсутствия точек спектра на каждом луче отдельно (см. разд. 3). Если определены значения $\omega (a)$ и $\omega (b)$ для обеих сторон угла, то в качестве критерия дихотомии углом можно взять их сумму:

Далее нужно вычислить проекторы ${{P}_{{{\text{in}}}}}$ и ${{P}_{{{\text{out}}}}} = {{I}_{n}} - {{P}_{{{\text{in}}}}}$ на инвариантные подпространства, соответствующие собственным значениям, лежащим внутри и вне угла $ab$ соответственно, а также привести подобным преобразованием матрицу A к клеточно-диагональному виду

Здесь собственные значения подматрицы ${{A}_{{{\text{in}}}}}$ лежат внутри угла $ab$, а собственные значения ${{B}_{{{\text{out}}}}}$ – вне угла.

4.2. Дихотомия прямыми $aa{\text{'}}$ и $bb{\text{'}}$

Для выделения части спектра, лежащей внутри угла $ab$, можно использовать алгоритм дихотомии прямыми, являющимися продолжением сторон угла. Обозначим соответствующие прямые через $aa{\text{'}}$ и $bb{\text{'}}$.

Случай 1. На обеих прямых $aa{\text{'}}$ и $bb{\text{'}}$ отсутствуют собственные значения матрицы $A$. В этой ситуации алгоритм решения задачи дихотомии максимально прост, но, очевидно, неприменим в общем случае. Решение сводится к последовательному выполнению дихотомии спектра A относительно обеих прямых с вычислением проекторов $P_{{bb'}}^{l}$ и $P_{{aa{\text{'}}}}^{r}$. Тогда проектор на инвариантное подпространство, соответствующее собственным значениям, находящимся внутри угла $ab$, представляет собой произведение ${{P}_{{{\text{in}}}}} = P_{{bb{\text{'}}}}^{l}P_{{aa{\text{'}}}}^{r}$. Сингулярное разложение проекторов ${{P}_{{{\text{in}}}}}$ и ${{I}_{n}} - {{P}_{{{\text{in}}}}}$ позволяет получить базисы инвариантных подпространств и привести матрицу $A$ к клеточно-диагональному виду.

Случай 2. На прямой $aa{\text{'}}$ нет собственных значений матрицы A. Произведем дихотомию спектра односительно этой прямой (см. разд. 3). В результате чего множество всех собственных значений разделится на два непересекающихся подмножества, лежащих в полуплоскостях $\mathbb{C}_{{aa'}}^{l}$ и $\mathbb{C}_{{aa'}}^{r}$, разделенных прямой $aa{\text{'}}$. При этом будут вычислены проекторы

на инвариантные подпространства матрицы A размерности ${{n}_{1}} = {\text{tr}}\;{{P}_{1}}$ и $n - {{n}_{1}}$, соответствующих собственным значениям, лежащим в полуплоскостях, разделенных прямой $aa{\kern 1pt} '$, базисы в данных инвариантных подпространствах $U_{{aa'}}^{l}$, $U_{{aa{\text{'}}}}^{r}$ и матрицу перехода приводящую матрицу A к клеточно-диагональному виду:

Здесь спектр подматрицы ${{A}_{1}}$ полностью лежит в левой полуплоскости $\mathbb{C}_{{aa'}}^{l}$, а спектр ${{B}_{1}}$ – в правой $\mathbb{C}_{{aa'}}^{r}$

Луч $b$ лежит в левой полуплоскости $\mathbb{C}_{{aa'}}^{l}$ прямой $aa{\text{'}}$, тогда как его продолжение $b{\text{'}}$ лежит в правой полуплоскости. А внутренняя часть угла $ab$ представляет собой пересечение полуплоскостей $\mathbb{C}_{{aa'}}^{l} \cap \mathbb{C}_{{bb'}}^{r}$. Следовательно, для завершения решения задачи необходимо разделить спектр матрицы ${{A}_{1}}$ относительно прямой $bb{\text{'}}$. В результате применения метода дихотомии матричного спектра прямой будут получены проекторы

базисы $U_{{bb{\text{'}}}}^{r}$, $U_{{bb{\text{'}}}}^{l}$ и матрица перехода приводящую матрицу ${{A}_{1}}$ к клеточно-диагональному виду: где

Случай 3. На прямой $bb{\text{'}}$ отсутствуют собственные значения матрицы A. Если при этом критерий дихотомии спектра прямой $aa{\text{'}}$ бесконечен, то порядок действий поменяется. Сначала нужно провести дихотомию прямой $bb{\text{'}}$. Результат этого действия: проекторы

на инвариантные подпространства матрицы $A$, соответствующих собственным значениям, лежащим в полуплоскостях, разделенных прямой $bb{\text{'}}$, базисы в данных инвариантных подпространствах $U_{{bb'}}^{r}$, $U_{{bb'}}^{l}$ и матрицу перехода приводящую матрицу A к клеточно-диагональному виду (4.4), причем

Далее производится дихотомия спектра подматрицы ${{A}_{1}}$ относительно прямой $aa{\text{'}}$. Будут получены: проекторы

базис $U_{{aa'}}^{l}$, $U_{{aa'}}^{r}$ и матрица перехода такая, что имеет место представление (4.7), где

Итоговые матрицы проекторов и матричные разложения будут иметь следующий вид:

Здесь символом ${{0}_{{n - {{n}_{1}}}}}$ обозначена квадратная нулевая матрица размера $(n - {{n}_{1}}) \times (n - {{n}_{1}})$.

4.3. Случай невозможности дихотомии прямыми $aa{\text{'}}$ и $bb{\text{'}}$

Более сложной является ситуация, когда на обоих лучах $a{\text{'}}$ и $b{\text{'}}$, продолжающих стороны угла $ab$, находятся собственные значения матрицы A. В этом случае в алгоритм вводится дополнительный этап, целью которого является приведение исходной матрицы к такому клеточно-диагональному виду

где у матрицы ${{A}_{0}}$ нет собственных значений на прямых $aa{\text{'}}$ и $bb{\text{'}}$. Описанный ниже алгоритм является одним из возможных методов построения разложения (4.14).

Разделим угол $b{\text{'}}a$ на равные части лучами ${{c}_{1}},{{c}_{2}},\; \ldots ,\;{{c}_{{n - 1}}}$. Так как у матрицы $A$ всего $n$ собственных значений с учетом кратности, и на прямых $aa{\text{'}}$ и $bb{\text{'}}$ по предположению находятся собственные значения, то найдется прямая ${{c}_{k}}c_{k}^{'}$, на которой точки спектра матрицы A будут отсутствовать. Здесь $c_{k}^{'}$ – луч, продолжающий луч ${{c}_{k}}$, при этом внутренность угла $ab$ лежит в левой полуплоскости $\mathbb{C}_{{{{c}_{k}}c_{k}^{'}}}^{l}$. Таким образом, дополнительным этапом алгоритма является дихотомия спектра матрицы $A$ прямой ${{c}_{k}}c_{k}^{'}$. Результатом этого этапа являются матрицы проекторов

базисы $U_{{{{c}_{k}}c_{k}^{'}}}^{r}$, $U_{{{{c}_{k}}c_{k}^{'}}}^{l}$ и матрица перехода которая позволяет представить исходную матрицу в виде (4.14).

Далее к матрице ${{A}_{0}}$ применяется упрощенный алгоритм. В итоге будут получены следующие матрицы:

где ${{n}_{0}} = {\text{tr}}\;{{P}_{0}}$, ${{n}_{1}} = {\text{tr}}\;{{P}_{1}}$.

5.1. Определение сектора дифференциального оператора

В качестве первого примера мы рассмотрим спектральную задачу для уравнения Орра–Зоммерфельда для плоскопараллельного течения Пуазейля вязкой несжимаемой жидкости. Изучению спектральных характеристик оператора Орра–Зоммерфельда посвящено множество работ, например, уже упоминавшиеся [4], [6], [17]. В статье [26] реализована идея разложения собственных функций в ряды в граничных точках, а обширный список литературы, приведенный в этой работе, содержит ссылки на разносторонние исследования оператора Орра–Зоммерфельда. Среди них можно отдельно отметить [27]. Изложенный в данной работе подход перекликается с задачей дихотомии, так как алгоритм применяется к точкам границы области с некоторым шагом, а позволяет сделать вывод о числе собственных значений, лежащих внутри области.

С одной стороны, спектральная задача для уравнения Орра–Зоммерфельда имеет важное прикладное значение с точки зрения исследования развития турбулентности. С другой – это содержательный пример несамосопряженного дифференциального оператора, и исследование его спектра представляет собой нетривиальную, но при этом хорошо изученную задачу. Таким образом, оператор Орра–Зоммерфельда является идеальным объектом для апробации новых вычислительных методов.

Рассмотрим матричный пучок $A - \lambda B$, являющийся дискретной аппроксимацией оператора Орра–Зоммефельда [28], [29]

Здесь ${{D}_{2}}$, ${{D}_{4}}$ – матрицы коллокационных производных второго и четвертого порядков, учитывающие однородные условия Дирихле и Неймана, I – единичная матрица, α и ${{\kappa }^{2}} = {{\alpha }^{2}} + {{\beta }^{2}}$ – спектральные параметры, Re – число Рейнольдса, $U(y) = 1 - {{y}^{2}}$, $ - 1 \leqslant y \leqslant 1$ – профиль основного течения Пуазейля, ${\text{diag}}(U)$ – диагональная матрица, на диагонали которой стоят значения функции U в точках Гаусса-Лобатто ${{y}_{j}} = \cos \pi {\text{/}}n$. Далее будем полагать $n = 100$, ${\text{Re}} = 6000$, $\alpha = 1.02$, $\beta = 0$.

Часть $\varepsilon $-спектра матрицы $\mathcal{A} = {{B}^{{ - 1}}}A$ изображена на фиг. 2 в виде линий уровня функции $f(\lambda ) = \mathop {\log }\nolimits_{10} {{\sigma }_{{\min }}}(\mathcal{A} - \lambda I)$. Цифрами обозначены значения, которые принимает функция $f(\lambda )$ на данном контуре.

Пусть стороны угла $ab$ образуют с положительной вещественной полуосью углы $5\pi {\text{/}}4$ и $7\pi {\text{/}}4$. Рассмотрим семейство ломаных ${{a}_{t}}{{b}_{t}}$, полученных из угла $ab$ переносом его вершины вдоль мнимой оси в точку $(0,\;t)$ (см. фиг. 2). Очевидно, что дихотомия спектра матрицы $\mathcal{A}$ относительно угла ${{a}_{t}}{{b}_{t}}$ равносильна дихотомии спектра матрицы $\mathcal{A} - itI$ относительно угла $ab$. Для каждого t из заданного интервала вычислим критерий дихотомии $\omega ({{a}_{t}}{{b}_{t}})$ и проектор ${{P}_{{{\text{in}}}}}({{a}_{t}}{{b}_{t}})$ на инвариантное подпространство, соответствующее собственным значениям матрицы $\mathcal{A}$, находящимся внутри угла ${{a}_{t}}{{b}_{t}}$.

На фиг. 3a изображен график функции $g(t) = \mathop {\log }\nolimits_{10} \omega ({{a}_{t}}{{b}_{t}})$, значения которой конечны, если на сторонах угла ${{a}_{t}}{{b}_{t}}$ нет точек спектра матрицы $\mathcal{A}$, и бесконечны в обратном случае. То есть “пики” графика $g(t)$ указывают на значения t, при которых собственные значения $\mathcal{A}$ оказываются на стороне угла ${{a}_{t}}{{b}_{t}}$.

На фиг. 3б изображен график кусочно-постоянной функции $h(t) = {\text{tr}}\;{{P}_{{{\text{in}}}}}({{a}_{t}}{{b}_{t}})$. Значения этой функции совпадают с числом собственных значений матрицы $\mathcal{A}$, находящихся внутри угла ${{a}_{t}}{{b}_{t}}$. При прохождении стороны угла ${{a}_{t}}{{b}_{t}}$ через точки спектра функция меняет значение на величину, равную числу собственных значений, оказавшихся в этот момент на стороне угла. На представленном графике величина каждого такого изменения равна двум. Это согласуется с тем фактом (см., например, [28]), что кратность собственных значений на ветви P (группа собственных значений, для которых ${\text{Re}}\lambda \to 1$) равна двум, каждому из них соответствуют две собственные функции – четная и нечетная. Кроме того, при $t > {{t}_{0}} = 0.9284$ функция $h(t)$ тождественно равна размерности матрицы $n = 100$. Это означает, что сектор с вершиной в точке $(0,\;t)$, $t > {{t}_{0}}$, и углом полураствора $\pi {\text{/}}4$ содержит все собственные значения матрицы $\mathcal{A}$.

5.2. Случай дугообразного спектрального пятна

Описанный в предыдущем разделе алгоритм даст ложноотрицательный результат, если пятно ε-спектра матрицы для достаточно малого ε имеет вид дуги, охватывающей вершину угла. В этом случае может оказаться, что при конечном значении критерия дихотомии $\omega (ab)$ предварительное разделение спектра любой прямой $cc{\text{'}}$ невозможно, а значит, невозможно выполнение алгоритма в целом.

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

Рассмотрим пример двухдиагональной матрицы A размера $n \times n$, на главной диагонали которой стоят числа

и $A(n + 1,n + 1) = - 2$, а на побочной диагонали $A(j,j + 1) = 2$. Изображение ε-спектра этой матрицы при $n = 10$ и $n = 40$ приведено на фиг. 4.

Зафиксируем угол $ab$, стороны которого образуют с положительной вещественной полуосью углы $3\pi {\text{/}}4$ и $5\pi {\text{/}}4$, и применим алгоритм из разд. 4 дихотомии относительно угла $ab$ к этим матрицам. При этом положим ${{\omega }_{{\max }}} = {{10}^{{16}}}$. Результаты вычислений представлены в строках 1–5 табл. 1.

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

В строке 5 табл. 1 мы видим, что при $n = 40$ $\omega (ab) < {{\omega }_{{\max }}}$, т.е. спектральные пятна достаточно хорошо разделены сторонами угла. Однако вспомогательная дихотомия относительно одной из прямых ${{c}_{k}}c_{k}^{'}$ невозможна, так как $\omega ({{c}_{k}}c_{k}^{'}) \geqslant {{\omega }_{{\max }}}$. Пробелы в последних двух ячейках этой строки говорят о том, что при $\omega \geqslant {{\omega }_{{\max }}}$ дальнейшие вычисления не проводились.

Однако дихотомия спектра при $n = 40$ может быть осуществлена, если для дополнительной дихотоми вместо прямой выбрать окружность. В данном конкретном случае это окружность радиуса $R = 3$ с центром в точке (–3, 0). Результаты вычислений см. в табл. 1, строка 6.