Программирование, 2022, № 2, стр. 53-62

ПОСТРОЕНИЕ ЭВОЛЮЦИОННЫХ УРАВНЕНИЙ В ЗАДАЧЕ МНОГИХ ТЕЛ С ИЗОТРОПНО ИЗМЕНЯЮЩИМИСЯ МАССАМИ С ПРИМЕНЕНИЕМ КОМПЬЮТЕРНОЙ АЛГЕБРЫ

А. Н. Прокопеня^a, *, М. Дж. Минглибаев^b, c, **, А. Б. Кошербаева^b, ***

^a Варшавский университет естественных наук – SGGW
02-776 Варшава, ул. Новоурсыновска, 159, Польша

^b Казахский национальный университет им. аль-Фараби
050040 Алматы, пр. аль-Фараби, 71, Казахстан

^c Астрофизический институт им. В.Г. Фесенкова
050020 Алматы, Обсерватория, 23, Казахстан

^* E-mail: alexander_prokopenya@sggw.pl
^** E-mail: minglibayev@gmail.com
^*** E-mail: kosherbaevaayken@gmail.com

Поступила в редакцию 28.07.2021
После доработки 09.08.2021
Принята к публикации 16.08.2021

EDN: PDPJIL
DOI: 10.31857/S0132347422020108

Полный текст (PDF)

Аннотация

Исследуется влияние переменности масс небесных тел на динамическую эволюцию планетных систем в случае, когда массы тел изменяются изотропно с различными скоростями, причем законы изменения масс считаются произвольными заданными функциями времени. В качестве модели многопланетной системы используется классическая задача n + 1 тел переменной массы, когда $n$ тел движутся вокруг центральной звезды по квазиэллиптическим непересекающимся орбитам и взаимодействуют друг с другом в соответствии с законом всемирного тяготения. Получены дифференциальные уравнения движения тел в оскулирующих элементах апериодического движения по квазиконическим сечениям. Обсуждается алгоритм вычисления возмущающей функции в виде степенных рядов по малым параметрам и получение дифференциальных уравнений, определяющих вековые возмущения орбитальных элементов. Все необходимые символьные вычисления выполняются с использованием системы компьютерной алгебры Wolfram Mathematica.

1. ВВЕДЕНИЕ

Движение сферически симмметричных тел, взаимодействующих друг с другом в соответствии с законом всемирного тяготения, обычно исследуется в рамках классической задачи многих тел, которая является базовой моделью небесной механики (см., напр., [1, 2]). Дифференциальные уравнения движения тел могут быть легко выписаны на основе второго закона Ньютона, но их общее решение известно только в случае двух взаимодействующих тел. Некоторые классы решений можно построить в виде степенных и тригонометрических рядов (см., например, [3–6]). Построение и исследование таких решений обычно связано с выполнением весьма громоздких символьных вычислений, что стимулировало разработку соответствующих аналитических и численно-аналитических методов и алгоритмов [7–14].

Открытие и изучение экзопланетных систем показало (см., например, [15–18]), что реальные небесные тела не являются стационарными, их масса, размеры, форма и внутренняя структура могут с течением времени изменяться. В нестационарном случае задача многих тел существенно усложняется, даже в простейшем случае двух тел переменной массы найти аналитическое решение уравнений движения удается только в специальных случаях (см. [19]). Поэтому для исследования нестационарных систем обычно применяют теорию возмущений, причем в качестве нулевого приближения используют точное решение нестационарной задачи двух тел (см. [20–26]).

Целью данной работы является исследование динамической эволюции планетных систем в случае, когда ведущим фактором эволюции является переменность масс небесных тел. Предполагается, что массы тел ${{P}_{0}},{{P}_{1}},\; \ldots ,\;{{P}_{n}}$ являются заданными функциями времени

${{m}_{0}} = {{m}_{0}}(t),{{m}_{1}} = {{m}_{1}}(t),\; \ldots ,\;{{m}_{n}} = {{m}_{n}}(t)$

и изменяются изотропно с различными скоростями, так что выполняются условия

$\frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}} \ne \frac{{{{{\dot {m}}}_{i}}}}{{{{m}_{i}}}},\quad \frac{{{{{\dot {m}}}_{i}}}}{{{{m}_{i}}}} \ne \frac{{{{{\dot {m}}}_{j}}}}{{{{m}_{j}}}},\quad i,j = 1,\;2,\; \ldots ,\quad n,i \ne j,$

а сами тела притягивают друг друга в соответствии с законом всемирного тяготения. Мы ограничимся наиболее интересным случаем планетарной проблемы n + 1 тел, когда $n$ планет ${{P}_{1}},\; \ldots ,\;{{P}_{n}}$ движутся вокруг центрального тела-звезды ${{P}_{0}}$ таким образом, что их орбиты не пересекаются [20, 27]. Основное внимание уделяется обсуждению вычислительных задач, возникающих при получении разложения возмущающей функции в степенной ряд по малым параметрам и определении эволюционных уравнений, для решения которых лучше всего использовать системы компьютерной алгебры. В данной работе все символьные вычисления выполняются с помощью системы компьютерной алгебры Wolfram Mathematica [28], которая имеет удобный интерфейс и позволяет легко комбинировать различные виды вычислений. Поскольку для использования системы Wolfram Mathematica требуется лицензия, для выполнения описанных в работе вычислений можно воспользоваться и другой свободно распространяемой системой компьютерной алгебы, в которой реализованы функции работы со степенными рядами.

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

Предполагая, что наиболее массивное тело ${{P}_{0}}$ находится в начале координат, и используя относительные декартовы координаты ${{\vec {R}}_{i}} = ({{X}_{i}},{{Y}_{i}},{{Z}_{i}})$, уравнения движения тел ${{P}_{1}},\; \ldots ,\;{{P}_{n}}$ можно записать в виде (см. [20, 26])

(1)

$\begin{gathered} {{{\ddot {\vec {R}}}}_{i}} + f({{m}_{0}} + {{m}_{i}})\frac{{{{{\vec {R}}}_{i}}}}{{R_{i}^{3}}} = \\ = \;f\sum\limits_{j = 1(j \ne i)}^n {{{m}_{j}}} \left( {\frac{{{{{\vec {R}}}_{j}} - {{{\vec {R}}}_{i}}}}{{R_{{ij}}^{3}}} - \frac{{{{{\vec {R}}}_{j}}}}{{R_{j}^{3}}}} \right),\quad i = 1,\;2,\; \ldots ,\;n, \\ \end{gathered} $

где f – гравитационная постоянная,

${{R}_{{ij}}} = \sqrt {{{{({{X}_{j}} - {{X}_{i}})}}^{2}} + {{{({{Y}_{j}} - {{Y}_{i}})}}^{2}} + {{{({{Z}_{j}} - {{Z}_{i}})}}^{2}}} ,$

${{R}_{i}} = \sqrt {X_{i}^{2} + Y_{i}^{2} + Z_{i}^{2}} ,\quad i,j = 1,\;2, \ldots ,\;n,$

а точка над символом в уравнениях (1) означает полную производную соответствующей функции по времени. Отметим, что при изотропном изменении масс тел реактивные силы не возникают и уравнения (1) с переменными массами выглядят точно так же, как и в случае постоянных масс (см. [1. 2]).

Поскольку найти общее решение системы (1) не представляется возможным даже при постоянных массах тел, для ее исследования воспользуемся теорией возмущений. Если пренебречь взаимодействием тел ${{P}_{1}},{{P}_{2}},\; \ldots ,\;{{P}_{n}}$ между собой и обнулить правую часть в (1), в нулевом приближении получим $n$ независимых задач двух тел, которые при постоянных массах легко интегрируются. Поскольку при переменных массах классическая задача двух тел в общем случае не интегрируется, перепишем уравнения (1) в виде

(2)

${{\ddot {\vec {R}}}_{i}} + f({{m}_{0}} + {{m}_{i}})\frac{{{{{\vec {R}}}_{i}}}}{{R_{i}^{3}}} - \frac{{{{{\ddot {\gamma }}}_{i}}}}{{{{\gamma }_{i}}}}{{\vec {R}}_{i}} = {{\vec {\nabla }}_{{{{{\vec {R}}}_{i}}}}}{{W}_{i}},$

где функции ${{W}_{i}}$ в правой части уравнений (2) имеют вид (см. [20, 26])

(3)

${{W}_{i}} = f\sum\limits_{j = 1(j \ne i)}^n {{{m}_{j}}} \left( {\frac{1}{{{{R}_{{ij}}}}} - \frac{{{{{\vec {R}}}_{i}} \cdot {{{\vec {R}}}_{j}}}}{{R_{j}^{3}}}} \right) - \frac{{{{{\ddot {\gamma }}}_{i}}}}{{2{{\gamma }_{i}}}}\vec {R}_{i}^{2}.$

Дважды непрерывно дифференцируемые функции ${{\gamma }_{i}}(t)$ в (2), (3) определяются соотношениями

${{\gamma }_{i}}(t) = \frac{{{{m}_{{00}}} + {{m}_{{i0}}}}}{{{{m}_{0}}(t) + {{m}_{i}}(t)}},$

где ${{m}_{{00}}} = {{m}_{0}}({{t}_{0}})$, ${{m}_{{i0}}} = {{m}_{i}}({{t}_{0}})$ – значения масс тел в начальный момент времени ${{t}_{0}}$. Отметим, что законы изменения масс во времени ${{m}_{j}}(t)$, ($j = 0,1$, ..., n) определяются на основе наблюдений за движением небесных тел и считаются известными. При постоянных массах производные функций ${{\gamma }_{i}}(t)$ в (2), (3) исчезают и уравнения (2) сводятся к классической задаче многих тел (см. [2]). Добавление в левой части уравнений (2) слагаемого, пропорционального второй производной функции ${{\gamma }_{i}}(t)$, обеспечивает интегрируемость уравнений (2) при ${{W}_{i}} = 0$, $i = 1,\;2,\; \ldots ,\;n$, что позволяет получить в нулевом приближении $n$ независимых уравнений, точные решения которых определяют движение тел ${{P}_{1}},{{P}_{2}},\; \ldots ,\;{{P}_{n}}$ вокруг тела ${{P}_{0}}$ по квазиконическим сечениям (см. [20]). Соответствующие точные решения, которые будем использовать в качестве нулевого приближения, можно представить в виде

${{X}_{j}} = {{\gamma }_{j}}{{\rho }_{j}}(\cos {{u}_{j}}\cos {{\Omega }_{j}} - \sin {{u}_{j}}\sin {{\Omega }_{j}}\cos {{i}_{j}}),$

${{Y}_{j}} = {{\gamma }_{j}}{{\rho }_{j}}(\cos {{u}_{j}}\sin {{\Omega }_{j}} + \sin {{u}_{j}}\cos {{\Omega }_{j}}\cos {{i}_{j}}),$

(4)

${{Z}_{j}} = {{\gamma }_{j}}{{\rho }_{j}}\sin {{u}_{j}}\sin {{i}_{j}},$

где ${{u}_{j}} = {{\theta }_{j}} + {{\omega }_{j}}$,

(5)

${{\rho }_{j}} = \frac{{{{a}_{j}}(1 - e_{j}^{2})}}{{1 + {{e}_{j}}\cos {{\theta }_{j}}}},\quad (j = 1,\;2,\; \ldots ,\;n),$

переменная ${{\theta }_{j}}$ является аналогом истинной аномалии, а параметры ${{a}_{j}}$, ${{e}_{j}}$, ${{i}_{j}}$, ${{\omega }_{j}}$, ${{\Omega }_{j}}$, определяемые из начальных условий движения, соответствуют известным из классической задачи двух тел кеплеровским орбитальным элементам и являются аналогами большой полуоси, эксцентриситета, наклонения, долготы перицентра и долготы восходящего узла невозмущенной квазиэллиптической орбиты каждого из тел ${{P}_{j}}$ (см. [1, 2, 20]).

Отметим, что в случае постоянных масс, когда ${{\gamma }_{j}}(t) \equiv 1$, уравнения (4)–(5) определяют движение тел ${{P}_{i}}$, ($i = 1,\;2,\; \ldots ,\;n$) вокруг тела ${{P}_{0}}$ по коническим сечениям. Наличие в выражениях (4) масштабного множителя ${{\gamma }_{j}}(t)$, зависящего от времени, приводит к деформации конических сечений и непериодичности движения. Поэтому говорят, что решения уравнений движения (2) в случае ${{W}_{i}} = 0$ описывают апериодическое движение тел по квазиконическим сечениям (см. [20]).

Аналог истинной аномалии ${{\theta }_{j}}$ характеризует положение тела на орбите и определяется уравнением

$\begin{gathered} \int\limits_0^{{{\theta }_{j}}} {\frac{{d{{\theta }_{j}}}}{{{{{(1 + {{e}_{j}}\cos {{\theta }_{j}})}}^{2}}}}} = \\ = \;\frac{1}{{{{{(1 - e_{j}^{2})}}^{{3/2}}}}}({{E}_{j}} - {{e}_{j}}\sin {{E}_{j}}) = \frac{{{{M}_{j}}}}{{{{{(1 - e_{j}^{2})}}^{{3/2}}}}}, \\ \end{gathered} $

где M_j – аналог средней аномалии, определяемой выражением

(6)

${{M}_{j}} = \frac{{\sqrt {{{\mu }_{{j0}}}} }}{{a_{j}^{{3/2}}}}({{\Phi }_{j}}(t) - {{\Phi }_{j}}({{\tau }_{j}})),$

${{\Phi }_{j}}(t) = \int\limits_0^t {\frac{{dt}}{{\gamma _{j}^{2}(t)}}} ,$

${{\mu }_{{j0}}} = f({{m}_{{00}}} + {{m}_{{j0}}}),\quad j = 1,\;2,\; \ldots ,\;n,$

${{\tau }_{j}}$ – аналог времени прохождения тела ${{P}_{j}}$ через перицентр, а эксцентрическая аномалия ${{E}_{j}}$ связана с истинной аномалией ${{\theta }_{j}}$ соотношением

(7)

${\text{tg}}\frac{{{{\theta }_{j}}}}{2} = \sqrt {\frac{{1 + {{e}_{j}}}}{{1 - {{e}_{j}}}}} {\text{tg}}\frac{{{{E}_{j}}}}{2}.$

При заданных орбитальных параметрах ${{a}_{j}}$, ${{e}_{j}}$, ${{i}_{j}}$, ${{\omega }_{j}}$, ${{\Omega }_{j}}$, ${{\tau }_{j}}$ каждого из тел ${{P}_{j}}$, а также известных функциях ${{\gamma }_{j}}(t)$, которые зависят от законов изменения масс всех тел, уравнение (6) позволяет найти средние аномалии M_j как функции времени. Решая затем уравнение Кеплера,

(8)

${{E}_{j}} - {{e}_{j}}\sin {{E}_{j}} = {{M}_{j}},$

находим эксцентрические аномалии ${{E}_{j}}$ и, используя соотношение (7), вычисляем истинные аномалии ${{\theta }_{j}}$. В результате выражения (4), (5) позволяют вычислить относительные декартовы координаты тел ${{P}_{j}}$ и полностью описать их невозмущенное движение.

В рассматриваемом случае многопланетной задачи предполагается, что масса центрального тела ${{P}_{0}}$ значительно превышает массы тел ${{P}_{j}}$ ($j = 1,\;2,\; \ldots ,\;n$). Потому их орбиты будут квазиконическими сечениями, определяемыми уравнениями (4)–(8), с малыми отклонениями, которые возникают вследствие взаимного гравитационного притяжения тел ${{P}_{i}}$ и ${{P}_{j}}$, ($i,j = 1,\;2,\; \ldots ,\;n$), а также изменения их масс. Эти факторы определяются возмущающими функциями ${{W}_{i}} \ne 0$, градиенты которых находятся в правых частях уравнений (2). Чтобы получить дифференциальные уравнения, определяющие временную эволюцию орбитальных параметров, удобно перейти к новым каноническим переменным, известным как вторая система канонических элементов Пуанкаре, и переписать уравнения (2) в канонической форме (см. [20, 27]). Используя решения (4) и выполняя стандартные, но довольно громоздкие символьные преобразования с применением встроенных функций $D$, $Expand$, $Replace$, $Simplify$, $Collect$, $Coefficient$, $Solve$ (см. [22, 28]), определяем три пары канонически сопряженных координат и импульсов $({{\lambda }_{j}},{{\Lambda }_{j}})$, $({{\eta }_{j}},{{\xi }_{j}})$, $({{q}_{j}},{{p}_{j}})$ для каждого из тел ${{P}_{j}}$, которые связаны с аналогами кеплеровских элементов орбиты соотношениями

(9)

${{\lambda }_{j}} = {{M}_{j}} + {{\pi }_{j}},\quad {{\Lambda }_{j}} = \sqrt {{{\mu }_{{j0}}}{{a}_{j}}} ,$

${{\eta }_{j}} = - \sqrt {2\sqrt {{{\mu }_{{j0}}}{{a}_{j}}} (1 - \sqrt {1 - e_{j}^{2}} )} \sin {{\pi }_{j}},$

(10)

${{\xi }_{j}} = \sqrt {2\sqrt {{{\mu }_{{j0}}}{{a}_{j}}} (1 - \sqrt {1 - e_{j}^{2}} )} \cos {{\pi }_{j}},$

${{q}_{j}} = - \sqrt {2\sqrt {{{\mu }_{{j0}}}{{a}_{j}}} \sqrt {1 - e_{j}^{2}} \left( {1 - \cos {{i}_{j}}} \right)} \sin {{\Omega }_{j}},$

(11)

${{p}_{j}} = \sqrt {2\sqrt {{{\mu }_{{j0}}}{{a}_{j}}} \sqrt {1 - e_{j}^{2}} \left( {1 - \cos {{i}_{j}}} \right)} \cos {{\Omega }_{j}},$

где

(12)

${{\pi }_{j}} = {{\Omega }_{j}} + {{\omega }_{j}}.$

Дифференциальные уравнения движения $n$ тел в оскулирующих аналогах второй системы переменных Пуанкаре (9)–(11) имеют канонический вид (см. [20])

${{\dot {\lambda }}_{j}} = \frac{{\partial R_{j}^{*}}}{{\partial {{\Lambda }_{j}}}},\quad {{\dot {\eta }}_{j}} = \frac{{\partial R_{j}^{*}}}{{\partial {{\xi }_{j}}}},\quad {{\dot {q}}_{j}} = \frac{{\partial R_{j}^{*}}}{{\partial {{p}_{j}}}},$

(13)

${{\dot {\Lambda }}_{j}} = - \frac{{\partial R_{j}^{*}}}{{\partial {{\lambda }_{j}}}},\quad {{\dot {\xi }}_{j}} = - \frac{{\partial R_{j}^{*}}}{{\partial {{\eta }_{j}}}},\quad {{\dot {p}}_{j}} = - \frac{{\partial R_{j}^{*}}}{{\partial {{q}_{j}}}},$

где функции Гамильтона

(14)

$R_{j}^{*} = - \frac{{\mu _{{j0}}^{2}}}{{2\gamma _{j}^{2}\Lambda _{j}^{2}}} - {{W}_{j}}(t,{{\Lambda }_{j}},{{\xi }_{j}},{{p}_{j}},{{\lambda }_{j}},{{\eta }_{j}},{{q}_{j}}),$

а возмущающие функции W_j, которые определяются соотношениями (3), должны быть выражены через новые канонические переменные (${{\lambda }_{j}}$, ${{\eta }_{j}}$, ${{q}_{j}}$, ${{\Lambda }_{j}}$, ${{\xi }_{j}}$, ${{p}_{j}}$). Отметим, что канонические уравнения возмущенного движения (13) наиболее удобны для описания динамической эволюции планетных систем в случае, когда аналоги эксцентриситетов ${{e}_{j}}$ и аналоги наклонений ${{i}_{j}}$ орбитальной плоскости достаточно малы.

При подстановке функции Гамильтона (14) в уравнения (13) получаем

${{\dot {\lambda }}_{j}} = \frac{{\mu _{{j0}}^{2}}}{{\gamma _{j}^{2}\Lambda _{j}^{3}}} - \frac{{\partial {{W}_{j}}}}{{\partial {{\Lambda }_{j}}}},\quad {{\dot {\Lambda }}_{j}} = \frac{{\partial {{W}_{j}}}}{{\partial {{\lambda }_{j}}}},$

(15)

${{\dot {\eta }}_{j}} = - \frac{{\partial {{W}_{j}}}}{{\partial {{\xi }_{j}}}},\quad {{\dot {\xi }}_{j}} = \frac{{\partial {{W}_{j}}}}{{\partial {{\eta }_{j}}}},$

${{\dot {q}}_{j}} = - \frac{{\partial {{W}_{j}}}}{{\partial {{p}_{j}}}},\quad {{\dot {p}}_{j}} = \frac{{\partial {{W}_{j}}}}{{\partial {{q}_{j}}}}.$

Из уравнений (15) легко видеть, что при ${{W}_{j}}$ = 0 канонические переменные (${{\eta }_{j}}$, ${{q}_{j}}$, ${{\Lambda }_{j}}$, ${{\xi }_{j}}$, ${{p}_{j}}$), а также орбитальные элементы ${{a}_{j}}$, ${{e}_{j}}$, ${{i}_{j}}$, ${{\omega }_{j}}$, ${{\Omega }_{j}}$, остаются постоянными, а средняя долгота ${{\lambda }_{j}}$ является возрастающей функцией времени.

3. ВЫРАЖЕНИЕ ВОЗМУЩАЮЩИХ ФУНКЦИЙ ЧЕРЕЗ ОСКУЛИРУЮЩИЕ ЭЛЕМЕНТЫ

Далее будем предполагать, что во время движения выполняются условия ${{i}_{j}} \ll 1$ и ${{e}_{j}} \ll 1$, т.е. тела движутся вблизи экваториальной плоскости по возмущенным квазиконическим сечениям с малым эксцентриситетом. Как следует из (10)–(11), в этом приближении канонические переменные (${{\eta }_{j}}$, ${{q}_{j}}$, ${{\xi }_{j}}$, ${{p}_{j}}$) являются малыми величинами. Действительно, из выражений (9)–(11) получаем точные соотношения

(16)

$e_{j}^{2} = \frac{{\eta _{j}^{2} + \xi _{j}^{2}}}{{{{\Lambda }_{j}}}}\left( {1 - \frac{{\xi _{j}^{2} + \xi _{j}^{2}}}{{4{{\Lambda }_{j}}}}} \right),$

(17)

$\mathop {\sin }\nolimits^2 {{i}_{j}} = \frac{{p_{j}^{2} + q_{j}^{2}}}{{{{\Lambda }_{j}}\sqrt {1 - e_{j}^{2}} }}\left( {1 - \frac{{p_{j}^{2} + q_{j}^{2}}}{{4{{\Lambda }_{j}}\sqrt {1 - e_{j}^{2}} }}} \right).$

Из (16), (17) следует, что при ${{e}_{j}} \ll 1$, ${{i}_{j}} \ll 1$ переменные ${{\eta }_{j}}$, ${{\xi }_{j}}$ и эксцентриситет ${{e}_{j}}$, а также переменные ${{q}_{j}}$, ${{p}_{j}}$ и наклонение i_j являются величинами одного порядка. Поэтому уравнения возмущенного движения (15) можно упростить, заменяя возмущающие функции W_j их соответствующими разложениями в степенные ряды по переменным ${{\eta }_{j}}$, ${{\xi }_{j}}$, ${{q}_{j}}$, ${{p}_{j}}$ (см. [2]). Следует отметить, что разложение возмущающих функций W_j в ряды требует выполнения довольно громоздких символьных вычислений и представляет собой нетривиальную задачу. Далее приведем весьма эффективный метод ее решения с помощью системы компьютерной алгебры Wolfram Mathematica.

Сначала выразим из соотношений (10) переменные ${{e}_{j}}$ и ${{\pi }_{j}}$ через ${{\eta }_{j}}$ и ${{\xi }_{j}}$. Используя разложение

$\sqrt {1 - \sqrt {1 - e_{j}^{2}} } = \frac{{{{e}_{j}}}}{{\sqrt 2 }}\left( {1 + \frac{{e_{j}^{2}}}{8} + \frac{{7e_{j}^{4}}}{{128}} + \; \ldots } \right),$

подставляем вместо $e_{j}^{2}$ его точное выражение из (16) и, выполняя в (10) разложение в степенной ряд по переменным ${{\eta }_{j}}$, ${{\xi }_{j}}$, находим:

${{e}_{j}}\sin {{\pi }_{j}} = - \frac{{{{\eta }_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }}\left( {1 - \frac{{\eta _{j}^{2} + \xi _{j}^{2}}}{{8{{\Lambda }_{j}}}} + \; \ldots } \right),$

(18)

${{e}_{j}}\cos {{\pi }_{j}} = \frac{{{{\xi }_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }}\left( {1 - \frac{{\eta _{j}^{2} + \xi _{j}^{2}}}{{8{{\Lambda }_{j}}}} + \; \ldots } \right).$

Разложения (18) получены с помощью встроенной функции $Series$ (см. [28]) с точностью до четвертого порядка. Поскольку разложения производятся независимо по двум переменным ${{\eta }_{j}}$ и ${{\xi }_{j}}$, в получаемом выражении могут содержаться произведения $\eta _{j}^{k}\xi _{j}^{m}$, ($k,m = 0,\;1,\; \ldots ,\;4$) вплоть до восьмого порядка включительно, например, $\eta _{j}^{4}\xi _{j}^{4}$. Чтобы ограничиться только членами четвертого порядка, в разлагаемой функции удобно сделать подстановку ${{\eta }_{j}} \to \varepsilon {{\eta }_{j}}$ и ${{\xi }_{j}} \to \varepsilon {{\xi }_{j}}$ и выполнить разложение функции по переменной $\varepsilon $ в окрестности точки $\varepsilon = 0$ с точностью до четвертого порядка, применяя функцию $Series$. В результате получим многочлен четвертого порядка по переменной $\varepsilon $, причем коэффициенты при ${{\varepsilon }^{j}}$, ($j = 0,1$, ..., 4) будут содержать только произведения $\eta _{j}^{k}\xi _{j}^{m}$, ($k,m = 0,1$, ..., 4) j-го порядка, удовлетворяющие условию $k + m = j$. Поскольку переменная $\varepsilon $ была введена только для упрощения вычислений, в полученном выражении достаточно положить $\varepsilon = 1$.

Аналогичным образом из уравнений (11) с учетом (17) находим разложения по степеням ${{\xi }_{j}}$, ${{\eta }_{j}}$, ${{p}_{j}}$ и ${{q}_{j}}$ для следующих выражений:

$\sin {{i}_{j}}\sin {{\Omega }_{j}} = - \frac{{{{q}_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }}\left( {1 - \frac{{p_{j}^{2} + q_{j}^{2}}}{{8{{\Lambda }_{j}}}} + \frac{{\eta _{j}^{2} + \xi _{j}^{2}}}{{4{{\Lambda }_{j}}}}} \right),$

(19)

$\sin {{i}_{j}}\cos {{\Omega }_{j}} = \frac{{{{p}_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }}\left( {1 - \frac{{p_{j}^{2} + q_{j}^{2}}}{{8{{\Lambda }_{j}}}} + \frac{{\eta _{j}^{2} + \xi _{j}^{2}}}{{4{{\Lambda }_{j}}}}} \right).$

Отметим, что соотношения (18), (19) представляют собой разложения в степенные ряды по ${{\xi }_{j}}$, ${{\eta }_{j}}$, ${{p}_{j}}$ и ${{q}_{j}}$ не самих переменых ${{e}_{j}}$, ${{\pi }_{j}}$, ${{i}_{j}}$, ${{\Omega }_{j}}$, а произведений ${{e}_{j}}\sin {{\pi }_{j}}$, ${{e}_{j}}\cos {{\pi }_{j}}$, $\sin {{i}_{j}}\sin {{\Omega }_{j}}$, $\sin {{i}_{j}}\cos {{\Omega }_{j}}$, которые входят в выражения (4) для декартовых координат тел ${{X}_{j}}$, ${{Y}_{j}}$, ${{Z}_{j}}$.

Кроме переменых ${{e}_{j}}$, ${{\omega }_{j}}$, ${{i}_{j}}$, ${{\Omega }_{j}}$ выражения (4) содержат аналог истинной аномалии ${{\theta }_{j}}$, которую также следует выразить через элементы Пуанкаре. Для этого можно воспользоваться известным решением уравнения Кеплера (8), которое определяет эксцентрическую аномалию E_j в виде сходящегося при ${{e}_{j}} \ll 1$ степенного ряда (см. [2])

(20)

${{E}_{j}} = {{M}_{j}} + {{e}_{j}}\sin {{M}_{j}} + \frac{{e_{j}^{2}}}{2}\sin (2{{M}_{j}}) + \; \ldots $

Подчеркнем, что встроенная функция $Series$ существенно упрощает вычисления и позволяет производить разложения выражений в степенные ряды с любой требуемой точностью (см. [28]). Для получения коэффициентов разложения (20), например, с точностью до $n$-го порядка по ${{e}_{j}}$ можно использовать функцию $solE[n]$, записанную на языке Wolfram Mathematica:

$solE[n\_]: = w - {{e}_{j}}Sin[w] - {{M}_{j}}{\text{/}}{\text{.}}w \to $

$Sum[e_{j}^{k}{{w}_{k}},\{ k,0,n\} ]//Series[\# ,\{ {{e}_{j}},0,n\} ]\& //$

$Normal//Collect\left[ {\# ,{{e}_{j}},Simplify} \right]\& ;$

Выполнение команды $solE[3]$, например, приводит к результату

$\begin{gathered} - {{M}_{j}} + {{w}_{0}} + {{e}_{j}}( - Sin[{{w}_{0}}] + {{w}_{1}}) + \\ + \;e_{j}^{2}( - Cos[{{w}_{0}}]{{w}_{1}} + {{w}_{2}}) + \\ + \;e_{j}^{3}\left( {\frac{1}{2}Sin[{{w}_{0}}]w_{1}^{2} - Cos[{{w}_{0}}]{{w}_{2}} + {{w}_{3}}} \right) \\ \end{gathered} $

Выделяя в полученном выражении с помощью функции $Coefficient$ коэффициенты при $e_{j}^{k}$, $k = 0$, 1, 2, ..., и приравнивая их к нулю, получаем систему уравнений, из которой последовательно находим коэффициенты ${{w}_{k}}$ разложения (20).

Поскольку увеличение точности вычислений приводит ко все более громоздким символьным выражениям, далее ограничимся вычислениями с точностью до второго порядка по переменным ${{\xi }_{j}}$, ${{\eta }_{j}}$, ${{p}_{j}}$ и ${{q}_{j}}$. Используя (20), находим

$\cos {{E}_{j}} = \cos {{M}_{j}} - {{e}_{j}}{\text{si}}{{{\text{n}}}^{2}}{{M}_{j}} - \frac{{3e_{j}^{2}}}{2}\cos {{M}_{j}}{\text{si}}{{{\text{n}}}^{2}}{{M}_{j}},$

(21)

$\begin{gathered} \sin {{E}_{j}} = \sin {{M}_{j}} + \frac{{{{e}_{j}}}}{2}\sin (2{{M}_{j}}) + \\ + \;\frac{{e_{j}^{2}}}{8}\left( {3\sin (3{{M}_{j}}) - \sin {{M}_{j}}} \right). \\ \end{gathered} $

Учитывая (7), получаем

$\cos {{\theta }_{j}} = \frac{{\cos {{E}_{j}} - {{e}_{j}}}}{{1 - {{e}_{j}}\cos {{E}_{j}}}} = \cos {{M}_{j}} - $

$ - \;{{e}_{j}}(1 - \cos (2{{M}_{j}})) + \frac{9}{8}e_{j}^{2}\left( {\cos (3{{M}_{j}}) - \cos {{M}_{j}}} \right),$

(22)

$\begin{gathered} \sin {{\theta }_{j}} = \frac{{\sin {{E}_{j}}\sqrt {1 - e_{j}^{2}} }}{{1 - {{e}_{j}}\cos {{E}_{j}}}} = \sin {{M}_{j}} + \\ + \;{{e}_{j}}\sin (2{{M}_{j}}) + \frac{{e_{j}^{2}}}{8}(9\sin (3{{M}_{j}}) - 7\sin {{M}_{j}}). \\ \end{gathered} $

Принимая во внимание (7), (9), (12), (16), (18), (21) и выполняя необходимые подстановки и разложения, перепишем (5) в виде

(23)

$\begin{gathered} \frac{{{{\rho }_{j}}}}{{{{a}_{j}}}} = 1 - {{e}_{j}}\cos {{E}_{j}} = 1 + \frac{{{{\eta }_{j}}\sin {{\lambda }_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }} - \frac{{{{\xi }_{j}}\cos {{\lambda }_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }} + \\ + \;\frac{{\eta _{j}^{2}}}{{2{{\Lambda }_{j}}}}\left( {1 + \cos (2{{\lambda }_{j}})} \right) + \frac{{\xi _{j}^{2}}}{{2{{\Lambda }_{j}}}}\left( {1 - \cos (2{{\lambda }_{j}})} \right) + \\ + \;\frac{{{{\eta }_{j}}{{\xi }_{j}}}}{{{{\Lambda }_{j}}}}\sin (2{{\lambda }_{j}}). \\ \end{gathered} $

Наконец, используя разложения (16)–(23), выражения (4) для декартовых координат тел ${{P}_{j}}$ с точностью до второго порядка по переменным Пуанкаре ${{\eta }_{j}}$, ${{\xi }_{j}}$, ${{q}_{j}}$, ${{p}_{j}}$ получаем в виде

$\frac{{{{X}_{j}}}}{{{{\gamma }_{j}}{{a}_{j}}}} = \cos {{\lambda }_{j}} - \frac{{{{\xi }_{j}}}}{{2\sqrt {{{\Lambda }_{j}}} }}(3 - \cos (2{{\lambda }_{j}})) - $

$ - \;\frac{{{{\eta }_{j}}}}{{2\sqrt {{{\Lambda }_{j}}} }}\sin (2{{\lambda }_{j}}) + \frac{{3\xi _{j}^{2}}}{{8{{\Lambda }_{j}}}}\left( {\cos (3{{\lambda }_{j}}) - \cos {{\lambda }_{j}}} \right) - $

$ - \;\frac{{{{\xi }_{j}}{{\eta }_{j}}}}{{4{{\Lambda }_{j}}}}\left( {3\sin (3{{\lambda }_{j}}) + \sin {{\lambda }_{j}}} \right) - \frac{{\eta _{j}^{2}}}{{8{{\Lambda }_{j}}}}\left( {3\cos (3{{\lambda }_{j}})} \right. + $

$\left. { - \;5\cos {{\lambda }_{j}}} \right) - \frac{{{{p}_{j}}{{q}_{j}}}}{{2{{\Lambda }_{j}}}}\sin {{\lambda }_{j}} - \frac{{q_{j}^{2}}}{{2{{\Lambda }_{j}}}}\cos {{\lambda }_{j}},$

$\frac{{{{Y}_{j}}}}{{{{\gamma }_{j}}{{a}_{j}}}} = \sin {{\lambda }_{j}} + \frac{{{{\eta }_{j}}}}{{2\sqrt {{{\Lambda }_{j}}} }}(3 + \cos (2{{\lambda }_{j}})) + $

$ + \;\frac{{{{\xi }_{j}}}}{{2\sqrt {{{\Lambda }_{j}}} }}\sin (2{{\lambda }_{j}}) - \frac{{3\eta _{j}^{2}}}{{8{{\Lambda }_{j}}}}(\sin (3{{\lambda }_{j}}) + \sin {{\lambda }_{j}}) + $

$ + \;\frac{{{{\xi }_{j}}{{\eta }_{j}}}}{{4{{\Lambda }_{j}}}}(3\cos (3{{\lambda }_{j}}) - \cos {{\lambda }_{j}}) + \frac{{\xi _{j}^{2}}}{{8{{\Lambda }_{j}}}}(3\sin (3{{\lambda }_{j}}) - $

$ - \;5\sin {{\lambda }_{j}}) - \frac{{{{p}_{j}}{{q}_{j}}}}{{2{{\Lambda }_{j}}}}\cos {{\lambda }_{j}} - \frac{{p_{j}^{2}}}{{2{{\Lambda }_{j}}}}\sin {{\lambda }_{j}},$

$\frac{{{{Z}_{j}}}}{{{{\gamma }_{j}}{{a}_{j}}}} = \frac{{{{p}_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }}\sin {{\lambda }_{j}} + \frac{{{{q}_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }}\cos {{\lambda }_{j}} + $

$ + \;\frac{{{{\xi }_{j}}{{p}_{j}}}}{{2{{\Lambda }_{j}}}}\sin (2{{\lambda }_{j}}) - \frac{{{{\xi }_{j}}{{q}_{j}}}}{{2{{\Lambda }_{j}}}}(3 - \cos (2{{\lambda }_{j}})) + $

(24)

$ + \frac{{{{\eta }_{j}}{{p}_{j}}}}{{2{{\Lambda }_{j}}}}(3 + \cos (2{{\lambda }_{j}})) - \frac{{{{\eta }_{j}}{{q}_{j}}}}{{2{{\Lambda }_{j}}}}\sin (2{{\lambda }_{j}}).$

Используя (4), (23), (24) и применяя функцию Series, получаем следующие выражения, которые будут использованы далее при вычислении разложений возмущающих функций ${{W}_{j}}$ (см. (3)):

$\frac{{R_{j}^{2}}}{{\gamma _{j}^{2}a_{j}^{2}}} = 1 + \frac{2}{{\sqrt {{{\Lambda }_{j}}} }}({{\eta }_{j}}\sin {{\lambda }_{j}} - {{\xi }_{j}}\cos {{\lambda }_{j}}) + $

$ + \;\frac{{\eta _{j}^{2}}}{{2{{\Lambda }_{j}}}}(3 + \cos (2{{\lambda }_{j}})) + \frac{{{{\eta }_{j}}{{\xi }_{j}}}}{{{{\Lambda }_{j}}}}\sin (2{{\lambda }_{j}}) + $

(25)

$ + \;\frac{{\xi _{j}^{2}}}{{2{{\Lambda }_{j}}}}(3 - \cos (2{{\lambda }_{j}})),$

$\frac{{\gamma _{j}^{3}a_{j}^{3}}}{{R_{j}^{3}}} = 1 - \frac{{3{{\eta }_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }}\sin {{\lambda }_{j}} + \frac{{3{{\xi }_{j}}}}{{\sqrt {{{\Lambda }_{j}}} }}\cos {{\lambda }_{j}} + $

$ + \;\frac{{3\eta _{j}^{2}}}{{2{{\Lambda }_{j}}}}(1 - 3\cos (2{{\lambda }_{j}})) - \frac{{9{{\eta }_{j}}{{\xi }_{j}}}}{{{{\Lambda }_{j}}}}\sin (2{{\lambda }_{j}}) + $

(26)

$ + \;\frac{{3\xi _{j}^{2}}}{{2{{\Lambda }_{j}}}}(1 + 3\cos (2{{\lambda }_{j}})).$

$\frac{{{{{\vec {R}}}_{i}} \cdot {{{\vec {R}}}_{j}}}}{{{{a}_{i}}{{a}_{j}}{{\gamma }_{i}}{{\gamma }_{j}}}} = \cos ({{\lambda }_{i}} - {{\lambda }_{j}}) + \frac{{{{\eta }_{i}}}}{{2\sqrt {{{\Lambda }_{i}}} }}(3\sin {{\lambda }_{j}} - $

$ - \;\sin (2{{\lambda }_{i}} - {{\lambda }_{j}})) + \frac{{{{\eta }_{j}}}}{{2\sqrt {{{\Lambda }_{j}}} }}(3\sin {{\lambda }_{i}} + \sin ({{\lambda }_{i}} - 2{{\lambda }_{j}})) - $

$ - \;\frac{{{{\xi }_{i}}}}{{2\sqrt {{{\Lambda }_{i}}} }}(3\cos {{\lambda }_{j}} - \cos (2{{\lambda }_{i}} - {{\lambda }_{j}})) - \frac{{{{\xi }_{j}}}}{{2\sqrt {{{\Lambda }_{j}}} }}(3\cos {{\lambda }_{i}} - $

$\begin{gathered} - \;{\text{cos}}({{\lambda }_{i}} - 2{{\lambda }_{j}})) - \\ - \frac{{\eta _{i}^{2}}}{{4{{\Lambda }_{i}}}}{\text{cos}}{{\lambda }_{i}}(3{\text{cos}}(2{{\lambda }_{i}} - {{\lambda }_{j}}) + {\text{cos}}{{\lambda }_{j}}) - \\ \end{gathered} $

$ - \;\frac{{\eta _{j}^{2}}}{{4{{\Lambda }_{j}}}}\cos {{\lambda }_{j}}(3\cos ({{\lambda }_{i}} - 2{{\lambda }_{j}}) + \cos {{\lambda }_{i}}) - $

$ - \;\frac{{\xi _{i}^{2}}}{{4{{\Lambda }_{i}}}}\sin {{\lambda }_{i}}(3\sin (2{{\lambda }_{i}} - {{\lambda }_{j}}) + \sin {{\lambda }_{j}}) + $

$ + \;\frac{{\xi _{j}^{2}}}{{4{{\Lambda }_{j}}}}\sin {{\lambda }_{j}}(3\sin ({{\lambda }_{i}} - 2{{\lambda }_{j}}) - \sin {{\lambda }_{i}}) + $

$\begin{gathered} + \;\frac{{{{\eta }_{i}}{{\eta }_{j}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{j}}} }}(9 + 3\cos (2{{\lambda }_{i}}) + \\ + \;\cos (2({{\lambda }_{i}} - {{\lambda }_{j}})) + 3\cos (2{{\lambda }_{j}})) + \\ \end{gathered} $

$\begin{gathered} + \;\frac{{{{\xi }_{i}}{{\xi }_{j}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{j}}} }}(9 - 3\cos (2{{\lambda }_{i}}) + \\ + \;\cos (2({{\lambda }_{i}} - {{\lambda }_{j}})) - 3\cos (2{{\lambda }_{j}})) - \\ \end{gathered} $

$ - \;\frac{{{{\eta }_{i}}{{\xi }_{i}}}}{{4{{\Lambda }_{i}}}}(3\sin (3{{\lambda }_{i}} - {{\lambda }_{j}}) + \sin ({{\lambda }_{i}} + {{\lambda }_{j}})) + $

$ + \;\frac{{{{\eta }_{j}}{{\xi }_{j}}}}{{4{{\Lambda }_{j}}}}(3\sin ({{\lambda }_{i}} - 3{{\lambda }_{j}}) - \sin ({{\lambda }_{i}} + {{\lambda }_{j}})) + $

$ + \;\frac{{{{\eta }_{i}}{{\xi }_{j}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{j}}} }}(3\sin (2{{\lambda }_{i}}) + 3\sin (2{{\lambda }_{j}}) - \sin (2({{\lambda }_{i}} - {{\lambda }_{j}}))) + $

$ + \;\frac{{{{\eta }_{j}}{{\xi }_{i}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{j}}} }}(3\sin (2{{\lambda }_{i}}) + 3\sin (2{{\lambda }_{j}}) + \sin (2({{\lambda }_{i}} - {{\lambda }_{j}}))) - $

$ - \;\frac{1}{{2{{\Lambda }_{i}}}}(q_{i}^{2}\cos {{\lambda }_{i}}\cos {{\lambda }_{j}} + p_{i}^{2}\sin {{\lambda }_{i}}\sin {{\lambda }_{j}}) - $

$ - \;\frac{1}{{2{{\Lambda }_{j}}}}(q_{j}^{2}\cos {{\lambda }_{i}}\cos {{\lambda }_{j}} + p_{j}^{2}\sin {{\lambda }_{i}}\sin {{\lambda }_{j}}) - $

$ - \;\frac{{{{q}_{i}}{{p}_{i}}}}{{2{{\Lambda }_{i}}}}\sin ({{\lambda }_{i}} + {{\lambda }_{j}}) - \frac{{{{q}_{j}}{{p}_{j}}}}{{2{{\Lambda }_{j}}}}\sin ({{\lambda }_{i}} + {{\lambda }_{j}}) + $

$ + \;\frac{{{{q}_{i}}{{p}_{j}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{j}}} }}\cos {{\lambda }_{i}}\sin {{\lambda }_{j}} + \frac{{{{q}_{j}}{{p}_{i}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{j}}} }}\cos {{\lambda }_{j}}\sin {{\lambda }_{i}} + $

(27)

$ + \;\frac{{{{q}_{i}}{{q}_{j}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{j}}} }}\cos {{\lambda }_{i}}\cos {{\lambda }_{j}} + \frac{{{{p}_{i}}{{p}_{j}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{j}}} }}\sin {{\lambda }_{i}}\sin {{\lambda }_{j}}.$

Используя (25), (27) и соотношение

${{R}_{{ij}}} = \sqrt {R_{i}^{2} + R_{j}^{2} - 2{{{\vec {R}}}_{i}} \cdot {{{\vec {R}}}_{j}}} ,$

вычисляем ${{R}_{{ij}}}$ и далее находим разложение выражения $1{\text{/}}{{R}_{{ij}}}$, которое описывает возмущения, связанные с гравитационным взаимодействием тел ${{P}_{i}}$ и ${{P}_{j}}$ (см. (3)). Поскольку получаемое выражение более громоздко, мы его не приводим. Отметим только, что в результате получается многочлен второго порядка относительно переменных Пуанкаре ${{\eta }_{j}}$, ${{\xi }_{j}}$, ${{q}_{j}}$, ${{p}_{j}}$, коэффициенты которого являются периодическими функциями средних долгот ${{\lambda }_{j}}$ и представляют собой рациональные выражения, в числителях которых содержатся тригонометрические функции $\sin (k{{\lambda }_{j}})$, $\cos (k{{\lambda }_{j}})$, $\sin (k{{\lambda }_{i}} \pm n{{\lambda }_{j}})$, $\cos (k{{\lambda }_{i}} \pm n{{\lambda }_{j}})$, ($k,n = 1,\;2,\;3,\; \ldots $), а знаменатели содержат выражения ${{\Delta }_{{ij}}}$, $\Delta _{{ij}}^{3}$, $\Delta _{{ij}}^{5}$, где

${{\Delta }_{{ij}}} = {{(a_{i}^{2}\gamma _{i}^{2} + a_{j}^{2}\gamma _{j}^{2} - 2{{a}_{i}}{{a}_{j}}{{\gamma }_{i}}{{\gamma }_{j}}\cos ({{\lambda }_{i}} - {{\lambda }_{j}}))}^{{1/2}}}.$

4. ВЕКОВАЯ ЧАСТЬ ВОЗМУЩАЮЩИХ ФУНКЦИЙ

Подставляя в (3) разложение выражения $1{\text{/}}{{R}_{{ij}}}$ и выражения (25)–(27) и выполняя необходимые символьные вычисления, получаем довольно громоздкое выражение для возмущающих функций W_j в виде многочлена второго порядка относительно переменных Пуанкаре ${{\eta }_{j}},{{\xi }_{j}},{{q}_{j}},{{p}_{j}}$, которое мы здесь не приводим. Поскольку нас интересует эволюция орбитальных параметров тел ${{P}_{j}}$ на больших интервалах времени, короткопериодические возмущения, связанные с орбитальным движением тел, следует устранить путем усреднения возмущающих функций W_j по средним долготам ${{\lambda }_{j}}$ (см. [2, 20]). Усреднение многочленов W_j, приводящее к вековой части $W_{j}^{{(sec)}}$ возмущающих функций, сводится к вычислению интегралов

(28)

$W_{j}^{{(sec)}} = \frac{1}{{{{{(2\pi )}}^{n}}}}\int\limits_0^{2\pi } d {{\lambda }_{1}}\; \ldots \;\int\limits_0^{2\pi } d {{\lambda }_{n}}{{W}_{j}}.$

Поскольку возмущающие функции W_j представляют собой суммы большого числа слагаемых, каждое из которых может зависеть не более чем от двух средних долгот ${{\lambda }_{i}}$ и ${{\lambda }_{j}}$ (см. (3), (25)–(27)), интеграл (28) сводится к двукратному интегралу по переменным ${{\lambda }_{i}}$ и ${{\lambda }_{j}}$.

В результате получим следующие вековые части возмущающих функций:

$W_{i}^{{(sec)}} = f\sum\limits_{s = 1}^{i - 1} {{{m}_{s}}} \left( {\frac{{A_{0}^{{is}}}}{2} + \Pi _{{ii}}^{{is}}\frac{{\eta _{i}^{2} + \xi _{i}^{2}}}{{2{{\Lambda }_{i}}}} + } \right.$

$ + \;\Pi _{{is}}^{{is}}\frac{{{{\eta }_{i}}{{\eta }_{s}} + {{\xi }_{i}}{{\xi }_{s}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{s}}} }} + \Pi _{{ss}}^{{is}}\frac{{\eta _{s}^{2} + \xi _{s}^{2}}}{{2{{\Lambda }_{s}}}} - $

$ - \;\left. {B_{1}^{{is}}\left( {\frac{{p_{i}^{2} + q_{i}^{2}}}{{8{{\Lambda }_{i}}}} - \frac{{{{p}_{i}}{{p}_{s}} + {{q}_{i}}{{q}_{s}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{s}}} }} + \frac{{p_{s}^{2} + q_{s}^{2}}}{{8{{\Lambda }_{s}}}}} \right)} \right) + $

$ + \;f\sum\limits_{k = i + 1}^n {{m}_{k}}\left( {\frac{{A_{0}^{{ik}}}}{2} + \Pi _{{ii}}^{{ik}}\frac{{\eta _{i}^{2} + \xi _{i}^{2}}}{{2{{\Lambda }_{i}}}} + } \right.$

$ + \;\Pi _{{ik}}^{{ik}}\frac{{{{\eta }_{i}}{{\eta }_{k}} + {{\xi }_{i}}{{\xi }_{k}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{k}}} }} + \Pi _{{kk}}^{{ik}}\frac{{\eta _{k}^{2} + \xi _{k}^{2}}}{{2{{\Lambda }_{k}}}} - $

$\left. { - \;B_{1}^{{ik}}\left( {\frac{{p_{i}^{2} + q_{i}^{2}}}{{8{{\Lambda }_{i}}}} - \frac{{{{p}_{i}}{{p}_{k}} + {{q}_{i}}{{q}_{k}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{k}}} }} + \frac{{p_{k}^{2} + q_{k}^{2}}}{{8{{\Lambda }_{k}}}}} \right)} \right) - $

(29)

$ - \;\frac{{{{{\ddot {\gamma }}}_{i}}\Lambda _{i}^{4}}}{{2{{\gamma }_{i}}\mu _{{i0}}^{2}}}\left( {1 + \frac{3}{{2{{\Lambda }_{i}}}}(\xi _{i}^{2} + \eta _{i}^{2})} \right),$

где приняты следующие обозначения для внутренних планет ($1 \leqslant s < i$):

$\Pi _{{ii}}^{{is}} = - \frac{{3{{\alpha }_{{is}}}}}{4}B_{0}^{{is}} - \frac{1}{2}B_{1}^{{is}} + \frac{{15 + 6\alpha _{{is}}^{2}}}{8}C_{0}^{{is}} - $

$ - \;\frac{{3{{\alpha }_{{is}}}}}{2}C_{1}^{{is}} - \frac{9}{8}C_{2}^{{is}},$

$\Pi _{{is}}^{{is}} = \frac{1}{8}(9B_{0}^{{is}} + B_{2}^{{is}}) - \frac{{9(1 + \alpha _{{is}}^{2})}}{{8{{\alpha }_{{is}}}}}C_{0}^{{is}} + $

$ + \;\frac{{21}}{{16}}C_{1}^{{is}} + \frac{{3(1 + \alpha _{{is}}^{2})}}{{8{{\alpha }_{{is}}}}}C_{2}^{{is}} + \frac{3}{{16}}C_{3}^{{is}},$

$\Pi _{{ss}}^{{is}} = - \frac{3}{{4{{\alpha }_{{is}}}}}B_{0}^{{is}} - \frac{1}{2}B_{1}^{{is}} + \frac{{15\alpha _{{is}}^{2} + 6}}{{8\alpha _{{is}}^{2}}}C_{0}^{{is}} - $

(30)

$ - \;\frac{3}{{2{{\alpha }_{{is}}}}}C_{1}^{{is}} - \frac{9}{8}C_{2}^{{is}},$

и введен параметр

${{\alpha }_{{is}}} = \frac{{{{\gamma }_{s}}{{a}_{s}}}}{{{{\gamma }_{i}}{{a}_{i}}}} < 1\;.$

Ограничение на величину параметра ${{\alpha }_{{is}}}$ является следствием предположения, что при $1 \leqslant s < i$ траектория планеты ${{P}_{s}}$ располагается внутри траектории планеты ${{P}_{i}}$ и потому планеты ${{P}_{s}}$ являются внутренними. Выражения (30) содержат коэффициенты Лапласа (см. [2]), определяемые соотношениями

$A_{0}^{{is}} = \frac{2}{{\pi {{a}_{i}}{{\gamma }_{i}}}}\int\limits_0^\pi {\frac{{d\lambda }}{{{{{(1 + \alpha _{{is}}^{2} - 2{{\alpha }_{{is}}}\cos \lambda )}}^{{1/2}}}}}} ,$

$B_{p}^{{is}} = \frac{{2{{a}_{s}}{{\gamma }_{s}}}}{{\pi {{{({{a}_{i}}{{\gamma }_{i}})}}^{2}}}}\int\limits_0^\pi {\frac{{\cos (p\lambda )d\lambda }}{{{{{(1 + \alpha _{{is}}^{2} - 2{{\alpha }_{{is}}}\cos \lambda )}}^{{3/2}}}}}} ,$

$C_{p}^{{is}} = \frac{{2({{a}_{s}}{{\gamma }_{s}}{{)}^{2}}}}{{\pi {{{({{a}_{i}}{{\gamma }_{i}})}}^{3}}}}\int\limits_0^\pi {\frac{{\cos (p\lambda )d\lambda }}{{{{{(1 + \alpha _{{is}}^{2} - 2{{\alpha }_{{is}}}\cos \lambda )}}^{{5/2}}}}}} ,$

где $p = 0,\;1,\;2,\;3$. Заметим, что с помощью системы Wolfram Mathematica приведенные выше интегралы вычисляются точно и выражаются через эллиптические интегралы первого и второго рода.

Если траектория планеты ${{P}_{i}}$ располагается внутри траекторий планет ${{P}_{k}}$, ($i < k \leqslant n$), то планеты ${{P}_{k}}$ называются внешними. Обозначения, используемые в (29) для внешних планет, аналогичны обозначениям (30):

$\Pi _{{ii}}^{{ik}} = - \frac{{3{{\alpha }_{{ik}}}}}{4}B_{0}^{{ik}} - \frac{1}{2}B_{1}^{{ik}} + \frac{{15 + 6\alpha _{{ik}}^{2}}}{8}C_{0}^{{ik}} - $

$ - \;\frac{{3{{\alpha }_{{ik}}}}}{2}C_{1}^{{ik}} - \frac{9}{8}C_{2}^{{ik}},$

$\Pi _{{ik}}^{{ik}} = \frac{1}{8}(9B_{0}^{{ik}} + B_{2}^{{ik}}) - \frac{{9(1 + \alpha _{{ik}}^{2})}}{{8{{\alpha }_{{ik}}}}}C_{0}^{{ik}} + $

$ + \;\frac{{21}}{{16}}C_{1}^{{ik}} + \frac{{3(1 + \alpha _{{ik}}^{2})}}{{8{{\alpha }_{{ik}}}}}C_{2}^{{ik}} + \frac{3}{{16}}C_{3}^{{ik}},$

$\Pi _{{kk}}^{{ik}} = - \frac{3}{{4{{\alpha }_{{ik}}}}}B_{0}^{{ik}} - \frac{1}{2}B_{1}^{{ik}} + \frac{{15\alpha _{{ik}}^{2} + 6}}{{8\alpha _{{ik}}^{2}}}C_{0}^{{ik}} - $

(31)

$ - \;\frac{3}{{2{{\alpha }_{{ik}}}}}C_{1}^{{ik}} - \frac{9}{8}C_{2}^{{ik}},$

где

${{\alpha }_{{ik}}} = \frac{{{{\gamma }_{i}}{{a}_{i}}}}{{{{\gamma }_{k}}{{a}_{k}}}} < 1,$

и коэффициенты Лапласа

$A_{0}^{{ik}} = \frac{2}{{\pi {{a}_{k}}{{\gamma }_{k}}}}\int\limits_0^\pi {\frac{{d\lambda }}{{{{{(1 + \alpha _{{ik}}^{2} - 2{{\alpha }_{{ik}}}\cos \lambda )}}^{{1/2}}}}}} ,$

$B_{p}^{{ik}} = \frac{{2{{a}_{i}}{{\gamma }_{i}}}}{{\pi {{{({{a}_{k}}{{\gamma }_{k}})}}^{2}}}}\int\limits_0^\pi {\frac{{\cos (p\lambda )d\lambda }}{{{{{(1 + \alpha _{{ik}}^{2} - 2{{\alpha }_{{ik}}}\cos \lambda )}}^{{3/2}}}}}} ,$

$C_{p}^{{ik}} = \frac{{2({{a}_{i}}{{\gamma }_{i}}{{)}^{2}}}}{{\pi {{{({{a}_{k}}{{\gamma }_{k}})}}^{3}}}}\int\limits_0^\pi {\frac{{\cos (p\lambda )d\lambda }}{{{{{(1 + \alpha _{{ik}}^{2} - 2{{\alpha }_{{ik}}}\cos \lambda )}}^{{5/2}}}}}} .$

Легко видеть, что выражения (30) и (31) отличаются только определением параметров ${{\alpha }_{{is}}}$ и ${{\alpha }_{{ik}}}$, которые должны быть меньше единицы. Отметим, что оба параметра зависят от времени.

5. ЭВОЛЮЦИОННЫЕ УРАВНЕНИЯ

Эволюционные уравнения, определяющие поведение орбитальных параметров на больших интервалах времени, получаются из уравнений движения (15), если вместо возмущающих функций W_i подставить их усредненные разложения (29). Поскольку $W_{i}^{{(sec)}}$ не зависят от средних долгот ${{\lambda }_{i}}$, соответствующие им канонические импульсы ${{\Lambda }_{i}}$, а также большие полуоси ${{a}_{i}} = \Lambda _{i}^{2}{\text{/}}{{\mu }_{{i0}}}$, ($i = 1,2$, ..., n), не зависят от времени.

Вековые возмущения элементов Пуанкаре ${{\eta }_{j}}$, ${{\xi }_{j}}$, ${{q}_{j}}$, ${{p}_{j}}$ определяются как решения следующей системы $4n$ дифференциальных уравнений:

${{\dot {\eta }}_{j}} = - f\sum\limits_{s = 1}^{i - 1} {{{m}_{s}}} \left( {\frac{{\Pi _{{ii}}^{{is}}}}{{{{\Lambda }_{i}}}}{{\xi }_{i}} + \frac{{\Pi _{{is}}^{{is}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{s}}} }}{{\xi }_{s}}} \right) - $

$ - \;f\sum\limits_{k = i + 1}^n {{{m}_{k}}} \left( {\frac{{\Pi _{{ii}}^{{ik}}}}{{{{\Lambda }_{i}}}}{{\xi }_{i}} + \frac{{\Pi _{{ik}}^{{ik}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{k}}} }}{{\xi }_{k}}} \right) + \frac{{3{{{\ddot {\gamma }}}_{i}}\Lambda _{i}^{3}}}{{2{{\gamma }_{i}}\mu _{{i0}}^{2}}}{{\xi }_{i}},$

${{\dot {\xi }}_{1}} = f\sum\limits_{s = 1}^{i - 1} {{{m}_{s}}} \left( {\frac{{\Pi _{{ii}}^{{is}}}}{{{{\Lambda }_{i}}}}{{\eta }_{i}} + \frac{{\Pi _{{is}}^{{is}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{s}}} }}{{\eta }_{s}}} \right) + $

$ + \;f\sum\limits_{k = i + 1}^n {{{m}_{k}}} \left( {\frac{{\Pi _{{ii}}^{{ik}}}}{{{{\Lambda }_{i}}}}{{\eta }_{i}} + \frac{{\Pi _{{ik}}^{{ik}}}}{{\sqrt {{{\Lambda }_{i}}{{\Lambda }_{k}}} }}{{\eta }_{k}}} \right) - \frac{{3{{{\ddot {\gamma }}}_{i}}\Lambda _{i}^{3}}}{{2{{\gamma }_{i}}\mu _{{i0}}^{2}}}{{\eta }_{i}},$

${{\dot {q}}_{i}} = f\sum\limits_{s = 1}^{i - 1} {{{m}_{s}}} B_{1}^{{is}}\left( {\frac{{{{p}_{i}}}}{{4{{\Lambda }_{i}}}} - \frac{{{{p}_{s}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{s}}} }}} \right) + $

$ + \;f\sum\limits_{k = i + 1}^n {{{m}_{k}}} B_{1}^{{ik}}\left( {\frac{{{{p}_{i}}}}{{4{{\Lambda }_{i}}}} - \frac{{{{p}_{k}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{k}}} }}} \right),$

${{\dot {p}}_{i}} = - f\sum\limits_{s = 1}^{i - 1} {{{m}_{s}}} B_{1}^{{is}}\left( {\frac{{{{q}_{i}}}}{{4{{\Lambda }_{i}}}} - \frac{{{{q}_{s}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{s}}} }}} \right) - $

(32)

$ - \;f\sum\limits_{k = i + 1}^n {{{m}_{k}}} B_{1}^{{ik}}\left( {\frac{{{{q}_{i}}}}{{4{{\Lambda }_{i}}}} - \frac{{{{q}_{k}}}}{{4\sqrt {{{\Lambda }_{i}}{{\Lambda }_{k}}} }}} \right).$

Отметим, что в рассматриваемом приближении правые части уравнений (32) являются линейными функциями переменных ${{\eta }_{j}}$, ${{\xi }_{j}}$, ${{q}_{j}}$, ${{p}_{j}}$ ($j = 1,\;2,\; \ldots ,\;n$). Однако массы тел ${{m}_{j}}(t)$, а также отношения масс ${{\gamma }_{j}}(t)$ и параметры ${{\alpha }_{{is}}}$, ($i,s = 1,2$, ..., $n;\;i \ne s$) являются заданными функциями времени, что приводит к зависимости от времени коэффициентов Лапласа $B_{0}^{{is}}$, $B_{1}^{{is}}$, $B_{2}^{{is}}$, $C_{0}^{{is}}$, $C_{1}^{{is}}$, $C_{2}^{{is}}$, $C_{3}^{{is}}$. Поэтому коэффициенты полученной системы из $4n$ линейных дифференциальных уравнений (32) являются довольно сложными функциями времени, и записать общее решение этой системы в символьной форме не представляется возможным. Однако эту систему можно решать численно, задавая различные законы изменения масс тел и исследуя их влияние на вековые возмущения орбитальных элементов. После решения системы (32) можно также найти численные решения дифференциальных уравнений для средних долгот ${{\lambda }_{j}}$:

${{\dot {\lambda }}_{j}} = \frac{{\mu _{{j0}}^{2}}}{{\gamma _{j}^{2}\Lambda _{j}^{3}}} - \frac{{\partial W_{j}^{{(sec)}}}}{{\partial {{\Lambda }_{j}}}}$

с учетом ${{\Lambda }_{j}} = {\text{const}}$.

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

В настоящей работе рассматривается классическая планетарная задача $n + 1$ тел переменной массы в случае, когда массы всех тел изменяются изотропно с различными скоростями, а сами тела взаимодействуют друг с другом в соответствии с законом всемирного тяготения. Предполагается, что $n$ тел ${{P}_{1}},\; \ldots ,\;{{P}_{n}}$ движутся вокруг центрального тела-звезды ${{P}_{0}}$ по квазиэллиптическим орбитам таким образом, что их орбиты не пересекаются. Поскольку дифференциальные уравнения движения являются неинтегрируемыми, проблема исследуется в рамках теории возмущений, причем в качестве нулевого приближения используется точное решение задачи двух тел с переменными массами, описывающее апериодическое движение тел по квазиконическим сечениям. Подробно описаны основные типы символьных вычислений, связанных с разложением возмущающих функций в степенные ряды, и получены выражения для возмущающих функций с точностью до второго порядка по каноническим переменным второй системы элементов Пуанкаре. Получены эволюционные уравнения, определяющие вековые возмущения элементов Пуанкаре, которые допускают численное интегрирование при заданных законах изменения масс тел.

Отметим, что все описанные символьные вычисления выполнены с использованием системы компьютерной алгебры Wolfram Mathematica.

Список литературы

Рой А.Э. Движение по орбитам. М.: Мир, 1981. 544 с.
Мюррей К., Дермотт С. Динамика Солнечной системы / Пер. с англ. под ред. И.И. Шевченко. М.: Физматлит, 2010. 588 с.
Мерман Г.А. О представлении общего решения задачи трех тел сходящимися рядами // Бюлл. ин-та теор. астрономии. 1958. Т. 6. С. 713–732.
Брумберг В.А. Аналитические алгоритмы в небесной механике. М.: Наука, 1980. 208 с.
Борунов В.П., Рябов Ю.А. Построение численно-аналитического тригонометрического решения задачи трех тел в ССВ Mathematica и Maple // Применение систем Mathematica и Maple в научных исследованиях. М.: ВЦ РАН, 2001. С. 63–77.
Прокопеня А.Н. Некоторые алгоритмы символьных вычислений в исследованиях проблем космической динамики // Программирование. 2006. Т. 32. № 2. С. 16–22.
Маркеев А.П. Точки либрации в небесной механике и космодинамике. М.: Наука, 1978. 312 с.
Абрамов С.А., Зима Е.Б., Ростовцев В.А. Компьютерная алгебра // Программирование. 1992. С. 4–25.
Prokopenya A.N. Determination of the stability boundaries for the Hamiltonian systems with periodic coeffficients // Mathematical Modelling and Analysis. 2005. V. 10. № 2. P. 191–204.
Prokopenya A.N. Computing the stability boundaries for the Lagrange triangular solutions in the elliptic restricted three-body problem // Mathematical Modelling and Analysis. 2006. V. 11. № 1. P. 95–104.
Прокопеня А.Н. Символьные вычисления в исследованиях устойчивости решений линейных систем дифференциальных уравнений с периодическими коэффициентами // Программирование. 2007. Т. 33. № 2. С. 9–16.
Прокопеня А.Н. Нормализация гамильтониана в ограниченной задаче многих тел методами компьютерной алгебры // Программирование. 2012. Т. 38. № 3. С. 65–78.
Гутник С.А., Сарычев В.А. Применение методов компьютерной алгебры для исследования стационарных движений спутника-гиростата // Программирование. 2017. Т. 43. № 2. С. 35–44.
Гутник С.А., Сарычев В.А. Применение методов компьютерной алгебры для исследования динамики системы двух связанных тел на круговой орбите // Программирование. 2019. Т. 45. № 2. С. 32–40.
Omarov T.B. Non-Stationary Dynamical Problems in Astronomy. N.Y.: Nova Science Publ., 2002. 260 p.
Bekov A.A., Omarov T.B. The theory of orbits in non-stationary stellar systems // Astronomical and Astrophysical Transactions. 2003. V. 22. № 2. P. 145–153.
Eggleton P. Evolutionary processes in binary and multiple stars. Cambridge: Cambridge University Press, 2006. 332 p.
Rapoport I., Bear E., Soker N. The future influence of six exoplanets on the envelope properties of their parent stars on the giant branches // Monthly Notices of the Royal Astronomical Society. 2021. V. 506. № 1. P. 468–472.
Berkovič L.M. Gylden-Meščerski problem // Celestial Mechanics. 1981. V. 24. P. 407–429.
Минглибаев М.Дж. Динамика гравитирующих тел с переменными массами и размерами. LAP LAMBERT Academic Publ., 2012. 224 с.
Minglibayev M.Zh., Mayemerova G.M. Evolution of the orbital-plane orientations in the two-protoplanet three-body problem with variable masses // Astronomy Reports. 2014. V. 58. № 9. P. 667–677.
Прокопеня А.Н., Минглибаев М.Дж., Маемерова Г.М. Символьные вычисления в исследованиях проблемы трех тел с переменными массами // Программирование. 2014. Т. 40. № 2. С. 51–59.
Prokopenya A.N., Minglibayev M.Zh., Beketauov B.A. Secular perturbations of quasi-elliptic orbits in the restricted three-body problem with variable masses // International Journal of Non-Linear Mechanics. 2015. V. 73. P. 58–63.
Прокопеня А.Н., Минглибаев М.Дж., Маемерова Г.М., Иманова Ж.У. Исследование ограниченной задачи трех тел с переменными массами методами компьютерной алгебры // Программирование. 2017. Т. 43. № 5. С. 18–23.
Minglibayev M.Zh., Prokopenya A.N., Mayemerova G.M., Imanova Zh.U. Three-body problem with variable masses that change anisotropically at different rates // Mathematics in Computer Science. 2017. V. 11. № 3–4. P. 383–391.
Minglibayev M.Zh., Kosherbayeva A.B. Differential equations of planetary systems // Reports of the National Academy of Sciences of the Republic of Kazakhstan. 2020. V. 2. № 330. P. 14–20.
Субботин М.Ф. Введение в теоретическую астрономию. М.: Наука, 1968. 800 с.
Wolfram S. An elementary introduction to the Wolfram Language. Wolfram Media, Inc., 2016.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Программирование

Архивы выпусков Информация о журнале Отправить рукопись в журнал