В этом параграфе, который носит вспомогательный характер, мы опишем методику исследования устойчивости стационарных решений по линейному приближению. Как известно, устойчивость решения зависит от собственных чисел матрицы линеаризованной системы. Мы коротко рассмотрим методы нахождения коэффициентов характеристического многочлена матрицы и приведем критерий устойчивости Рауса – Гурвица. В дальнейшем мы используем характеристический многочлен при нахождении бифуркационных точек ( $\$ \$ 5.4,5.5,5.8$ ), а также для определения устойчивости периодических решений ( $\$ 5.8$ ). Наконец, мы напомним о методе нахождения собственных чисел матрицы непосредственно как корней характеристического многочлена, что является более предпочтительным при небольших размерах системы.

5.3.1. Линейные уравнения

Система линейных дифференциальных уравнений с постоянными коэффициентами
\[
\frac{d \mathbf{x}}{d t}=\mathbf{A} \mathbf{x}
\]

имеет (в случае невырожденной матрицы А) единственное стационарное решение $\mathbf{x}=0$. Если у матрицы А нет кратных

собственных значений, то решение системы (5.3.1) можно представить в виде
\[
\mathbf{x}(t)=\sum_{i=1}^{n} C_{i}[\ldots]_{i} e^{\operatorname{Re} \lambda_{i} t}
\]

где в квадратных скобках стоит либо собственный вектор матрицы $\mathbf{A}$, соответствующий вещественному собственному числу $\lambda_{i}$, либо ограниченная векторная функция переменной $t$, если $\lambda_{i}-$ мнимое число. Постоянные $C_{i}$ определяются начальным условием $\mathbf{x}(0)$. Из разложения (5.3.2) следует, что устойчивость (в данном случае глобальная) стационарного решения системы (5.3.1) определяется вещественными частями собственных чисел $\lambda_{i}$. Если для любого собственного числа матрицы А имеет место условие
\[
\operatorname{Re} \lambda_{i}<0, \quad i=1, \ldots, n,
\]

то нулевое стационарное решение является асимптотически устойчивым; в частности, для всякого решения системы (5.3.1) выполняется условие
\[
\lim _{t \rightarrow \infty}\|\mathbf{x}(t)\|=0 .
\]

Если какое-нибудь собственное число имеет положительную вещественную часть, то нулевое решение, очевидно, неустойчиво: существуют решения $\mathbf{x}(t)$, для которых $\|\mathbf{x}(t)\| \rightarrow \infty$ при $t \rightarrow \infty$ (и сколь угодно малых $\mathbf{x}(0)$ ). Эти утверждения справедливы и для случая кратных собственных чисел. При этом в квадратных скобках в разложении (5.3.2) появляются многочлены от $t$, однако, как известно, экспоненциальная функция убывает (или возрастает) быстрее, чем многочлен любого порядка.

5.3.2. Характеристический многочлен

Собственные числа матрицы А являются корнями характеристического многочлена
\[
\begin{aligned}
P(\lambda)=(-1)^{n} \operatorname{det}(\mathbf{A}-\lambda \mathbf{I})=\lambda^{n}+a_{1} \lambda^{n-1} & +\ldots \\
& \ldots+a_{n-1} \lambda+a_{n} .
\end{aligned}
\]

Для $n=2$ или 3 мы можем получить характеристический многочлен прямо из определения детерминанта. Для больших же значений $n$ этот процесс оказывается очень трудоемким и неудобным даже при использовании ЭВМ.

Существует целый ряд методов для вычисления коэффициентов характеристического многочлена (см., например, ра-

боты [5.6, 5.7, 5.8]). Ниже мы рассмотрим только два из них, наиболее простые с вычислительной точки зрения – это метод Крылова и метод интерполяции.

Метод Крылова основывается на тождестве Гамильтона Кэли
\[
P(\mathbf{A})=\mathbf{A}^{n}+a_{1} \mathbf{A}^{n-1}+\ldots+a_{n-1} \mathbf{A}+a_{n} \mathbf{I}=0 .
\]

Умножим тождество (5.3.6) на некоторый вектор $\mathbf{h}^{0}$; в результате мы получим соотношение
\[
a_{1} \mathbf{h}^{n-1}+a_{2} \mathbf{h}^{n-2}+a_{n-1} \mathbf{h}^{1}+a_{n} \mathbf{h}^{0}=-\mathbf{h}^{n},
\]

где через $\mathbf{h}^{k}$ обозначены произведения $\mathbf{h}^{k}=\mathbf{A}^{k} \mathbf{h}^{0}$. Эти векторы могут быть вычислены рекуррентно:
\[
\mathbf{h}^{k+1}=\mathbf{A h}^{k}, \quad k=0,1, \ldots, n-1 .
\]

Соотношение (5.3.7) представляет собой систему $n$ линейных алгебраических уравнений относительно $n$ неизвестных $a_{1}, \ldots, a_{n}$. Мы можем ее решить, например, методом исключения Гаусса. Если матрица системы (5.3.7) окажется вырожденной, что представляется маловероятным, то необходимо выбрать другой вектор $\mathbf{h}^{0}$. ${ }^{19]}$

Метод интерполяции основан на том факте, что два многочлена степени $n$, значения которых совпадают в $n+1$ различных точках, тождественно равны друг другу. Выберем $n+1$ различных вещественных чисел $s_{0}, s_{1}, \ldots, s_{n}$ и вычислим значение $P\left(s_{i}\right)$ :
\[
P\left(s_{i}\right)=(-1)^{n} \operatorname{det}\left(\mathbf{A}-s_{i} \mathbf{I}\right) .
\]

Отметим, что в формулу (5.3.9) входит определитель числовой матрицы ( $s_{i}$-заданное число), и для его вычисления можно использовать программное обеспечение, основанное на методе исключения Гаусса. Построив по точкам $\left(s_{i}, P\left(s_{i}\right)\right.$ ), $i=0, \ldots, n$, интерполяционный многочлен Лагранжа, мы получим выражение для характеристического многочлена $P(\lambda)$.

Для построения характеристического многочлена можно воспользоваться также методом неопределенных коэффициентов, записав (5.3.9) в виде системы линейных алгебраических уравнений для коэффициентов $a_{0}, a_{1}, \ldots, a_{n}$ :
\[
\begin{array}{c}
a_{0} s_{0}^{n}+a_{1} s_{0}^{n-1}+a_{2} s_{0}^{n-2}+\ldots+a_{n-1} s_{0}+a_{n}=P\left(s_{0}\right), \\
a_{0} s_{1}^{n}+a_{1} s_{1}^{n-1}+a_{2} s_{1}^{n-2}+\ldots+a_{n-1} s_{1}+a_{n}=P\left(s_{1}\right), \\
\cdots \cdots \cdot . \cdots \cdots \\
a_{0} s_{n}^{n}+a_{1} s_{n}^{n-1}+a_{2} s_{n}^{n-2}+\ldots+a_{n-1} \varsigma_{n}+a_{n}=P\left(s_{n}\right) .
\end{array}
\]

Коэффициент $a_{0}$ должен равняться единице (см. (5.3.5)). Поэтому отклонение $a_{0}$ от единицы можно рассматривать в качестве меры погрешностей (округления), возникающих в процессе вычислений. (Можно включить равенство $a_{0}=1$ в алгоритм вычислений и вычислять определитель (5.3.9) только $n$ раз.) ${ }^{[10]}$

5.3.3. Критерий Рауса – Гурвица

Для выяснения вопроса о том, все ли корни характеристического многочлена $P(\lambda)$ имеют отрицательные вещественные части, можно воспользоваться критерием Рауса – Гурвица [5.8]. Обозначим $D_{1}=a_{1}$ и
\[
D_{k}=\operatorname{det}\left[\begin{array}{ccccc}
a_{1} & a_{3} & a_{5} & \ldots & a_{2 k-1} \\
1 & a_{2} & a_{4} & \ldots & a_{2 k-2} \\
0 & a_{1} & a_{3} & \ldots & a_{2 k-3} \\
0 & 1 & a_{2} & \ldots & a_{2 k-4} \\
\ldots & \ldots & \ldots & \ldots & \ldots \\
0 & 0 & \ldots & \ldots & a_{k}
\end{array}\right]
\]

для $k=2,3, \ldots, n$ (если $j>n$, то считается, что $a_{i}=0$ ). Если теперь для всех $k$ выполнены неравенства
\[
D_{k}>0, \quad k=1,2, \ldots, n,
\]

то все корни многочлена $P(\lambda)$ имеют отрицательную вещественную часть.

Отметим, что определители $D_{i}$ представляют собой многочлены от коэффициентов $P(\lambda)$ и, следовательно, являются многочленами относительно элементов матрицы А. При небольших значениях $n$ эти многочлены можно выписать в явном виде.

5.3.4. Нахождение собственных чисел матрицы

Часто оказывается необходимым знать все собственные числа матрицы $\mathbf{A}$, а стало быть, и все (в том числе комплексные) корни характеристического многочлена. Вещественные корни можно найти, например, с помощью метода Ньютона
\[
\lambda^{(k+1)}=\lambda^{(k)}-\frac{P\left(\lambda^{(k)}\right)}{P^{\prime}\left(\lambda^{(k)}\right)}, \quad k=0,1, \ldots,
\]

задаваясь при этом подходящим образом выбранными начальными приближениями $\lambda^{(0)}$. Если мы таким способом найдем все $n$ вещественных корней, то задача решена. Если у нас по-

лучится $n-1$ вещественных корней, то оставшийся корень также должен быть вещественным. Если найдены $n-2$ вещественных корней $\lambda_{1}, \lambda_{2}, \ldots, \lambda_{n-2}$, то можно разделить многочлен $P(\lambda)$ на произведение корневых сомножителей вида $\left(\lambda-\lambda_{1}\right)\left(\lambda-\lambda_{2}\right) \ldots\left(\lambda-\lambda_{n-2}\right)$, после чего квадратное уравнение даст нам оставшиеся два корня. Если, наконец, число вещественных корней оказывается меньше, чем $n-2$, то после деления на соответствующее произведение корневых сомножителей мы получаем многочлен более высокого порядка, чем 2 , решение которого уже нельзя построить в явном виде. Существует целый ряд методов нахождения набора корней вещественного многочлена. Одним из наиболее популярных среди них является метод Лина-Бэрстоу. Здесь мы рассмотрим лишь оснсвную идею этого метода.

Выделим из произведения корневых сомножителей многочлена $P(\lambda)$ два множителя, соответствующие либо двум вещественным, либо двум комплексно-сопряженным корням $\lambda_{i}$ и $\lambda_{k}$ :
\[
P(\lambda)=\left(\lambda^{2}-\left(\lambda_{j}+\lambda_{k}\right) \lambda+\lambda_{j} \lambda_{k}\right) \prod_{\substack{i
eq 1 \\ i
eq k}}\left(\lambda-\lambda_{i}\right)
\]

Коэффициенты $\left(\lambda_{j}+\lambda_{k}\right)$ и $\lambda_{j} \lambda_{k}$ при таком выборе будут вещественными.

Основная идея метода Лина-Бэрстоу состоит в последовательном нахождении этих «квадратичных корневых множителей». Представим многочлен $P(\lambda)$ в виде
\[
P(\lambda)=\left(\lambda^{2}-\alpha \lambda-\beta\right)\left(\lambda^{n-2}+b_{1} \lambda^{n-3}+\ldots+b_{n-2}\right)+A \lambda+B,
\]

где $\alpha$ и $\beta$ – параметры, которые мы хотим подобрать так, чтобы $A=B=0$. Значения $\alpha$ и $\beta$ определяют с помощью итерационной процедуры методом Ньютона.

Самостоятельной проблемой при этом остается выбор начальных значений $\alpha$ и $\beta$. Предположим, что $P(\lambda)$ имеет только простые корни, а именно, $2 k$ комплексных и $m$ вещественных корней $(2 k+m=n)$. Тогда существует $k+\left(\begin{array}{l}m \\ 2\end{array}\right)$ решений $(\alpha, \beta)$ для уравнений $A=0, B=0$. Таким образом, решений может оказаться относительно много, и поэтому вполне вероятно, что метод Ньютона со случайно выбранным начальным приближением $\alpha$ и $\beta$ будет сходиться. Произвольное решение системы $A(\alpha, \beta)=B(\alpha, \beta)=0$ позволяет понизить степень многочлена на два. Однако при таком последовательном понижении степени многочлена могут накапливаться погрешности (в частности, неточности вычисления $\alpha$ и $\beta$, погрешности округления

и т. п.), так что последние из найденных корней могут оказаться уже недостаточно точными. Эти погрешности можно скорректировать, применяя, например, метод Ньютона или метод Лина Бэрстоу для исходного многочлена.

Вычисление набора собственных чисел матрицы с помощью характеристического многочлена оказывается эффективным при небольших $n$. При больших значениях $n$, например при $n>20$ (в зависимости от конкретной задачи), происходит накопление погрешностей округления, и результаты могут оказаться очень неточными, а часто вовсе неприемлемыми. Поэтому для вычисления всех собственных чисел матрицы обычно используются прямые итерационные методы, являющиеся стандартной составной частью программного обеспечения ЭВМ (например, программы EISPACK, MATLAB и др.). В настоящее время чаще всего используется $Q R$-алгоритм (см., например, [5.34]). Используя это программное обеспечение, можно эффективно находить все собственные числа матрицы с числом строк $n$, измеряемым десятками.

5.3.5. Нелинейные уравнения – исследование устойчивости в линейном приближении

Обратимся вновь к нелинейной системе (5.1.1)
\[
x_{i}^{\prime}=f_{i}\left(x_{1}, \ldots, x_{n}, \alpha\right), \quad i=1, \ldots, n .
\]

Будем исследовать поведение решения в окрестности стационарного состояния $\overline{\mathbf{x}}$ при заданном $\alpha$. Введем вектор отклонений решения от $\overline{\mathbf{x}}$
\[
\delta=\mathbf{x}-\overline{\mathbf{x}} .
\]

Исходное нелинейное уравнение в окрестности точки $\overline{\mathbf{x}}$ аппроксимируем с помощью формулы Тейлора, сохраняя члены до первого порядка включительно (‘ $=d / d t)$ :
\[
\begin{aligned}
x_{i}^{\prime} & =\vec{x}_{i}^{\prime}+\delta_{i}^{\prime}=f_{i}(\mathbf{x}, \alpha)=f_{i}(\overline{\mathbf{x}}+\delta, \alpha) \approx \\
& \approx f_{i}(\overline{\mathbf{x}}, \alpha)+\sum_{j=1}^{n} \frac{\partial f_{i}(\overline{\mathbf{x}}, \alpha)}{\partial x_{j}} \delta_{j}, \quad i=1,2, \ldots, n .
\end{aligned}
\]

Для стационарного решения $\overline{\mathbf{x}}$ имеем $\overline{\mathbf{x}}_{i}^{\prime}=0, f_{i}(\overline{\mathbf{x}}, \alpha)=0$, и из разложения (5.3.17) мы получаем
\[
\delta_{i}^{\prime}=\sum_{j=1}^{n} a_{i j} \delta_{j}, \quad i=1,2, \ldots, n, a_{i j}=\frac{\partial f_{i}}{\partial x_{j}}(\overline{\mathbf{x}}, \alpha) .
\]

Система (5.3.18) представляет собой линейную систему с постоянными коэффициентами, для ее исследования можно воспользоваться результатами предыдущих пунктов. Так, если все собственные числа матрицы $\mathbf{A}=\left[a_{i j}\right]$ имеют отрицательные вещественные части, то $\boldsymbol{\delta} \rightarrow \mathbf{0}$ при $t \rightarrow \infty$ и, стало быть, принимая во внимание формулу (5.3.16), $\mathbf{x} \rightarrow \overline{\mathbf{x}}$ при $t \rightarrow \infty^{1)}$. Если некоторые собственные числа матрицы А имеют положительные вещественные части, то стационарное решение $\overline{\mathbf{x}}$ не будет устойчивым. Если, наконец, некоторые из собственных чисел имеют вещественные части, равные нулю, а остальные $\operatorname{Re} \lambda_{i}<0$, то линейное приближение не позволяет решить вопрос об устойчивости стационарного решения.

Результаты, полученные с помощью указанного метода линеаризации, были использованы для выделения областей устойчивости на диаграммах стационарных решений, приведенных в §5.2. Покажем, как подсчитываются значения $\lambda_{1}$ и $\lambda_{2}$ в табл. 5.2. Дифференцируя по соответствующим переменным правые части дифференциальных уравнений задачи 1, мы получаем (ср. формулу (5.2.17))
\[
\mathbf{A}=\left[\begin{array}{l}
-\Lambda-\mathrm{Da} E, \quad \mathrm{Da}(1-x) E /(1+\Theta / \gamma)^{2} \\
-\mathrm{Da} B E, \quad-\Lambda+\mathrm{Da} B(1-x) E /(1+\Theta / \gamma)^{2}-\beta
\end{array}\right],
\]

где $E=\exp (\Theta /(1+\Theta / \gamma))$. Соответствующий характеристический многочлен имеет вид
\[
P(\lambda)=\lambda^{2}-\left(a_{11}+a_{22}\right) \lambda+a_{11} a_{22}-a_{12} a_{21} .
\]

В соответствии с терминологией, принятой в теории динамических систем на плоскости, будем называть стационарное решение (состояние равновесия) узлом, седлом или фокусом в зависимости от значений собственных чисел $\lambda_{1}, \lambda_{2}$ :
седло: $\lambda_{1}, \lambda_{2}$ вещественные,
\[
\begin{array}{l}
\lambda_{1} \cdot \lambda_{2}<0 ; \\
\lambda_{1} \cdot \lambda_{2}>0 ;
\end{array}
\]

узел: $\lambda_{1}, \lambda_{2}$ вещественные,
фокус: $\lambda_{1}, \lambda_{2}$ комплексно-сопряженные ( $\lambda_{1}=\bar{\lambda}_{2}$ ).
Фокус и узел могут быть как устойчивыми, так и неустойчивыми, седло же всегда является неустойчивым. Соответствующие типы решений и значения собственных чисел $\lambda_{1}, \lambda_{2}$ представлены в табл. 5.3.

В заключение отметим, что результаты, найденные с помощью метода линеаризации, имеют локальный характер. Это означает, что они говорят лишь о поведении траекторий, которые начинаются «вблизи» стационарного состояния ${ }^{11}$. В нелинейном случае вокруг устойчивого стационарного состояния существует так называемая «область притяжения», т. е. область начальных условий, которые дают траекторию, приближающуюся к указанному состоянию. В случае одного устойчивого состояния такой областью притяжения может оказаться все пространство $\mathbb{R}^{n}$. Если одновременно существует несколько устойчивых стационарных состояний, то отдельные области притяжения разделяются гиперповерхностями, называемыми сепаратрисами. В случае $n=2$ они находятся сравнительно просто (это кривые на плоскости); при $n>2$ указанная задача становится гораздо более трудной.