Напомним (см. рис. 5.2) определение характеристического показателя Ляпунова б для траектории $\boldsymbol{x}(t)$ и близкой к ней $\boldsymbol{x}+\boldsymbol{w}$, где $w$ – касательный вектор:
\[
\begin{array}{c}
\sigma\left(x_{0}, w_{0}\right)= \\
=\lim _{\substack{t \rightarrow \infty \\
d(0) \rightarrow 0}} \frac{1}{t} \ln \frac{d(t)}{d(0)} .
\end{array}
\]

В системах, близких к интегрируемым, аналитический и численный расчет показателей $\sigma$ (особенно максимального $\sigma=\sigma_{1}$ ) широко используется в качестве критерия стохастичности. Отметим прежде всего [57], что в интегрируемых гамильтоновых системах все $\sigma$ равны нулю. Рассмотрим в качестве примера движение на торе
Рис. 5.5. Линейный рост расстояния $d(t)$ между близкими траекториями в интегрируемой системе на примере отображения поворота.
\[
\begin{array}{c}
\sigma\left(x_{0}, w_{0}\right)= \\
=\lim _{\substack{t \rightarrow \infty \\
d(0) \rightarrow 0}} \frac{1}{t} \ln \frac{d(t)}{d(0)} .
\end{array}
\]
[в проекции на плоскость $(J, \theta)$ см. рис. 3.1, $a$ ]:
\[
\begin{array}{l}
J(t)=J_{0}, \\
\theta(t)=\theta_{0} \div \omega(J) t .
\end{array}
\]

Как показано на рис. 5.5, траектории в этом случае представляют собой окружности, а расходимость близких траекторий оказывается максимальной, когда вектор $w=(\Delta J, 0)$, т. е. направлен вдоль радиуса. Для двух близких траекторий имеем
\[
\begin{array}{l}
\theta_{1}(t)=\theta_{0}+\omega_{0} t, \\
\theta_{2}(t)=\theta_{0}+\omega_{0} t-\omega_{0}^{\prime} d_{0} t,
\end{array}
\]

где $\omega_{0}^{\prime}$ – производная по $J$. Расстояние между траекториями равно
\[
d^{2}(t)=J_{0}^{2}\left(\theta_{2}-\theta_{1}\right)^{2}+d_{0}^{2}=d_{0}^{2}\left[\left(\omega_{0}^{\prime} J_{0} t\right)^{2}+1\right]
\]

и при больших $t$ растет линейно со временем. Тогда из (5.2.8) имеем
\[
\sigma_{1}=\lim _{t \rightarrow \infty} \frac{1}{2 t} \ln \left[\left(\omega_{0}^{\prime} J_{0} t\right)^{2}-1\right]=0 .
\]

С другой стороны, в хаотической области
\[
\begin{array}{c}
d(t) \sim d_{0} e^{\sigma_{1}(x) t}, \\
\lim _{t \rightarrow \infty} \frac{1}{t} \ln \frac{d}{d_{0}}=\sigma_{1}>0 .
\end{array}
\]
5.3а. Аналитические оценки

Для достаточно простой динамической системы, в которой области устойчивости занимают пренебрежимо малую часть фазового пространства, максимальный показатель Ляпунова $\sigma_{1}$ можно получить аналитически. Для стандартного отображения при больших $K$ это было сделано Чириковым [70]. Линеаризовав стандартное отображение вблизи некоторой точки, не обязательно неподвижной, получим из характеристического уравнения при больших $K$ наибольшее из двух собственных значений в виде
\[
\lambda^{+}=K|\cos \theta| .
\]

Так как вся фазовая плоскость считается стохастической, то величина $\sigma_{1}$ (и $h$ ) определяется просто усреднением $\ln \lambda^{+}$. В данном случае это эквивалентно усреднению по $\theta$ и дает
\[
\sigma_{1}=\frac{1}{2 \pi} \int_{0}^{2 \pi} \ln |K \cos \theta| d \theta=\ln \frac{K}{2} .
\]

Этот аналитический результат сравнивался с численным значением $\sigma_{n}$, полученным для начальных условий на стохастической компоненте. При $K \approx 6$, когда еще можно ожидать количественного согласия с (5.3.8), Чириков получил ${ }^{1}$ ) $\sigma_{n} / \sigma_{1}=1,02$.

Как мы видели в § 4.1, в стандартном отображении островки устойчивости существуют при сколь угодно больших $K$. С увеличением $K$ эти островки уменьшаются, но поскольку они могут существовать вокруг периодических точек произвольно большого периода, то возникает вопрос, не занимают ли они значительную часть фазовой плоскости, что привело бы к неправильному значению КС-энтропии (5.3.8) даже для больших $K$. Чириков исследовал этот вопрос двумя методами. В первом квадрат $2 \pi \times 2 \pi$ разделялся на $100 \times 100$ ячеек и вычислялась доля $\mu_{s}$ ячеек, в которые попадала траектория с начальными условиями на стохастической компоненте. Ясно, что такое «огрубление» может давать правильные результаты только для относительно малых значений $K$, когда
1) Более полные аналитические и численные данные по КС-энтропии для различных отображений приведены в работе [73].- Прим. ред.

островки устойчивости достаточно велики. Во втором методе выбиралось некоторое число, скажем $N=100$, траекторий со случайными начальными условиями и определялась доля $\mu$, тех из них, для которых численное значение КС-энтропии близко к нулю. Естественно ожидать, что $\mu_{s}+\mu_{r} \approx 1$. Результаты этих численных экспериментов представлены в табл. 5.1. Чириков пришел к заключению, что мера устойчивых областей с увеличением $K$ стремится к нулю ${ }^{1}$ ), так что при больших $K$ возможно простое статистическое описание движения.

Таблица 5.1. Численное определение меры стохастической $\left(\mu_{\mathrm{s}}\right)$ и регулярной ( $\left.\mu_{r}\right)$ компонент стандартного отображения (по данным работы [70], табл. 5.3)
5.3б. Численные методы

Удобным практическим критерием стохастичности в заданной области фазового пространства служит численное определение показателей Ляпунова. Максимальный показатель $\sigma_{1}$ широко использовался в качестве критерия стохастичности в ряде работ (см., например, [19, 41-43, 73, 133, 371$]$ ).
Вычисление $\sigma_{1}$. Ранее было показано (п. 5.2б), что $\boldsymbol{\sigma}\left(\boldsymbol{x}, \boldsymbol{w}_{0}\right)=$ $=\sigma_{1}(\boldsymbol{x})$ для почти всех касательных векторов $\boldsymbol{w}_{0}$. Поэтому для вычисления $\sigma_{1}$ начальный вектор $w_{0}$ можно выбирать произвольно. Интегрируя совместно уравнения (5.2.6) и (5.2.7), находим величину
\[
d(t)=|w(t)|,
\]

где для удобства принято $d_{0}=d(0)=1$. Трудности возникают в том случае, когда $|w|$ растет экспоненциально, что приводит к переполнению и другим численным ошибкам. Чтобы избавиться от этих неприятностей, выберем определенный интервал времени $\tau$, как показано ${ }^{2}$ ) на рис. 5.6 , и будем перенормировывать $\qquad$
1) Согласно данным работы [73], $\mu_{r} \sim \exp (-C \sqrt{\bar{K}})$, где $C \sim$ const, за исключением специальных значений $K$, близких к (4.1.14).- Прим. ред.
2) Фактически на рис. 5.6 показан другой метод вычисления $\sigma_{1}$ – по двум близким траекториям системы (5.2.6), а не по одной траектории и линеаризованному уравнению (5.2.7), как в тексте. Использование линеаризованного уравнения имеет то преимущество, что его решение не зависит от модуля |w|. Однако в некоторых специальных случаях удобнее интегрировать две траектории (см., например, [470]).- Прим. ред.

$|w|$ на единицу в конце каждого такого интервала ${ }^{1}$ ). Таким способом будем последовательно вычислять величины
\[
\begin{aligned}
d_{k} & =\left|w_{k-1}(\tau)\right|, \\
w_{k}(0) & =\frac{w_{k-1}(\tau)}{d_{k}},
\end{aligned}
\]

где $w_{k}(\tau)$ получается интегрированием (5.2.6), (5.2.7) на интервале $\tau$ с начальными значениями $\boldsymbol{x}(k \tau)$, $\boldsymbol{w}_{k}(0)$. Если ввести величину
\[
\sigma_{n}=\frac{1}{n \tau} \sum_{i=1}^{n} \ln d_{i},
\]

Рис. 5.6. Схема вычисления максимального показателя Ляпунова (по данным работы [19]).
$y=x+w, \tau-$ конечный интервал времени.

то из (5.2.8) получаем
\[
\sigma_{1}=\lim _{n \rightarrow \infty} \sigma_{n}
\]

Для регулярного движения $\sigma_{1}=0$, в то время как на стохастической компоненте $\sigma_{1}>0$ и не зависит от $\boldsymbol{x}$. Описанный метод применим как для непрерывных траекторий, так и для отображений.

Следуя работе Бенеттина и др. [19], проиллюстрируем использование этого метода на примере модели Хенона-Хейлеса (п. 1.4а). На рис. 5.7 приведена зависимость $\sigma_{n}$ от $n \tau$ для шести траекторий с энергией $E=0,125$ (см. рис. 1.13). Три траектории (1-3) находятся, очевидно, в области устойчивого движения, тогда как три другие (4-6) – в широкой стохастической области, которая хорошо видна на рис. 1.13. Как и следовало ожидать, для первых траекторий с увеличением $n$ величина $\sigma_{n}$, очевидно, стремится к $\qquad$
1) То есть сокращать длину вектора $\boldsymbol{w}\left(|\boldsymbol{w}| \rightarrow d_{0}=1\right)$, не изменяя его направления (см. [73]).- Прим. ред.

нулю $^{\mathbf{1}}$ ), а для вторых, по-видимому,-к одному и тому же значению $\left.{ }^{2}\right) \sigma_{1}$.

Используя этот метод, Бенеттин с соавторами нашли зависимость $\sigma_{1}$ от энергии $E$ в областях стохастического и регулярного движения (рис. 5.8). В последнем случае $\sigma_{n} \rightarrow 0$ (черточки на оси абсцисс). Сплошная кривая соответствует подгонке значений $\sigma_{1}>0$ к экспоненциальной зависимости. Используя соотно-

Рис. 5.7. Зависимость показателя $\sigma_{n}$ от времени для траекторий на регулярной ( $1-3$ ) и стохастической (4-6) компонентах (по данным работы [19]),

шение (5.2.27) и зависимость $1-\mu_{s}$ от $E$, полученную Хеноном и Хейлесом [188], Бенеттин с соавторами вычислили КС-энтропию $h(E)$ (пунктирная кривая на рис. 5.8). Они высказывают предположение, что $h(E)>0$ для всех $E>0$, что согласуется с нашими общими представлениями о свойствах движения, хотя и не видно из рис. 5.8 из-за различных приближений при вычислении $\sigma_{1}$. Это указывает на некоторый присущий методу недостаток. Возможно, например, что одна из траекторий на рис. 5.7, которая кажется регулярной, на самом деле лежит в узком стохастическом слое. Такую траекторию одинаково трудно обнаружить как непосредственно, так и при помощи настоящего метода, поскольку выход на очень низкое плато $\sigma_{1}$ может оказаться за пределами фактической $\qquad$
1) Грубо говоря, как $1 / n \tau$ [см. (5.3.4)].- Прим. ред.
2) Вообще говоря, это не всегда так (см., например, [73], рис. 3). В данном примере близкие значения $\sigma_{n}$ указывают на принадлежность всех трех траекторий к одной и той же стохастической компоненте.- Прим. ред.

длительности счета ${ }^{1}$ ). Мы вернемся к этим вопросам в гл. 6 при обсуждении более сложного движения в многомерных системах.

Вычисление всех показателей Ляпунова ${ }^{2}$ ). Следуя Бенеттину и др. $[20]^{3}$ ), покажем, каким образом вычисляется полный набор показателей Ляпунова в $M$-мерном фазовом пространстве. Ясно, что любая попытка определить $\sigma_{2}, \sigma_{3}$ и т. д., выбирая касательный вектор $w$ вдоль векторов $\boldsymbol{e}_{2}, \boldsymbol{e}_{3}$ и т. д. (см. рис. 5.2), обречена на неудачу из-за неустойчивости этих направлений, так что любые ошибки повернут в конце концов $w(t)$ вдоль $e_{1}$. Вместо этого выберем начальный базис из $p$ ортонормированных касательных векторов и численно определим $p$-мерный объем $V_{p}(t)$, заданный этими векторами. Отсюда можно найти показатель Ляпунова $\sigma_{1}^{(p)}$ порядка $p$ (5.2.14). Проделав эту процедуру для $p=1,2, \ldots, M$, из (5.2.15) определим все показатели $\sigma_{1}$, $\sigma_{2}, \ldots, \sigma_{M}$. Здесь, однако, возникает следующая трудность. В процессе движения углы между касательными векторами, вообще говоря, экспоненциально уменьшаются и численные ошибки резко возрастают. Поэтому в дополнение к перенормировке длины векторов $\boldsymbol{w}$ необходимо также периодически ортогонализовать их. При этом новые векторы !должны стве, что и старые.

Оказывается, что такую процедуру одновременного вычисления всех $p$-мерных объемов можно свести к расчету эволюции $M$ векторов и их специальной ортогонализации по методу Грама-
1) Заметим, что хотя действительно таким методом нельзя установить регулярность движения (интегрируемость) в некоторой области фазового пространства, однако вполне можно обнаружить хаотическую (случайную) компоненту движения, т. е. неинтегрируемость. В некотором смысле этот численный метод дополняет аналитический метод обратной задачи рассеяния (см. примечание редактора на с. 56), который позволяет, наоборот, доказать только интегрируемость.- Прим. ред.
2) Заметим, что для критерия стохастичности нужен по существу только максимальный показатель $\sigma_{1}$ ввиду неравенства $h \geqslant \sigma_{1}>0$. Иначе говоря, точное значение КС-энтропии $h$, которое требует знания всех показателей, несущественно.- Прим. ред.
$\left.{ }^{3}\right)$ См. также [502].- Прим. ред.

Шмидта. Вспомним, что $w_{k-1}(\tau)$ – касательный вектор $w_{k-1}(0)$ через время $\tau$. Вычислим сначала для каждого временного интервала $\tau$ следующие величины:
\[
\begin{aligned}
d_{k}^{(1)} & =\left|w_{k-1}^{(1)}(\tau)\right|, \\
w_{k}^{(1)}(0) & =\frac{w_{k-1}^{(1)}(\tau)}{d_{k}^{(1)}} .
\end{aligned}
\]

Pис. 5.9. Все три положительных показателя Ляпунова для модели Контопулоса (см. текст) (промежуточная область) (по данным работы [20]).

Затем для $j=2$, . ., $M$ найдем последовательно величины
\[
\begin{aligned}
\boldsymbol{u}_{k-1}^{(i)}(\tau) & =\boldsymbol{w}_{k-1}^{(j)}(\tau)-\sum_{i=1}^{j-1}\left(\boldsymbol{w}_{k}^{(i)}(0) \cdot \boldsymbol{w}_{k-1}^{(j)}(\tau)\right) \boldsymbol{w}_{k}^{(i)}(0),(5.3 .12 \text { в }) \\
d_{k}^{(j)} & =\left|\boldsymbol{u}_{k-1}^{(j)}(\tau)\right|, \\
\boldsymbol{w}_{k}^{(j)}(0) & =\frac{\boldsymbol{u}_{k-1}^{(i)}(\tau)}{d_{k}^{(j)}} .
\end{aligned}
\]

Тогда в течение ( $k-1$ )-го интервала времени $\tau$ объем $V_{p}$ возрастает в $d_{k}^{(1)} d_{k}^{(2)} \ldots d_{k}^{(p)}$ раз. Отсюда, согласно (5.2.14), получаем
\[
\sigma_{1}^{(p)}=\lim _{n \rightarrow \infty} \frac{1}{n \tau} \sum_{i=1}^{n} \ln \left(d_{i}^{(1)} d_{i}^{(2)} \ldots d_{i}^{(p)}\right) .
\]

Вычитая $\sigma_{1}^{(p-1)}$ из $\sigma_{1}^{(p)}$ и используя (5.2.15), находим $p$-й показатель Ляпунова
\[
\sigma_{p}=\lim _{n \rightarrow \infty} \frac{1}{n \tau} \sum_{i=1}^{n} \ln d_{i}^{(p)} .
\]

Это соотношение фактически и используется при численном моделировании.

Бенеттин и др. [20] нашли таким способом все показатели Ляпунова для нескольких гамильтоновых систем, включая 4- и 6-мерные отображения. Мы приведем их результаты для системы с тремя степенями свободы, которая исследовалась Контопулосом и др. [93]:
\[
H=\frac{1}{2}\left(p_{1}^{2}+q_{1}^{2}\right)+\frac{\sqrt{2}}{2}\left(p_{2}^{2} \div q_{2}^{2}\right)+\frac{\sqrt{3}}{2}\left(p_{3}^{2}+q_{3}^{2}\right)+q_{1}^{2} q_{2}+q_{1}^{2} q_{3} .
\]

Движение этой системы ограничено при $H<0,097$. Для $H=0,09$ фазовое пространство разделяется, грубо говоря, на три области: большая область стохастичности с $\sigma_{1} \approx 0,03 ; \sigma_{2} \approx 0,008$ и $\sigma_{3} \approx 0$; область регулярного движения ( $\sigma_{1}=\sigma_{2}=\sigma_{3}=0$ ) и промежуточная область. На рис. 5.9 представлены результаты вычисления первых трех показателей Ляпунова для начальных условий в промежуточной области. Видно, что сходимость имеет место для $\sigma_{1} \approx 3 \times 10^{-3}$ и $\sigma_{3}=0$, соответствующего направлению вдоль траектории. Поведение показателя $\sigma_{2}$ не вполне ясно.

Қак мы увидим в § 6.2, эти результаты на самом деле обманчивы. Действительно, в системе с тремя степенями свободы первая и третья области должны быть связаны слабой диффузией Арнольда, благодаря которой траектория переходит из одной области в другую. Поэтому, по-видимому, и для промежуточной области $\sigma_{1} \approx 0,03$, а $\sigma_{2} \approx 0,008$, что противоречит данным на рис. 5.9. Это еще раз указывает на основную трудность численного определения показателей Ляпунова: не существует априорного условия для определения достаточного числа итераций $n$. Поэтому при численном моделировании необходимо использовать и другие методы, такие, например, как метод сечения Пуанкаре ${ }^{1}$ ).