Лабораторная работа No. 7
практикума "Анализ временных рядов":
Расчет старшего показателя Ляпунова
по временному ряду

Задачи лабораторной работы:

$\bullet$ Ознакомление с методом вычисления старшего ляпуновского показателя по скалярному временному ряду.

$\bullet$ Изучение влияния параметров алгоритма на результат вычислений, оценка точности расчетов для случая регулярной динамики.

$\bullet$ Расчет старшего показателя Ляпунова хаотического режима колебаний предложенной модельной системы.

Теоретические сведения:

Интерпретация ляпуновских показателей

В данной работе будет рассмотрен наиболее часто используемый на практике метод расчета ляпуновских показателей по временному ряду, предложенный в статье Вольфа с соавторами (Physica D, 1985). Этот метод немного напоминает стандартный подход к вычислению спектра ЛХП потоковых систем по известной математической модели (метод Бенеттина). В целом, техника расчета ляпуновских показателей базируется на следующей идеологии.

Предположим, что задана $n$-мерная динамическая система с непрерывным временем. Чтобы охарактеризовать устойчивость ее решения, анализируется временная эволюция бесконечно малой $n$-мерной сферы начальных условий; с течением времени эта сфера преобразуется в эллипсоид. Если говорить о спектре ЛХП, то $i$-й показатель Ляпунова может быть определен в терминах длин осей эллипсоида $p_i(t)$:

$$\lambda_i = \lim_{t \rightarrow \infty} {1 \over t} \ln {p_i (t) \over p_i (0)},$$

где $\lambda_i$ упорядочены от наибольшего к наименьшему. Таким образом, показатели Ляпунова определяются расширением либо сжатием сферы по различным направлениям в фазовом пространстве. Поскольку ориентация эллипсоида непрерывно меняется со временем, меняются и направления, ассоциирующиеся с тем или иным показателем. Поэтому нельзя говорить о каком-то одном направлении для каждой ляпуновской экспоненты.

Отметим, что если задается только одно возмущение (это значит, что мы следим лишь за главной осью эллипсоида), то оно будет в линейном приближении расти по закону $e^{\lambda_1 t}$. Для двух независимых возмущений площадь образуемого ими квадрата (или прямоугольника) эволюционирует по закону $e^{(\lambda_1 + \lambda_2) t}$; для трех - эволюция объема описывается законом $e^{(\lambda_1 + \lambda_2 + \lambda_3) t}$ и т.д. Данное свойство приводит к несколько иному определению спектра ляпуновских экспонент: сумма первых $j$ показателей определяется скоростью экспоненциального роста $j$-мерного элемента объема. Такая интерпретация обеспечивает основу для техники анализа экспериментальных данных (когда требуется вычислить несколько показателей).

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

Иногда ляпуновские показатели анализируют в терминах информации, в частности, как скорость потери информации о начальном состоянии. Поэтому их порой определяют в таких единицах как bits/s (количество бит информации в секунду), bits/period (для потоковых систем) или bits/iteration (для систем с дискретным временем). Соответственно, вычисляют скорость разбегания соседних траекторий не по закону $e^{\lambda_1 t}$, а по закону $2^{\lambda_1 t}$. Например, если старший показатель Ляпунова аттрактора Лоренца равен 2.16 bits в единицу времени, а начальные условия определены с точностью до одной миллионной части ($\approx 20$ bits), то дальнейшее поведение не может быть предсказано спустя $\approx 9$ единиц безразмерного времени (20/2.16), что соответствует $\approx$ 20 условным периодам (оборотам вокруг одного из состояний равновесия). Спустя это время малая начальная неопределенность будет покрывать весь аттрактор, т.е. фазовая точка может оказаться где угодно.

Поскольку расчет ляпуновских экспонент по экспериментальным данным имеет одну общую идеологию с расчетом показателей по системе ОДУ, кратко проанализируем эту процедуру.

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

Как известно, для режима динамического хаоса характерно наличие экспоненциальной неустойчивости траекторий, количественной мерой которой является положительный ляпуновский показатель, характеризующий степень чувствительности системы к выбору начальных условий. Число положительных экспонент в спектре ЛХП определяется количеством неустойчивых направлений периодических орбит, встроенных в хаотических аттрактор, хотя в принципе возможны и более сложные ситуации (состоящие в сосуществовании периодических орбит с различным числом неустойчивых направлений). Мы ограничимся рассмотрением динамических систем с одним положительным показателем Ляпунова $\lambda_1$. При вычислении $\lambda_1$ делается предположение о "типичности" фазовой траектории, являющейся решением ДС при выбранных начальных условиях. В противном случае величина показателя, определенная на интервале времени $T$, может отличаться от предельного значения $\lambda_1$ (соответствующего $T \rightarrow \infty$). В частности, известны случаи, когда дискретные модели демонстрируют очень длительные переходные процессы (до 1,500,000 итераций динамика является "хаотической", после чего сменяется регулярной). С точки зрения вычисления ляпуновского показателя это соответствует тому, что величина $\lambda_1$ в течение переходного процесса сходится к некоторому положительному значению, и только на очень длительных временах спадает до нуля. Чтобы охарактеризовать поведение типичной фазовой траектории, иногда используют понятие "ляпуновских показателей на конечном времени" (finite time Lyapunov exponents), которые характеризуют скорость разбегания или сжатия по различным направлениям в течение конечного интервала времени $T$.

Если уравнения ДС, генерирующей фазовую траекторию, известны, то определить величину максимального показателя (или полный спектр ЛХП) можно с помощью хорошо известного алгоритма, предложенного в работе Бенеттина с соавторами (независимо похожую технику предложили Шимада и Нагашима). Данный метод часто называют стандартным алгоритмом расчета ляпуновских экспонент.

Рассмотрим систему ОДУ:

$${d \vec{x} \over dt} = \vec{f} (\vec{x}, \vec{\mu}), \hspace{1cm} \vec{x} \in R^n, \vec{\mu} \in R^m,$$

в которой $\vec{x}$ - вектор состояния, $\vec{f}$ - нелинейная вектор-функция, $\vec{\mu}$ - вектор управляющих параметров. Исследование данной системы на устойчивость ее частного решения $\vec{x}_p(t)$ сводится к анализу уравнений в вариациях и введению понятия $k$-мерного ляпуновского показателя. Однако, если ограничиться только вычислением старшего ЛХП, данный алгоритм существенно упрощается. Искомая величина $\lambda_1$ будет определять эволюцию во времени вектора возмущения

$$r(t) = r_0 e^{\lambda_1 t},$$

где $r_0$ - величина начального возмущения (в момент времени $t_0=0$), $r_0=\mid \vec{x}(t_0) - \vec{x}_p(t_0) \mid$. Данная формула является очень приближенной, поскольку скорость разбегания траекторий не является постоянной, а зависит от выбора точки на аттракторе. Строго говоря, ляпуновские показатели определяют путем решения уравнений в вариациях, рассматривая бесконечно малые возмущения. Однако существуют работы, в которых обосновано, что экспоненциальный закон разбегания траекторий справедлив не только для бесконечно малых возмущений, но и для малых возмущений конечной величины, поэтому последняя формула вполне может быть использована при практическом анализе локальной неустойчивости. Поскольку большие значения $r(t)$ принимать не может (в противном случае не будет выполняться условие линейного приближения), при расчете старшего ЛХП проводят перенормировки, в ходе которых задаются новые возмущения - малые по величине, но выбранные в направлении, которое соответствует вектору возмущения непосредственно перед перенормировкой. Затем вновь оценивают скорость экспоненциального разбегания близлежащих траекторий. В результате вычисляется усредненная вдоль фазовой траектории $\vec{x}_p(t)$ количественная характеристика степени хаотичности. Схематично процедура перенормировки изображена на рис. 1.

Расчет старшего ЛХП по одномерной реализации (временному ряду)

Метод расчета $\lambda_1$ по реализации использует ту же идеологию, но теперь возникает одна проблема: мы не можем произвольным образом задавать вектор возмущения. Поскольку анализируется экспоненциальная неустойчивость траекторий в фазовом пространстве, на первом этапе процедуры расчета показателя нужно осуществить реконструкцию траектории в фазовом пространстве, т.е. сформировать множество векторов $\vec{z}(t)$ (см. работу 6).

Мы не можем выбирать вектор возмущения при перенормировках строго в заданном направлении, и должны брать новую точку там, где мы ее можем найти (в некотором конусе) - рис. 2. Невозможность выбора вектора возмущения строго в заданном направлении приводит к появлению ошибки ориентации (угол $\alpha$). Поскольку работа с одной единственной траекторией ограничивает возможности выбора вектора возмущения, на практике приходится искать компромисс между уменьшением ошибки его ориентации в фазовом пространстве и минимизацией его длины. Ограничения на длину могут быть сформулированы следующим образом:

$$l_1 < r(t) < l_2.$$

Необходимо вводить минимальное значение $l_1$, поскольку при $r(t) < l_1$ сказывается влияние шума, нарушающего экспоненциальный характер разбегания траекторий. Величина $l_2$ задает условие линейного приближения и часто может быть введена в процентном отношении от размера аттрактора (скажем, 5-10%). По аналогии со стандартным алгоритмом проводятся перенормировки вектора возмущения, если расстояние между траекториями перестает удовлетворять условию линейного приближения ($r(t) > l_2$). Поскольку вычисление $\lambda_1$ предполагает реконструкцию аттрактора, результат расчета данной величины будет зависеть от качества реконструкции, что приводит к появлению дополнительных параметров численной схемы - размерности пространства вложения, задержки и т.п.

Несколько слов о перенормировках вектора. Здесь возможны варианты:

• перенормировки проводятся после достижения определенного расстояния между траекториями ($r(t) > l_2$);
• перенормировки проводятся через фиксированные интервалы времени (например, через характерный период колебаний).

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

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

Рисунок 1

Рисунок 2

color=red> Ссылки на технические руководства:

$\bullet$ Как написать, откомпилировать и запустить простую С-программу
$\bullet$ Как выводить результаты расчета на экран в виде графиков
$\bullet$ Руководство по языку C
$\bullet$ Использование программы gnuplot для оперативного вывода графической информации.

Порядок выполнения лабораторной работы:

1) Ознакомление с теоретическим материалом.
2) Расчет старшего показателя Ляпунова периодического режима колебаний и оценка точности вычислений. Исследование влияния параметров алгоритма (задержки, величины вектора возмущения и т.д.) на погрешность оценки показателя (задание 1).
3) Расчет старшего ляпуновского показателя хаотической динамики по реализации предложенной преподавателем модели (задание 2).
4) Оформление отчета по полученным результатам.

Содержание и оформление отчета по лабораторной работе:

Отчет по лабораторной работе представляется в виде latex или html документа по указанию преподавателя. Он должен содержать:

1) Заголовок, с указанием названия лабораторной работы, Ф.И.О. выполнявших ее студентов, номер учебной группы, номер задания.

2) Сведения (уравнения и комментарий) о модельной системе, смысл параметров, диапазон их изменения.

3) Построенную двупараметрическую бифуркационную диаграмму.

4) Информацию о количестве и типе аттракторов либо особых точек/траекторий в исследованной области и о смысле бифуркационных линий

5) Наброски фазовых портретов либо копии экрана, характеризующие поведение исследуемой системы в различных областях управляющих параметров.

6) Краткое резюме - заключение по пунктам 1-5.

Список вариантов заданий на выполнение лабораторной работы:

Задание 1. С помощью программы ``lhp.x'' провести расчеты старшего ляпуновского показателя периодического режима колебаний. Построить зависимости величины $\lambda_1$ от:
(а) времени задержки (при фиксированных остальных параметрах);
(б) размерности пространства вложения;
(в) минимального расстояния между точками в фазовом пространстве ($l_1$);
(г) максимального расстояния между точками в фазовом пространстве ($l_2$);
(д) времени между перенормировками.

N варианта	динамическая система	примечание
1	система Ресслера	-
2	система Лоренца	-
3	генератор с инерционной нелинейностью	-
4	модель нефрона	-
5	модель нейрона (Хиндмарш-Розе)	-
6	модель бета-клетки	-

Задание 2. Провести те же самые расчеты для случая хаотического режима колебаний. Оценить погрешность вычислений.