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

5.4.1. Нахождение точек поворота и точек ветвления

Пусть снова уравнения
\[
f_{i}\left(x_{1}, x_{2}, \ldots, x_{n} ; \alpha\right)=0, \quad i=1,2, \ldots, n,
\]
1) В случае линейной системы (см. п. 5.3.1) выводы были справедливы для траекторий, начинающихся в любой точке пространства $\mathrm{R}^{n}$.

определяют стационарные решения системы дифференциальных уравнений. Для того, чтобы существовали точки (вещественной) бифуркации стационарных решений, необходимо, чтобы (где-то в изучаемой области ( $\mathbf{x}, \alpha)$ ) выполнялось условие
\[
f_{n+1}\left(x_{1}, x_{2}, \ldots, x_{n}, \boldsymbol{\alpha}\right)=\operatorname{det} \mathbf{J}\left(x_{1}, x_{2}, \ldots, x_{n}, \boldsymbol{\alpha}\right)=0,
\]

тде $\mathbf{J}$-матрица Якоби; $\mathbf{J}=\left[\partial f_{i} / \partial x_{j}\right]$. Равенство (5.4.1) есть отрицание основного условия теоремы о неявных функциях. Оно формально эквивалентно утверждению о том, что у матрицы $\mathbf{J}$ имеется собственное число, равное нулю.

Для того, чтобы в изучаемой области ( $\mathbf{x}, \alpha)$ существовали точки ветвления, необходимо, чтобы дополнительно к (5.4.1) при некотором $k$ выполнялось условие
\[
f_{n+2}\left(x_{1}, x_{2}, \ldots, x_{n}, \alpha\right)=\operatorname{det} \mathrm{J}_{k}\left(x_{1}, x_{2}, \ldots, x_{n}, \alpha\right)=0,
\]

тде $\mathbf{J}_{k}$ – матрица, определяемая формулой (5.2.12) и получаемая из матрицы $\mathbf{J}$ вычеркиванием $k$-го столбца, [ $\left.\partial f_{i} / \partial f_{n+1}\right]=$ $=\left[\partial f_{i} / \partial \alpha\right]$. В большинстве практических задач для точки поворота $\left(\mathbf{x}^{*}, \alpha^{*}\right)$ при всех $k$ выполняется условие
\[
f_{n+2}\left(\mathbf{x}^{*}, \alpha^{*}\right)
eq 0 .
\]

Наоборот, для точки ветвления ( $\mathbf{x}^{* *}, \alpha^{* *}$ ) из теоремы о неявных функциях следует, что для всех $k$
\[
f_{n+2}\left(\mathbf{x}^{* *}, \alpha^{* *}\right)=0 .
\]

Таким образом, соотношения (5.4.3) и (5.4.4) можно использовать для установления различий между точкой поворота и точкой ветвления.

Остановимся теперь на определении так называемых точек первичной бифуркации, т. е. точек бифуркации на тривиальном решении ${ }^{1)}: \mathbf{x}(\alpha) \equiv \overline{\mathbf{x}}$ для всех значений параметра $\alpha$. При этом уравнения системы (5.4.0) выполняются автоматически, а уравнение (5.4.1) определяет собой соотношение для определения значения параметра $\alpha^{* *}$ в точке бифуркации. Решение этого (одного) нелинейного уравнения может быть найдено, например, с помощью метода Ньютона.

В общем случае уравнения (5.4.0) и (5.4.1) образуют вместе систему из $n+1$ уравнений относительно $(n+1)$ неизвестных координат $x_{1}^{*}, x_{2}^{*}, \ldots, x_{n}^{*}, \alpha^{*}$ критической точки. Для нахождения решения этой системы мы вновь можем воспользоваться

методом Ньютона:
\[
\begin{array}{c}
\mathbf{F}^{\prime}\left(\mathbf{X}^{k}\right) \Delta \mathbf{X}^{k}=-\mathbf{F}\left(\mathbf{X}^{k}\right), \\
\mathbf{X}^{k+1}=\mathbf{X}^{k}+\Delta \mathbf{X}^{k} .
\end{array}
\]

Здесь использованы обозначения $\mathbf{X}=\left(x_{1}, x_{2}, \ldots, x_{n}, \boldsymbol{\alpha}\right), \mathbf{F}=$ $=\left(f_{1}, f_{2}, \ldots, f_{n+1}\right)$, причем $\mathbf{F}^{\prime}(\mathbf{X})=\left[\partial f_{i} / \partial x_{j}\right]$-это матрица Якоби. Последняя строка этой матрицы подсчитывается обычно с помощью соответствующих разностных формул. Формулы для их аналитического вычисления (полезные при небольшом $n$ ). представлены в работах $[5.9,5.10]$.

Рассмотрим теперь точку ( $x, \alpha$ ), которая удовлетворяет соотношениям (5.1.3), (5.4.1) и (5.4.2), т. е. условиям для точки ветвления. Здесь мы имеем $n+2$ уравнений для $n+1$ неизвестных составляющих вектора $\mathbf{X}$. Положим теперь $\overline{\mathbf{F}}=\left(f_{1}\right.$, $f_{2}, \ldots, f_{n+2}$ ). Для решения переопределенной системы
\[
\overline{\mathbf{F}}(\mathbf{X})=0
\]

воспользуемся методом Гаусса-Ньютона (см. [5.1]):
\[
\left[\overline{\mathbf{F}}^{\prime T}\left(\mathbf{X}^{k}\right) \overline{\mathbf{F}}^{\prime}\left(\mathbf{X}^{k}\right)\right] \Delta \mathbf{X}^{k}=-\overline{\mathbf{F}}^{\prime T}\left(\mathbf{X}^{k}\right) \overline{\mathbf{F}}\left(\mathbf{X}^{k}\right) .
\]

Қаждое новое приближение мы подсчитываем вновь по формуле (5.4.6). Указанный метод сходится к стационарной точке (минимуму) функции $\varphi=\overline{\mathbf{F}}^{T} \overline{\mathbf{F}}$; при этом в точке ветвления функция $\varphi$ имеет минимум, равный нулю (и только в этом случае метод дает нужное нам решение переопределенной системы). Если матрица $\overline{\mathbf{F}}^{\prime}\left(\mathbf{X}^{* *}\right) \overline{\mathbf{F}}^{\prime}\left(\mathbf{X}^{* *}\right)$ оказывается невырожденной, то сходимость метода, даваемого формулой (5.4.8), будет квадратической.

Наиболее часто встречающиеся ситуации представлены в табл. 5.5. Конечно, можно найти и соответствующие контрпримеры, однако результаты табл. 5.5, судя по всему, будут верны для большинства практически важных случаев.
Таблица 5.5. Наиболее вероятные ситуации
Рассмотрим простой пример для $n=1$. Пусть
\[
f_{1}=x\left(x^{2}-\alpha\right)\left[(x-1)^{2}+\alpha-4\right](1-k+k(x-\alpha))=0 .
\]

Диаграмма решений для этого случая изображена на рис. 5.7. Параметр $k$ может принимать два значения $k=0$ и $k=1$ (для
$k=1$ мы имеем двойную точку бифуркации № 5). Конкретная форма уравнений (5.4.1) и (5.4.2) в данном случае очевидна:
\[
\begin{aligned}
f_{2}= & 6 k x^{5}+5(1-3 k-k \alpha) x^{4}+4(2 k \alpha-k-2) x^{3}+ \\
& +3(5 k \alpha-3+3 k) x^{2}+2\left(-3 k \alpha^{2}+\alpha(k+2)\right) x+ \\
& +k \alpha^{3}-(2 k+1) \alpha^{2}+3(1-k) \alpha \\
f_{3}= & -k x^{5}+2 k x^{4}+5 k x^{3}+(k+2-6 k \alpha) x^{2}+ \\
& +\left(3 k \alpha^{2}-2(2 k+1) \alpha+3(1-k)\right) x .
\end{aligned}
\]

В табл. 5.6 указано число итераций, необходимых для обеспечения неравенства
\[
\max \left(\left|x^{k}-x^{*}\right|,\left|\alpha^{k}-\alpha^{*}\right|\right) \leqslant \varepsilon
\]

Таблица 5.6. Исследование примера (5.4.9). Показано число итераций, необходимое для выполнения неравенства (5.4.1).
1) Точки, отмеченные цифрами на рис. 5.7. ${ }^{2)}$ Стационарная точка $\varphi$.

при различных начальных приближениях. Из таблицы видно, что итерационный процесс (5.4.5), (5.4.6) сходится квадратично $к$ точке поворота и линейно (если он вообще сходится) к точке ветвления ${ }^{1)}$. Итерационный процесс (5.4.8) сходится квадратично к точке ветвления. Точка № 5 на рис. 5.7 при $k=1$ представляет собой точку ветвления (в ней пересекаются три ветви решения). Характер обоих итерационных процессов можно проследить в нижней части табл. 5.6. Отметим, что процесс (5.4.8) сходится к упомянутой точке лишь линейно.

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

$\operatorname{det} \mathbf{J}=0$ можно записать иначе: должен существовать не равный нулю вектор $\mathbf{v}$, такой, что
\[
\mathbf{J}(\mathbf{x}, \boldsymbol{\alpha}) \mathbf{v}=0 .
\]

Этот вектор определен с точностью до постоянного множителя, и поэтому к соотношению (5.4.12) добавляется еще условие «нормировки» вектора v, например,
\[
\sum_{i=1}^{n} v_{i}^{2}=1
\]

или
\[
v_{k}=1
\]

при заданном $k$. Уравнения (5.4.0), (5.4.12) и (5.4.13) (или (5.4.14)) образуют систему $2 n+1$ нелинейных уравнений относительно $2 n+1$ неизвестных составляющих $\mathbf{x}, \mathbf{v}, \alpha$. Решение этой системы можно искать, например, с помощью метода Ньютона. Построение матрицы Якоби для применения метода Ньютона требует вычисления вторых производных функций $f_{i}$ или использования соответствующих разностных формул. Заметим, что в случае $n=1$ метод, описываемый формулами (5.4.12), (5.4.14), оказывается тождественным подходу (5.4.5), поскольку условие для функции $f_{n+1}$ (5.4.1) идентично соотношению (5.4.12).

Другой метод нахождения точек поворота основан на использовании условия
\[
\frac{\partial \alpha}{\partial x_{k}}=0
\]

при соответствующим образом выбранном $k$. Условие (5.4.15) вместе с системой (5.4.0) вновь образует систему $n+1$ нелинейных уравнений относительно неизвестных $x_{1}, x_{2}, \ldots, x_{n}, \alpha$. Левую часть формулы (5.4.15) можно найти как последнюю составляющую вектора v, т. е. решая систему линейных алгебраических уравнений
\[
\mathbf{J}_{k}(\mathbf{x}, \alpha) \cdot \mathbf{v}=-\frac{\partial \mathrm{f}}{\partial x_{k}},
\]

где матрица $\mathbf{J}_{k}$ определяется формулой (5.2.12). При $n=1$ условие (5.4.15) вновь оказывается эквивалентным соотношению (5.4.5).

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

одна переменная $\Theta$ и один «живой» параметр Da.-Peд.) В п. 5.2.1 мы вывели соотношение (5.2.4), из которого находим
\[
\begin{array}{c}
\mathrm{Da}=\left[B-\Theta-\frac{\beta\left(\Theta-\Theta_{c}\right)}{\Lambda}\right]^{-1}\left[\Lambda \Theta+\beta\left(\Theta-\Theta_{c}\right)\right] \cdot \exp \frac{-\Theta}{1+\Theta / \gamma}= \\
=g\left(\Theta, \Lambda, \gamma, B, \beta, \Theta_{c}\right) .
\end{array}
\]

В точке поворота имеет место условие $d \mathrm{Da} / d \Theta=0$, или
\[
\frac{\partial g\left(\boldsymbol{\theta}, \boldsymbol{\Lambda}, \boldsymbol{\gamma}, B, \boldsymbol{\beta}, \Theta_{c}\right)}{\partial \Theta}=0 .
\]

Если мы найдем корень $\Theta^{*}$ уравнения (5.4.17) при фиксированных значениях параметров $\Lambda, \gamma, B, \beta$ и $\Theta_{c}$, то из формулы (5.4.16) получим значение $\mathrm{Da}^{*}$ в точке поворота. Проведя дифференцирование в (5.4.17), после некоторых преобразований мы получаем квадратное уравнение относительно неизвестной $\Theta$. Тем самым можно найти максимум два положительных корня $\Theta$, которые соответствуют точкам поворота на диаграмме решений – зависимости от Da (см. также данные табл. 5.3).

5.4.2. Продолжение решения из точки ветвления

В предыдущем пункте мы рассмотрели некоторые методы нахождения точек ветвления на диаграмме решений. В $\S 3.4$ были получены соотношения для определения касательных векторов ветвей решения, отходящих из данной точки ветвления. В то время как процедуры для определения точек ветвления носили итерационный характер, угловые коэффициенты ветвей разыскивались с помощью различных финитных подходов, основанных на методе исключения Гаусса (точка ветвления предполагалась известной заранее). В настоящем пункте мы опишем, как используются эти сведения для процесса продолжения по параметру, т. е. построения диаграммы стационарных решений (см. [5.12]).

При этом мы будем предполагать, что имеет место случай «общего положения»: ранг матрицы Якоби $\mathbf{J}$ и расширенной матрицы $\mathbf{J}$ понижается только на 1:
\[
\operatorname{rank} \mathbf{J}=\operatorname{rank} \tilde{\mathbf{J}}=n-1 .
\]

Указанный алгоритм состоит из 9 шагов:

1. Нахождение точки ветвления ( $\left.x^{* *}, \alpha^{* *}\right)$, см. п. 5.4.1.
2. Вычисление частных производных функций $f_{i}$ в этой точке, построение матрицы $\tilde{\mathbf{J}}$ по формуле (3.3.3).
3. Решение двух систем линейных алгебраических уравнений (3.4.2) и (3.4.3) с целью нахождения частных производных
\[
\frac{\partial \varphi_{i}}{\partial x_{1}}, \quad \frac{\partial \varphi_{i}}{\partial \alpha}, \quad i=2,3, \ldots, n .
\]
Это решение можно получить, применяя метод исключения Гаусса для двух различных правых частей одновременно, поскольку системы (3.4.2) и (3.4.3) имеют одну и ту же матрицу.
4. Нам потребуются частные производные второго порядка от функции $F\left(x_{1}, \alpha\right)$, даваемой формулой (3.3.11a). С этой целью вычислим производные функции $f_{i}$ в точке ( $\left.x^{* *}, \alpha^{*}\right)$ :
\[
f_{i j}^{\prime}=\frac{\partial f_{i}}{\partial x_{j}}, \quad f_{i \alpha}^{\prime}=\frac{\partial f_{i}}{\partial \alpha}, \quad f_{i j k}^{\prime}=\frac{\partial^{2} f_{i}}{\partial x_{j} \partial x_{k}}, \quad f_{i j \alpha}^{\prime}=\frac{\partial^{2} f_{i}}{\partial x_{j} \partial \alpha} .
\]
Значения частных производных первого порядка мы уже нашли на шаге 2. Обозначим аналогично производные функций $\varphi_{i}$ :
\[
\begin{array}{c}
\varphi_{i 1}^{\prime}=\frac{\partial \varphi_{i}}{\partial x_{1}}, \quad \varphi_{i \alpha}^{\prime}=\frac{\partial \varphi_{i}}{\partial \alpha}, \quad \varphi_{i 11}^{\prime}=\frac{\partial^{2} \varphi_{i}}{\partial x_{1}^{2}}, \\
\varphi_{i 1 \alpha}^{\prime}=\frac{\partial^{2} \varphi_{i}}{\partial x_{1} \partial \alpha}, \quad \varphi_{i \alpha \alpha}^{\prime}=\frac{\partial^{2} \varphi_{i}}{\partial \alpha^{2}} .
\end{array}
\]
В этих обозначениях системы (3.4.2) и (3.4.3) переписываются в виде
\[
\begin{array}{c}
\sum_{j=2}^{n} f_{i j}^{\prime} \varphi_{j 1}^{\prime}=-f_{i 1}^{\prime}, \quad i=1,2, \ldots, n-1, \\
\sum_{j=2}^{n} f_{i j}^{\prime} \varphi_{j \alpha}^{\prime}=-f_{i \alpha}^{\prime}, \quad i=1,2, \ldots, n-1 .
\end{array}
\]
Дифференцируя соотношение (3.4.2a) по $x_{1}$, находим
\[
\begin{array}{c}
\sum_{j=2}^{n} f_{i j}^{\prime} \varphi_{j 11}^{\prime}=-\sum_{j=2}^{n}\left[2 f_{i j 1}^{\prime} \varphi_{j 1}^{\prime}+\sum_{k=2}^{n}\left(f_{i j k}^{\prime} \varphi_{k 1}^{\prime}\right) \varphi_{j 1}^{\prime}\right]-f_{i 11}^{\prime}, \\
i=1,2, \ldots, n-1 .
\end{array}
\]
Аналогичным образом, дифференцируя формулу (3.4.2a) по $\alpha$, получаем
\[
\begin{array}{c}
\sum_{j=2}^{n} f_{i j}^{\prime} \varphi_{j 1 \alpha}^{\prime}=-\sum_{j=2}^{n}\left[f_{i j \alpha}^{\prime} \varphi_{j 1}^{\prime}+f_{i j 1}^{\prime} \varphi_{j \alpha}^{\prime}+\sum_{k=2}^{n}\left(f_{i j k}^{\prime} \varphi_{k \alpha}^{\prime}\right) \varphi_{j 1}^{\prime}\right]-f_{i 1 \alpha}^{\prime}, \\
i=1,2, \ldots, n-1 .
\end{array}
\]

Наконец, дифференцирование (3.4.3а) по $\alpha$ дает
\[
\begin{array}{c}
\sum_{j=2}^{n} f_{i j}^{\prime} \varphi_{i \alpha \alpha}^{\prime}=-\sum_{j=2}^{n}\left[2 f_{i j \alpha}^{\prime} \varphi_{j \alpha}^{\prime}+\sum_{k=2}^{n}\left(f_{i j k}^{\prime} \varphi_{k \alpha}^{\prime}\right) \varphi_{j \alpha}^{\prime}\right]-f_{i \alpha \alpha}^{\prime}, \\
i=1,2, \ldots, n-1 .
\end{array}
\]
Соотношения (5.4.21), (5.4.22) и (5.4.23) представляют собой 3 системы линейных алгебраических уравнений относительно неизвестных $\varphi_{i 11}^{\prime}, \varphi_{i 1 \alpha}^{\prime}, \varphi_{i \alpha a}^{\prime}$, обладающие одинаковой матрицей. Тем самым, с помощью метода исключения Гаусса их можно решать одновременно. (Правые части систем составлены из величин, численные значения которых уже определены на предыдущих шагах алгоритма.)
5. Находим частные производные второго порядка для функции $F$, задаваемой формулой (3.3.11a):
\[
\begin{aligned}
A= & \frac{\partial^{2} F}{\partial x_{1}^{2}}=f_{n 11}^{\prime}+\sum_{j=2}^{n}\left[2 f_{n 1 j}^{\prime} \varphi_{j 1}^{\prime}+f_{n j}^{\prime} \varphi_{j 11}^{\prime}+\sum_{k=2}^{n}\left(f_{n j k}^{\prime} \varphi_{k 1}^{\prime}\right) \varphi_{j 1}^{\prime}\right] \\
B= & \frac{\partial^{2} F}{\partial x_{1} \partial \alpha}=f_{n 1 \alpha}^{\prime}+\sum_{j=2}^{n}\left[f_{n 1 j}^{\prime} \varphi_{j \alpha}^{\prime}+f_{n j \alpha}^{\prime} \varphi_{j 1}^{\prime}+f_{n j}^{\prime} \varphi_{j 1 \alpha}^{\prime}+\right. \\
& \left.+\sum_{k=2}^{n}\left(f_{n j k}^{\prime} \varphi_{k \alpha}^{\prime}\right) \varphi_{j 1}^{\prime}\right] \\
C= & \frac{\partial^{2} F}{\partial \alpha^{2}}=f_{n \alpha \alpha}^{\prime}+\sum_{j=2}^{n}\left[2 f_{n j \alpha}^{\prime} \varphi_{j \alpha}^{\prime}+f_{n j}^{\prime} \varphi_{j \alpha \alpha}^{\prime}+\sum_{k=2}^{n}\left(f_{n j k}^{\prime} \varphi_{k \alpha}^{\prime}\right) \varphi_{j \alpha}^{\prime}\right] .
\end{aligned}
\]
6. Строим стартовые точки для алгоритма продолжения, например, алгоритма DERPAR (см. п. 5.2.3). Выбираем две малые числовые постоянные
\[
h_{1}>0, \quad h_{\alpha}>0 .
\]
7. По крайней мере одно из двух искомых направлений ответвления характеризуется конечным значением производной $d x_{1} / d \alpha$ (см. (3.2.4)). Если оба значения производной оказываются конечными, то, выбирая меньшее из этих двух значений, полагаем
направление 1:
\[
\frac{d x_{1}}{d \alpha}=S_{1},
\]
где $S_{1}$ вычисляется согласно формуле (3.2.6). Другое направление характеризуется конечным значением производной $d \alpha / d x_{1}$ и, следовательно,

направление 2:
\[
\frac{d \alpha}{d x_{1}}=S_{2} .
\]
8. С помощью полученных таким образом двух характеристических направлений мы находим четыре исходных точки для

Таблица 5.7. Стартовые точки для процесса продолжения. Для каждого направления мы получаем две точки, одну для $P=1$ и другую для $P=-1$. Здесь $N_{i}$ – параметры направления для процесса продолжения, $k$ – индекс переменной, которая остается фиксированной в ходе ньютоновских итераций для уточнения стартовой точки.

процесса продолжения решения (см. табл. 5.7). Здесь использовано то, что изменение остальных переменных подчиняется соотношениям
\[
\frac{d x_{i}}{d \alpha}=\varphi_{i 1}^{\prime} \frac{d x_{1}}{d \alpha}+\varphi_{i \alpha}^{\prime}
\]
или
\[
\frac{d x_{i}}{d x_{1}}=\varphi_{i 1}^{\prime}+\varphi_{i \alpha}^{\prime} \frac{d \alpha}{d x_{1}} .
\]
В таблице приведены также направляющие параметры $N_{i}$ для каждой из переменных (см. алгоритм в п. 5.2.3) и индекс $k$ переменной, которая не изменяется в ходе итераций метода. Ньютона для уточнения стартовой точки.
9. По окончании проведения ньютоновских итераций мы можем проконтролировать, соответствуют ли найденные стартовые

точки ( $\left.x_{1}, \ldots, x_{n}, \alpha\right)$ предсказанным точкам; например, для направления 1 это можно сделать с помощью оценок
\[
\begin{array}{l}
\left|\frac{x_{i}-x_{i}^{* *}}{\alpha-\alpha^{* *}}-\varphi_{i 1}^{\prime} S_{1}-\varphi_{i \alpha}^{\prime}\right|<\varepsilon, \quad i=2, \ldots, n, \\
\left|\frac{x_{1}-x_{1}^{* *}}{\alpha-\alpha^{* *}}-S_{1}\right|<\varepsilon .
\end{array}
\]

Для направления 2 тестовая оценка производится аналогично.

Рис. 5.8. Диаграмма решений задачи (5.4.31) в окрестности точки бифуркаций-
Проиллюстрируем применение описанного алгоритма на следующем примере:
\[
\begin{array}{l}
f_{1}=x_{1}-2 x_{2}-\alpha-x_{1} x_{2}-3 x_{2} \alpha+x_{1}^{2}-2 x_{2}^{2}-\alpha^{2}=0, \\
f_{2}=x_{2}^{2}-x_{2}^{2} x_{1}^{2}+\alpha^{2} x_{1}^{2}-\alpha^{2}=0 .
\end{array}
\]

Сходимость метода Гаусса-Ньютона (5.4.8) к точке ветвления $\mathbf{x}^{* *}=(1,1 / 3), \alpha^{* *}=1 / 3$ показана в табл. 5.8. Ход вычислений представлен в табл. 5.9, а поведение ветвей решения изображено графически на рис. 5.8.

Таблица 5.8. Сходимость метода Гаусса – Ньютона к точке ветвления для системы (5.4.31).

Исследуем теперь характер разветвления в точках первичной бифуркации для задачи о двух реакторах с модельной реакцией типа «брюсселятор» (см. задачу 8). В качестве параметра $\alpha$ выберем параметр $D_{1}$, полагая $D_{2}=D_{1} / \rho, \rho=0,1$. Независимо от $D_{1}$ в задаче существует тривиальное стационарное решение $X_{1}=X_{2}=A, Y_{1}=Y_{2}=B / A$. Матрица Якоби для правых частей оказывается равной (если расположить переменные в порядке $X_{1}, Y_{1}, X_{2}, Y_{2}$ )
\[
\mathbf{J}=\left[\begin{array}{cccc}
-(B+1)+2 X_{1} Y_{1}-\alpha & X_{1}^{2} & \alpha & 0 \\
B-2 X_{1} Y_{1} & -X_{1}^{2}-10 \alpha & 0 & 10 \alpha \\
\alpha & 0 & -(B+1)+2 X_{2} Y_{2}-\alpha & X_{2}^{2} \\
0 & 10 \alpha & B-2 X_{2} Y_{2} & -X_{2}^{2}-10 \alpha
\end{array}\right]
\]

и, следовательно, в случае тривиального решения мы имеем
\[
\mathbf{J}=\left[\begin{array}{cccc}
B-1-\alpha & A^{2} & \alpha & 0 \\
-B & -A^{2}-10 \alpha & 0 & 10 \alpha \\
\alpha & 0 & B-1-\alpha & A^{2} \\
0 & 10 \alpha & -B & -A^{2}-10 \alpha
\end{array}\right] .
\]

Положим $A=2, B=6$. После соответствующих вычислений находим, что значения параметра $\alpha$, при которых матрица Якоби имеет нулевое собственное число, равны 0,0443 и 2,256 . Таким образом, мы имеем две точки бифуркации $(2 ; 3 ; 2 ; 3$; $0,0443)$ и $(2 ; 3 ; 2 ; 3 ; 2,256)$. Направления ветвей решения в первой из них, полученные с помощью описанного выше алгоритма, приведены в табл. 5.9. Эти результаты хорошо согласуются с диаграммой решений, представленной на рис. 5.9; точка

a) Ход вычислений для точки ветвления $(1 ; 0,3333 ; 0,3333$ ) из примера (5.4.31).
\[
\tilde{\mathbf{J}}=\left[\begin{array}{ccc}
2,6667 & -5,3333 & -2,6667 \\
0 & 0 & 0
\end{array}\right], \quad \varphi_{21}^{\prime}=0,5, \quad \varphi_{2 \alpha}^{\prime}=0,5 .
\]
\[
\begin{aligned}
& A=-1,3333, \quad B=2, \quad C=0 ; \quad S_{1}=d x_{1} / d \alpha=0, \quad d x_{2} / d \alpha=-0,5 ; \quad S_{2}= \\
=d \alpha / d x_{1}=0,3333, & d x_{2} / d x_{1}=0,3333 .
\end{aligned}
\]

Стартовые точки ( $\left.h_{1}=h_{\alpha}=0,01\right)$ :
b) Направления разветвления в точке первичной бифуркации $(2,3,2,3$. 0,04433 ) для задачи 8. $A=2, B=6, \rho=0,1$.

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

Рис. 5.9. Диаграмма стационарных решений задачи $8, N=2 ; A=2, B=6$, $\rho=0,1$.
При построении диаграммы решений важной проблемой является отыскание всех ветвей решения. Наряду со случайным выбором начальных приближений для метода Ньютона (см. §5.1) мы можем теперь воспользоваться еще одной возможностью.

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

5.4.3. Возникновение изол на диаграмме решений

Возникновение и существование изол на диаграмме решений мы попытаемся объяснить с помощью двухпараметрической системы уравнений
\[
f_{i}\left(x_{1}, x_{2}, \ldots, x_{n}, \alpha_{1}, \alpha_{2}\right)=0, \quad i=1,2, \ldots, n .
\]

Принимая во внимание сложности графического представления такой системы, будем полагать в ней $n=1$, т. е. рассматривать уравнение вида
\[
f_{1}\left(x_{1}, \alpha_{1}, \alpha_{2}\right)=0 .
\]

Множество $S\left(f_{1}\right)=\left\{\left(x_{1}, \alpha_{1}, \alpha_{2}\right) \in \mathrm{R}^{3}, f_{1}\left(x_{1}, \alpha_{1}, \alpha_{2}\right)=0\right\}$ представляет собой в общем случае некоторую поверхность в $R^{3}$. Предположим, что эта поверхность имеет вид параболоида вращения с вершиной в точке P, ось которого параллельна оси $\alpha_{2}$ (см. рис. 5.10). Зафиксируем теперь параметр $\alpha_{2}$ в формуле (5.4.32a), полагая $\alpha_{2}=\tilde{\alpha}_{2}$. Мы получаем уравнение с одним параметром, диаграмма решений которого представляет собой пересечение множества $\mathbf{S}\left(f_{1}\right)$ (параболоида вращения) с плоскостью $\alpha_{2}=\tilde{\alpha}_{2}$. На рис. 5.10 видно, что при $\tilde{\alpha}_{2}<\alpha_{2}^{*}$ это пересечение пусто, при $\tilde{\alpha}_{2}=\alpha_{2}^{*}$ оно вырождается в точку $\mathrm{P}$, а при $\tilde{\alpha}_{2}>\alpha_{2}^{*}$ указанное пересечение представляет собой некоторую замкнутую кривую. Ортогонально проектируя эту замкнутую кривую на плоскость. $x_{1}-\alpha_{1}$, мы получаем в этой плоскости искомую изолу. Описанная ситуация схематически представлена на рис. 5.11.

Сформулируем теперь необходимые условия, позволяющие определить точку $\mathrm{P}^{11}$ для системы общего вида (5.4.32). В окрестности этой точки переменную $\alpha_{2}$ можно рассматривать.

1) $\mathrm{P}=\left(\mathrm{x}^{*}, \alpha_{1}^{*}, \alpha_{2}^{*}\right)$ – точка на двумерной поверхности $\mathrm{S}=\left\{\mathbf{x}, \alpha_{1}, \alpha_{2} \in\right.$ $\left.\in \mathrm{R}^{n+2} ; \mathrm{f}\left(\mathrm{x}, \alpha_{1}, \alpha_{2}\right)=0\right\}$, являющаяся изолированной точкой в пересечении S с плоскостью $\alpha_{2}=\alpha_{2}^{*}$.

Рис. 5.10. Возникновение изол в двухпараметрической проекции.

Рис. 5.11. Схематическое
изображение возникновения
изол на диаграмме решений.

как функцию переменных $x_{k}, \alpha_{1}$ (при некотором фиксированном $k$ )
\[
\alpha_{2}=\varphi\left(x_{k}, \alpha_{1}\right)
\]

Действительно, выбрав значения $x_{k}$ и $\alpha_{1}$, из уравнений (5.4.32) можно определить оставшиеся неизвестные $x_{1}, \ldots, x_{k-1}, x_{k+1}, \ldots$ $\ldots, x_{n}$ и $\alpha_{2}$ и таким образом получить функциональные зависимости $x_{j}=x_{j}\left(x_{k}, \alpha_{1}\right)$ и $\alpha_{2}=\alpha_{2}\left(x_{k}, \alpha_{1}\right) \quad\left(=\varphi\left(x_{k}, \alpha_{1}\right)\right)$. Функция $\varphi$ имеет экстремум в точке $\mathrm{P}$ и, следовательно, для нее должнығ выполняться соотношения
\[
\frac{\partial \varphi}{\partial x_{k}}=0, \quad \frac{\partial \varphi}{\partial \alpha_{1}}=0 .
\]

Положим $x_{n+1}=\alpha_{2}$, считая, что матрица $\mathbf{J}_{k}$ задана выражением (5.2.12). Дифференцируя систему (5.4.32) по переменной $x_{k}$, мы получаем систему линейных алгебраических уравнений
\[
\mathbf{J}_{k} \cdot \mathbf{r}=-\frac{\partial \mathrm{f}}{\partial x_{k}}
\]

относительно составляющих вектора
\[
\mathbf{r}=\left(r_{1}, \ldots, r_{n}\right)=\left(\frac{\partial x_{1}}{\partial x_{k}}, \ldots, \frac{\partial x_{k-1}}{\partial x_{k}}, \frac{\partial x_{k+1}}{\partial x_{k}}, \ldots, \frac{\partial x_{n}}{\partial x_{k}}, \frac{\partial \alpha_{2}}{\partial x_{k}}\right) .
\]

Первое из условий (5.4.34) можно при этом записать в виде
\[
r_{n}\left(x_{1}, x_{2}, \ldots, x_{n}, \boldsymbol{\alpha}_{1}, \alpha_{2}\right)=0^{1)} .
\]

Дифференцируя (5.4.32) по $\alpha_{1}$, получим систему
\[
\mathbf{J}_{k} \cdot \mathbf{s}=-\frac{\partial \mathrm{f}}{\partial \alpha_{1}},
\]

где использовано обозначение
\[
\mathbf{s}=\left(s_{1}, \ldots, s_{n}\right)=\left(\frac{\partial x_{1}}{\partial \alpha_{1}}, \ldots, \frac{\partial x_{k-1}}{\partial \alpha_{1}}, \frac{\partial x_{k+1}}{\partial \alpha_{1}}, \ldots, \frac{\partial x_{n}}{\partial \alpha_{1}}, \frac{\partial \alpha_{2}}{\partial \alpha_{1}}\right)
\]

и, следовательно, второе из условий (5.4.34) записывается в виде
\[
s_{n}\left(x_{1}, x_{2}, \ldots, x_{n}, \alpha_{1}, \alpha_{2}\right)=0 .
\]
1) Компоненты $\mathrm{r}$ имеют смысл производных $\frac{\partial x_{j}}{\partial x_{k}}$ (или $\frac{\partial \alpha_{2}}{\partial x_{k}}$ ) на поверхности $\mathbf{f}=0$. Однако систему (5.4.35) можно рассматривать и вне этой: поверхности; тогда $\mathbf{r}=\mathbf{r}\left(\mathbf{x}, \alpha_{1}, \alpha_{2}\right)$. Именно так понимается $r_{\text {п }}$ в (5.4.37).-.

Отметим, что системы (5.4.35) и (5.4.38) имеют одну и ту же матрицу и, следовательно, решение их может быть найдено с помощью метода исключения Гаусса для обеих правых частей одновременно.

Итак, для определения координат точки Р мы имеем в самом общем случае $n+2$ уравнений (5.4.32), (5.4.37) и (5.4.40) относительно $n+2$ неизвестных $x_{1}, x_{2}, \ldots, x_{n}, \alpha_{1}, \alpha_{2}$ [5.13]. Задавая значения этих неизвестных, мы можем найти левые части указанных $n+2$ уравнений и затем использовать для нахождения решений этой системы соответствующие итерационные процедуры, например обобщенный метод секущих или метод Ньютона с матрицей Якоби, подсчитываемой с помощью разностных формул (при аналитическом вычислении производных для соотношений (5.4.37) и (5.4.40) нам потребовалось бы вычислять частные производные второго порядка от функций $f_{i}$ ).

Изолы могут существовать либо при $\alpha_{2}>\alpha_{2}^{*}$, либо при $\alpha_{2}<$ $<\alpha_{2}^{*}$. Различить между собой эти два случая можно опытным путем, используя метод Ньютона для $\alpha_{2}>\alpha_{2}^{*}$ и $\alpha_{2}<\alpha_{2}^{*}$, причем существование решения доказывается сходимостью к этому решению на изоле ${ }^{1 \text { ). }}$. Другая возможность заключается в том, чтобы вычислить частные производные второго порядка $\partial^{2} \varphi / \partial x_{k}^{2}, \partial^{2} \varphi / \partial x_{k} \partial \alpha_{1}, \partial^{2} \varphi / \partial \alpha_{1}^{2}$ и использовать условия Сильвестра для минимума и максимума функции $\varphi\left(x_{k}, \alpha\right)$. Так, если
\[
\frac{\partial^{2} \varphi}{\partial x_{k}^{2}}>0 \quad \text { и } \quad \frac{\partial^{2} \varphi}{\partial x_{k}^{2}} \cdot \frac{\partial^{2} \varphi}{\partial \alpha_{1}^{2}}-\left[\frac{\partial^{2} \varphi}{\partial x_{k} \partial \alpha_{1}}\right]^{2}>0,
\]

то функция $\varphi$ имеет минимум и, следовательно, изолы существуют при $\boldsymbol{\alpha}_{2}>\boldsymbol{\alpha}_{2}^{*}$. Если же первое из неравенств (5.4.41) имеет обратный знак, то изолы существуют при $\alpha_{2}<\alpha_{2}^{*}$.

Рассмотрим пример использования описанного выше алгоритма для нахождения точки возникновения изол. Обратимся вновь к задаче 1 в формулировке с использованием времени задержки $\tau$, которое в данном случае играет роль параметра $\alpha_{1}$. Обращаясь к рис. 5.2, мы видим, что критическое значение параметра $a$ (выступающего здесь в роли параметра $\alpha_{2}$ ) приблизительно равно 5. Воспользуемся соотношением (5.2.3). Стационарные состояния при этом описываются уравнением (5.2.4).

Используя формулы $\mathrm{Da}=k_{0} \tau, \beta=\alpha \tau$, имеем из (5.2.4)
\[
\begin{aligned}
f_{1}(\Theta, \tau, a)=-\Lambda \Theta+k_{0} \tau & \left(B-\Theta-\frac{a \tau\left(\Theta-\Theta_{c}\right)}{\Lambda}\right) \times \\
\times & \exp \frac{\Theta}{1+\theta / \gamma}-a \tau\left(\Theta-\Theta_{c}\right)=0 .
\end{aligned}
\]

Поскольку $n=1$, мы получаем $x_{k}=x_{1}=\Theta$. Уравнение принимает вид
\[
\frac{\partial f_{1}}{\partial a} \cdot r_{1}=-\frac{\partial f_{1}}{\partial \Theta}
\]

и, следовательно, соотношение (5.4.37) эквивалентно условию
\[
\begin{aligned}
\frac{\partial f_{1}}{\partial \Theta}=-\Lambda+\left[-\left(1+\frac{a \tau}{\Lambda}\right)\right. & \left.+\left(B-\Theta-\frac{a \tau\left(\Theta-\Theta_{c}\right)}{\Lambda}\right) \frac{1}{(1+\Theta / \gamma)^{2}}\right] \times \\
& \times k_{0} \tau \exp \frac{\Theta}{1+\Theta / \gamma}-a \tau=0 .
\end{aligned}
\]

Аналогичным образом, уравнение (5.4.38) записывается в виде
\[
\frac{\partial f_{1}}{\partial a} \cdot s_{1}=-\frac{\partial f_{1}}{\partial \tau},
\]

а условие (5.4.40) дает
\[
\frac{\partial f_{1}}{\partial \tau}=k_{0}\left[B-\Theta-\frac{2 a \tau\left(\Theta-\Theta_{c}\right)}{\Lambda}\right] \exp \frac{\Theta}{1+\Theta / \gamma}-a\left(\Theta-\Theta_{c}\right)=0 .
\]

Таким образом, мы имеем три уравнения (5.4.42), (5.4.44) и (5.4.46) для трех неизвестных $\Theta, \tau$ и $a$. Ход итераций, полученных при использовании метода Ньютона с матрицей Якоби, которая вычислялась с помощью соответствующих разностных формул, представлен в табл. 5.10. Полученные результаты соответствуют диаграмме решений, изображенной на рис. 5.2.

Таблица 5.10. Расчет точки возникновения изол для задачи 1. (Метод Ньютона для системы уравнений (5.4.42), (5.4.44) и (5.4.46); $\gamma=20, B=10, \Lambda=1$, $\left.\boldsymbol{\Theta}_{c}=-5, k_{0}=1\right)$.

В работе [5.13] был рассмотрен еще один алгоритм для определения точки возникновения изол. Так, на зависимости

Рис. 5.12. Диаграмма стационарных решений задачи $3 ; \mu=8,4 \cdot 10^{-6}, f=2$, $\varepsilon=6,6667 \cdot 10^{-4}, \varepsilon^{\prime}=1,7778 \cdot 10^{-5}$. Точка P-см. рис. 5.10, точки P’ и $\mathrm{Q}-$ см. рис. $5.20 b$.
$\mathbf{x}\left(\alpha_{2}\right)$ при $\alpha_{1}=\alpha_{1}^{*}$ точка $\mathrm{P}$ представляет собой точку поворота (см. также рис. 5.10). Необходимое условие для точки Р имест

при этом вид
\[
f_{n+1}\left(x_{1}, \ldots, x_{n}, \alpha_{1}, \alpha_{2}\right)=\operatorname{det} \mathbf{J}\left(x_{1}, \ldots, x_{n}, \alpha_{1}, \alpha_{2}\right)=0 \text {, }
\]

где $\mathbf{J}=\left[\partial f_{i} / \partial x_{j}\right]$ – матрица Якоби. При фиксирөванном $x_{k}=x_{k}^{*}$ зависимость $\left(x_{1}, \ldots, x_{k-1}, x_{k+1}, \ldots, x_{n}, \alpha_{1}\right.$ ) от $\alpha_{2}$ также имеет точку поворота, и поэтому должно быть выполнено условие
\[
f_{n+2}\left(x_{1}, \ldots, x_{n}, \alpha_{1}, \alpha_{2}\right)=\operatorname{det} \mathbf{J}_{k}\left(x_{1}, \ldots, x_{n}, \alpha_{1}, \alpha_{2}\right)=0,
\]

где матрица $\mathbf{J}_{k}$ определяется по формуле (5.2.12) при $x_{n+1}=\alpha_{1}$. Объединив (5.4.32), (5.4.47) и (5.4.48), мы вновь получаем $n+2$ нелинейных уравнений относительно $n+2$ неизвестных. Отметим, что уравнения (5.4.47) и (5.4.48) формально тождественны соотношениям (5.4.1) и (5.4.2) для нахождения точек ветвления. Однако здесь мы имеем на одну переменную больше (добавится $\alpha_{2}$ ), и для решения этой системы мы используем метод Ньютона (в отличие от п. 5.4.1, где применялся метод Гаусса Ньютона). Читатель, который воспользуется вторым методом для нахождения точек возникновения изол в уравнении (5.4.42), получит те же результаты, что и в случае применения первого метода, поскольку здесь $n=1$.

Применим оба описанных подхода для анализа задачи 3, где роль параметра $\alpha_{1}$ будет играть параметр $\beta$, а роль параметра $\alpha_{2}$ – параметр $\alpha$. Пример расчета итераций по методу Ньютона (матрица Якоби вычисляется с помощью соответ-

Таблица 5.11. Расчет точек возникновения изол в задаче 3 ( $\mu=8,4 \times 10^{-6}$, $\left.f=2, \quad \varepsilon=6,6667 \times 10^{-4}, \quad \varepsilon^{\prime}=1,7778 \times 10^{-5}\right)$. Ход итераций по методу’ Ньютона.

ствующих разностных формул) приведен в табл. 5.11. Соответствующая диаграмма решений для различных значений параметра $\alpha$ изображена на рис. 5.12.