Разложение гамильтониана планетной задачи в ряд Пуассона в гелиоцентрической системе отсчета

Микрюков Денис Викторович

Разложение гамильтониана планетной задачи в ряд Пуассона в гелиоцентрической системе отсчета

Строится разложение гамильтониана N-планетной задачи в ряд Пуассона по системе модифицированных (комплексных) элементов Пуанкаре в гелиоцентрической системе координат. Используется лагранжев и гамильтонов формализм. Приведены первые члены разложения главной и дополнительной части возмущающей функции. Численными экспериментами получены оценки числа слагаемых в представленных разложениях. Проведено сравнение с результатами других авторов.

Математика

Диссертации

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d368e5f1be77c40d59258

UUID: 11a515d7-cf1b-4c76-9ed1-67fbd1532aae

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 35

Микрюков Денис Викторович

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 537307 bytes

Поделиться работой

Санкт-Петербургский государственный университет Математико-механический факультет Кафедра небесной механики Микрюков Денис Викторович Разложение гамильтониана планетной задачи в ряд Пуассона в гелиоцентрической системе отсчета Выпускная квалификационная работа Научный руководитель: д.ф.-м.н., заведующий Кафедрой небесной механики, профессор Холшевников К.В. Рецензент: д.ф.-м.н., доцент УрФУ Кузнецов Э.Д. Санкт-Петербург 2016

SAINT-PETERSBURG STATE UNIVERSITY Mathematics & Mechanics Faculty Celestial Mechanics Chair Mikryukov Denis Viktorovich Expansion of the planetary Hamiltonian into a Poisson series in a heliocentric reference frame Final qualifying work Scientific supervisor: Head of Celestial Mechanics Chair, Professor K.V. Kholshevnikov Reviewer: Associated Professor E.D. Kuznetsov Saint-Petersburg 2016

Содержание 1 ВВЕДЕНИЕ 4 2 ОСНОВНАЯ ЧАСТЬ 6 2.1 Гелиоцентрические координаты . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Кинетическая и потенциальная энергия . . . . . . . . . . . . . . . . . . . . . 8 2.3 Уравнения Лагранжа и функция Гамильтона . . . . . . . . . . . . . . . . . . 10 2.4 Системы оскулирующих элементов . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5 Разложение гамильтониана . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.6 2.5.1 Разложение дополнительной части возмущающей функции . . . . . . 14 2.5.2 Разложение главной части возмущающей функции . . . . . . . . . . . 14 2.5.3 Вычислительный алгоритм разложения функции a0 /∆ . . . . . . . . . 16 Программное обеспечение для построения разложений . . . . . . . . . . . . . 20 2.6.1 Модуль M_KEPLER.py . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6.2 Модуль M_POISSON.py . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.6.3 Модуль M_LAPLACE.py . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.7 Результаты . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.8 Обсуждение результатов . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3 ЗАКЛЮЧЕНИЕ 26 4 ПОЛОЖЕНИЯ, ВЫНОСИМЫЕ НА ЗАЩИТУ 26 Список литературы 27 Приложение А 29 Приложение Б 30 Приложение В 31 Приложение Г 32 Приложение Д 33 3

1 ВВЕДЕНИЕ При решении задачи нескольких тел основными используемыми системами отсчета являются система координат Якоби и гелиоцентрическая система координат [14]. Первая является наиболее удобной для теоретического изучения задачи трех тел, например, для описания системы Солнце — Земля — Луна [14]. Однако в случае большого числа планет якобиевы координаты громоздки, так как положение каждого тела определяется положением всей предыдущей подсистемы тел [6]. Это приводит к необходимости разложения возмущающего потенциала системы в ряд по степеням малого параметра [19]. При изучении планетного движения гелиоцентрические координаты оказываются более удобными в практическом отношении. При любом числе планет главная часть возмущающей функции представляется суммой слагаемых, каждое из которых зависит от положения лишь пары планет [15]. Терминологическое замечание. Во многих исследованиях, посвященных построению планетных теорий [6, 9, 17], авторы используют при составлении исходных уравнений движения гамильтонов подход и вводят систему обобщенных координат и импульсов, которую называют каноническими гелиоцентрическими координатами1 . Эта система канонических координат была впервые рассмотрена Пуанкаре. Она определяется с помощью классических декартовых координат и ньютоновских импульсов (произведение массы тела на его скорость) в фиксированной инерциальной (абсолютной) системе отсчета. Именно, если в системе, состоящей из Солнца и N планет, обозначить радиус–вектор Солнца, отнесенный к абсолютным осям, как X0 , а абсолютные радиусы–векторы планет как Xk (k = 1, . . . , N ), то за обобщенные координаты uk принимаются координаты всех планет, отнесенные к Солнцу: uk = Xk − X0 , k = 1, . . . , N. e k по построению берут N соответствующих ньюЗа сопряженные обобщенные импульсы u тоновских импульсов этих планет e k = Pk , u k = 1, . . . , N, отнесенных к абсолютным осям. Для доопределения системы координат (система содержит N + 1 тело) вводится еще одна пара координата–импульс, за которые по построению 1 В литературе также используются названия канонические относительные координаты [6], а также канонические астроцентрические координаты [13]. 4

полагают соответственно радиус–вектор Солнца и полный ньютоновский импульс материальной системы в абсолютной системе отсчета: u0 = X0 , e0 = u N X Pk . k=0 При указанном выборе обобщенных координат и импульсов уравнения движения системы могут быть записаны в гамильтоновой форме и положены в основу всего дальнейшего построения планетной теории. В настоящей работе предлагается метод исследования планетных орбит, в котором движение каждого тела рассматривается также в гелиоцентрической системе отсчета, но в основу которого положен лагранжев подход при составлении исходных уравнений. А именно, для описания материальной системы используются гелиоцентрические координаты и производные этих координат по времени. Фундаментальным этапом построения теории планетного движения (в любой системе координат) является разложение гамильтониана системы в ряд при помощи базовых функций. Чаще всего — это ряд Пуассона. Выбор метода разложения определяется особенностями и условиями рассматриваемой задачей. Наиболее трудоемким шагом (в любом методе) является разложение величины обратного расстояния между двумя планетами, входящей в главную часть возмущающей функции [8, 12]. Со времен Эйлера было предложено множество способов разложения возмущающей функции системы. Одним из наиболее известных способов представления главной части является разложение с помощью полиномов Лежандра [11]. Метод эффективен, если отношение расстояний возмущаемого и возмущающего тел до центрального является малым [17]. В другом классическом методе применяются коэффициенты Лапласа (частный случай функций Пуассона [21]). Как отмечается в работе [17], использующий коэффициенты Лапласа метод можно употреблять в большем числе случаев по сравнению с методом, опирающимся на полиномы Лежандра. В любом способе, использующем ряды по степеням эксцентриситетов и наклонов (или по степеням величин порядка эксцентриситета и наклона), последние должны быть малы. При достижении значений 0.2 ÷ 0.3 соответствующий ряд может расходиться [16, 20]. В настоящей работе для разложения главной части возмущающей функции мы используем эффективный метод, опирающийся на коэффициенты Лапласа [8]. Этот метод может быть применен для исследования эволюции орбит, близких к круговым и компланарным. 5

2 2.1 ОСНОВНАЯ ЧАСТЬ Гелиоцентрические координаты Следуя Уинтнеру [15], введем в употребление гелиоцентрическую систему координат. Будем придерживаться системы обозначений, принятой в работе [19]. В этой работе используется система координат Якоби. В одинаковой системе обозначений удобно сравнивать выражения одной и той же физической величины в разных системах отсчета. Пусть в трехмерном евклидовом пространстве R3 находятся материальные точки Q0 (“Солнце”) с массой m0 и Qk (“планеты”) с массами µmk m0 , k = 1, . . . , N . Величина m0 имеет размерность массы, mk и µ полагаются безразмерными. “Массовый множитель” mk имеет порядок единицы, µ — малый параметр. На практике за µ обычно берут отношение массы самой массивной планеты к массе центральной звезды или отношение суммы масс всех планет к массе звезды. В последнем случае mk < 1 для всех k. При исследовании Солнечной системы обычно принимают µ = 10−3 . Такой выбор обозначений планетных масс связан с тем, что в большинстве известных планетных систем массы планет существенно меньше массы центральной звезды и в то же время сравнимы между собой. Q0 Q0 r1 r1 Q1 r2 Q1 r2 r3 Q2 Q2 r3 r0 r0 Q3 Q3 O O Гелиоцентрические (слева) и якобиевы (справа) координаты для одной и той же конфигурации, состоящей из четырех тел Введем декартову инерциальную систему координат с началом в точке O. Положение тел Q0 и Qk определим в ней соответственно радиусами–векторами ρ0 и ρk . Векторы 6

r1 , r2 , . . . , rN , определяющие соответственно положения “планет” Q1 , Q2 , . . . , QN относительно “Солнца” Q0 , называются по определению гелиоцентрическими координатами. Если мы задаем систему гелиоцентрических координат таким образом, то N + 1 вектор ρ0 , ρ1 , . . . , ρN заменяется лишь на N векторов r1 , . . . , rN . Для невырожденности преобразования, следуя монографии [15], добавим к координатам планет относительно Солнца абсолютные координаты центра масс системы2 . В результате система гелиоцентрических координат определяется соотношениями N µ X 1 mj ρj , r0 = ρ0 + m e m e j=1 rk = ρk − ρ0 , (1) k = 1, . . . , N, где m e = 1 + µm1 + . . . + µmN , а r0 — вектор, определяющий положение центра масс системы относительно точки O. Из определения используемых в работе [19] координат Якоби (см. формулы (П1) в Приложении А, а также общее определение координат Якоби в руководствах [5, 14]) и определения гелиоцентрических координат (1) следует, что векторы r0 и r1 в обех системах имеют одинаковый смысл (см. рисунок). Преобразование (1) является линейным с постоянными зависящими лишь от масс коэффициентами. Определитель DH системы (1) не изменится, если к первому столбцу прибавить сумму всех остальных столбцов:

1 µm1 µm2 µmN −1 µmN

µmN −1 µmN

µm1 µm2

... ...

m m e m e m e m e

m e m e m e m e

−1

1 0 . . . 0 0 1 0 . . . 0 0

−1

0 1 . . . 0 0 0 1 . . . 0 0

. DH =

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 0 ... 1 0

0 0 0 ... 1 0

−1

0 0 ... 0 1 0 0 0 ... 0 1

Очевидно, DH = 1 и преобразование (1) неособое. 2 В работе [19] при переходе от абсолютных координат к якобиевым возникает такое же затруднение: N + 1 вектор заменяется лишь на N векторов. Оно также разрешается путем учета радиуса–вектора барицентра системы в абсолютной системе отсчета. 7

Замечание. В работе [19] рассматривается аналогичное преобразование абсолютных координат ρ0 , ρ1 , . . . , ρN к якобиевым и утверждается, что его определитель DJ также равен единице. Доказать это утверждение можно также прибавлением в определителе к первому столбцу всех остальных столбцов. Доказательство мы приводим в Приложении А. Обратное к (1) преобразование выглядит следующим образом: N ρ 0 = r0 − µ X mj rj , m e j=1 (2) N µ X mj rj , ρk = rk + r0 − m e j=1 k = 1, . . . , N. Элементарный вывод соотношений (2) из (1) приведен в [15]. Важно, что разности гелиоцентрических (как и якобиевых [19]) координат не зависят от r0 . Более того, разности абсолютных и гелиоцентрических координат для произвольной пары планет совпадают: 2.2 ρk − ρ0 = rk , 1 6 k 6 N, ρj − ρk = rj − rk , 1 6 k < j 6 N. (3) Кинетическая и потенциальная энергия Для записи уравнений Лагранжа и функции Гамильтона необходимо определить кинетическую и потенциальную энергию системы в гелиоцентрических координатах. Будем использовать следующие удобные обозначения знаков сумм: X 1 = N X X , = 2 k=1 X . 16j<k6N Считаем, что пробегаемые индексы стоящих внутри этих сумм величин могут обозначаться различными буквами. Для кинетической энергии E1 по определению имеем 2E1 = m0 ρ̇20 + µm0 X mk ρ̇2k . 1 Подстановка (2) после преобразований приводит к представлению 2E1 = mm e 0 ṙ20 + µm0 X mk ṙ2k − 1 8 2 µ2 m0 X mk ṙk . m e 1 (4)

Квадрат суммы можно раскрыть следующим образом: X mk ṙk 2 = 1 X m2k ṙ2k + 2 1 X mj mk ṙj ṙk . (5) 2 Приведем для сравнения выражение кинетической энергии в координатах Якоби [19]: 2E1 = mm e 0 ṙ20 + µm0 X mk m e k−1 1 m ek ṙ2k , где m e 0 = 1, m e k = 1 + µm1 + . . . + µmk . Как и в координатах Якоби, в гелиоцентрической системе отсчета ṙ0 входит в E1 только через слагаемое, пропорциональное ṙ20 . Однако остальные скорости входят не только квадратами, но и попарными произведениями ṙj ṙk , что является недостатком гелиоцентрических координат. Потенциальную энергию E2 удобно представить в виде суммы двух частей. Первая определяется взаимодействием планет с Солнцем, а вторая — взаимодействием планет друг с другом: − X X mj mk mk E2 2 = µ + µ , Gm20 |ρk − ρ0 | |ρj − ρk | 1 2 где G — постоянная тяготения. Если массы планет малы по сравнению с массой Солнца (а мы только такие системы и рассматриваем), то второе слагаемое в этом выражении мало́ по сравнению с первым. Согласно (3), переход от координат ρ к r дает − X mk X mj mk E2 2 = µ + µ , Gm20 rk |rj − rk | 1 2 rk = |rk | . (6) Представление потенциальной энергии в координатах Якоби [19] имеет более сложный вид по сравнению с (6): разность радиусов–векторов ρj −ρk в общем случае зависит (линейно) также от положения промежуточных тел (см. формулы (3) в работе [19]). Например, в выражение ρ5 − ρ2 , помимо r2 и r5 , входят векторы r3 и r4 : ρ5 − ρ2 = r5 + µ m4 m3 m e1 r4 + µ r3 − r2 . m e4 m e3 m e2 Итак, можно сделать вывод, что кинетическая энергия имеет более простой и удобный вид в якобиевой системе координат, а потенциальная — в гелиоцентрической. 9

2.3 Уравнения Лагранжа и функция Гамильтона Действующие на планеты силы обусловлены стационарным потенциальным полем. Кинетическая энергия E1 представляется однородной квадратичной формой скоростей (4) с постоянными (зависящими лишь от масс) коэффициентами. Согласно [2], функции Лагранжа и Гамильтона в этом случае определяются соответственно равенствами L = E1 − E2 и H = E1 + E2 , а классические уравнения Лагранжа второго рода имеют вид ∂L d ∂L − = 0, s = 1, . . . , 3N + 3. dt ∂ q̇s ∂qs (7) В этой системе второго порядка с 3N + 3 степенями свободы за обобщенные координаты qs принимаются координаты всех N + 1 гелиоцентрических векторов r0 , r1 , . . . , rN . Компоненты вектора r0 являются циклическими координатами. Кинетическая энергия E1 не зависит от координат, а потенциальная энергия E2 — от скоростей. Поэтому уравнения Лагранжа (7) расщепляются на систему d ∂E1 ∂E2 , s = 1, . . . , 3N, =− dt ∂ q̇s ∂qs (8) записанную для координат векторов r1 , . . . , rN , и на одно векторное уравнение r̈0 = 0. Последнее равенство означает равномерное и прямолинейное движение центра масс системы в инерциальном пространстве. Поведение планет описывается замкнутой системой (8), причем первое слагаемое справа в (4) можно опустить как не влияющее на уравнения движения. Далее, соотношения (8) в силу их однородности можно сократить на постоянный множитель µm0 , что равносильно замене E1 и E2 соответственно функциями E3 и E4 : 2 X µ X mk ṙ2k − 2E3 = mk ṙk , m e 1 1 X mk X mj mk E4 − = +µ . Gm0 rk |rj − rk | 1 2 (9) При этом гамильтониан системы принимает вид h = E3 + E4 . Представим функцию Гамильтона h стандартной суммой невозмущенного и возмущающего гамильтониана. Как показывают соотношения (9), такое представление можно осуществить двумя способами. Во-первых, можно положить h=e h0 + µe h1 , 10 (10)

где невозмущенная часть e h0 и возмущающая µe h1 определяются равенствами X ṙ2k κ 2 e h0 = mk − , κ 2 = Gm0 , κ > 0 , 2 r k 1 X X mj mk 2 1 e mk ṙk − κ 2 . h1 = − 2m e 1 |rj − rk | 2 (11) Во-вторых, можно раскрыть с помощью (5) квадрат суммы в (9): X µmk 2µ X 2E3 = 1− mk ṙ2k − mj mk ṙj ṙk , m e m e 2 1 после чего представить h в виде h = h0 + µh1 , где невозмущенный гамильтониан h0 и возмущающий µh1 имеют вид X ṙ2 κ 2 k − k , h0 = βk 2 rk 1 κ2 X ṙj ṙk mj mk h1 = − + . |rj − rk | m e 2 (12) (13) В последних равенствах положено µmk βk = mk 1 − , m e m e µmk −1 = κ2 , κk2 = κ 2 1 − m e m e − µmk (14) причем величину µ в (14) следует считать данным числом, а не параметром. В уравнениях (11) гравитационный параметр κ 2 не зависит от массы k-ого тела. В уравнениях же (13) параметр κk2 для каждой планеты имеет свое определенное значение. Оба представления гамильтониана (10) и (12) можно использовать на практике, однако далее мы будем работать с более удобным с физической точки зрения представлением (12). 2.4 Системы оскулирующих элементов Если массами планет можно пренебречь по сравнению с массой Солнца, а взаимные расстояния между планетами имеют положительную нижнюю границу, то в (12) можно положить µ = 0. Задача планетного гелиоцентрического движения сводится к задаче о движении каждой планеты в статическом гравитационном поле Солнца: r̈s + κs2 rs = 0, rs3 s = 1, . . . , N. 11 (15)

Уравнения (15) соответствуют N независимым задачам одного притягивающего центра, каждая со своим гравитационным параметром κs2 . Решение каждой задачи дает шесть постоянных интегрирования, в качестве которых обычно выбирают кеплеровы элементы соответствующей орбиты или различные системы канонических элементов. Для большой полуоси, эксцентриситета и наклона орбиты будем использовать стандартные обозначения a, e, i. Среднюю аномалию, аргумент перицентра и долготу восходящего узла обозначим через M, g, Ω. Для планетных задач характерна малость эксцентриситетов и наклонов. В этом случае подходящей канонической системой является вторая система элементов Пуанкаре [5, 14] √ L = βκ a , p p √ ξ1 = 2L 1 − 1 − e2 cos ge , √ √ √ 4 ξ2 = 2L 1 − e2 1 − cos i cos Ω , λ = M + ge , p p √ η1 = − 2L 1 − 1 − e2 sin ge , √ √ √ 4 η2 = − 2L 1 − e2 1 − cos i sin Ω , (16) где через ge = g + Ω обозначена долгота перицентра. Величины ξ1 и η1 называются эксцентрическими элементами, а ξ2 и η2 — облическими (связанными с наклоном) элементами. Вместо элементов (16) далее будем использовать систему модифицированных (комплексных) элементов Пуанкаре [8, 10], определяемую равенствами3 в которых L0 = L , λ0 = λ , ξ10 = X , η10 = −iX̄ , ξ20 = Y , η20 = −iȲ , (17) p √ 1 − 1 − e2 eieg , √ √ √ 4 Y = L0 1 − e2 1 − cos i eiΩ , X = p L0 e — неперово число и i — мнимая единица. Элементы (17) можно назвать комплексной формой канонических элементов Пуанкаре (16). Однако их следует отличать от используемой в работе [6] системы комплексных переменных Пуанкаре — Харцера 3 L00 = L , λ00 = iλ , ξ100 = X , η100 = X̄ , ξ200 = Y , η200 = Ȳ , Чертой сверху мы везде обозначаем комплексное сопряжение. 12

в которой за сопряженный к L элемент вместо λ принимается iλ. Замечание. В любой системе элементов оскулирующие элементы каждой s-й планеты (включая величины β и κ) снабжаются соответствующим индексом s. Согласно работам [8, 10], основные разложения кеплеровского движения удобно строить по степеням безразмерных величин r 2 , X = ξ10 L0 ξ0 Y =√2 0, 2L r X̄ = iη10 2 , L0 iη 0 Ȳ = √ 2 0 . 2L (18) Пара канонически сопряженных переменных (L0 , λ0 ) тождественно совпадает с парой (L, λ), поэтому в системе элементов (17) штрихи над L и λ далее опустим. Система (17) имеет определенные преимущества перед классическими элементами (16). А именно, элементы (17) позволяют получить более экономные (по числу членов ряда) основные разложения кеплеровского движения при заданной точности разложения. Подробное изложение этого вопроса приведено в работе [10]. 2.5 Разложение гамильтониана Перейдем к разложению гамильтониана h планетной системы в ряд Пуассона по комплексным элементам Пуанкаре (17). Согласно (13), невозмущенная часть гамильтониана h0 элементарно выражается через оскулирующие элементы: h0 = − X βs κ 2 s 1 2as =− X β 3κ4 s 1 s 2L2s . Величина h1 является сложной функцией всех элементов, имеющей физическую размерность квадрата скорости. Следуя [19], представим ее в виде h1 = Gm0 h2 , a0 где множитель Gm0 /a0 имеет размерность квадрата скорости, а h2 безразмерно. За a0 можно принять любую характерную для планетной системы имеющую размерность длины постоянную. Для Солнечной системы a0 можно положить равной астрономической единице длины. Безразмерная возмущающая функция системы h2 принимает вид a X a0 0 h2 = − mj mk + ṙj ṙk . |rj − rk | Gm0 m e 2 13

Представим h2 суммой h3 + h4 , где X mj mk , |r − r | j k 2 a0 X h4 = − mj mk ṙj ṙk . Gm0 m e 2 h3 = −a0 (19) Величины h3 и h4 называются соответственно главной и дополнительной частью возмущающей функции. Главная часть обусловлена прямым взаимодействием планет, а дополнительная определяется выбором системы координат. Разложение главной части представляет основную сложность при разложении h2 в ряд Пуассона. Разложение h4 выполнить намного проще и оно в своем представлении (при одной и той же точности разложения) содержит значительно меньше слагаемых. 2.5.1 Разложение дополнительной части возмущающей функции Разложение дополнительной части h4 сводится к разложению скалярного произведения гелиоцентрических скоростей: ṙj ṙk = ẋj ẋk + ẏj ẏk + żj żk . Один из возможных методов разложения функций ẋ, ẏ, ż по комплексным элементам Пуанкаре (17) приведен в [10]. В работе [8] для этой цели используются промежуточные разложения по модифицированным (комплексным) элементам Лагранжа. 2.5.2 Разложение главной части возмущающей функции Разложение главной части h3 представляет собой трудоемкую операцию, так как обратное расстояние |rj − rk |−1 имеет наиболее сложную структуру в буквенных теориях классического типа. Представим h3 в виде h3 = X Rjk , Rjk = −a0 mj mk /∆jk , ∆jk = |rj − rk | . (20) 2 Число N функций Rjk равно числу всевозможных планетных пар в системе, то есть N (N − 1)/2. Для двупланетной системы N = 1, для восьмипланетной имеем N = 28. Следуя работе Ласкара и Робютеля [8], приведем метод разложения функции Rjk в ряд Пуассона по элементам (17). Величина Rjk определяется положением лишь двух планет, поэтому условимся обозначать внутреннюю и внешнюю планеты соответственно как Q 14

и Q0 и вместо Rjk писать R. Относящиеся к внешней планете величины будем снабжать штрихами. Запишем функцию R в виде R=− a0 mm0 a0 , a0 ∆ (21) где ∆ = |r − r0 |. Обозначим гелиоцентрический угол между планетами Q и Q0 через ϕ. По теореме косинусов ∆2 = r2 + r02 − 2rr0 cos ϕ, где r = |r|. Введем обозначения σ = r/r0 , α = a/a0 . За λ и λ0 примем соответственно средние долготы планет Q и Q0 в орбите. Легко убедиться, что a0 a0 a0 = 0 (1 + σ 2 − 2σ cos ϕ)−1/2 = 0 (A + P )−1/2 , ∆ r r где A = 1 + α2 − 2α cos(λ − λ0 ) , P1 = cos(λ − λ0 ) − P2 = σ 2 α σ cos ϕ , α (22) − 1, P = 2αP1 + α2 P2 . Величина A зависит только от двух аргументов: α и λ − λ0 . Из первого равенства (22) следует, что A > (1 − α)2 при любых вещественных значениях λ − λ0 . Для всякой рассматриваемой нами планетной пары α < 1, так что A отделена от нуля и можно записать −1/2 a0 a0 −1/2 P = 0A . 1+ ∆ r A Если эксцентриситеты и наклоны орбит достаточно малы и если α не слишком близко к единице, справедливо P < A. В самом деле, согласно (22), величина P разлагается в ряд по степеням эксцентрических и облических (позиционных) элементов4 . Непосредственным получением первых членов можно убедиться, что P имеет первый порядок малости по эксцентрическим элементам (см. также формулы (47, 48) в работе [8]). В то же время при не слишком близком к единице значении α величина A имеет порядок единицы. 4 Далее для краткости под позиционными будем везде подразумевать совокупность эксцентрических и облических переменных. 15

Итак, P/A < 1. Поэтому величина (1 + P/A)−1/2 разлагается в ряд по степеням отношения P/A, что приводит к представлению a0 a0 1 a0 3 a0 2 −5/2 5 a0 3 −7/2 35 a0 4 −9/2 −3/2 = 0 A−1/2 − P A + P A − P A + P A − ... ∆ r 2 r0 8 r0 16 r0 128 r0 Положим Uk = a0 k P , r0 k = 0, 1, 2, . . . (23) Тогда 1 3 5 35 a0 = U0 A−1/2 − U1 A−3/2 + U2 A−5/2 − U3 A−7/2 + U4 A−9/2 − . . . (24) ∆ 2 8 16 128 Для всякого полуцелого положительного s функция A−s разлагается в ряд Фурье по косинусам углов, кратных λ−λ0 , который сходится для всех вещественных значений λ−λ0 . Коэффициенты Uk являются степенными рядами относительно эксцентрических и облических элементов. В результате получаем ряд Пуассона. Покажем, каким образом строится разложение в ряд функций A−s и Uk . Величина A−s является производящей функцией коэффициентов Лапласа. Ее представление в виде ряда Лорана выглядит следующим образом: A−s = 1 X (k) ik(λ−λ0 ) b e . 2 k∈Z s (25) (k) Коэффициенты Лапласа bs (α) являются функциями одного непрерывного аргумента α. Разложение величин Uk , согласно (22, 23), сводится к разложению основных функций кеплеровского движения, поскольку σ r a0 = , α a r0 cos ϕ = xx0 + yy 0 + zz 0 . rr0 Каждая из безразмерных функций r/a, a/r, x/r, y/r, z/r разлагается в степенной ряд относительно соответствующих планете позиционных элементов. Метод разложения этих функций по безразмерным комплексным параметрам (18) приведен в работе [10]. 2.5.3 Вычислительный алгоритм разложения функции a0 /∆ На практике работа с a0 /∆ означает использование отрезка ряда (24). Количество удерживаемых слагаемых определяется точностью, с которой решается задача. Ряд (24) является рядом Пуассона, поэтому разложение строится до определенной степени d позиционных и до определенной кратности w угловых элементов. 16

Функции A−s не зависят от позиционных элементов, поэтому в (24) до степени d разлагаются величины U (то есть U0 , U1 , U2 , . . .). Разложение всех U в степенной ряд строится с помощью (23). Функция P имеет первый порядок малости по эксцентрическим элементам, а a0 /r0 — нулевой, поэтому Uk является величиной k-го порядка малости относительно позиционных элементов. Отсюда следует, что величины Uk разлагаются лишь для k 6 d. Для вычисления Ud достаточно взять разложение P до первой степени. Количество элементов в ряде P больше по сравнению с числом элементов в обычных кеплеровских разложениях (P зависит от положения двух планет). При заданной точности разложения это приводит к увеличению длины используемого отрезка ряда. Численные эксперименты показывают, что число членов в отрезках рядов U возрастает при этом на 2-3 порядка. Например, разложения функций x/r и y/r до 12-го порядка малости по степеням параметров (18) содержат по 446 слагаемых [10]. В разложении функции U3 до 12-го порядка малости по тем же элементам содержится 256 401 слагаемое. Эти оценки размеров используемых частных сумм U означают, что замечание о порядках малости величин U может быть эффективно использовано на практике. Именно, для разложения рядов Uk (k = 0, 1, 2, . . . , d) до заданной степени d не обязательно каждый раз брать ряд P (в соответствии с (23)), разложенный до степени d (разложение P до степени d потребуется лишь при вычислении U1 ). Будем использовать метод, основанный на следующем свойстве однородных многочленов (форм): их произведение есть также форма [7]. Этот метод использует лишь операции сложения и умножения рядов. Прежде всего, из (23) вытекает рекуррентное соотношение для коэффициентов U : Uk = Uk−1 P , k = 1, 2, . . . , d. Отсюда следует, что разложение рядов U можно строить по очереди: сначала U0 = a0 /r0 , затем U1 , U2 , . . . , Ud путем домножения каждый раз на P . Сначала ряд P разлагаем до степени d (по позиционным элементам). Полученный многочлен Pd представляем суммой однородных многочленов pj степени j относительно позиционных элементов: Pd = p1 + p2 + . . . + pd . Построение рядов U производится также путем получения и сложения однородных многочленов относительно эксцентрических и облических переменных: Uk = d X j=k 17 ukj .

Каждая форма ukj имеет двойной индекс. Первый означает индекс величины Uk , а второй j = k, k + 1, . . . , d равен степени формы ukj . После получения Pd разложение U0 = a0 /r0 до степени d представляем суммой U0 = u00 + u01 + . . . + u0d . Свободный член в разложении a0 /r0 равен 1, так что u00 = 1. После этого шага получаем последовательно остальные величины U1 , U2 , . . . , Ud . Ряд U1 складывается из форм u1j , где 1 6 j 6 d. Каждая форма u1j складывается из произведений некоторых форм, входящих в состав U0 , на соответствующую форму из многочлена Pd . Из U0 и Pd берутся все такие формы, степени которых в сумме дают j: u1j = u00 pj + u01 pj−1 + u02 pj−2 + . . . + u0,j−1 p1 = j−1 X u0s pj−s . s=0 Поясним сказанное примерами: u11 = u00 p1 , u12 = u00 p2 + u01 p1 , u14 = u00 p4 + u01 p3 + u02 p2 + u03 p1 , u1d = u00 pd + u01 pd−1 + u02 pd−2 + . . . + u0,d−1 p1 . Получение остальных Uk (k = 2, 3, . . . , d) совершается с помощью форм, входящих в состав предыдущего Uk−1 , и форм из состава Pd . Дадим несколько примеров: u22 = u11 p1 , u48 = u33 p5 + u34 p4 + u35 p3 + u36 p2 + u37 p1 , u9,11 = u88 p3 + u89 p2 + u8,10 p1 , 12 u12,12 = u11,11 p1 = u10,10 p1 p1 = u99 p31 = . . . = u00 p12 1 = p1 . Общие формулы получения ukj и Uk для k = 1, 2, . . . , d имеют вид ukj = j−1 X uk−1,s pj−s , Uk = s=k−1 j−1 d X X uk−1,s pj−s . (26) j=k s=k−1 На практике — при работе с системами компьютерной алгебры — такой метод построения рядов U оказывается эффективным, так как размеры используемых отрезков рядов U 18

являются относительно большими и удобнее работать с их более мелкими частями, которыми являются рассмотренные формы ukj . Рассмотрим теперь, каким образом останавливается разложение a0 /∆ на заданной кратности w по λ и λ0 . Средние долготы λ и λ0 входят в слагаемые полученных отрезков 0 0 рядов U тригонометрическими многочленами Λv Λ0v , где Λ = eiλ , Λ0 = eiλ , а v и v 0 — целые числа. Поэтому максимальная кратность w разложения a0 /∆ относительно λ и λ0 определяется кратностью разложения функций A−s . Необходимо узнать, какие пределы следует выбрать для k в (25), чтобы после перемножения рядов U и A−s у всех полученных членов разложения a0 /∆ степени элементов Λ и Λ0 не превышали по модулю w. Поскольку 0 0 Λv Λ0v Λk Λ0−k = Λv+k Λ0v −k , задача сводится к нахождению всех целых k, удовлетворяющих системе   −w 6 v + k 6 w ,  −w 6 v 0 − k 6 w . (27) 0 Из (27) заключаем, что при заданном w член из рядов U , имеющий множителем Λv Λ0v , домножается лишь на те члены ряда A−s = 1 X (k) k 0−k 1 −1 0 (0) (1) 0−1 b Λ Λ , . . . + b(−1) Λ Λ + b + b ΛΛ + . . . = s s s 2 2 k∈Z s (28) у которых max{−v − w, v 0 − w} 6 k 6 min{w − v, w + v 0 }. При изучении долговременных изменений планетных орбит может возникнуть необходимость извлечения из a0 /∆ лишь тех слагаемых, угловые аргументы которых являются заданными линейными комбинациями средних долгот nλ + n0 λ0 , где n и n0 — целые числа (например, 2λ − 5λ0 ). Согласно работе [8], для получения слагаемых из a0 /∆, ассоциированных с угловым аргументом nλ + n0 λ0 , удобно использовать определение долготной характеристики n + n0 члена разложения. Долготная характеристика каждого слагаемого ряда (28) равна нулю, поэтому для получения членов ряда a0 /∆ с аргументом nλ + n0 λ0 достаточно из разложения величин U взять только те слагаемые, у которых степени элементов Λ и Λ0 в сумме дают n + n0 (с последующим домножением этих слагаемых на соответствующий член из ряда (28)). 19

2.6 Программное обеспечение для построения разложений Практическое разложение главной h3 и дополнительной h4 части возмущающей функции выполняется с помощью специальных алгебраических пакетов. Нами были написаны вычислительные программы, позволяющие получать необходимые разложения до заданной степени точности. Все разложения строились в системе компьютерной алгебры Piranha [3]. На некоторых этапах работы для контроля правильности получаемых результатов мы применяли систему общего назначения Maple. В системе Piranha поддерживаются различные форматы представления чисел. С целью избежать накопления ошибок, коэффициенты всех получаемых разложений представлялись в виде рациональных чисел. Пользовательский интерфейс Piranha реализован на языке Python, поэтому библиотеку всех написанных программ мы составили в виде нескольких модулей (файлов с расширением .py), расположенных в корневой директории системы. 2.6.1 Модуль M_KEPLER.py Этот модуль содержит процедуры5 , реализующие разложения основных формул кеплеровского движения. Алгоритм получения необходимых разложений r/a, a/r, x/r, y/r, z/r, а также ẋ, ẏ, ż по параметрам (18) подробно описан в работе [10]. Следуя [10], в основу кеплеровского процессора мы положили функции z1 = e sin M и z2 = e cos M . Аналитическое представление z1 и z2 через параметры (18) имеет следующий вид (см. формулы (13) в работе [10]): r i ih −1 z1 = XΛ − X̄Λ 1 − 2 r i 1h −1 z2 = XΛ + X̄Λ 1 − 2 X X̄ , 4 (29) X X̄ . 4 Разложение правых частей (29) в степенной ряд выполняется системой Piranha автоматически. Единственным аргументом dg строящей эти разложения процедуры z1z2(dg) является степень, до которой требуется получить разложение. Полученный отрезок ряда записывается системой в отдельный файл и сохраняется. В приложении Б приведена 5 Напомним, что процедурой в языке Python называют заданную последовательность команд, которая в точку вызова не возвращает никакого значения. Последовательность команд, возвращающая в точку вызова какое-либо значение, называется функцией. 20

машинная запись разложения z1 до 12 степени включительно в файл z1_12.qpc. Приведем первые члены разложения функции z1 : i i i i z1 = XΛ−1 − X̄Λ − X 2 X̄Λ−1 + X X̄ 2 Λ− 2 2 16 16 i i i i − X 3 X̄ 2 Λ−1 + X 2 X̄ 3 Λ − X 4 X̄ 3 Λ−1 + X 3 X̄ 4 Λ + . . . 256 256 2048 2048 После получения отрезков рядов z1 и z2 строятся последовательно все остальные необходимые разложения. Программная реализация разложения каждой функции в данном модуле по существу одна и та же: дается аналитическое выражение функции, указывается степень разложения, указываются переменные, по которым выполняется разложение, полученный многочлен записывается в файл с определенным именем в указанную директорию. В качестве примера в Приложении В приводится программа (процедура), реализующая разложение функции z3 = E −M по степеням z1 и z2 (через E обозначена эксцентрическая аномалия). Разложение можно выполнять с помощью формулы Бурмана–Лагранжа [21], однако мы предпочли использовать для этой цели не требующий операций дифференцирования метод итераций [21]. 2.6.2 Модуль M_POISSON.py В этом модуле находятся процедуры, выполняющие разложения пуассоновского процессора σ/α, cos ϕ, P , а также процедуры, строящие в соответствии с (23, 24) разложения рядов U и a0 /∆. Предложенный способ разложения коэффициентов U кажется громоздким, однако алгоритм работы с формулами (26) программируется без затруднений. Кроме того, разложение величин U происходит лишь один раз, поскольку полученные ряды записываются и сохраняются в соответствующих файлах. Метод получения отрезка ряда a0 /∆ при заданных d и w и метод получения членов a0 /∆, ассоциированных с заданным угловым аргументом nλ + n0 λ0 , программируются с помощью встроенных функций системы Piranha и стандартных средств языка Python. 2.6.3 Модуль M_LAPLACE.py Для сравнения полученного нами разложения a0 /∆ с приведенным в работе Ласкара и Робютеля [8] разложением a0 /∆ (см. следующие два параграфа) был написан мо21

дуль M_LAPLACE.py. В этом модуле находятся функции, позволяющие с помощью рекуррентых соотношений [1, 4] выполнять преобразования коэффициентов Лапласа между собой (возвращаемыми значениями этих функций являются коэффициенты Лапласа). В Приложении Г для примера мы приводим содержимое функции UP_IND(n,k) (смысл аргументов n и k приведен там же), получающей при помощи рекуррентного соотношения = b(k+1) s k k + s − 1 (k−1) (α + α−1 )b(k) b s − k−s+1 k−s+1 s (k) (0) (1) выражение произвольного коэффициента Лапласа bs через bs и bs (с коэффициентами, зависящими от целых степеней α). Из этого соотношения мы выводим следующие нужные нам представления (см. формулы (30, 31) и пояснения к ним): (2) (1) (0) b3/2 = 2(α + α−1 )b3/2 − 3b3/2 , (2) (0) (1) b5/2 = −2(α + α−1 )b5/2 + 5b5/2 , (0) (1) (3) b5/2 = −(8α−2 + 23 + 8α2 )b5/2 + 20(α + α−1 )b5/2 . Функция DOWN_IND(n,k,d) при помощи рекуррентных соотношений s s (k) (k+1) b(k) (1 + α2 )bs+1 − 2αbs+1 , s = k+s k+s s s (k) (k+1) 2αbs+1 − (1 + α2 )bs+1 , b(k+1) = s k−s+1 k−s+1 (k) bs+1 = (k+1) bs+1 = s + k 1 + α2 (k) k − s + 1 2α bs − b(k+1) , 2 2 s (1 − α ) s (1 − α2 )2 s s+k 2α k − s + 1 1 + α2 (k+1) (k) b b − , s (1 − α2 )2 s s (1 − α2 )2 s (k) а также при помощи вызова функции UP_IND(n,k) выражает коэффициент Лапласа bs (0) (1) для четного d (нуль тоже считаем четным) через bd/2+1/2 и bd/2+1/2 , а для нечетного — через6 (0) (1) bbd/2c+3/2 и bbd/2c+3/2 . Неотрицательное целое d (аргумент d) обозначает степень одночлена (k) по позиционным элементам, имеющего множителем bs в разложении a0 /∆. Приведем для примера два результата применения этой функции: (0) b5/2 = (3) b5/2 = 6 1 + α2 (0) 2α (1) b3/2 + b3/2 , 2 2 2 2 (1 − α ) 3(1 − α ) 4α−1 − 6α + 4α3 (0) −8α−2 + 9 + 9α2 − 8α4 (1) b + b3/2 . 3/2 (1 − α2 )2 3(1 − α2 )2 Скобками b c мы обозначаем целую часть числа. 22

2.7 Результаты Использование коэффициентов Лапласа означает представление a0 /∆ в виде одиннадцатиаргументного разложения. Участвуют 8 позиционных элементов X, X̄, Y, Ȳ , X 0 , X̄ 0 , Y 0 , Ȳ 0 , два угловых элемента Λ, Λ0 (при работе с вычислительными системами средние долготы λ, λ0 удобнее представлять элементами Λ, Λ0 ) и единственный аргумент коэффициентов Лапласа α. Разложение безразмерного a0 /∆ происходит по степеням одиннадцати указанных безразмерных элементов. Приведем первые члены полученного разложения a0 /∆. При построении всех разложений коэффициенты Лапласа сохранялись в символьном виде. Таблица 1. Разложение вековой части a0 /∆ по позиционным элементам до второй степени включительно Функция от α X X̄ X0 X̄ 0 Y Ȳ Y0 Ȳ 0 c0 (α) 0 0 0 0 0 0 0 0 c1 (α) 1 1 0 0 0 0 0 0 c1 (α) 0 0 1 1 0 0 0 0 c2 (α) 1 0 0 1 0 0 0 0 c2 (α) 0 1 1 0 0 0 0 0 c3 (α) 0 0 0 0 1 1 0 0 c3 (α) 0 0 0 0 0 0 1 1 c4 (α) 0 0 0 0 1 0 0 1 c4 (α) 0 0 0 0 0 1 1 0 В табл. 1 представлено разложение вековой части a0 /∆ до второй степени включительно по эксцентрическим и облическим элементам. В этом отрезке ряда содержится 31 слагаемое. Для большей наглядности и сокращения длины таблицы в первую колонку мы поместили функции от одного аргумента α, которые домножаются на соответствующий одно23

член по позиционным элементам. Степени позиционных элементов в этом одночлене приводятся в остальных колонках. Указанные функции от α имеют следующий вид: 1 (0) c0 (α) = b1/2 , 2 3 3 2 (0) 1 (1) 15 2 3 4 (0) 9 (1) (2) c1 (α) = − α b3/2 − αb3/2 + α + α b5/2 − α3 b5/2 − α2 b5/2 , 8 4 16 8 4 16 3 1 (2) 3 3 (0) 21 2 3 4 (1) 9 3 2 (1) (2) (3) α + α b5/2 + α3 b5/2 + α2 b5/2 , c2 (α) = α b3/2 + αb3/2 + α b5/2 − 8 4 8 32 8 8 32 (30) 1 (1) −c3 (α) = c4 (α) = αb3/2 . 2 Как видим из табл. 1, функция c0 (α) является свободным членом разложения a0 /∆ как ряда Пуассона (24). Можно показать, что разложение a0 /∆ в виде ряда (24) обладает даламберовской характеристикой [18], поэтому в разложении вековой части a0 /∆ содержатся одночлены только четной степени (по эксцентрическим и облическим элементам). Таблица 2. Члены разложения a0 /∆ до второй степени включительно по позиционным элементам, ассоциированные с угловыми аргументами λ − 2λ0 и −λ + 2λ0 X X0 X̄ X̄ 0 α b1/2 (1) b3/2 (0) b3/2 (1) b3/2 (2) b3/2 (3) Λ Λ0 1/4 1 0 0 0 2 0 0 0 1 0 1 −2 −3/8 1 0 0 0 1 0 0 1 0 0 1 −2 1/8 1 0 0 0 1 0 0 0 0 1 1 −2 −1/4 0 1 0 0 2 0 0 1 0 0 1 −2 3/8 0 1 0 0 1 0 1 0 0 0 1 −2 −1/8 0 1 0 0 1 0 0 0 1 0 1 −2 1/4 0 1 0 0 0 1 0 0 0 0 1 −2 1/4 0 0 1 0 2 0 0 0 1 0 −1 2 −3/8 0 0 1 0 1 0 0 1 0 0 −1 2 1/8 0 0 1 0 1 0 0 0 0 1 −1 2 −1/4 0 0 0 1 2 0 0 1 0 0 −1 2 3/8 0 0 0 1 1 0 1 0 0 0 −1 2 −1/8 0 0 0 1 1 0 0 0 1 0 −1 2 1/4 0 0 0 1 0 1 0 0 0 0 −1 2 24

В табл. 2 приведено разложение по эксцентрическим и облическим переменным долготной (периодической) части a0 /∆ до второй степени включительно. В этой таблице мы приводим члены, ассоциированные с аргументами λ − 2λ0 и −λ + 2λ0 . Итак, в таблицах 1 и 2 представлены некоторые первые члены разложения функции a0 /∆, входящей множителем в каждое слагаемое R главной части h3 возмущающей функции h2 (см. формулы (20, 21)). В Приложении Д приводятся (аналогично в виде таблицы) первые члены разложения дополнительной части h4 возмущающей функции h2 . В таблице дается разложение безразмерного множителя W 0 , входящего в скалярное произведение гелиоцентрических скоростей ṙṙ0 = ωω 0 aa0 W 0 (через ω обозначено среднее движение планеты). Согласно (19), величина ṙṙ0 входит в каждое слагаемое дополнительной части h4 . Согласно [10], для имеющего размерность квадрата скорости множителя ωω 0 aa0 справедливо равенство ωω 0 aa0 = κ 2 κ 02 ββ 0 , LL0 c помощью которого мы получаем окончательное разложение дополнительной части h4 по комплексным элементам (17). 2.8 Обсуждение результатов В основу всех разложений мы положили метод, разработанный в [10]. Этот метод опирается на изложенные в монографии Субботина [14] идеи. Представляется целесообразным сравнение с полученным в работе [8] разложением a0 /∆, поскольку мы используем одинаковые элементы: параметры (18) и коэффициенты Лапласа. Ласкар и Робютель [8] в разложении a0 /∆ используют преобразование коэффициентов Лапласа, указанное в пункте 2.6.3. Его применение существенно сокращает количество членов (см. подробности в работе [8]). Проиллюстрируем выигрыш в числе слагаемых на примере представленного таблицей 1 разложения вековой части a0 /∆. После необходимых преобразований в (30) функции c1 (α), c2 (α) приобретают вид 1 (1) c1 (α) = αb3/2 , 8 3 (0) 1 1 2 (1) c2 (α) = αb3/2 − + α b3/2 , 8 4 4 25 (31)

а функции c0 (α), c3 (α), c4 (α) не изменяются. Формулы (31) показывают, что наше разложение по позиционным элементам вековой части a0 /∆ до второй степени включительно совпадает с разложением a0 /∆, представленным в работе [8]. Представление (31) уменьшает число слагаемых в этом отрезке ряда с 31 (таблица 1) до 13 [8, стр. 215]. Согласно [8], при увеличении степени разложения a0 /∆ экономия слагаемых должна быть все более значительной. Например, в построенном нами разложении по позиционным элементам вековой части a0 /∆ до четвертой степени включительно содержится 577 слагаемых. Согласно [8], этот отрезок ряда сокращается до 165 слагаемых [8, стр. 215]. 3 ЗАКЛЮЧЕНИЕ Основанное на подходящей замене коэффициентов Лапласа сокращение числа слагаемых в разложении a0 /∆ вида (24) проиллюстрировано Ласкаром и Робютелем [8] лишь на примере разложения вековой части. Представляется целесообразным применить указанный метод также для периодической части разложения. Достоинством указанного метода является возможность избежать использования производных от коэффициентов Лапласа, а также то, что после его применения количество коэффициентов Лапласа с различными нижними и верхними индексами в разложении (24) значительно снижается. Результаты настоящей работы предполагается использовать для построения осредненных уравнений движения планетной системы. Автор выражает благодарность научному руководителю К. В. Холшевникову за руководство работой и помощь в подготовке рукописи, а также ресурсному центру «Вычислительного центра Санкт-Петербургского государственного университета» за предоставленную возможность использования вычислительного кластера при построении всех разложений. Работа выполнена при финансовой поддержке СПбГУ (грант 6.37.341.2015). 4 ПОЛОЖЕНИЯ, ВЫНОСИМЫЕ НА ЗАЩИТУ На защиту выносятся следующие положения: 1. Исследована гелиоцентрическая система координат. В отличие от многих современных работ, для вывода всех уравнений использован лагранжев формализм. Проведено элементарное сравнение гелиоцентрических координат с системой координат Якоби [19]. Выявлены сравнительные преимущества и недостатки указанных систем. 26

2. Предложена специальная система оскулирующих элементов. По сравнению с классическими системами, используемая в работе система комплексных канонических элементов Пуанкаре (17) позволяет получать более компактные разложения кеплерова движения, что оказывается выгодным при построении разложения возмущающей функции. 3. Построено разложение гамильтониана N -планетной задачи в ряд Пуассона. Предложен аглоритм разложения главной и дополнительной части возмущающей функции. Представлены первые члены разложения обеих частей. Проведено сравнение с результатами работы [8]. 4. Протестирована и изучена специальная система компьютерной алгебры [3]. Написан комплекс программ, реализующих все рассмотренные в работе разложения. Список литературы [1] Аксенов Е.П. Специальные функции в небесной механике (М.: Наука, 1986). [2] Арнольд В.И. Математические методы классической механики (М.: Наука, 1989). [3] Бискани (F. Biscani). The Piranha algebraic manipulator // arXiv:0907.2076v1. [4] Брауэр Д., Клеменс Дж. Методы небесной механики (М.: Мир, 1964). [5] Дубошин Г.Н. Небесная механика. Основные задачи и методы (М.: Наука, 1975). [6] Красинский Г.А. Основные уравнения планетной теории (М.: Наука, сб. Малые Планеты, под ред. Н.С. Самойловой-Яхонтовой, 1973), с. 81. [7] Курош А.Г. Курс высшей алгебры (М.: Наука, 1975). [8] Ласкар, Робютель (J. Laskar and P. Robutel). Stability of the Planetary Three-Body Problem. I. Expansion of the Planetary Hamiltonian // Celest. Mech. Dynam. Astron. 62, 193 (1995). [9] Либер, Сансотера (A.-S. Libert and M. Sansottera). On the extension of the LaplaceLagrange secular theory to order two in the masses for extrasolar systems // Celest. Mech. Dynam. Astron. 117, 149 (2013). 27

[10] Микрюков Д.В., Холшевников К.В. Разложение основных функций кеплеровского движения с использованием комплексных переменных // Письма в Астрон. журн. 42, 302 (2016). [11] Мюррей К., Дермотт С. Динамика Солнечной системы (М.: Физматлит, 2009). [12] Петровская М.С., Иванова Т.В. О построении разложений планетной возмущающей функции // Бюлл. ИТА АН СССР XIV, 288 (1978). [13] Родригес, Галлардо (A. Rodrı́guez and T. Gallardo). The dynamics of the HD12661 extrasolar planetary system // Astrophys. J. 628, 1006 (2005). [14] Субботин М.Ф. Введение в теоретическую астрономию (М.: Наука, 1968). [15] Уинтнер А. Аналитические основы небесной механики (М.: Наука, 1967). [16] Ферраз-Меллу (S. Ferraz-Mello). The convergence domain of the Laplacian expansion of the Disturbing Function // Celest. Mech. Dynam. Astron. 58, 37 (1994). [17] Ферраз-Меллу и др. (S. Ferraz-Mello, T.A. Michtchenko and C.Beaugé). Regular Motions in Extra-Solar Planetary Systems // arXiv:astro-ph/0402335v2. [18] Холшевников К.В. Даламберовские функции в небесной механике // Астрон. журн. 74, 146 (1997). [19] Холшевников К.В., Греб А.В., Кузнецов Э.Д. Разложение гамильтониана планетной задачи в ряд Пуассона по всем элементам (теория) // Астрон. вестник 35, 267 (2001). [20] Холшевников К.В., Греб А.В., Кузнецов Э.Д. Разложение гамильтониана планетной задачи в ряд Пуассона по всем элементам: оценка и прямое вычисление коэффициентов // Астрон. вестник 36, 75 (2002). [21] Холшевников К.В., Титов В.Б. Задача двух тел: Учебное пособие (СПб.: СПбГУ, 2007). 28

Приложение А Определитель преобразования абсолютных координат к якобиевым Пусть планетная система состоит из центральной звезды и N планет. Следуя работе [19], систему координат Якоби определим соотношениями N 1 µ X r0 = ρ0 + mj ρj , m e m e j=1 k−1 1 µ X rk = ρk − ρ0 − mj ρj , m e k−1 m e k−1 j=1 (П1) k = 1, . . . , N, где m e 0 = 1, m e k = 1 + µm1 + . . . + µmk , а остальные обозначения в правых частях имеют с (1) одинаковый смысл. Пустая сумма (для k = 1) считается равной нулю. Чтобы показать, что определитель DJ преобразования (П1) равен единице, прибавим к первому столбцу DJ все остальные столбцы:

µm1 µm2 µmN −1 µmN

1 ... ...

m e m e m e m e

m e m e m e m e m e

−1 1 0 ... 0 0

1 0 ... 0 0

−µm1

−1

−µm1

0 1 ... 0 0

1 ... 0 0

m e e m e 1 1 1

DJ =

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

−1 −µm −µm −µm −µm 1 2 1 2

... 1 0

m m e m e e N −2 m e N −2 N −2 N −2

e N −2 m

−1

−µm1 −µm2 −µmN −1 −µm1 −µm2 −µmN −1

... 1

0 ... 1

m e N −1 m e N −1 m e N −1 m e N −1 m e N −1 m e N −1 m e N −1 Дополнительный минор стоящей в левом верхнем углу единицы равен единице, поэтому разложение преобразованного определителя по первому столбцу дает DJ = 1. 29

Приложение Б Разложение функции z1 = e sin M до 12 степени включительно по параметрам X, X̄, Λ = eiλ в системе Piranha [poly_arg] name=Xc time_eval= [poly_arg] name=X time_eval= [poly_arg] name=LMB time_eval= [terms] (0,1/16)|2;1;1 (0,-1/16)|1;2;-1 (0,-1/2)|1;0;1 (0,-7/524288)|5;6;-1 (0,-1/2048)|3;4;-1 (0,-5/65536)|4;5;-1 (0,1/256)|3;2;1 (0,1/2048)|4;3;1 (0,1/2)|0;1;-1 (0,-1/256)|2;3;-1 (0,5/65536)|5;4;1 (0,7/524288)|6;5;1 30

Приложение В Процедура, выполняющая разложение функции z3 = E − M по степеням z1 = e sin M и z2 = e cos M методом итераций def z3_z1z2(dg): #аргумент dg - степень разложения z3 по z1 и z2 import pyranha Z1=pyranha.Qpoly.qpolyc(pyranha.Core.psym(’z1’)) Z2=pyranha.Qpoly.qpolyc(pyranha.Core.psym(’z2’)) W=pyranha.Qpoly.qpolyc(pyranha.Core.psym(’w’)) wz=Z1 SIN=0 COS=W**0 term=1 n=1 while n<dg: n=n+1 print "-"*80 print "n=",n term=term*W/(n-1) print "term = ",term if 1j**n==-1: # n=2,6,10,14,... SIN=SIN+term elif 1j**n==-1j: # n=3,7,11,15,... COS=COS-term elif 1j**n==1: # n=4,8,12,16,... SIN=SIN-term elif 1j**n==1j: # n=5,9,13,17,... COS=COS+term SINz=SIN.sub(’w’,wz) COSz=COS.sub(’w’,wz) pyranha.Truncators.truncators.degree.set(n+1) wz=Z1*COSz+Z2*SINz print "wz = ",wz pyranha.Truncators.truncators.degree.unset() wz.save_to(’z3_z1z2_’+str(dg)+’.qpc’) 31

Приложение Г (k) Функция, получающая выражение коэффициента Лапласа bs (0) (1) через bs и bs #Lp cf = Laplace coefficient #this function expreses b_n^k in terms of b_n^0 and b_n^1 #the lower index n is not transformed #this function is applied to all Lp cfs in the series #n=2*s, where s - half integer #k - the upper index def UP_IND(n,k): # аргумент n=2s - нижний индекс; аргумент k - верхний индекс import pyranha alp=pyranha.Qpoly.qpoly(pyranha.Core.psym(’ALP’)) #initial values of Lp cfs b_f=pyranha.Qpoly.qpoly(pyranha.Core.psym(’b’+str(n)+’^0’)) b_s=pyranha.Qpoly.qpoly(pyranha.Core.psym(’b’+str(n)+’^1’)) if k == 0: return b_f if k == 1: return b_s if k > 1: i=0 #i runs from 1 to k-1 #tne number of the loops is k-1 while i < k-1: i=i+1 b_nxt=2*i*(alp+alp**-1)*b_s/(2*i-n+2)-(2*i+n-2)*b_f/(2*i-n+2) b_f=b_s b_s=b_nxt return b_nxt 32

Приложение Д Разложение по эксцентрическим и облическим элементам функции W 0 до второй степени включительно 1/2 1/2 1/2 1/2 1/2 1/2 −1/16 9/16 −1/4 −1/4 1/2 9/16 −1/16 1/2 9/16 −1/16 −1/4 −1/4 −1/16 9/16 −1/2 −1/2 −1/2 1 1 −1/2 1 1 −1/2 −1/2 −1/2 −1/2 X 0 0 1 0 0 0 2 2 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X̄ 0 0 0 1 0 0 0 0 1 1 0 2 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 2 2 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X̄ 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 1 2 2 0 0 0 0 0 0 0 0 0 0 0 0 Y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 1 1 1 0 0 0 0 0 0 0 33 Ȳ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 2 1 1 0 0 0 0 Y0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 2 1 1 0 Ȳ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 1 2 Λ 1 −1 −2 2 1 −1 −1 −3 1 −1 −2 3 1 2 1 −1 1 −1 1 −1 −1 1 −1 −1 −1 1 1 1 −1 1 −1 1 Λ0 −1 1 1 −1 −2 2 −1 1 −1 1 2 −1 1 −2 −3 −1 −1 1 1 3 −1 −1 1 −1 1 1 −1 1 −1 −1 1 1

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв