Рассмотрим автономную систему $n$ обыкновенных дифференциальных уравнений
\[
\frac{d \mathbf{x}}{d t}=\mathbf{f}(\mathbf{x}, \boldsymbol{\alpha})
\]

с одним параметром $\alpha$ (значение которого мы пока будем считать фиксированным). Периодическое решение $\mathbf{x}(t)$ системы (5.8.1), имеющее период $T$, при любом $t \in \mathrm{R}$ удовлетворяет условию
\[
\mathbf{x}(t+T)=\mathbf{x}(t) .
\]

Устойчивые периодические решения (т. е. притягивающие, или асимптотически орбитально устойчивые решения. – Ред.) (см. § 2.3) можно найти посредством процесса установления, имитируя (численно) динамику системы (5.8.1) с помощью методов, описанных в $\S 5.7^{11}$. После установления колебаний нетрудно оценить (найти с помощью интерполяции) и величину периода $T$. Неустойчивое периодическое решение этим способом найти нельзя, за исключением случаев, когда оно устойчиво при $t \rightarrow-\infty$ (т. е. тогда, когда ни один мультипликатор не лежит внутри единичного круга, при $n=2$ для неустойчивого периодического решения это всегда выполняется). Такие решения можно найти, интегрируя уравнение (5.8.1) в отрицательном направлении по оси времени.

Примером приложения указанного подхода служат периодические решения задачи 1, представленные в форме фазовых портретов на рис. 5.19. На рис. 5.22 изображены два вложенных один в другой предельных цикла. Внешний из них является устойчивым и был найден с помощью установления при $t \rightarrow \infty$. Внутренний цикл неустойчив и был найден процессом установления в обращенном времени (при $t \rightarrow-\infty$ ).

Переходя к общему случаю, преобразуем независимую переменную $t$ следующим образом:
\[
t=T z \text {. }
\]
1) Разумеется, для этого нужно знать точку, достаточно близкую к искомому решению.

Тогда соотношения (5.8.1) и (5.8.2) принимают вид
\[
\begin{array}{c}
\frac{d \mathbf{x}}{d z}=T \mathbf{f}(\mathbf{x}, \alpha), \\
\mathbf{x}(1)=\mathbf{x}(0),
\end{array}
\]
т. е. решение $\mathbf{x}(z)$ имеет период, равный единице (учитывая автономность системы, значение $t$ в формуле (5.8.2) можно выбрать равным нулю).

Рис. 5.22. Периодические решения задачи 1. Устойчивый (сплошная линия) и неустойчивый (штриховая линия) предельные циклы. $\Lambda=1, B=14, \gamma \rightarrow$ $\rightarrow \infty, \mathrm{Da}=0,1618, \beta=3, \Theta_{c}=0 ;+$ устойчивое стационарное решение.

Соотношения (5.8.4), (5.8.5) определяют собой нелинейную краевую задачу со смешанными граничными условиями (более подробно нелинейные краевые задачи рассматриваются в § 6.1). Для решения такого рода задач чаще всего используются разностные методы и метод стрельбы.

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

Идея разностных методов заключается в замене производных конечными разностями на достаточно густой сетке узловых точек $z_{i}=i h, i=0,1, \ldots, N, h=1 / N$. Простейшая такая замена дает вместо (5.8.4)
\[
\frac{\mathrm{x}^{i+1}-\mathrm{x}^{i}}{h}=T \mathrm{f}\left(\frac{1}{2}\left(\mathrm{x}^{i+1}+\mathrm{x}^{i}\right), \alpha\right), \quad i=0,1, \ldots, N-1 ;
\]

здесь $\mathbf{x}^{i} \sim \mathbf{x}\left(z_{i}\right)$. Из граничных условий (5.8.5) имеем
\[
\mathbf{x}^{0}=\mathbf{x}^{N} .
\]

Соотношения (5.8.6), (5.8.7) представляют собой систему $n(N+1)$ нелинейных уравнений относительно $n(N+1)+1$ неизвестных $\mathbf{x}^{0}, \mathbf{x}^{\prime}, \ldots, \mathbf{x}^{N}, T$. Одну из неизвестных мы можем фиксировать заранее (это не может быть $T$, поскольку периодическое решение существует только при дискретных, заранее неизвестных значениях $T^{11}$. Однако для того, чтобы такой подход был успешным, это фиксированное значение должно соответствовать существующему периодическому решению. Ясно,
Рис. 5.23. Схема размещения для системы (5.8.6), (5.8.7).

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

Для решения системы нелинейных уравнений (5.8.6)(5.8.7) можно воспользоваться, например, итерационным методом Ньютона. Матрица размещения для этих уравнений является почти ленточной (ленточный характер нарушается лишь вхождением неизвестной $T$ и условиями (5.8.7), см. рис. 5.23). Для решения линейных систем в случае ленточной матрицы обычно используется одна из модификаций метода исключения Гаусса.

Для уменьшения числа решаемых уравнений можно использовать неравномерную сетку узловых точек на интервале $[0,1]$. Можно также вместо схемы $O\left(h^{2}\right)$ в формуле (5.8.6) использовать схемы высших порядков.

В табл. 5.22 приведены полученные в результате профили периодического решения (при различных значениях шага сетки h) для случая задачи 10 ( $\left.x=x_{1}, y=x_{2}, z=x_{3}\right)$. Сравнивая

Таблица 5.22. Профили периодического решения, найденные разностным методом (5.8.6), (5.8.7) (задача 10: $\sigma=16, b=4, r=200,0$ ). Значение $x(t=0)=0$ в ходе итераций фиксировано $(\tau \in[0,1], t=T \tau)$.

значения периода $T$, полученные для различных величин шага $h$, мы видим, что при уменьшении шага $h$ период $T$ стремится к некоторому предельному значению, приближенно равному $T_{\operatorname{Re}}=0,8226$. Этот предел получается экстраполяцией по Ричардсону из значений $T$ при $h=0,005$ и $h=0,0025$. Сравнивая теперь найденные значения периода $T$ с результатами, полученными с помощью метода стрельбы (см. табл. 5.24), мы видим, что точность, обеспечиваемая разностным методом, является не очень высокой.

5.8.2. Метод стрельбы

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

и определенное значение периода $T$. Тогда систему дифференциальных уравнений (5.8.4) можно проинтегрировать от точки $z=0$, где мы имеем полные начальные условия (5.8.8), до точки $z=1$, используя при этом методы, описанные в §5.7. (В случаях, когда задача Коши «более устойчива» в отрицательном направлении, используется интегрирование от $z=0$ до $z=-1$. Читатель может легко внести в алгоритм соответствующие изменения.) В результате интегрирования мы получаем значения неизвестных $x_{i}$ в точке $z=1$ :
\[
x_{i}(1)=\varphi_{i}\left(\eta_{1}, \ldots, \eta_{n}, T\right), \quad i=1,2, \ldots, n .
\]

Для выполнения условия (5.8.5) необходимо, чтобы удовлетворялась система $n$ нелинейных уравнений
\[
F_{i}\left(\eta_{1}, \ldots, \eta_{n}, T\right)=\varphi_{i}\left(\eta_{1}, \ldots, \eta_{n}, T\right)-\eta_{i}=0, \quad i=1,2, \ldots, n
\]

относительно $n+1$ неизвестных $\eta_{1}, \ldots, \eta_{n}, T$. Так же, как и в случае разностных методов, нужно задать произвольно значение одной из неизвестных. Выберем в качестве такой неизвестной переменную $\eta_{k}$ и положим $\eta_{k}=\bar{\eta}_{k}$. Этот выбор будет успешным лишь в случае, когда $\bar{\eta}_{k}$ «лежит» на искомой зависимости $x_{k}(z)$, т. е. когда при некотором $\bar{z} \in[0,1]$ выполнено соотношение $\bar{\eta}_{k}=x_{k}(\bar{z})$ (см. рис. 5.24). Для решения системы нелинейных уравнений (5.8.10) относительно неизвестных $\eta_{1}, \ldots$ $\ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}, T$ воспользуемся методом Ньютона. При этом соответствующая матрица Якоби будет включать в себя

элементы $\partial F_{i} / \partial \eta_{j}(j
eq k)$ и $\partial F_{i} / \partial T$. Эти элементы можно вычислить либо с помощью разностных формул, когда мы решаем задачу Коши (т. е. вычисляем $F_{i}$ ) для одной итерации в общей сложности $n+1$ раз (см. §5.1), либо с помощью соответствующих вариационных переменных и системы дифференциальных уравнений в вариациях. Покажем, как это можно сделать в по-
Рис. 5.24. Схематическое изображение адаптивного выбора $\eta_{k}$.

следнем случае. Для этого, рассматривая $\mathbf{x}$ как функцию $\eta$, введем обозначения
\[
p_{i j}(z)=\frac{\partial x_{i}}{\partial \eta_{j}}, \quad i, j=1,2, \ldots, n .
\]

Дифференцируя уравнение (5.8.4) по переменной $\eta_{j}$, находим
\[
\frac{d p_{i j}}{d z}=T \sum_{s=1}^{n} \frac{\partial f_{i}}{\partial x_{s}} p_{s j}, \quad i, j=1,2, \ldots, n .
\]

Далее, дифференцируя (5.8.4) по $T$, для производных $q_{i}=$ $=\partial x_{i} / \partial T$ получаем
\[
\frac{d q_{i}}{d z}=f_{i}+T \sum_{s=1}^{n} \frac{\partial f_{i}}{\partial x_{s}} q_{s}, \quad i=1,2, \ldots, n .
\]

Начальные условия для этих уравнений в вариациях имеют вид
\[
p_{i j}(0)=\delta_{i j}, \quad q_{i}(0)=0,
\]

где $\delta_{i j}$ – символ Кронекера (равный 0 при $i
eq j$ и 1 при $i=j$ ).

Для элементов матрицы Якоби системы (5.8.10) имеют место соотношения ${ }^{1)}$
\[
\begin{array}{c}
\frac{\partial F_{i}}{\partial \eta_{j}}=p_{i j}(1)-\delta_{i j}, \quad i, j=1,2, \ldots, n, \\
\frac{\partial F_{i}}{\partial T}=q_{i}(1)=f_{i}(\mathbf{x}(1), \alpha), \quad i=1,2, \ldots, n .
\end{array}
\]

Из последней формулы следует, что уравнения (5.8.13) интегрировать не нужно, поскольку достаточно знать лишь концевое значение $\mathbf{x}(1)$.

Заметим, что столбец производных $\partial F_{i} / \partial \eta_{k}$ матрицы Якоби также оказывается ненужным, поскольку значение $\eta_{k}$ уже зафиксировано. Выбор этого значения $\eta_{k}$ является довольно непростым делом, если только не использовать какую-либо априорную информацию об искомых периодических решениях (вытекающую, например, из неких физических соображений).

Если никакой предварительной информации у нас нет, то можно воспользоваться «поисковым» подходом, описанным в $\$ 5.1$, т. е. случайным выбором начальных приближений (всех $\eta_{1}, \ldots, \eta_{n}, T$ ) для метода Ньютона в сочетании со случайным выбором индекса $k$. В настоящее время, когда требуемый объем машинного времени перестает быть ограничивающим фактором, указанный подход становится, по-видимому, вполне реалистичным.

Пример сходимости метода Ньютона для задачи 10 показан в табл. 5.23. При этом, хотя начальное приближение выбиралось вблизи устойчивого периодического решения (ср. четвертое решение в табл. 5.24), метод привел к неустойчивому периодическому решению («дважды повторенному»- его минимальный период вдвое меньше найденного). Это обстоятельство, по-видимому, вызывается существованием в рассматриваемой области большого числа периодических решений (см. рис. 5.26). Оно подчеркивает преимущества продолжения периодических решений (см. п. 5.8.4) по сравнению с поиском решений для отдельных значений параметра.

Метод стрельбы является эффективным в тех случаях, когда соответствующие задачи Коши легко интегрируются. Если же функция $\mathbf{x}(z)$ сильно зависит от $\boldsymbol{\eta}$, то задача Коши становится «плохой» и сам метод уже не работает; в частности, так обстоит дело в случае сильно неустойчивых периодических решений. При этом в качестве альтернативного метода можно

Таблица 5.23. Метод стрельбы для нахождения периодического решения задачи 10: $\sigma=16, b=4, r=197, \eta_{k}=-5,5, k=1$. Найденное неустойчивое решение «составлено» из двух идентичных периодических решений с периодом $T / 2$

использовать либо разностный метод, либо метод многократной стрельбы (см. п. 5.8.6).

5.8.3. Устойчивость периодических решений

Исследуем теперь устойчивость уже известного периодического решения, характеризующегося значениями $\bar{\eta}_{1}, \ldots, \bar{\eta}_{n}, \bar{T}$ (см. п. 5.8.2).

Орбитальная устойчивость исследуемого периодического решения определяется собственными числами матрицы монодромии (см. (5.8.9) и §2.3)
\[
\mathbf{B}=\left[\frac{\partial \varphi_{i}}{\partial \eta_{j}}\right]=\left[p_{i j}(1)\right] .
\]

При этом вариационные переменные $p_{i j}(z)$ (5.8.11) находятся для значений $\eta$ и $T$, отвечающих изучаемому периодическому решению. Напомним, что для нахождения собственных чисел матрицы В можно использовать методику, описанную в § 5.3. В случае автономной системы одно собственное число всегда равняется единице ( $\lambda_{1}=1$ ). Остальные $n-1$ собственных чисел $\lambda_{2}, \ldots, \lambda_{n}$ (которые мы в гл. 2 называли мультипликаторами) соответствуют собственным числам линейной части отображения Пуанкаре и определяют устойчивость периодического решения. При этом периодическое решение является устойчивым (точнее, орбитально асимптотически устойчивым), если выполняются неравенства
\[
\left|\lambda_{i}\right|<1, \quad i=2, \ldots, n .
\]

Периодическое решение является неустойчивым, если хотя бы для одного мультипликатора имеет место условие $\left|\lambda_{i}\right|>1$. Если

$\left|\lambda_{i}\right| \geqslant 10^{-4}-10^{-5}$, можно говорить о сильно неустойчивом периодическом решении. В п. 5.8.5 мы рассмотрим различные случаи изменения устойчивости периодических решений при изменении параметров задачи.

5.8.4. Продолжение периодических решений по параметру

Займемся теперь исследованием решения задачи (5.8.1), (5.8.2) (или, что то же самое, (5.8.4), (5.8.5)) в зависнмости от параметра $\alpha$. Опишем алгоритм для нахождения указанной зависимости периодических решений, который основывается на алгоритме DERPAR из § 5.2, использовавшемся для продолжения решений нелинейных (алгебраических) уравнений.

При фиксированном $\eta_{k}$ (значение индекса $k$ также не изменяется) соотношения (5.8.10) принимают вид
\[
F_{i}\left(\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}, T, \alpha\right)=0, i=1,2, \ldots, n \text {. }
\]

Решение этой системы $n$ уравнений относительно $n+1$ неизвестных $\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}, T, \alpha$ можно исследовать с помощью алгоритма DERPAR.

Для процесса продолжения нам потребуется вычислять частные производные функций $F_{i}$ по $\boldsymbol{\eta}, T$ и $\alpha$. Вычисление $\partial F_{i} / \partial \eta_{j}$, $\partial F_{i} / \partial T$ описано в п.5.8.2. (Непосредственно для самого процесса продолжения производные $\partial F_{i} / \partial \eta_{k}$ не нужны; для определения устойчивости по матрице В (5.8.17) они необходимы.) Покажем, как вычислить производную $\partial F / \partial \alpha$. Обозначим $r_{i}=\partial x_{i} / \partial \alpha$; дифференцируя уравнение (5.8.4) по $\alpha$, получим уравнения в вариациях
\[
r_{i}^{\prime}=T \sum_{s=1}^{n} \underset{\partial x_{s}}{{ }_{2}} r_{s}+T \frac{\partial f_{i}}{\partial \alpha}, \quad i=1,2, \ldots, n
\]

с начальными условиями
\[
r_{i}(0)=0 .
\]

К формулам (5.8.15) и (5.8.16) следует теперь добавить соотношения
\[
\frac{\partial F_{i}}{\partial \alpha}=r_{i}(1), \quad i=1,2, \ldots, n .
\]

Успешно продолжать решение системы (5.8.19) (т. е. двигаться вдоль (связной компоненты) кривой (5.8.19).- Ред.) можно до тех пор, пока $\bar{\eta}_{k}$ будет лежать на профиле $x_{k}(z)$ (см. рис. 5.24). Если же $\bar{\eta}_{k}$ «покинет» профиль $x_{k}(z)$, т. е. при данном выборе $\bar{\eta}_{k}$ периодическое решение не будет существовать,

то данный алгоритм не сработает. Поэтому в процессе продолжения мы будем адаптивно «приспосабливать» значение $\bar{\eta}_{k}$ к функции $x_{k}(z)$.

Рассмотрим этот вопрос более подробно. Предположим, что нам известно некоторое решение системы (5.8.19) и найдено соответствующее периодическое решение системы (5.8.4). Опишем, как выбирается значение $\bar{\eta}_{k}$ при следующем значении $\alpha$ (см. рис. 5.24).

1. Возьмем монотонный участок периодического профиля ${ }^{1)}$ на промежутке $\left[0, z_{1}\right] \cup\left[z_{2}, 1\right]$ (где лежит значение $\bar{\eta}_{k}$ ). Через $x_{k}^{\max }$ (или соответственно $x_{k}^{\min }$ ) обозначим максимум (или соответственно минимум) на этом участке. Введем следующие обозначения ( $0<\omega_{1}<1$ ):
\[
\begin{array}{l}
x_{k}^{+}=\frac{1}{2}\left[\left(1-\omega_{1}\right) x_{k}^{\min }+\left(1+\omega_{1}\right) x_{k}^{\max }\right], \\
x_{k}^{-}=\frac{1}{2}\left[\left(1+\omega_{1}\right) x_{k}^{\min }+\left(1-\omega_{1}\right) x_{k}^{\max }\right] .
\end{array}
\]

При этом для $\bar{\eta}_{k}$ потребуем выполнения неравенства
\[
x_{k}^{-}<\bar{\eta}_{k}<x_{k}^{+} .
\]

Это условие означает, что $\bar{\eta}_{k}$ располагается достаточно далеко от локального максимума и локального минимума. Поэтому можно надеяться, что для следующего значения параметра $\alpha$ значение $\bar{\eta}_{k}$ «не пропадает» с профиля $x_{k}$.
2. Обозначим через $H_{\max }$ максимальную разность между максимумом и минимумом на монотонном участке профиля $x_{k}(z)$, если этот монотонный участок выбирать из всего промежутка $z \in[0,1]$ :
\[
H_{\max }=x_{k}\left(z^{+}\right)-x_{k}\left(z^{-}\right) .
\]

Потребуем теперь, чтобы
\[
x_{k}^{\max }-x_{k}^{\min }>\left(1-\omega_{2}\right) H_{\max } .
\]

Это второе условие означает, что $\bar{\eta}_{k}$ лежит на достаточно широком монотонном участке.

Опыт показывает, что обычно разумно выбирать значения $\omega_{1}$ и $\omega_{2}$ в диапазоне $[0.5,0.8]$.

Если хотя бы одно из условий (5.8.24) и (5.8.25) не выполняется, то значение $\bar{\eta}_{k}$ выбирается заново по формуле
\[
\bar{\eta}_{k}^{\text {новое }}=\frac{1}{2}\left[x_{k}\left(z^{+}\right)+x_{k}\left(z^{-}\right)\right] .
\]
1) Речь идет о монотонности функции $x_{k}(z),-$

Остальные составляющие $\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}$ должны пересчитываться с использованием профилей $x_{1}(z), \ldots, x_{n}(z)$, которые были найдены для последнего значения параметра $\alpha$. Для этого пересчета выберем координату $z^{*}$ между $z^{-}$и $z^{+}$таким образом, чтобы выполнялось условие
\[
x_{k}\left(z^{*}\right) \approx \bar{\eta}_{k}^{\text {новое }} \text {. }
\]

Тогда значения $\eta_{i}^{\text {новое }}=x_{i}\left(z^{*}\right), i=1, \ldots, n$, будут представлять собой стартовые значения для дальнейшего процесса продолжения. При этом многошаговая схема интегрирования (например, схема Адамса-Бэшфорта), используемая в качестве предиктора в алгоритме DERPAR, вновь начинается со схемы первого порядка (т. е. схемы Эйлера). При такой организации вычислений должны также пересчитываться и некоторые управляющие параметры в алгоритме DERPAR. Прежде всего это относится к параметрам направления $N_{i}$, поскольку нам необходимо сохранять направление движения вдоль ветви периодических решений. Все детали указанного подхода читатель может найти в работе [5.17]. На рис. 5.25 схематически изображена схема описанного алгоритма нахождения периодических решений, который далее мы будем называть алгоритмом DERPER.

Отметим, что в процессе продолжения для каждого найденного периодического решения легко вычисляются его мультипликаторы. Действительно, элементы матрицы монодромии В = $=\left(p_{i j}(1)\right)$ получаются попутно (см. формулы (5.8.12), (5.8.17) ), после чего остается лишь найти собственные значения матрицы В.

Трудности с нахождением параметрической зависимости возникают иногда в окрестности точек бифуркации (например, точки бифуркации с потерей симметрии), однако в большинстве случаев описываемый алгоритм отслеживает (гладкое) продолжение исходной ветви решений. Алгоритм беспрепятственно проходит точки рождения инвариантного тора и точки ответвления решений с двукратным периодом и прослеживает продолжение исходной ветви решений; в случае точки поворота (отвечающей бифуркации «+1» (слияние двух периодических решений).-Pед.) алгоритм, пройдя эту точку, переходит на другую ветвь. Окончание ветви решений в точке бифуркации Андронова-Хопфа и в точке, где периодическое решение превращается в гомоклиническую траекторию, также контролируется этим алгоритмом (см. рис. 5.25). Проиллюстрируем использование алгоритма DERPER на нескольких примерах.

Модель Лоренца (задача 10, см. (5.5.1)) обладает симметрией в следующем смысле. Введем отображение $\mathrm{g}: \mathbf{R}^{3} \rightarrow \mathbf{R}^{3}$
Рис. 5.25. Схема применения алгоритма DERPER.
(симметрию относительно оси $z$ ) посредством соотношения $\mathrm{g}(x, y, z) \equiv(-x,-y, z)$. Если $\mathbf{P}(t)=(x(t), y(t), z(t))$ – некото-

рое решение уравнений (5.5.1) с траекторией $\gamma$, то $\overline{\mathbf{P}}(t)=$ $=(-x(t),-y(t), z(t))$ – также решение (с траекторией $\mathbf{g}(\gamma))$. В частностн, еслн $\mathbf{P}(t)$ – периодическое решение, то $\overline{\mathbf{P}}(t)$, очевидно, тоже. При этом могут иметь место два случая:
1. $\mathbf{g}(\gamma)=\gamma$, т. е. решения $\mathbf{P}$ и $\overline{\mathbf{P}}$ имеют одну и ту же траекторию, которая симметрична относительно оси $z$; при этом $\mathbf{P}(t)=\overline{\mathbf{P}}(t+T / 2)$. Такие решения мы будем обозначать символом S; примером служит решенне, представленное на рис. 5.26b.
2. $\boldsymbol{g}(\gamma)
eq \gamma$ и, следовательк, функции $\mathbf{P}$ и $\overline{\mathbf{P}}$ задают два существенно различных периодических решения, траектории которых $\gamma$ и $g(\gamma)$ взаимно симметричны относительно оси $z$ (такие решения мы будем обозначать символом А). Примером может служить решение на рис. 5.26c.

Другие модели также обладают симметрией, в соответствии с которои́ с данным решением сосуществуют несколько других симметричных решений.

Опишем теперь зависимость периодических решений от параметра $r$ для задачи 10, показанную на рис. 5.26а.

Из точки бифуркации Андронова-Хопфа субкритически отделяется ветвь неустойчивых периодических решений. С учетом упомянутой симметрии таких точек бифуркации оказывается две ${ }^{1)}$ (на рис. 5.26 а им отвечает одна точка), и из них ответвляются две ветви периодических решений, траектории которых взаимно симметричны по отношению к оси $z$. При уменьшении значения параметра $r$ уменьшается расстояние периодических траекторий до состояния равновесия $x=y=z=0$. Эти периодические траектории сходятся к двум гомоклиническим траекториям, которые входят и выходят из состояния равновесия $x=y=z=0$ (мы обозначим их как $\gamma_{1}$ и $\gamma_{2}$ ). Некоторые ветви периодических решений сходятся при изменении $r$ к множеству, представляющему собой объединение этих гомоклинических траекторий, т. е. к множеству $\gamma_{1} \cup \gamma_{2}$. Амплитуда изменения переменной $y$ на этом множестве вдвое больше, чем на одной гомоклинической траектории. На рис. 5.26а это отображено в виде скачка кривой.
$\mathrm{Ha}$ этом же рисунке в диапазоне параметра $r \in(50,250)$ показано несколько точек поворота (где встречаются две ветви периодических решений). Вблизи каждой точки поворота существует узкий интервал изменения параметра $r$, в котором периодическое решение оказывается устойчивым. Отметим, что рис. 5.26 не показывает все существующие ветви периодических
1) Имеются в виду точки на такой диаграмме периодических решений, которая «различает» два взаимно симметричных решения.

решений [4.40]. Описание того, как периодическое решение теряет устойчивость, мы дадим в следующем пункте.

Рис. 5.26. Периодические решения задачи $10, \sigma=16, b=4$. a) Диаграмма периодических решений. $A_{y}$-амплитуда 2 -й составляющей вектора решения, НВВ – точка бифуркации Хопфа, НО – гомоклиническая орбита, сплошная линия – устойчивое решение, штриховые линии – неустойчивые решения, – -точка бифуркации с потерей симметрии. b) Симметричное устойчивое решение для значения $r=137 ; \mathrm{T}=1,48629, A_{y}=154$. Показана проекция на плоскость $x-y, \eta=(-0,02624,-36,5315,135,279)$, + стационарное решение. c) Несимметричное устойчивое решение для значения $r=87,487$; $T=1,54563 ; \quad A_{y}=101$. Показана проекция на плоскость $x-y, \quad \eta=$ $=(7,7979,-26,7059,99,7729) ;+$ стационарное решение
Верхние ветви решений, идущие из точек поворота, оканчиваются на разных гомоклинических траекториях, нижние же ветви на концах стремятся к паре гомоклинических траекторий («восьмерка»). При $r>500$ существует только одна ветвь пе-

риодических решений, которые являются устойчивыми и симметричными, при $r \sim 470$ эта ветвь теряет устойчивость в точке бифуркации с потерей симметрии. В табл. 5.24 приведены устойчивые периодические решения из всех устойчивых участков ветвей, показанных на рис. 5.26a, а также два неустойчивых периодических решения. Напомним, что устойчивые периодические решения можно искать также, имитируя (численно) динамику уравнений (5.8.1) (процесс установления периодического режима).

Таблица 5.24. Некоторые периодические решения для задачи 101)

Скажем теперь несколько слов о графическом представлении периодических решений и зависимости периодических решений от параметра. Периодические решения (некоторые их составляющие) мы либо представляем графически в зависимости от времени (см. рис. $5.27 \mathrm{~d}$ ), либо изображаем проекцию траектории на фазовую плоскость $x_{i}-x_{j}, i
eq j$ (см. рис. $5.26 \mathrm{~b}, \mathrm{c}$, а также рис. 5.27 c). В последнем случае мы получаем представление о траектории, но теряем информацию о зависимости решения от времени. Для заданного периодического решения часто используется также проекция траектории отображения Пуанкаре на некоторую плоскость (см. § 5.9). При изображении диаграммы периодических решений (см. рис. 5.26a, 5.27a, b и 5.28) на оси ординат мы откладываем какую-либо величину, которая характеризует периодические решения в зависимости от параметра. Это может быть, например, координата какойлибо из неподвижных точек отображения Пуанкаре. Более наглядно изображение значений амплитуды (какой-либо пере-

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

Рис. 5.27. Периодические решения задачи 8, $N=2, A=2, B=5,9, \rho=$ $=D_{1} / D_{2}=0,1$ a) Диаграмма периодических решений. $A_{1}$ – амплитуда переменной $X_{1}$, сплошная линия – изолированное семейство периодических решений с периодом $T \sim 11-13$, штриховая линия – ветви периодических решений с удвоенным периодом $T \sim 22-26$. $O-1$ точки бифуркации удвоения периода. $b$ ) Диаграмма периодических решений. $A_{1}$-амплитуда переменной $X_{1}$; сплошная линия – устойчивое решение, штриховая линия – неустойчивые решения, НВВ – бифуркация Андронова – Хопфа, О-1 – точки бифуркации удвоения периода. c) Траектория периодического решения для значения параметра $D_{1}=1,26 . T=12,31918, A_{1}=6,20, \eta=(3,7335,2,0840$, $1,2159,2,5000)$. Показана прсекция на плоскость $X_{1}-X_{2}$. d) Периодическое решение с рис. с).

ный период (для ветви решения с периодом $T_{2} \sim 2 T_{1}$, которая отделилась от ветви решений с периодом $T_{1}$, модифицированный период определяется как $T_{2}=T_{2} / 2$ ), в результате чего диаграмма остается в этих точках непрерывной.

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

связанных реакторов, в которых имеет место реакция типа «брюсселятор» (задача 8). На рис. $5.27 \mathrm{a}$, b приведены две небольшие части диаграммы периодических решений для этой задачи. Изола, изображенная на рисунке $5.27 \mathrm{a}$, содержит четыре точки поворота и очень близко к ним четыре точки бифуркации удвоения периода (последовательность дальнейших бифуркаций удвоения на рисунке не представлена). На части диаграммы,

Рис. 5.28. Диаграмма периодических решений задачи $2 ; \gamma=1000, B=12$, $\beta_{1}=\beta_{2}=2, \boldsymbol{\theta}_{c}=0, \Lambda=0,8, \mathrm{Da}_{2}=0,2$. Сплошные линии – устойчивые решения, штриховые линии – неустойчивые решения, $A_{4}$ – амплитуда переменной $\Theta_{2}, \mathrm{HBB}$ – точка бифуркации Андронова – Хопфа.

представленной на рис. $5.27 \mathrm{~b}$, изображены две ветви, которые выходят из точки бифуркации Андронова-Хопфа, а затем вновь возникает последовательность бифуркаций удвоения периода (на рисунке показаны только точки первого удвоения). В целом диаграмма периодических решений оказывается гораздо более сложной, чем показано на рис. $5.27 \mathrm{a}, \mathrm{b}$ (см. работу [5.29]). Точки бифуркации, изображенные на рис. $5.27 \mathrm{a}$, мы рассмотрим в следующем пункте.

Диаграмма периодических решений для задачи 2 изображена на рис. $5.28 \mathrm{a}, \mathrm{b}$.

Для сильно неустойчивых периодических решений (мультипликаторы порядка $10^{5}$ и больше) метод простой стрельбы не

работает. В таких случаях для процесса продолжения можно использовать разностные методы или метод многократной стрельбы. Метод многократной стрельбы можно применять также для продолжения периодических решений вблизи гомоклинических орбит. Таким способом можно находить периодические решения с бо́льшими значениями периода $T$, чем при использовании метода простой стрельбы.

При запуске процесса продолжения периодических решений нам необходимо определить начальную точку на ветви, т. е. периодическое решение для одного значения параметра $\alpha$. Относительно легко такое решение можно найти в окрестности точки бифуркации Андронова-Хопфа $\left(x_{1}^{+}, x_{2}^{+}, \ldots, x_{n}^{+}, \boldsymbol{\alpha}^{+}\right)$, методы определения которой изложены в § 5.5. Здесь используются, по существу, два подхода.

Первый из них основан на асимптотическом анализе характера ветвления в данной точке путем учета высших производных. Описание этого метода можно найти в работе [5.33], где приведен также собственно алгоритм вычислений. В результате мы получаем:
a) оценку периода $T$;
b) возможность узнать, является ли ответвившееся решение устойчивым;
c) возможность определить направление ответвления, т. е. узнать, отходит ветвь периодических решений при $\alpha<\alpha^{+}$или при $\alpha>\alpha^{+}$.

Другой подход к данной задаче является эвристическим. Именно, делается попытка найти периодическое решение при $\alpha$, близком к $\alpha^{+}$(см. п. 5.8.2). В качестве начального приближения для величины периода $T$ обычно годится величина $\widetilde{T}=$ $=2 \pi / \sqrt{\omega^{+}}$, где $\lambda_{1,2}= \pm i \sqrt{\omega^{+}}-$собственные числа матрицы Якоби правых частей в точке бифуркации Андронова-Хопфа. Начальное приближение $\eta$ можно выбрать, например, в виде
\[
\begin{array}{c}
\eta_{1}=x_{1}^{+}, \ldots, \eta_{k-1}=x_{k-1}^{+}, \quad \eta_{k}=x_{k}^{+}+\varepsilon, \\
\eta_{k+1}=x_{k+1}^{+}, \ldots, \eta_{n}=x_{n}^{+},
\end{array}
\]

где $\varepsilon
eq 0$ есть некоторое фиксированное малое число, а $\eta_{k}$ в процессе ньютоновских итераций остается неизменным. В зависимости от того, будет ли этот выбор успешным для $\alpha=$ $=\alpha^{+}+\delta$ или для $\alpha=\alpha^{+}-\delta$, выясняется обычно, в какую сторону ответвляются периодические решения. Здесь $\delta
eq 0$ есть опять подходящим образом выбранное малое число. В целом указанный подход дает вполне удовлетворительные результаты

Оба этих подхода оказываются слишком сложными в случае, когда разыскивается периодическое решение на устойчивой ветви. В этом случае достаточно просто интегрировать дифференциальные уравнения вплоть до установления колебаний. ${ }^{[13]}$

Направления ветвлений в точках бифуркации удвоения периода или в точках бифуркации с потерей симметрии приближенно можно подсчитать с помощью методики, описанной в п. 5.4.2.

5.8.5. Бифуркации периодических решений

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

Если одно собственное число матрицы монодромии В проходит по вещественной оси через точку -1 , то от ветви периодических решений отходит ветвь решений с двойным периодом (см. § 2.3). На рис. 5.29c схематически представлены некоторые возможные варианты такого ветвления, при этом если исходная ветвь была устойчивой, она эту устойчивость теряет, а устойчивость новой ветви определяется направлением ветвления. Бифуркация такого типа называется в литературе бифуркацией удвоения периода, или бифуркацией типа «-1». Указанная бифуркация обычно последовательно повторяется в достаточно узком интервале значений параметра $\alpha$, образуя при этом так называемую последовательность Фейгенбаума (см. табл. 5.31). Описанная структура бифуркаций изображена на рис. $5.29 \mathrm{~d}$.

Другие явления возникают в случае, когда одно собственное число проходит по вещественной оси через точку +1 ; как правило, это соответствует точке поворота на диаграмме решений (рис. 5.29a). Если одна из ветвей отвечает устойчивым периодическим решениям, то другая отвечает неустойчивым. Более сложная картина наблюдается в системах, обладающих симметрией (см. задачу 10, которая рассматривается в п. 5.8.4). Схематически эта ситуация представлена на рис. $5.29 \mathrm{~b}$. От ветви устойчивых симметричных решений при бифуркации типа «+1»

отходит ветвь устойчивых несимметричных решений. При этом исходная ветвь симметричных решений теряет устойчивость.

Рис. 5.29. Схематическое изображение бифуркаций периодических решений; сплошные линии – устойчивые решения, штриховые линии – неустойчивые решения, $A_{p}$ – амплитуда или какая-либо другая величина, характеризующая периодические решения, +1 – точка бифуркации с потерей симметрии, – 1точка бифуркации удвоения периода, LB – предельная точка (точка поворота), S – симметричное решение, А-асимметричное решение; $a$ ) точка поворота, b) точка бифуркации с потерей симметрии, c) точка бифуркации удвоения периода, $d$ ) последовательность бифуркаций, удваивающих период, которая следует за бифуркацией с потерей симметрии, $e$ ) последовательность бифуркаций, удваивающих период.
Следующий тип бифуркации возникает при переходе двух комплексно сопряженных собственных чисел через единичную окружность. Здесь от исходной ветви, которая утрачивает устойчивость, отходит ветвь двухчастотных колебаний. Более точно, рождается инвариантный двумерный тор, устойчивый или не-

устойчивый в зависимости от суперкритического или субкритического характера бифуркации (т. е. от того, что больше, $\alpha$ или $\alpha^{*}$; см. п. 2.3.1).

В качестве примера описанных бифуркаций обратимся вновь к рис. 5.26а для задачи 10, приведенному в п. 5.8.4. В окрестности каждой точки поворота на узком интервале изменения $r$ существует устойчивое периодическое решение. Это решение теряет устойчивость либо через бифуркацию удвоения периода (в случае ветви несимметричных решений, см. рис. 5.29e), либо в точках бифуркации «+1» с потерей симметрии (в случае ветви симметричных решений, см. выделенные точки на рис. 5.26a). За этой бифуркацией возникает последовательность бифуркаций типа «-1» (см. рис. 5.29d). Кроме того, на основной ветви на рис. 5.26а при $r \sim 350$ имеется последовательность бифуркаций удвоения периода. На рисунке для последней указаны только ветви решения с периодом $\sim 2 T$ и $\sim 4 T$.

В последнее время внимание исследователей уделялось, в основном, разработке новых методов нахождения точек бифуркации на ветвях диаграммы периодических решений. Опишем здесь два итерационных метода для нахождения точек бифуркации удвоения периода. Оба метода основаны на нахождении периодического решения, для которого один из мультипликаторов строго равен -1 [5.26].
Метод 1 (бифуркация «-1»)
Пусть характеристический многочлен матрицы монодромии В есть $P(\lambda)=(-1)^{n} \operatorname{det}(\mathbf{B}-\lambda \mathbf{I})$.
Соответствующее характеристическое уравнение записывается в виде
\[
P(\lambda)=\lambda^{n}+a_{1} \lambda^{n-1}+\ldots+a_{n-1} \lambda+a_{n}=0 .
\]

Вычисление коэффициентов может производиться произвольным образом, например, с помощью методов, рассмотренных в § 5.3.

Для того чтобы матрица В соответствовала точке бифуркации типа «-1», значение $\lambda=-1$ должно быть корнем уравнения (5.8.28), т. е. должно выполняться следующее соотношение:
\[
F_{n+1}\left(\eta_{1}, \ldots, \eta_{n}, T, \alpha\right)=1+\sum_{i=1}^{n}(-1)^{t} a_{i}=0 .
\]

Уравнения (5.8.19) и (5.8.29) представляют собой систему $n+1$ нелинейных (алгебраических) уравнений относительно $n+1$ неизвестных $\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}, T$, $\alpha$ (неизвестную $\eta_{k}$ мы вновь считаем фиксированной и «лежащей» на искомом

профиле). Для решения этой системы воспользуемся методом Ньютона. Первые $n$ строк матрицы Якоби найдем с помощью соотношений (5.8.15), (5.8.16) и (5.8.22). Элементы последней строки $\partial F_{n+1} / \partial \eta_{j}, \partial F_{n+1} / \partial T, \partial F_{n+1} / \partial \alpha$ можно вычислить, используя, например, соответствующие разностные формулы.
Метод 2 (бифуркация типа «-1»)
Матрица монодромии В имеет в точке бифуркации удвоения периода собственное число, равное – 1 . Это означает, что существует не равный нулю вектор $\mathbf{v}=\left(v_{1}, \ldots, v_{n}\right)$, для которого
\[
(\mathbf{B}+\mathbf{I}) \mathbf{v}=\mathbf{0} .
\]

Вектор $\mathbf{v}$, задаваемый системой (5.8.30), определен с точностью до некоторого множителя. Поэтому одну составляющую вектора v мы можем считать фиксированной в ходе всего процесса вычислений, положив, например,
\[
v_{s}=1 \text {, }
\]

где $1 \leqslant s \leqslant n$. Тем самым мы получаем $2 n$ уравнений (5.8.19), (5.8.30) относительно $2 n$ неизвестных $\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}$, $T, \alpha, v_{1}, \ldots, v_{s-1}, v_{s+1}, \ldots, v_{n}$. Для решения этой системы вновь, как и в методе 1 , воспользуемся схемой Ньютона. При этом по аналогии с методом 1 часть матрицы Якоби вычислим с помощью соотношений (5.8.15), (5.8.16), (5.8.22), а оставшуюся часть с помощью соответствующих разностных формул (существует также возможность их аналитического вычисления).

Проиллюстрируем описанные методы на примере задачи 8, в которой рассматривалась модель двух связанных между собой реакторов. На рис. 5.27 а была представлена диаграмма периодических решений для этой задачи, на которой появились четыре точки бифуркации типа «-1». Координаты этих точек ${ }^{1)}$ приведены в табл. 5.25. Для точного определения критического значения параметра $\alpha=\mathrm{Da}_{1}$ и соответствующего периодического решения в одной из бифуркационных точек были применены оба описанных выше метода. Результаты расчетов представлены в табл. 5.26 и на рис. 5.30. Четыре строчки таблицы отвечают одному и тому же периодическому решению, более точно – четырем периодическим решениям, отличающимся лишь сдвигом во времени. То, что таких решений при $\eta_{1}=2$ должно быть ровно 4 , очевидно из рисунка.

В табл. 5.27 показан ход итераций по схеме Ньютона в случае метода 1. Метод сходится к первой точке из табл. 5.25, а более точно – к четвертому решению из табл. 5.26. При проведении расчетов оказалось, что области начальных приближений, для которых метод Ньютона сходился к выбранным точкам бифуркации «-1», относительно невелики. Это объясняется

Рис. 5.30. Составляющая $X_{1}(z)$ периодической траектории задачи 8, $N=2$, в точке бифуркации типа «-1». Параметры взяты из табл. 5.25. Первая точка из этой таблицы выбрана в качестве начальной точки траектории. Нумерация тех точек, где $X_{1}(z)=2$, используется в табл. 5.26.

тем, что в окрестности выбранных значений параметров существует много точек бифуркации (более подробные результаты можно найти в работе [5.26]).

Таблица 5.25. Точки бифуркации удвоения периода для задачи 8 (рис. 5.27a): $A=2, B=5,9, \rho=D_{1} / D_{2}=0,1, k=1, \eta_{k}=2$.

Можно разработать аналогичные методы для нахождения точек бифуркации «+1» и бифуркации рождения тора. Так, в случае бифуркации «+1» (в автономной системе дифференци-

Таблица 5.26. Четыре различных решения, соответствующие первой точке в табл. 5.25 (см. рис. 5.30): $T=11,59421, \alpha=1,22382$

Таблица 5.27. Вычисление точки бифуркации удвоения периода для задачи 8. Пример итераций в случае метода $1: A=2, B=5,9, \quad \rho=D_{1} / D_{2}=0,1$, $k=1, \eta_{k}=2$. Начальное приближение для метода Ньютона выбиралось случайным образом.

альных уравнений) характеристическое уравнение (5.8.28) имеет двукратный корень +1 , т. е. должно выполняться соотношение
\[
F_{n+1}\left(\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}, T, \alpha\right) \stackrel{\text { def }}{=} \frac{d P}{d \lambda}(1)=0
\]
(условие $P(1)=0$ выполняется автоматически). Уравнения (5.8.19) и (5.8.32) вновь представляют собой систему $n+1$ нелинейных (алгебраических) уравнений относительно $n+1$ неизвестных $\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}, T$, $\alpha$ (переменную $\eta_{k}$ мы опять считаем фиксированной). Далее действуем точно так же, как в методе 1 для точек бифуркации типа «-1».

В случае бифуркации типа тора указанный подход несколько усложняется. Эта бифуркация (см. гл. 2) характеризуется наличием пары комплексно сопряженных собственных чисел $\lambda_{2,3}$ матрицы монодромии, которые переходят через единичную окружность и для которых в точке бифуркации выполняются соотношения
\[
\left|\lambda_{2,3}\right|=1, \quad \lambda_{2,3}^{m}
eq 1, \quad m=1,2,3,4 .
\]

Будем искать периодические решения, у которых два мультипликатора имеют вид
\[
\lambda_{2,3}=a \pm i b, \quad a^{2}+b^{2}=1, \quad b
eq 0 .
\]

Оба последующих метода основаны на использовании подходящего разложения характеристического многочлена.
Метод 3 (бифуркация типа тора)
Выделив из $P(\lambda)$ множитель $\lambda^{2}-\omega \lambda+1$, напишем
\[
P(\lambda)=\left(\lambda^{2}-\omega \lambda+1\right)\left(\lambda^{n-2}+b_{1} \lambda^{n-3}+\ldots+b_{n-2}\right)+C \lambda+D .
\]

Здесь $\omega=2 a$ и $\lambda^{2}-\omega \lambda+1=\left(\lambda-\lambda_{2}\right)\left(\lambda-\lambda_{3}\right)$. Коэффициенты $b_{1}, \ldots, b_{n-2}, C, D$ находятся с помощью рекуррентных формул $\left(b_{0}=1, b_{-1}=0\right)$ :
\[
\begin{array}{c}
b_{m}=a_{m}+\omega b_{m-1}-b_{m-2}, \quad m=1,2, \ldots, n-2, \\
C=a_{n-1}+\omega b_{n-2}-b_{n-3}, \quad D=a_{n}-b_{n-2} .
\end{array}
\]

При этом в точке бифуркации типа тора имеют место соотношения
\[
\begin{array}{l}
F_{n+1}\left(\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}, T, \alpha, \omega\right) \stackrel{\text { def }}{=} C=0, \\
F_{n+2}\left(\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}, T, \alpha, \omega\right) \stackrel{\text { def }}{=} D=0 .
\end{array}
\]

Уравнения (5.8.19) и (5.8.37) представляют собой систему $n+2$ нелинейных (алгебраических) уравнений относительно $n+2$ неизвестных $\eta_{1}, \ldots, \eta_{k-1}, \eta_{k+1}, \ldots, \eta_{n}, T, \alpha, \omega$ (неизвестная $\eta_{k}$ снова считается фиксированной).
Метод 4 (бифуркация типа тора)
При разложении характеристического многочлена $P(\lambda)$ можно использовать то обстоятельство, что (в автономном случае) всегда есть корень $\lambda_{1}=1$. Разложение при этом принимает вид
\[
\begin{array}{r}
P(\lambda)=\left(\lambda^{3}-(2 a+1) \lambda^{2}+(2 a+1) \lambda-1\right) \times \\
\times\left(\lambda^{n-3}+b_{1} \lambda^{n-4}+\ldots+b_{n-3}\right)+C \lambda^{2}+D \lambda+E .
\end{array}
\]

Коэффициенты $b_{1}, \ldots, b_{n-3}, C, D, E$ вновь подсчитываются по рекуррентным формулам ( $b_{0}=1, b_{-1}=0, b_{-2}=0 ; \omega=2 a$ )
\[
\begin{array}{c}
b_{m}=a_{m}+(\omega+1)\left(b_{m-1}-b_{m-2}\right)+b_{m-3}, \quad m=1,2, \ldots, n-3, \\
C=a_{n-2}+b_{n-5}+(\omega+1)\left(b_{n-3}-b_{n-4}\right), \\
D=a_{n-1}+b_{n-4}-(\omega+1) b_{n-3}, \\
E=a_{n}+b_{n-3} .
\end{array}
\]

Если коэффициенты $C$ и $D$ в соотношениях (5.8.37) вычислить по формулам (5.8.39), то уравнения (5.8.19) и (5.8.37) будут представлять собой систему $n+2$ нелинейных уравнений относительно $n+2$ неизвестных, как и в методе 3 (в найденной точке бифуркации типа тора условие $E=0$ выполняется автоматически). Вместо соотношений (5.8.37) мы могли бы взять условия $C=E=0$ или $D=E=0$ (поскольку $P(1)=0$, то $C+$ $+D+E=0$, см. (5.8.38)). Для решения получаемых систем нелинейных уравнений можно вновь воспользоваться методом Ньютона. Можно поступить и так: с помощью метода Гаусса Ньютона решить систему $n+3$ уравнений (5.8.19), (5.8.37) и $E=0$ относительно $n+2$ неизвестных.

В табл. 5.28 приведены точки бифуркации типа тора для задачи 8, найденные описанными выше методами. Отметим, что метод Ньютона является очень чувствительным к выбору начального приближения. При этом подходящее начальное приближение находится обычно с помощью продолжения соответствующей ветви периодических решений. В табл. 5.29 указана одна точка бифуркации типа тора из табл. 5.28 в четырех различных представлениях. Эти различные представления отвечают наличию на траектории периодического решения нескольких точек, в которых $x_{k}(z)$ принимает заданное значение. Для точки бифуркации при $B=5,5$ из табл. 5.28 график функции $x_{1}(z)$ пересекает прямую $\eta_{1}=2$ четыре раза. В табл. 5.29 приведены только значения координат $\eta_{2}, \eta_{3}$ и $\eta_{4}$ (значения $T$ и $\alpha$, конечно, одинаковы для всех представлений).
Таблица 5.28. Точки бифуркации типа тора для задачи 8 ( $A=2$, $\rho=D_{1} / D_{2}=0,1, k=1, \eta_{k}=2$ ).

—————————————————————-
0035_fiz_kol_vol_366_no_photo_page-0236.jpg.txt

5.8. Вычисление периодических решений в автономном случае
235
Таблица 5.29. Точка бифуркации типа тора в четырех представлениях; задача 8 ( $\left.A=2, B=5,5, \quad \rho=D_{1} / D_{2}=0,1, \quad k=1, \quad \eta_{k}=2\right) ; T=8,32408$, $\alpha=D_{1}=0,052495, \omega=-1,08277$.

На рис. 5.31 представлена бифуркационная диаграмма для точек бифуркации типа тора в плоскости параметров $B-D_{1}$. Точки, в которых $\lambda^{3}=1$ или $\lambda^{4}=1$, обозначены как $3 T$ и $4 T$

Рис. 5.31. Бифуркационная диаграмма периодических решений. Бифуркация типа тора. Задача $8, A=2, \rho=D_{1} / D_{2}=0,1$.

соответственно. Изображенная на рисунке кривая заканчивается в обозначенных кружком точках, для которых $\lambda_{1}=\lambda_{2}=$ $=\lambda_{3}=1$.

В заключение данного пункта отметим, что описанные выше методы можно без затруднений использовать для вычисления точек бифуркации и в случае неавтономных систем обыкновенных дифференциальных уравнений (см. § 5.11). Кроме того, на основе результатов п. 5.4.2 в точке бифуркации типа «-1» можно определить направление ветви решений с двойным периодом.

5.8.6. Метод многократной стрельбы

Метод стрельбы, называемый также иногда методом простой стрельбы (см. п. 5.8.2), в случае сильно неустойчивых периодических решений практически не работает. В таких случаях можно использовать конечно-разностный метод, который, однако, не дает нам информации об устойчивости найденного периодического решения; кроме того, размерность решаемых задач оказывается достаточно высокой. Oпределенным компромиссом между методом стрельбы и разностным методом является метод многократной стрельбы. Он позволяет не только находить сильно неустойчивые периодические решения, но и достаточно просто получать информацию об их устойчивости.

Приступая к описанию метода, разобьем интервал $[0,1]$ на $m-1$ подынтервалов с помощью точек $z_{j}$ :
\[
\begin{array}{c}
0=z_{1}<z_{2}<\ldots<z_{m-1}<z_{m}=1, \\
\Delta z_{j}=z_{j+1}-z_{l}, \quad j=1,2, \ldots, m-1 .
\end{array}
\]

В точках $z_{j}, j=1,2, \ldots, m-1$, выберем начальные условия вида
\[
\mathbf{x}\left(z_{j}\right)=\eta_{j}=\left(\eta_{j 1}, \ldots, \eta_{j n}\right) .
\]

Кроме того, зададимся определенным значением периода $T$. Интегрируя уравнения (5.8.4) на интервалах $\left[z_{j}, z_{j+1}\right]$ с начальными условиями (5.8.40), можно найти решение в точке $z_{j+1}$; обозначим его как $\mathbf{x}\left(z_{j+1} \mid \boldsymbol{\eta}_{i}, T, \alpha\right)$. Для того чтобы решение в точке $z_{j+1}$ было непрерывным, потребуем выполнения условий
\[
F_{j}=x\left(z_{j+1}^{\prime} \mid \eta_{j}, T, \alpha\right)-\eta_{j+1}=0, \quad j=1,2, \ldots, m-1 .
\]

Кроме того, из условий периодичности (5.8.5) следует, что $\boldsymbol{\eta}_{1}=\boldsymbol{\eta}_{m}$. Уравнения (5.8.41) представляют собой систему $n(m-1)$ нелинейных (алгебраических) уравнений относительно $n(m-1)+1$ неизвестных $\left(\eta_{1}, \boldsymbol{\eta}_{2}, \ldots, \boldsymbol{\eta}_{m-1}, T\right)$. Так же, как в случае простой стрельбы, одну составляющую вектора $\eta_{1}$, а именно $\eta_{1 k}$, зафиксируем. Обозначим $\tilde{\eta}_{1}=\left(\eta_{1,1}, \ldots, \eta_{1, k-1}\right.$, $\eta_{1, k+1}, \ldots, \eta_{1 n}$ ). Тогда уравнения (5.8.41) можно рассматривать как систему $n(m-1)$ нелинейных уравнений относительно $n(m-1)$ неизвестных $\left(\tilde{\boldsymbol{\eta}}_{1}, \boldsymbol{\eta}_{2}, \ldots, \boldsymbol{\eta}_{m-1}, T\right)$ и одного параметра $\alpha$.

Для продолжения решения уравнений (5.8.41) в зависимости от параметра $\alpha$ можно вновь воспользоваться алгоритмом DERPER, который был описан в п. 5.8.4. Матрица Якоби

системы (5.8.41) имеет вид
\[
\left[\begin{array}{ccccccc}
\tilde{G}_{1} & -I & 0 & \ldots & & \partial F_{1} / \partial T & \partial F_{1} / \partial \alpha \\
0 & G_{2} & -I & \ldots & & \partial F_{2} / \partial T & \partial F_{2} / \partial \alpha \\
\cdots & \cdots & \cdots & \cdots & & \cdots & \cdots \\
& & 0 & G_{m-2} & -I & \cdots & \cdots \\
-\tilde{I} & 0 & & 0 & G_{m-1} & \partial F_{m-1} / \partial T & \partial F_{m-1} / \partial \boldsymbol{\alpha}
\end{array}\right]
\]

Здесь $G_{j}$ – матрицы размера $n \times n$,
$G_{j}=\left\{\frac{\partial x_{i}\left(z_{l+1} \mid \eta_{j}, T, \alpha\right)}{\partial \eta_{j p}}\right\}, i, p=1, \ldots, n, I$ – единичная матрица;

матрица $\sigma_{1}$ получается из матрицы $G_{1}$ вычеркиванием $k$-го столбца; аналогично из $I$ получается $I$.

Частные производные $\frac{\partial x}{\partial \eta}$ можно найти, используя уравнения в вариациях
\[
\frac{\partial U}{\partial z}=T \frac{\partial f(x(z), \alpha)}{\partial x} U,
\]

где $U$-матрица размера $n \times n$. Если выбрать начальное условие в виде $U\left(z_{j}\right)=I$, то, интегрируя уравнения (5.8.44) на интервале $\left[z_{j}, z_{j+1}\right]$, мы получим $G_{j}=U\left(z_{j+1}\right)$. Для частных производных $\partial F_{j} / \partial T$ имеют место соотношения
\[
\partial F_{j} / \partial T=\Delta z_{j} \cdot f\left(x\left(z_{j+1}\right) \mid \boldsymbol{\eta}_{j}, T, \boldsymbol{\alpha}\right) .
\]

Частные производные $\partial F_{j} / \partial \alpha$ находятся интегрированием уравнений в вариациях, аналогичных уравнениям (5.8.44).

В п. 5.8.3 мы говорили о структуре матрицы монодромии в случае метода простой стрельбы. В методе многократной стрельбы матрица монодромии имеет вид
\[
B=G_{m-1} G_{m-2} \ldots G_{2} G_{1},
\]

где матрицы $G_{j}, j=1,2, \ldots, m-1$, определяются формулой (5.8.43). Отметим, что столь простую структуру матрицы монодромии невозможно получить при использовании разностных методов.

Метод многократной стрельбы был нами использован для нахождения сильно неустойчивых периодических решений в так называемой модели Ходжкина-Хаксли передачи возмущения по нервному волокну [5.20], а также для продолжения по параметру периодических решений задачи 8 [5.29].

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