Санкт-Петербургский государственный университет
Математико-механический факультет
Кафедра небесной механики
Микрюков Денис Викторович
Разложение гамильтониана планетной задачи в ряд Пуассона
в гелиоцентрической системе отсчета
Выпускная квалификационная работа
Научный руководитель:
д.ф.-м.н., заведующий Кафедрой небесной
механики, профессор Холшевников К.В.
Рецензент:
д.ф.-м.н., доцент УрФУ Кузнецов Э.Д.
Санкт-Петербург
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
0
.
.
.
0
0
1
0
.
.
.
0
0
0
1
.
.
.
0
0
0
1
.
.
.
0
0
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
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 ṙ20 + µm0
X
mk ṙ2k −
1
8
2
µ2 m0 X
mk ṙk .
m
e
1
(4)
Квадрат суммы можно раскрыть следующим образом:
X
mk ṙk
2
=
1
X
m2k ṙ2k + 2
1
X
mj mk ṙj ṙk .
(5)
2
Приведем для сравнения выражение кинетической энергии в координатах Якоби [19]:
2E1 = mm
e 0 ṙ20 + µm0
X mk m
e k−1
1
m
ek
ṙ2k ,
где m
e 0 = 1, m
e k = 1 + µm1 + . . . + µmk .
Как и в координатах Якоби, в гелиоцентрической системе отсчета ṙ0 входит в E1 только через слагаемое, пропорциональное ṙ20 . Однако остальные скорости входят не только
квадратами, но и попарными произведениями ṙj ṙ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 ṙ2k −
2E3 =
mk ṙ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
ṙ2k κ 2
e
h0 =
mk
−
, κ 2 = Gm0 , κ > 0 ,
2
r
k
1
X
X mj mk
2
1
e
mk ṙk − κ 2
.
h1 = −
2m
e 1
|rj − rk |
2
(11)
Во-вторых, можно раскрыть с помощью (5) квадрат суммы в (9):
X
µmk
2µ X
2E3 =
1−
mk ṙ2k −
mj mk ṙj ṙk ,
m
e
m
e 2
1
после чего представить h в виде
h = h0 + µh1 ,
где невозмущенный гамильтониан h0 и возмущающий µh1 имеют вид
X ṙ2 κ 2
k
− k ,
h0 =
βk
2
rk
1
κ2
X
ṙj ṙ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 = −iȲ ,
(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 = Ȳ ,
Чертой сверху мы везде обозначаем комплексное сопряжение.
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
Ȳ = √ 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
+
ṙj ṙk .
|rj − rk | Gm0 m
e
2
13
Представим h2 суммой h3 + h4 , где
X mj mk
,
|r
−
r
|
j
k
2
a0 X
h4 = −
mj mk ṙj ṙk .
Gm0 m
e 2
h3 = −a0
(19)
Величины h3 и h4 называются соответственно главной и дополнительной частью возмущающей функции. Главная часть обусловлена прямым взаимодействием планет, а дополнительная определяется выбором системы координат. Разложение главной части представляет основную сложность при разложении h2 в ряд Пуассона. Разложение h4 выполнить
намного проще и оно в своем представлении (при одной и той же точности разложения)
содержит значительно меньше слагаемых.
2.5.1
Разложение дополнительной части возмущающей функции
Разложение дополнительной части h4 сводится к разложению скалярного произведения
гелиоцентрических скоростей: ṙj ṙk = ẋj ẋk + ẏj ẏk + żj żk .
Один из возможных методов разложения функций ẋ, ẏ, ż по комплексным элементам
Пуанкаре (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,
а также ẋ, ẏ, ż по параметрам (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, Ȳ , X 0 , X̄ 0 , Y 0 , Ȳ 0 ,
два угловых элемента Λ, Λ0 (при работе с вычислительными системами средние долготы
λ, λ0 удобнее представлять элементами Λ, Λ0 ) и единственный аргумент коэффициентов
Лапласа α. Разложение безразмерного a0 /∆ происходит по степеням одиннадцати указанных безразмерных элементов.
Приведем первые члены полученного разложения a0 /∆. При построении всех разложений коэффициенты Лапласа сохранялись в символьном виде.
Таблица 1. Разложение вековой части a0 /∆ по позиционным элементам
до второй степени включительно
Функция от α
X
X̄
X0
X̄ 0
Y
Ȳ
Y0
Ȳ 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 , входящего в скалярное произведение гелиоцентрических
скоростей ṙṙ0 = ωω 0 aa0 W 0 (через ω обозначено среднее движение планеты). Согласно (19),
величина ṙṙ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.Beaugé). 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 все остальные столбцы:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
−1
−µm
−µm
−µm
−µm
1
2
1
2
m
m
e
m
e
e N −2 m
e N −2
N −2
N −2
−µm1 −µm2
−µmN −1
−µm1 −µm2
−µmN −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
Ȳ
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
Ȳ 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
Отзывы:
Авторизуйтесь, чтобы оставить отзыв