Для вычисления таких статистических средних, как размерность, энтропия, спектр показателей Ляпунова, и других характеристик аттрактора, необходимо иметь множество точек, определенных в фазовом пространстве размерности n и принадлежащих аттрактору. Число точек M в расчетах конечно, но обязано быть достаточно большим. Согласно формуле, предложенной в [34]
(1.13)
где D - размерность аттрактора. В случае, когда динамическая система задана дискретным оператором отображения, точки находятся автоматически после задания начальных условий. Если динамическая система задана системой дифференциальных уравнений, то в общем случае решение может быть найдено только численным интегрированием системы на компьютере. Обычно используют метод Рунге-Кутта 4-го порядка, погрешность задают 10-4-10-8, шаг счета определяется конкретной системой и должен быть выбран в сравнении с наименьшим из ее характерных времен.
Однако часто требуется вычислить характеристики аттрактора некоторой реальной системы, математическая модель которой неизвестна. При этом, как правило, неизвестна и размерность ее фазового пространства. В этой ситуации мы располагаем информацией о поведении во времени какой-либо одной из динамических переменных. К тому же и интервал времени экспериментальной реализации естественно ограничен. Можно ли в таких условиях получить характеристики аттрактора?
Путь к решению этой проблемы был предложен Такенсом. В [41] доказано, что почти для всех гладких динамических систем по имеющейся временной реализации одной наблюдаемой динамической переменной можно сконструировать новый аттрактор, основные свойства которого будут такими же, как у исходного.
Пусть имеется временной ряд экспериментальных данных, представляющий собой отсчеты некоторой физической величины: . Если известен шаг по времени , то время t = k*. Предполагается, что физическая величина s является одной из переменных динамической системы. Система находится в стационарном режиме, т.е. фазовая траектория проходит внутри аттрактора. Для восстановления аттрактора Такенсом предложен метод временной задержки координат. В n-мерном фазовом пространстве строится последовательность точек вида:
(1.14)
здесь - временная задержка, n - размерность вложения.
Основной результат Такенса состоит в следующем. Если , то множество точек задает вложение исходного аттрактора почти при любом выборе наблюдаемой переменной, если n не меньше удвоенной размерности исходного аттрактора. Для оценки характеристик реального исследуемого аттрактора можно вычислять характеристики восстановленного аттрактора. С целью уменьшения ошибки, обусловленной конечностью набора экспериментальных точек , необходимо проводить расчеты при нескольких различных значениях M и n и добиваться независимости получаемых оценок характеристик от M и n в пределах заданной точности.
Для малых шагов по времени значения sk и sk+1 будут близкими, поэтому большое значение приобретает правильный выбор временной задержки . Необходимо стремиться выбрать так, чтобы корреляция между sk и была по возможности минимальной. Традиционный способ выбора временной задержки состоит в вычислении автокорреляционной функции временного ряда:
(1.15)
Задержка выбирается равной времени первого пересечения нуля автокорреляционной функции. Второй способ [21] требует вычисления спектра мощности временного ряда, т.е. быстрого преобразования Фурье автокорреляционной функции. Если в спектре мощности присутствуют кратные пики, то задержка выбирается равной четверти периода самой высокой из доминирующих частот. Третий способ [23] основан на вычислении средней взаимной информации между двумя измерениями. Пусть даны два множества измерений A и B. Взаимная информация между элементом ai множества A и элементом bj множества B определяется как количество информации, которое имеют измерения ai и bj по отношению к друг другу:
(1.16)
Если измерения независимы, то взаимная информация равна нулю. Усредняя по всем измерениям, получаем:
(1.17)
Заменяя ai и bj на sk и соответственно, получаем среднюю взаимную информацию как функцию временной задержки . Задержка выбирается равной времени первого минимума во взаимной информации.
В случае модельных данных, когда нам известна размерность n фазового пространства динамической системы и все n координат каждой точки на аттракторе, корреляционную размерность D2 аттрактора находят следующим образом [26, 28].
Рассмотрим корреляционный интеграл C(r), показывающий относительное число пар точек аттрактора, находящихся на расстоянии, не большем r:
(1.18)
здесь - функция Хевисайда:
- расстояние в n-мерном фазовом пространстве,
m - число точек xi на аттракторе.
Если выполняется условие
(1.19)
то D2 считают корреляционной размерностью аттрактора.
Справедливость приведенного степенного закона ограничена значениями r, достаточно малыми по сравнению с размером аттрактора. При увеличении r величина C(r) достигает насыщения (при r, сравнимых с размером аттрактора). С другой стороны, при очень малых значениях r число пар точек xi,xj , расстояние между которыми не превышает r, становится малым (из-за конечности числа точек на аттракторе) и статистика становится бедной. Кроме того, приобретает решающее значение влияние инструментальных ошибок измерения сигнала. Следовательно, на практике степенной закон выполняется только в ограниченном диапазоне значений r (скейлинговом диапазоне), который и может быть использован для определения размерности аттрактора.
Учитывая, что из (1.19) следует
(1.20)
получаем оценку размерности аттрактора как тангенс угла наклона прямой, аппроксимирующей график корреляционного интеграла C(r) в двойном логарифмическом масштабе.
В случае экспериментальных данных мы обычно не знаем размерность фазового пространства системы и располагаем информацией только об одной координате точек на аттракторе. Поэтому все расчеты проводятся для нескольких размерностей фазового пространства n=1,2,3,... Для восстановления аттрактора используется метод Такенса. При этом корреляционная размерность аттрактора D2(n) сначала возрастает, но затем обычно выходит на постоянный уровень . Таким образом, получают искомую корреляционную размерность D2 аттрактора и оценку размерности фазового пространства системы . Если же выходной сигнал динамической системы сильно зашумлен, то размерность аттрактора постоянно растет.
Корреляционная энтропия K2 может быть вычислена достаточно просто [27]. Для этого также вычисляют корреляционный интеграл (1.18), но рассматривают не только его зависимость от расстояния r, но и от размерности фазового пространства n. При этом полагают, что
(1.21)
откуда
(1.22)
Энтропия K2 аппроксимируется в приемлемом диапазоне значений r и n.
Пусть динамическая система задана дискретным дифференцируемым отображением:
xk+1 = F(xk) , , DF(x) - матрица Якоби.
Бенеттин [14] придумал численный алгоритм вычисления показателей Ляпунова аттрактора такого отображения. Он начинается с ортонормированного базиса:
.
На k-м шаге вычисляют:
(1.23)
Затем применяют процедуру ортогонализации Грама-Шмидта:
, где (1.24)
На последующих шагах полагают:
(1.25)
где (u,v) - обычное скалярное произведение (индексы (k+1) опущены для простоты). Здесь - длина m-го вектора после ортогонализации, но перед нормализацией.
Показатели Ляпунова численно выражаются как
(1.26)
для больших k. Подобный алгоритм применим и к потоку , за исключением того, что вариационные уравнения
(1.27)
интегрируют от t до t+, начиная с ортонормированного базиса в t. Ортонормализационная процедура применяется как и выше, но множитель 1/k в (1.26) заменяется на 1/k.
Если математическая модель динамической системы неизвестна, то требуются другие алгоритмы нахождения показателей Ляпунова, поскольку в этом случае мы не знаем ни матрицы Якоби, ни размерности фазового пространства системы.
В 1985 году Вольф [46] описал алгоритм нахождения наибольшего показателя Ляпунова по временному ряду данных. Он следит за парой точек на аттракторе для оценки на каждом шаге по времени, по которому наибольший показатель Ляпунова может быть вычислен используя (1.26). Схематическая иллюстрация алгоритма приведена на рисунке 2. Он начинается с первой точки данных
(1.28)
снова меньше, чем и такой, что
Рисунок 2.
Процедура продолжается до тех пор, пока принятая за основу сравнения траектория y не дойдет до конца временного ряда. Наибольший показатель Ляпунова аттрактора оценивается как
(1.29)
где L - число шагов замены и m - общее число шагов по времени, в течение которых движется траектория y.
Каждая заменяющая точка должна лежать в том же направлении, что и старая, но из-за ограниченного числа отсчетов данных необходимы компромиссы. Первоначально поиск точки для замены ограничен в конусе угловой ширины и высоты около y(t). Обычно устанавливают равным . Значение увеличивают по мере необходимости пока сосед y не будет найден. Если это не приводит к успеху, то берут ближайшего соседа y безотносительно . Как показано в [46], ориентационные ошибки обычно приводят к малым ошибкам в оценке . Таким образом, результаты обычно нечувствительны к выбору .
Если не очень маленькое, то алгоритм Вольфа может быть полезен для определения того, хаотичный ли наблюдаемый временной ряд. Часто, однако, интересуются более чем одним первым показателем Ляпунова. Алгоритм может расширяться для измерения , но становится заметно сложнее. Он должен следить за треугольником точек и сохранять ориентацию относительно треугольника, которому выбирают замену. Этот подход становится плохим для более чем двух положительных показателей Ляпунова.В том же 1985 году Экманн и Рюэль [20] предложили метод оценки вариационных уравнений (1.27), который в принципе может дать все показатели Ляпунова аттрактора. Они следят, по крайней мере, за n+1 точкой, по которым методом наименьших квадратов оценивают матрицу Якоби DF(x). Показатели Ляпунова затем можно получить, проинтегрировав вариационные уравнения (1.27), используя алгоритм Бенеттина. Как показано в [39, 45], основной недостаток такого подхода - в зависимости результатов от размерности вложения аттрактора.