Материал данного параграфа основывается на сведениях, изложенных в п. 2.5.2, где были введены понятия аттрактора и хаотического инвариантного множества.

Хаотическое инвариантное множество, представляющее собой аттрактор, называется хаотическим (странным) аттрактором.

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

О существовании сложных траекторий в автономных системах обыкновенных дифференциальных уравнений с размерностью, больше или равной трем, было известно еще Пуанкаре. Лоренц в работе [4.40] на простом численном примере (задача 10 при значениях параметров $\sigma=10, b=8 / 3, r=28$ ) продемонстрировал наличие хаотического поведения решений.

Один из механизмов возникновения хаотического множества связан с последовательностью бифуркаций удвоения периода (или бифуркаций типа «-1», см. §5.8,2.3 и 2.5). Этой последовательности бифуркаций часто соответствует сходящаяся последовательность $\left\{\alpha_{n}\right\}_{n=1}^{\infty}$ бифуркационных значений параметра $\alpha$. Положим $\lim _{n \rightarrow \infty} \alpha_{n}=\alpha_{\infty}$. В случае если $\alpha_{n} \bigwedge \alpha_{\infty}$ (или соответственно $\alpha_{n} \backslash \alpha_{\infty}$ ) при $\alpha>\alpha_{\infty}$ (соответственно $\alpha<\alpha_{\infty}$ ) в фазовом пространстве системы $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x}, \boldsymbol{\alpha})$ существует хаотическое множество (см. рис. 5.34f).

Ниже мы рассмотрим некоторые методы анализа хаотических аттракторов. Для описания хаотического аттрактора мы будем использовать показатели Ляпунова (см. § 2.5) или будем изучать его структуру с помощью отображения Пуанкаре. Хаотический аттрактор характеризуется наличием хотя бы одного положительного одномерного показателя Ляпунова или тем, что инвариантное множество соответствующего отображе-

ния Пуанкаре имеет характер канторова множества 1). Более подробно вопросы возникновения и анализа хаотического пове-

Рис. 5.32. Фазовый портрет задачи 8, $N=2, A=2, B=5,9, D_{1}=1,21$, $D_{2}=12,1$. Часть хаотической траектории в проекции на плоскость $X_{1}-Y_{1}$.

дения систем с соответствующими приложениями рассматриваются в работе [5.18].

5.9.1. Вычисление показателей Ляпунова

Рассмотрим систему $n$ дифференциальных уравнений
\[
\frac{d \mathbf{x}}{d t}=\mathbf{f}(\mathbf{x}, \alpha), \quad\left(\mathbf{x} \in \mathrm{R}^{n}\right), \quad \alpha \in \mathrm{R}^{1} .
\]
1) Понимание смысла последнего высказывания для дальнейшего не обязательно. – Прим. ред.

В п. 2.5.1 мы ввели определение показателей Ляпунова. Напомним, что для траектории системы (5.9.1), которая в момент $t=0$ проходит через точку $\mathbf{x}_{0}$ (при фиксированном значении $\alpha$ ), одномерный показатель $\lambda^{1}$ определяется как
\[
\lambda^{1}\left(\mathbf{x}_{0} ; \mathbf{v}_{1}^{0}\right)=\lim _{t \rightarrow \infty} \frac{1}{t} \ln \frac{\left\|\mathbf{v}_{1}^{t}\right\|}{\left\|\mathbf{v}_{1}^{0}\right\|} .
\]

Соответствующий $k$-мерный показатель задается соотношением
\[
\lambda^{k}\left(\mathbf{x}_{0} ; \mathbf{v}_{1}^{0}, \ldots, \mathbf{v}_{k}^{0}\right)=\lim _{t \rightarrow \infty} \frac{1}{t} \ln \frac{\left\|\mathbf{v}_{1}^{t} \wedge \mathbf{v}_{2}^{t} \ldots \wedge \mathbf{v}_{k}^{t}\right\|}{\left\|\mathbf{v}_{1}^{0} \wedge \mathbf{v}_{2}^{0} \ldots \wedge \mathbf{v}_{k}^{0}\right\|} .
\]

Здесь $\left\|\mathbf{w}_{1} \wedge \mathbf{w}_{2} \ldots \wedge \mathbf{w}_{k}\right\|$ есть объем $k$-мерного параллелепипеда, образуемого векторами $\mathbf{w}_{1}, \ldots, \mathbf{w}_{k}$. В формулах (5.9.2) и (5.9.3) $\mathbf{v}_{i}^{t}$ есть решение линеаризованного уравнения
\[
\frac{d \mathbf{y}}{d t}=\frac{\partial \mathbf{f}}{\partial \mathbf{x}}(\varphi(t)) \mathbf{y},
\]

для которого $\mathbf{v}_{i}(0)=\mathbf{v}_{i}^{0}$ (матрица $\partial \mathrm{f} / \partial \mathbf{x}$ вычисляется в точках выбранной траектории $\mathbf{x}=\varphi(t)$ системы (5.9.1)).
Замечания
1. Предел в формуле (5.9.2) существует для почти всех $\mathbf{x}_{0}$ и векторов $\mathbf{v}_{1}^{0}$.
2. При случайном выборе линейно независимых векторов $\mathbf{v}_{1}^{0}, \ldots, \mathbf{v}_{k}^{0}$ формула (5.9.3) с вероятностью 1 дает максимальный $k$-мерный показатель $\lambda_{\max }^{k} \quad$ (см. теорему в п. 2.5.1) ${ }^{1)}$.

Вычисление одномерных показателей Ляпунова обычно проводится следующим образом. С помощью формулы (5.9.3) подсчитываем $n$ максимальных показателей $\lambda_{\max }^{1}, \lambda_{\max }^{2}, \ldots, \lambda_{\max }^{n}$ и затем находим одномерные показатели $\lambda_{1}^{1} \geqslant \lambda_{2}^{1} \geqslant \ldots \geqslant \lambda_{n}^{1}$ с помощью соотношений
\[
\lambda_{1}^{1}=\lambda_{\max }^{1}, \quad \lambda_{k}^{1}=\lambda_{\max }^{k}-\lambda_{\max }^{k-1}, \quad k=2, \ldots, n,
\]

поскольку каждый $k$-мерный показатель представляет собой сумму $k$ одномерных показателей (см. п. 2.5.1).

Если подсчитывать показатели Ляпунова с помощью выражения (5.9.3), то у нас могут возникнуть трудности вычислительного порядка. Так, например, если точка $\mathbf{x}_{0}$ лежит на хаотической траектории, то с ростом $t$ векторы $\mathbf{v}_{i}^{t}$ увеличиваются, а углы между ними уменьшаются, в результате чего вычисление

объема становится все более неточным. Чтобы преодолеть эти затруднения, заметим, что отношение объемов $v_{k}$ (дробь, стоящая в формуле (5.9.3)) не изменится, если вместо векторов $\mathbf{v}_{i}^{0}$ взятьих (независимые) линейные комбинации $\mathbf{w}_{i}$, поскольку $v_{k}$ зависит только от подпространства, натянутого на $\mathbf{v}_{i}^{0}$.

При вычислении показателей Ляпунова мы поступаем так. Задав небольшое $\tau$, вычисляем величину
\[
v_{k}(0 ; \tau)=\frac{\left\|\mathbf{v}_{1}^{\tau} \wedge \ldots \wedge \mathbf{v}_{k}^{\tau}\right\|}{\left\|\mathbf{v}_{1}^{0} \wedge \ldots \wedge \mathbf{v}_{k}^{0}\right\|}
\]
– «коэффициент изменения» $k$-мерного объема за время $0 \leqslant$ $\leqslant t \leqslant \tau$. Затем выбираем ортонормированный базис $\tilde{\mathbf{w}}_{i}(i=$ $=1, \ldots, k$ ) в подпространстве, натянутом на $\mathbf{v}_{i}^{\tau}$, например, ортогонализируя $\mathbf{v}_{i}^{\tau}$ по Граму-Шмидту. Затем на промежутке $\tau \leqslant t \leqslant 2 \tau$ вычисляем $v_{k}(\tau ; 2 \tau)$, используя решения $\mathbf{w}_{i}(t)$ линеаризованного уравнения (5.9.4) с начальными данными $\mathbf{w}_{i}(\tau)=\tilde{\mathbf{w}}_{i}$. Теперь коэффициент изменения объема за время $2 \tau$ есть, очевидно, произведение: $v_{k}(0 ; 2 \tau)=v_{k}(0 ; \tau) \cdot v_{k}(\tau ; 2 \tau)$. Проводя ортогонализацию и нормировку $\mathbf{w}_{i}$ в конце каждого интервала длины $\tau$, полагаем при $t=L \tau$
\[
\lambda_{\max }^{k} \approx \lambda^{k}(L)=\frac{1}{L \tau} \ln v_{k}(0 ; L \tau) .
\]

В реальных вычислениях целое число $L$ подбирается настолько большим, чтобы $\lambda^{k}(L+1)$ мало отличалось от $\lambda^{k}(L)$. Разумеется, для разных $k$ могут понадобиться разные $L$.

Мы ограничимся здесь двумя примерами. На рис. 5.33 изображена зависимость максимального показателя Ляпунова $\lambda_{\max }^{1}$ от параметра $\alpha$ для притягивающей траектории (в задаче 10). Учитывая большой объем необходимого машинного времени, вычисления проводились на относительно редкой сетке значений параметра $r$, выбранной на основе данных рис. 5.26. На рис. 5.33 можно видеть несколько интервалов изменения параметра (не-
Таблица 5.30. Одномерные показатели Ляпунова для задачи 8: $N=2, A=2$, $B=5,9, D_{1} / D_{2}=0,1$.

которые из них весьма малы), в которых существуют устойчивые периодические решения $\left(\lambda_{\max }^{1}=0\right.$ ) (см. рис. 5.26). При этом на интервалах, где $\lambda_{\max }>0$, аттрактор является хаотическим.

В табл. 5.30 приведены численные значения показателей Ляпунова для задачи о двух связанных реакторах (задача 8)

Рис. 5.33. Зависимость максимального одномерного піказателя Ляпунова $\lambda_{\max }^{1}$ от параметра $r$ в задаче $10, \sigma=16, b=4$.

при трех значениях параметра $D_{1}$. Положительная величина $\lambda_{1}^{1}$ указывает на существование хаотических траекторий; при этом $\lambda_{2}^{1}$ всегда равно нулю (см. замечание 2 в п. 2.5.1).

5.9.2. Отображение Пуанкаре

В п. 2.3.2 мы ввели определение отображения Пуанкаре, создаваемого потоком в окрестности замкнутой траектории $\gamma$ системы (5.9.1). Отображение Пуанкаре можно рассматривать как дискретную динамическую систему, определенную на сечении $\Sigma$. По характеру поведения орбит этой динамической системы можно судить о поведении траекторий системы (5.9.1). Аналогичный подход мы будем использовать для исследований траекторий на хаотическом аттракторе $\mathscr{A}$. Выберем подходящее сечение $\Sigma$ – участок гиперповерхности $S$ :
\[
S\left(x_{1}, \ldots, x_{n}\right)=0 .
\]

Лежащая в аттракторе траектория $\Gamma$, которая проходит через точку $\mathbf{x}_{0} \in \Sigma$, через определенное время пересекает сечение $\Sigma$ в точке $\mathbf{x}_{1}=\mathbf{P}\left(\mathbf{x}_{0}\right)$, затем в точке $\mathbf{x}_{2}=\mathbf{P}\left(\mathbf{x}_{1}\right)=\mathbf{P}^{2}\left(\mathbf{x}_{0}\right)$ и т. д. Таким образом, мы получаем орбиту $\mathbf{O}\left(\mathbf{x}_{0}\right)=\left\{\mathbf{x}_{0}, \mathbf{x}_{1}, \mathbf{x}_{2}, \ldots\right\}$ динамической системы, которая порождается отображением $\mathbf{P}$. По характеру орбиты $\mathbf{0}\left(\mathbf{x}_{0}\right)$ мы можем до некоторой степени судить о геометрической структуре аттрактора $\mathscr{A}$. Так, например, если орбита $\mathbf{O}\left(\mathbf{x}_{0}\right)$ плотно заполняет некоторую дугу $\mathbf{C}$, то дуга $\mathbf{C}$ является пересечением аттрактора $\mathscr{A}$ с поверхностью $\Sigma$, и можно думать, что аттрактор $\mathscr{A}$ является двумерным. В частности, если точки орбиты О плотно заполняют замкнутую кривую (типа окружности), то можно считать, что аттрактор имеет вид двумерного тора, который плотно заполняется одной траекторией системы (5.9.1).

Обычный численный подход к нахождению орбиты отображения Пуанкаре заключается в том, что мы интегрируем систему (5.9.1) и на каждом шаге интегрирования оцениваем знак функции $S$ в формуле (5.9.9). При изменении знака $S$ мы находим точку пересечения $\Gamma$ с поверхностью $\Sigma$ на основе интерполяции между двумя последними точками, найденными в процессе интегрирования. Очевидно, что грубая интерполяция будет источником численных погрешностей; для интерполяции более высокого порядка нам потребуется сохранять в ходе интегрирования большее число точек.

Нетрудно преобразовать процесс вычислений так, чтобы одна из точек, найденных численным интегрированием, оказалась прямо на $\Sigma$ (см. [5.21]).

Вместо общего соотношения (5.9.9) рассмотрим сначала случай, чаще всего используемый на практике, а именно $x_{i}$ $-a=0$, где $a$ – некоторая постоянная, а $i$ фиксировано $(1 \leqslant$ $\leqslant i \leqslant n)$.

Преобразуем систему (5.9.1) в эквивалентную систему, в которой роль независимой переменной вместо $t$ играет $x_{i}$. Разделим каждое уравнение системы (5.9.1) на $i$-е уравнение, причем будем предполагать, что $f_{i}(\mathbf{x}, \alpha)
eq 0$ в окрестности сечения $\Sigma$. Тогда наша система преобразуется к виду
\[
\begin{array}{l}
\frac{d x_{k}}{d x_{i}}=\frac{f_{k}(\mathbf{x}, \alpha)}{f_{i}(\mathbf{x}, \alpha)} ; \quad k=1, \ldots, n \quad(k
eq i) \\
\frac{d t}{d x_{i}}=\frac{1}{f_{i}(\mathbf{x}, \alpha)} .
\end{array}
\]

Теперь будем интегрировать систему (5.9.1) до момента изменения знака $x_{i}-a$. Далее, перейдем к интегрированию системы (5.9.9) с шагом $\Delta x_{i}=a-x_{i}$, причем начальные условия мы

будем брать в последней или в предпоследней точке, найденной интегрированием уравнения (5.9.1) ${ }^{1}$. Найденная таким образом точка лежит прямо на $\Sigma$ (с точностью до погрешностей аппроксимации метода интегрирования). Тем самым мы находим следующую точку орбиты отображения Пуанкаре и можем продолжить интегрирование системы (5.9.1). Удобно использовать при этом какой-нибудь одношаговый метод, например, метод Рунге-Кутты с автоматическим изменением длины шага. Отметим, что если нам не нужно знать моменты времени, когда траектория системы проходит через сечение $\Sigma$, то последнее уравнение в (5.9.9) можно опустить.

В общем случае гиперповерхности (5.9.8) мы вводим еще одну переменную
\[
x_{n+1}=S\left(x_{1}, \ldots, x_{n}\right)
\]

и добавляем к системе (5.9.1) дифференциальное уравнение вида
\[
\frac{d x_{n+1}}{d t}=f_{n+1}(\mathbf{x}, \alpha),
\]

где
\[
f_{n+1}=\sum_{j=1}^{n} f_{j}(\mathbf{x}, \alpha) \frac{\partial S}{\partial x_{j}} .
\]

Тем самым мы получаем новую систему из $n+1$ дифференциальных уравнений, причем начальное условие для неизвестной $x_{n+1}$ задается в соответствии с формулой (5.9.10).

Соотношение (5.9.8), описывающее гиперповерхность, имеет теперь вид $x_{n+1}=0$. Далее мы поступаем точно так же, как это делалось ранее для плоскости $x_{i}-a=0$, полагая $i=n+1$.

На рис. $5.34 \mathrm{a}$ – $\mathrm{f}$ приведены несколько периодических и одна хаотическая орбита отображения Пуанкаре для задачи 10. Гиперплоскость $\Sigma$ определялась при этом уравнением $y=0$. На рис. 5.34а изображена двухточечная орбита. Эта орбита возникла после бифуркации удвоения периода от основной ветви устойчивых периодических решений (см. рис. 5.26 из §5.8). На рисунках $5.34 \mathrm{~b}$, c, d, е приведены орбиты отображения Пуанкаре, отвечающие периодическим решениям задачи 10 с периодами 4T (возникающим после двух бифуркаций удвоения периода. – Ред.), $8 T, 16 T$ для разных значений параметра $r$. На этих рисунках хорошо прослеживается эволюция, которую претерпевает орбита рис. 5.34а в ходе последовательных бифур-

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

Рис. 5.34. Орбиты отображения Пуанкаре. Задача $10(a-f), \sigma=16, b=4$, $y=0$. a) $r=339,0, b) r=338,0$, c) $r=334,5, d) r=334,2, e) r=333,25$, f) $r=332,5$. Задача $8, N=2(g, h), A=2, B=5,9, D_{2}=10 D_{1}$; показана проекция орбиты с гиперплоскости $X_{1}-Y_{1}+X_{2}-Y_{2}+0,9=0$ на плоскость $X_{1}-X_{2}$ g) $\left.D_{1}=1,194, h\right) D_{1}=1,21$.

рис. 5.34f, где изображена «хаотическая» орбита отображения Пуанкаре для значения параметра $r=332,5$. Множество точек пересечения траектории системы (5.9.1) с гиперплоскостью $\Sigma$ плотно заполняет дугу кривой, которая входит в пересечение

(общую часть) хаотического аттрактора с гиперплоскостью $\Sigma$.
Значения параметра $r$ в бифуркационных точках удвоения периода образуют так называемую последовательность Фейгенбаума, стремящуюся к некоторому пределу (см. п. 2.5.3). Эти значения параметра $r$ можно подсчитать с помощью алгоритмов, описанных в п. 5.8.5, и на их основе вычислить величины
\[
\delta_{j}=\frac{r_{j}-r_{j-1}}{r_{j+1}-r_{j}} .
\]

Координаты точек бифуркации и значения $\delta_{i}$ приведены в табл. 5.31 [5.26]. Из таблицы видно, что найденные значения $\delta_{j}$ стремятся к пределу $\delta^{*} \sim 4,6692$ [5.27].

Таблица 5.31. Каскад бифуркационных точек с удвоением периода для модели Лоренца, задача $10\left(\sigma=16, b=4, k=1, x_{k}=3,82038\right)$.

На рис. $5.34 \mathrm{~g}$, h представлены две хаотические орбиты отображения Пуанкаре для задачи о двух связанных между собой реакторах (задача 8). Трехмерная гиперплоскость $\Sigma$ задается при этом соотношением $X_{1}-Y_{1}+X_{2}-Y_{2}+0,9=0$. Для более наглядного геометрического представления орбита $\mathbf{O}$ отображения Пуанкаре, лежащая в плоскости $\Sigma$, обычно проектируется на какую-либо из координатных плоскостей. В данном случае орбита 0 спроектирована на плоскость $X_{1}-X_{2}$ (см. рис. $5.34 \mathrm{~g}, \mathrm{~h}$ ). Из рисунка ясно, что точки орбиты 0 располагаются на некоторой гладкой кривой (точнее, в ее очень малой плоскости). Показатели Ляпунова, которые подсчитывались с помощью подхода, описанного в п. 5.9.1, приведены в табл. 5.30 для значений параметров, указанных на рис. $5.34 \mathrm{~g}$, h. При этом положительное значение $\lambda_{1}^{1}$ указывает на наличие хаотического аттрактора.