«Численно-экспериментальное исследование вращательной динамики спутников планет ...»

Главная (Пулковская) астрономическая обсерватория РАН

На правах рукописи

Куприянов Владимир Викторович

Численно-экспериментальное исследование

вращательной динамики спутников планет

01.03.01 – Астрометрия и небесная механика

ДИССЕРТАЦИЯ

на соискание ученой степени

кандидата физико-математических наук

Научный руководитель д. ф.-м. н.

Шевченко Иван Иванович Санкт-Петербург – 2014 Оглавление Введение................................... 4 Глава 1. Исторический обзор..................... 1.1. Численное моделирование в задаче тел............ 1.2. Численный эксперимент и динамический хаос.......... 1.3. Методы компьютерной алгебры в небесной механике и динами­ ческой астрономии.......................... 1.4. Выводы к первой главе....................... Глава 2. Хаотическое вращение спутников планет: ляпуновские спектры и максимальные характеристические показатели Ля­ пунова................................... 2.1. Введение............................... 2.2. Основные определения и алгоритмы................ 2.3. Аналитическое оценивание МХПЛ................. 2.4. Численные методы определения полных ляпуновских спектров и МХПЛ............................... 2.5. ХПЛ хаотического вращения спутников планет......... 2.6. Сравнение случаев плоского и пространственного вращения.. 2.7. Выводы ко второй главе...................... Глава 3. Хаотическое вращение спутников планет: ляпуновские спектры и константа Якоби..................... 3.1. Введение............................... 3.2. Вычисление ляпуновских спектров................ 3.3. Орбитальное движение и ХПЛ вращения............. 3.4. Аналитическое оценивание ХПЛ: границы применимости... 3.5. Зависимость ХПЛ от значения константы Якоби........ 3.6. Точность вычисления компонент ляпуновского спектра..... 3.7. Выводы к третьей главе....................... Глава 4. Вращательная динамика спутников планет: обзор ре­ гулярного и хаотического поведения................ 4.1. Введение............................... 4.2. Постановка численного эксперимента............... 4.3. Угловые скорости и ляпуновские времена............ 4.4. Устойчивость движения в синхронном резонансе........ 4.5. Приливное замедление вращения................. 4.6. Выводы к четвертой главе..................... Заключение.................................. Список литературы............................ Приложение А. Программный комплекс для вычисления ля­ пуновских спектров динамических систем на основе HQRB­ метода................................... А.1. Обозначения............................. А.2. Структура комплекса........................ А.3. Использование комплекса...................... А.4. Подключение интегратора, отличного от DOP853........ А.5. Вопросы переносимости....................... Введение Актуальность работы К настоящему времени в Солнечной системе открыто уже более 170 спут­ ников планет [71]. Значительный материал об их орбитальных и физических характеристиках накоплен в результате как наземных наблюдений, так и кос­ мических миссий («Вояджер-1», «Вояджер-2», «Галилео», «Кассини»). Зада­ чи динамики спутников и спутниковых систем, формирования их современ­ ных динамических состояний являются одними из актуальнейших в совре­ менной небесной механике и космогонии Солнечной системы. Орбитальная и вращательная динамика спутников связана с их физическими свойствами — массой, размерами, формой, составом и внутренним строением — и, таким образом, имеет важное значение в планетологии.

В результате исследований, выполненных в последние три десятилетия, стало ясно, насколько существенную роль в динамике Солнечной системы — и, в частности, в динамике спутниковых систем — играют резонансные явле­ ния (см. напр. книгу Мюррея и Дермотта [11]). Во многих случаях резонансы определяют пространственную конфигурацию орбит спутников и структуру колец планет. Многие из известных естественных спутников находятся в на­ стоящее время в состоянии синхронного спин-орбитального резонанса, про­ цесс захвата в который является важным событием в динамической истории спутника. Детали этого процесса, так же как и многих других эффектов, связанных с резонансами, все еще остаются мало изученными.

Взаимодействие резонансов порождает фундаментальный динамический эффект — хаотическое поведение. Уиздом и др. [70] в 1984 г. на основе ана­ лиза возможности существования основных резонансных спин-орбитальных состояний и их устойчивости у известных спутников планет сделали вывод, что вращение 7-го спутника Сатурна Гипериона должно быть хаотическим.

Позднее этот вывод был подтвержден в наблюдениях Клаветтером [42], Бл­ эком и др. [26], А. В. Девяткиным и др. [2] и — строгим образом — путем моделирования кривых блеска А. В. Мельниковым [46]. Недавно обработка наблюдений с КА «Кассини» позволила Харбисон и др. [40] сделать вывод о неоднородности распределения вещества внутри Гипериона и несовпадении геометрических осей его фигуры с осями инерции.

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

Цель диссертационной работы. В работе были поставлены и решены следующие задачи:

1. Численно-экспериментальное исследование резонансных и хаотических режимов вращательной динамики спутников планет и анализ наблюда­ тельных проявлений этих режимов.

2. Развитие методов и программных средств для исследования вращатель­ ной динамики спутников, основанное на массовом вычислении значений характеристических показателей Ляпунова путем численного интегри­ рования уравнений движения.

3. Построение диаграмм устойчивости вращательных режимов спутников планет, сравнение результатов численного моделирования с аналитиче­ ской теорией; выявление качественных закономерностей в хаотическом вращении с целью определения границ применимости теории.

Научная новизна. В процессе выполнения работы был получен ряд новых результатов:

1. Создан новый программный комплекс для численного интегрирования и вычисления ляпуновских спектров динамических систем с непрерыв­ ным временем, ориентированный на анализ вращательной динамики спутников планет.

2. Впервые численно-экспериментально подтверждены выводы теории се­ паратрисных отображений о свойствах хаотической вращательной ди­ намики спутников планет.

3. Впервые выявлены наиболее вероятные кандидатуры (помимо Гипери­ она) — спутники Сатурна Прометей и Пандора — для наблюдательного поиска проявлений хаоса во вращательной динамике спутников планет.

Научная и практическая значимость работы Созданные в рамках данной диссертационной работы методика и про­ граммный комплекс для расчета ляпуновских спектров динамических систем с непрерывным временем могут быть использованы как инструмент модели­ рования для выявления различных качественных закономерностей во вра­ щательной динамике спутников планет. Универсальность методов и их про­ граммной реализации позволяет распространить их использование также на более широкий круг задач динамики тел Солнечной системы — как враща­ тельной, так и орбитальной.

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

Следует отметить, что с использованием развитых в настоящей диссер­ тационной работе программных средств и методик был получен результат о режимах вращения Гипериона и Фебы, вошедший в перечень НСА РАН важнейших достижений астрономических исследований в России в 2008 г.

На защиту выносятся следующие основные результаты и поло­ жения:

1. Программная реализация алгоритмов расчета ляпуновских спектров ди­ намических систем с непрерывным временем. Создание программного комплекса, ориентированного на анализ вращательной динамики спут­ ников планет.

2. Численно-экспериментальное подтверждение выводов теории сепаратрис­ ных отображений о свойствах хаотической вращательной динамики спут­ ников планет.

3. Эмпирические зависимости компонент ляпуновского спектра от инер­ ционных параметров в задаче о пространственном вращении спутника.

4. Выявление наиболее вероятных кандидатур планетных спутников (по­ мимо Гипериона), которые могут находиться в хаотическом вращении, — а именно, 16-го и 17-го спутников Сатурна Прометея и Пандоры.

Апробация работы. Результаты работы докладывались на семинарах научных подразделений ГАО РАН и на следующих конференциях:

1. Всероссийская астрономическая конференция «ВАК–2001», С.-Петербург, АИ СПбГУ, 6–11 августа 2001 г.;

2. «Небесная механика — 2002. Результаты и перспективы», С.-Петербург, ИПА РАН, 10–14 сентября 2002 г.;

3. Всероссийская астрономическая конференция «ВАК–2004», Москва, ГАИШ МГУ, 2004 г.;

4. «Астрономия–2005 — современное состояние и перспективы», Москва, ГАИШ МГУ, 1–6 июня 2005 г.

Публикации. Материалы диссертации опубликованы в 7 печатных ра­ ботах, из них 4 статьи в журналах, входящих в перечень ВАК рецензируемых научных изданий, и 3 статьи в других изданиях.

Личный вклад автора. Содержание диссертации и положения, выно­ симые на защиту, отражают персональный вклад автора в опубликованные работы. В опубликованных по теме диссертации работах подготовка к публи­ кации полученных результатов проводилась совместно с соавторами, причем вклад диссертанта был равнозначным с соавторами.

Структура и объем диссертации. Диссертация состоит из введения, 4 глав, заключения, приложения и библиографии. Общий объем диссертации 148 страниц c Приложением, включая 26 рисунков и 8 таблиц. Библиография включает 71 наименование.

Возникновение небесной механики как науки непосредственно связано с созданием Ньютоном теории тяготения. С ее открытием один из важнейших разделов практической астрономии – расчет эфемерид (положений небесных тел в заданные моменты времени) – обрел твердую математическую базу.

Родившись из крайне насущных целей измерения времени и навигации, тео­ рия движения небесных тел явилась в то же время одним из наиболее ярких успехов в познании человеком законов строения Вселенной и, по мере роста точности предсказания положений Луны и планет и с открытием Адамсом и Леверье Нептуна, стала подлинным триумфом механистической картины мира. Это предопределило безграничную веру Лапласа в детерминизм, гос­ подствовавшую в науке в течение всего XIX столетия.

Но постепенно небесномеханические расчеты все усложнялись. Несмот­ ря на простоту исходных уравнений, отыскать их точное аналитическое ре­ шение в применении к задачам движения тел Солнечной системы удалось только в нескольких простейших случаях. Поэтому, начиная с работ Лагран­ жа и Лапласа в конце XVIII и начале XIX века, небесная механика пошла по пути использования теории возмущений. Однако ряды, даваемые аналитиче­ ской теорией возмущений, требуют большого объема алгебраических выкла­ док, которые приходилось проводить вручную. Так, опубликованная Делоне в 1867 году в результате кропотливых 20-летних расчетов теория движения Луны состоит их трех формул, каждая из которых занимает 200 страниц. Подчеркнем, что теория Делоне при этом являлась не просто неким узко­ Интересно отметить, что эта теория не имеет ошибок вплоть до 9-го порядка, за исключением одного члена 7-го порядка, несущественного для конечного результата. Это было показано в 1970 году Депри с помощью средств компьютерной алгебры.

специальным ограниченным научным результатом, а служила еще и весьма насущной практической цели вычисления морских навигационных таблиц.

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

Развитие классической небесной механики достигло апогея в работах Пу­ анкаре в конце XIX века, в особенности в его важнейшей работе по этой теме «Новые методы небесной механики». Но эта же работа выявила и границы аналитических средств. Пуанкаре показал, что большинство небесномехани­ ческих рядов расходится, так что с их помощью невозможно достичь сколь угодно точного решения. Также стало ясно, что в общем случае невозмож­ но найти аналитическое решение важнейшей задачи небесной механики — задачи тел. Вместе с тем, эта работа Пуанкаре содержала уже основные идеи современной теории динамических систем и, в частности, теории хаоса.

Эти идеи определили ход развития небесной механики, начиная с середины XX века.

Вплоть до этого времени в науке господствовали два взаимно дополни­ тельных взгляда на теорию динамических систем. С одной стороны, основан­ ная на ньютоновой теории, а затем и на релятивистской теории Эйнштейна небесная механика была средоточием регулярности и детерминизма. С дру­ гой, статистическая механика, основы которой были заложены в XIX веке в работах Больцмана и Гиббса, рассматривала статистические свойства всей со­ вокупности частиц в целом, игнорируя индивидуальные траектории частиц и полагая их случайными. Видимое противоречие между детерминистским характером уравнений, которым подчиняется динамика отдельных частиц, и их случайным поведением было преодолено в XX веке с созданием эргодиче­ ской теории в работах Биркгофа, Синая и других; эта теория дала строгое математическое обоснование статистической механике. В астрономии первый подход целиком господствовал в небесной механике, а второй получил распро­ странение в звездной динамике. Так или иначе, было широко распространено мнение, что динамика реальных систем является либо регулярной (в случае небесномеханических систем — систем с малым числом тел), либо эргоди­ ческой (в случае звезднодинамических систем — систем с большим числом тел).

Как было отмечено, точное решение большинства задач, связанных с взаимодействием трех и более тел, невозможно получить аналитическими методами. Численное же их решение вручную сопряжено с колоссальными вычислительными трудностями. Поэтому вполне естественно, что появление в середине XX века первых электронно-вычислительных машин сразу же при­ влекло внимание исследователей. Таким образом, история развития числен­ ных методов и применения их в научных исследованиях — и, в частности, в небесной механике и динамической астрономии — насчитывает всего полвека, и за это время компьютеры успели оказать огромное влияние на большинство областей науки.

Однако сама потребность в механизации процесса вычислений имеет, возможно, столь же древнюю историю, как и начало широкого использова­ ния математики в экономической и хозяйственной деятельности человека во­ обще. Древнейшее известное приспособление для счета — абак — достоверно упоминается с V–IV века до н. э. как «саламинская доска»; оно было известно также в Китае и Японии и в видоизмененном виде — как, например, создан­ ные в XVI веке русские счеты — дошло до наших дней.

В первой половине XIX века прогресс прецизионной механики привел к появлению «бытовых» механических вычислительных устройств – арифмо­ метров2, которые, наряду с логарифмической линейкой, широко применялись Патент на первый арифмометр получил в 1820 году К. Томас из Германии, который занялся промышленным производством этих устройств и изготовил за 50 лет 1500 экземпляров.

в научных исследованиях в течение, по крайней мере, трех четвертей XX ве­ ка. Однако эти устройства всего лишь облегчали выполнение ручных вы­ числений. В тот же период Чарлз Бэббедж сконструировал «аналитическую машину» совершенно нового типа, способную хранить данные и выполнять различные программы. Идеи Бэббеджа и его ученицы, первого программи­ ста Ады Лавлейс, оказали большое влияние на кибернетику3 XX века, но сама машина, ввиду ее экзотичности и трудоемкости изготовления, распро­ странения не получила. Качественный прорыв произошел лишь после созда­ ния электронно-вычислительных машин, основные элементы архитектуры ко­ торых — использование двоичной арифметики, процессор, работающий под управлением программы, наличие устройств хранения данных и устройств ввода-вывода — сохраняются и сейчас и, по мере развития электронных тех­ нологий и связанного с ним роста скорости вычислений, позволяют решать все более сложные задачи.

1.1. Численное моделирование в задаче тел Первая серьезная попытка использования цифровой электронно-вычис­ лительной машины в задаче тел была предпринята в 1953 году в Лос-Ала­ мосе Ферми, Паста и Уламом на одном из первых компьютеров MANIAC.

Рассмотренная этими исследователями система состояла из осцилляторов, расположенных вдоль прямой и моделирующих колеблющуюся струну; си­ ла взаимодействия между соседними осцилляторами полагалась линейной, с малой нелинейной добавкой. Если бы сила была в точности линейной, энер­ гия каждой из колебательных мод, заданных начальными условиями, сохра­ нялась бы. Наличие же нелинейного возмущения, как предполагал Ферми, Термин «кибернетика» был введен в 1834 году Ампером для гипотетической науки об управлении обществом и государством.

приведет к тому, что энергия со временем равномерно распределится меж­ ду всеми модами, то есть, в соответствии с предсказаниями статистической механики, система придет к тепловому равновесию. Вопреки этому, реальное поведение системы оказалось далеко от ожидаемого эргодического поведения и гораздо сложнее его — энергия каждой моды квазипериодически возвраща­ лась к своему первоначальному значению. Позднее было показано, что пара­ докс Ферми–Пасты–Улама можно объяснить наличием двух различных режи­ мов — квазипериодического, являющегося следствием наличия формального третьего интеграла движения системы, и хаотического, объясняемого пере­ крытием резонансов в фазовом пространстве. Этот результат явился полной неожиданностью для научного сообщества, воспитанного на представлении о том, что сложные системы подчиняются законам статистической механики.

Он показал, насколько нетривиальным может быть поведение нелинейных систем, и продемонстрировал важность численных экспериментов для их ис­ следования. С этим же результатом можно связать рождение нелинейной динамики как полноправного научного направления.

В 1956 году компьютер был использован П.-О. Линдбладом в Стокголь­ ме для выявления механизма образования спиральной структуры галактик.

Моделируя траектории движения звезд в плоской галактике, Линдблад пока­ зал, что спиральные рукава закручиваются в направлении против вращения галактики. В том же году к этой работе подключился Дж. Контопулос. Рас­ считав трехмерные траектории звезд в галактике, он обнаружил, что, вопре­ ки предсказаниям господствовавшей тогда эргодической теории, траектории не заполняют все пространство, а образуют ограниченные области, в про­ екции представляющие собой криволинейные параллелограммы, напоминаю­ щие деформированные фигуры Лиссажу. В 1960 году Контопулос доказал, что этот численный результат также можно объяснить посредством третьего интеграла движения, получившего впоследствии название «интеграла Конто­ пулоса».

Примерно в те же годы компьютерное моделирование применялось в за­ дачах звездной динамики, для исследования звездных скоплений и галактик как систем тел. Такие исследования были начаты фон Хорнером, исполь­ зовавшим модели всего с несколькими десятками тел. Позднее аналогичные вычисления были распространены на системы, состоящие из тысяч и мил­ лионов тел, и это позволило объяснить многие детали эволюции звездных систем.

Тем не менее, в эти годы использование компьютеров в научных прило­ жениях встречало и суровое противодействие, особенно в среде математиков.

Так, «после смерти Джона фон Неймана его бывшие коллеги из Института перспективных исследований в Принстоне на много лет очистили свои здания от компьютеров» [3, с. 127]. Однако общая тенденция была, все же, противопо­ ложной, и компьютеры в эти годы стали прочно входить в практику научных исследований.

В 1966 году Виктор Себехей предпринял попытку исследовать с помо­ щью компьютера «пифагорейскую» задачу трех тел, в которой тела с масса­ ми, относящимися как 3:4:5, помещены в вершинах пифагорейского треуголь­ ника со сторонами, находящимися в том же отношении, имеют первоначаль­ но нулевые скорости и движутся под действием взаимного ньютоновского притяжения. Эта задача была исследована численными методами и ранее, математиками Мейсселем в конце XIX в. и Бурро в первой четверти XX в., однако трудоемкость ручных вычислений не позволила им достичь успеха и сделать какие-либо качественные выводы о динамике такой системы. Себехей выполнил расчеты первоначально с помощью М. Стендиша на компьютере Йельского университета, а затем, в соавторстве со Спинелли и Лекаром, в Нью-Йорке, в Институте космических наук NASA. Одновременно аналогич­ ная работа была проведена Л. Станеком в Цюрихе. Полученный этими иссле­ дователями результат оказался крайне неожиданным: через некоторое время два тела из трех образуют связанную систему, а третье на огромной скорости выбрасывается из системы благодаря своего рода «эффекту рогатки». Сейчас этот результат широко известен, подтвержден на многочисленных примерах и доказан аналитически. Он проливает свет на образование двойных звезд­ ных систем и позволяет указать источник происхождения «звезд-странни­ ков», движущихся в Галактике с огромными скоростями и даже покидающих ее пределы.

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

Характерным примером служат работы Тоомре 70-х годов, в которых моде­ лируется поведение взаимодействующих галактик как больших ансамблей ча­ стиц. И, наоборот, специалисты в области звездной динамики осознали необ­ ходимость обратиться к методам небесной механики и рассмотрению индиви­ дуальных траекторий звезд для разрешения таких парадоксов, как вычислен­ ное Чандрасекаром время релаксации для звездной системы, превышающее возраст Вселенной — так что звезды не могли бы к настоящему времени до­ стичь наблюдаемого состояния статистического равновесия.

В применении к динамике Солнечной системы один из наиболее мас­ штабных численных экспериментов был выполнен в 1988 году Зюссманом и Уиздомом, которые использовали специально сконструированный для этого компьютер, названный «Цифровым планетарием», для вычисления орбит пя­ ти внешних планет на интервале 1/5 возраста Солнечной системы. В ходе этого эксперимента было, в частности, показано, что орбитальное движение Плутона является хаотическим. Примерно тогда же был осуществлен проект LONGSTOP (“LONg-term Gravitational Stability Test of the Outer Planets” — долговременный тест гравитационной устойчивости внешних планет), заклю­ чавшийся в численном интегрировании движения пяти планет на интервале в 100 млн лет. Кроме ответа на вопрос об устойчивости Солнечной систе­ мы и о ее будущем, такие эксперименты помогают, в частности, уточнить значения частот, амплитуд и фаз для возмущений, что позволяет вывести бо­ лее точную вековую теорию движения планет. В дополнение к эксперименту LONGSTOP на компьютере Cray-1S Лондонского университета Нобили и ее коллеги провели моделирование движения известных тогда пяти спутников Урана, позволившее, в частности, уточнить значения масс этих спутников.

В свою очередь, потребности численного эксперимента привели к созда­ нию новых эффективных вычислительных методов — таких, как симплекти­ ческие интеграторы, широко используемые сейчас в задачах моделирования динамики гамильтоновых систем.

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

1.2. Численный эксперимент и динамический хаос Как отмечалось выше, до середины XX века считалось, что динамиче­ ские системы являются либо регулярными, либо эргодическими. Развитие теории динамического хаоса, основы которой были заложены еще в работах Пуанкаре, показало, что практически любая нелинейная система может де­ монстрировать хаотическое поведение, определяющим свойством которого яв­ ляется непредсказуемость движения, то есть существенная зависимость его от малых изменений начальных условий, несмотря на детерминистский ха­ рактер уравнений движения. В частности, Хенон и Хейлес в 1964 году в чис­ ленном эксперименте впервые продемонстрировали хаотическое поведение простой неинтегрируемой гамильтоновой системы, названной впоследствии в их честь системой Хенона–Хейлеса. При этом, однако, сохраняются остро­ ва устойчивого движения. Таким образом, полностью упорядоченные — так же, как и полностью хаотические — системы являются достаточно исклю­ чительным случаем в природе; абсолютное большинство реальных систем в динамической астрономии может характеризоваться обоими указанными ти­ пами движения.

Важность роли хаоса в динамике Солнечной системы была осознана в 80-е годы XX века. С этим связано понимание таких явлений, как наличие люков Кирквуда в поясе астероидов, в которых практически нет вещества, и нерегулярное пространственное вращение 7-го спутника Сатурна, Гипериона.

Широко известное как «эффект бабочки» — непредсказуемость отдаленных последствий даже самых незначительных воздействий — явление хаоса имеет огромное значение и для эволюции всей Солнечной системы [20]. Принципи­ ально непредсказуемый характер некоторых явлений заставляет, ни в коей мере не умаляя прогностической ценности научной теории, внести корректи­ вы в интерпретацию ее результатов и связь их с наблюдательными данны­ ми — так же, как это произошло с квантовой теорией в 20-х годах прошлого века.

Наряду с аналитическими методами, обеспечивавшими «фундамент» и строгое обоснование теории хаоса, важнейшую роль в ее развитии играли чис­ ленные методы. Дополняя данные наблюдений, численно-экспериментальные результаты позволяют проверить справедливость качественных оценок раз­ личных параметров хаотического движения — таких, как, например, характе­ ристический показатель Ляпунова, величина, обратная которому, дает время предсказуемости движения — и определить границы применимости этих оце­ нок. Например, критическая величина возмущения, при которой происходит скачкообразный переход к крупномасштабному хаосу, определяется, как пра­ вило, численно. Кроме того, численное моделирование, ставшее с появлением компьютеров полноправным инструментом исследования, позволяет нагляд­ но проиллюстрировать проявления хаотического поведения и поставляет бо­ гатый материал для выявления новых динамических закономерностей.

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

В качестве другой иллюстрации может служить численное моделиро­ вание возможной динамической эволюции орбиты астероида Хирон, выпол­ ненное Оикавой и Эверхардтом в 1979 году. Орбита этого астероида сильно вытянута, и перигелий ее лежит внутри орбиты Сатурна, а афелий — вбли­ зи орбиты Урана. Численные эксперименты показали, что в будущем Хирон испытает несколько тесных сближений с планетами; при этом незначитель­ ная разница в начальных условиях, в пределах той точности, с которой была известна орбита Хирона, приводит к совершенно различным сценариям его дальнейшей судьбы после сближения с Сатурном — он может как перейти во внутреннюю часть Солнечной системы, так и полностью покинуть Сол­ нечную систему. Это одно из наиболее наглядных проявлений динамического хаоса. Аналогичная картина имеет место в случае кометы Шумейкер–Леви 9, упавшей на Юпитер в июле 1994 года. Численное моделирование показыва­ ет, что первоначально эта комета имела, по-видимому, орбиту с малым экс­ центриситетом, лежащую внутри орбиты Юпитера, но перешла на орбиту, пересекающуюся с орбитой Юпитера, примерно в первой половине XX века.

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

Таким образом, развитие компьютерных методов оказало огромное вли­ яние на понимание важности хаотического поведения в динамике тел Сол­ нечной системы и ее эволюции.

1.3. Методы компьютерной алгебры в небесной механике и динамической астрономии Как было уже отмечено, аналитические выкладки в небесной механи­ ке, связанные с использованием разложений возмущающей функции, крайне длинны и трудоемки. Однако сами по себе они достаточно рутинны, сводятся к набору хорошо формализуемых правил и легко поддаются алгоритмизации.

По этой причине представляется вполне естественным поручить эту задачу компьютерным программам. Первые такие программы, обеспечивающие ав­ томатизацию символьных вычислений, — системы компьютерной алгебры — начали создаваться с 60-х годов XX века. Наиболее известны среди них паке­ ты MAO, TRIGMAN, CAMAL и некоторые другие. Они предназначались, как правило, для узкоспециализированных задач. Наиболее заметным успехом таких систем явилась отмеченная выше проверка теории Луны Делоне. Пер­ вые универсальные системы компьютерной алгебры, предназначенные для работы на персональных компьютерах, — такие, как REDUCE, MACSYMA, Maple, DERIVE, — стали появляться с конца 70-х годов. Однако первона­ чально их возможностей было недостаточно для серьезного применения в небесной механике.

Во второй половине 80-х годов XX века Ласкар предложил комбини­ рованный подход, сочетающий численное интегрирование с аналитическими разложениями и позволяющий эффективно моделировать движение планет на интервалах времени, сравнимых с возрастом Солнечной системы. В этом подходе алгоритмы компьютерной алгебры используются для вывода усред­ ненных уравнений (усредняются короткопериодические эффекты, не влияю­ щие на долговременную эволюцию),4 а интегрирование этих уравнений про­ Ласкар разработал для вывода усредненных уравнений собственные специализированные проце­ дуры на языке FORTRAN.

водится численно с гораздо бльшим, чем при прямом численном интегри­ ровании исходных уравнений движения, шагом по времени. Применив такой подход, Ласкар в 1988 году показал, что движение внутренних планет явля­ ется хаотическим, а характерное время экспоненциальной расходимости тра­ екторий (то есть характерное время предсказуемости движения) составляет порядка 5 млн лет. Близкие результаты были получены и в других иссле­ дованиях, причем было показано, что орбиты всех четырех планет-гигантов (Юпитера, Сатурна, Урана и Нептуна) также являются хаотическими.

В настоящее время рост мощности персональных компьютеров и разви­ тие универсальных систем компьютерной алгебры способствуют применению таких систем в широком классе небесномеханических задач. Они используют­ ся, в частности, для расчета эфемерид тел Солнечной системы с учетом возму­ щений от планет и крупных астероидов и релятивистских эффектов. Напри­ мер, точность определения положения Луны достигает при этом нескольких сантиметров. Использование компьютерных средств для символьных вычис­ лений помогает также, например, получить аналитические выражения для разложения в ряд планетной возмущающей функции до высоких порядков, вычислить нормальные формы различных гамильтоновых систем, решить аналитически обобщенное уравнение Кеплера и некоторые другие уравнения классической небесной механики, разработать эффективные теории движе­ ния искусственных спутников Земли и многое другое.

1.4. Выводы к первой главе Использование компьютерных методов в современной динамике Солнеч­ ной системы связано с двумя основными группами задач. В первой из них численное интегрирование уравнений движения реальных или модельных си­ стем используется для выявления качественных закономерностей регулярной и хаотической динамики различных групп тел, составляющих Солнечную си­ стему — планет, их спутников, астероидов, комет, транснептуновых объектов.

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

Поскольку пока не существует аналитического решения проблемы устой­ чивости Солнечной системы, большой интерес представляют исследования этого вопроса численными методами. Численное интегрирование уравнений, моделирующих Солнечную систему, на больших временных масштабах поз­ волило выявить хаотический характер орбитального движения планет. Од­ нако, несмотря на это, нет никаких указаний на то, что Солнечная система неустойчива и распадется в будущем: результаты всех проведенных до сих пор численных экспериментов говорят о том, что все планеты (кроме, воз­ можно, Меркурия) сохраняют орбиты, близкие к существующим сейчас, на временах в миллиарды лет. Наличие же хаоса означает ограниченность того срока, на который мы можем предсказать положение планет.

Помимо «классических» приложений компьютерного моделирования си­ стем частиц в динамике Солнечной системы и в звездной динамике, сейчас оно широко применяется и в «новых» областях астрономии — в таких, напри­ мер, как динамика экзопланетных систем, динамика аккреционных дисков звезд и черных дыр, эволюция крупномасштабной структуры Вселенной.

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

Хаотическое вращение спутников планет:

ляпуновские спектры и максимальные характеристические показатели Ляпунова 2.1. Введение Вычисление характеристических показателей Ляпунова (ХПЛ) служит одним из важнейших инструментов качественного анализа движения, в част­ ности, в небесной механике. ХПЛ характеризуют скорость экспоненциальной расходимости траекторий, близких друг к другу в фазовом пространстве.

Величина, обратная максимальному ХПЛ (МХПЛ), дает время предсказуе­ мости движения. ХПЛ тесно связаны с динамической энтропией [12, 16, 25, 30, 45].

В настоящей главе вычисляются ХПЛ хаотического вращения малых спутников планет. Полагается, что динамически асимметричный спутник дви­ жется по фиксированной эллиптической орбите. Используются имеющиеся данные по инерционным и орбитальным параметрам для выборки спутников Марса, Юпитера, Сатурна и Нептуна. В настоящее время большинство из них, вероятнее всего, вращается регулярно в синхронном спин-орбитальном резонансе. Однако, несмотря на это, начальные значения динамических пере­ менных всегда выбираются в хаотической области фазового пространства для того, чтобы исследовать хаотический режим вращения, в котором каждый из спутников мог находиться до захвата в синхронный резонанс. Данное исследо­ вание охватывает и случай плоского (ось вращения ортогональна плоскости орбиты), и случай пространственного вращения. Полные ляпуновские спек­ тры вычисляются при помощи HQRB-метода фон Бремена и др. [66]. Для расчета МХПЛ используется также традиционный метод «теневой траекто­ рии» [7].

Для лучшего качественного понимания результатов численного модели­ рования описывается также методика аналитического оценивания МХПЛ в рамках модели нелинейного резонанса как возмущенного нелинейного маят­ ника, разработанная И. И. Шевченко [18, 19] в теории сепаратрисных отоб­ ражений, для получения аналитических оценок МХПЛ в задаче плоского вращения. Кроме того, исследуется вопрос о применимости данных оценок к случаю пространственного вращения.

2.2. Основные определения и алгоритмы ХПЛ траектории имеет физический смысл средней скорости расходимо­ сти траекторий, близких к данной. Ненулевое значение ХПЛ указывает на хаотический характер движения, в то время как нулевое служит проявлени­ ем регулярного (периодического или квазипериодического) движения.

Рассмотрим две траектории, близкие друг к другу в фазовом простран­ стве. Одну из них назовем ведущей, а другую — теневой. Пусть (0 ) — длина вектора смещения, направленного от ведущей к теневой траектории, в началь­ ный момент 0. ХПЛ определяется следующей формулой [7]:

Для гамильтоновых систем величина может принимать 2 различных зна­ чений (в зависимости от направления первоначального смещения), где — число степеней свободы; ХПЛ разделяются на пары: для каждого > В практических вычислениях всегда получается лишь приближение к значению ХПЛ за конечный промежуток времени. Например, так называе­ мый метод теневой траектории дает оценку ХПЛ по формуле [7] где = /1, обозначает расстояние между ведущей и теневой фазовы­ ми точками на -м шаге итерации, — длина шага, а — число шагов, на которых измеряется расстояние. При вычислении МХПЛ согласно (2.1) необходимо проводить периодическую перенормировку положения теневой фазовой точки по отношению к ведущей так, чтобы расстояние оставалось бы малым.

Заметим, что в практических вычислениях с использованием форму­ лы (2.1) получается оценка только лишь МХПЛ, поскольку начальные дан­ ные теневых траекторий, которые могут породить остальные элементы ляпу­ новского спектра, принадлежат множеству меры нуль [7]. Поэтому начальные данные для теневой траектории можно выбирать произвольным образом — однако так, чтобы расстояние от ведущей траектории было достаточно ма­ лым.

Рассмотрим алгоритм вычисления полного ляпуновского спектра для случая динамической системы с дискретным временем. Для того, чтобы вы­ числить спектр, необходимо знать матрицу касательного отображения [66].

Пусть есть исходное отображение. Тогда матрица касательного отображения опре­ деляется соотношением (x обозначает положение фазовой точки на -м шаге итерации, а x — со­ ответствующий касательный вектор). Обозначим матричное представление касательного отображения () на -м шаге в стандартном базисе как ().

Тогда приближение для каждого ХПЛ дается выражением где =, = 1,..., 2, а верхнетреугольная матрица () получается итерациями посредством QR-факторизации произведения () (1) : то есть где — ортогональная, а (0) — единичная матрица.

Истинное значение ХПЛ является пределом (или в случае, ес­ ли вычисляется полный ляпуновский спектр) при. Здесь и далее значения (или ) при больших принимаются в качестве истинных значений (соответственно, ) при численном оценивании ХПЛ. Следует, однако, отметить, что формально эти величины различны. Традиционный численный метод нахождения предельных значений ХПЛ (см. напр. [68]) за­ ключается в следующем: строя зависимость log от log, мы ищем зна­ чение log, при котором эта зависимость «достигает насыщения», то есть выходит на горизонтальное плато. Значение на этом плато принимается в качестве истинного значения.

2.3. Аналитическое оценивание МХПЛ В работах И. И. Шевченко [18, 19] в рамках теории сепаратрисных отоб­ ражений был предложен метод оценивания МХПЛ движения в хаотической области в окрестности сепаратрис возмущенного нелинейного резонанса. Ре­ зонанс моделируется следующим гамильтонианом:

где = +0. Переменная — резонансная фаза, — сопряженный импульс;

— угол фазы возмущения, а 0 — его начальное значение. Величина обозначает постоянную частоту возмущения;,,, — константы.

Сепаратрисное отображение в форме Б. В. Чирикова [7, 15, 16, 30] описы­ вает движение в окрестности сепаратрис нелинейного резонанса под действи­ ем симметричного периодического возмущения. Различные аспекты общей теории сепаратрисных отображений были рассмотрены в работах [1, 14, 17, 21–23, 54, 55, 65] и многих других. Эта теория имеет ключевое значение для изучения хаотической динамики небесных тел. Построение и анализ сепара­ трисных отображений является мощным инструментом современной нелиней­ ной динамики — прежде всего, для изучения хаотического поведения.

В рамках модели возмущенного нелинейного маятника И. И. Шевчен­ ко [17, 52] была разработана процедура приведения сепаратрисного отображе­ ния к единой поверхности сечения фазового пространства исходной гамильто­ новой системы (процедура синхронизации отображения). При помощи этой процедуры было показано [17], что сепаратрисные отображения адекватно описывают фазовые портреты движения около сепаратрис как в случае вы­ соких, так и в случае низких частот возмущения. В применении к задачам небесной механики были выведены общие сепаратрисные отображения, под­ ходящие к случаям асимметричных возмущений и возмущений, характер­ ных для задач с орбитальными резонансами [18, 19, 54, 55]. В частности, И. И. Шевченко [18, 19, 54] применил теорию сепаратрисных отображений к построению фазовых портретов хаотического движения, к анализу резо­ нансной структуры фазового пространства и к аналитическому оцениванию МХПЛ в задаче плоской вращательной динамики асимметричного спутника на эллиптической орбите.

В случае симметричного возмущения ( = ) сепаратрисное отображе­ ние имеет два параметра. Первый,, есть отношение, частоты возмущения, к 0 = ()1/2, частоте малых колебаний на резонансе. Второй дается фор­ мулой где есть интеграл Мельникова–Арнольда, как определено в [15, 17, 30, 52].

Для асимметричного возмущения ( = ) имеем два параметра возму­ щения вместо единственного. Две эти величины, + и, являются значениями для прямого и обратного движения модельного маятника, со­ ответственно:

Следуя И. И. Шевченко [18, 19], примем зависимость МХПЛ симметрич­ ного сепаратрисного отображения от в виде где 0.8 — постоянная Чирикова.

Если =, то средний период хаотического вращения (или, что то же самое для движения, близкого к сепаратрисе, средний полупериод хаотиче­ ской либрации) модельного маятника различается для прямого и обратного движений:

МХПЛ в компонентах хаотического слоя, отвечающих прямому и обратному вращениям модельного маятника, на единицу времени исходной гамильтоно­ вой системы равен [18, 19]. Вычисление среднего (среднего по всему слою) представляет со­ бой весьма сложную задачу. В частности, необходимо знать относительные средние времена нахождения системы в трех различных компонентах хаоти­ ческого слоя, отвечающих прямому вращению, обратному вращению и либра­ ции маятника. Эти времена зависят от параметров системы — среди других, от асимметрии возмущения. Поэтому здесь принимается просто в качестве оценки МХПЛ, усредненной по всему хаотическому слою. Такое усреднение является лишь грубым приближением, но вывод более точной формулы потребовал бы гораздо более сложной и строгой теории, основыва­ ющейся на многих численных экспериментах с сепаратрисными отображени­ ями.

исходной гамильтоновой системы. Для задач данной главы этот период есть ни что иное, как период обращения спутника orb. Тогда МХПЛ в расчете на единицу времени (скажем, на сутки, если orb выражен в сутках) равен [18, 19]. Имеющая размерность времени величина 1 — это так назы­ ваемое ляпуновское время. Она характеризует время предсказуемости хаоти­ ческого движения.

Таким образом, если параметры сепаратрисного отображения известны, формулы (2.10) – (2.14) позволяют получить теоретическую оценку МХПЛ.

2.4. Численные методы определения полных ляпуновских спектров и МХПЛ Для вычисления ХПЛ в данной главе используется два метода. Тра­ диционный метод «теневой траектории» (иначе «двухчастичный» метод) непосредственно следует из определения (2.1); он эффективен и прост в реализации. Однако, как было отмечено выше, он позволяет получить только МХПЛ, но не весь спектр.

К тому же, метод теневой траектории может давать ошибочные оценки ХПЛ, если начальная теневая траектория выбрана неправильно (см. напр.

[61]). А именно, оценки ХПЛ, полученные этим методом, могут зависеть от величины начального сдвига теневой траектории 0. Если этот сдвиг слиш­ ком велик, векторы смещения в фазовом пространстве могут перестать быть аппроксимацией касательных векторов. С другой стороны, нельзя сделать этот сдвиг произвольно малым из-за определяющего влияния ошибок округ­ ления при приближении к машинной точности представления чисел. В обоих случаях ошибки накапливаются в точках ренормировки, что может привести к ошибочным оценкам ХПЛ. Танкреди и др. [61] рекомендуют проверять, нет ли зависимости оценок ХПЛ от 0. Для вычислений, приведенных ниже, был выполнен такой тест, и было найдено, что значение сдвига 0 = удовлетворяет целям данной главы. Как отмечено в [61], необходимо также правильно выбрать интервал времени ренормировки. Слишком малое его значение может привести к быстрому накоплению ошибок округления при ренормировке, а слишком большое — к арифметическому переполнению. В описанных ниже экспериментах использовались различные значения шага ре­ нормировки; показано, что во всех случаях подходит величина, равная шагу итерации.

Значительно более эффективный метод вычисления ХПЛ предложен фон Бременом и др. [66]. Он основан на QR-разложении матрицы касательно­ го отображения (2.3) с использованием преобразований Хаусхолдера и изве­ стен поэтому как метод HQRB (Householder QR-based). Суть метода состоит в вычислении рефлекторов Хаусхолдера — набора из 1 ( — размерность динамической системы) матриц () ( = 1, 2,..., 1), с помощью которых можно явным образом осуществить факторизацию (2.5):

Матрицы-рефлекторы имеют следующую структуру:

где — единичная матрица, и первые ( 1) элементов -векторов () равны нулю. В выражение (2.4) для расчета ХПЛ входят только диагональ­ ные элементы матриц (), а матрицы () используются там только в виде произведения () (1). Это позволило фон Бремену и др. [66] оптимизи­ ровать алгоритм расчета ХПЛ на основе вычисления рефлекторов Хаусхол­ дера, а также показать, что метод HQRB обладает более высокой числен­ ной устойчивостью, чем метод QR-разложения, основанный на ортогонали­ зация Грама–Шмидта (последний использовался, в частности, Бенеттином и др. [25]). Метод HQRB также менее чувствителен к величине шага итера­ ции, чем другие алгоритмы факторизации, поскольку не использует ренор­ мировки. Далее все вычисления полных ляпуновских спектров выполнены с помощью программного комплекса HQRB, описанного в Приложении и осно­ ванного на методе фон Бремена и др. [66].

В случае, когда матрица касательного отображения () не задана явным образом, существует три возможности для ее вычисления. Во-первых, мож­ но заменить касательные векторы x в (2.3) малыми векторами смещения:

x = x x, где x и x — ведущая и теневая фазовые точки, соответствен­ но, на шаге. Тогда обе точки интегрируются независимо, согласно (2.2). Это приводит к x+1 = x+1 x+1. Повторив эту процедуру для 2 линейно независимых векторов x, можно решить уравнение относительно матрицы касательного отображения (). Этот метод требует знания только исходного отображения (2.2). Зависимость результата от на­ чального сдвига теневой траектории здесь не так критична, как в методе Бенеттина и др. [25]. Причина в том, что нет накопления ошибок округления при ренормировке. Этот алгоритм используется практически во всех вычис­ лениях в данной главе.

На производительности данного метода, однако, негативно сказывается необходимость выполнять 2 дополнительных итераций исходного отобра­ жения на каждом шаге. Тем не менее, существует еще один путь вычисления матрицы касательного отображения. Он применим, когда отображение (2.2) порождается непосредственно исходной динамической системой Тогда матрица касательного отображения может быть приближенно выраже­ на как где x ( ) — матрица Якоби системы при = =, а — единичная матрица. Это приближение, однако, выполняется только для достаточно ма­ лых величин шага итерации. По этой причине данный метод используется в этой главе только с целью тестирования.

Матрицу касательного отображения можно также вычислить непосред­ ственно с более высокой точностью, хотя и ценой несколько бльших вычис­ лительных затрат, путем одновременного интегрирования исходной и линеа­ ризованной систем. Такой подход не требует введения дополнительного ма­ лого параметра — смещения 0. Сопоставление его с первым методом будет проведено в следующей главе.

Как один из необходимых этапов исследования точности и надежности численного метода, в ходе подготовки к вычислениям были выполнены те­ сты с целью определить оптимальное значение. Использовались значения = 0.01 · 2, = 0.1 · 2 и = 1. Было показано, что величина ша­ га практически не отражается на результате вычислений; при этом меньшие значения оказываются чуть более предпочтительными.

Еще одним тестом надежности вычисления полного ляпуновского спек­ тра служит контроль суммы всех ХПЛ 2 (). В данной главе = 3, и указанную сумму можно для краткости обозначить как (6). Эта сумма должна равняться нулю (см. раздел 2.2), что позволяет контролировать внут­ реннюю точность метода. На фактическом значении (6) может сказаться как выбор начального сдвига теневой точки или величина шага итерации, так и точность используемого интегратора.

2.5. ХПЛ хаотического вращения спутников планет Значительная часть естественных спутников диаметром менее 500 кило­ метров имеет явно асимметричную форму (см. напр. таблицы в Ephmrides Astronomiques [32]). В настоящее время львиная доля информации о форме спутников получается из снимков с различных межпланетных космических аппаратов.

Многие известные спутники планет вращаются в синхронном резонан­ се с их орбитальным движением. Плоское вращение (в плоскости орбиты) в синхронном резонансе, как следует из теории приливных спин-орбитальных взаимодействий (см. напр. [24, 69]), является естественной конечной стадией долговременной динамической эволюции спутника. На этой стадии ось вра­ щения совпадает с осью наибольшего момента инерции.

Однако когда фазовая траектория в ходе долговременной динамической эволюции сближается с сепаратрисами синхронного резонанса и входит в ха­ отический слой, образованный сепаратрисами, спутник начинает хаотически кувыркаться, поскольку плоское вращение в этой области фазового простран­ ства неустойчиво по отношению к наклону оси вращения [9, 10, 69, 70]. Поэто­ му спутник вращается хаотически по всем углам Эйлера. Чем больше дина­ мическая асимметрия и эксцентриситет орбиты спутника, тем больше размер области хаотического движения в фазовом пространстве; правда, роли этих двух параметров совершенно различны.

Как отметил Уиздом [69], для того, чтобы спутник был захвачен в син­ хронный резонанс, прохождение его через область неустойчивости должно быть достаточно быстрым; соответственно, ширина хаотического слоя долж­ на быть достаточно малой. Конечно, движение в центре самог синхронного резонанса должно быть устойчиво по отношению к наклону оси вращения.

В ходе динамической эволюции спутник проходит через различные резонанс­ ные спин-орбитальные состояния и хаотические слои вблизи их сепаратрис.

Оценки МХПЛ в этих слоях могут дать важную информацию о характере динамической эволюции.

До сих пор только для Гипериона, седьмого спутника Сатурна, в на­ блюдениях был установлен хаотический характер вращения [2, 26, 42, 46].

Большинство других спутников, вероятно, находится в регулярном состоя­ нии, близком к синхронному [49]. Однако, в любом случае, захват в синхрон­ ный резонанс не может произойти без прохождения через главный хаотиче­ ский слой, поэтому знание величин МХПЛ в этом слое может иметь большое значение для изучения динамической истории спутников. В зависимости от сценария, можно перевычислить аналитические оценки МХПЛ для значений эксцентриситета, соответствующих различным этапам динамической эволю­ ции.

Для спутников с хаотическим или неизвестным характером вращения теоретические оценки ХПЛ дают информацию, полезную для планирования наблюдений этих спутников. Основным методом для наземных наблюдений является построение и анализ кривых блеска спутников [2, 42]. Величина, обратная МХПЛ (ляпуновское время), дает время предсказуемости хаотиче­ ского вращения. Ясно, что если целью наблюдений является изучение враща­ тельной динамики (возможно, хаотической), характерный интервал между наблюдениями должен быть меньше ляпуновского времени. Это накладыва­ ет важные ограничения на стратегию наблюдений.

В данном разделе рассмотрен хаотический режим вращения 11 избран­ ных спутников. В настоящее время из этих спутников только Гиперион изве­ стен как вращающийся хаотически (с очень большой вероятностью; см. [2, 26, 42, 46]). Для пяти спутников (Елена, Атлас, Прометей, Пандора и Протей) характер вращения неизвестен. Вращение Фобоса, Деймоса, Амальтеи, Эпи­ метея и Януса регулярно, и все они находятся в синхронном резонансе (см.

напр. [49, 69]). Следует отметить, что в рассмотрение включены спутники с из­ вестным регулярным характером вращения, но набор начальных данных для интегрирования для них, тем не менее, выбирается в хаотической области фа­ зового пространства (а именно, вблизи сепаратрис синхронного резонанса).

Иными словами, изучается возможный режим вращения спутника, который осуществлялся бы, если бы спутник не был захвачен в синхронный резонанс.

Полный набор спутников рассматривается для того, чтобы получить по воз­ можности более статистически значимый результат для сравнения с теорией и сопоставления различных динамических случаев. Кроме того, спутники, находящиеся в настоящее время в состоянии регулярного вращения, должны были оказаться в хаотическом состоянии на некоторой стадии их динами­ ческой эволюции, до их захвата в синхронный резонанс. Прохождение че­ рез хаотическую область вблизи сепаратрис синхронного резонанса в ходе приливной эволюции неизбежно. Эксцентриситеты орбит спутников, правда, могли во время этого прохождения иметь другие значения — но, поскольку невозможно судить об их значениях в то время, принимаются современные величины.

Истинные значения ХПЛ хаотического вращения могут отличаться от того, что мы ожидаем, еще по одной причине. Спутники предполагаются дви­ жущимися по невозмущенным эллиптическим орбитам. Хорошо известно, од­ нако, что в реальности элементы их орбит испытывают возмущения, вызван­ ные сплюснутостью центральной планеты и/или наличием других массивных спутников. Например, в случае Гипериона эксцентриситет орбиты колеблет­ ся от 0.08 до 0.12 с периодом в 18.8 лет [26] благодаря возмущениям от Титана. В приводимых ниже расчетах ляпуновские времена в большинстве случаев составляют порядка нескольких дней. Это гораздо меньше харак­ терных периодов возмущений. Следовательно, принимать во внимание такие возмущения представляется излишним.

2.5.1. Случай плоского вращения Рассмотрим движение асимметричного твердого тела («спутника») во­ круг его центра масс, который обращается по невозмущенной эллиптической орбите вокруг «планеты» (неподвижной массивной точки). Предполагается, что размеры спутника много меньше радиуса его орбиты, а масса пренебрежи­ мо мала по сравнению с массой планеты. Главные моменты инерции спутника по отношению к его осям инерции,1, обозначим как,,, соответ­ Обозначение обычно используется также в качестве большой полуоси орбиты спутника. Однако в модельных расчетах, проводимых в настоящей диссертационной работе, принята система единиц, в ственно. Оси выбираются таким образом, что.

В данном разделе сравниваются численные и теоретические оценки МХПЛ в задаче плоского хаотического вращения спутника, движущегося по фиксированной эллиптической орбите. Ось вращения спутника полагается совпадающей с осью его наибольшего момента инерции; она ортогональна плоскости орбиты. Вращательное движение описывается уравнением Белец­ кого [4] для плоских либраций и вращений спутника на эллиптической орбите.

В случае малых эксцентриситетов уравнения движения даются гамильтони­ аном (2.6), получаемым в пренебрежении всеми степенями эксцентриситета выше первой (см. напр. [29, 70]). Согласно И. И. Шевченко [18, 19], для пара­ инерционный параметр где — главные центральные моменты инерции спутника; — экс­ обращения.

Отображение (2.2) порождается уравнениями движения. Используется интегратор Хайрера и др. [39]. Это явный метод Рунге–Кутты 8-го порядка в реализации Дормана и Принса (см. [39]) с управляемой длиной шага. Инте­ грирование проводилось на временнм интервале [0, 105 ]. Такого времени оказалось достаточно для достижения насыщения (выхода на плато) зави­ симости от. Начальные условия для интегрирования были выбраны в хаотической области фазового пространства.

На рис. 2.1 показаны численные оценки МХПЛ для = 0.1 — среднего эксцентриситета орбиты Гипериона (см. напр. [70]). Параметр динамической которой значение большой полуоси орбиты равно 1, а то, какая из этих двух величин имеется в виду, всегда ясно из контекста.

асимметрии 0 меняется в диапазоне 0.2 0 1. Для сравнения с данны­ ми численного моделирования приводится соответствующая теоретическая кривая согласно (2.10–2.13).

Сравним теперь численные оценки МХПЛ с теоретическими для 11 ука­ занных выше спутников. Сравнение приведено в табл. 2.1. Значения инерци­ онных параметров /, / (или геометрические размеры, и ; тогда /, / вычисляются согласно модели эллипсоида однородной плотности) и эксцентриситеты орбит взяты из работ [59] (Фобос), [69] (Деймос, Амальтея, Янус и Эпиметей), [64] (Гиперион), [60] (Елена), [27] (Атлас), [36] (Прометей и Пандора), [32] (Протей). Значения /, /, необходимые для моделиро­ вания пространственного вращения и представленные ниже в табл. 2.3, дают по формуле (2.18) значение параметра 0, приведенного в табл. 2.1.

В графическом виде сопоставление результатов численного моделиро­ вания и теоретического оценивания МХПЛ для данной выборки спутников представлено на рис. 2.2. Теоретические оценки МХПЛ вычисляются на ос­ нове теории сепаратрисных отображений И. И. Шевченко [18, 19] (см. раз­ дел 2.3). Численные оценки являются средними от текущих значений ХПЛ по интервалу [5 · 104, 105 ], соответствующему плато насыщения. Для рас­ четов, результаты которых представлены в табл. 2.1, использованы два значе­ ния. В обоих случаях наблюдается практически одинаково хорошее согла­ сие с теорией. Для графического представления использованы результаты, полученные с = 0.01 · 2.

Рис. 2.1 и 2.2 оба наглядно показывают, что численные оценки находятся в хорошем согласии с теорией в случае плоского вращения.

В табл. 2.2 представлены результаты вычисления полных ляпуновских спектров плоского вращения методом HQRB. Следует отметить, что «вирту­ альные» теневые траектории в методе HQRB не ограничены плоскостью орби­ ты. Следовательно, (2) и (3) не обязательно равны нулю. Вместо уравнений Таблица 2.1. Случай плоского вращения. Инерционный параметр, эксцентриситет орбиты и численные оценки МХПЛ методом теневой траектории Рис. 2.1. Сравнение численных (кружки) и теоретических (кривая) оценок МХПЛ для случая эксцентриситета орбиты Гипериона движения, определяемых гамильтонианом (2.6), здесь использованы уравне­ ния Эйлера (см. раздел 2.5.2). Как определено выше, длина шага итерации = 0.01·2, а начальное смещение теневой частицы 0 = 107. Как и ранее, ХПЛ усредняются по интервалу времени [5 · 104, 105 ]. Здесь представлены средние от положительных ХПЛ и модулей соответствующих им отрицатель­ ных, то есть под () понимается 1 (() + |(+ ) |). Сумма (6), вычисляв­ шаяся для контроля внутренней согласованности метода HQRB, оказалась равной нулю для всех спутников с точностью до 6-го знака после запятой.

Для того, чтобы наглядно показать достижение насыщения при вычис­ лении ХПЛ, на рис. 2.3 приводится зависимость мгновенных оценок ХПЛ от времени для случая Гипериона.

Как можно видеть, сравнив табл. 2.1 и 2.2, значения МХПЛ, полученные методом теневой траектории, и значения (1), вычисленные методом HQRB, согласуются.

Таблица 2.2. Случай плоского вращения. Численные оценки компонент ляпуновского спек­ тра методом HQRB Рис. 2.2. Сравнение численных (кружки) и теоретических (квадраты) оценок МХПЛ для 11 спутников 2.5.2. Случай пространственного вращения Определим систему координат,, с началом в перицентре орбиты спутника следующим образом: ось направлена по прямой «перицентр – пла­ нета», ось параллельна вектору орбитальной скорости, ось ортогональна плоскости орбиты и дополняет систему до правой. Ориентация спутника в этой системе определяется последовательностью воображаемых поворотов на углы Эйлера,, из исходного положения до тех пор, пока спутник не ока­ жется в своем истинном положении. В исходном положении оси инерции,, направлены по осям,,, соответственно. Последовательность вообра­ жаемых поворотов, принятая в данной диссертационной работе, следующая:

во-первых, осуществляется поворот на угол вокруг оси, затем — на угол вокруг оси и, окончательно, на угол вокруг оси. Эти повороты наглядно представлены на рис. 2.4.

Такое определение углов Эйлера идентично используемому Уиздомом и др. в [70] и отличается от обычного. Причина в том, что стандартная си­ стема сингулярна, когда ось спутника (ось наибольшего момента инерции) Рис. 2.3. Мгновенное значение МХПЛ (полученное методом теневой траектории; пунктир­ ная линия) и мгновенные значения компонент ляпуновского спектра (полученные методом HQRB; сплошные линии) для Гипериона в зависимости от времени; плоский случай ортогональна плоскости орбиты. В принятой здесь системе это положение отвечает условию = 0, а сингулярность смещена в точку = ±/2. По­ следнее значение отвечает ситуации, когда ось вращения спутника лежит в плоскости орбиты. Как показано в [44], используемая здесь система связана со стандартной, [0, 2), [0, ] соотношениями Примем систему единиц = = = = 1, где — постоянная всемирного тяготения, — масса планеты, — большая полуось орбиты спутника, а — среднее движение. Тогда время измеряется в единицах периода обращения — так же, как в разделе 2.5.1. Вектор угловой скорости спутника, вращающегося в поле тяготения планеты, определяется динамиче­ Рис. 2.4. Определение углов Эйлера в инерциальной системе отсчета скими уравнениями Эйлера (см. напр. [4, 70]), которые в выбранных единицах приводятся к виду (см. [43]), где — расстояние «спутник – планета»;,, — компонен­ ты вектора угловой скорости в системе,, ;,, — направляющие косинусы главных осей инерции по отношению к направлению на планету.

Согласно [70], в принятой нами системе координат выражения для направля­ ющих косинусов следующие:

Следует отметить, что выражения для моментов сил тяготения (правые части уравнений (2.20)) верны в предположении, что размеры спутника малы по сравнению с, а его масса = 1.

Кинематические уравнения Эйлера в той же системе координат, соглас­ но [70], сводятся к Расстояние между спутником и планетой равно где — эксцентриситет орбиты; эксцентрическая аномалия получается ре­ шением уравнения Кеплера где — время. Истинная аномалия для подстановки в уравнения (2.21) выражается соотношением Уравнения Эйлера (2.20), (2.22) и уравнения орбитального движения (2.23)–(2.25) совместно образуют замкнутую систему, определяющую враща­ тельное движение спутника. Совместное интегрирование этой системы поз­ воляет определить эволюцию ориентации спутника, то есть углы,, как функции времени.

Численные оценки ХПЛ в данном разделе получаются интегрированием указанной системы уравнений при следующих начальных условиях: |=0 = 1.5, |=0 = 0.001, |=0 = 0.001, |=0 = 1, |=0 = 0, |=0 = 0. Таким обра­ зом, для того, чтобы вывести движение из плоскости орбиты, к начальным значениям |=0 и |=0 добавляется малое смещение. Эксцентриситет спут­ ника и отношения /, / его моментов инерции являются параметрам задачи.

Мы используем и метод теневой траектории, и метод HQRB. Интегра­ тор и величина шага итерации те же самые, что и в плоском случае. В про­ странственном случае, однако, требуется большее число итераций, поскольку структура хаотических областей фазового пространства в этом случае гораз­ до сложнее, и поэтому требуется гораздо большее время, чтобы () достиг насыщения. Необходимо проводить усреднение на бльших временных мас­ штабах, чтобы получить истинную оценку ХПЛ. С другой стороны, сумма (6), служащая индикатором накопления ошибок, растет при этом быстрее.

Согласно нашему опыту, зависимость достигает насыщения после = 5 · для всех спутников. Поэтому результирующие оценки ХПЛ усредняются на интервале [5 · 105, 106 ].

В табл. 2.3 представлены оценки МХПЛ пространственного вращения для указанных в разделе 2.5.1 одиннадцати спутников, полученные числен­ ным моделированием методом теневой траектории.

В табл. 2.4 представлены результаты вычисления полного ляпуновского спектра для случая пространственного вращения. Параметры интегрирова­ ния те же, что и в разделе 2.5.1. Последний столбец содержит значения сум­ мы (6) (индикатора накопления ошибок) на последнем шаге итерации.

Отметим, что вклады отдельных парных сумм () + (+ ) ( = 1,..., 3) в полную сумму (6) оказались одного порядка, так что внутреннюю ошибку Таблица 2.3. Случай пространственного вращения. Инерционные параметры, эксцентри­ ситет орбиты и численные оценки МХПЛ методом теневой траектории Таблица 2.4. Случай пространственного вращения. Численные оценки компонент ляпунов­ ского спектра методом HQRB и сумма всех ХПЛ каждого можно приблизительно принять равной /3.

Как и в разделе 2.5.1, здесь приведены зависимости мгновенных оценок ХПЛ от времени для Гипериона (рис. 2.5). На рис. 2.6 показан график зависи­ мости суммы (6) от времени; можно видеть, как ошибки накапливаются на больших временах интегрирования.

Полученный здесь спектр ХПЛ для Гипериона можно сравнить с резуль­ татами Уиздома и др. [70], хотя нами использованы несколько другие исход­ ные данные для расчета моментов инерции этого спутника [64]. На рис. 9 рабо­ ты [70] показана полученная Уиздомом и др. зависимость текущих значений ХПЛ от времени интегрирования для случая пространственного вращения Гипериона. Из этого рисунка можно заключить, что имеется близкое соответ­ ствие результатам нашего численного моделирования, хотя детали графиков, конечно, различаются — главным образом, из-за различного времени инте­ Рис. 2.5. Мгновенное значение МХПЛ (полученное методом теневой траектории; пунктир­ ная линия) и мгновенные значения компонент ляпуновского спектра (полученные методом HQRB; сплошные линии) для Гипериона в зависимости от времени; пространственный случай грирования.

2.5.3. Параметры Родрига–Гамильтона Задача пространственного вращения спутника, обращающегося по фик­ сированной орбите, может быть решена также с использованием параметров Родрига–Гамильтона (РГ) 0, 1, 2, 3 [5]. В отличие от используемой нами системы углов Эйлера, содержащей координатную сингулярность при = ± (см. выше), параметры РГ регулярны во всей области задания ( [1, 1]).

Размерность динамической системы = 7, и ее состояние определяется на­ бором переменных (,,, 0, 1, 2, 3 ).

Рис. 2.6. Сумма всех ХПЛ для Гипериона в зависимости от времени; пространственный случай Параметры РГ определяются через углы Эйлера следующим образом:

При этом обратное преобразование задается следующими соотношениями:

Соотношения (2.26) полезны для задания начальных данных интегрирования в форме углов Эйлера, имеющих наглядную геометрическую интерпретацию, а соотношения (2.27) — для контроля траектории в ходе интегрирования.

Согласно [5], система (2.20), (2.22) в параметрах РГ записывается следу­ ющим образом:

с направляющими косинусами где истинная аномалия определяется уравнением (2.25).

Для ХПЛ системы (2.28) выполняются соотношения (8) = () ( = 1,..., 4). Отсюда следует, в частности, что (4) = 0, что, наряду с условием нормировки 3 2 = 1, может служить для контроля вычислений.

Все проведенные нами численные эксперименты показывают хорошее согласие результатов интегрирования, выполненного в системе углов Эйле­ ра и в системе параметров РГ. В качестве иллюстрации рассмотрим задачу свободного вращения тела (случай Эйлера), когда члены с 1/3 в правых ча­ стях уравнений (2.20) и (2.28) равны нулю. В первом случае на интервале [0, 108 ] интегрируется система (2.20), (2.22). Во втором случае интегри­ руется система (2.28), а значения параметров РГ преобразуются для сопо­ ставления к соответствующим им значениям углов Эйлера при помощи фор­ мул (2.27). На рис. 2.7 приводятся графики разностей компонент вектора уг­ ловой скорости,, и углов Эйлера,,, вычисленных непосредствен­ но и через параметры РГ. Видно, что на временах до 3.5 · 107 различие пренебрежимо мало, на уровне относительной машинной точности. Однако и на бльших временах различие имеет малую величину; постепенный его рост объясняется, по-видимому, накоплением ошибок округления.

Рис. 2.7. Сравнение результатов интегрирования для случая свободного вращения с ис­ пользованием углов Эйлера и параметров РГ (помечены индексом «RH »); показаны за­ висимости от времени разностей компонент вектора угловой скорости и углов Эйлера, которые в случае использования параметров РГ вычисляются по формулам (2.27) На рис. 2.8 приведено сравнение МХПЛ, вычисленных с использованием углов Эйлера и параметров РГ для нескольких значений параметров /, / и при одних и тех же (в терминах углов Эйлера) начальных данных.

Видно хорошее согласие значений МХПЛ, определенных обоими методами, которое подтверждает, что численная устойчивость метода в случае регуляр­ ного вращения (рис. 2.7) имеет место также для хаотического движения. Это служит дополнительным указанием на надежность полученных в рамках дан­ ной диссертационной работы оценок ХПЛ.

Рис. 2.8. Сравнение МХПЛ, вычисленных с использованием углов Эйлера и параметров РГ (помечены индексом «RH »), для модельных значений инерционных и орбитальных па­ раметров 2.6. Сравнение случаев плоского и пространственного вращения Сравним, прежде всего, результаты вычислений МХПЛ методом теневой траектории и методом HQRB.

На рис. 2.9а приведен график зависимости оценки (1), полученной ме­ тодом HQRB, от оценки МХПЛ, полученной методом теневой траектории, для выбранных нами 11 спутников в случае плоского вращения. Значения (1) взяты из табл. 2.2, а значения — из табл. 2.1. «Усы» у каждой точки (на этом и всех последующих графиках) отражают максимальные отклоне­ ния мгновенных значений ХПЛ на плато от среднего значения по интервалу усреднения ( [5 · 104, 105 ] для плоского случая и [5 · 105, 106 ] — для пространственного).

Прямая, проведенная на рисунке, является линейной аппроксимацией (1) = ·+. Можно сделать заключение о хорошем согласии результатов, полученных обоими методами: наклон прямой = 0.99 ± 0.07, а сдвиг = 0.002±0.004. Интервалы ошибок параметров приводятся для доверительного уровня 95%.

На рис. 2.9б показан аналогичный график для случая пространственно­ го вращения. Значения (1) взяты из табл. 2.4, а — из табл. 2.3. Наклон аппроксимирующей прямой составляет 0.75 ± 0.11, а сдвиг — 0.02 ± 0.01.

Видно, что статистическое согласие в пространственном случае хуже, чем в плоском. Этого и следовало ожидать, так как структура хаотической обла­ сти здесь гораздо сложнее, и простой метод теневой траектории оказывается, очевидно, менее надежным.

По сравнению с плоским случаем, получение корректных оценок ХПЛ теперь требует гораздо большего времени интегрирования. Траектория в про­ странственном случае много сложнее, она проходит через разнообразные об­ ласти фазового пространства с различными локальными значениями ХПЛ.

В обоих методах полное время интегрирования ограничено накоплением оши­ бок. Если в случае метода теневой траектории оно вызвано ренормировками, то в методе HQRB его определяет отличие точной матрицы касательного отображения от ее численной аппроксимации по формуле (2.15) или (2.17).

Обратимся теперь к важному вопросу о том, может ли описанная вы­ ше (и, очевидно, верная в плоском случае) теория аналитического оценива­ ния МХПЛ быть применена также и к пространственному случаю. На этот Рис. 2.9. МХПЛ, вычисленные двумя различными методами, для 11 спутников: (а) плос­ кий случай, (б) пространственный случай вопрос помогает ответить график сравнения МХПЛ (1) для двух этих слу­ чаев, приведенный на рис. 2.10а. Легко видеть, что какая-либо корреляция отсутствует. При этом интересно отметить, что для (2) и (3) глобальная корреляция между плоским и пространственным случаем, по-видимому, су­ ществует (см. рис. 2.10б,в, хотя зависимости и нелинейны.

Сопоставляя значения параметров спутников в нашей выборке, можно очень приблизительно определить границы применимости теории как > 0.02 и / > 0.8. Если характеристики спутника лежат вне этого диапазо­ на, внешние резонансы с явной зависимостью от времени не играют большой роли во вращательной динамике; поскольку спутник существенно асиммет­ ричен, важнее всего взаимодействие резонансов связи, проявляющихся уже в круговой задаче. Соответствующая гамильтонова система тогда близка к автономной. Для того, чтобы вывести аналитические оценки МХПЛ и в этом случае, требуется выяснить, какой из внутренних резонансов связи является ведущим (определение ведущего резонанса см. в [30]). Другой способ добить­ ся успеха в получении теоретических оценок — это исследовать зависимость МХПЛ от константы Якоби системы в антиинтегрируемом пределе. Такое исследование проводится далее в главе 3.

Рассмотрим теперь статистические соотношения между компонентами ляпуновского спектра по данным таблиц 2.2 и 2.4. Соответствующие зависи­ мости показаны на рис. 2.11. На этом рисунке квадраты обозначают оценки (2) в зависимости от (1), а треугольники — (3) от (1). Как можно видеть из рис. 2.11, значения компонент приблизительно пропорциональны друг другу.

Пропорциональность наблюдается и в плоском, и в пространственном слу­ чае. Средние по всем спутникам отношения компонент ляпуновского спектра в плоском случае составляют (1) /(2) 1.8, (1) /(3) 24 и (2) /(3) 13.

Доверительные интервалы (на уровне 95%) для них [1.53, 2.06], [9.68, 38.4] и [6.74, 19.2], соответственно. Средние отношения ХПЛ в пространственном случае (1) /(2) 3.7, (1) /(3) 22 и (2) /(3) 5.4, 95%-ные доверитель­ ные интервалы для них: [2.31, 5.18], [8.61, 35.3] и [4.56, 6.15], соответственно.

В заключение приведем результаты регрессионного анализа для ком­ понент ляпуновского спектра (1), (2), (3) в зависимости от инерционных параметров /, / и 0. Найденные отношения могут быть полезными, пока не построена теория для аналитического оценивания ХПЛ в данной задаче. Коэффициенты линейной аппроксимации даются вместе с их 95% доверительными интервалами; — коэффициент корреляции.

Рис. 2.10. Сравнение компонент ляпуновского спектра в плоском и пространственном слу­ чаях: (а) (1), (б) (2), (в) (3) Рис. 2.11. Отношения компонент ляпуновского спектра: (а) плоский случай, (б) простран­ ственный случай (1) = (0.065 ± 0.091)0 + (0.044 ± 0.073), = 0.48, (3) = (0.008 ± 0.004)0 + (0.000 ± 0.004), = 0.82.

Наилучшие корреляции наблюдаются для /. Это понятно, поскольку мо­ менты инерции определены таким образом, что ; поэтому / характеризует максимальную асимметрию. Корреляция для / гораздо ху­ же. Зависимости для 0 важнее всего в плоском случае. Поскольку 0 зави­ сит и от /, и от /, корреляция носит промежуточный характер между зависимостями для / и для /. Отметим, что в случаях / и 0 кор­ реляция лучше для компонент спектра более высокого порядка. И снова это вполне объяснимо: для некоторых спутников роль эксцентриситета значима, а это сказывается на МХПЛ.

2.7. Выводы ко второй главе В данной главе были рассмотрены методы численного и аналитического оценивания ХПЛ вращения спутников планет. Величина, обратная МХПЛ (ляпуновское время), составляет характерное время предсказуемости дина­ мики. В приложении к изучению вращательной динамики спутников планет определение этой величины может оказаться полезным для планирования наблюдений кривых блеска спутников, характер вращательной динамики ко­ торых пока не известен.

Используемая модель включает асимметричное (трехосное) твердое те­ ло, вращающееся вокруг своего центра масс и обращающееся по фиксирован­ ной эллиптической орбите. Рассматривается и плоское, и пространственное вращение. В плоском случае ось вращения спутника совпадает с осью наи­ большего момента инерции и ортогональна плоскости орбиты. В простран­ ственном случае спутник может вращаться произвольно.

В рамках данной модели для 11 избранных спутников с известными из наблюдений инерционными и орбитальными параметрами вычислены полные ляпуновские спектры хаотического движения. Использован алгоритм, осно­ ванный на QR-факторизации матрицы касательного отображения [66]. Для сравнения использован также более традиционный метод «теневой траекто­ рии», позволяющий вычислить только МХПЛ. Значения МХПЛ, получен­ ные обоими методами, находятся в хорошем согласии. Выполнено сравнение результатов численного интегрирования с использованием углов Эйлера с результатами, полученными на основе параметров Родрига–Гамильтона, сво­ бодных от координатных сингулярностей; совпадение результатов для обоих методов служит дальнейшим подтверждением надежности применяемых ал­ горитмов численного расчета ХПЛ.

Можно ли для рассматриваемой задачи предсказать значения ХПЛ ана­ литически? Для того, чтобы исследовать этот вопрос, было рассмотрено при­ ложение аналитического метода, основанного на теории сепаратрисных отоб­ ражений [18, 19]. Приведены дополнительные свидетельства того, что данный метод дает хорошие результаты при оценивании значений МХПЛ в плоской задаче; это иллюстрируют рис. 2.1 и 2.2.

Однако метод аналитического оценивания компонент ляпуновского спек­ тра в пространственной задаче во всем диапазоне значений инерционных и орбитальных параметров по-прежнему остается предметом для дальнейших исследований. В настоящее время можно аналитически предсказать только значения МХПЛ в некоторых интервалах этих параметров. Показано, что ме­ тод И. И. Шевченко [18, 19] применим для наименьших (среди рассмотренных нами) значений динамической асимметрии или/и наибольших эксцентрисите­ тов орбиты. В этом случае ведущим является синхронный резонанс, и теория сепаратрисных отображений, разработанная в применении к движению вбли­ зи сепаратрис этого резонанса, обеспечивает оценки МХПЛ, которые согла­ суются с результатами численного моделирования.

Получены также простые эмпирические зависимости компонент ляпу­ новского спектра от параметров задачи (/, /, 0 ).

Хаотическое вращение спутников планет:

ляпуновские спектры и константа Якоби 3.1. Введение Даже в простейшем реалистичном случае спутника — трехосного твердо­ го тела на фиксированной эллиптической орбите — вращательная динамика спутника представляет множество интересных динамических явлений, как регулярных, так и хаотических [9, 69, 70]. Пример регулярного движения да­ ет синхронный резонанс, в котором находятся многие спутники. Существует два основных типа синхронного резонанса — - и -резонансы [10]. Редким режимом движения в синхронном резонансе является бифуркационный ре­ жим [8]. Все это — регулярные типы движения; хаотическое поведение про­ является в окрестностях сепаратрис резонансов. В реальности выявляются и регулярные, и хаотические режимы вращения — прежде всего, из кривых блеска [2, 26, 42, 46]. Новые открытия неизвестных ранее малых спутников планет (см. напр. [34, 35]) служат обширным источником новой информации для интерпретации существующих теорий и одновременно создают потреб­ ность в новых теоретических исследованиях вращательной динамики спут­ ников.

В главе 2 были вычислены полные спектры ХПЛ хаотического враще­ ния для набора из 11 спутников при помощи метода HQRB, предложенного фон Бременом и др. [66]. Для вычисления МХПЛ использовался также более традиционный метод теневой траектории. Рассматривались пространствен­ ная (трехмерная) и плоская (с осью вращения, ортогональной плоскости орби­ ты) модели вращения. Численные оценки ХПЛ, полученные в плоском и про­ странственном случаях хаотического вращения, сопоставлялись с аналитиче­ скими оценками, даваемыми теорией сепаратрисных отображений [18, 19] в модели нелинейного резонанса (в данном случае — синхронного спин-орби­ тального резонанса) как возмущенного нелинейного маятника. Были полу­ чены свидетельства хорошего согласия данных численного моделирования с теорией сепаратрисных отображений в плоском случае. Более того, было показано, что разработанная для плоского случая теория, вероятнее всего, сохраняет свою применимость и в случае пространственного вращения, если динамическая асимметрия спутника достаточно мала и/или эксцентриситет его орбиты достаточно велик (но не настолько, чтобы динамическая модель стала неприменимой).

Теория какого типа будет применима во всем диапазоне значений пара­ метров, в настоящее время неизвестно. Неясна также роль начальных усло­ вий в пространственной задаче, которая имеет три с половиной степени свобо­ ды — то есть на две степени свободы больше, чем в плоском случае. В данной главе делается попытка дальнейшего прояснения перспектив аналитического оценивания ХПЛ хаотического вращения спутников. Объектом исследования служит та же модельная задача, что и в главе 2: спутник представлен трехос­ ным твердым телом на фиксированной эллиптической или круговой орбите.

Метод исследования заключается в вычислении ляпуновских спектров в за­ даче пространственного хаотического вращения и анализе наблюдаемых за­ висимостей ХПЛ от инерционных и орбитальных параметров спутника и от значения константы Якоби динамической системы.

Рассматриваемый в данной главе набор реальных спутников включает 12 спутников Марса, Юпитера, Сатурна и Нептуна. По сравнению с преды­ дущей главой, в набор включен один дополнительный спутник, Адрастея.

Все это — по-прежнему спутники с известными значениями инерционных и орбитальных параметров.

В данной главе, во-первых, вычисленные МХПЛ сравниваются с МХПЛ, полученными для эксцентриситета орбиты, формально положенного равным нулю. Это позволяет получить информацию о границах применимости тео­ рии, основанной на методе сепаратрисного отображения [18, 19], в простран­ ственном случае. Во-вторых, изучается вопрос о применимости метода се­ паратрисных отображений для аналитического оценивания МХПЛ в случае вытянутого осесимметричного спутника. Этот случай важен, поскольку неко­ торые спутники с известной формой близки к вытянутому осесимметричному случаю (см. данные в табл. 3.1). В третьих, исследуется зависимость компо­ нент ляпуновского спектра от значения интеграла Якоби системы для спут­ ника произвольной формы на круговой орбите. В результате будет предло­ жен дополнительный метод аналитического оценивания ХПЛ, работающий в тех случаях, когда начальные условия играют в определении ХПЛ замет­ ную роль. Как будет показано, эта роль является определяющей при высоких значениях константы Якоби системы.

3.2. Вычисление ляпуновских спектров Моделью спутника в данной главе служит трехосный эллипсоид одно­ родной плотности, обращающийся по фиксированной эллиптической или кру­ говой орбите, движение которого подчиняется системе уравнений (2.20)–(2.25), определенной в разделе 2.5.2. Спектры ХПЛ вычисляются посредством мето­ да HQRB, описанного в разделе 2.4, с применением того же самого интегра­ тора [39].

Имея вектор начального положения x0 = (0, 0, 0, 0, 0, 0 ), мы вы­ числяем итеративно вектор положения x = (,,,,, ) на равномер­ ной сетке = · ( = 1,..., ) путем интегрирования уравнений движе­ ния (2.16), являющихся уравнениями Эйлера (2.20, 2.22).

В данной главе используется метод вычисления матрицы касательного отображения (, = 1,..., 6), отличный от применявшегося в главе 2. А именно, вместо использования прямого численного приближения для она вычисляется путем интегрирования вариационных уравнений где — матрица Якоби исходной нелинейной системы (2.16). Начальным условием для в (3.1) служит = (единичная матрица) на каждой ите­ рации. Уравнения (2.16) и (3.1) интегрируются совместно. Таким образом, полная система суть ложения», (X, ) = (x, ) для первых 6 переменных и (X, ) = (x, ) для остальных. Система (3.2) имеет порядок 6 + 36 = 42.

Как было обнаружено в ходе сопоставления с методом вычисления мат­ рицы касательного отображения, использованным в главе 2, метод, основан­ ный на уравнениях (3.1), (3.2), дает значительное преимущество в точности, поскольку он не требует введения дополнительного произвольного малого параметра, отвечающего длине векторов смещения, аппроксимирующих ис­ тинные касательные векторы. Тем самым повышается возможная длитель­ ность интервала интегрирования. Более того, несмотря на существенно более высокий порядок интегрируемой системы, данный метод также приводит к возрастанию скорости расчетов приблизительно вдвое — благодаря тому, что в нем не нужно выполнять 2 дополнительных итераций исходного отобра­ жения (см. раздел 2.4).

Описанная процедура позволяет получить оценки ХПЛ на каждой ите­ рации. Эти оценки испытывают сильные колебания в начале интервала интегрирования. Поэтому, как и ранее, значения ХПЛ можно считать окон­ чательными, если они оставались стабильными в течение достаточно долго­ го времени («вышли на плато»). Для повышения точности в данной главе осуществляется усреднение по большой части времени интегрирования.

Как оказалось, для всех вычислений в данной главе достаточно интервала = 106. В качестве «истинного» значения ХПЛ принимается среднее по второй половине этого интервала, то есть [5 · 105, 106 ].

Аналогично, в качестве теста внутренней устойчивости метода исполь­ зуется сумма всех ХПЛ () + (7) ( = 1, 2, 3), которая должна равняться нулю для гамильтоновой системы. Данный тест показал, что, начиная с мо­ мента достижения плато, эти суммы как правило не превышают 0.01–0.02% от среднего значения на плато. Этот результат указывает на гораздо более высокую (примерно на два порядка) точность этого метода по сравнению с тем, который был использован в главе 2. Некоторые другие тесты надежности вычислений представлены ниже в разделе 3.6.

орбитального периода. МХПЛ в расчете на сутки получается умножением этих значений на orb, где период обращения orb выражен в сутках.

3.3. Орбитальное движение и ХПЛ вращения Методом, описанным в разделe 3.2, вычислены ляпуновские спектры для выборки из 12 малых спутников планет с известными орбитальными и инер­ ционными параметрами. Поскольку целью данной главы является исследо­ вание возможности хаотического вращения спутника в ходе его динамиче­ ской истории, а не его фактическое состояние в настоящее время, начальные данные для интегрирования выбираются внутри хаотического слоя фазового пространства: 0 = 1.5, 0 = 0 = 0.001; 0 = 1, 0 = 0 = 0.1 Эти данные при­ близительно соответствуют неустойчивому верхнему положению маятника в модели синхронного резонанса [18, 19] (см. также главу 2). Начальная точка берется вне плоскости орбиты (0 = 0 = 0.001) для того, чтобы обеспечить пространственное вращение.

Результаты расчетов представлены вместе с параметрами спутников в табл. 3.1. Отношения /, / и эксцентриситеты орбит взяты из следую­ щих источников: [59] (Фобос); [69] (Деймос, Амальтея, Янус и Эпиметей); [37] (Адрастея); [64] (Гиперион); [60] (Елена); [27] (Атлас); [36] (Прометей и Пан­ дора); [32] (Протей). В случае, когда в источнике приведены геометрические размеры,,, отношения /, / вычисляются в предположении, что спутник является трехосным эллипсоидом однородной плотности, то есть (см. напр. [5]).

Влияние внешних (спин-орбитальных) резонансов — и, в частности, син­ хронного резонанса — на ХПЛ должно возрастать с увеличением эксцен­ триситета орбиты и с уменьшением геометрической асимметрии спутника (см. [56]). И наоборот, естественно предположить, что хаотическая динамика пространственного вращения существенно асимметричного спутника на по­ чти круговой орбите определяется, главным образом, взаимодействием внут­ ренних резонансов связи. В настоящее время отсутствует какая-либо теория аналитического оценивания ХПЛ в этом случае.

В данной главе вычислены спектры ХПЛ для выборки из 12 спутников с эксцентриситетом орбиты, формально положенным равным нулю. Резуль­ таты представлены в табл. 3.1. Заметим, что наблюдаемый эксцентриситет Для двух спутников, Атласа и Протея, использованы несколько иные начальные условия: 0 = 1.51 и 0 = 1.52, так как точка 0 = 1.5, которая подходит для остальных спутников, для этих двух лежит вне хаотической области при = 0.

Таблица 3.1. Инерционные параметры, эксцентриситеты орбит и спектры ХПЛ Фобос (M1) 0.7234 0.8504 0.01500 0.085 0.027 0.0057 0.124 0. Деймос(M2) 0.6612 0.8808 0.00050 0.161 0.048 0.0019 0.149 0. Амальтея (J5) 0.5049 0.9371 0.00300 0.102 0.032 0.0028 0.117 0. Адрастея (J14) 0.6098 0.8293 0.00000 0.172 0.030 0.0000 0.172 0. Гиперион (S7) 0.6220 0.8840 0.10000 0.054 0.019 0.0035 0.154 0. Янус (S10) 0.7302 0.8757 0.00900 0.082 0.025 0.0040 0.127 0. Эпиметей (S11) 0.7055 0.9081 0.00700 0.067 0.019 0.0021 0.124 0. Елена (S12) 0.8293 0.9466 0.00500 0.046 0.013 0.0015 0.072 0. Атлас (S15) 0.7493 0.8220 0.00200 0.137 0.028 0.0027 0.073 0. Прометей (S16) 0.4124 0.8723 0.00400 0.148 0.056 0.0043 0.182 0. Пандора (S17) 0.5717 0.8622 0.00400 0.143 0.050 0.0056 0.176 0. Протей (N8) 0.9215 0.9685 0.00046 0.043 0.004 0.0002 0.038 0. орбиты Адрастеи и так нулевой. В случае = 0 наша гамильтонова система автономна, и компонента =0 = 0.

В графическом виде результаты для ненулевого и нулевого эксцентри­ ситетов сопоставлены на рис. 3.1. «Усы» на рис. 3.1 — так же, как и всю­ ду в этой главе — отражают максимальные отклонения мгновенной оценки ХПЛ от среднего значения на интервале [5 · 105, 106 ]. Прямые линии представляют собой линейные аппроксимации =0 = () + ( = 1, 2), где = 0.571 ± 0.650, = 0.066 ± 0.069 для = 1 и = 0.543 ± 0.338, = 0.009 ± 0.011 для = 2. Статистические погрешности приведены для доверительного уровня 95%, как и всюду в настоящей работе.

Какая-либо корреляция на рис. 3.1а отсутствует (коэффициент корре­ ляции = 0.552). Это свидетельствует о том, что эксцентриситет орбиты и, тем самым, взаимодействие спин-орбитальных резонансов играют важней­ шую роль в определении (1) для рассмотренной выборки спутников. С дру­ гой стороны, некоторая корреляция (линейная зависимость) присутствует на рис. 3.1б ( = 0.771). Это говорит о том, что в определении (2) решающая роль принадлежит внутренним резонансам связи, чье доминирование обу­ словлено наличием дополнительных степеней свободы по сравнению с плос­ ким случаем.

По данным табл 3.1 можно получить простые эмпирические зависимо­ сти, аналогичные приведенным в разделе 2.6. Отличие от оценок (2.30)–(2.32) в численном методе определения ХПЛ и небольшом расширении набора спут­ ников.

Рис. 3.1. Сопоставление численных оценок компонент ляпуновского спектра для = 0 с оценками для реального значения эксцентриситета; (а) (1), (б) (2) (2) = (0.271 ± 0.133)(1) + (0.001 ± 0.015), = 0.82;

Статистическая значимость большинства этих оценок низка, в согласии с выводами главы 2. В разделе 3.5 будет показано, что корреляция (1) с вели­ чиной константы Якоби системы для данных значений / и / в авто­ номном случае гораздо лучше, так что по сравнению с приведенными выше зависимостями (1) от / и / зависимость (1) от значения константы Якоби представляется гораздо более предпочтительным средством оценива­ ния МХПЛ.

Следует отметить, однако, достаточно явную корреляцию ( = 0.82) (2) с (1). При этом (3) не коррелирует с (1) ( = 0.08), что может сви­ детельствовать о больших относительных ошибках определения (3) в силу малой его величины.

3.4. Аналитическое оценивание ХПЛ: границы применимости В численных экспериментах, проведенных в главе 2, получены предвари­ тельные указания на применимость теории сепаратрисных отображений для аналитического оценивания МХПЛ хаотического пространственного враще­ ния спутников: 0.02 и / 0.8.