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

6.3.1. Первичные бифуркации

Некоторые задачи обладают так называемыми тривиальными решениями, не зависящими от части параметров. K ним, в частности, относятся системы типа «реакция – диффузия», описанные в гл. 4. Тривиальное стационарное решение для такой двухкомпонентной системы (4.3.5), (4.3.6) однородно по пространству, т. е. имеет вид
\[
x(z)=\bar{x}, \quad y(z) \equiv \bar{y} .
\]

Оно не зависит от величины параметра $L$, задающего размеры системы, а также от величины коэффициентов диффузии $D_{x}$ и $D_{y}$.

Значения $\bar{x}$ и $\bar{y}$ определяются уравнениями $f(\bar{x}, \bar{y})=0$, $g(\bar{x}, \bar{y})=0$ и могут, естественно, зависеть от параметров, входящих в функции $f$ и $g$. При изменении этих параметров могут происходить бифуркации, не нарушающие пространственной однородности. Их можно исследовать методами, описанными в гл. 5.

В этом пункте мы рассмотрим бифуркации тривиальных решений, нарушающие пространственную однородность и связан-

ные с изменением параметра $L$-размера системы. Такие бифуркации мы будем называть первичными.

Устойчивость тривиального решения (6.3.1) определяется собственными числами ${ }^{1)}$ линеаризованного оператора [2.32]

при соответствующих (однородных) граничных условиях (типа ГУ2 или ГУ1). Здесь
\[
\begin{array}{ll}
a_{11}=\frac{\partial f(\bar{x}, \bar{y})}{\partial x}, & a_{12}=\frac{\partial f(\bar{x}, \bar{y})}{\partial y}, \\
a_{21}=\frac{\partial g(\bar{x}, \bar{y})}{\partial x}, & a_{22}=\frac{\partial g(\bar{x}, \bar{y})}{\partial y} .
\end{array}
\]

Тривиальное решение ( $\bar{x}, \bar{y}$ ) устойчиво, если все собственные числа оператора $\mathscr{L}$ имеют отрицательную вещественную часть и неустойчиво, если хотя бы одно собственное число имеет положительную вещественную часть.

Обсудим теперь устойчивость тривиального решения в зависимости от величины параметра $L$ и типа граничных условий. Рассмотрим сначала очень малые $L$. В случае ГУ1 при $L \rightarrow+0$ (точнее, при $L<L_{\min }$ ) всякое тривиальное решение будет устойчивым [6.17]. В случае же ГУ2 вопрос об устойчивости при $L \rightarrow 0$ решается в зависимости от поведения системы при отсутствии пространственных градиентов. (Более формально – от поведения системы обыкновенных дифференциальных уравнений, получающейся при $D_{\mathrm{x}}=D_{\mathrm{y}}=0$ ). Именно, если все собственные числа матрицы
\[
\mathbf{A}=\left[\begin{array}{ll}
a_{11}, & a_{12} \\
a_{21}, & a_{22}
\end{array}\right]
\]

лежат в левой (комплексной) полуплоскости, тривиальное решение будет устойчивым в случае ГУ 2 при $L \rightarrow+0$. Наоборот, наличие собственного числа матрицы А с положительной вещественной частью влечет за собой неустойчивость тривиального решения в случае ГУ 2.

С ростом $L$ какое-нибудь из собственных чисел оператора $\mathscr{L}$ может пересечь мнимую ось. При этом в результате веществен-

ной или комплексной бифуркации (бифуркации Андронова Хопфа) может появиться новое решение (нелинейной) задачи. Значение параметра $L$, соответствующее первичной бифуркации, называют первичной бифуркационной длиной. Это значение находится следующим образом.

В рассматриваемом случае (тривиального стационарного решения) оператор $\mathscr{L}$ имеет постоянные коэффициенты и его собственные функции выписываются в явном виде. А именно:
в случае ГУ1: $\mathbf{u}=\left[c_{n} \sin (n \pi z), d_{n} \sin (n \pi z)\right], n=1,2, \ldots$,
в случае ГУ2: $\mathbf{u}=\left[c_{n} \cos (n \pi z), d_{n} \cos (n \pi z)\right], n=0,1,2, \ldots$.

Подставив эти выражения в равенство $\mathscr{L} \mathbf{u}=\lambda \mathbf{u}$, мы получаем: 1) собственное число $\lambda$ оператора $\mathscr{L}$, соответствующее собственной функции (6.3.4) или (6.3.5) при определенном значении $n$, одновременно является собственным числом матрицы
\[
\mathbf{A}_{n}=\left[\begin{array}{cc}
-D_{x} E_{(n)}+a_{11}, & a_{12} \\
a_{21}, & -D_{\mathrm{y}} E_{(n)}+a_{22}
\end{array}\right], \quad E_{(n)}=\left(\frac{n \pi}{L}\right)^{2},
\]

которая имеет один и тот же вид для ГУ1 и ГУ2;2) вектор $\left(\begin{array}{l}c_{n} \\ d_{n}\end{array}\right)$ есть собственный вектор матрицы $\mathbf{A}_{n}$. Собственныечисла $\lambda$ матрицы $\mathbf{A}_{n}$ – корни характеристического уравнения
\[
\lambda^{2}-P_{n}(L) \lambda+Q_{n}(L)=0,
\]

где
\[
\begin{array}{c}
P_{n}(L)=a_{11}+a_{22}-\left(D_{x}+D_{y}\right) E_{(n)}, \\
Q_{n}(L)=D_{x} D_{y} E_{(n)}^{2}-\left(D_{x} a_{22}+D_{y} a_{11}\right) E_{(n)}+a_{11} a_{22}-a_{12} a_{21} .
\end{array}
\]

Для того, чтобы $\lambda=0$ было собственным числом оператора $\mathscr{L}$, т. е. чтобы имела место вещественная бифуркация, необходимо ${ }^{1)}$, чтобы при некотором $n$
\[
Q_{n}(L)=0 .
\]

Для точки первичной комплексной бифуркации должно выполняться условие
\[
P_{n}(L)=0, \quad Q_{n}(L)>0 .
\]

Значение $E_{(n)}^{*}$, которое удовлетворяет соотношению (6.3.10), не зависит от $n$ :
\[
\begin{aligned}
E_{(n)}^{*}=E_{1,2}^{*}=\frac{1}{2 D_{\mathrm{x}} D_{\mathrm{y}}}\left\{D_{\mathrm{x}} a_{22}+D_{\mathrm{y}} a_{11} \pm\left[\left(D_{\mathrm{x}} a_{22}+D_{\mathrm{y}} a_{11}\right)^{2}-\right.\right. \\
\left.\left.-4 D_{\mathrm{x}} D_{\mathrm{y}}\left(a_{11} a_{22}-a_{12} a_{21}\right)\right]^{1 / 2}\right\} .
\end{aligned}
\]

Мы можем получить два, одно или ни одного значения $E^{*}$, для которых имеет место вещественная бифуркация. После этого находим бифуркационную длину
\[
L_{(n) 1,2}^{*}=\frac{n \pi}{\sqrt{E_{1,2}^{*}}},
\]

если $E^{*}>0$. Из последней формулы также следует, что
\[
L_{(n)}^{*}=n L_{(1)}^{*}
\]

для обеих групп бифуркационных длин (определяемых индексами 1 и 2 в выражении (6.3.13)). Величину $L_{(1)}^{*}$ называют эле-

Рис. 6.10. Бифуркационная диаграмма первичных бифуркаций; а) общая схема, $b$ ) для задачи $11, A=2, D_{\mathrm{x}}=0,0016, D_{\mathrm{y}}=0,008$; сплошные линии $Q=0$, прерывистые $-P=0$.

ментарной бифуркационной длиной; остальные бифуркационные длины оказываются кратными элементарным бифуркационным длинам.

Для значения $L^{+}$, при котором имеет место комплексная бифуркация, из соотношений (6.3.11) находим
\[
L_{(1)}^{+}=\pi\left[\frac{D_{\mathrm{x}}+D_{\mathrm{y}}}{a_{11}+a_{22}}\right]^{1 / 2}, \quad Q_{1}\left(L_{(1)}^{+}\right)>0 .
\]

Формула (6.3.14) остается справедливой и в этом случае, т. е.
\[
L_{(n)}^{+}=n L_{(1)}^{+} .
\]

Условия существования бифуркационных длин в зависимости от величины отношения $D_{\mathrm{x}} / D_{\text {у }}$ и коэффициентов $a_{i j}$ приведены, например, в работе [5.10]. Если построить графики зависимостей величин $L^{*}$ и $L^{+}$от другого параметра задачи $\alpha$, то мы получим бифуркационную диаграмму первичных бифуркаций. Такого рода бифуркационная диаграмма представлена схематически на рис. 6.10 a. При $\alpha<\alpha_{1}$ не существует какой-либо бифуркационной длины, при $\alpha \in\left(\alpha_{1}, \alpha_{2}\right) \cup\left(\alpha_{3}, \alpha_{4}\right)$ имеются два значения $L_{(1)}^{*}$ а при $\alpha \in\left(\alpha_{2}, \alpha_{3}\right) \cup\left(\alpha_{4}, \infty\right)$ имеются два значения $L_{(1)}^{*}$ и одно значение $L_{(1)}^{+}$. На рис. 6.10 b эта диаграмма изображена для задачи 11 (показаны также значения некоторых «кратных» бифуркационных длин $\left.L_{(n)}\right)$. Читатель может сравнить точки первичных бифуркаций на рис. 6.3 и 6.4 со значениями $L$ на рис. $6.10 \mathrm{~b}$.
6.3.2. Нахождение точек вещественной бифуркации

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

6.3.2.1. Разностные методы

Рассмотрим для простоты краевую задачу (6.2.1), (6.2.2). Необходимым условием существования вещественной бифуркации является требование, чтобы линеаризованное уравнение
\[
\delta^{\prime \prime}-\frac{\partial f\left(z, y, y^{\prime}, \alpha\right)}{\partial y} \delta-\frac{\partial f\left(z, y, y^{\prime}, \alpha\right)}{\partial y^{\prime}} \delta^{\prime}=0
\]

с граничными условиями
\[
a_{0} \delta(0)+b_{0} \delta^{\prime}(0)=0, \quad a_{1} \delta(1)+b_{1} \delta^{\prime}(1)=0
\]

имело ненулевое решение. Если такое решение $\delta(z)$ существует, то $c \delta(z)$ также будет решением. Чтобы выделить одно конкретное решение, нужно некое условие нормировки. Одна из возможностей заключается в том, чтобы выбрать
\[
\delta(0)=1
\]

в случае, когда $b_{0}
eq 0$. (Если $b_{0}=0$, то можно положить $\boldsymbol{\delta}^{\prime}(0)=1$.) Уравнения (6.2.1), (6.3.17) вместе с граничными условиями (6.2.2), (6.3.18), (6.3.19) представляют собой нелинейную краевую задачу, число условий в которой на единицу превышает порядок системы. Это, однако, компенсируется тем, что параметр $\alpha$ мы рассматриваем в качестве неизвестной; решив краевую задачу, мы получим значение $\alpha$, отвечающее вещественной бифуркации. Замена данной задачи соответствующими разностными формулами приводит к системе нелинейных (алгебраических) уравнений с почти ленточной схемой размещения (ленточный характер матрицы нарущает столбец, который соответствует неизвестной $\alpha$ ). Результаты так проведенных расчетов для задачи 16 представлены в работе [6.10].

Другую возможность применения разностных методов можно продемонстрировать при нахождении точек поворота в задаче $17^{11}$. Введем следующие обозначения:
\[
f=\frac{\partial F}{\partial k}, \quad g=\frac{\partial G}{\partial k}, \quad h=\frac{\partial H}{\partial k} .
\]

Продифференцировав уравнения (Р17-16)-(P17-18) по переменной $k$, найдем
\[
\begin{array}{l}
f^{\prime \prime}=\sqrt{\operatorname{Re}}\left(h F^{\prime}+H f^{\prime}\right)+\operatorname{Re}(2 F f-2 G g+1), \\
g^{\prime \prime}=2 \operatorname{Re}(F g+f G)+\sqrt{\operatorname{Re}}\left(g^{\prime} H+G^{\prime} h\right), \\
h^{\prime}=-2 \sqrt{\operatorname{Re}} f, \\
h(0)=f(0)=h(1)=f(1)=g(0)=0, \quad g(1)=S^{\prime}(k) .
\end{array}
\]

В точке поворота (по параметру $S$ ) должно быть
\[
\frac{d S}{d k}=0, \quad \text { т. е. } g(1)=0 .
\]

Таким образом мы получаем систему шести дифференциальных уравнений (P17-16), (P17-17), (P17-18), (6.3.21a), (6.3.21b), (6.3.21c) (имеющую суммарно десятый порядок) и двенадцать граничных условий (Р17-19), (P17-20), (6.3.21d), (6.3.22). На

первый взгляд задача переопределена; на самом деле это не так, поскольку кроме $f(z), g(z), h(z)$ ищутся значения двух неизвестных $k$ и $S$.

Для решения указанной нелинейной краевой задачи использовались стандартные разностные замены на сетке с узлами $z_{i}=i h, i=0,1, \ldots, n, h=1 / n$ (аналогичные заменам вида
Рис. 6.11. Бифуркационная диаграмма задачи 17.
(6.1.23)). В частности, уравнению (6.3.21b) отвечает разностная система
\[
\begin{aligned}
\frac{g_{i-1}-2 g_{i}+g_{i+1}}{h^{2}}=2 \operatorname{Re} & \left(F_{i} g_{i}+f_{i} G_{i}\right)+ \\
& +\sqrt{\operatorname{Re}}\left[\frac{g_{i+1}-g_{i-1}}{2 h} H_{i}+\frac{G_{i+1}-G_{i-1}}{2 h} h_{i}\right] .
\end{aligned}
\]

а уравнению (6.3.21c) – соотношения вида
\[
\frac{h_{i}-h_{i-1}}{h}=-\sqrt{\operatorname{Re}}\left(f_{i}+f_{i-1}\right), \quad i=1,2, \ldots, n .
\]

Если ввести вектор неизвестных величин
$\mathbf{X}=\left(H_{0}, h_{0}, F_{n}, f_{0}, G_{0}, g_{0}, H_{1}, h_{1}, F_{1}, \ldots\right.$
\[
\left.\ldots, H_{n}, h_{n}, F_{n}, f_{n}, G_{n}, g_{n}, k, S\right)
\]

и подходящим образом упорядочить разностные уравнения, то мы получим систему $6(n-1)+2$ нелинейных уравнений
\[
\varphi(\mathbf{X})=\mathbf{0}
\]

с матрицей размещения, близкой к 15-диагональной. Для решения системы (6.3.24) можно снова воспользоваться методом

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

На рис. 6.11 приведена бифуркационная диаграмма в плоскости параметров $\operatorname{Re}$ и $S$, построенная с помощью описанного выше подхода. Эта бифуркационная диаграмма не является полной. В некоторых ее областях указано число решений задачи 17 при соответствующих значениях параметров $\operatorname{Re}$ и $S$. На рис. 6.11 при $\mathrm{Re}=625$ легко найти все шесть точек поворота, указанных на рис. 6.2 .
6.3.2.2. Метод стрельбы

Опишем кратко этот метод на примере системы «реакция диффузия» (6.1.1), (6.1.2) с граничными условиями типа 2 , т. е. с условиями (6.1.4). Выбирая $\boldsymbol{\eta}$ в (6.1.24) и используя условия (6.1.4a), мы можем решить задачу Коши для уравнений (6.1.1), (6.1.2) на промежутке от точки $z=0$ до точки $z=1$, где нам нужно удовлетворить условию $(6.1 .4 \mathrm{~b})$, т. е. соотношению (6.2.32). Если ввести матрицу Якоби
\[
\mathbf{J}(\boldsymbol{\eta}, L)=\left[\frac{\partial F_{i}}{\partial \eta_{k}}\right],
\]

то, в соответствии с п. 5.4.1, необходимым условием вещественной бифуркации будет выполнение соотношения
\[
\left.F_{3}\left(\eta_{1}, \eta_{2}, L\right)=\operatorname{det} \mathbf{J}(\boldsymbol{\eta}, L)=0{ }^{1}\right) .
\]

Преимуществом метода стрельбы является то обстоятельство, что бифуркационное условие (6.3.26) записывается в пространстве меньшей размерности. Используя теперь формулу (6.1.29), перепишем (6.3.26) в виде
\[
F_{3}\left(\eta_{1}, \eta_{2}, L\right)=p_{x 1}^{\prime}(1) p_{y 2}^{\prime}(1)-p_{x 2}^{\prime}(1) p_{y 1}^{\prime}(1)=0 .
\]
Точно так же, как и в п. 5.4.1, мы будем рассматривать уравнения (6.2.32), (6.3.27) как систему трех нелинейных уравнений относительно трех неизвестных $\eta_{1}, \eta_{2}$ и $L$ («координат» точки вещественной бифуркации). Для решения этой системы

можно воспользоваться, например, методом Ньютона. Матрица Якоби для метода Ньютона вычисляется следующим образом. Из соотношений (6.1.29), (6.2.37) найдем значения производных $\partial F_{i} / \partial \eta_{1}, \partial F_{i} / \partial \eta_{2}, \partial F_{i} / \partial L, i=1,2$. Производные для $i=3$ можно получить, либо используя разностные формулы, либо с помощью вариационных переменных второго порядка
\[
\begin{array}{ll}
p_{x i j}=\frac{\partial^{2} x}{\partial \eta_{i} \partial \eta_{j}}, & p_{x i L}=\frac{\partial^{2} x}{\partial \eta_{i} \partial L}, \\
p_{y i j}=\frac{\partial^{2} y}{\partial \eta_{i} \partial \eta_{j}}, \quad p_{y i L}=\frac{\partial^{2} y}{\partial \eta_{i} \partial L} .
\end{array}
\]

Дифференцируя уравнения в вариациях (6.1.27) по $\eta_{1}, \eta_{2}$ и $L$, можно получить дифференциальные уравнения для этих переменных и соответствующие начальные условия.

Рис. 6.12. Бифуркационная диаграмма для задачи $14 ; \gamma=20, P e_{M}=10$, $P e_{H}=5, B=15, \Theta_{c}=0$. В отдельных областях указано число решений.

Для нахождения точки ветвления (см. п. 5.4.1) потребуем, помимо выполнения соотношений (6.2.32), (6.3.27), чтобы было выполнено также уравнение
\[
\begin{array}{l}
F_{4}\left(\eta_{1}, \eta_{2}, L\right)=\operatorname{det} \overline{\mathbf{G}}(\boldsymbol{\eta}, L)=p_{x 1}^{\prime}(1) p_{y L}^{\prime}(1)- \\
-p_{x L}^{\prime}(1) p_{y 1}^{\prime}(1)=0,
\end{array}
\]

где матрица $\overline{\mathbf{G}}$ получается из матрицы (6.2.33) вычеркиванием среднего столбца. Для решения системы четырех уравнений

(6.2.32), (6.3.27), (6.3.28), так же как в п. 5.4.1, используется метод Гаусса – Ньютона. Приложение указанного подхода к задаче 11 читатель может найти в работе [5.10].

Продемонстрируем использование этого метода для нахождения точек поворота в задаче 14 в случае, когда параметром является число Da. Таким образом, нас будет интересовать нахождение $y(z), \Theta(z)$ – решения краевой задачи (6.1.35), (6.1.36). Недостающие начальные условия в точке $z=1$ выберем как в $(6.1 .37): y(1)=\eta_{1}, \Theta(1)=\eta_{2}$. Требуя выполнения краевых условий в точке $z=0$ (см. 6.1.36a), мы получаем два уравнения:
\[
\begin{array}{l}
F_{1}\left(\eta_{1}, \eta_{2}, \mathrm{Da}\right) \stackrel{\text { def }}{=} \mathrm{Pe}_{\mathrm{M}} y(0)-y^{\prime}(0)=0, \\
F_{2}\left(\eta_{1}, \eta_{2}, \mathrm{Da}\right) \stackrel{\text { def }}{=} \mathrm{Pe}_{\mathrm{H}} \Theta(0)-\Theta^{\prime}(0)=0
\end{array}
\]
(см. 6.2.41). Третье, бифуркационное уравнение здесь имеет вид (ср. с (6.2.42)):
\[
\begin{array}{r}
\boldsymbol{F}_{3}(\boldsymbol{\eta}, \mathrm{Da})=\left(\mathrm{Pe}_{\mathrm{M}} p_{y 1}(0)-p_{y 1}^{\prime}(0)\right)\left(\mathrm{Pe}_{\mathrm{H}} p_{\boldsymbol{\theta} 2}(0)-p_{\boldsymbol{\theta} 2}^{\prime}(0)\right)- \\
-\left(\mathrm{Pe}_{\mathrm{H}} p_{\boldsymbol{\theta} 1}(0)-p_{\boldsymbol{\theta} 1}^{\prime}(0)\right)\left(\mathrm{Pe}_{\mathrm{M}} p_{y 2}(0)-p_{y 2}^{\prime}(0)\right)=0 .
\end{array}
\]

Результаты, полученные при решении этих 3 уравнений методом Ньютона (с использованием уравнений второго порядка в вариациях), приведены в табл. 6.10. Читатель может сравнить координаты полученных точек поворота с данными рис. 6.8. Если продолжить теперь полученные координаты точек поворота в зависимости от какого-либо другого параметра задачи, то мы получим соответствующую бифуркационную диаграмму. Пример такого продолжения по параметру представлен на рис. 6.12 в плоскости параметров $\mathrm{Da}-\beta[6.18]^{11}$

6.3.3. Нахождение точек комплексной бифуркации (бифуркации Андронова – Хопфа)
Нахождение точек комплексной бифуркации для тривиального пространственно-однородного решения рассматривалось нами в п. 6.3.1. В этом пункте мы займемся нахождением точек комплексной бифуркации на нетривиальных ветвях стационарных решений ${ }^{2}$. Наиболее простым методом является здесь
1) Точный смысл сказанного таков. Разрешив меняться $\beta$, мы получаем систему $F_{i}\left(\eta_{1}, \eta_{2}, \mathrm{Da}, \beta\right)=0(i=1,2,3)$. Ее решение определяет кривую. в пространстве ( $\left.\eta_{1}, \eta_{2}, \mathrm{Da}, \beta\right)$, связные компоненты которой можно находить,. используя какой-либо «алгоритм продолжения». Проекция найденной кривой: на плоскость ( $\mathrm{Da}, \beta$ ) дает часть линий бифуркационной диаграммы.
2) Мы называем эти точки точками вторичной комплексной бифуркации.

Таблица 6.10. Метод Ньютона для определения точек поворота в задаче 14 $\left(\gamma=20, B=15, \beta=2, \Theta_{c}=0, \mathrm{Pe}_{\mathrm{M}}=10, \mathrm{Pe}_{\mathrm{H}}=5\right)$

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

альных уравнений, методы анализа которых рассматривались в гл. 5.

Рассмотрим, к примеру, систему типа «реакция – диффузия» (4.3.7) с граничными условиями типа 1 (см. (4.3.9)). Выберем равномерную сетку узловых точек ( $z_{i}=i h, i=0,1, \ldots, n$, $h=1 / n$ ) и обозначим приближенные значения решения в этих точках как
\[
x_{i}(t) \sim x\left(z_{i}, t\right), \quad y_{i}(t) \sim y\left(z_{i}, t\right) .
\]

Апроксимируем теперь производные по координате $z$ с помощью простейших (трехточечных) разностных формул. Мы получим систему обыкновенных дифференциальных уравнений для $i=1,2, \ldots, n-1$ :
\[
\begin{aligned}
\frac{d x_{i}}{d t} & =\frac{D_{\mathrm{x}}}{L^{2} h^{2}}\left(x_{i-1}-2 x_{i}+x_{i+1}\right)+f\left(x_{i}, y_{i}\right), \\
\frac{d y_{i}}{d t} & =\frac{D_{\mathrm{y}}}{L^{2} y^{2}}\left(y_{i-1}-2 y_{i}+y_{i+1}\right)+g\left(x_{i}, y_{i}\right) .
\end{aligned}
\]

Граничные условия (4.3.9) дают
\[
x_{0}(t) \equiv \bar{x}, \quad x_{n}(t) \equiv \bar{x}, \quad y_{0}(t) \equiv \bar{y}, \quad y_{n}(t) \equiv \bar{y} .
\]

Таким образом, мы имеем систему $2(n-1)$ обыкновенных дифференциальных уравнений, для решения которой можно использовать методы, рассмотренные в § 5.5. Система (6.3.31a,b) представляет собой некоторую аппроксимацию исходной «распределенной» системы. Погрешность аппроксимации зависит от $h$ и имеет порядок $O\left(h^{2}\right)$. При уменьшении $h$ (т. е. с увеличением точности аппроксимации) возрастают размеры системы (6.3.31) и возникают трудности с вычислением собственных чисел соответствующих матриц. Поэтому описываемый подход применим лишь для не очень точного нахождения точек комплексной бифуркации, хотя в некоторых случаях даже при сравнительно большом шаге $h$ он дает вполне удовлетворительные результаты. Дискретизацию по координате $z$ можно проводить и с помощью более точных формул, например посредством ортогональной коллокации. Читателя, который более глубоко заинтересуется применением метода прямых для нахождения точек комплексной бифуркации, мы отсылаем к оригинальным работам $[6.19,6.20,6.21,6.22]$.

6.3.3.2. Прямой итерационный метод

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

зации) позволяет находить точку комплексной бифуркации гораздо точнее. Результаты, полученные с помощью первых двух подходов, удобно использовать в качестве начального приближения для прямого итерационного метода. Обратимся вновь к системе (4.3.7) с граничными условиями типа 1, т. е. условиями $\quad x(0, t)=x(1, t) \equiv \bar{x} ; \quad y(0, t)=y(1, t) \equiv \bar{y} \quad$ (см. (4.3.9)). Обозначим через $\mathscr{L}$ линеаризованный оператор правых частей уравнений (4.3.7) и вычислим его для стационарного решения $x_{0}(z), y_{0}(z)$, т. е. для решения уравнений (6.1.1), (6.1.2). Вид оператора $\mathscr{L}$ задается выражением (6.3.2) $\left(a_{11}=\frac{\partial f}{\partial x}\left(x_{0}(z), y_{0}(z)\right)\right.$ и т. д.).

В точке комплексной бифуркации оператор $\mathscr{L}$ имеет чисто мнимое собственное число $\lambda=$ is и, следовательно, имеется ненулевое решение и уравнения
\[
\mathscr{L} \mathbf{u}=i s \mathbf{u} ; \quad \mathbf{u}=\left(\begin{array}{c}
v \\
w
\end{array}\right) .
\]

Положим $\mathbf{u}=\mathbf{u}_{\mathrm{R}}+i \mathbf{u}_{\mathrm{I}}, v=v_{\mathrm{R}}+i v_{\mathrm{I}}$ и $w=w_{\mathrm{R}}+i w_{\mathrm{I}}$. Отделяя теперь вещественную и мнимую части в соотношении (6.3.32), находим
\[
\begin{array}{c}
\mathscr{L} \mathbf{u}_{\mathrm{R}}=-s \mathbf{u}_{\mathrm{I}}, \\
\mathscr{L} \mathbf{u}_{\mathrm{I}}=s \mathbf{u}_{\mathrm{R}},
\end{array}
\]

или, более подробно,
\[
\begin{array}{l}
v_{\mathrm{R}}^{\prime \prime}+\frac{L^{2}}{D_{\mathrm{x}}}\left[\frac{\partial f}{\partial x} v_{\mathrm{R}}+\frac{\partial f}{\partial y} w_{\mathrm{R}}\right]=-s v_{\mathrm{I}}, \\
v_{\mathrm{I}}^{\prime \prime}+\frac{L^{2}}{D_{\mathrm{x}}}\left[\frac{\partial f}{\partial x} v_{\mathrm{I}}+\frac{\partial f}{\partial y} w_{\mathrm{I}}\right]=s v_{\mathrm{R}}, \\
w_{\mathrm{R}}^{\prime \prime}+\frac{L^{2}}{D_{\mathrm{y}}}\left[\frac{\partial g}{\partial x} v_{\mathrm{R}}+\frac{\partial g}{\partial y} w_{\mathrm{R}}\right]=-s w_{\mathrm{I}}, \\
w_{\mathrm{I}}^{\prime \prime}+\frac{L^{2}}{D_{\mathrm{y}}}\left[\frac{\partial g}{\partial x} v_{\mathrm{I}}+\frac{\partial g}{\partial y} w_{\mathrm{I}}\right]=s w_{\mathrm{R}} .
\end{array}
\]

Граничные условия для функций $v_{\mathrm{R}}, v_{\mathrm{I}}, w_{\mathrm{R}}$, $w_{\mathrm{I}}$ с учетом формул (4.3.9) имеют вид
\[
\begin{array}{l}
v_{R}(0)=v_{\mathrm{I}}(0)=w_{\mathrm{R}}(0)=w_{\mathrm{I}}(0)=0, \\
v_{\mathrm{R}}(1)=v_{\mathrm{I}}(1)=w_{\mathrm{R}}(1)=w_{\mathrm{I}}(1)=0 .
\end{array}
\]

Собственная функция $\mathbf{u}=(v, w)$ определена с точностью до ненулевого комплексного множителя. Это означает, что для функций $v(z)$ и $w(z)$ мы можем достаточно произвольно задать два фиксированных ненулевых значения. Двумя возмож-

ными примерами такого выбора могут служить условия
\[
v_{\mathrm{R}}\left(\frac{1}{2}\right)=1, \quad w_{\mathrm{I}}\left(\frac{1}{2}\right)=1
\]

или
\[
v_{\mathrm{R}}^{\prime}(0)=1, \quad w_{\mathrm{I}}^{\prime}(0)=1 .
\]

Читатель сам может предложить другие комбинации.
Таким образом, мы получили систему обыкновенных дифференциальных уравнений (6.1.1), (6.1.2), (6.3.34) (суммарно 12 -го порядка) с 14 -ю граничными условиями (4.3.9), (6.3.35)(6.3.37) и с двумя параметрами $s$, $L$, которые нам предстоит определить в точке комплексной бифуркации. Число уравнений «согласовано» с числом неизвестных, и описанную задачу можно пытаться решать либо с помощью разностных методов, либо методом стрельбы. Рассмотрим здесь оба этих подхода.

Разностный подход. Заменим вторые производные простейшими трехточечными центрально-разностными соотношениями. В результате вместо уравнений (6.1.1), (6.1.2) и (6.3.34) мы получим следующую систему (ср. с (6.3.31)):
\[
\begin{array}{c}
x_{i-1}-2 x_{i}+x_{i+1}+\frac{L^{2} h^{2}}{D_{\mathrm{x}}} f_{i}=0, \quad f_{i}=f\left(x_{i}, y_{i}\right), \\
y_{i-1}-2 y_{i}+y_{i+1}+\frac{L^{2} h^{2}}{D_{\mathrm{y}}} g_{i}=0, \quad g_{i}=g\left(x_{i}, y_{i}\right), \\
v_{\mathrm{R}, i-1}-2 v_{\mathrm{R}, i}+v_{\mathrm{R}, i+1}+\frac{L^{2} h^{2}}{D_{\mathrm{x}}}\left(f_{x, i} v_{\mathrm{R}, i}+f_{y, i} w_{\mathrm{R}, i}\right)=-s v_{\mathrm{I}, i}, \\
w_{1, i-1}-2 w_{\mathrm{I}, i}+w_{\mathrm{I}, i+1}+\frac{L^{2} h^{2}}{D_{\mathrm{y}}}\left(g_{x, i} v_{\mathrm{I}, i}+g_{y, i} w_{\mathrm{I}, i}\right)=s w_{\mathrm{R}, i},
\end{array}
\]

где $i=1, \ldots, n-1$.
Дополним эту систему очевидными граничными условиями:
\[
\begin{array}{cl}
x_{0}=x_{n}=\bar{x}, & y_{0}=y_{n}=\bar{y}, \quad v_{\mathrm{R}, 0}=v_{\mathrm{I}, 0}=w_{\mathrm{R}, 0}=w_{\mathrm{I}, 0}=0 ; \\
v_{\mathrm{R}, n}=v_{\mathrm{I}, n}=w_{\mathrm{R}, n}=w_{\mathrm{I}, n}=0 .
\end{array}
\]

Тем самым мы получим систему $6(n+1)$ нелинейных (алгебраических) уравнений относительно $6(n+1)+2$ неизвестных
\[
x_{i}, y_{i}, v_{\mathrm{R}, i}, v_{\mathrm{I}, i}, w_{\mathrm{R}, i} w_{\mathrm{I}, i}, \quad(i=0, \ldots, n), s, L .
\]

В качестве двух недостающих условий можно, например, взять соотношения (6.3.37a). Полученную в результате систему $6(n+1)+2$ нелинейных алгебраических уравнений можно решать с помощью метода Ньютона.

В п. 6.3.3.1 мы рассмотрели аппроксимацию исходных дифференциальных уравнений в частных производных системой

обыкновенных дифференциальных уравнений с помощью метода прямых. Точку комплексной бифуркации для этой аппроксимирующей системы (6.3.31) можно найти с помощью одного из методов, описанных в параграфе 5.5. При этом, если воспользоваться для решения указанной системы прямым итерационным методом в его несокращенной версии, которая была описана в п. 5.5.3, то нетрудно показать, что получающаяся нелинейная (алгебраическая) система имеет тот же самый вид, что и система (6.3.38), если подходящим образом выбрать для нее граничные условия ${ }^{1)}$. Это утверждение читатель может легко проверить самостоятельно.

Метод стрельбы. Этот подход мы продемонстрируем на примере задачи 14, т. е. уравнений (Р14-7), (Р14-8) с граничными условиями (P14-9), (P14-10). Правые части этих уравнений отличаются от правых частей уравнений (4.3.7) членами с первой производной; кроме того, имеются различия и в граничных условиях (представляющих собой комбинацию граничных условий 1 -го и 2 -го рода). В качестве бифуркационного параметра выберем параметр $\mathrm{Da}$ и положим $\mathrm{Le}=1$.

Обозначим $E=(1-y) \exp (\theta /(1+\Theta / \gamma))$ и запишем уравнения, которым должны удовлетворять стационарные решения
\[
\begin{array}{c}
\frac{1}{\mathrm{Pe}_{M}} y^{\prime \prime}-y^{\prime}+\mathrm{Da} E=0, \\
\frac{1}{\mathrm{Pe}_{\mathrm{H}}} \Theta^{\prime \prime}-\Theta^{\prime}+B \mathrm{Da} E-\beta\left(\Theta-\Theta_{c}\right)=0 .
\end{array}
\]

Граничные условия при этом имеют вид
\[
\begin{array}{l}
y^{\prime}(0)-\mathrm{Pe}_{\mathrm{M}} y(0)=0, \quad \Theta^{\prime}(0)-\mathrm{P}_{\mathrm{H}} \Theta(0)=0, \\
y^{\prime}(1)=0, \quad \Theta^{\prime}(1)=0 . \\
\end{array}
\]

Уравнения (6.3.33) для этого случая (в компонентах $v, w$ ) записываются так:
\[
\begin{array}{c}
\frac{1}{\mathrm{Pe}_{\mathrm{M}}} v_{\mathrm{R}}^{\prime \prime}-v_{\mathrm{R}}^{\prime}+\mathrm{Da}\left[\frac{\partial E}{\partial y} v_{\mathrm{R}}+\frac{\partial E}{\partial \Theta} w_{\mathrm{R}}\right]=-s v_{\mathrm{I}}, \\
\frac{1}{\mathrm{Pe}_{\mathrm{M}}} v_{\mathrm{I}}^{\prime \prime}-v_{\mathrm{I}}^{\prime}+\mathrm{Da}\left[\frac{\partial E}{\partial y} v_{\mathrm{I}}+\frac{\partial E}{\partial \Theta} w_{\mathrm{I}}\right]=s v_{\mathrm{R}}, \\
\frac{1}{\mathrm{Pe}_{\mathrm{H}}} w_{\mathrm{R}}^{\prime \prime}-w_{\mathrm{R}}^{\prime}+B \mathrm{Da}\left[\frac{\partial E}{\partial y} v_{\mathrm{R}}+\frac{\partial E}{\partial \Theta} w_{\mathrm{R}}\right]-\beta w_{\mathrm{R}}=-s w_{\mathrm{I}}, \\
\frac{1}{\mathrm{Pe}_{\mathrm{H}}} w_{\mathrm{I}}^{\prime \prime}-w_{\mathrm{I}}^{\prime}+B \mathrm{Da}\left[\frac{\partial E}{\partial y} v_{\mathrm{I}}+\frac{\partial E}{\partial \Theta} w_{\mathrm{I}}\right]-\beta w_{\mathrm{I}}=s w_{\mathrm{R}} .
\end{array}
\]
1) Имеются в виду дополнительные уравнения типа (6.3.37).- Прим. ред.

Граничные условия для функций $v$ и ш получаются из условий (6.3.40) в виде
\[
\begin{array}{cc}
v_{\mathrm{R}}^{\prime}(0)-\mathrm{Pe}_{\mathrm{M}} v_{\mathrm{R}}(0)=0, & v_{\mathrm{I}}^{\prime}(0)-\mathrm{Pe}_{\mathrm{M}} v_{\mathrm{I}}^{\prime}(0)=0, \\
w_{\mathrm{R}}^{\prime}(0)-\mathrm{Pe}_{\mathrm{H}} w_{\mathrm{R}}(0)=0, & w_{\mathrm{I}}^{\prime}(0)-\mathrm{Pe}_{\mathrm{H}} w_{\mathrm{I}}^{\prime}(0)=0, \\
v_{\mathrm{R}}^{\prime}(1)=0, & v_{\mathrm{I}}^{\prime}(1)=0, \\
w_{\mathrm{R}}^{\prime}(1)=0, & w_{\mathrm{I}}^{\prime}(1)=0 .
\end{array}
\]

Для решения системы шести дифференциальных уравнений второго порядка (6.3.39), (6.3.41) воспользуемся методом стрельбы. Недостающие начальные условия выберем в точке $z=1$ (это связано с тем, что соответствующие задачи Коши неустойчивы, если двигаться от точки $z=0$ к точке $z=1$, см. п. 6.1.2):
\[
\begin{array}{cl}
y(1)=\eta_{1}, & \Theta(1)=\eta_{2}, \quad v_{\mathrm{R}}(1)=\eta_{3}, \quad v_{\mathrm{I}}(1)=\eta_{4}, \\
w_{\mathrm{R}}(1)=\eta_{5}, \quad w_{\mathrm{I}}(1)=\eta_{6} .
\end{array}
\]

Вместе с условиями (6.3.40b) и (6.3.42b) мы получаем 12 начальных условий, необходимых для интегрирования уравнений (6.3.39) и (6.3.41). Поскольку собственная функция
\[
\mathbf{u}=\mathbf{u}_{\mathrm{R}}+i \mathbf{u}_{\mathrm{I}}=\left(v_{\mathrm{R}}, w_{\mathrm{R}}\right)+i\left(v_{\mathrm{I}}, w_{\mathrm{I}}\right)
\]

определена с точностью до ненулевого комплексного множителя, для функций $\mathbf{u}_{\mathrm{R}}$ и $\mathbf{u}_{\mathrm{I}}$ можно зафиксировать два каких-либо конкретных значения; например, можно задать $v_{\mathrm{R}}(1)=v_{\mathrm{I}}(1)=1$ (т. е. положить в (6.3.43) $\eta_{3}=\eta_{4}=1$ ). В точке $z=0$ мы получаем шесть нелинейных уравнений (6.3.40a) и (6.3.42а) относительно шести неизвестных $\eta_{1}, \eta_{2}, \eta_{5}, \eta_{6}, s$ и Dа. Эту систему можно решать методом Ньютона, причем матрица Якоби вычисляется с помощью уравнений в вариациях, либо исходя из соответствующих разностных формул. Пример сходимости метода Ньютона к точке комплексной бифуркации представлен в табл. 6.11. При этом соответствующие итерации сходятся к точке, которая лежит на ветви решений с высокой выходной степенью конверсии $y$ (1) (ср. качественно близкую диаграмму решений на рис. 6.8 для аналогичных значений параметров). В указанной точке стационарное решение (при убывании Da) теряет устойчивость, и от него отходит ветвь устойчивых периодических решений (см. п. 6.5.1). Описанный здесь прямой итерационный метод был рассмотрен на примере двух уравнений (4.3.7) (или, соответственно, (P14-7), (Р14-8)). Метод легко обобщается на задачи с бо́льшим числом уравнений, но размерность возникающих при этом (конечномерных) задач, очевидно, возрастает.

Таблица 6.11. Пример сходимости метода Ньютона к точке комплексной бифуркации для задачи $14\left(\gamma \rightarrow \infty, B=12, \beta=2, \Theta_{c}=0, \mathrm{Pe}_{\mathrm{H}}=2\right.$, $\mathrm{Pe}_{\mathrm{M}}=2, \mathrm{Le}=1$ ).