Содержание к диссертации
Введение
1. Моделирование регулярных и хаотических режимов динамических систем 18
1.1. Современные принципы исследования динамических систем 18
1.2. Методы построения отображения Пуанкаре 26
1.3. Поиск периодических решений и анализ их устойчивости 28
1.4. Продолжение семейств периодических решений 31
1.5. Симметричные периодические решения 42
1.6. Исследование сценариев перехода к динамическому хаосу 48
1.7. Особенности небесно-механических задач 50
2. Описание комплекса программ 56
2.1. Общая структура комплекса программ 56
2.2. Программы численного интегрирования 64
2.3. Программы поиска и продолжения периодических решений 69
2.4. Программы исследования перехода к динамическому хаосу 72
3. Результаты исследования плоской задачи Хилла 75
3.1. Общие свойства задачи Хилла 75
3.2. Основные семейства периодических решений задачи Хилла 81
3.3. Классификация периодических решений второго рода . 90
3.4. Каскады бифуркаций удвоения периода в задаче Хилла 108
3.5. Расщепление инвариантных многообразий 111
3.6. Глобальная динамика 112
Заключение 114
Литература 116
Приложения
- Поиск периодических решений и анализ их устойчивости
- Исследование сценариев перехода к динамическому хаосу
- Программы поиска и продолжения периодических решений
- Классификация периодических решений второго рода
Введение к работе
Формулируя основные фундаментальные направления развития математической науки в XXI столетии, С. Смейл [1] отнес к ним систематическое развитие компьютерной динамики. Вместе с тем далеко не все задачи численного исследования динамических систем решены. Этому отчасти мешает сложившаяся традиция отдельно рассматривать системы с непрерывным временем и порождаемые ими потоки с дискретным временем. В каждом из этих направлений разработаны методы обнаружения, продолжения и анализа перестроек инвариантных структур. Известно значительное число пакетов прикладных программ, реализующих различные методы исследования динамических систем общего вида. Существует сайт , содержащий обзор программного обеспечения для исследования динамических систем. В этом обзоре имеется информация о двух десятках программных пакетов, относящихся к данной теме.
На начальном этапе исследования гамильтоновых систем, возникающих при изучении небесно-механических задач, автору пришлось применять различные пакеты прикладных программ. Однако, ни один из таких программных комплексов, ни даже их комбинация не могли обеспечить необходимого качества численных исследований. Дадим описание комплексов программ для исследования динамических систем. Выбор комплексов определялся на основании следующих факторов:
доступность пакета в виде исполняемого модуля или исходных текстов;
наличие в пакете средств исследования инвариантных структур динамической системы;
адаптивность пакета.
Этим требованиям удовлетворяют пакеты WinSet, DESIR, CONTENT, AUTO2000.
Программа WinSet предназначена для исследования и демонстрации инвариантных множеств целого ряда динамических систем. Она является приложением к книге [2]. Заявленная в описании к программе
возможность построения отображений Пуанкаре реализована либо для аналитически заданных отображений, либо для систем обыкновенных дифференциальных уравнений (ОДУ) с полутора степенями свободы. В качестве встроенного интегратора систему ОДУ используется явный метод Рунге-Кутты 4-го порядка с автоматическим выбором шага интегрирования в соответствии с априорно заданной локальной погрешностью на шаге. Использование этого метода для сильно неустойчивых фазовых траекторий приводит к быстрому накоплению ошибки и, как следствие, к неадекватному представлению траекторий и сечений Пуанкаре. Таким образом, применение этого программного средства возможно исключительно с целью первоначального ознакомления со структурой фазового пространства исследуемой системы, но не для ее количественного изучения.
Пакет DESIR разработан на механико-математическом факультете Ростовского государственного университета Говорухиным В.Н. Этот программный комплекс представляет собой набор исполняемых в операционной системе DOS модулей, позволяющих в интерактивном режиме исследовать систему ОДУ в нормальной форме. Возможно также исследование отображений. Следует отметить, что после определения исходных данных пользователь получает готовый исполняемый модуль с откомпилированными процедурами интегрирования системы ОДУ и ее уравнения в вариациях. Пакет предлагает широкий выбор современных интеграторов (в частности, хорошо зарекомендовавшие себя алгоритмы Дормана-Принса и Грэгга-Булирша-Штера), включает необходимые для анализа положений равновесия и предельных циклов алгоритмы определения собственных чисел. Пакет позволяет исследовать и консервативные системы ОДУ, используя в процессе численного интегрирования коррекцию Накози [3]. Но основным недостатком пакета является то, что практически все методы численного исследования инвариантных структур фазового пространства динамической системы доступны только для систем общего вида, и, следовательно, совершенно не используют сущностных свойств фазового потока канонических уравнений гамиль-тоновых систем. Продолжение семейств положений равновесия или периодических решений возможно только при условии, что параметр явно присутствует в правой части уравнений. Была предпринята попытка применения данного пакета к исследованию периодических решений задачи Хилла, но пакет смог лишь визуализировать проекции фазовой траектории на двумерные подпространства.
Пакет CONTENT был разработан Ю. Кузнецовым [4] и его группой в 90-х годах. В отличие от описанных выше средств он является открытым проектом и реализован на различных платформах, таких, как UNIX и Windows. Пакет предназначен для работы с динамическими система-
ми, задаваемыми в виде систем ОДУ, дифференциально-алгебраических уравнений, отображений и уравнений в частных производных первого порядка. Он содержит средства визуализации инвариантных структур и бифуркационных диаграмм, а также большой набор библиотек, реализующих алгоритмы численного интегрирования систем ОДУ, поиска и продолжения по параметру положений равновесия и предельных циклов, алгоритмы бифуркационного анализа. Работая совместно со средствами компиляции, этот пакет позволяет создавать пользователю высокоэффективные исполняемые модули, вычисляющие правые части систем ОДУ, не используя при этом навыки программирования последних. Предлагая оконный пользовательский интерфейс, пакет CONTENT обеспечивает интерактивный режим работы пользователя. Существенным недостатком данного пакета является то, что все алгоритмы поиска и продолжения инвариантных структур предназначены исключительно для работы с системами, правые части которых явно содержат внешние параметры. Пакет не использует первые интегралы системы ОДУ. Поэтому применить реализованные в нем алгоритмы для исследования, например, задачи Хилла не удалось. Построение отображений Пуанкаре средствами CONTENT возможно только для систем, заданных аналитически, но не для отображений, индуцированных непрерывным фазовым потоком.
Наиболее близким по архитектуре является пакет AUTO2000, разработанный Е. Доделем и др. Этот пакет является дальнейшим развитием пакета AUT094 [5], написанного на языке FORTRAN. Он поставляется в исходных текстах, предназначен для работы в среде UNIX, хотя удалось его без каких-либо изменений перенести на платформу Windows. Он позволяет исследовать алгебраические системы и системы автономных ОДУ с внешними параметрами, осуществлять продолжение периодических решений и применять к ним различные бифуркационные алгоритмы. Основной техникой исследования периодических решений является сведение последней к двухточечной краевой задаче на единичном интервале. Имеется возможность использовать программы пакета в среде кластерных вычислений. К сожалению, из-за отсутствия поддержки пакетом гамильтоновых систем и возможности использования библиотек высокоточной арифметики возникают трудности с применением AUTO2000 к исследованию интересующих автора небесно-механических моделей.
Возникла необходимость создания алгоритмов и программ для комплексного исследования гамильтоновых систем, позволяющих изучать разные объекты динамической системы, которые играют ключевую роль в ее динамике. Традиционные алгоритмы поиска и продолжения инвариантных динамических структур должны быть дополнены алгоритмами, позволяющими исследовать явления динамического хаоса и получать некоторые количественные характеристики последнего. Различные
проекты космонавтики, такие как низкоэнергетические переходы, движение специализированных космических аппаратов, требуют информации об инвариантных структурах небесно-механических задач, обладающих специфическими свойствами.
Целью работы является разработка пакета прикладных программ для комплексного исследования небесно-механических задач, использующего современные методы компьютерного исследования гамильто-новых динамических систем, и проведение средствами пакета изучения основных регулярных и хаотических структур конкретной небесно-механической модели — плоской круговой модели Хилла.
Первая глава носит методологический характер и содержит описание современных методов исследования автономных гамильтоновых систем.
В первом параграфе излагаются современные принципы исследования динамических систем в соответствии с работами [6,7,8,9,10, И, 12, 13,14,15].
В динамической системе есть объекты, играющие ключевую роль в ее динамике. Они образуют в некотором роде «скелет» поведения системы, и ее динамика будет достаточно понятна, если мы знаем эти объекты, знаем, как они связаны, знаем возможные переходы от одного к другому и как долго эти переходы могут реализовываться. К таким объектам относятся: неподвижные точки и периодические орбиты, инвариантные кривые, торы и другие инвариантные многообразия, устойчивые и неустойчивые инвариантные многообразия гиперболических объектов, пересечения устойчивых и неустойчивых многообразий (гомоклинические и гетероклинические).
Исследуя эти объекты, будем отдавать предпочтение параметрическому подходу к динамическим системам, который предполагает изучение зависимости всех найденных объектов от параметров, возможное продолжение по параметру и поиск бифуркаций. В диссертационной работе математический аппарат и инструментальные средства изучения выбирались соответствующими указанному подходу [16].
Сечение Пуанкаре устанавливает связь между динамическими системами с непрерывным и дискретным временем. Это особенно удобно для гамильтоновых систем с двумя степенями свободы на заданном уровне энергии. Пусть некоторая динамическая система задана гамильтонианом Н(х), х Є М4. В качестве «плоскости сечения» выберем подмногообразие Г С R4, полученное путём пересечения гиперплоскости с многообразием Н(х) = С для фиксированного значения параметра энергии С:Г = Н~г(С) ҐЇ7- Фазовый поток Ф* : R4 —> Е4, определённый гамильтонианом Н(х), корректно индуцирует отображение Пуанкаре Р : Г ^> Г при условии, что векторное поле гамильтониана Я трансверсально в каждой точке из Г. Обозначим через gr : Г — R4 отоб-
ражение вложения. Выберем z Є Г, тогда через t(z) будем обозначать наименьшее время, необходимое для возвращения фазовой траектории Фі (&г(і)) на плоскость сечения с тем же направлением вектора фазовой скорости. Тем самым отображение Р задаётся формулой
P(z) = g? (Фт(я) (gr(z))).
Определение отображения Пуанкаре гарантирует, что его предельные множества соответствуют предельным множествам исходного фазового потока [17].
Во втором параграфе рассматриваются основные методы численного определения отображения Пуанкаре для автономных гамильтоновых систем. Дискретная реализация непрерывной фазовой кривой требует уточнения точки пересечения последней с секущей гиперплоскостью. Такое уточнение может быть реализовано стандартной итерационной ньютоновой схемой [10]. Если в качестве секущей гиперплоскости выбрана одна из координатных плоскостей, например, *,-, то применяется метод, предложенный М. Непоп (см., например, [10]). В окрестности плоскости сечения дифференциальные уравнения движения заменяются уравнениями с новой независимой переменной — координатой х\. Тогда условие пересечения фазовой кривой с секущей гиперплоскостью выглядит наиболее просто: xi = щ. Однако сходимость вышеуказанных методов резко ухудшается, если условие трансверсальности становится плохо обусловленным. В этом случае можно применять вариант метода дихотомии. Поскольку на каждом шаге интегрирования известны фазовые координаты вместе с фазовыми скоростями, то время пересечения фазовой кривой с секущей гиперплоскостью эффективно определяется методом кубической интерполяции.
В третьем параграфе рассматриваются два подхода к поиску периодических решений (ПР) автономной гамильтоновой системы. Оба подхода сводят задачу поиска ПР к задаче поиска нулей некоторого векторного поля, то есть позволяют применять методы поиска положения равновесия, использующие ньютоновы итерации. Первый подход является разновидностью метода пристрелки, адаптированного для автономных систем ДУ. Периодическим решениям соответствуют нули векторного поля невязок .F(X), где X = (х,Т)
Т : Шп+1 -> Ш" : X н-> фт(х) - х.
Пусть известно начальное приближение Хо периодического решения, тогда вектор поправок ЛХ должен удовлетворять в линейном приближении матричному уравнению
(М-Е | f(tf>T(xo))).4X = -;F(Xo)
В этой системе число уравнений меньше числа неизвестных, поэтому, либо стандартная процедура ньютоновых интераций должна быть модифицирована (например, с использованием псевдообратной матрицы), либо должно быть добавлено дополнительное условие (например, условие изоэнергетической редукции или условие Мея [10]).
Основной недостаток метода связан с тем, что матрица монодромии М автономной гамильтоновой системы имеет два собственных числа равных +1, что делает матрицу М — Е особенной, поэтому эффективность численного решения системы, определяющей поправки, резко снижается.
Второй подход сводит проблему поиска периодического решения к проблеме поиска неподвижной точки к-й итерации отображения Пуанкаре. Условие трансверсальности фазового потока секущей гиперплоскости вместе с условием изоэнергетической редукции гарантируют невырожденность матрицы DPk — Е для всех точек сечения, кроме бифуркационных. Использование отображения Пуанкаре снижает размерность системы для определения поправок АХ, а отсутствие периода Т решения, как параметра, избавляет от необходимости накладывать на систему дополнительные условия. Тем не менее, неглобальность сечения Пуанкаре и наличие петель периодических орбит приводит либо к нарушению условия трансверсальности, либо к скачкообразному изменению порядка итерации к. Заметим, что этот параметр априорно не задан и должен выбираться из некоторого диапазона [18]. Более того, вследствие отсутствия для конкретных гамильтоновых систем аналитического способа задания отображения Пуанкаре, приходится и точки, и производные отображения строить численно путем интегрирования уравнения в вариациях.
Поиск периодического решения автономной системы общего вида без привлечения дополнительной информации достаточно сложен, поэтому обычно ему предшествует определенная аналитическая работа позволяющая установить наличие таких решений и определить области их существования. Возможен и так называемый метод «грубой силы», позволяющий путем интенсивных вычислений получить представление о структуре фазового потока в некоторых областях фазового пространства.
В четвертом параграфе обсуждаются методы продолжения семейств периодических решений гамильтоновых систем по параметру. Пусть известно некоторое периодическое (с периодом Т) решение исходной системы ОДУ х(Т,хо) = хо. Для нахождения близкого к нему периодического решения x(t) = хо(0 + у(0 с периодом Т + dT нужно решить систему уравнений
хо(Г + dT) + у(Г + dT) = хо(0) + у(0),
или, в линейном приближении
(М - Е)у(Г) + drjVH(xo(0)) = 0.
(*)
Ранг матрицы этой системы линейных уравнений всегда меньше числа неизвестных, так как для гамильтоновых систем матрица монодро-мии имеет по крайней мере два собственных числа, равных +1. Справиться с этой проблемой позволяет метод, предложенный Б.Б. Крей-сманом [19]. Этот метод предлагает переход в так называемый сопутствующий базис с помощью ортогональной и симплектической матрицы, элементами которой являются компоненты нормированного вектора VH. Структура матрицы (М — Е) в новом базисе позволяет описывать периодические решения вектором, лежащим в ортогональном дополнении к линейному подпространству, задаваемому касательным к решению хо вектором. Это позволяет избавиться от одного собственного значения +1, связанного с движением по циклу. Решая полученную систему уравнений, находим вектор, касательный к продолжаемому семейству, и поправки к периоду.
Зная нормированный вектор, касательный к ветви продолжения семейства, можем найти следующее периодическое решение этого семейства, пользуясь методом Келлера [4] продолжения семейства ПР по длине дуги. Особенностью метода продолжения по длине дуги является то, что матрица Якоби, используемая в ньютоновых итерациях, является невырожденной во всех регулярных точках бифуркационной диаграммы, в том числе в точках поворота («складках»).
Кроме того, метод Б.Б. Крейсмана позволяет определять направление продолжения не только данного семейства, но и всех взаимодействующих с ним семейств. Продолжая последние, находим направление ветвления взаимодействующих с ними семейств и так далее. Метод оказался чрезвычайно эффективным, чем и обусловлен его выбор для последующего построения алгоритмов исследования периодических решений.
Пятый параграф посвящен методам поиска и продолжения симметричных периодических решений. Для симметричных периодических решений стандартное условие периодичности вида х(0) = х(Г) может быть заменено условием ортогонального пересечения орбиты с осью симметрии через половину (для орбит с одной симметрией) или четверть (для орбит с двумя симметриями) периода [20]. Поэтому симметричное периодическое решение однозначно определяется трехмерным вектором X = (xi(0),Xj(0),T), і = 1,2, j = 3,4, и все формулы 5 значительно упрощаются [21]. Анализ решений системы ( * на предшествующей странице) позволяет сделать выводы о наличии симметрии у продолжаемых решений. Установлено, что окрестность дважды симметричного р/^-резонансного периодического решения устроена следующим образом: если числа р и q нечетные, то через него проходит одно дважды симметричное семейство (/-кратных периодических решений; если же хотя бы одно из чисел р или q четное, то на нем заканчиваются два семейства д-кратных
периодических решений с одной симметрией.
Методы исследования сценариев перехода к динамическому хаосу обсуждаются в шестом параграфе. Последовательность усложняющих структуру фазового пространства бифуркаций может приводить к возникновению хаотических режимов. Методы, описанные в двух предыдущих параграфах, позволяют определять наличие у семейства ПР бифуркации удвоения периода, а также каскадов этих бифуркаций и их количественных характеристик, таких, как постоянная Фейгенбаума [22,23] и масштабные константы [24]. Аналогично можно проанализировать каскады ^-кратного увеличения периода, например, для q = 3.
Возникновение расщепленных асимптотических поверхностей в окрестности седловых точек может быть причиной появления сложных динамических режимов. Определим, как это сделано в [17] устойчивые и неустойчивые инвариантные многообразия (сепаратрисы):
Ws'"(z0) = (z Є Г : lim Tn(z0) = z
где z0 - седловая точка отображения P. Сепаратрису можно реализовать как вложение вещественной прямой М в многообразие Г с помощью функций Yu,s: Ш — Г, удовлетворяющих условию
^(Л^О = Р(^()), Y"*(0) = z0,
где AM/S - собственные числа линейной части отображения Р в точке zo.
Представляя сепаратрисы в виде ряда по степеням малого параметра , который есть расстояние вдоль дуги сепаратрисы до гиперболической точки, получим бесконечную систему реккурентных линейных уравнений относительно коэффициентов этого разложения. Последовательное определение компонентов z"^, і = 1,2,... этих коэффициентов возможно при условии, что матрица (Е — DP(zJ's)) обратима, то есть AM/S ф ±1. Реккурентные уравнения, начиная со второго, содержат значения производных отображения Р, вычисленных в седловой точке. Сложность выражений для вычисления производных с увеличением их порядка значительно возрастает [17], поэтому в работе используется вычисление производных отображения через аппроксимацию последнего с помощью полиномов Чебышева.
В седьмом параграфе отмечены особенности небесно-механических задач как динамических систем. Поскольку гамильтонианы небесно-механических моделей содержат особенности в окрестности тяготеющих тел, то при интегрировании уравнений движения вблизи этих особенностей накапливаются большие ошибки. Наличие в семействе периодических решений столкновительных траекторий приводит к невозможности
продолжения этого семейства. Для преодоления этих трудностей предлагается использовать регуляризацию уравнений движения. Если глобальная регуляризация или невозможна, или сильно усложняет понимание динамики фазового пространства, то применяется локальная регуляризация Леви-Чивита [20] для уравнений движения и для уравнений в вариациях. В этом случае используется техника производящих функций, описанная в работе М. Лидова [25]. Также гамильтонианы большого числа небесно-механических моделей допускают некоторые виды симметрии, что делает возможным применение к их исследованию методов параграфа 1.5 на с. 42. Например, ограниченная задача трех тел (ОЗТТ) обладает одной симметрией относительно оси абсцисс [20].
Вторая глава содержит описание разработанного автором комплекса программ по исследованию динамических систем. Общая структура комплекса программ приведена в первом параграфе. Комплекс программ представляет собой набор исходных текстов программ и сценариев сборки исполняемых файлов для операционных систем, поддерживающих стандарт POSIX.1. Весь программный код комплекса и большая часть программ написаны на языках С и C++, а сценарии сборки комплекса написаны с использованием средств make. Наличие свободно распространяемых средств разработки открытых приложений GNU таких как эффективный компилятор С и C++, отладчик, загрузчик и другие, делает программный код комплекса мобильным. Объектно-ориентированные свойства языка программирования C++ позволяют перегружать стандартные арифметические операции, операции ввода-вывода и использовать один и тот же код алгоритма для различных классов данных, поддерживающих высокоразрядную арифметику. Для нормального функционирования программного комплекса программная среда должна удовлетворять стандартным требованиям, предъявляемым к разработке открытых программ в среде Linux.
В состав комплекса входят различные свободно распространяемые библиотеки и вспомогательные программные средства. Это позволяет свободно переносить комплекс с одной платформы на другую и дополнять его новыми компонентами по мере необходимости. Программы комплекса были испытаны на платформах Linux и Win32. Основной особенностью комплекса является то, что он ориентирован на использование высокоточной арифметики.
В структуре комплекса выделяются три уровня: библиотеки алгоритмов низкого уровня, которые содержат описание базовых типов данных и различных функций, используемых алгоритмами высокого уровня; вспомогательные программы и модули, которые обеспечивают получение начальных данных, их хранение и визуализацию; программы верхнего уровня, специализированные следующим образом: программы построе-
ния траекторий (орбит) и отображений Пуанкаре; программа поиска и продолжения периодических решений; программа поиска и продолжения неподвижных точек отображения; программа исследования сценариев перехода к динамическому хаосу.
В первом параграфе обсуждаются взаимосвязи различных компонент комплекса, а также описана общая последовательность использования программ комплекса с точки зрения пользователя.
Второй параграф посвящен программам численного интегрирования динамических систем. В качестве базового алгоритма интегрирования уравнений движения комплекс использует метод рядов Тейлора, реализованный в пакете Taylor версии 1.4. Выбор данного метода интегрирования связан со следующими особенностями: метод рядов Тейлора является методом интегрирования с автоматическим выбором шага и порядка разложения, для чего используется эффективный алгоритм вычисления производных любого порядка правых частей уравнений движения; метод поддерживает как стандартные типы данных, так и типы из библиотек высокоточной арифметики; сравнение этого метода с другими показало высокую его эффективность [26]. Для программ поиска и продолжения периодических решений необходимо вычислять матрицу монодромии или матрицу производных отображения Пуанкаре, а для исследования эффекта расщепления необходимо знание производных отображения высоких порядков. Определение этих объектов осуществляется либо путем интегрирования уравнений в вариациях совместно с уравнениями движения, либо средствами интерполяции. В первом случае алгоритм является строго последовательным, а во втором допускает несложное распараллеливание, что приводит существенному сокращению времени счета даже на высокоточной разрядной сетке.
Программы поиска и продолжения семейств периодических решений описаны в третьем параграфе. Поиск периодических решений может осуществляться или с применением метода «грубой силы», или с помощью ньютоновых итераций. В первом случае пользователю представляется возможность исследовать структуру фазового пространства путем построения и визуализации отображения Пуанкаре с помощью программы POINCAREMAP. Особенностью этой программы является то, что вычисления могут осуществляться как последовательно, так и параллельно (средствами вычислительного кластера). Если в качестве секущей гиперплоскости выступает одна из координатных плоскостей, то для вычисления отображения применяется метод М. Эно, иначе «работает» итерационная схема уточнений с использованием кубической интерполяции. Визулиазация отображений Пуанкаре позволяет определить начальные приближения для программ уточнения и продолжения семейств периодических решений DERPAR или CONTPER.
Для работы программ пользователь должен задать начальное приближение, погрешности интегратора и определения ПР, начальный шаг смещения вдоль семейства, максимальное число ньютоновых итераций. Программа DERPAR основана на известном алгоритме продолжения по параметру нулей векторного поля [12] и может применяться к системам общего вида, но обладает рядом недостатков. Во-первых, в процессе продолжения семейства может происходить скачок числа пересечений ПР с секущей плоскостью, что приводит к отказам в работе, а, во-вторых, алгоритм в силу своей универсальности не использует структуру матрицы производных отображения для определения направления ветвления ПР второго рода. Рекомендуется использование программы DERPAR для поиска и продолжения несимметричных периодических решений. Вторая из программ свободна от указанных выше недостатков, но является специализированной для поиска и продолжения семейств симметричных ПР. Следует отметить, что в обеих программах применяется адаптивный метод выбора шага продолжения, учитывающий число итераций, потраченных на уточнение ПР.
Программы исследования перехода к динамическому хаосу описаны в 2.4 на с. 72. Здесь дано описание трех программ: CASCADE - для поиска каскадов бифуркации удвоения периода, LYAPUNOV - для определения показателей Ляпунова и SEPSPLIT - для исследования явления расщепления сепаратрисы. Первая из программ основана на алгоритме CONTPER и позволяет, используя анализ структуры матрицы монодро-мии М, находить бифуркации удвоения периода, а также определять количественные показатели - постоянную Фейгенбаума и масштабные константы. Характерной особенностью алгоритма является использование высокоточной арифметики, которая обеспечивает необходимую точность вычисления ПР Алгоритм вычисления показателей Ляпунова использует стандартную технику ортогонализации Грамма-Шмидта [27] в процессе интегрирования уравнения в вариациях методом рядов Тейлора. Для изучения эффекта расщепления сепаратрис применяется несколько описанных выше алгоритмов, однако существенным является то, что для разложения сепаратрисы в окрестности гиперболической точки, требуется вычисление производных отображения Пуанкаре высоких порядков. В отличие от методики работы [17] в программе SEPSPLIT применяется аппроксимация отображения с помощью двумерных многочленов Чебы-шева. Для вычисления производных отображения Пуанкаре применяется дискретный аналог уравнения в вариациях, с помощью которого строятся интерполяционные многочлены Чебышева, обеспечивающие равномерную аппроксимацию на интервале и простую оценку погрешности аппроксимации. Для интенсификации вычислений может использоваться параллельный компьютинг. Узлы прямоугольной сетки рассматриваются
как независимые начальные условия задачи Коши, и последняя решается численным интегрированием на подчиненных узлах вычислительного кластера. Таким образом, параллелизм затрагивает самую трудоемкую часть алгоритма.
В третьей главе излагаются результаты изучения регулярных и хаотических режимов задачи Хилла, полученные с использованием комплекса программ, описанного во второй главе.
Задача Хилла имеет чрезвычайно интересную историю, связанную с именами Л. Эйлера, А. Пуанкаре, Дж. Хилла, A.M. Ляпунова. Ею занимались советские ученые. Эта задача явилась источником целого ряда ярких идей и замечательных достижений, составивших эпоху в развитии небесной механики.
Четкую математическую формулировку рассматриваемой задачи дал Хилл в своем знаменитом мемуаре «Исследования по теории Луны» [28]. Однако, ранее, Эйлер в «Новой теории Луны» составляет уравнения движения Луны в прямоугольных геоцентрических эклиптических координатах, равномерно вращающихся с угловой скоростью, равной среднему движению Луны. Он разлагает правые части уравнений в ряды, используя малые параметры (параллакс Солнца, наклон и эксцентриситет лунной орбиты), и последовательно «перебирая» все слагаемые правых частей, находит решение в периодических функциях. Эйлер показывает, как возмущения каждого класса могут быть найдены путем быстро сходящихся последовательных приближений [29].
Многие идеи Эйлера были развиты Хиллом далее. Вместе с тем он внес в рассматриваемую проблему много своих весьма оригинальных идей. Хилл пользуется прямоугольной геоцентрической эклиптической системой координат, равномерно вращающейся с угловой скоростью, равной среднему движению Солнца. Ось абсцисс он направляет по прямой, соединяющей Землю и Солнце. В этих координатах дифференциальные уравнения задачи имеют вид:
(х = 2у + 3х-^, \у=-2х-
где г — у/х2 + у2. Они имеют первый интеграл — интеграл Якоби
2 x2 + if = 3x2 + - -С,
где С — постоянная Якоби.
Модель Хилла представляет собой модель нулевого уровня в исследованиях сложных небесно-механических динамических систем (см.,
например, [20,30,31]). Занимая в некотором роде «промежуточное» положение между интегрируемой кеплеровои задачей и неинтегрируемои ограниченной задачей трех тел, гамильтониан задачи Хилла обладает лишь одной особенностью, как и гамильтониан первой задачи, но при этом задача Хилла является неинтегрируемои, как вторая (см. [32], [33]). Поэтому модель Хилла можно рассматривать как возмущение кеплеровои задачи, а ограниченную задачу трех тел — как возмущение задачи Хилла, представляющей начальный уровень сложности в иерархии моделей ограниченной задачи трех тел. Модель Хилла — это реальная динамическая система, которую можно использовать не только в теории Луны, но и в других спутниковых задачах, в которых отношение планетоцентрических расстояний спутника и внешнего тела является малой величиной, а орбита внешнего тела близка к круговой. Её можно применять и в исследовании тройных звездных систем, у которых одно из взаимных расстояний мало по сравнению с двумя другими. С помощью этой модели исследовалось движение звезд шарового скопления в поле галактики [34].
В первом параграфе обсуждаются общие свойства задачи Хилла. Показано, как уравнения этой задачи получаются из уравнений плоской ограниченной задачи трех тел (ОЗТТ) при стремлении массового параметра к нулю; записаны гамильтониан и канонические уравнения; доказано, что дифференциальные уравнения задачи Хилла эквивариант-ны относительно некоторой дискретной группы симметрии. Проведена, в соответствии с [20] и [25], регуляризация уравнений движения и гамильтониана. Рассмотрены положения равновесия задачи Хилла и условия их существования.
Второй параграф посвящен основным семействам периодических решений задачи Хилла. Эти семейства изучались ранее численно Дж. Хил-лом, Кельвином, Дж. Джексоном, Т. Матукумой, М. Эно [35]. Автором выполнены аналогичные численные исследования с более высокой точностью для подтверждения достоверности результатов, получаемых с помощью предложенного комплекса программ. В работе приведены графики устойчивости, характеристики для основных семейств и описана эволюция орбит каждого семейства при изменении константы Якоби. Ранее проводились аналитические исследования основных семейств [35,36,37] и др. Автором применен метод Депри-Хори [27] для аналитического исследования однооборотного семейства g и его бифуркаций [38,39].
В третьем параграфе предлагается некоторая классификация и подробно описаны свойства периодических решений второго рода задачи Хилла. Большинство из исследованных семейств являются новыми, в том смысле, что они не описаны ранее в работах других авторов. Для каждого из семейств построены графики устойчивости, характеристики,
показаны основные перестройки орбит семейства. Установлены различия между взаимодействием периодического решения второго рода с порождающим семейством с одной и двумя симметриями.
Обнаружены замкнутые периодические решения, для которых кривые устойчивости и характеристики являются кривыми, определенными на отрезке. Эти семейства либо дважды взаимодействуют с семейством /, либо взаимодействуют с семействами fug. Все они являются высокооборотными.
Выявлены возможные бифуркации [40] периодических решений второго рода и их роль в формировании динамического хаоса. К ним относятся: бифуркация рождения-гибели, бифуркация потери симметрии, бифуркация удвоения периода (кратного увеличения периода), транскритическая бифуркация. Для таких семейств определены бифуркационные значения параметра С с высокой точностью. Семейства прямооборотных периодических решений описаны в работах [41,42].
Четвертый параграф посвящен описанию каскадов бифуркаций удвоения периода некоторых периодических решений задачи Хилла, которые были впервые обнаружены автором в ходе численных экспериментов. Была построена бифуркационная последовательность для орбит с 1 _ 2 - 4 - 8 - 16 - 32 - 64 - 128 - 256 - 512 оборотами. Бифуркационные значения параметра С образуют сходящуюся геометрическую прогрессию, знаменатель которой достаточно быстро стремится к универсальной постоянной Фейгенбаума. Аналогичные расчеты проводились и для других орбит. Так, например, для семейства #^(1/5)(). неподвижная точка отображения которого лежит на оси симметрии, был построен каскад бифуркаций удвоения 5 — 10 — 20 — 40 — 80 —160 и вычислены не только постоянная 5, но и масштабные константы ос и р. Значения всех констант оказались близкими к теоретическим. При предельном значении С = Соо каждый каскад завершается появлением цикла бесконечного периода. Для семейства gf и его периодических решений второго рода эти циклы располагаются вблизи сепаратрисной поверхности, соответствующей семейству g. Наличие каскадов является причиной появления хаоса в окрестности неподвижных точек, соответствующих семейству g. Результаты этого параграфа изложены в работах [43,44].
Еще одной причиной возникновения локального хаоса является расщепление инвариантных многообразий гиперболических точек отображения Пуанкаре. Исследованию этого явления в задаче Хилла посвящен 5-ый параграф 3-ей главы. Сначала было получено расщепление сепаратрисы неподвижной гиперболической точки отображения, соответствующей семейству g, а затем и точек, соответствующих семействам /(1/9)(). /(1/7)(). /(1/8)(). /(1/6)(). Были определены численно такие количественные характеристики расщепления, как пло-
щадь луночки и гомоклинический инвариант, а также исследована их зависимость от собственных чисел матрицы производной отображения в неподвижной гиперболической точке. Проведенное исследование позволяет утверждать, что при значениях параметра, близких к бифуркационному, наблюдается эффект экспоненциальной малости расщепления сепаратрис. Результаты этого параграфа опубликованы в работах [45], [46]. Гомоклинические и гетероклинические пересечения инвариантных многообразий, соответствующих ляпуновским периодическим орбитам а и с, изучены подробно К.Симо в работе [47], поэтому они в данной работе не рассматривались.
Последний параграф настоящей диссертации содержит описание глобальной динамики задачи Хилла при различных значениях интеграла энергии, полученное на основе изложенных выше результатов. Изложение, в основном, следует работам автора [48,49,50].
В конце диссертации приводится список основных обозначений и список литературы, содержащий 67 наименований.
По теме диссертации опубликовано 19 работ [48,49,51,52,53,43,38, 54,44,39,45,55,46,50,41,21,56,57,42].
Результаты диссертации докладывались и обсуждались на международной конференции «Астрометрия, геодинамика и небесная механика на пороге XXI века» (Санкт-Петербург, 2000 г.); на всероссийской научной конференции «Новые результаты аналитической и качественной небесной механики» (Москва, 2000 г.); на международной конференции «Celestial Mechanics» (Санкт-Петербург, 2002 г.); на международной конференции «Колмогоров и современная математика» (Москва, 2003 г.); на международном конгрессе «New Geometry of Nature» (Казань, 2003 г.); на всероссийской астрономической конференции «Горизонты Вселенной» (Москва, 2004 г.); на семинарах отдела небесной механики Государственного астрономического института им. П.К. Штернберга МГУ (Москва) в 2003 и 2005 г.г.; на семинаре «Механика. Управление. Информатика.» Института космических исследований РАН в 2005 г.; на семинаре им. В.А. Егорова по динамике космического полета при кафедре теоретической механики МГУ в 2005 г.; на семинарах кафедры математического анализа и теории функций Вол ГУ в 2003-2004 г.г.; на семинаре кафедры экспериментальной математики и информатики Вол ГУ в 2004 г.; на научных конференциях профессорско-преподавательского состава ВГИ Вол ГУ в 1997-2004 г.; на других международных и всероссийских конференциях.
Пользуясь случаем, автор выражает глубокую благодарность к.ф.-м.н. Сумарокову СИ. за руководство работой, а также к.т.н. Крейсма-ну Б.Б. за постоянное внимание и поддержку в работе.
Поиск периодических решений и анализ их устойчивости
Поэтому модель Хилла можно рассматривать как возмущение кеплеровои задачи, а ограниченную задачу трех тел — как возмущение задачи Хилла, представляющей начальный уровень сложности в иерархии моделей ограниченной задачи трех тел. Модель Хилла — это реальная динамическая система, которую можно использовать не только в теории Луны, но и в других спутниковых задачах, в которых отношение планетоцентрических расстояний спутника и внешнего тела является малой величиной, а орбита внешнего тела близка к круговой. Её можно применять и в исследовании тройных звездных систем, у которых одно из взаимных расстояний мало по сравнению с двумя другими. С помощью этой модели исследовалось движение звезд шарового скопления в поле галактики [34].
В первом параграфе обсуждаются общие свойства задачи Хилла. Показано, как уравнения этой задачи получаются из уравнений плоской ограниченной задачи трех тел (ОЗТТ) при стремлении массового параметра к нулю; записаны гамильтониан и канонические уравнения; доказано, что дифференциальные уравнения задачи Хилла эквивариант-ны относительно некоторой дискретной группы симметрии. Проведена, в соответствии с [20] и [25], регуляризация уравнений движения и гамильтониана. Рассмотрены положения равновесия задачи Хилла и условия их существования.
Второй параграф посвящен основным семействам периодических решений задачи Хилла. Эти семейства изучались ранее численно Дж. Хил-лом, Кельвином, Дж. Джексоном, Т. Матукумой, М. Эно [35]. Автором выполнены аналогичные численные исследования с более высокой точностью для подтверждения достоверности результатов, получаемых с помощью предложенного комплекса программ. В работе приведены графики устойчивости, характеристики для основных семейств и описана эволюция орбит каждого семейства при изменении константы Якоби. Ранее проводились аналитические исследования основных семейств [35,36,37] и др. Автором применен метод Депри-Хори [27] для аналитического исследования однооборотного семейства g и его бифуркаций [38,39].
В третьем параграфе предлагается некоторая классификация и подробно описаны свойства периодических решений второго рода задачи Хилла. Большинство из исследованных семейств являются новыми, в том смысле, что они не описаны ранее в работах других авторов. Для каждого из семейств построены графики устойчивости, характеристики, показаны основные перестройки орбит семейства. Установлены различия между взаимодействием периодического решения второго рода с порождающим семейством с одной и двумя симметриями.
Обнаружены замкнутые периодические решения, для которых кривые устойчивости и характеристики являются кривыми, определенными на отрезке. Эти семейства либо дважды взаимодействуют с семейством /, либо взаимодействуют с семействами fug. Все они являются высокооборотными.
Выявлены возможные бифуркации [40] периодических решений второго рода и их роль в формировании динамического хаоса. К ним относятся: бифуркация рождения-гибели, бифуркация потери симметрии, бифуркация удвоения периода (кратного увеличения периода), транскритическая бифуркация. Для таких семейств определены бифуркационные значения параметра С с высокой точностью. Семейства прямооборотных периодических решений описаны в работах [41,42].
Четвертый параграф посвящен описанию каскадов бифуркаций удвоения периода некоторых периодических решений задачи Хилла, которые были впервые обнаружены автором в ходе численных экспериментов. Была построена бифуркационная последовательность для орбит с 1 _ 2 - 4 - 8 - 16 - 32 - 64 - 128 - 256 - 512 оборотами. Бифуркационные значения параметра С образуют сходящуюся геометрическую прогрессию, знаменатель которой достаточно быстро стремится к универсальной постоянной Фейгенбаума. Аналогичные расчеты проводились и для других орбит. Так, например, для семейства # (1/5)(). неподвижная точка отображения которого лежит на оси симметрии, был построен каскад бифуркаций удвоения 5 — 10 — 20 — 40 — 80 —160 и вычислены не только постоянная 5, но и масштабные константы ос и р. Значения всех констант оказались близкими к теоретическим. При предельном значении С = Соо каждый каскад завершается появлением цикла бесконечного периода. Для семейства gf и его периодических решений второго рода эти циклы располагаются вблизи сепаратрисной поверхности, соответствующей семейству g. Наличие каскадов является причиной появления хаоса в окрестности неподвижных точек, соответствующих семейству g. Результаты этого параграфа изложены в работах [43,44].
Еще одной причиной возникновения локального хаоса является расщепление инвариантных многообразий гиперболических точек отображения Пуанкаре. Исследованию этого явления в задаче Хилла посвящен 5-ый параграф 3-ей главы. Сначала было получено расщепление сепаратрисы неподвижной гиперболической точки отображения, соответствующей семейству g, а затем и точек, соответствующих семействам /(1/9)(). /(1/7)(). /(1/8)(). /(1/6)(). Были определены численно такие количественные характеристики расщепления, как площадь луночки и гомоклинический инвариант, а также исследована их зависимость от собственных чисел матрицы производной отображения в неподвижной гиперболической точке. Проведенное исследование позволяет утверждать, что при значениях параметра, близких к бифуркационному, наблюдается эффект экспоненциальной малости расщепления сепаратрис. Результаты этого параграфа опубликованы в работах [45], [46]. Гомоклинические и гетероклинические пересечения инвариантных многообразий, соответствующих ляпуновским периодическим орбитам а и с, изучены подробно К.Симо в работе [47], поэтому они в данной работе не рассматривались.
Последний параграф настоящей диссертации содержит описание глобальной динамики задачи Хилла при различных значениях интеграла энергии, полученное на основе изложенных выше результатов. Изложение, в основном, следует работам автора [48,49,50].
Результаты диссертации докладывались и обсуждались на международной конференции «Астрометрия, геодинамика и небесная механика на пороге XXI века» (Санкт-Петербург, 2000 г.); на всероссийской научной конференции «Новые результаты аналитической и качественной небесной механики» (Москва, 2000 г.); на международной конференции «Celestial Mechanics» (Санкт-Петербург, 2002 г.); на международной конференции «Колмогоров и современная математика» (Москва, 2003 г.); на международном конгрессе «New Geometry of Nature» (Казань, 2003 г.); на всероссийской астрономической конференции «Горизонты Вселенной» (Москва, 2004 г.); на семинарах отдела небесной механики Государственного астрономического института им. П.К. Штернберга МГУ (Москва) в 2003 и 2005 г.г.; на семинаре «Механика. Управление. Информатика.» Института космических исследований РАН в 2005 г.; на семинаре им. В.А. Егорова по динамике космического полета при кафедре теоретической механики МГУ в 2005 г.; на семинарах кафедры математического анализа и теории функций Вол ГУ в 2003-2004 г.г.; на семинаре кафедры экспериментальной математики и информатики Вол ГУ в 2004 г.; на научных конференциях профессорско-преподавательского состава ВГИ Вол ГУ в 1997-2004 г.; на других международных и всероссийских конференциях.
Исследование сценариев перехода к динамическому хаосу
Данная формула показывает, что вектор bi имеет нулевые проекции на оси Х2 и хз, а вектор Ьг на оси х\ и х±, соответственно, а, значит, если хотя бы один из миноров d2 з отличен от нуля, то продолжение симметричного относительно оси х\ или Xi решения дает симметричное периодическое решение.
Для симметричных периодических решений стандартное условие периодичности вида х(0) = х(Т) может быть заменено в зависимости от типа симметрии одним из следующих условий. 1. В случае орбит симметричных относительно оси х\ с начальным условием х2(0) = з(0) = 0 имеем х2(Г/2) = Хз(Т/2) = 0. 2. В случае орбит симметричных относительно оси х2 с начальным условием i(0) = х4(0) = 0 имеем Xi(T/2) = х4(Г/2) = 0. 3. В случае орбит симметричных относительно осей Х\ и Х2 с начальным условием 2(0) = з(0) = 0 имеем Xi(T/4) = 3:4(Г/4) = 0. 4. В случае орбит симметричных относительно осей i и x2 с начальным условием Xi(0) = 4(0) = 0 имеем л:2(Г/4) = Хз(Т/4) = 0. Эти условия означают, что периодическая орбита ортогонально пересекает соответствующую ось симметрии. В случае одной симметрии это одна и та же ось, а в случае двух симметрии это разные оси. Таким образом, симметричное периодическое решение однозначно определяется трехмерным вектором X = (хі(0),Х4(0),Т) в случаях 1 и 3, или X = (x2(0),X3(0), Г) в случаях 2 и 4. Будем рассматривать двояко симметричные периодические орбиты, которые последовательно ортогонально пересекают оси симметрии через четверть периода. Пусть есть порождающая двояко симметричная периодическая орбита с периодом Г, в окрестности которой появляется периодическая орбита 2 рода с периодом qT. Для того, чтобы эта орбита имела также две симметрии, необходимо выполнение указанного выше условия ортогональности. Пусть начальная точка порожденной орбиты лежит на оси ОХ, тогда через четверть периода следующая ортогональная точка будет лежать на оси ОУ, при этом точка орбиты должна совершить (2к —1)/4 оборота вокруг начала координат, где к Є N. Поэтому за весь период точка орбиты совершит 2к — 1 оборотов, то есть число q всегда нечетное. С другой стороны, орбита, соответствующая резонансу p/q, может быть представлена в фазовом пространстве, как замкнутая винтовая линия, лежащая на некотором торе [22], которая делает за один период q оборотов по одной образующей тора и р оборотов по другой образующей. Лоэтому для дважды симметричной орбиты число р тоже должно быть «четным. Вывод: окрестность дважды симметричного / -резонансного периодического решения устроена следующим образом: если числа р и q нечетные, то через него проходит одно дважды симметричное семейство -кратных периодических решений; если же хотя бы одно из чисел р или q четное, то на нем заканчиваются два семейства -кратных одно-симметричных периодических решений. Выясним, как наличие симметрии может быть использовано для упрощения процедуры продолжения неподвижных точек отображения. Преобразование S : R — Ш называется преобразованием симметрии отображения Р, если [24] S2 = (Р о S)2 = I, где I — тождественное преобразование, и detS(z) = 1. Неподвижные точки преобразования S образуют линию, называемую линией симметрии. Если Р обладает симметрией, то Р — генератор циклической группы {Рк}, к Є Z, причем Р-1 = SoPoS и S — преобразование симметрии для любого Рк. Если неподвижная точка /с-той итерации отображения Р лежит не на линии симметрии, то можно провести некоторый содержательный анализ. Не нарушая общности, будем считать, что линией симметрии является одна из координатных осей, для определенности — ось абсцисс, то есть S : (z\,Z2) — (zi, —22). Если точка z лежит на оси симметрии и является неподвижной точкой отображения Рк, то, как показано в [66], матрица DzPk(z ) имеет вид ( 1, а условие продолжения в линейном приближении имеет вид Если а ф 1, то be ф 0, и из второго и четвертого уравнений системы ( 1.59) следует, что Az2 = 0, то есть при продолжении по параметру симметрия сохраняется. Устойчивость точки 2 эквивалентна условию \а\ 1. Пусть точка z устойчива, тогда a = coso , а, в силу того, что отображение Р сохраняет площадь, be = cos2 со — 1 = — sin2 со. Непосредственными вычисленнями проверяется, что Поскольку отображение Пуанкаре индуцировано фазовым потоком, то наличие симметрии у Р связано с симметрией канонических уравнений. Если неподвижная точка отображения лежит на оси симметрии, то орбита ортогонально пересекает эту ось. Как показано в [20], существует вторая точка ортогонального пересечения орбиты с осью, то есть на оси симметрии есть еще одна неподвижная точка. Будем обозначать эти точки z , а элементы матрицы DzPk(z ) — а,-, Ь,-, с;, dt, і = 1,2. Рассмотрим случай а = 1, соответствующий некоторой бифуркации. Тип бифуркации может быть определен лишь с учетом дополнительной информации. В силу условия сохранения площади при а = 1 либо Ъ = 0, либо с = 0. 1) Если / ф 0, то вектор, касательный к характеристике, имеет направление (1,0,0), то есть в этой точке наблюдается экстремум по параметру а и одновременно смена устойчивости. Такая бифуркация имеет несколько названий: бифуркация рождения-гибели; седло-узел; бифуркация типа складки. При этом ранг матрицы DP(z ,# ) по прежнему равен 2, поэтому система () может быть разрешима, если в качестве свободной переменной выбрать переменную Z\. 2) Если / = 0, то rangDP(z ,a ) = 1, что соответствует пересечению двух семейств. При этом возможны следующие ситуации: либо одно из семейств имеет экстремум по ос (е ф 0, бифуркация типа «вилки»), либо происходит обмен устойчивостями между этими семействами (транскритическая бифуркация).
Пусть теперь Ъ = 0, а с ф 0. Тогда Агг ф 0, то есть происходит уход с оси симметрии. Появляется пара несимметричных относительно оси 22 = 0 решений, которые взаимно симметричны относительно этой оси.
Используя информацию о структуре матрицы монодромии [21], можно показать, что матрицы DzP(z] /a ) и DzP(z2,a ) устроены идентично, поэтому результаты бифуркационного исследования не зависят от того, в какой из двух точек ортогонального пересечения с осью симметрии он проводится.
Программы поиска и продолжения периодических решений
Метод продолжения семейства симметричных периодических решений требует вычисления матрицы монодромии путем интегрирования уравнения в вариациях ( 1.24 на с. 32). Для орбит соударения удобно заменить уравнение в вариациях его дискретным аналогом. Для этого зададим вариацию начального условия периодического решения в виде набора одномерных сеток с к узлами вдоль каждой из координатных осей. Поскольку вариация изохронная, то каждый из 4к узлов сетки задает начальное условие задачи Коши для уравнения ( 1.1 на с. 18). Решая эти задачи на интервале, соответствующем типу симметрии (Т/2 для однократно симметричных, Т/4 - для двукратно симметричных орбит), и, применяя соответствующую разностную схему, получим значения элементов матрицанта Y(), и, следовательно, матрицу монодромии М. Заметим, что такой способ вычисления может быть легко распараллелен на 4к независимых процессов, что существенно повышает скорость счета без изменения процедуры интегрирования системы уравнений [55].
Если в процессе продолжения семейства одна из точек ортогонального пересечения оси орбитой попадает в окрестность особой точки гамильтониана, то даже регуляризация не позволяет корректно определить матрицу монодромии. Была опробована адаптивная схема, позволяющая выбирать такое условие периодичности, чтобы и начальная и конечная точки находились вне окрестности особой точки. Продемонстрируем идею на следующем примере. Пусть периодическое решение обладает симметрией ( 1.47 на с. 42), то есть в качестве начальной точки Л может быть выбрана точка с координатами (xf, 0,0, х). Пусть для этого условия выполнено i Rsafe- В силу условия симметрии через интервал Т/2 в точке В орбита будет иметь координаты ( f,0,0,#4)-Если xf Rsafe, то при определении элементов матрицы М возможны большие погрешности, которые не позволят корректно продолжить семейство. Заменив условие периодичности на полупериоде Т/2 условием xf{T) = з СП = 0 на периоде, обеспечим определение элементов матрицы М, а, значит, и параметров продолжения семейства с необходимой точностью. При этом почти вдвое увеличиваются вычислительные затраты. Если по мере продолжения семейства точка, соответствующая половине периода, покидает окрестность особенности, то можно вернуться к более экономной схеме, основанной на условии периодичности на половине периода.
Одним из признаков сложной структуры фазового потока является наличие бесконечных каскадов кратного увеличения периода семейств ПР. Доказательство наличия бесконечного каскада численными методами невозможно. Но, если характерные константы оказываются достаточно близкими к теоретическим, то можно сделать вывод о существовании в изучаемой динамической системе такого явления.
Программа CASCAD в своей работе существенно использует алгоритм продолжения семейства ПР с одновременным мониторингом ветвления ПР второго рода с периодом qT. Для функционирования программы необходимо задать начальные условия порождающего симметричного ПР, начальный шаг смещения вдоль характеристики семейства и шаг смещения при поиске ответвившегося решения. На первом этапе алгоритм продолжает порождающее ПР до тех пор, пока его индекс устойчивости не станет равным cos(2np/q). По уже вычисленным точкам семейства с помощью рациональной аппроксимации определяются координаты точки ветвления периодического решения с соизмеримостью p/q и направление ветвления, то есть формируются начальные условия для продолжения ответвившегося семейства. Одновременно фиксируются параметры, необходимые для численного определения постоянной Фейгенбаума и скейлинговых констант. Согласно формулы ( 1.61 на с. 48), необходимо обнаружить не менее трех бифуркаций кратного увеличения периода. С учетом того, что при каждой новой бифуркации ветвление происходит с измельчением масштаба по параметру и с увеличением периода ПР, приходится пропорционально увеличивать точность интегрирования и уменьшать шаг смещения вдоль характеристики. Полностью процесс определения каскада сложно автоматизировать, работа с программой выглядит, как метод «проб и ошибок». Когда в процессе продолжения каскада достигнута предельная точность соответствующей разрядной сетки, то можно воспользоваться вариантом программы, откомпилированной для двойной разрядной сетки. Если же и такой точности недостаточно, то можно откомпилировать программу для «четверной» разрядной сетки.
Следует отметить, что численное определение скейлинговых констант а и /5 реализовано только для каскадов удвоения периода, а вычисление константы Фейгенбаума 5 и предельного значения параметра Соо возможно для любых каскадов, в том числе и сложных, например 1/2/3/2/3/2....
Для исследования каскадов удвоения периода можно использовать совместно две другие программы комплекса DERPAR и POINMAP. Первая позволяет продолжать по параметру неподвижную точку отображения соответствующей кратности, а вторая позволяет визуализировать сечение Пуанкаре в окрестности каскада, чтобы наблюдать появление новых неподвижных точек удвоенной кратности. С помощью последней программы возможна визуализация некоторого приближения нестранного хаотического аттрактора Фейгенбаума — цикла бесконечного периода, соответствующего значению С ».
Программа исследования эффекта расщепления сепаратрисы SEPSPLIT. Численное изучение эффекта расщепления инвариантных многообразий, связанных с неподвижной гиперболической точкой отображения, реализуется с применением техники отображений Пуанкаре. Существенным условием успешного применения программы SEPSPLIT является возможно более точное определение координат неподвижной гиперболической точки отображения, лежащей на оси симметрии. Это позволяет с минимальными затратами определить положение гомокли-нической точки, в которой трансверсально пересекаются устойчивая и неустойчивая ветви сепаратрисной поверхности, а также вычислить такие инварианты расщепления, как площадь луночки и гомоклиническии инвариант.
Алгоритм работы программы предполагает несколько этапов. На первом этапе происходит определение с высокой точностью неподвижной гиперболической точки отображения (то есть точки, в которой собственные числа матрицы производных отображения DP вещественны). Затем в этой точке вычисляются производные отображения первого, второго и третьего порядков, причем эти вычисления могут быть выполнены двумя способами: либо путем интегрирования уравнений в вариациях соответствующих порядков ( 1.10 на с. 25), ( 1.11 на с. 25), ( 1.12 на с. 25), либо путем аппроксимации отображения с помощью многочленов Чебышева. Такая аппроксимация осуществлялась путем расчета коэффициентов полиномов по образам узлов некоторой разностной сетки, построенной в окрестности неподвижной точки отображения. Поскольку итерации соседних узлов сетки могут выполняться независимо, то этот вычислительный процесс легко распараллеливается. Полиномы Чебышева дают на интервале равномерную аппроксимацию и позволяют легко определять точность приближения. Если точность приближения недостаточна, то увеличивается количество узлов сетки.
Классификация периодических решений второго рода
Оказалось, что семейства g[r(l/2)r(Zi), gJr(l/3)(Zi), также, как и g7, вновь становятся устойчивыми при отрицательных значениях С. Индекс устойчивости семейства r(l/3)(Zi) имеет минимум « —2-Ю12, но затем начинает монотонно расти и при С s=s —4,69 семейство вновь становится устойчивым.
Итак, при С « —4,69217 индекс устойчивости семейства % вновь становится равным —1, и происходит бифуркация удвоения периода, однако по типу II (см. [69]). То есть, само семейство g в точке бифуркации «сменило» неустойчивость на устойчивость, а в его окрестности появляется неустойчивое решение вдвое большего периода. Поэтому при С О не наблюдается каскадов бифуркаций удвоения периода. Появление и поведение других резонансов при отрицательных С отличается от их появления и поведения при положительных значениях константы. Прежде всего, все существенные изменения индекса устойчивости происходят в очень узкой полосе изменения постоянной Якоби. Резонанс 1/4 уже не является сильным, его появление при С = —4,6984690 совпадает с достижением максимума по С вдоль семейства. Вдоль неустойчивой ветви при уменьшении С индекс устойчивости монотонно и быстро растет. Вдоль устойчивой ветви индекс устойчивости сначала уменьшается монотонно, затем достигает минимума s » —1,057, затем максимума s « —0,935, а затем еще раз — минимума s : —8,869, после чего индекс устойчивости монотонно и быстро возрастает. Таким образом, при изменении С от —4,702 до —4,7 индекс устойчивости четырежды проходит значение —1, и каждый раз происходит бифуркация удвоения периода. Резонансы 1/5, 1/7 и т. д. ведут себя аналогично, с той лишь разницей, что вдоль устойчивой ветви индекс устойчивости достигает минимума один раз, после чего резко растет. Резонанс 1/3 сохраняет свойства сильных ре-зонансов, его появление и поведение совпадает с поведением семейства g[r(l/3)(Zi) при положительных С.
Далее, при С = —4,70476 индекс устойчивости семейства gf проходит значение +1, само семейство становится неустойчивым, а в его окрестности появляется семейство gflt (1/1), уже не обладающее симметрией относительно оси ОХ. Впервые это семейство обнаружил Эно [18].
Поскольку / обладает двумя симметриями, появление резонансов p/q происходит по-разному в зависимости от четности чисел р и q. Для случая, когда либо р, либо q четное, все периодические решения второго рода обладают только одной симметрией. При этом, от пары ортогональных точек на оси ОХ ответвляется пара неустойчивых -оборотных орбит, сохранивших симметрию относительно ОХ, и симметричных друг другу относительно оси OY. А от ортогональных точек пересечения орбиты / с осью OY ответвляются две устойчивые {/-оборотные орбиты, сохранившие симметрию относительно оси ОУ, и взаимно симметричные относительно оси ОХ. (Слова «устойчивые» и «неустойчивые» для некоторых орбит поменяются местами). Так как для взаимно симметричных семейств все изменения происходят синхронно, то при их описании ограничимся одной ветвью с симметрией L\ и одной ветвью с симметрией Т.2 Для случая, когда и р и q нечетные, все порожденные решения двукратно симметричны и могут быть продолжены от ортогональных точек, расположенных на любой из осей. При этом от одной из ортогональных точек каждой оси ответвляется устойчивая ветвь -оборотного решения, а от другой — неустойчивая. Кроме того для интервалов i\ и 1 2 ветвление локально одинаково, только для орбит, ответвившихся в интервале І2, в момент появления достигается минимум по С, а для орбит, ответвившихся в интервале i\ — с достижением максимума по С.
Рассмотрим периодические решения второго рода, ответвившиеся от / в интервале 2 с p/q 1/3. Можно утверждать, что неустойчивые ветви слабых резонансов имеют много общего в изменении их индекса устойчивости. Прежде всего, ветвление устойчивой и неустойчивой ветвей происходит с достижением минимума по С. Затем вдоль неустойчивой ветви индекс устойчивости растет, достигает максимума, уменьшается, проходит значение +1 с одновременным достижением максимума по С, проходит значение —1, уменьшается до некоторого минимума, затем монотонно возрастает, опять проходя все значения от —1 до +1. Причем, когда индекс становится равен —1, происходит бифуркация удвоения периода по типу I, то есть с возможным образованием каскадов. В момент прохождения индекса устойчивости каждого семейства через +1 от него ответвляется семейство, потерявшее симметрию относительно оси ОХ. Это ответвившееся семейство во всех исследованных случаях является порождающим для несимметричных орбит. Все существенные изменения индекса устойчивости происходят после того как нечетнооборотные орбиты испытали соударения и стали прямооборотными. Поэтому в сечении Пуанкаре мы наблюдаем очень сложное поведение в окрестности этих резонансов в области, соответствующей прямооборотным орбитам. В поведении устойчивых ветвей нечетнооборотных семейств есть существенные различия, мы укажем на них при описании каждого отдельного семейства. К общим свойствам этих семейств следует отнести то, что после прохождения соударений количество их оборотов вокруг начала координат уменьшается на два.
Для слабых резонансов интервала i\ характерно другое поведение. Если р и q нечетные, то от / ответвляются двоякосимметричные q-оборотные орбиты, причем локально это происходит также, как в интервале І2, только теперь в момент ветвления по С достигается минимум. Устойчивая ветвь каждого из резонансов имеет минимум индекса устойчивости, равный —1, то есть мультипликаторы каждой орбиты проходят по всей единичной окружности, и происходит локальное ветвление резонансов со всеми соизмеримостями p/q. Если р и q нечетные, то ветвятся двоякосимметричные семейства, если хотя бы одно из чисел р или q четное, то ветвление происходит с потерей симметрии. При S = +1 происходит бифуркация потери симметрии, а именно теряется симметрия относительно оси ОУ. Устойчивая ветвь ответвившегося семейства с симметрией Ei также имеет минимум индекса устойчивочти s — 1, поэтому кривая устойчивости дважды проходит —1, причем бифуркация удвоения периода происходит в одном случае с образованием каскада, а в другом — без. При прохождении индекса устойчивости через +1 происходит ответвление несимметричного семейства.
