Устойчивость стационарного решения системы (5.1.1) характеризуется собственными числами матрицы Якоби, определяемой правыми частями соответствующих уравнений. Если мы изменяем один из параметров системы, то вдоль ветви решения на соответствующей диаграмме решений характер устойчивости может изменяться лишь в точках, где собственное число переходит из левой половины комплексной плоскости в правую. Переход вещественного собственного числа через нуль обсуждался нами в §5.2. Если пара комплексно-сопряженных собственных чисел пересекает мнимую ось, то матрица Якоби все время остается невырожденной, и на диаграмме стационарных решений мы имеем регулярную точку данной ветви. Однако характер стационарного решения при этом переходе может измениться: устойчивое решение может стать неустойчивым (или наоборот). Точки диаграммы стационарных решений, в которых пара комплексно-сопряженных собственных чисел пересекает мнимую ось, называются точками комплексной бифуркации ${ }^{1)}$ или точками бифуркации Хопфа, по имени математика, опубликовавшего одну из основополагающих работ о характере решений в окрестности таких точек. Следующий существенный факт мотивирует разработку алгоритмов для нахождения точек комплексной бифуркации: в указанных точках (при выполнении определенных условий) от ветви стационарных решений отходит ветвь периодических решений. Этот параграф мы посвятим методам определения этих точек; алгоритмам для расчета и продолжения периодических решений будет посвящен $\S 5.8$.

5.5.1. Аналитические подходы

В случае некоторых задач невысокой размерности для нахождения точек комплексной бифуркации нет необходимости обращаться к численным методам. Можно воспользоваться аналитическими подходами, основанными на том, что мы легко раскрываем определители матриц порядка 2 и 3 и находим корни многочленов такой же степени ${ }^{2}$.
1) Этот термин не является общепринятым. – Прим. ред.
2) Корни многочлена третьей степени легко находятся, если среди них есть чисто мнимые.

Проиллюстрируем сказанное на примере двух задач: модели Лоренца (см. задачу 10) и модели реактора с перемешиванием для реакции 1-го порядка (см. задачу 1).
Модель Лоренца описывается системой уравнений
\[
\begin{array}{l}
\dot{x}=-\sigma x+\sigma y, \\
\dot{y}=-x z+r x-y, \\
\dot{z}=x y-b z,
\end{array}
\]

правые части которых имеют матрицу Якоби вида
\[
\mathbf{J}=\left[\begin{array}{crr}
-\sigma & \sigma & 0 \\
r-z & -1 & -x \\
y & x & -b
\end{array}\right] .
\]

Для нетривиального стационарного решения (см. формулу (P10-5))
\[
x=y= \pm[b(r-1)]^{1 / 2}, \quad z=r-1, \quad r>1
\]

построим характеристический многочлен матрицы $\mathbf{J}$
\[
P(\lambda)=\lambda^{3}+(\sigma+b+1) \lambda^{2}+b(r+\sigma) \lambda+2 \sigma b(r-1) .
\]

В точке бифуркации Хопфа этот многочлен должен иметь два взаимно сопряженных чисто мнимых корня, которые мы обозначим как $\pm i \lambda_{0}$. Третий (вещественный) корень обозначим $\lambda_{1}$. Тогда
\[
P(\lambda)=\left(\lambda+i \lambda_{0}\right)\left(\lambda-i \lambda_{0}\right)\left(\lambda-\lambda_{1}\right)=\lambda^{3}-\lambda_{1} \lambda^{2}+\lambda_{0}^{2} \lambda-\lambda_{1} \lambda_{0}^{2} .
\]

Сравнивая формулы (5.5.4) и (5.5.5), получим
\[
\begin{array}{c}
\sigma+b+1=-\lambda_{1}, \\
b(r+\sigma)=\lambda_{0}^{2}, \\
2 \sigma b(r-1)=-\lambda_{1} \lambda_{0}^{2} .
\end{array}
\]

Считая $\sigma$ и $b$ фиксированными, имеем систему трех уравнений относительно неизвестных $\lambda_{0}$, $\lambda_{1}$ и $r$. Если учесть, что $r>1$, то из уравнений (5.5.6) следует, что $\lambda_{1}<0$ и $\lambda_{0}
eq 0$. Перемножая: левые и правые части первых двух уравнений системы (5.5.6) и вычитая полученный результат из последнего уравнения (5.5.6), мы получаем соотношение ${ }^{1)}$
\[
(\sigma+b+1) b(r+\sigma)=2 \sigma b(r-1),
\]
1) Если вещественный многочлен $\lambda^{3}+A_{2} \lambda^{2}+A_{1} \lambda+A_{0}$ имеет чисто мнимый корень, то $A_{0}=A_{1} A_{2}$.

в которое входит только неизвестная $r$. Отсюда
\[
r^{+}=\frac{\sigma(\sigma+b+3)}{(\sigma-b-1)},
\]

тде $r^{+}$-значение параметра $r$ в точке комплексной бифуркации; значения переменных $x^{+}, y^{+}, z^{+}$подсчитываются по формулам (5.5.3).

Для системы двух дифференциальных уравнений с матрицей Якоби
\[
\mathbf{J}=\left[\begin{array}{ll}
a_{11} & a_{12} \\
a_{21} & a_{22}
\end{array}\right]
\]

характеристический многочлен $\mathbf{J}$ имеет вид
\[
P(\lambda)=\lambda^{2}-\left(a_{11}+a_{22}\right) \lambda+\left(a_{11} a_{22}-a_{21} a_{12}\right) .
\]

В точке комплексной бифуркации этот многочлен имеет два qисто мнимых корня:
\[
P(\lambda)=\left(\lambda-i \lambda_{0}\right)\left(\lambda+i \lambda_{0}\right)=\lambda^{2}+\lambda_{0}^{2} .
\]

Сравнение (5.5.9) и (5.5.10) дает
\[
\begin{aligned}
a_{11}+a_{22} & =0, \\
a_{11} a_{22}-a_{21} a_{12} & =\lambda_{0}^{2}>0 .
\end{aligned}
\]

Используем теперь описанный подход для отыскания точек комплексной бифуркации в задаче 1. Матрица Якоби для правых частей уравнений (P1-6) и (P1-7) имеет вид
\[
\mathbf{J}=\left[\begin{array}{lc}
-\Lambda-\mathrm{Da} E(\Theta) & \mathrm{Da}(1-x) E_{1}(\Theta) \\
-\mathrm{Da} B E(\Theta) & -\Lambda+\mathrm{Da} B(1-x) E_{1}(\Theta)-\beta
\end{array}\right],
\]

где введены обозначения $E(\Theta)=\exp (\Theta /(1+\Theta / \gamma))$ и $E_{1}(\Theta)=$ $=E(\Theta) /(1+\Theta / \gamma)^{2}$.
Первое из условий (5.5.11) приводит нас к уравнению
\[
-2 \Lambda-\mathrm{Da} E(\Theta)+\mathrm{Da}\left(B-\Theta-\frac{\beta}{\Lambda}\left(\Theta-\Theta_{c}\right)\right) E_{1}(\Theta)-\beta=0,
\]

где мы подставили $x=\Theta / B+\beta\left(\Theta-\Theta_{c}\right) / B \Lambda$ (см. формулу (5.2.3)). Уравнения, описывающие стационарное состояние, дают еще соотношение между $\Theta$ и $\mathrm{Da}$ (см. (5.4.16):
\[
\mathrm{Da}=\frac{\beta\left(\Theta-\Theta_{c}\right)+\Lambda \Theta}{\left(B-\Theta-\frac{\beta\left(\Theta-\Theta_{c}\right)}{\Lambda}\right) E(\Theta)} .
\]

Подставляя выражение для Dа в уравнение (5.5.13), мы получаем кубическое уравнение относительно $\Theta$ (при фиксированных $B, \Lambda, \beta, \Theta_{c}$ ) следующего вида:
\[
\begin{aligned}
\left(\beta\left(\Theta-\Theta_{c}\right)+\Lambda \Theta-\right. & \left.(2 \Lambda+\beta)(1+\Theta / \gamma)^{2}\right)\left(B-\Theta-\frac{\beta\left(\Theta-\Theta_{c}\right)}{\Lambda}\right)- \\
– & \left(\beta\left(\Theta-\Theta_{c}\right)+\Lambda \Theta\right)(1+\Theta / \gamma)^{2}=0 .
\end{aligned}
\]

Решение этого уравнения дает нам значения переменной $\Theta^{+}$, после чего по формуле (5.5.14) нетрудно найти значения параметра $\mathrm{Da}^{+}$, отвечающие точке комплексной бифуркации. В табл. 5.12 приведен пример такого расчета для двух значений параметра $B$ в предельном случае $\gamma \rightarrow \infty$ (когда (5.5.15) сводится к квадратному уравнению).
Таблица 5.12. Нахождение точек комплексной бифуркации для задачи 1 $\left(\gamma \rightarrow \infty, \Lambda=0,5, \beta=0,8, \Theta_{c}=0\right)$.

В приведенных примерах мы определяли точки бифуркации Хопфа для задач, размерность которых не превышала 3. Для задач, описываемых системами более чем из 3 дифференциальных уравнений, исследователям приходится обычно использовать различного рода численные подходы. Часть из них будет рассмотрена в следующих двух пунктах $[5.14,5.15,5.16]$.

5.5.2. Методы декомпозиции

В точке комплексной бифуркации (бифуркации Хопфа) $\left(x_{1}^{+}, \ldots, x_{n}^{+}, \alpha^{+}\right)$матрица Якоби $\mathbf{J}$ имеет два (взаимно) сопряженных чисто мнимых собственных числа $\lambda_{1}$ и $\lambda_{2}$, т. е.
\[
\lambda_{1,2}= \pm i \sqrt{\omega^{+}}, \quad \omega^{+}>0 .
\]

Запишем характеристический многочлен матрицы $\mathbf{J}$,
\[
P(\lambda)=\lambda^{n}+a_{1} \lambda^{n-1}+\ldots+a_{n-1} \lambda+a_{n},
\]

и представим его в виде
\[
P(\lambda)=\left(\lambda^{2}+\omega^{+}\right) P_{n-2}(\lambda),
\]

где $P_{n-2}(\lambda)$ – многочлен степени $n-2$. Для нахождения величины $\omega^{+}$(заранее неизвестной) перепишем $P(\lambda)$ в форме
\[
P(\lambda)=\left(\lambda^{2}+\omega\right)\left(\lambda^{n-2}+p_{1} \lambda^{n-3}+\ldots+p_{n-3} \lambda+p_{n-2}\right)+A \lambda+B .
\]

Коэффициенты $p_{1}, \ldots, p_{n}, A$ и $B$ могут быть вычислены по рекуррентным формулам
\[
\begin{array}{c}
p_{-1}=0, \quad p_{0}=1, \quad p_{k}=a_{k}-\omega p_{k-2}, \quad k=1,2, \ldots, n-2 \\
A=a_{n-1}-\omega p_{n-3}, \quad B=a_{n}-\omega p_{n-2} .
\end{array}
\]

Величины $p_{i}, A$ и $B$ зависят от $\mathbf{x}, \alpha$ и $\omega$. При этом для точки комплексной бифуркации выполняются соотношения
\[
\begin{array}{l}
f_{n+1}\left(x_{1}, \ldots, x_{n}, \alpha, \omega\right)=A\left(x_{1}, \ldots, x_{n}, \boldsymbol{\alpha}, \omega\right)=0, \\
f_{n+2}\left(x_{1}, \ldots, x_{n}, \boldsymbol{\alpha}, \omega\right)=B\left(x_{1}, \ldots, x_{n}, \boldsymbol{\alpha}, \omega\right)=0 .
\end{array}
\]

Таким образом, в общем случае мы имеем $(n+2)$ нелинейных уравнений (5.1.3), (5.5.20) относительно $(n+2)$ неизвестных $x_{1}, \ldots, x_{n}, \alpha$, $\omega$. Решение этих уравнений позволяет нам найти точку комплексной бифуркации. Для нахождения решения можно воспользоваться методом Ньютона, а для вычисления элементов матрицы Якоби в методе Ньютона можно частично (для нахождения производных от $f_{n+1}$ и $f_{n+2}$ ) использовать разностные формулы. Можно использовать и аналитические выражения, однако в случае больших $n$ это может оказаться весьма трудоемким.

Отметим, что коэффициенты $p_{1}, \ldots, p_{n-2}$ несут информацию об устойчивости стационарного решения в окрестности точки комплексной бифуркации (бифуркации Хопфа). Если все корни многочлена с указанными коэффициентами располагаются в левой полуплоскости (см. § 5.3), то стационарное решение будет устойчивым либо при $\alpha<\alpha^{+}$, либо при $\alpha>\alpha^{+}$.

Примеры применения описанного метода для трех различных точек комплексной бифуркации в случае задачи 2 представлены в табл. 5.13. При этом в качестве параметра $\alpha$ использовалось число Дамкелера в первом реакторе каскада $\mathrm{Da}_{1}$.

Если матрица $\mathbf{J}$ имеет собственное число $\lambda$, то матрица $\mathbf{J}^{2}$ обладает собственным числом $\lambda^{2}$. Для собственных чисел (5.5.16) мы имеем $\lambda_{1,2}^{2}=-\omega^{+}$. Итак, матрица $\mathbf{J}^{2}$ имеет отрицательное вещественное собственное число – $\omega^{+}$кратности два. Обозначим характеристический многочлен матрицы $\mathbf{J}^{2}$ через $Q(\omega)$. Тогда в точке бифуркации Хопфа имеют место условия

Таблица 5.13. Сходимость метода Ньютона для системы (5.1.3), (5.5.20) к точкам бифуркации Хопфа в случае задачи 2 ( $\gamma=1000, B=12$, $\beta_{1}=\beta_{2}=2, \Theta_{c 1}=\Theta_{c 2}=0, \Lambda=0,8, \mathrm{Da}_{2}=0,2, \alpha=\mathrm{Da}_{1}$ )
$Q\left(-\omega^{+}\right)=0, Q^{\prime}\left(-\omega^{+}\right)=0$. Соотношения (5.5.20) заменяются при этом соотношениями
\[
\begin{array}{l}
f_{n+1}\left(x_{1}, \ldots, x_{n}, \alpha, \omega\right)=Q(\omega)=0, \\
f_{n+2}\left(x_{1}, \ldots, x_{n}, \alpha, \omega\right)=\frac{d Q(\omega)}{d \omega}=0 .
\end{array}
\]

Систему $f_{j}(\mathbf{x}, \alpha, \omega)=0(j=1, \ldots, n+2$; см. (5.1.3) и (5.5.21) ) можно решать методом Ньютона. Если этот метод дает сходимость к точке $\left(\mathrm{x}^{+}, \alpha^{+},-\omega^{+}\right)$, где $\omega^{+}<0$, то она не является точкой бифуркации Хопфа, поскольку тогда матрица $\mathbf{J}$ имеет два вещественных собственных числа
\[
\lambda_{1,2}= \pm \sqrt{-\omega^{+}} .
\]

5.5.3. Прямые итерационные процедуры

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

Пусть $\mathbf{u}$-собственный вектор матрицы Якоби J, соответствующий собственному числу $\lambda$ :
\[
\mathbf{J} \mathbf{u}=\lambda \mathbf{u} .
\]

В точке комплексной бифуркации мы имеем $\lambda=i s, s \in \mathrm{R}^{1}$ и $\mathbf{u}=\mathbf{v}+i \mathbf{w}$, где $\mathbf{v}$ и $\mathbf{w}$ – вещественные векторы. Разделяя в (5.5.22) вещественные и мнимые части, получаем
\[
\begin{array}{l}
\mathbf{J} \mathbf{v}+s \mathbf{w}=\mathbf{0}, \\
\mathbf{J} \mathbf{w}-s \mathbf{v}=\mathbf{0} .
\end{array}
\]

Запишем теперь уравнения (5.5.23) и (5.5.24) покомпонентно в виде
\[
\begin{array}{c}
f_{i}\left(x_{1}, \ldots x_{n}, \alpha, s, v_{1}, \ldots, v_{n}, w_{1}, \ldots, w_{n}\right)=0, \\
i=n+1, n+2, \ldots, 3 n,
\end{array}
\]

где $f_{i}, i=1,2, \ldots, n,-$ правые части исходной системы (5.1.1).
Произведение вектора и на произвольное комплексное число также является собственным вектором матрицы $\mathbf{J}$, соответствую-
Рис. 5.13. Матрица размещения для уравнений (5.1.3), (5.5.23), (5.5.24).

щим собственному числу $\lambda=i s$. Это означает, что, вообще говоря, мы можем выбрать произвольно две составляющих векторов $\mathbf{v}$ и w. Мы не будем обсуждать здесь условия успешности подобного выбора. Так, могут иметь место ситуации, когда этот выбор окажется неудачным, однако при случайном выборе это маловероятно.

Система уравнений (5.1.3) и (5.5.25) представляет собой систему $3 n$ нелинейных уравнений относительно $(3 n+2)$ неизвестных $\left(x_{1}, x_{2}, \ldots, x_{n}, \alpha, s, v_{1}, v_{2}, \ldots, v_{n}, w_{1}, w_{2}, \ldots, w_{n}\right)$. Из последних $2 n$ неизвестных мы, однако, считаем две переменных фиксированными, в результате чего у нас получается система $3 n$ уравнений относительно $3 n$ неизвестных, решение которой можно вновь строить с помощью метода Ньютона. Часть элементов матрицы Якоби нетрудно найти аналитически, а оставшиеся элементы определить, аппроксимируя частные производные конечными разностями. Таблица вхождения переменных в уравнения решаемой системы представлена на рис. 5.13. Она определяет также отличные от нуля элементы матрицы Якоби (порядка $3 n$ ). Те из них, которые можно определить.

аналитически ${ }^{1)}$, обозначены на рисунке кружком. С учетом выбора двух компонент векторов $\mathbf{v}$ и $\mathbf{w}$ (которые остаются фиксированными в ходе итераций по методу Ньютона), векторы v и w могут иметь размерность в соответствии с одним из следующих трех вариантов:
1. $\operatorname{dim} \tilde{\mathbf{v}}=n-2, \operatorname{dim} \tilde{\mathbf{w}}=n, \quad$ фиксированы 2 компоненты вектора $\mathbf{v}$.
2. $\operatorname{dim} \tilde{\mathbf{v}}=n-1, \operatorname{dim} \tilde{\mathbf{w}}=n-1$, фиксированы по одной компоненте векторов $\mathbf{v}$ и $\mathbf{w}$.
3. $\operatorname{dim} \tilde{\mathbf{v}}=n, \operatorname{dim} \tilde{\mathbf{w}}=n-2$, фиксированы 2 компоненты вектора $\mathbf{w}$.

Пример сходимости метода Ньютона при решении системы уравнений (5.1.3), (5.5.23), (5.5.24) приведен в табл. 5.14. В таблице указано, какие две компоненты векторов $\mathbf{v}$ и $\mathbf{w}$ выбирались фиксированными. Предлагаем читателю сравнить найденную точку комплексной бифуркации с результатами, представленными в табл. 5.13.

Еще один метод, который мы здесь рассмотрим, основан на аналогичном принципе, однако размерность решаемой нелинейной системы оказывается более низкой ( $2 n$ вместо $3 n$ ).
Из уравнений (5.5.23) и (5.5.24) легко находим
\[
\mathbf{J}^{2} \mathbf{v}+s^{2} \mathbf{v}=\mathbf{0} .
\]

Запишем (5.5.26) в виде
\[
f_{i}\left(x_{1}, \ldots, x_{n}, \alpha, \omega, v_{1}, \ldots, v_{n}\right)=0, \quad i=n+1, \ldots, 2 n,
\]

где использовано обозначение $\omega=s^{2}$. Таким образом, мы имеем систему $2 n$ нелинейных уравнений (5.1.3), (5.5.27) относительно $2 n+2$ неизвестных $x_{1}, \ldots, x_{n}, \alpha, \omega, v_{1}, \ldots, v_{n}$. В векторе $\mathbf{v}$ мы снова можем произвольно задать две его составляющие. На рис. 5.14 представлена «таблица вхождения» неизвестных для решаемой системы. Элементы матрицы Якоби (порядка 2n), которые можно легко найти аналитически, обозначены кружками. Характер сходимости метода Ньютона хорошо виден из табл. 5.15.

В заключение параграфа сравним особенности методов, описанных в пп. 5.5.2 и 5.5.3. Так, если сравнивать ограничения на объем памяти и число операций для каждого из этих мето-

дов (приведены только главные члены в оцениваемых величинах), то из данных табл. 5.16 следует, что метод, основанный на решении системы (5.1.3), (5.5.23), (5.5.24), оказывается более предпочтительным для тех ЭВМ, у которых объем памяти достаточно велик.

Дальнейшие выводы мы можем сделать из рассмотрения данных табл. 5.17. Исходные значения всех переменных формировались случайным образом (с помощью генератора случайных чисел) из следующих промежутков: $x_{1}, x_{2} \in(0,1), \Theta_{1}$, $\Theta_{2} \in(0,4), \alpha \in(0,1), \mathrm{s} \in(0,10), v_{i}, w_{i} \in(-20,20)$, где $s^{2}=\omega$.

Рис. 5.14. Матрица размещения для уравнений (5.1.3), (5.5.26), $\operatorname{dim} \tilde{\mathbf{v}}
eq n-2$.

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

Таблица 5.16. Скорость роста объема памяти и числа операций при увеличении $n$ для разных методов отыскания точки комплексной бифуркации.

Таблица 5.17. Эффективность различных методов для задачи 2. Здесь же указано число начальных приближений (из общего числа 1000 приближений, выбранных случайным образом), при которых метод Ньютона сходится к какому-либо из трех решений. Значения параметров: $\gamma=1000, B=12$, $\beta_{1}=\beta_{2}=2, \Theta_{c 1}=\theta_{c 2}=0, \Lambda=0,8, \mathrm{Da}_{1}=0,2, \alpha=\mathrm{Da}_{1}$.

Методы, рассмотренные в данном пункте, легко обобщаются на случай нахождения точек комплексной бифуркации для дифференциальных уравнений с частными производными параболического типа. Это обобщение будет проведено в гл. 6.