CC BY
6ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ И
ПРОЦЕССЫ УПРАВЛЕНИЯ N 3, 2006 Электронный журнал, рег. N П23275 от 07.03.97
http://www.neva.ru/journal e-mail: [email protected]
Моделирование динамических систем
РАЗРАБОТКА И РЕАЛИЗАЦИЯ АЛГОРИТМОВ ПОСТРОЕНИЯ СИМВОЛИЧЕСКОГО ОБРАЗА
Е.И.Петренко
Россия, 198904, Санкт-Петербург, Петродворец, Университетский пр., 28, Санкт-Петербургский Государственный Университет, математико-механический факультет, e-mail: [email protected]
Реферат.
Главной причиной интереса к изучению динамических систем послужило их все более возрастающее значение при исследовании процессов связанных с окружающим миром. Развитие теории динамических систем стимулировало активную разработку алгоритмов их исследования, а бурный рост возможностей вычислительной техники позволил широко применять методы компьютерного моделирования.
Для компьютерного моделирования и исследования динамической системы мы предлагаем использовать символический образ [5], который представляет собой ориентированный граф и является дискретизацией исходной динамической системы. Он строится по заданному покрытию фазового пространства ячейками {О,}. Вершины графа соответствуют ячейкам, дуги — связям между ними, а именно: вершины г и ] соединяются дугой, если образ ячейки О.1 при действии динамической системы пересекается с ячейкой О^. Многие задачи исследования динамических систем могут быть сведены к задачам исследования построенного ориентированного графа — символического образа.
В данной работе раскрываются аспекты реализации построения символического образа. Приводится способ введения координат, представления ячеек. Описаны четыре метода построения образа ячейки. Приводится их анализ и сравнение производительности.
1 Определения
Пусть М — подмножество т-мерного пространства Кт. Как правило, М является замкнутым ограниченным множеством (компактом) или гладким многообразием в Кт. Пусть Ъ — множество целых чисел и К — множество вещественных. Динамической системой называется непрерывное отображение Ф(х,£), где х Е М, £ Е Ъ (или £ Е К), такое, что
Ф : М х Ъ — М (или Ф : М х К — М),
Ф(х, 0) = х, (1)
Ф(Ф(х,£),в) = Ф(х, £ + 5), (2)
где £, в принадлежат Ъ (или К). Переменная £ трактуется как время, а многообразие М называется фазовым пространством. Если £ Е Ъ, то мы получаем динамическую систему с дискретным временем, которая называется каскадом или дискретной динамической системой. Часто дискретные динамические системы порождены итерационными процессами вида хп+1 = I(хп) или разностными уравнениями. В случае £ Е К, мы имеем дело с системой с непрерывным временем, которая называется непрерывной динамической системой. Иногда непрерывные динамические системы называются потоками. Как правило, непрерывные динамические системы порождены автономными системами дифференциальных уравнений х = I(х).
Мы будем рассматривать дискретные динамические системы, порожденные итерационными процессами вида хп+1 = I(хп), где I — диффеоморфизм, О — компакт, I : О С Кт — О задает динамическую систему с дискретным временем.
Определение 1 Степень отображения будем обозначать, и определять следующим образом.
/°(х) = х,
Гр(х) = 1р-1(1(х)),р > 2, (3)
I-р(х) = (I-1)р(х),р> 0.
Определение 2 Пусть хо £ О. Положительной траекторией или орбитой точки хо будем называть множество
{хо, X! = /(Х0),Х2 = /(хх) = /2(жо),...} . (4)
Отрицательная орбита или траектория точки х0 — это множество
{жо, хх = /-1(хо), Х2 = /-1(хх) = /-2(жо),...} . (5)
Определение 3 Если /(хо) = хо, то такая точка называется неподвижной точкой динамической системы.
Например, рассмотрим динамическую систему, заданную отображением
f : х —► 2006х. (6)
Неподвижной точкой этой динамической системы является точка хо = 0.
Определение 4 Если /к (хо) = хо, при к > 0 и для любого 0 < р < к /р(хо) = хо, то тогда точка хо называется к-периодической.
При численном моделировании динамических систем оказывается удобным рассматривать ^-траектории.
Определение 5 Пусть задан набор точек {хо,хх,...}. Тогда будем говорить, что задана е-траектория отображения /, если
\^+х - /(х01 <е, (7)
для всех г > 0
При помощи е-траектории мы можем задать погрешность определения орбиты или погрешность вычисления значения функции. Аналогично можно определить (р, е)-периодическую точку, е-неподвижную точку. Например, при е = 2, е-траекторией системы (6) будет являться траектория {1, 2005, 4022031,...}
2 Символический образ динамической системы
Пусть С = {Ох,..., Ок} — конечное покрытие компакта О замкнутыми множествами. Множества О.1 будем называть ячейками и / : О ^ О — диффеоморфизм. По множеству О, диффеоморфизму / и множеству ячеек С можно
построить ориентированный граф следующим образом: каждой ячейке О^ соответствует узел г графа. Между узлами г и ] есть ориентированное ребро тогда и только тогда, когда I(О^ П О^ = 0. Таким образом построенный ориентированный граф будем называть символическим образом динамической системы, заданной диффеоморфизмом I.
Практически любую задачу исследования динамической системы можно свести к задаче исследования ориентированного графа — символического образа [5]. Например, можно локализовать неподвижные точки, периодические траектории, оценивать спектр Морса динамической системы [4]. Между исходной системой и символическим образом существуют следующие связи:
• траектории системы образуют пути на символическом образе;
• символический образ отражает глобальную структуру динамической системы;
• символический образ является конечным приближением динамической системы;
• точность этого приближения зависит от максимального размера ячейки. Введем несколько параметров символического образа. Будем обозначать:
• ё — максимальный диаметр ячейки;
• г — нижняя граница символического образа. Определяется как минимальное расстояние между I(О¡) и множеством ячеек разбиения, которые не пересекаются с образом ячейки Оi.
Символический образ можно рассматривать как конечную дискретную аппроксимацию динамической системы. При этом более мелкое покрытие порождает более точную аппроксимацию. С помощью процесса последовательного подразбиения элементов покрытия можно строить последовательность символических образов и, тем самым, уточнять структурные характеристики системы.
Описанный метод может быть успешно применен к решению следующих задач [5]:
• локализация периодических траекторий заданного периода;
• построение периодической траектории;
• локализация цепно-рекуррентного множества;
• построение (положительно, отрицательно) инвариантного множества;
• построение аттрактора и его области притяжения;
• построение фильтраций и точной последовательности фильтраций;
• определение структурного графа динамической системы;
• оценка топологической энтропии;
• оценка показателей Ляпунова;
• оценка спектра Морса;
• построение изолирующих окрестностей инвариантных множеств.
Определение 6 Последовательность гк вершин графа С называется допустимым путем (или просто путем), если для любого к граф С содержит ориентированное ребро гк ^ zk+х. Путь называется периодическим, если последовательность {гк} является периодической.
Существует естественная связь между допустимыми путями на символическом образе С и е-траекториями отображения /. Можно сказать, что допустимый путь является следом е-траектории, причем обратное также верно. Справедлива следующая теорема [5]:
Теорема 1 [5] Справедливы следующие утверждения для символического образа С:
1. Если последовательность {гк} есть путь на символическом образе С и хк £ ОХк, тогда последовательность {хк} есть е-траектория гомеоморфизма / для любого е > д + ¿, где д и й — наибольшие диаметры ячеек О.1 и их образов /(О^ соответственно. В частности, если последовательность {гк} есть периодический путь на символическом образе, то последовательность {хк} есть е-периодическая траектория.
2. Если последовательность {гк} есть путь на символическом образе С, тогда существует последовательность {хк}, хк £ М(гк), которая есть е-траектория гомеоморфизма / для любого е > й.
3. Пусть последовательность {хк} есть е-траектория гомеоморфизма /, е < г и хк £ ОХк, где г нижняя граница символического образа. Тогда последовательность {гк} есть допустимый путь на символическом образе С. В частности, если последовательность {хк} есть е-периодическая траектория, то последовательность {гк} есть периодический путь на символическом образе С.
3 Анализ символического образа
При исследовании динамических систем наибольший интерес представляют неподвижные точки, периодические орбиты, а также цепно-рекуррентные множества.
Определение 7 Точка х называется цепно-рекуррентной, если для любого е > 0 существует периодическая е-траектория, проходящая через х. Цепно-рекуррентным множеством называется множество всех цепно-рекуррентных точек.
Цепно-рекуррентное множество инвариантно, замкнуто и содержит возвращающиеся траектории всех типов [5]. В частности цепно-рекуррентное множество содержит периодические орбиты.
Пусть Q — цепно-рекуррентное множество. Q(е) — множество точек е-траекторий. Тогда из определения цепно-рекуррентного множества и множества е-траекторий следует Q С Q(е) и Q(ех) С Q(е2), при ех < е2.
Определение 8 Вершина символического образа называется возвратной, если существует периодический путь, проходящий через нее. Две вершины символического образа называются эквивалентными, если существует периодический путь, проходящий через них.
Согласно определению, множество возвратных вершин разбивается на несколько классов эквивалентности. Ясно, что каждый периодический путь £ находится в некотором классе, который однозначно определяется по
Рассмотрим те ячейки О,, для которых вершина г является возвратной. Такое множество зависит от выбранного покрытия С и от наибольшего диаметра ячейки й. Поскольку зависимость от й для нас наиболее важно, обозначим
Р(й) = {х £ О,, : г — возвратные}.
(8)
Т(й) = {х £ Ок : к — невозвратные}.
Теорема 2 [5] Справедливы следующие утверждения:
1. Множество Р(() является замкнутой окрестностью цепно-рекуррентного множества. Кроме того, Р(() состоит из е-периодических точек для любого е > д + (, т. е.
Р(() С Я(е), е > д + (. (9)
2. Цепно-рекуррентное множество Q совпадает с пересечением множеств Р(() для всех положительных (:
Q = П Р (()• (10)
¿>о
3. Точки из множества Т(() не являются цепно-рекуррентными. Кроме того, если е < г, то не существует е-периодической траектории, проходящей через точку х множества Т((), т.е.
Q(е^ Т(() = 0, е < г. (11)
Множество Т(( ) является замкнутым по построению и пара Р(( ), Т(( ) образует замкнутое покрытие О. Следовательно, множество Р(ё)\Т(() является окрестностью цепно-рекуррентного множества Q. Из теоремы 2 следует включение:
Q С Q(еl) С О \ Т(() = Р(() \ Т(() С Р(() С Q(е2), (12)
где е1 < г < ( < д + ( < е2.
Отметим, что множества Р(() не являются монотонными по (, т. е. из условия (1 > (2 не обязательно следует, что множество Р(ё1) содержит Р((2). Однако, если С2 является подразбиением покрытия С1, то Р((2) С Р((1). Это свойство лежит в основе следующего алгоритма.
3.1 Алгоритм локализации цепно-рекуррентного множества
1. Строим исходное покрытие С компакта О. Находим символический образ О отображения I. Заметим, что ячейки исходного покрытия могут иметь произвольный диаметр
2. Выделяем на графе С возвратные вершины {%к}. Если множество таких вершин пустое, значит локализуемое цепно-рекуррентное множество является пустым и процесс его локализации прекращается. Иначе, используя их, находим замкнутую окрестность Р = {х £ О,к : г1к — возвратные} цепно-рекуррентного множества Q.
3. Разбиваем ячейки {О,к}, соответствующие возвратным вершинам символического образа и, таким образом, определяем новое покрытие.
4. Строим символический образ С для нового покрытия.
5. Переходим ко второму пункту, если размеры ячеек построенного символического образа достаточно велики.
Повторяя процесс последовательного измельчения покрытия, мы получаем последовательность окрестностей P0, P1, P2,... цепно-рекуррентного множества Q и последовательность наибольших диаметров do, di, , . . . ячеек, соответствующих возвратным вершинам символического образа для покрытия Ck. Следующая теорема обосновывает описанный алгоритм локализации множества Q.
Теорема 3 [5] Последовательность множеств P0,P1,P2,... обладает следующими свойствами:
1. Окрестности Pk вложены друг в друга, т. е.
Po D Pi D P2 D ... D Q, (13)
2. Если наибольшие диаметры dk ^ 0 при k ^ ж, то
lim Pk = П Pk = Q. (14)
k
k
Таким образом, описанный алгоритм дает монотонно убывающую последовательность окрестностей, сходящуюся к цепно-рекуррентному множеству.
Пусть G(V, E) — ориентированный граф, V — множество вершин, E — множество ребер.
Определение 9 Вершины v1 и v2 сильно связаны в G, если существует путь из v1 в v2 и из v2 в v1. Если все вершины в ориентированном графе сильно связаны, то G называется сильно связным.
Как было показано в [5], задача о локализации цепно-рекуррентного множества заданной динамической системы сводится к исследованию соответствующего символического образа и выделению на нем классов возвратных вершин. Из определения возвратной вершины следует, что компоненте сильной связности соответствует объединение классов возвратных вершин. Таким образом, выделение таких классов на графе эквивалентно нахождению компонент сильной связности. Для этого использован алгоритм Тарьяна [7], который обладает достаточно хорошей оценкой сложности: 0(п + т), где п — количество узлов, т — количество ребер.
3.2 Алгоритм выделения компонент сильной связности
Алгоритм основан на обходе графа в глубину и использует два стека "стек" и "маршрут". Стек "маршрут" содержит путь от начальной вершины до текущей. Каждая новая исследуемая вершина опускается в стек "маршрут", а при возвратах — извлекается. В "стек" добавляются все просмотренные вершины. Все элементы найденной компоненты сильной связности удаляются после ее окончательного формирования.
Заведем счетчик вершин с некоторым начальным значением и припишем к каждой вершине 2 числовых параметра: "номер" и "связка". Поле "номер" определяется простой последовательной нумерацией вершин по мере их обхода алгоритмом. Поле "связка" для произвольной вершины хранит номер другой вершины, которая была нумерована раньше. Если рассматриваемая вершина является корнем дерева поиска компоненты сильной связности, то значения полей "номер" и "связка" совпадут. Заметим, что значение поля "связка" всегда меньше или равно значению поля "номер". Схему работы алгоритма можно представить следующим образом.
• Шаг 1. Выбрать произвольную нерассмотренную вершину V, т.е. вершину для которой значение поля "номер" не было инициализировано.
• Шаг 2.
— Положить вершину V в стеки "стек" и "маршрут".
— Увеличить счетчик вершин на 1.
— Присвоить полям "номер" и "связка" этой вершины значение счетчика.
• Шаг 3. Выбрать некоторое нерассмотренное ребро, выходящее из вершины V и рассмотреть вершину (ш), в которую оно ведет.
- Если ребро идет в не рассмотренную ранее вершину, положить V = ш и перейти на Шаг 2.
- Если ребро идет в уже рассмотренную вершину, перейти на Шаг 4.
- Если вершина V не имеет неисследованных выходов и значение поля "связка" меньше значения поля "номер", перейти на Шаг 5.
- Если вершина V не имеет неисследованных выходов и значение поля "связка" равно значению поля "номер", перейти на Шаг 6.
• Шаг 4.
- Если значение поля "номер" вершины ш меньше значения поля "номер" вершины V и ш находится в стеке "стек", тогда положить значение поля "связка" вершины V равным минимуму из значений полей "связка" вершин V и ш.
- Перейти на Шаг 3.
• Шаг 5.
- Извлечь вершину V из стека "маршрут" и рассмотреть вершину, которая оказалась на вершине стека (и).
- Положить значение поля "связка" вершины и равным минимуму из значений полей "связка" вершин V и и. Положить V равным и и перейти на Шаг 3.
• Шаг 6.
- Взять все вершины с вершины стека "стек" до вершины V включительно и поместить их в новую компоненту сильной связности. Больше эта компонента изменяться не будет.
- Извлечь вершину V из стека "маршрут". Если в итоге "маршрут" окажется пустым, перейти на Шаг 1, иначе положить V равным вершине на вершине "маршрута" и перейти на Шаг 3.
Заметим, что мы можем получить компоненту сильной связности, состоящую из одной вершины, при этом у нее нет ребра, ведущего в нее саму. Для того чтобы не рассматривать такие "компоненты сильной связности", было добавлено дополнительное условие на шаге формирования компоненты
сильной связности (Шаг 6). А именно, если на верхушке стека "маршрут" лежит текущая вершина V, то она единственная в своей компоненте. Согласно описанию алгоритма она становится компонентой сильной связности. Если у этой вершины есть ребро идущее в неё, то это компонента сильной связности, иначе это просто проходящая вершина.
На Шаге 3 алгоритма требуется находить необработанные ребра. Для этого в каждой вершине графа будем хранить индекс последнего обработанного ребра.
На Шаге 4 требуется проверять, находится ли данная вершина в стеке "стек". Для ускорения этой проверки введем в каждой вершине флаг. Флаг будет установлен, если вершина находится в стеке. Для того, чтобы выяснить, содержится ли вершина в стеке будет достаточно посмотреть на значение этого флага.
На Шаге 6 для формирования компоненты сильной связности будем создавать новую структуру данных графа для сохранения выделенных вершин. Если для дальнейшего исследования требуются ребра, то после завершения работы алгоритма локализации компонент сильной связности, можно восстановить ребра в выделенной компоненте сильной связности.
Обозначим С исходный граф. Пусть для каждой компоненты связности графа С, задан набор узлов Ж,, входящих в нее. Тогда для построения графов компонент сильной связности СN графа С следует выполнить следующий алгоритм:
• Для каждой вершины V компоненты сильной связности N найти соответствующую ей вершину V в исходном графе С.
• Для всех ребер вида V ^ и', и' вершина С, если вершина и' £ N добавить это ребро в граф СN.
После завершения процесса локализации компонент сильной связности, мы получаем набор из графов-компонент сильной связности исходного графа. Каждый такой граф соответствует цепно-рекуррентному множеству исходной системы [5]. Теперь исходный граф можно удалить из памяти.
4 Представление ячейки
Ячейка представляется координатами ее вершин и центра. В реализации алгоритма удобно рассматривать одинаковые ячейки, тогда информацию о раз-
мерах ячейки можно хранить отдельно и всего один раз. Каждая ячейка в таком случае представляется точкой ее "верхнего левого угла" единственным образом.
Для избежания ошибок при работе с плавающей арифметикой на компьютере, рассматривается целочисленная система координат, за единицу длины, в которой принимается размер ячейки.
Введем обозначение [а, Ь] = [а1,Ь1] х • • • х [ат, Ьт].
Возьмем а,Ь Е Кт, такие, что а = Ь, тогда
О = [а, Ь] С Кт (15)
Зафиксируем т положительных чисел р1,р2,... ,рт и разобьем каждый отрезок [а-1,Ь-1] на р-1 равных частей:
( = Ьь—г,
\ Рг
^ = а + а - \)(1,а1+ а Ь] = Wi1 и Wi2 и •••и wpг, Wi1 П Wj2 = 0,л = п.
Определим множество всех ячеек разбиения.
W = ^^ х Wj2 х^^х Wjm |1 < ^ < pk, 1 < к < т}. (17)
Элемент множества W будем называть ячейкой. Целочисленной координатой ячейки с Е W будет набор значений (]1,]2,... ,]т), с которыми эта ячейка входит в множество W.
Пусть а, Ь Е К*, (р$= 1 — набор положительных чисел. Будем обозначать через W([а,Ь],р) множество ячеек, построенное описанным выше способом.
Использование представления ячеек в целочисленной системе координат позволяет уменьшить объем памяти, требуемой для хранения ячейки. Для вычисления координат вершин ячейки в исходной системе координат проводится линейное преобразование координат.
При построении последовательности символических образов ячейки, соответствующие вершинам символического образа разбиваются на равное количество частей. Удобно задавать количество частей по каждой координате отдельно. Пусть задан вектор разбиения ячеек г = (г1,г2,... ,гт), Г{ > 0. Тогда каждой ячейке с с координатами (х1, х2,..., хт) соответствуем множество ячеек (х1г1 + х^2 + ]2,..., хтгт + ]т), где ji Е [0, г^ П Ж.
5 Структура данных символического образа
Процесс построения символического образа представляет собой последовательную обработку ячеек. Для экономии памяти оказывается удобным хранить в графе только вершины, соответствующие обработанным ячейкам.
Для того чтобы обеспечить быстрое построение символического образа требуется, чтобы структура данных, представляющая ориентированный граф, была устроена так, чтобы оптимально выполнялись следующие операции:
• быстрый поиск вершин. При добавлении новой вершины, нужно проверить, не была ли эта вершина добавлена в граф раньше. И если была, то нужно просто вернуть указатель на нее и не производить никаких дополнительных действий;
• быстрый поиск ребер. При добавлении нового ребра, следует проверить, не было ли добавлено это ребро раньше.
• оптимальное использование памяти. Построенная нами структура данных должна работать с большими графами.
Граф представляется с помощью списка инциденции, т.е. списка вершин, каждая из которых содержит список исходящих из нее ребер. Для ускорения поиска списки ребер каждой вершины и список узлов хешируются при помощи хэш-таблиц.
6 Методы построения образа ячейки
Построение символического образа основано на построении образа ячейки. На рис. 1 показан образ ячейки (заштрихован темным цветом). Понятно, что в зависимости от выбора ячеек покрытия, попадающих в этот образ, мы будем получать различные графы.
6.1 Линейный метод
В этом методе мы оцениваем возможные значения функции на заданной ячейке. Для этого производятся следующие действия (рис. 2):
• берем ячейку разбиения О,;
/Ш
\
\
I
/ /
Рис. 1: Разбиение множества О, образ ячейки.
• вычисляем значения функции системы / на вершинах ячейки О,;
• строим минимальный прямоугольник Е, ориентированный по осям координат, который содержит образы вершин ячейки О,;
• образом ячейки будем считать множество ячеек покрытия, которые пересекаются с построенным прямоугольником Е.
Для ускорения работы этого метода значения функции системы в вершинах ячейки вычисляются с точностью до членов первого порядка при разложении этой функции в ряд Тейлора в окрестности центра ячейки х0, т.е.
/(х) = /(хо) + О](хо)(х - хо), (18)
где х — вершина ячейки.
Чтобы увеличить вероятность того, что построенный этим методом образ ячейки содержит все точки образа ячейки, можно ввести коэффициент расширения полученного прямоугольника.
6.2 Точечный метод
Пусть дано число п £ Мт, т.е. п = (щ,п2, • • • ,пт), где п, £ N. Тогда точечный метод с дроблением п действует следующим образом. В ячейке равномерно берем п, точек по 1-й координате. Образом ячейки будем называть
Рис. 2: Построения образа ячейки. Линейный метод.
набор ячеек, в которые попали образы равномерно брошенных точек (рис. 3). Результат и время работы метода зависит от выбора значения вектора п.
Рис. 3: Построение образа ячейки при N = (2, 2). Точечный метод.
6.3 Улучшенный точечный метод
Пусть даны число п £ Nm и Ь £ [0,1]т, тогда улучшенный точечный метод с дроблением п и наложением Ь действует аналогично точечному методу с дроблением п. Образу равномерно брошенной точки х будем сопоставлять ячейку, которая содержит образ /(х) и все ячейки, которые находятся ближе чем ^ от этой точки, где % номер оси, в направлении которой идет отрезок-граница ячейки.
Если образ некоторой равномерно брошенной точки попадет на светлосерую область ячейки (рис. 4), то тогда добавляем к результату не только эту ячейку, но и соседнюю, как показано на рис. 4. Если точка лежит близко
от нескольких границ ячейки, то следует добавлять несколько ячеек (рис. 4).
■
Области ячейки Точка на краю ячейки
□Я
Точка, близкая к двум границам
Рис. 4: Точки на краю ячейки
Для поиска соседних ячеек был разработан специальный алгоритм, который довольно эффективно вычисляет целочисленные координаты всех ячеек, которые должны быть добавлены.
6.4 Адаптивный метод
Недостаток описанных выше методов построения образа ячейки заключается в том, что эти методы вычисляют функцию в заранее определенных точках. Предлагается новый метод построения образа ячейки, который будет вычислять значение функции только в тех точках, выбранных в процессе работы алгоритма, в которых это действительно необходимо, при этом количество
ш
таких точек, по которым строится образ ячейки тоже будет зависеть от поведения системы на данной ячейке. Такой способ построения образа ячейки будем называть адаптивным методом.
Адаптивный метод работает по принципу точечного метода (см. 6.2). Единственное отличие его заключено в том, что точки, в которых вычисляются значения функции берутся не равномерно, как раньше, а исходя из поведения системы.
Определение 10 Пусть есть множество точек Р С О и рефлексивное симметричное отношение ф на множестве Р. Будем называть такое отношение ф отношением соседства. Пусть задана функция правой части системы / = (/1} ¡2,..., /т). Фиксируем £ = (е!ь)'т=1, £1 > 0. Тогда будем называть множество Р £-множеством для функции / и отношения ф, если для любых двух соседних (ф(х,у)) точек Х,у Е Р выполнено соотношение
Ц\(Х) - ¡г(У)1 < £г. (19)
Будем обозначать такое множество через Р^.
Определение 11 Пусть задан набор т положительных чисел £ = (£,)'т=1. Тогда £-окрестностью точки х = (х1,х2,... ,хт) Е Кт будем называть множество
[Х1 - £1 ,Х1 + £1] х [Х2 - £2,Х2 + £2] х • • • х [Хт - £т, Хт + £т] С . (20)
Пусть задано множество Р^ ^, и ячейка Оз, все вершины которой принадлежат множеству Р. Тогда образ ячейки будет состоять из £-окрестностей образов точек, которые лежат в Р П Оз. При этом если элементы £ достаточно малы, т.е.
£ < | > (21)
то можно искать ячейки по принципу, описанному в улучшенном точечном методе (см. 6.3).
Введенное отношение соседства ф можно представить как неориентированный граф, где вершины соответствуют точкам множества Р, а ф(Х,у) эквивалентно существованию ребра между вершинами Х и у или совпадению этих вершин Х = у.
6.5 Алгоритм построения графа для Р^ф
Предположим, задано начальное множество точек и граф, определяющий соотношение соседства для точек множества. Зафиксируем £ > 0.
Определение 12 Вектором длины ребра т между узлами графа х и у будем называть вектор пространства компонентами которого будут расстояния между проекциями образов точек х и у под действием функции системы /.
¡впдЩт), = |/,(х) - /,(у)|. (22)
На множестве векторов длин можно ввести отношение порядка. Будем говорить, что вектор а больше вектора Ь, если это соотношение выполнено покомпонентно.
Определение 13 Обозначим через Е(п) множество всех вершин, которые соединены одним ребром с вершиной п.
Для того, чтобы заданное множество удовлетворяло требованиям определения, нужно проверить, что все ребра не длиннее £. Представим все ребра исходного графа в виде упорядоченного по убыванию длины ребра списка Ь. В качестве длины будем брать сумму модулей компонент вектора длин.
Пока список Ь не пуст, возьмем ребро т из начала списка и удалим его из списка. Если ребро т между точками х и у такое, что вектор его длины больше £, выполним следующее:
1. Добавим в граф разбиения новую вершину р с координатами (у )т=1.
2. Удалим ребро т из графа.
3. Добавим ребра (х,р) и (у,р) в граф и в список Ь, с сохранением порядка.
4. Добавим ребра к вершинам, принадлежащим Е(х)ПЕ(у) в граф и список Ь с сохранением порядка.
Основная сложность в данном алгоритме возникает на шаге 4. Нужно определить, какие ребра добавлять к новому узлу графа, т.е. следует определить расширение отношения ф на новую точку р. Следует наиболее полно расширить отношение ф, чтобы учесть все расстояния от новой точки р до остальных точек графа.
На практике возникают ситуации, когда описанный процесс работает очень медленно из-за большого количества добавлений новых вершин в граф. В этом случае мы вводим ограничение на количество вершин в графе разбиения. А именно: множество ребер графа делится на два класса: ребра с известной длиной и ребра, длина которых не определена. Алгоритм построения образа ячейки по графу разбиения работает аналогично приведенному выше, однако, при построении £-окрестности для вершины этого графа берем £ равным длине наибольшего ребра. Можно приписать каждому ребру его длину и хранить упорядоченный список ребер для каждой вершины. В вершинах можно хранить значение образа точки. Введенные оптимизации позволили увеличить скорость работы алгоритма.
Пример работы алгоритма для 2-х мерного случая
Рассмотрим граф для множества Р, которое состоит из 4-х точек (рис. 5). Предположим мы проверяем длину выделенного жирно ребра. Тогда, согласно алгоритму, следует ввести новую точку и некоторое число новых ребер.
Пометим двойной линией ребра, которые имеют одну общую вершину с рассматриваемым (выделенным жирно) ребром. Тогда, согласно алгоритму, следует добавить новые ребра только к тем вершинам, которые соединены ребрами с двумя вершинами рассматриваемого ребра (рис. 5).
Исходный граф Дробление ребра
Рис. 5: Разбиение длинного ребра.
Выделим такие вершины кругом. Тогда появляются два новых ребра, которые нарисованы серым цветом. Эти ребра будут добавлены в конец списка Ь. На следующих шагах алгоритм проверит длины оставшихся ребер.
6.6 Начальное разбиение ячейки
Вершины графа разбиения соответствуют точкам фазового пространства исследуемой системы. Результат и скорость работы адаптивного метода зависит от выбора начального разбиения ячейки. Мы предлагаем следующие графы разбиения ячейки для одномерного, двумерного и трехмерного случаев (рис. 6):
• для одномерного случая достаточно взять взять граф с одним ребром и двумя вершинами, которые соответствуют границам одномерной ячейки;
• в двумерном случае будем рассматривать граф из 5 вершин — центр ячейки и ее четыре вершины и 8 ребер, которые связывают соседние вершины;
• в трехмерном случае мы рассматриваем все вершины ячейки, все центры ее граней и центр ячейки. Соединяем ребром соседние вершины.
1-мерная ячейка
2-мерная ячейка
3-мерная ячейка
Рис. 6: Графы начального разбиения
7 Примеры
В этой части приведены примеры построения цепно-рекуррентных множеств динамических систем с помощью построения последовательности символических образов с использованием разных методов построения образа ячейки. Приведены рисунки полученных цепно-рекуррентных множеств и таблицы сравнений работы различных методов построения образа ячейки. В таблице также показано время работы алгоритма и количество компонент сильной связности графа.
Количество шагов процесса подразбиения символического образа определяется экспериментально для каждой системы. Следующий шаг построения символического образа проводится, если для этого хватит оперативной памяти компьютера.
На каждом шаге построения символического образа будем разбивать исходную ячейку на 2 части по каждой координате. Точечный метод будет строить образ ячейки по 2 точкам по каждой координате, улучшенный точечный метод использует наложение, равное 10% от размера ячейки. Ограничение на количество точек для построения образа ячейки адаптивным методом не устанавливается.
Все эксперименты были проведены на компьютере Intel Pentium 4 3Ghz, 1Gb оперативной памяти. Операционная система Microsoft Windows XP SP2.
7.1 Отображение Хенона [10]
Рассмотрим систему:
'Л ^ - аХ + Ьу), (23)
при а = 1.4, Ь = 0.3.
Известно [5], что при заданных параметрах у этой системы существует инвариантное множество в области [-1.5,1.5] х [-у, ]. Возьмем О = [-10,10] х [-10,10]. Проведем 13 шагов построения символического образа. Сведем результаты вычислений в таблицу. На каждом шаге будем выполнять построение символического образа и его анализ, согласно алгоритму, приведенному в 3.1, т.е. выделение цепно-рекуррентного множества с помощью
алгоритма выделения компонент сильной связности.
Метод Кол-во Кол-во Время Кол-во
узлов ребер работы компонент
Линейный метод 572 435 5 509 083 24 047 ШЙ 2
Точечный метод 298 397 853 915 40 656 ШЙ 5
Улучшенный точечный метод 505 152 3 877 289 465 563 ШЙ 3
Адаптивный метод 507 637 5 171 045 464 844 шй 2
Каждым методом было сделано по 13 шагов. Начальное разбиение — по 10 ячеек на каждую координату. На каждом шаге ячейка разбивается на 2 части по каждой координате.
Приведем иллюстрации построенного цепно-рекуррентного множества. Построенные цепно-рекуррентрые множества получились практически одинаковые, отличие видно только по количеству узлов. Самый маленький граф получается при использовании точечного метода.
Линейный метод
Адаптивный метод
Точечный метод
Улучшенный точечный метод
Рис. 7: Результаты вычислений. Отображение Хенона.
7.2 Отображение Икеда [11]
Отображение имеет вид:
C
1 + x2 + y
x\ id — C2(x cos t (x,y) — y sin т (x,y)) y) V C2(x sin T(x,y)+ y cos T(x,y))
т(x, y) = C — 1 , J 2 (24)
(25)
где А = 2, С = 0.4, С2 = 0.9, Сз = 6.
Пусть Б = [—10,10] х [—10,10]. Отображение Икеда возникает при моделировании оптических носителей (кристаллов) информации. Оно было рассмотрено и исследовано в [11, 13]; было показано, что при заданных параметрах это отображение имеет цепно-рекуррентное множество.
Метод Кол-во Кол-во Время Кол-во Кол-во
узлов ребер работы компонент шагов
Линейный метод 3 117 236 25 988 943 172 484 ms 2 11
Точечный метод 1 732 774 4 378 370 112 359 ms 4 11
Улучшенный точечный метод 2 949 314 21 666 315 512 969 ms 2 11
Адаптивный метод 1 017 624 11 740 358 640 609 ms 2 10
Результаты вычисления символического образа отображения Икеда приведены на рис. 8. Видно, что результат, полученный точечным методом немного более "дырявый", по сравнению с результатами, полученными другими методами. Согласно [13], наиболее хорошие результаты дали линейный и адаптивный методы.
Точечный метод Улучшенный точечный метод
Рис. 8: Результаты вычислений. Отображение Икеда.
7.3 Множество Жюлиа для г = -0.12 + 0.74г [6, 12]
В исследованиях Пайтгена и Рихтера [6] были приведены алгоритмы и результаты построения множеств Жюлиа и наполненных множеств Жюлиа для комплексного отображения
2 ^ 22 + Г. (26)
Система, эквивалентная (26), [12]:
Л ^ ( '-<2 - у2 + те _ (27)
у I \ 2ху + гт
где те = -0.12, гт = 0.74
Описанные выше методы построения инвариантных множеств с помощью символического образа дали похожие результаты.
Метод Кол-во Кол-во Время Кол-во Кол-во
узлов ребер работы компонент шагов
Линейный метод 2 412 630 20 337 231 1 174 000 ШЙ 4 13
Точечный метод 92 140 1 406 ШЙ 15 15
Улучшенный точечный метод 2 164 382 14 282 466 844 844 ШЙ 6 13
Адаптивный метод 1 068 206 12 145 512 2 058 797 шй 2 12
На этом примере точечный метод показал очень странный результат. Это можно объяснить тем, что количество точек, взятых для построения образа ячейки слишком мало. Изображения построенных цепно-рекуррентых множеств приведены на рис. 9.
Линейный метод
Адаптивный метод
Точечный метод
Улучшенный точечный метод
Рис. 9: Результаты вычислений. Отображение 2 ^ х2 + г, г = —0.14 + 0.741
7.4 Множество ^Кюлиа в случае параболической неподвижной точки [6]
Рассмотрим отображение на комплексной плоскости
f : х ^ х2 + exp2пг20 х.
;28)
Для ^-периодической точки 20 ее орбита представляет собой множество точек {20, Г(20), Г2(20),..., ¡к-1(го)} ( цикл периода к).
Определение 14 Пусть 20 периодическая точка периода п. Собственным значением точки х0 для отображения Г называется комплексное число
А = (ГУ (20)
По правилам дифференцирования сложной функции, мы получаем, что это число одинаково для любой точки цикла.
Периодическая или неподвижная точка 20 называется:
• притягивающей, если 0 < |А| < 1,
• отталкивающей, если |А| > 1,
• сверхпритягивающей, если А = 0,
• нейтральной, если |А| = 1.
Нейтральные точки в свою очередь тоже можно классифицировать. Пусть | А| = 1 и А — собственное значение точки 20 для отображения Г, тогда его можно переписать в виде
А = е2пга, а е [0,1]. (29)
Неподвижная точка 20 называется параболической, если а е Q.
При заданных параметрах, для исследуемой системы 20 = 0 является неподвижной параболической точкой. Обозначим А — собственное число функции Г в точке 20. Поскольку А20 = 1, множество Жюлиа подходит к точке 20 = 0 с 20 различных направлений (между 20 лепестками). Основная проблема моделирования поведения этой системы заключается в потере точности компьютерных вычислений, в результате чего возникают лишние периодические или неподвижные точки.
В этом примере техника символического образа для построения цепно-реккурентных множеств позволяет построить лишь наполненные множества Жюлиа, т.е. вышеописанная проблема остается.
Отметим, что в работах [2, 6] были приведены методы решения этой задачи.
Проведем 10 шагов построения символического образа (см. 3.1). Полученные результаты сведем в таблицу:
Метод Кол-во Кол-во Время Кол-во Кол-во
узлов ребер работы компонент шагов
Линейный метод 2 737 711 24 507 847 182 922 шэ 2 10
Точечный метод 2 684 661 14 941 697 1 145 969 шэ 1 10
Улучшенный точечный метод 2 801 326 32 547 551 4 795 844 шэ 1 10
Адаптивный метод 2 814 564 35 579 838 9 043 422 шэ 1 10
Все методы дали визуально схожие результаты (рис. 10). Полученные нами результаты согласуются с результатами, полученными в [2, 6].
Линейный метод
Адаптивный метод
Точечный метод
Улучшенный точечный метод
Рис. 10: Результаты вычислений в случае параболической неподвижной точки.
7.5 Двойное логистическое отображение [9]
Рассмотрим динамическую систему
х У
—>
(1 — а)х + аЬу(1 — у) (1 — а)у + аЬх(1 — х)
30)
где а = 0.38, Ь = 4.11.
Данная система, была исследована в работах ([1], [9]). Были получены оценки на область притяжения к бесконечно удаленной точке, в частности было показано, что в этой области лежит внешность круга
(х - г)2 + (у - г)2 < 2г2, (31)
где г = 1~ас+аЬ.
При заданных значениях параметров инвариантным множеством является объединением инвариантных кривых бифуркации Хопфа в окрестности двух неподвижных точек, расположенных симметрично относительно главной диагонали.
Проведем построение символического образа и сравним полученные численные результаты.
Метод Кол-во Кол-во Время Кол-во Кол-во
узлов ребер работы компонент шагов
Линейный метод 393 478 2 715 484 16 875 шв 15 13
Точечный метод 258 874 1 108 540 57 188 шв 13 13
Улучшенный точечный метод 492 115 4 868 381 320 266 шв 13 13
Адаптивный метод 503 102 5 275 606 517 078 шв 13 13
Цепно-рекуррентые множества показаны на рис. 11. На этой системе различия методов построения образа ячейки не видны. Однако по численным показателям отличия все же есть.
Линейный метод
Адаптивный метод
Точечный метод
Улучшенный точечный метод
Рис. 11: Двойное логистическое отображение. а = 0.38, Ь = 4.11.
7.6 Отображение с задержкой [8]
Рассмотрим динамическую систему с задержкой.
х у
—»
у
ау(1 — х)
32)
Система обладает двумя неподвижными точками 0\ = (0,0) и 02 = (1 — а, 1 — а)• Начало координат является седловой точкой при а > 1. Собственные числа в точке О1 равны 0 и а. Точка 02 является фокусом при а > |. Собственные числа в этой точке равны Х12 = 1 (1 ± у/5 — 4а). Фокус устойчив при а < 2, а при а > 2 фокус теряет устойчивость через бифуркацию Хопфа. Возникающая при этой бифуркации инвариантная кривая разрушается при а = 2.27 с появлением странного аттрактора [5]. Кроме того, при этом значении параметра неустойчивая сепаратриса седловой точки касается устойчивой, лежащей на ости ОХ.
Проведем построение цепно-рекуррентного множества при а = 2.21 (рис. 12). Разобьем каждый шаг построения символического образа с меньшим разбиением на два этапа. На первом этапе построения будем делить ячейки только по первой координате на две равные части. На втором — только вторые координаты. Каждый этап представляет собой отдельный шаг построения уточненного символического образа. Такой способ позволяет получить символической образ с более мелким разбиением и затратить на это меньше оперативной памяти, чем аналогичный процесс построения образа ячейки, на каждом шагу которого деление проводилось сразу по нескольким координатам.
Метод Кол-во Кол-во Время Кол-во Кол-во
узлов ребер работы компонент шагов
Линейный метод 1 489 037 12 863 578 357 859 ШЙ 5 17
Точечный метод 378 579 954 108 311 875 ШЙ 7 17
Улучшенный точечный метод 1 244 114 8 723 417 13 201 563 ШЙ 6 17
Адаптивный метод 1 390 661 13 941 695 18 473 141 шй 6 17
Линейный метод
Адаптивный метод
Точечный метод
Улучшенный точечный метод
Рис. 12: Отображение с задержкой, а = 2.21.
Проведем построение цепно-рекуррентного множества при а = 2.27.
Метод Кол-во узлов Кол-во ребер Время работы Кол-во компонент Кол-во шагов
Линейный метод 422 636 3 504 285 17 297 шэ 2 13
Точечный метод 185 861 450 270 28 406 шэ 3 13
Улучшенный точечный метод 367 115 2 501 061 252 625 шэ 2 13
Адаптивный метод 381 186 3 712 343 388 609 шэ 3 13
Изображение цепно-рекуррентых множеств дано на рис. 13. На рис. 14 показана область вблизи седловой точки, где возникает сложное поведение ее инвариантных многообразий. Визуально все методы построения образа ячейки работают примерно одинакого. Согласно таблице, самым быстрым методом является линейный метод, самым медленным — адаптивный.
Точечный метод Улучшенный точечный метод
Рис. 13: Отображение с задержкой при а = 2.27
Линейный метод
Адаптивный метод
Точечный метод
Улучшенный точечный метод
Рис. 14: Отображение с задержкой, а = 2.27. Область [0, 0.1] х [0, 0.1]
7.7 3х мерная модель Хищник-^Кертва-Суперхищник [14]
В работах [5, 14] рассматривается следующая динамическая система
(Л
у
—>
4.522х ехр(-у)
1+х тах(ехр(-у),к(г)к(у))
0.962ху ехр(—г)к(4уг )к(у) 4ух
\
/
(33)
_ 1— ехр(- Ь)
где к{Ь) =
Пусть В = [0.01,10] х [0.01,10] х [0.01,10]. В этом примере рассматривается трехмерная динамическая система. Для хранения ячейки требуется больше памяти, время вычисления образа одной точки тоже увеличивается. Наиболее хорошие результаты построения инвариантного множества для этой системы можно получить лишь используя все методы построения образа ячейки в совокупности. Эффективным по памяти оказывается на каждом шаге дробить ячейки только на две части (только по одной координате). В этом случае прирост ячеек в построенном графе оказывается минимальным, что позволяет уложиться в оперативную память компьютера и провести больше итераций построения символического образа, тем самым построить символический образ с меньшим разбиением.
Значение константы Липшица системы на множестве В оказывается очень большим, в результате чего адаптивный метод может работать очень медленно. Приходится вводить ограничения на количество узлов в графе представления Р^.
Результаты построения цепно-рекуррентного множества приведены в таблице и на рис. 15.
Результаты работы каждого из методов для этого примера сильно отличаются. Видно, что быстрее всего работает точечный метод, и самым долгим является адаптивный метод, несмотря на введенные ограничения на граф (не более 100 вершин в графе адаптивного метода для ячейки).
Метод Кол-во Кол-во Время Кол-во Кол-во
узлов ребер работы компонент шагов
Линейный метод 287 627 11 399 582 30 219 шэ 1 7
Точечный метод 6 744 33 956 5 297 шэ 7 7
Улучшенный точечный метод 110 361 3 771 889 43 359 шэ 5 7
Адаптивный метод 150 514 6 525 427 1 167 094 шэ 1 7
Линейный метод
Адаптивный метод
Точечный метод
Улучшенный точечный метод
Рис. 15: Результаты вычислений. 3-х мерная модель хищник-жертва-суперхищник.
8 Заключение
Сравним средние скорости работы описанных методов для исследованных двумерных систем. Для этого для каждой из двумерных систем для каждого метода вычислим среднее количество ячеек V, ребер е, построенных за 1 миллисекунду (шб), и значение отношения числа ребер к числу ячеек V •
Усредним эти значения по всем исследованным системам для каждого ме-
тода (V, ё и в/у). Где под V понимается среднее арифметическое. По каждому из устредняемых параметров вычислим среднее квадратичное отклонение по
формуле у г=1 N—, где N количество значений. Сведем полученные результаты в таблицу.
Из собранных данных видно, что самым производительным оказывается линейный метод. Самым быстрым оказывается точечный метод, но среднее количество ребер на один узел покрытия оказывается довольно маленьким. По этому показателю лидирует адаптивный метод, для которого в среднем на каждую ячейку приходится чуть более 10 ребер. Однако отношение числа ячеек к ребрам меньше всего зависит от системы, если мы применяем линейный метод.
На практике самым удобным оказывается линейный метод. Он позволяет получить достаточно хорошие результаты при минимальных затратах времени. Для приближенного рассмотрения системы удобным оказывается использовать точечный метод или улучшенный точечный метод. Они работают быстрее, но результат получается менее точны. Для некоторых систем оказывается удобным использовать адаптивный метод. Этот метод является наиболее сложным и наиболее медленным, он позволяет получить наибольшее значение отношения среднего числа ребер исходящих из каждого узла. Исходя из среднего квадратичного отклонения среднего числа исходящих ребер из каждого узла можно считать, что результат работы этого метода сильнее зависит от поведения системы, чем результат работы линейного метода.
Метод Ср. кол-во узлов за шй V Ср. квадр. отклонение Время работы
ЕЫ-*)2
V N
Линейный метод 15.802 8.6635 1 945 484 ШЙ
Точечный метод 5.4305 4.7180 1 696 593 ШЙ
Улучшенный точечный метод 2.0527 2.1462 29 589 571 ШЙ
Адаптивный метод 0.7916 0.4803 31 858 960 ШЙ
Метод Кол-во ребер за шй в Ср. квадр. отклонение Время работы
V N
Линейный метод 132.9332 73.6511 1 945 484 ШЙ
Точечный метод 16.0435 11.7427 1 696 593 ШЙ
Улучшенный точечный метод 9.8190 4.9137 29 589 571 ШЙ
Адаптивный метод 8.5444 5.2827 31 858 960 шй
Метод Ребер на Ср. квадр. Время
узел отклонение работы
в/у
V N
Линейный метод 8.4534 0.7656 1 945 484 ШЙ
Точечный метод 3.1109 1.2832 1 696 593 ШЙ
Улучшенный точечный метод 7.3512 2.9207 29 589 571 ШЙ
Адаптивный метод 10.8550 0.9591 31 858 960 шй
где N = 7 количество рассмотренных двумерных систем.
Список литературы
[1] Ампилова Н.Б., Осипов А.В. Локальные бифуркации для полного отображения Гардини. Деп. ВИНИТИ 14.06.96,N 1969-B96.
[2] Ампилова Н.Б., Петренко Е.И. Алгоритм компьютерного моделирования множества Жюлиа в случае неподвижной параболической точки. // Электронный Журнал Дифференциальные Уравнения и Процессы Управления (http://www.neva.ru/journal), 2,2002
[3] Осипенко Г.С. О символическом образе динамической системы // сб. Граничные задачи, Пермь, 1983, с.101-105.
[4] Осипенко Г.С., Романовский И.В., Ампилова Н.Б., Петренко Е.И., О вычисление спектра Морса // Проблемы Математического Анализа, Выпуск 27, Январь 2004, с. 151-169.
[5] Осипенко Г.С., Ампилова Н.Б. Введение в символический анализ динамических систем. // Изд. СПбГУ, 2005. ISBN: 5-288-03656-Х
[6] Пайтген Х.й, Рихтер П.Х. Красота фракталов. // Образы комплексных динамических систем. M, 1993.
[7] Писсанецки С. Технология разреженных матриц. // М., Мир, 1988.
[8] Aronson D.G.,Chory M.A., Hall G.R. et.all Bifurcation from an invariant circle for two-parameter families of maps of the plane: A computer-assisted study. // Commun.Math.Phys.83,3(1982), p.303-354.
[9] Gardini L., Abraham R., Record R.J., Fournier-Prunaret D. A double logistic map. // Int.J.Bif.and Chaos,4,1(1994), p.145-176.
[10] Henon M. A two-dimensional mapping with a strange attractor. // Comm. Math.Phys. v.50,69-77(1976).
[11] Ikeda K. Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system // Opt. Comm. 1979, Vol.30, p.257-261.
[12] Julia G. Sur l'iteration des fonctions rationneles. // Journal de Math.Pure at Appl. 8:47-245. 1918
[13] Osipenko G. Numerical Explorations of the Ikeda mapping dynamics // Electronic Journal of Differentional Equations and Control Processes (http://www.neva.ru/journal), Vol.2, 2004
[14] Linsdrom T. Dependencies between competition and predation — and their consequences for initial value sensitivity, // SIAM J. Appl. Math. Vol59, pp.1468-1486
