Доклады Российской академии наук. Математика, информатика, процессы управления, 2023, T. 512, № 1, стр. 33-41

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим круговую пространственную ограниченную задачу трех тел. Предположим, что пассивно гравитирующее тело (астероид) P находится в гравитационном поле двух массивных тел, движущихся друг относительно друга по круговой орбите радиуса ${{r}_{J}}$; центральное тело S (Солнце) имеет массу ${{m}_{S}}$ и воздействует на астероид с силой ${{{\mathbf{F}}}_{S}}$, второе тело J (Юпитер) массы ${{m}_{J}}$ оказывает возмущающее действие с силой ${{{\mathbf{F}}}_{J}}$. Считаем, что невозмущенная траектория спутника есть кеплеровский эллипс с фокусом в S, плоскость которого $\Pi $ образует с плоскостью движения ${{\Pi }_{0}}$ притягивающих тел угол i (рис. 1).

Известно, что в такой классической постановке задача исследована достаточно подробно. Так, в статьях [1, 2] получено первое (хилловское) приближение для дважды усредненной силовой функции задачи, подробно проведен количественный и качественный анализ уравнений движения в указанном приближении. В статье [3] получено усеченное значение усредненной силовой функции в четвертом приближении, однако анализ уравнений движения ограничен хилловским случаем. В работе [4] исследованы усредненные уравнения второго приближения, делается вывод о топологической эквивалентности фазового портрета хилловскому случаю. В статье [5] приводится подробный численный анализ редуцированных уравнений движения, когда усредненная силовая функция рассматривается в виде кратного интеграла по двум быстрым переменным.

Основная цель статьи – получить новые результаты в двукратно усредненной задаче Гаусса, используя явное компактное представление усредненной силовой функции $R{\kern 1pt} *{\kern 1pt} *$ в виде степенного ряда, либо в виде ряда Фурье. Постановка задачи и результаты исследований отличаются от численных исследований М.А. Вашковьяка [5], когда $R{\kern 1pt} *{\kern 1pt} *$ рассматривают в виде квадратуры, непредставимой в конечном виде. Численный метод, несмотря на свою универсальность, ограничен в силу сложности функции $R{\kern 1pt} *{\kern 1pt} *$: невозможно рассчитать все многообразие стационарных точек при произвольных значениях параметров, их бифуркации, трудно построить гетероклинические траектории.

Приближенно-аналитический метод, используемый в статье, заменяет силовую функцию $R{\kern 1pt} *{\kern 1pt} *$ усеченной, что позволяет строить бифуркационные диаграммы, точки бифуркации равновесий, сепаратрисы, фазовые портреты колебаний по эксцентриситету $e$ и аргументу перицентра $\omega $ в высоких приближениях по малому параметру. При этом решается важный вопрос о поведении силовой функции на кривых неаналитичности, препятствующих численному расчету квадратуры $R{\kern 1pt} *{\kern 1pt} *$: степенные ряды и ряды Фурье, описывающие $R{\kern 1pt} *{\kern 1pt} *$, являются асимптотическими рядами на этих кривых.

Введем неинерциальную гелиоцентрическую систему координат с центром в $S$. Ось $Sx$ направим в точку весеннего равноденствия, $SN$ – линия узлов орбиты спутника, остальные оси, образующие правую систему координат, на рис. 1 не указаны. Пусть ${\mathbf{r}}$ – радиус-вектор тела P, ${{{\mathbf{r}}}_{J}}$ – радиус вектор тела $J$.

Пусть $\Omega $ – долгота восходящего узла невозмущенной орбиты спутника на плоскости ${{\Pi }_{0}}$, $e,\;\omega $ – эксцентриситет и аргумент перицентра этой орбиты, ${{\lambda }_{J}}$ – долгота тела $J$, а $\lambda $ – долгота тела P в плоскости ${{\Pi }_{0}}$.

Рассмотрим “внутренний” вариант задачи трех тел, когда параметры движения тела P удовлетворяют условию

Здесь $a$ – большая полуось невозмущенной орбиты спутника.

Запишем выражение для возмущенной силовой функции задачи:

где f – постоянная тяготения, $\gamma $ – угол между ${{{\mathbf{r}}}_{J}}$ и ${\mathbf{r}}$, $\nu $ – истинная аномалия в движении спутника вдоль невозмущенной орбиты.

Функция (1) может быть разложена в ряд по полиномам Лежандра с точностью до членов, не зависящих от координат тела $P$:

2. УСРЕДНЕНИЕ СИЛОВОЙ ФУНКЦИИ

Считаем, что частота ω_J$({{\lambda }_{J}} = {{\omega }_{J}}{\kern 1pt} t + {{\lambda }_{{j0}}})$ не резонирует с частотой n невозмущенного движения спутника. Проводя усреднение (3) по долготе ${{\lambda }_{J}}$ тела $J$, получим, с учетом (2), теоремы сложения для полиномов Лежандра и равенства ${{P}_{{2n + 1}}}(0)$ = 0, однократно усредненную силовую функцию R* задачи:

Выражение (4) совпадает с силовой функцией кольца Гаусса [6] при r < r_J, что было впервые установлено в работе Аксенова [7].

Усредним теперь выражение (4) по средней долготе $\lambda $ тела $P$, учитывая, что dλ = = ${{r}^{2}}d\nu {\text{/}}({{a}^{2}}\sqrt {1\, - \,{{e}^{2}}} )$:

Для вычисления этого интеграла представим выражение для r через истинную аномалию $\nu $, используя  формулы невозмущенного движения. Вычисления показывают, что функция R** примет вид

где

Основная техническая проблема исследований – вычисление этого интеграла. Современные программные комплексы, такие как Wolfram Mathematica и Maple, с этой задачей не справляются. Однако, учитывая периодичность по $\nu $ сомножителей, входящих в подынтегральное выражение интеграла (7), можем применить формулу Парсеваля [8], описывающую главный (секулярный) член произведения двух рядов Фурье:

Здесь $\{ {{a}_{0}},\;{{a}_{m}},\;{{b}_{m}}\} $ – коэффициенты Фурье функции ${{\left( {1 + e\cos \nu } \right)}^{{ - 2n - 2}}}$, $\{ {{\alpha }_{0}},\;{{\alpha }_{m}},\;{{\beta }_{m}}\} $ – коэффициенты Фурье функции ${{P}_{{2n}}}(\sin i\cos \theta )$.

Вычисления показывают, что интеграл (7) имеет вид

Здесь ${{F}_{{2,1}}}$ – гипергеометрическая функция Гаусса, а $F_{{3,2}}^{{{\text{reg}}}}$ – обобщенная гипергеометрическая регуляризованная функция Клаузена, $P_{{2n}}^{{(m)}}$ – присоединенная функция Лежандра. Отметим, что функция Клаузена $F_{{3,2}}^{{{\text{reg}}}}$ является аналитическим продолжением гипергеометрического ряда, представляет собой решение некоторого линейного дифференциального уравнения третьего порядка, имеет особенность в точках $z = 0,1,\infty $.

Достоверность формулы (9) была подтверждена численным интегрированием выражения (1 + + ecosν)^{–2n – 2}${{P}_{{2n}}}\left( {\sin i\cos \theta } \right)$ – результаты вычислений левой и правой частей формулы совпадают по всем значащим цифрам для следующего декартового произведения значений параметров: $e$ от 0 до 1 с шагом 0.1, $\omega $ от 0 до $2\pi $ с шагом $\pi {\text{/}}6$, $i$ от $ - \pi {\text{/}}2$ до $\pi {\text{/}}2$ с шагом $\pi {\text{/}}6$ и $n$ от 1 до 15 с шагом 1.

Выражение для усредненной силовой функции можно разложить в ряд Фурье по $\cos 2\omega $:

Здесь

Представим также усредненную силовую функцию в виде степенного ряда по малому параметру ${{(a{\text{/}}{{r}_{J}})}^{{2k}}}$:

где

Полученные разложения (10), (11) удобны с аналитической и вычислительной точки зрения. Более того, с каждым из этих разложений можно работать как с единой однозначной аналитической функцией (вычислять производные любого порядка по аргументам $e,\;\omega ,\;i$, проводить усечение рядов до членов любого порядка малости), если использовать Wolfram Mathematica.

Отметим, что аналитическое представление дважды усредненной силовой функции задачи в виде ряда Фурье по $\omega $ впервые получено в работе Аксенова [7], в которой усреднение проводилось по истинной аномалии. Такое усреднение отлично от усреднения по Гауссу. Коэффициенты этого ряда автор статьи выразил через неизвестные специальные функции, связанные между собой рекуррентными соотношениями. Заметим также, что иное представление степенного ряда (11) через специальные функции Хансена $X_{0}^{{k, - 2p}}(e)$ и функции наклона ${{F}_{{k0p}}}(i)$ содержится в статье [9] (см. также [13]). К сожалению, эти разложения имеют более громоздкий вид и менее удобны для аналитических исследований. Отметим также, что формула (10) переходит в ряд Фурье статьи [7], если усреднение по $\lambda $ заменить усреднением по $\nu $ с переходом к новому времени $d\tau = (\sqrt {fmp} {\text{/}}{{r}^{2}})dt$.

3. ВЕРИФИКАЦИЯ ДВУКРАТНО УСРЕДНЕННОЙ СИЛОВОЙ ФУНКЦИИ

Сравним формулу (10) для вычисления силовой функции задачи с результатами, получаемыми при численном усреднении ${{R}^{{(int)}}}$:

С этой целью возьмем следующие параметры для численного счета: $e = 0.2$, $i = - \pi {\text{/}}6$, $\omega = 0$, $a{\text{/}}{{r}_{J}} = 0.85$, $\Omega = \pi {\text{/}}12$.

Подсчитаем значения силовой функции R**, удерживая по возрастающей число $n$ членов ряда Фурье (10), и сравним эти значения с результатами численного интегрирования выражения (12) при том же самом числе $n$ членов подынтегрального выражения.

Как можно видеть из рис. 2, функция R** выходит на свое значение в фиксированной точке начиная с частичной суммы, содержащей 12 членов ряда (коричневая кривая), в то время как результат численного интегрирования представлен в виде колебаний вплоть до $n = 20$ (синяя кривая), с последующим выходом кривой на значение функции в точке. Показано, что при $n \geqslant 20$ результаты расчетов по формуле (10) совпадают с результатами численного интегрирования ${{R}^{{(int)}}}$ на сетке значений $e \in \left[ {0,0.9} \right]$ с шагом 0.1, $(a{\text{/}}{{r}_{J}}) \in $ $ \in $ [0.5, 0.95] с шагом $0.05$, $\omega \in \left[ {0,\pi } \right]$ с шагом π/6, $i \in \left[ { - \pi {\text{/}}2,\pi {\text{/}}2} \right]$ с шагом π/6.

Для верификации ряда (11) запишем его хилловское приближение, в котором удерживается только первый член ряда:

Здесь

Подставляя вышеперечисленное в выражение для $R_{{hill}}^{{**}}$ и упрощая, получим:

Это выражение тождественно совпадает с хилловским приближением усредненной силовой функции, полученным ранее в работе [1] (см. также [2]).

4. АСИМПТОТИЧЕСКОЕ ПОВЕДЕНИЕ УСРЕДНЕННОЙ СИЛОВОЙ ФУНКЦИИ НА КРИВЫХ НЕАНАЛИТИЧНОСТИ

Исследуем поведение рядов (10), (11) на кривых неаналитичности усредненной силовой функции, впервые обнаруженных в статье [10]:

Отметим, что в статьях [5, 10] эти кривые часто называют “сеператрисами”. Название условное, так как кривые не являются интегральными. На рис. 3 можно видеть зависимость значения силовой функции (10) в фиксированной точке фазового пространства оскулирующих элементов от числа удерживаемых членов ряда. Параметры, использованные для построения рисунка, связаны функциональной зависимостью ${{f}_{1}} = 0$, их значения $e = 0.134237$, $\omega = \pi {\text{/}}3$, $a{\text{/}}{{r}_{J}} = 0.95$, i = $ - \pi {\text{/}}3$. Частичная сумма первых 40 членов ряда (10) хорошо приближает силовую функцию задачи R**, и только удержание бóльшего количества членов ряда приводит к его расходимости. Таким образом, ряд (10) является асимптотическим по Пуанкаре [11], при этом расходимость ряда сопровождается нарастанием амплитуды колебаний.

Аналогичные результаты наблюдаются и для ряда (11). Таким образом, ряды (10), (11) расходятся вдоль кривых неаналитичности, но ряды удовлетворяют условию асимптотичности по Пуанкаре, поэтому мы можем их аппроксимировать аналитической функцией, представляющей собой сумму первых 40 членов ряда.

5. УСРЕДНЕННЫЕ УРАВНЕНИЯ ДВИЖЕНИЙ. ФАЗОВЫЕ ПОРТРЕТЫ КОЛЕБАНИЙ

Как хорошо известно [2, 12], уравнения в оскулирующих элементах имеют два первых интеграла:

Второй из этих интегралов принято называть интегралом Лидова–Козаи. С его помощью исключаем угол i

из усредненной силовой функции R**, получаем редуцированную систему уравнений с одной степенью свободы [5]:

Силовая функция $\hat {R} = \hat {R}\left( {e,\omega } \right)$ – результат исключения угла i. Очевидно, что интеграл энергии уравнений (15) имеет вид $\hat {R}\left( {e,\omega } \right) = {{c}_{2}}$.

Из второго интеграла (14) следует ограничение на область изменения оскулирующего эксцентриситета орбиты: $0 \leqslant e \leqslant \sqrt {1 - {{c}_{1}}} $. Здесь $e = 0$, $e = \sqrt {1 - {{c}_{1}}} $ – интегральные многообразия.

Построим фазовый портрет колебаний в плоскости переменных $(\omega ,e)$ редуцированной системы (15) при следующих значениях параметров: $a{\text{/}}{{r}_{J}}$ = = 0.8, ${{c}_{1}} = 0.1$. Учитывая, что хилловское приближение подробно исследовано (см. [1–3]), рассмотрим второе и четвертое приближение силовой функции (11).

Для начала удержим в $\hat {R}$ члены до четвертого порядка малости по $(a{\text{/}}{{r}_{J}})$ включительно (второе приближение усредненной силовой функции редуцированной системы):

Здесь

На рис. 4 изображен фазовый портрет колебаний в оскулирующих элементах $e,\omega $. Топологически он подобен хилловскому приближению, что соответствует выводу статьи [4] о качественном соответствии первого и второго приближений в усредненной задаче. В частном случае полярной орбиты (${{c}_{1}} = 0$), для которой угол наклона $i$ ее плоскости к плоскости эклиптики равен π/2, вывод о падении астероида на Солнце сохраняет силу. Мы также видим область $A$, в которой линия апсид совершает либрационные движения в окрестности единственной стационарной точки, как и в случае Хилла. Отличие этой области от области первого приближения состоит в том, что фазовые кривые в окрестности неподвижной точки имеют больший размах колебаний по эксцентриситету, а сама стационарная точка сдвигается в область бóльших значений эксцентриситета.

Исследуем теперь четвертое приближение силовой функции (11):

Здесь $\alpha = {{e}^{2}} - 1 + {{c}_{1}}$. Прежде всего отметим, что для обоснования усечений ряда (11) были проведены вычисления, при которых частичная сумма $R_{n}^{{**}}$ ряда (11) сравнивалась со значением интеграла (12), в котором число удерживаемых членов в подынтегральном выражении фиксировано и равно 45. Для параметров $e = 0.2$, $\omega = 0$, $a{\text{/}}{{r}_{J}} = 0.8$, $i = - \pi {\text{/}}6$, $\Omega = \pi {\text{/}}12$ показано, что сумма $R_{2}^{{**}}$ дает относительную ошибку приближения 0.06, а $R_{4}^{{**}}$ – близкую к нулю, хотя для $R_{5}^{{**}}$, $R_{6}^{{**}}$ она несколько портится; $R_{n}^{{**}}$ $(n \geqslant 10)$ практически совпадает с интегралом (12). Это дает частичное обоснование удержанию первых четырех членов ряда (11) при $a{\text{/}}{{r}_{J}} = 0.8$.

Кривая равновесий для четвертого приближения представлена на рис. 5. Она отвечает параметру $a{\text{/}}{{r}_{J}} = 0.8$ и значению аргумента перицента $\omega = \pi {\text{/}}2$, состоит из трех ветвей: ниспадающей, петлеобразной и зарождающейся петли. Ниспадающая ветвь присутствует в первом – четвертом приближениях и пересекает ось ${{c}_{1}}$ в точке (0, 0.847). Область параметров, лежащая выше ниспадающей кривой, не содержит равновесий и отвечает ротационным движениям. Рассмотрим область ниже этой кривой. Петля начинает формироваться во втором приближении и полностью формируется только в третьем приближении, в четвертом же приближении можно наблюдать появление следующей петли. Для петель характерно наличие точек бифуркации. Для сформировавшейся петли точка имеет координаты $c_{1}^{{(1)}} = 0.382,\,{{e}^{{(1)}}}$ = = 0.447 (на рис. 5 точка бифуркации обозначена цифрой 1). Для формирующейся петли точка имеет координаты $c_{1}^{{(2)}} = 0.015,\,{{e}^{{(2)}}} = 0.672$ (на рис. 5 точка обозначена цифрой 2). При прохождении через точку бифуркации появляются дополнительно два устойчивых равновесия и два неустойчивых. Последние являются седлами и они не представлены на рис. 5, так как отвечают случаю $\omega \ne \pi {\text{/}}2$. Эти равновесия можно видеть на рис. 6. Устойчивые равновесия имеют координаты $(\pi {\text{/}}2,0.74)$, $(\pi {\text{/}}2,0.8)$ и $(\pi {\text{/}}2,0.938)$.

Фазовый портрет колебаний изображен на рис. 6. Он состоит из четырех либрационных зон $A,\;C,\;D,\;E$ и одной области ротационных движений $B$. Зоны отделены друг от друга пятью гетероклиническими траекториями. Одна из них соединяет между собой неустойчивые равновесия, принадлежащие интегральному многообразию $e = 0$. Эту кривую мы наблюдаем также в хилловском случае. Остальные гетероклинические траектории соединяют между собой два седла, проявивших себя только в третьем приближении.

Итак, в четвертом приближении по малому параметру $a{\text{/}}{{r}_{J}}$ имеем зависимость фазового портрета колебаний от значений константы ${{c}_{1}}$ интеграла Лидова–Козаи. Так, при $0.382 < {{c}_{1}} < 1$ фазовый портрет соответствует хилловскому приближению; в случае $0.015 < {{c}_{1}} < 0.382$ его топология меняется за счет появления дополнительных равновесий фазовой точки и гетероклинических траекторий (рис. 6), а при $0 < {{c}_{1}} < 0.015$ топология снова меняется – появляется пара дополнительных равновесий.

Достоверность четвертого приближения вытекает из сравнения фазовых портретов колебаний – при фиксированных значениях константы ${{c}_{1}}$ – с третьего по пятое приближение включительно. На рис. 7 представлен фазовый портрет пятого приближения.

Топология сохраняется при ${{c}_{1}} = 0.1$. Она сохраняется и при других значениях ${{c}_{1}}$, меняясь только в точках бифуркации равновесий (точки 1 и 2) на рис. 5).

Численный анализ показывает, что удержание членов ряда более высокого порядка малости ведет к появлению дополнительных ветвей кривых равновесия в малой окрестности точки $e = \sqrt {1 - {{c}_{1}}} $, $\omega = 0$ (верхний левый угол фазового портрета) при малых ${{c}_{1}}$ , как следствие, меняется топология портрета в этой окрестности. Чем выше порядок малости добавляемых членов, тем меньше указанная окрестность по $e$. Случай ${{c}_{1}} \approx 0,e \approx \sqrt {1 - {{c}_{1}}} $ отвечает полярным сильно вытяным орбитам. Отметим, что при $e \to 1$ точность расчетов падает.

6. ЗАКЛЮЧЕНИЕ

Исследованы, с использованием схемы Гаусса двукратного усреднения, пространственные оскулирующие эллиптические движения астероида под действием гравитационного притяжения от Солнца и гравитационного возмущения от внешней планеты (Юпитера). Получены явные аналитические выражения для дважды усредненной силовой функции в виде ряда Фурье и степенного ряда. Аналитическое выражение позволяет быстро и эффективно строить любое усечение усредненной силовой функции, исследовать поведение этой функции вдоль кривых неаналитичности. Численно показано, что на этих кривых расходящиеся ряды усредненной силовой функции являются асимптотическими по Пуанкаре. Поэтому они допускают аппроксимацию аналитическими функциями в виде суммы $n$ первых слагаемых вплоть до $n = 40$.

Построены фазовые портреты колебаний во втором и четвертом приближениях. Показано, что второе приближение сохраняет все качественные выводы первого (хилловского) приближения. Показано также, что в четвертом приближении качественные выводы о фазовых колебаниях по эксцентриситету $e$ и аргументу перицентра $\omega $ хилловского случая сохраняются лишь при значениях константы интеграла Лидова–Козаи ${{c}_{1}}$ из интервала $(0.382,1)$. В противном случае, когда ${{c}_{1}} \in (0.015,\,0.382)$, имеем существенное усложнение картины колебаний: в силу бифуркации равновесий появляются дополнительные устойчивые и неустойчивые равновесия и гетероклинные траектории, из которых формируются границы четырех областей либрации и одной области ротационных движений по $e$ и $\omega $. Дальнейшее уменьшение ${{c}_{1}}$ вплоть до нуля ведет к еще большей сложности фазового портрета: появляются новые положения равновесия в силу бифуркации равновесий, новые гетероклинические траектории.

Показано, что удержание членов ряда (11) более высокого порядка малости ведет к появлению дополнительных ветвей кривых равновесия в малой окрестности точки $e = \sqrt {1 - {{c}_{1}}} $, $\omega = 0$ при малых ${{c}_{1}}$. Следовательно, четвертое приближение усредненной силовой функции дает топологически обоснованную картину колебаний за исключением малой окрестности указанной точки при ${{c}_{1}} \approx 0$.

Результаты этих исследований могут быть нарушены в очень высоких приближениях вдоль кривых неаналитичности. Поведение фазовых траекторий в окрестности кривых неаналитичности подробно исследовано в статье [5]. Седловые точки, не принадлежащие многообразию $e = 0$, впервые были получены в препринте [14].