Журнал вычислительной математики и математической физики, 2023, T. 63, № 1, стр. 154-164

Исследование вековых возмущений в ограниченной задаче трех тел переменной массы с применением компьютерной алгебры

А. Т. Ибраимова^1, 2, *, М. Дж. Минглибаев^1, 2, **, А. Н. Прокопеня^3, ***

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

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

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

^* E-mail: ibraimova@aphi.kz
^** E-mail: minglibayev@gmail.com
^*** E-mail: alexander_prokopenya@sggw.edu.pl

Поступила в редакцию 04.08.2022
После доработки 04.08.2022
Принята к публикации 10.09.2022

EDN: LMGDHY
DOI: 10.31857/S004446692301009X

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

Аннотация

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

Ключевые слова: ограниченная задача трех тел, переменная масса, уравнения движения, вековые возмущения.

1. ВВЕДЕНИЕ

Наблюдения за движением небесных тел показывают, что реальные космические системы являются нестационарными, их массы, размеры, форма и ряд других физических характеристик изменяются с течением времени в процессе эволюции (см. [1]–[6]). Учет нестационарности небесных тел приводит к существенному усложнению математических моделей их движения, что требует разработки новых алгоритмов символьных вычислений для их исследования.

Целью настоящей работы является получение дифференциальных уравнений, определяющих вековые возмущения орбитальных элементов двух тел ${{P}_{1}}$, ${{P}_{2}}$ переменной массы ${{m}_{1}} = {{m}_{1}}(t)$, ${{m}_{2}} = {{m}_{2}}(t)$, движущихся вокруг нестационарного центрального тела ${{P}_{0}}$ массой ${{m}_{0}} = {{m}_{0}}(t)$, в рамках ограниченной задачи трех тел. Предполагается, что движение двух массивных тел ${{P}_{0}}$, ${{P}_{1}}$ определяется их взаимным гравитационным притяжением, а также реактивными силами, которые возникают при неизотропном изменении их масс. Третье тело ${{P}_{2}}$ пренебрежимо малой массы $\left( {{{m}_{2}} \ll {{m}_{0}},{{m}_{1}}} \right)$ не влияет на движение двух массивных тел и движется в их гравитационном поле. Хотя масса третьего тела предполагается малой, ее изменение также приводит к возникновению реактивных сил и возмущению его движения. Исследование этой проблемы связано с довольно громоздкими символьными вычислениями, которые лучше всего выполнять с применением систем компьютерной алгебры (см. [7]–[11]). В данной работе все необходимые вычисления выполнены с помощью системы Wolfram Mathematica (см. [12], [13]).

2. ОПИСАНИЕ МОДЕЛИ

Движение тел ${{P}_{1}}$, ${{P}_{2}}$ удобно рассматривать в относительной системе координат $Oxyz$, начало которой связано с центром масс тела ${{P}_{0}}$, а оси параллельны осям абсолютной инерциальной системы координат. Обозначая радиус-векторы тел через ${{{\mathbf{r}}}_{j}} = \left( {{{x}_{j}},{{y}_{j}},{{z}_{j}}} \right),$ $j = 1,~2$, и используя второй закон Ньютона, уравнения движения можно записать в виде (см. [6], [9], [14])

(1)

$\frac{{{{d}^{2}}{{{\mathbf{r}}}_{1}}}}{{d{{t}^{2}}}} + G\left( {{{m}_{0}} + {{m}_{1}}} \right)\frac{{{{{\mathbf{r}}}_{1}}}}{{r_{1}^{3}}} - \frac{{{{d}^{2}}{{\gamma }_{1}}}}{{d{{t}^{2}}}}\frac{{{{{\mathbf{r}}}_{1}}}}{{{{\gamma }_{1}}}} = {{{\mathbf{F}}}_{1}},$

(2)

$\frac{{{{d}^{2}}{{{\mathbf{r}}}_{2}}}}{{d{{t}^{2}}}} + G{{m}_{0}}\frac{{{{{\mathbf{r}}}_{2}}}}{{r_{2}^{3}}} - \frac{{{{d}^{2}}{{\gamma }_{2}}}}{{d{{t}^{2}}}}\frac{{{{{\mathbf{r}}}_{2}}}}{{{{\gamma }_{2}}}} = {{{\mathbf{F}}}_{2}}.$

Здесь $G$ – гравитационная постоянная, дважды непрерывно дифференцируемые функции ${{\gamma }_{1}}\left( t \right)$ и ${{\gamma }_{2}}(t)$ определяются соотношениями

${{\gamma }_{1}}(t) = \frac{{{{m}_{{00}}} + {{m}_{{10}}}}}{{{{m}_{0}}(t) + {{m}_{1}}(t)}},\quad {{\gamma }_{2}}(t) = \frac{{{{m}_{{00}}}}}{{{{m}_{0}}(t)}},$

где ${{m}_{{00}}} = {{m}_{0}}({{t}_{0}})$, ${{m}_{{10}}} = {{m}_{1}}({{t}_{0}})$ – значения масс тел в начальный момент времени. Силы ${{{\mathbf{F}}}_{1}}$, ${{{\mathbf{F}}}_{2}}$ в правых частях уравнений (1), (2) можно представить в виде

(3)

${{{\mathbf{F}}}_{1}} = {\mathbf{F}}_{1}^{{\left( r \right)}} - \frac{{{{{\ddot {\gamma }}}_{1}}}}{{{{\gamma }_{1}}}}{{{\mathbf{r}}}_{1}},\quad {{{\mathbf{F}}}_{2}} = G{{m}_{1}}\left( {\frac{{{{{\mathbf{r}}}_{1}} - {{{\mathbf{r}}}_{2}}}}{{r_{{12}}^{3}}} - \frac{{{{{\mathbf{r}}}_{1}}}}{{r_{1}^{3}}}} \right) + {\mathbf{F}}_{2}^{{\left( r \right)}} - \frac{{{{{\ddot {\gamma }}}_{2}}}}{{{{\gamma }_{2}}}}{{{\mathbf{r}}}_{2}},$

где

${{r}_{{12}}} = \sqrt {{{{\left( {{{x}_{2}} - {{x}_{1}}} \right)}}^{2}} + {{{\left( {{{y}_{2}} - {{y}_{1}}} \right)}}^{2}} + {{{\left( {{{z}_{2}} - {{z}_{1}}} \right)}}^{2}}} ,\quad {{r}_{1}} = \sqrt {x_{1}^{2} + y_{1}^{2} + z_{1}^{2}} ,$

а реактивные силы ${\mathbf{F}}_{1}^{{\left( r \right)}}$, ${\mathbf{F}}_{2}^{{\left( r \right)}}$ определяются выражениями (см. [15])

(4)

${\mathbf{F}}_{1}^{{\left( r \right)}} = \frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{{\mathbf{V}}}_{1}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{{\mathbf{V}}}_{0}}~,\quad {\mathbf{F}}_{2}^{{\left( r \right)}} = \frac{{{{{\dot {m}}}_{2}}}}{{{{m}_{2}}}}{{{\mathbf{V}}}_{2}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{{\mathbf{V}}}_{0}}.$

Точка над символом в выражениях (3), (4) означает полную производную соответствующей функции по времени, а V_j, $j = 0,1,2,$ – относительные скорости частиц, покидающих тело ${{P}_{j}}$ или осаждающихся на нем, которые являются заданными функциями времени. Следует отметить, что слагаемые $\frac{{{{d}^{2}}{{\gamma }_{i}}}}{{d{{t}^{2}}}}\frac{{{{{\mathbf{r}}}_{i}}}}{{{{\gamma }_{i}}}}$ в левой части уравнений (1), (2), а также в выражениях (3) для сил ${{{\mathbf{F}}}_{1}}$, ${{{\mathbf{F}}}_{2}}$, вводятся для того, чтобы получаемые при ${{{\mathbf{F}}}_{1}} = 0$, ${{{\mathbf{F}}}_{2}} = 0$ уравнения были интегрируемы и можно было записать их общее решение при произвольных законах изменения масс тел. Такое решение далее используется в качестве первого приближения при построении возмущенных решений уравнений (1), (2). Законы изменения масс тел во времени ${{m}_{j}}(t)$ определяются на основе наблюдений и также считаются известными, причем в общем случае массы изменяются в различных темпах, т.е. справедливо соотношение

$\frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}} \ne \frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}} \ne \frac{{{{{\dot {m}}}_{2}}}}{{{{m}_{2}}}}.$

Поэтому далее будем считать, что реактивные силы ${\mathbf{F}}_{1}^{{\left( r \right)}},~{\mathbf{F}}_{2}^{{\left( r \right)}}$ являются заданными функциями времени.

Поскольку получить общее решение уравнений (1), (2) не представляется возможным даже в случае постоянных масс, воспользуемся теорией возмущений. Заметим, что при ${{{\mathbf{F}}}_{1}} = 0$, ${{{\mathbf{F}}}_{2}} = 0$ уравнения (1), (2) являются независимыми, причем каждое из них имеет точное решение, которое описывает апериодическое движение по квазиконическому сечению и может быть представлено в виде

${{x}_{j}} = {{\gamma }_{j}}{{\rho }_{j}}\left( {\cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\cos {{{{\Omega }}}_{j}} - \sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{{{\Omega }}}_{j}}\cos {{i}_{j}}} \right),$

(5)

${{y}_{j}} = {{\gamma }_{j}}{{\rho }_{j}}\left( {\cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{{{\Omega }}}_{j}} + \sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\cos {{{{\Omega }}}_{j}}\cos {{i}_{j}}} \right),$

${{z}_{j}} = {{\gamma }_{j}}{{\rho }_{j}}~\sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{i}_{j}},\quad j = 1,2,~$

где ${{\nu }_{j}}$ – истинная аномалия,

(6)

${{\rho }_{j}} = \frac{{{{a}_{j}}\left( {1 - e_{j}^{2}} \right)}}{{1 + {{e}_{j}}{\text{\;cos}}{{\nu }_{j}}}}.$

Постоянные ${{a}_{j}},\;{{e}_{j}},\;{{i}_{j}},\;{{{{\Omega }}}_{j}},\;{{\omega }_{j}}$ в (5), (6) являются аналогами известных кеплеровских орбитальных элементов (см. [6]). Введение аналога эксцентрической аномалии ${{E}_{j}}$ согласно соотношению

(7)

$\operatorname{tg} \frac{{{{\nu }_{j}}}}{2} = \sqrt {\frac{{1 + {{e}_{j}}}}{{1 - {{e}_{j}}}}} \operatorname{tg} \frac{{{{E}_{j}}}}{2}$

приводит к известному уравнению Кеплера

(8)

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

При этом зависимости эксцентрической ${{E}_{j}}$ и средней ${{M}_{j}}$ аномалий в (8) от времени

(9)

${{E}_{j}} = {{E}_{j}}(t),\quad {{M}_{j}} = \frac{{\sqrt {{{\kappa }_{j}}} }}{{a_{j}^{{3/2}}}}\left( {{{\phi }_{j}}(t) - {{\phi }_{j}}({{\tau }_{j}})} \right),\quad {{\kappa }_{1}} = G\left( {{{m}_{{00}}} + {{m}_{{10}}}} \right),\quad {{\kappa }_{2}} = G{{m}_{{00}}}$

определяются законами изменения масс тел, поскольку функции ${{\phi }_{j}}\left( t \right)$ в (9) имеют вид

(10)

${{\phi }_{j}}(t) = \mathop \smallint \limits_{{{t}_{0}}}^t \frac{{dt}}{{\gamma _{j}^{2}(t)}}.$

Через ${{\tau }_{j}}$ в (9) обозначен аналог времени прохождения тела ${{P}_{j}}$ через перицентр.

Отметим, что при заданных орбитальных параметрах ${{a}_{j}},\;{{e}_{j}},\;{{i}_{j}},\;{{{{\Omega }}}_{j}},\;{{\omega }_{j}},\;{{\tau }_{j}}$ каждого из тел ${{P}_{1}}$, ${{P}_{2}}$, а также известных функциях ${{\gamma }_{1}}(t)$, ${{\gamma }_{2}}(t)$, которые определяются законами изменения масс всех трех тел, уравнения (9), (10) позволяют найти средние аномалии ${{M}_{j}}$ как функции времени. Решая затем уравнение Кеплера (8), находим эксцентрические аномалии ${{E}_{j}}$ и, используя соотношение (7), вычисляем истинные аномалии ${{\nu }_{j}}$ и декартовы координаты тел (5). Таким образом, соотношения (5)–(10) позволяют вычислить относительные декартовы координаты тел ${{P}_{1}}$, ${{P}_{2}}$ в любой момент времени и полностью описать их невозмущенное движение. Отметим также, что при постоянных массах тел $\left( {{{\gamma }_{1}}(t) = {{\gamma }_{2}}(t) \equiv 1} \right)$ решения (5), (6) сводятся к хорошо известным решениям задачи двух тел (см. [16], [17]).

3. ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ ВОЗМУЩЕННОГО ДВИЖЕНИЯ

Как и в случае постоянных масс, гравитационное взаимодействие тел ${{P}_{1}}$ и ${{P}_{2}}$ приводит к возмущению траектории тела ${{P}_{2}}$ пренебрежимо малой массы. В отличие от классической ограниченной задачи трех тел (см. [18]), в рассматриваемом случае переменных масс учет реактивных сил приводит не только к дополнительному возмущению траектории тела ${{P}_{2}}$, но и к возмущению траектории массивного тела ${{P}_{1}}$. Это означает, что учет сил ${{\vec {F}}_{1}} \ne 0$, ${{\vec {F}}_{2}} \ne 0$ в правых частях уравнений (1), (2) приводит к зависимости орбитальных параметров ${{a}_{j}},\;{{e}_{j}},\;{{i}_{j}},\;{{{{\Omega }}}_{j}},\;{{\omega }_{j}},\;{{\tau }_{j}}$ от времени. Используя метод вариации постоянных, решения уравнений движения (1), (2) можно представить в общем виде

${{x}_{j}} = {{x}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{i}_{j}},{{{{\Omega }}}_{j}},{{\omega }_{j}},{{M}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{\tau }_{j}}} \right),t} \right),$

(11)

${{y}_{j}} = {{y}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{i}_{j}},{{{{\Omega }}}_{j}},{{\omega }_{j}},{{M}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{\tau }_{j}}} \right),t} \right),$

${{z}_{j}} = {{z}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{i}_{j}},{{{{\Omega }}}_{j}},{{\omega }_{j}},{{M}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{\tau }_{j}}} \right),t} \right),$

где функции в правых частях определяются выражениями (5), в которых вместо истинных аномалий ${{\nu }_{j}}$ используются средние аномалии ${{M}_{j}}$, а орбитальные параметры ${{a}_{j}},\;{{e}_{j}},\;{{i}_{j}},\;{{{{\Omega }}}_{j}},\;{{\omega }_{j}},\;{{\tau }_{j}}$ ($j = 1,~2$) каждого из тел ${{P}_{1}}$, ${{P}_{2}}$ являются функциями времени, причем соотношения (8), (9) определяют неявную зависимость средних аномалий ${{M}_{j}}$ от параметров ${{a}_{j}},\;{{e}_{j}},\;{{\tau }_{j}}$.

Частные производные функций (11) по времени совпадают с производными по времени функций (5) в отсутствие возмущений $\left( {{{{\mathbf{F}}}_{1}} \ne 0,\;~{{{\mathbf{F}}}_{2}} \ne 0} \right)$ и определяют скорости их изменения. Используя уравнения (5)–(10), получаем

${{\dot {x}}_{j}} = \left( {\frac{{{{{\dot {\gamma }}}_{j}}}}{{{{\gamma }_{j}}}} + \frac{{{{{\dot {\rho }}}_{j}}}}{{{{\rho }_{j}}}}} \right){{x}_{j}} - {{\gamma }_{j}}{{\rho }_{j}}{{\dot {\nu }}_{j}}\left( {\sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\cos {{{{\Omega }}}_{j}} + \cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{{{\Omega }}}_{j}}\cos {{i}_{j}}} \right),$

(12)

${{\dot {y}}_{j}} = \left( {\frac{{{{{\dot {\gamma }}}_{j}}}}{{{{\gamma }_{j}}}} + \frac{{{{{\dot {\rho }}}_{j}}}}{{{{\rho }_{j}}}}} \right){{y}_{j}} - {{\gamma }_{j}}{{\rho }_{j}}{{\dot {\nu }}_{j}}\left( {\sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{{{\Omega }}}_{j}} - \cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\cos {{{{\Omega }}}_{j}}\cos {{i}_{j}}} \right),$

${{\dot {z}}_{j}} = \left( {\frac{{{{{\dot {\gamma }}}_{j}}}}{{{{\gamma }_{j}}}} + \frac{{{{{\dot {\rho }}}_{j}}}}{{{{\rho }_{j}}}}} \right){{z}_{j}} + {{\gamma }_{j}}{{\rho }_{j}}{{\dot {\nu }}_{j}}~\cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{i}_{j}},\quad j = 1,~2,~$

где

${{\dot {\nu }}_{j}} = \frac{{\sqrt {{{\kappa }_{j}}} }}{{a_{j}^{{3/2}}{{{\left( {1 - e_{j}^{2}} \right)}}^{{3/2}}}}}\frac{{{{{\left( {1 + {{e}_{j}}~{\text{cos}}{{\nu }_{j}}} \right)}}^{2}}}}{{\gamma _{j}^{2}(t)}}.$

При наличии возмущений орбитальные параметры ${{a}_{j}},\;{{e}_{j}},\;{{i}_{j}},\;{{{{\Omega }}}_{j}},\;{{\omega }_{j}},\;{{\tau }_{j}}$ ($j = 1,2$) каждого из тел ${{P}_{1}}$, ${{P}_{2}}$ в выражениях для производных (12) также являются функциями времени, формально их можно представить в виде

${{\dot {x}}_{j}} = {{\dot {x}}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{i}_{j}},{{{{\Omega }}}_{j}},{{\omega }_{j}},{{M}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{\tau }_{j}}} \right),t} \right),$

(13)

${{\dot {y}}_{j}} = {{\dot {y}}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{i}_{j}},{{{{\Omega }}}_{j}},{{\omega }_{j}},{{M}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{\tau }_{j}}} \right),t} \right),$

${{\dot {z}}_{j}} = {{\dot {z}}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{i}_{j}},{{{{\Omega }}}_{j}},{{\omega }_{j}},{{M}_{j}}\left( {{{a}_{j}},{{e}_{j}},{{\tau }_{j}}} \right),t} \right),$

где функции в правых частях определяются выражениями (12).

Для получения дифференциальных уравнений, определяющих зависимость орбитальных параметров от времени, следует подставить выражения (11), (13) в уравнения (1), (2), что приводит к шести уравнениям второго порядка относительно 12 искомых функций ${{a}_{j}}(t),\;{{e}_{j}}(t),\;{{i}_{j}}(t),\;{{{{\Omega }}}_{j}}(t),$ ${{\omega }_{j}}(t),\;{{\tau }_{j}}(t)$ ($j = 1,2$). Поскольку такая система имеет бесконечное число решений, требуется ввести дополнительные шесть уравнений для переменных ${{a}_{j}},\;{{e}_{j}},\;{{i}_{j}},\;{{{{\Omega }}}_{j}},\;{{\omega }_{j}},\;{{\tau }_{j}}$ ($j = 1,2$). Для получения таких уравнений естественным с физической точки зрения является условие равенства частных производных функций (11) по времени скоростям изменения возмущенных координат ${{\dot {x}}_{j}},\;{{\dot {y}}_{j}},\;{{\dot {z}}_{j}}$ (см. [17], [19]). Это условие означает, что в каждый момент времени выражения для координат (11) и скоростей (13) определяются точными решениями (5), (12) задачи двух тел с переменными массами, в которых орбитальные параметры являются функциями времени. В результате получаются шесть уравнений вида (ниже мы приводим только уравнения, записанные для функций ${{x}_{j}}$, заданных уравнениями (11); уравнения для функций ${{y}_{j}}$, ${{z}_{j}}$ записываются аналогичным образом)

(14)

$\frac{{\partial {{x}_{j}}}}{{\partial {{a}_{j}}}}\frac{{d{{a}_{j}}}}{{dt}} + \frac{{\partial {{x}_{j}}}}{{\partial {{e}_{j}}}}\frac{{d{{e}_{j}}}}{{dt}} + \frac{{\partial {{x}_{j}}}}{{\partial {{i}_{j}}}}\frac{{d{{i}_{j}}}}{{dt}} + \frac{{\partial {{x}_{j}}}}{{\partial {{{{\Omega }}}_{j}}}}\frac{{d{{{{\Omega }}}_{j}}}}{{dt}} + \frac{{\partial {{x}_{j}}}}{{\partial {{\omega }_{j}}}}\frac{{d{{\omega }_{j}}}}{{dt}} + \frac{{\partial {{x}_{j}}}}{{\partial {{\tau }_{j}}}}\frac{{d{{\tau }_{j}}}}{{dt}} = 0,\quad j = 1,2.$

Еще шесть уравнений получаем при подстановке решений (11) и (13) в уравнения движения (1), (2). Соответствующее уравнение для функций ${{\dot {x}}_{j}}$, заданных уравнениями (13), имеет вид (уравнения для функций ${{\dot {y}}_{j}}$, ${{\dot {z}}_{j}}$ записываются аналогичным образом)

(15)

$\frac{{\partial {{{\dot {x}}}_{j}}}}{{\partial {{a}_{j}}}}\frac{{d{{a}_{j}}}}{{dt}} + \frac{{\partial {{{\dot {x}}}_{j}}}}{{\partial {{e}_{j}}}}\frac{{d{{e}_{j}}}}{{dt}} + \frac{{\partial {{{\dot {x}}}_{j}}}}{{\partial {{i}_{j}}}}\frac{{d{{i}_{j}}}}{{dt}} + \frac{{\partial {{{\dot {x}}}_{j}}}}{{\partial {{{{\Omega }}}_{j}}}}\frac{{d{{{{\Omega }}}_{j}}}}{{dt}} + \frac{{\partial {{{\dot {x}}}_{j}}}}{{\partial {{\omega }_{j}}}}\frac{{d{{\omega }_{j}}}}{{dt}} + \frac{{\partial {{{\dot {x}}}_{j}}}}{{\partial {{\tau }_{j}}}}\frac{{d{{\tau }_{j}}}}{{dt}} = {{F}_{{jx}}},\quad j = 1,2.$

Решение системы (14), (15) позволяет найти выражения для производных орбитальных параметров $\frac{{d{{a}_{j}}}}{{dt}}$, $\frac{{d{{e}_{j}}}}{{dt}}$, $\frac{{d{{i}_{j}}}}{{dt}}$, $\frac{{d{{{{\Omega }}}_{j}}}}{{dt}}$, $\frac{{d{{\omega }_{j}}}}{{dt}}$, $\frac{{d{{\tau }_{j}}}}{{dt}}$ для каждого из тел ${{P}_{1}}$, ${{P}_{2}}$. Отметим, что реализация описанной процедуры получения дифференциальных уравнений для орбитальных параметров требует выполнения довольно громоздских символьных вычислений. Такие вычисления можно реализовать, используя, например, систему компьютерной алгебры Wolfram Mathematica (см. [12]) и ее встроенные функции для дифференцирования выражений, преобразования тригонометрических функций, решения алгебраических уравнений и упрощения символьных выражений. Выполняя соответствующие вычисления и преобразования, в результате получаем систему дифференциальных уравнений для определения зависимости орбитальных параметров от времени (для дальнейших вычислений вместо параметра ${{\tau }_{j}}$ более удобно использовать среднюю аномалию ${{M}_{j}}$):

(16)

$\frac{{d{{a}_{j}}}}{{dt}} = \frac{{2a_{j}^{{3/2}}}}{{\sqrt {{{\kappa }_{j}}} }}\frac{{{{\gamma }_{j}}(t)}}{{1 - {{e}_{j}}~{\text{cos}}{{E}_{j}}}}\left( {\sqrt {1 - e_{j}^{2}} ~{{F}_{{\tau j}}} + {{e}_{j}}\sin {{E}_{j}}{{F}_{{rj}}}} \right),$

(17)

$\frac{{d{{e}_{j}}}}{{dt}} = \frac{{\sqrt {{{a}_{j}}\left( {1 - e_{j}^{2}} \right)} }}{{\sqrt {{{\kappa }_{j}}} }}\frac{{{{\gamma }_{j}}(t)}}{{1 - {{e}_{j}}~{\text{cos}}{{E}_{j}}}}\left( {~{{F}_{{\tau j}}}\left( {2\cos {{E}_{j}} - {{e}_{j}} - {{e}_{j}}{{{\cos }}^{2}}{{E}_{j}}} \right) + \sqrt {1 - e_{j}^{2}} \sin {{E}_{j}}{{F}_{{rj}}}} \right),$

(18)

$\frac{{d{{i}_{j}}}}{{dt}} = \frac{{\sqrt {{{a}_{j}}} ~{{\gamma }_{j}}(t)}}{{\sqrt {{{\kappa }_{j}}\left( {1 - e_{j}^{2}} \right)} }}{{F}_{{nj}}}\left( {\left( {\cos {{E}_{j}} - {{e}_{j}}} \right)\cos {{\omega }_{j}} - \sqrt {1 - e_{j}^{2}} \sin {{E}_{j}}\sin {{\omega }_{j}}} \right),$

(19)

$\frac{{d~{{{{\Omega }}}_{j}}}}{{dt}} = \frac{{\sqrt {{{a}_{j}}} ~{{\gamma }_{j}}(t)}}{{\sqrt {{{\kappa }_{j}}\left( {1 - e_{j}^{2}} \right)} \sin {{i}_{j}}}}{{F}_{{nj}}}\left( {\left( {\cos {{E}_{j}} - {{e}_{j}}} \right)\sin {{\omega }_{j}} + \sqrt {1 - e_{j}^{2}} \sin {{E}_{j}}\cos {{\omega }_{j}}} \right),$

(20)

$\begin{gathered} \frac{{d{{\omega }_{j}}}}{{dt}} = \frac{{\sqrt {{{a}_{j}}} }}{{{{e}_{j}}\sqrt {{{\kappa }_{j}}} }}\frac{{{{\gamma }_{j}}(t)}}{{1 - {{e}_{j}}~{\text{cos}}{{E}_{j}}}}\left( {~\left( {2 - e_{j}^{2} - {{e}_{j}}\cos {{E}_{j}}} \right)\sin {{E}_{j}}{{F}_{{\tau j}}} - \sqrt {1 - e_{j}^{2}} \left( {\cos {{E}_{j}} - {{e}_{j}}} \right){{F}_{{rj}}}} \right) - \\ - \;\frac{{\sqrt {{{a}_{j}}} ~{{\gamma }_{j}}(t)\cot {{i}_{j}}}}{{\sqrt {{{\kappa }_{j}}\left( {1 - e_{j}^{2}} \right)} }}{{F}_{{nj}}}\left( {\left( {\cos {{E}_{j}} - {{e}_{j}}} \right)\sin {{\omega }_{j}} + \sqrt {1 - e_{j}^{2}} \sin {{E}_{j}}\cos {{\omega }_{j}}} \right), \\ \end{gathered} $

(21)

$\begin{gathered} \frac{{d{{M}_{j}}}}{{dt}} = \frac{{\sqrt {{{\kappa }_{j}}} }}{{a_{j}^{{3/2}}\gamma _{j}^{2}(t)~}} + \frac{{\sqrt {{{a}_{j}}} }}{{{{e}_{j}}\sqrt {{{\kappa }_{j}}} }}~\frac{{{{\gamma }_{j}}(t)}}{{1 - {{e}_{j}}{\text{\;cos}}{{E}_{j}}}}\left( {\sqrt {1 - e_{j}^{2}} ~{{F}_{{\tau j}}}\sin {{E}_{j}}\left( { - 2 + e_{j}^{2} + {{e}_{j}}\cos {{E}_{j}}} \right) + } \right. \\ \left. {\mathop + \limits_{_{{}}} \;{{F}_{{rj}}}\left( {\left( {1 + 3e_{j}^{2}} \right)\cos {{E}_{j}} - {{e}_{j}}\left( {3 + e_{j}^{2}\cos 2{{E}_{j}}} \right)} \right)} \right),\quad j = 1,2. \\ \end{gathered} $

Силы ${{F}_{{rj}}},\;{{F}_{{\tau j}}},\;{{F}_{{nj}}}$ в правых частях уравнений (16)–(21) представляют собой соответственно радиальные, трансверсальные и нормальные составляющие сил ${{{\mathbf{F}}}_{1}}$, ${{{\mathbf{F}}}_{2}}$, определяемых выражениями (3). Использование орбитальных систем координат тел ${{P}_{1}}$, ${{P}_{2}}$ связано с тем, что реактивные силы ${\mathbf{F}}_{1}^{{\left( r \right)}}$, ${\mathbf{F}}_{2}^{{\left( r \right)}}$ обычно определяются именно в этих системах, а направляющие косинусы единичных векторов ${{{\mathbf{e}}}_{{rj}}} = \left( {{{e}_{{xj}}},{{e}_{{yj}}},{{e}_{{zj}}}} \right)$, ${{{\mathbf{e}}}_{{\tau j}}} = \left( {{{\tau }_{{xj}}},{{\tau }_{{yj}}},{{\tau }_{{zj}}}} \right)$, ${{{\mathbf{e}}}_{{nj}}} = \left( {{{n}_{{xj}}},{{n}_{{yj}}},{{n}_{{zj}}}} \right)$ вдоль радиального, трансверсального и нормального направлений соответственно можно легко записать на основе решений (5):

${{e}_{{xj}}} = \cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\cos {{{{\Omega }}}_{j}} - \sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{{{\Omega }}}_{j}}\cos {{i}_{j}},$

(22)

${{e}_{{yj}}} = \cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{{{\Omega }}}_{j}} + \sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\cos {{{{\Omega }}}_{j}}\cos {{i}_{j}},$

${{e}_{{zj}}} = \sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{i}_{j}},$

${{\tau }_{{xj}}} = - \sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\cos {{{{\Omega }}}_{j}} - \cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{{{\Omega }}}_{j}}\cos {{i}_{j}},$

(23)

${{\tau }_{{yj}}} = - \sin \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{{{\Omega }}}_{j}} + \cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\cos {{{{\Omega }}}_{j}}\cos {{i}_{j}},$

${{\tau }_{{zj}}} = \cos \left( {{{\nu }_{j}} + {{\omega }_{j}}} \right)\sin {{i}_{j}},~~$

${{n}_{{xj}}} = \sin {{{{\Omega }}}_{j}}\sin {{i}_{j}},$

(24)

${{n}_{{yj}}} = - \cos {{{{\Omega }}}_{j}}\sin {{i}_{j}},$

${{n}_{{zj}}} = \cos {{i}_{j}},\quad j = 1,~2.~$

Используя (3), (22)–(24), получаем для первого тела

(25)

${{F}_{{r1}}} = \frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{V}_{{r1}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{r0}}}~ - \frac{{{{{\ddot {\gamma }}}_{1}}}}{{{{\gamma }_{1}}}}{{r}_{1}},\quad {{F}_{{\tau 1}}} = \frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{V}_{{\tau 1}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{\tau 0}}},\quad {{F}_{{n1}}} = \frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{V}_{{n1}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{n0}}},$

где ${{V}_{{r1}}}$, ${{V}_{{r0}}}$, ${{V}_{{\tau 1}}}$, ${{V}_{{\tau 0}}}$, ${{V}_{{n1}}}$, ${{V}_{{n0}}}$ – составляющие относительных скоростей частиц, покидающих тела ${{P}_{1}}$ и ${{P}_{0}}$ или осаждающихся на них, вдоль радиального, трансверсального и нормального направлений в орбитальной системе координат, связанной с телом ${{P}_{1}}$. Аналогично определяем составляющие относительных скоростей частиц, покидающих тело ${{P}_{2}}$ или осаждающихся на нем, вдоль радиального, трансверсального и нормального направлений ${{V}_{{r2}}}$, ${{V}_{{\tau 2}}}$, ${{V}_{{n2}}}$ в орбитальной системе координат, связанной с телом ${{P}_{2}}$. В результате радиальная, трансверсальная и нормальная составляющие силы ${{\vec {F}}_{2}}$ принимают вид

${{F}_{{r2}}} = - G{{m}_{1}}\frac{{{{r}_{2}}}}{{r_{{12}}^{3}}} + G{{m}_{1}}\left( {\frac{{{{r}_{1}}}}{{r_{{12}}^{3}}} - \frac{1}{{r_{1}^{2}}}} \right)\left( {{{{\mathbf{e}}}_{{r1}}} \cdot {{{\mathbf{e}}}_{{r2}}}} \right) - \frac{{{{{\ddot {\gamma }}}_{2}}}}{{{{\gamma }_{2}}}}{{r}_{2}} + {{Q}_{{r2}}},$

(26)

${{F}_{{\tau 2}}} = G{{m}_{1}}\left( {\frac{{{{r}_{1}}}}{{r_{{12}}^{3}}} - \frac{1}{{r_{1}^{2}}}} \right)\left( {{{{\mathbf{e}}}_{{r1}}} \cdot {{{\mathbf{e}}}_{{\tau 2}}}} \right) + {{Q}_{{\tau 2}}},$

${{F}_{{n2}}} = G{{m}_{1}}\left( {\frac{{{{r}_{1}}}}{{r_{{12}}^{3}}} - \frac{1}{{r_{1}^{2}}}} \right)\left( {{{{\mathbf{e}}}_{{r1}}} \cdot {{{\mathbf{e}}}_{{n2}}}} \right) + {{Q}_{{n2}}}.$

Поскольку относительные скорости ${{V}_{{r0}}}$, ${{V}_{{\tau 0}}}$, ${{V}_{{n0}}}$ заданы в орбитальной системе координат, связанной с телом ${{P}_{1}}$, реактивные силы, действующие на тело ${{P}_{2}}$, принимают вид

${{Q}_{{r2}}} = \frac{{{{{\dot {m}}}_{2}}}}{{{{m}_{2}}}}{{V}_{{r2}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}\left( {{{V}_{{r0}}}\left( {{{{\mathbf{e}}}_{{r1}}} \cdot {{{\mathbf{e}}}_{{r2}}}} \right) + {{V}_{{\tau 0}}}\left( {{{{\mathbf{e}}}_{{\tau 1}}} \cdot {{{\mathbf{e}}}_{{r2}}}} \right) + {{V}_{{n0}}}\left( {{{{\mathbf{e}}}_{{n1}}} \cdot {{{\mathbf{e}}}_{{r2}}}} \right)} \right),$

(27)

${{Q}_{{\tau 2}}} = \frac{{{{{\dot {m}}}_{2}}}}{{{{m}_{2}}}}{{V}_{{\tau 2}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}\left( {{{V}_{{r0}}}\left( {{{{\mathbf{e}}}_{{r1}}} \cdot {{{\mathbf{e}}}_{{\tau 2}}}} \right) + {{V}_{{\tau 0}}}\left( {{{{\mathbf{e}}}_{{\tau 1}}} \cdot {{{\mathbf{e}}}_{{\tau 2}}}} \right) + {{V}_{{n0}}}\left( {{{{\mathbf{e}}}_{{n1}}} \cdot {{{\mathbf{e}}}_{{\tau 2}}}} \right)} \right),$

${{Q}_{{n2}}} = \frac{{{{{\dot {m}}}_{2}}}}{{{{m}_{2}}}}{{V}_{{n2}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}\left( {{{V}_{{r0}}}\left( {{{{\mathbf{e}}}_{{r1}}} \cdot {{{\mathbf{e}}}_{{n2}}}} \right) + {{V}_{{\tau 0}}}\left( {{{{\mathbf{e}}}_{{\tau 1}}} \cdot {{{\mathbf{e}}}_{{n2}}}} \right) + {{V}_{{n0}}}\left( {{{{\mathbf{e}}}_{{n1}}} \cdot {{{\mathbf{e}}}_{{n2}}}} \right)} \right).$

При заданных реактивных силах и законах изменения масс тел уравнения (16)–(27) полностью определяют возмущенное движение тел ${{P}_{1}}$, ${{P}_{2}}$.

4. ЭВОЛЮЦИОННЫЕ УРАВНЕНИЯ ДЛЯ ОРБИТАЛЬНЫХ ПАРАМЕТРОВ

Уравнения возмущенного движения (16)–(21) показывают, что при условиях ${{{\mathbf{F}}}_{1}} = 0$, ${{{\mathbf{F}}}_{2}} = 0$ орбитальные элементы ${{a}_{j}},\;{{e}_{j}},\;{{i}_{j}},\;{{{{\Omega }}}_{j}},\;{{\omega }_{j}}$ остаются постоянными, а средние аномалии ${{M}_{j}}$ являются возрастающими функциями времени. При наличии возмущений орбитальные параметры изменяются с течением времени, но получить точные решения дифференциальных уравнений (16)–(21), определяющих эти изменения, не представляется возможным. Поскольку во многих практически важных случаях эксцентриситеты e_j и наклонения i_j орбит являются малыми величинами (${{e}_{j}} \ll 1$, ${{i}_{j}} \ll 1$), уравнения (16)–(21) можно упростить, заменяя их правые части соответствующими разложениями в ряды по малым параметрам. Получаемые при этом уравнения описывают динамику исследуемых систем на длительных временных интервалах и позволяют моделировать движение тел с требуемой точностью (см. [16]–[19]). Далее рассмотрим случай малых эксцентриситетов ${{e}_{j}} \ll 1$ и наклонений ${{i}_{j}} \ll 1$ и разложим правые части уравнений (16)–(21) в ряды по этим параметрам с точностью до второго порядка. Хотя такие разложения требует выполнения довольно громоздких символьных вычислений, использование систем компьютерной алгебры позволяет весьма эффективно решать эту задачу.

При ${{e}_{j}} \ll 1$ уравнение Кеплера (8) позволяет получить выражение для эксцентрической аномалии в виде сходящегося степенного ряда по ${{e}_{j}}$ (см. [17])

${{E}_{j}} = {{M}_{j}} + {{e}_{j}}\sin {{M}_{j}} + \frac{1}{2}e_{j}^{2}\sin \left( {2{{M}_{j}}} \right) + \cdots .$

Используя это разложение и встроенную функцию Series (см. [12]), получаем

(28)

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

Разлагая в ряды по малым параметрам решения (5) с учетом (28), получаем

$\begin{gathered} {{x}_{j}} = {{\gamma }_{j}}{{a}_{j}}\left( {\cos \left( {{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right) + \frac{1}{2}{{e}_{j}}\left( {\cos \left( {2{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right) - 3\cos \left( {{{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right)} \right) + } \right. \\ + \;\frac{1}{8}e_{j}^{2}\left( {\cos \left( {{{M}_{j}} - {{\omega }_{j}} - {{{{\Omega }}}_{j}}} \right) - 4\cos \left( {{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right) + 3\cos \left( {3{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right)} \right) + \\ \left. {\mathop + \limits_{_{{_{{_{{}}}}}}} \;s_{j}^{2}\left( {\cos \left( {{{M}_{j}} + {{\omega }_{j}} - {{{{\Omega }}}_{j}}} \right) - \cos \left( {{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right)} \right)} \right), \\ \end{gathered} $

(29)

${{y}_{j}} = {{\gamma }_{j}}{{a}_{j}}\left( {\sin \left( {{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right) + \frac{1}{2}{{e}_{j}}\left( {\sin \left( {2{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right) - 3\sin \left( {{{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right)} \right) - } \right.$

$\begin{gathered} - \;\frac{1}{8}e_{j}^{2}\left( {\sin \left( {{{M}_{j}} - {{\omega }_{j}} - {{{{\Omega }}}_{j}}} \right) + 4\sin \left( {{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right) - 3\sin \left( {3{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right)} \right) - \\ \left. {\mathop - \limits_{_{{_{{_{{}}}}}}} \;s_{j}^{2}\left( {\sin \left( {{{M}_{j}} + {{\omega }_{j}} - {{{{\Omega }}}_{j}}} \right) + \sin \left( {{{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}} \right)} \right)} \right), \\ {{z}_{j}} = {{\gamma }_{j}}{{\rho }_{j}}{{s}_{j}}\left( {2\sin \left( {{{M}_{j}} + {{\omega }_{j}}} \right) + {{e}_{j}}\left( {\sin \left( {2{{M}_{j}} + {{\omega }_{j}}} \right) - 3\sin {{\omega }_{j}}} \right)} \right),\quad j = 1,2, \\ \end{gathered} $

где введенный для удобства вычислений параметр ${{s}_{j}} = \sin \frac{{{{i}_{j}}}}{2}$ является малым, так как наклонения ${{i}_{j}} \ll 1$. Используя (29), находим

(30)

${{r}_{j}} = {{\gamma }_{j}}{{\rho }_{j}}\left( {1 - {{e}_{j}}\cos {{M}_{j}} + e_{j}^{2}{{{\sin }}^{2}}{{M}_{j}}} \right),\quad j = 1,2.$

Далее вычисляем разложение выражения для $r_{{12}}^{2} = r_{1}^{2} + r_{2}^{2} - 2{{{\mathbf{r}}}_{1}}{{{\mathbf{r}}}_{2}}$, которое определяет силу ${{{\mathbf{F}}}_{2}}$:

$\begin{gathered} r_{{12}}^{2} = a_{1}^{2}\gamma _{1}^{2} + a_{2}^{2}\gamma _{2}^{2} - 2{{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right) + {{e}_{1}}\left( { - 2a_{1}^{2}\gamma _{1}^{2}\cos \left( {{{\lambda }_{1}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right) - } \right. \\ \left. { - \;{{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\cos \left( {2{{\lambda }_{1}} - {{\lambda }_{2}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right) + 3{{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\cos \left( {{{\lambda }_{2}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right)} \right) + \\ \end{gathered} $

$\begin{gathered} + \;{{e}_{2}}\left( {3{{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\cos \left( {{{\lambda }_{1}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right) - {{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\cos \left( {{{\lambda }_{1}} - 2{{\lambda }_{2}} + {{\omega }_{2}} + {{{{\Omega }}}_{2}}} \right){{ - }^{{\kern 1pt} }}} \right. \\ \left. { - \;2a_{2}^{2}\gamma _{2}^{2}\cos \left( {{{\lambda }_{2}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right)} \right) + \frac{{e_{1}^{2}}}{4}\left( {2a_{1}^{2}\gamma _{1}^{2}\left( {3 - \cos \left( {2{{\lambda }_{1}} - 2{{\omega }_{1}} - 2{{{{\Omega }}}_{1}}} \right)} \right) + } \right. \\ + \;{{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\left( {4\cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right) - 3\cos \left( {3{{\lambda }_{1}} - {{\lambda }_{2}} - 2{{\omega }_{1}} - 2{{{{\Omega }}}_{1}}} \right){{ - }^{{\kern 1pt} }}} \right. \\ \end{gathered} $

(31)

$\begin{gathered} \left. {\left. {{}^{{\kern 1pt} } - \;\cos \left( {{{\lambda }_{1}} + {{\lambda }_{2}} - 2{{\omega }_{1}} - 2{{{{\Omega }}}_{1}}} \right)} \right)} \right) + \frac{{e_{2}^{2}}}{4}\left( {2a_{2}^{2}\gamma _{2}^{2}\left( {3 - \cos \left( {2{{\lambda }_{2}} - 2{{\omega }_{2}} - 2{{{{\Omega }}}_{2}}} \right)} \right) + } \right. \\ + \;{{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\left( {4\cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right) - 3\cos \left( {{{\lambda }_{1}} - 3{{\lambda }_{2}} + 2{{\omega }_{2}} + 2{{{{\Omega }}}_{2}}} \right){{ - }^{{\kern 1pt} }}} \right. \\ \end{gathered} $

$\begin{gathered} \left. {\left. {{}^{{\kern 1pt} } - \;\cos \left( {{{\lambda }_{1}} + {{\lambda }_{2}} - 2{{\omega }_{2}} - 2{{{{\Omega }}}_{2}}} \right)} \right)} \right) - \frac{{{{e}_{1}}{{e}_{2}}}}{2}{{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\left( {9\cos \left( {{{\omega }_{1}} - {{\omega }_{2}} + {{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right){{ + }^{{\kern 1pt} }}} \right. \\ + \;\cos \left( {2{{\lambda }_{1}} - 2{{\lambda }_{2}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}} + {{\omega }_{2}} + {{{{\Omega }}}_{2}}} \right) - 3\cos \left( {2{{\lambda }_{1}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right) - \\ \left. {{}^{{\kern 1pt} } - \;3\cos \left( {2{{\lambda }_{2}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right)} \right) + {{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\left( {s_{1}^{2}\left( {2\cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right) - } \right.} \right. \\ \end{gathered} $

$\begin{gathered} \left. { - \;\cos \left( {{{\lambda }_{1}} + {{\lambda }_{2}} - 2{{{{\Omega }}}_{1}}} \right)} \right) + s_{2}^{2}\left( {2\cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right) - \cos \left( {{{\lambda }_{1}} + {{\lambda }_{2}} - 2{{{{\Omega }}}_{2}}} \right)} \right) + \\ \left. {{}^{{\kern 1pt} } + \;4{{s}_{1}}{{s}_{2}}\left( {4\cos \left( {{{\lambda }_{1}} + {{\lambda }_{2}} - {{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right) - \cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}} - {{{{\Omega }}}_{1}} + {{{{\Omega }}}_{2}}} \right)} \right)} \right), \\ \end{gathered} $

где для удобства дальнейших вычислений вместо средней аномалии ${{M}_{j}}$ введена средняя долгота ${{\lambda }_{j}} = {{M}_{j}} + {{\omega }_{j}} + {{{{\Omega }}}_{j}}$ (см. [16–19]). Используя разложения (29), (30) и выражения для направляющих косинусов (22)–(24), находим скалярные произведения единичных векторов, входящих в выражения для составляющих силы ${{\vec {F}}_{2}}$ вдоль радиального, трансверсального и нормального направлений (см. (26), (27)). Например,

$\begin{gathered} \left( {{{{\mathbf{e}}}_{{r1}}} \cdot {{{\mathbf{e}}}_{{r2}}}} \right) = \cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right) - 2{{e}_{1}}\sin \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right)\sin \left( {{{\lambda }_{1}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right) + 2{{e}_{2}}\sin \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right)\sin \left( {{{\lambda }_{2}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right) + \\ + \;\frac{{e_{1}^{2}}}{4}\sin \left( {{{\lambda }_{1}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right)\left( {\sin \left( {{{\lambda }_{2}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right) - 9\sin \left( {2{{\lambda }_{1}} - {{\lambda }_{2}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right)} \right) + \\ \end{gathered} $

(32)

$ + \;\frac{{e_{2}^{2}}}{4}\sin \left( {{{\lambda }_{2}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right)\left( {\sin \left( {{{\lambda }_{1}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right) + 9\sin \left( {{{\lambda }_{1}} - 2{{\lambda }_{2}} + {{\omega }_{2}} + {{{{\Omega }}}_{2}}} \right)} \right) + $

$\begin{gathered} + \;4{{e}_{1}}{{e}_{2}}\cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right)\sin \left( {{{\lambda }_{1}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right)\sin \left( {{{\lambda }_{2}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right) + 2s_{1}^{2}\sin \left( {{{\lambda }_{1}} - {{{{\Omega }}}_{1}}} \right)\sin \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{1}}} \right) + \\ + \;2s_{2}^{2}\sin \left( {{{\lambda }_{1}} - {{{{\Omega }}}_{2}}} \right)\sin \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{2}}} \right) + 4{{s}_{1}}{{s}_{2}}\sin \left( {{{\lambda }_{1}} - {{{{\Omega }}}_{1}}} \right)\sin \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{2}}} \right), \\ \end{gathered} $

$\begin{gathered} \left( {{{{\vec {e}}}_{{\tau 1}}} \cdot {{{\vec {e}}}_{{r2}}}} \right) = - \sin \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right) - 2{{e}_{1}}\cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right)\sin \left( {{{\lambda }_{1}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right) + 2{{e}_{2}}\cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right)\sin \left( {{{\lambda }_{2}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right) - \\ - \frac{{e_{1}^{2}}}{4}\sin \left( {{{\lambda }_{1}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right)\left( {\cos \left( {{{\lambda }_{2}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right) + 9\cos \left( {2{{\lambda }_{1}} - {{\lambda }_{2}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right)} \right) + \\ \end{gathered} $

(33)

$ + \;\frac{{e_{2}^{2}}}{4}\sin \left( {{{\lambda }_{2}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right)\left( {\cos \left( {{{\lambda }_{1}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right) + 9\cos \left( {{{\lambda }_{1}} - 2{{\lambda }_{2}} + {{\omega }_{2}} + {{{{\Omega }}}_{2}}} \right)} \right) - $

$\begin{gathered} - \;4{{e}_{1}}{{e}_{2}}\sin \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right)\sin \left( {{{\lambda }_{1}} - {{\omega }_{1}} - {{{{\Omega }}}_{1}}} \right)\sin \left( {{{\lambda }_{2}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right) - 2s_{1}^{2}\cos \left( {{{\lambda }_{1}} - {{{{\Omega }}}_{1}}} \right)\sin \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{1}}} \right) \\ - \;2s_{2}^{2}\cos \left( {{{\lambda }_{1}} - {{{{\Omega }}}_{2}}} \right)\sin \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{2}}} \right) + 4{{s}_{1}}{{s}_{2}}\cos \left( {{{\lambda }_{1}} - {{{{\Omega }}}_{1}}} \right)\sin \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{2}}} \right), \\ \end{gathered} $

(34)

$\begin{gathered} \left( {{{{\mathbf{e}}}_{{n1}}} \cdot {{{\mathbf{e}}}_{{r2}}}} \right) = - 2{{s}_{1}}\sin \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{1}}} \right) + 2{{s}_{2}}\sin \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{2}}} \right) - \\ - \;4{{e}_{2}}\left( {{{s}_{1}}\cos \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{1}}} \right) - {{s}_{2}}\cos \left( {{{\lambda }_{2}} - {{{{\Omega }}}_{2}}} \right)} \right)\sin \left( {{{\lambda }_{2}} - {{\omega }_{2}} - {{{{\Omega }}}_{2}}} \right). \\ \end{gathered} $

Подобным образом вычисляются и остальные скалярные произведения $\left( {{{{\mathbf{e}}}_{{r1}}} \cdot {{{\mathbf{e}}}_{{\tau 2}}}} \right)$, $\left( {{{{\mathbf{e}}}_{{\tau 1}}} \cdot {{{\mathbf{e}}}_{{\tau 2}}}} \right)$, $\left( {{{{\mathbf{e}}}_{{n1}}} \cdot {{{\mathbf{e}}}_{{\tau 2}}}} \right)$, $\left( {{{{\mathbf{e}}}_{{r1}}} \cdot {{{\mathbf{e}}}_{{n2}}}} \right)$, $\left( {{{{\mathbf{e}}}_{{\tau 1}}} \cdot {{{\mathbf{e}}}_{{n2}}}} \right)$, $\left( {{{{\mathbf{e}}}_{{n1}}} \cdot {{{\mathbf{e}}}_{{n2}}}} \right)$. В результате правые части дифференциальных уравнений (16)–(21) можно представить в виде разложений по степеням малых параметров ${{e}_{1}}$, ${{e}_{2}}$, ${{s}_{1}}$, ${{s}_{2}}$. Коэффициенты этих разложений являются периодическими функциями средних долгот λ₁, λ₂ и представляют собой рациональные выражения, в числителях которых содержатся тригонометрические функции cos(k λ_j), sin(k λ_j), cos(k λ ± nλ₂), sin(k λ₁ ± nλ₂) (k, n = 1, 2, ...). Поскольку получаемые выражения весьма громоздки, мы их не приводим. Отметим только, что при разложении выражения $r_{{12}}^{{ - 3}}$ с точностью до второго порядка по малым параметрам появляются рациональные выражения вида $1{\text{/}}{{\rho }_{0}}$, $1{\text{/}}\rho _{0}^{3}$, $1{\text{/}}\rho _{0}^{5}$, где

(35)

${{\rho }_{0}} = {{\left( {a_{1}^{2}\gamma _{1}^{2} - 2{{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}\cos \left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right) + a_{2}^{2}\gamma _{2}^{2}} \right)}^{{1/2}}}.$

Так как ρ₀ является периодической функцией переменных λ₁, λ₂, выражения $1{\text{/}}{{\rho }_{0}}$, $1{\text{/}}\rho _{0}^{3}$, $1{\text{/}}\rho _{0}^{5}$ следует заменить их разложениями в ряды Фурье вида

$\frac{1}{{{{\rho }_{0}}}} = \frac{1}{2}\mathop \sum \limits_{k = - \infty }^{ + \infty } {\kern 1pt} {{A}_{k}}\cos \left( {k\left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right)} \right),$

(36)

$\frac{{{{a}_{1}}{{a}_{2}}{{\gamma }_{1}}{{\gamma }_{2}}}}{{\rho _{0}^{3}}} = \frac{1}{2}\mathop \sum \limits_{k = - \infty }^{ + \infty } {\kern 1pt} {{B}_{k}}\cos \left( {k\left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right)} \right),$

$\frac{{a_{1}^{2}a_{2}^{2}\gamma _{1}^{2}\gamma _{2}^{2}}}{{\rho _{0}^{5}}} = \frac{1}{2}\mathop \sum \limits_{k = - \infty }^{ + \infty } {\kern 1pt} {{C}_{k}}\cos \left( {k\left( {{{\lambda }_{1}} - {{\lambda }_{2}}} \right)} \right),$

где ${{A}_{k}}$, ${{B}_{k}}$, ${{C}_{k}}$ – коэффициенты Лапласа, удовлетворяющие следующим рекуррентным соотношениям (см. [20]):

${{A}_{k}} = \frac{{2\left( {k - 1} \right)}}{{2k - 1}}\left( {\alpha + \frac{1}{\alpha }} \right){{A}_{{k - 1}}} - \frac{{2k - 3}}{{2k - 1}}{{A}_{{k - 2}}},\quad k \geqslant 2,$

${{B}_{k}} = \frac{{\left( {2k + 1} \right)\alpha \left( {1 + {{\alpha }^{2}}} \right)}}{{{{{\left( {1 - {{\alpha }^{2}}} \right)}}^{2}}}}{{A}_{k}} - 2{{\alpha }^{2}}\frac{{2k + 1}}{{{{{\left( {1 - {{\alpha }^{2}}} \right)}}^{2}}}}{{A}_{{k + 1}}},\quad k \geqslant 0,$

${{C}_{k}} = \frac{{\left( {2k + 3} \right)\alpha \left( {1 + {{\alpha }^{2}}} \right)}}{{{{{\left( {1 - {{\alpha }^{2}}} \right)}}^{2}}}}{{B}_{k}} - 2{{\alpha }^{2}}\frac{{2k - 1}}{{3{{{\left( {1 - {{\alpha }^{2}}} \right)}}^{2}}}}{{B}_{{k + 1}}},\quad k \geqslant 0.$

Все коэффициенты Лапласа можно вычислить, используя приведенные рекуррентные соотношения и следующие выражения для ${{A}_{0}}$ и ${{A}_{1}}$:

${{A}_{0}} = \frac{2}{{\pi {{a}_{2}}{{\gamma }_{2}}}}\mathop \smallint \limits_0^\pi \frac{{d\lambda }}{{{{{\left( {1 + {{\alpha }^{2}} - 2\alpha \cos \lambda } \right)}}^{{1/2}}}}} = \frac{4}{{\pi {{a}_{2}}{{\gamma }_{2}}\left( {1 + \alpha } \right)}}~K\left( {\frac{{4\alpha }}{{{{{\left( {1 + \alpha } \right)}}^{2}}}}} \right),$

${{A}_{1}} = \frac{2}{{\pi {{a}_{2}}{{\gamma }_{2}}}}\mathop \smallint \limits_0^\pi \frac{{\cos \lambda d\lambda }}{{{{{\left( {1 + {{\alpha }^{2}} - 2\alpha \cos \lambda } \right)}}^{{1/2}}}}} = \frac{2}{{\pi {{a}_{2}}{{\gamma }_{2}}\alpha \left( {1 + \alpha } \right)}}\left( {\left( {1 + {{\alpha }^{2}}} \right)~K\left( {\frac{{4\alpha }}{{{{{\left( {1 + \alpha } \right)}}^{2}}}}} \right) - {{{\left( {1 + \alpha } \right)}}^{2}}E\left( {\frac{{4\alpha }}{{{{{\left( {1 + \alpha } \right)}}^{2}}}}} \right)} \right).$

Здесь функции $K\left( {\frac{{4\alpha }}{{{{{\left( {1 + \alpha } \right)}}^{2}}}}} \right)$ и $E\left( {\frac{{4\alpha }}{{{{{\left( {1 + \alpha } \right)}}^{2}}}}} \right)$ обозначают полный эллиптический интеграл первого и второго рода соответственно, а параметр $\alpha = \frac{{{{a}_{1}}{{\gamma }_{1}}}}{{{{a}_{2}}{{\gamma }_{2}}}} < 1$ в случае, если тело ${{P}_{2}}$ является внешней планетой. Если тело ${{P}_{2}}$ является внутренней планетой, то ${{a}_{1}}{{\gamma }_{1}}$ и ${{a}_{2}}{{\gamma }_{2}}$ в выражениях для ${{A}_{0}}$ и ${{A}_{1}}$ нужно поменять местами и параметр $\alpha $ определить как $\alpha = \frac{{{{a}_{2}}{{\gamma }_{2}}}}{{{{a}_{1}}{{\gamma }_{1}}}} < 1$.

Выполняя описанные выше вычисления, получаем для правых частей уравнений (16)–(21) выражения в виде рядов по малым параметрам e₁, e₂, i₁, i₂, коэффициенты которых являются периодическими функциями средних долгот λ₁, λ₂. Наличие таких функций приводит к колебаниям орбитальных параметров и значительно усложняет получение численных решений уравнений (16)–(21) на больших временных интервалах. Поскольку нас интересует поведение орбитальных параметров именно на больших промежутках времени, при отсутствии соизмеримостей средних движений тел правые части уравнений (16)–(21) можно заменить их осредненными выражениями относительно средних долгот λ₁, λ₂, что приводит к уравнениям, определяющим вековые возмущения орбитальных параметров (см. [16], [19]–[21]).

Напомним, что осреднение функции $W\left( {{{\lambda }_{1}},{{\lambda }_{2}}} \right)$ по переменным ${{\lambda }_{1}},\;{{\lambda }_{2}}$ и переход к вековым возмущениям сводится к вычислению интеграла

(37)

${{W}^{{\left( {{\text{sec}}} \right)}}} = \frac{1}{{{{{\left( {2\pi } \right)}}^{2}}}}\mathop \smallint \limits_0^{2\pi } \,W\left( {{{\lambda }_{1}},{{\lambda }_{2}}} \right)d{{\lambda }_{1}}d{{\lambda }_{2}}.$

Поскольку мы рассматриваем ограниченную задачу трех тел, правые части уравнений (16)–(21), записанные для тела ${{P}_{1}}$, зависят только от средней долготы ${{\lambda }_{1}}$, поэтому интегрирование по ${{\lambda }_{2}}$ в выражении (37) является тривиальным. При подстановке вместо функции $W\left( {{{\lambda }_{1}},{{\lambda }_{2}}} \right)$ правых частей уравнений (16)–(21), в которых выполнены разложения по малым параметрам с учетом выражений для сил (25)–(27) и разложений (28)–(36), в случае тела ${{P}_{1}}$ получаем уравнения

(38)

$\frac{{d{{a}_{1}}}}{{dt}} = \frac{{a_{1}^{{3/2}}{{\gamma }_{1}}\left( {2 - e_{1}^{2}} \right)}}{{\sqrt {{{\kappa }_{1}}} }}\left( {\frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{V}_{{\tau 1}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{\tau 0}}}} \right),$

(39)

$\frac{{d{{e}_{1}}}}{{dt}} = - \frac{{3\sqrt {{{a}_{1}}} {{\gamma }_{1}}{{e}_{1}}}}{{2\sqrt {{{\kappa }_{1}}} }}\left( {\frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{V}_{{\tau 1}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{\tau 0}}}} \right),$

(40)

$\frac{{d{{i}_{1}}}}{{dt}} = - \frac{{3\sqrt {{{a}_{1}}} {{\gamma }_{1}}{{e}_{1}}}}{{2\sqrt {{{\kappa }_{1}}} }}\cos {{\omega }_{1}}\left( {\frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{V}_{{n1}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{n0}}}} \right),$

(41)

$\frac{{d~{{{{\Omega }}}_{1}}}}{{dt}} = - \frac{{3\sqrt {{{a}_{1}}} {{\gamma }_{1}}{{e}_{1}}}}{{2{{i}_{1}}\sqrt {{{\kappa }_{1}}} }}\sin {{\omega }_{1}}\left( {\frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{V}_{{n1}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{n0}}}} \right),$

(42)

$\frac{{d{{\omega }_{1}}}}{{dt}} = \frac{{3\sqrt {{{a}_{1}}} {{\gamma }_{1}}{{e}_{1}}}}{{2{{i}_{1}}\sqrt {{{\kappa }_{1}}} }}\sin {{\omega }_{1}}\left( {\frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{V}_{{n1}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{n0}}}} \right) + \frac{{\sqrt {{{a}_{1}}} {{\gamma }_{1}}}}{{\sqrt {{{\kappa }_{1}}} }}\left( {\frac{{{{{\dot {m}}}_{1}}}}{{{{m}_{1}}}}{{V}_{{r1}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{r0}}}} \right) - \frac{{3a_{1}^{{3/2}}{{\gamma }_{1}}{{{\ddot {\gamma }}}_{1}}}}{{~2\sqrt {{{\kappa }_{1}}} }}.$

Следует отметить, что интегрирование уравнений (16)–(21) по средним долготам приводит к утрате информации о положении тел на орбитах и позволяет предсказать только медленные изменения во времени орбитальных параметров a_j, e_j, i_j, Ω_j, ω_j. Поэтому усредненное уравнение (21) мы не приводим.

Получение подобных уравнений для орбитальных элементов тела ${{P}_{2}}$ требует выполнения существенно более громоздких символьных вычислений. Однако использование системы компьютерной алгебры Wolfram Mathematica позволяет успешно выполнить все вычисления. В результате получаем уравнения для вековых возмущений орбитальных элементов тела ${{P}_{2}}$:

(43)

$\begin{gathered} \frac{{d{{a}_{2}}}}{{dt}} = - \frac{{G{{m}_{1}}\sqrt {{{a}_{2}}} }}{{8\sqrt {{{\kappa }_{2}}} }}~\left( {4{{B}_{2}} - 3{{C}_{1}} + 3{{C}_{3}}} \right){{i}_{1}}{{i}_{2}}\sin \left( {{{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right) + \frac{{a_{2}^{{3/2}}{{\gamma }_{2}}}}{{\sqrt {{{\kappa }_{2}}} }}\frac{{{{{\dot {m}}}_{2}}}}{{{{m}_{2}}}}{{V}_{{\tau 2}}}\left( {2 - e_{2}^{2}} \right) + \frac{{G{{m}_{1}}\sqrt {{{a}_{2}}} }}{{32\sqrt {{{\kappa }_{2}}} }}{{e}_{1}}{{e}_{2}} \times \\ \times \;\sin \left( {{{\omega }_{1}} - {{\omega }_{2}} + {{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right)\left( {6\alpha \left( {4{{C}_{0}} - 6{{C}_{1}} + 4{{C}_{2}} - 10{{C}_{3}} + 5{{D}_{0}} - 5{{D}_{4}}} \right) - 24{{B}_{0}} + 60{{B}_{1}} + 8{{B}_{2}} - 60{{B}_{3}}\mathop - \limits_{_{{_{{_{{}}}}}}} } \right. \\ \left. { - \;42{{C}_{0}} - 108{{C}_{1}} + 288{{C}_{2}} + 12{{C}_{3}} - 54{{C}_{4}} - 150{{D}_{1}} + 165{{D}_{3}} + \frac{6}{\alpha }\left( {12{{C}_{0}} - 10{{C}_{1}} - 4{{C}_{2}} - 6{{C}_{3}} + 5{{D}_{0}} - 5{{D}_{4}}} \right)} \right){\kern 1pt} , \\ \end{gathered} $

(44)

$\begin{gathered} \frac{{d{{e}_{2}}}}{{dt}} = \frac{{3\sqrt {{{a}_{2}}} {{\gamma }_{2}}~}}{{2\sqrt {{{\kappa }_{2}}} }}\frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}\left( {{{V}_{{\tau 0}}}{{e}_{1}}\cos \left( {{{\omega }_{1}} - {{\omega }_{2}} + {{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right) - {{V}_{{n0}}}{{i}_{2}}\cos {{\omega }_{2}} + {{V}_{{r0}}}{{e}_{1}}\sin \left( {{{\omega }_{1}} - {{\omega }_{2}} + {{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right)\mathop + \limits_{} } \right. \\ \left. {\mathop + \limits_{} \;{{V}_{{n0}}}{{i}_{1}}\cos \left( {{{\omega }_{2}} - {{{{\Omega }}}_{1}} + {{{{\Omega }}}_{2}}} \right)} \right) - \frac{{3\sqrt {{{a}_{2}}} {{\gamma }_{2}}~}}{{2\sqrt {{{\kappa }_{2}}} }}\frac{{{{{\dot {m}}}_{2}}}}{{{{m}_{2}}}}{{V}_{{\tau 2}}}{{e}_{2}} + \frac{{G{{m}_{1}}~{{e}_{1}}}}{{16\sqrt {{{a}_{2}}{{\kappa }_{2}}} }}\sin \left( {{{\omega }_{1}} - {{\omega }_{2}} + {{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right) \times \\ \times \;\left( {6\alpha \left( {{{C}_{0}} + {{C}_{1}} + {{C}_{2}} - {{C}_{3}}} \right) - 6{{B}_{0}} - 10{{B}_{1}} + 2{{B}_{2}} - 6{{B}_{3}} + 3{{C}_{0}} - 27{{C}_{1}} + 3{{C}_{3}} - 3{{C}_{4}} + \frac{6}{\alpha }\left( {3{{C}_{0}} - {{C}_{2}}} \right)} \right), \\ \end{gathered} $

(45)

$\frac{{d{{i}_{2}}}}{{dt}} = - \frac{{G{{m}_{1}}~{{s}_{1}}}}{{4\sqrt {{{a}_{2}}{{\kappa }_{2}}} }}\left( {{{B}_{0}} + {{B}_{2}} - \frac{{2{{a}_{2}}{{\gamma }_{2}}}}{{a_{1}^{2}\gamma _{1}^{2}}}} \right)\sin \left( {{{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right) - \frac{{3\sqrt {{{a}_{2}}} {{\gamma }_{2}}{{e}_{2}}}}{{2\sqrt {{{\kappa }_{2}}} }}\cos {{\omega }_{2}}\left( {\frac{{{{{\dot {m}}}_{2}}}}{{{{m}_{2}}}}{{V}_{{n2}}} - \frac{{{{{\dot {m}}}_{0}}}}{{{{m}_{0}}}}{{V}_{{n0}}}} \right),$

(46)

$\frac{{d~{{{{\Omega }}}_{2}}}}{{dt}} = - \frac{{G{{m}_{1}}}}{{8\sqrt {{{a}_{2}}{{\kappa }_{2}}} }}\left( {{{B}_{0}} + {{B}_{2}}} \right)\left( {1 - \frac{{{{s}_{1}}}}{{{{s}_{2}}}}\cos \left( {{{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right)} \right)~,~~$

(47)

$\begin{gathered} \frac{{d{{\omega }_{2}}}}{{dt}} = - \frac{{G{{m}_{1}}}}{{8\sqrt {{{a}_{2}}{{\kappa }_{2}}} }}\left( {{{B}_{0}} + {{B}_{2}}} \right)\frac{{{{s}_{1}}}}{{{{s}_{2}}}}\cos \left( {{{{{\Omega }}}_{1}} - {{{{\Omega }}}_{2}}} \right) + \frac{{G{{m}_{1}}}}{{8\sqrt {{{a}_{2}}{{\kappa }_{2}}} }}\left( {{{B}_{0}} + 4{{B}_{1}} - 7{{B}_{2}} + 3{{C}_{0}}\mathop + \limits_{_{{_{{_{{}}}}}}} } \right. \\ \left. { + \;6{{C}_{1}} + 3{{C}_{2}} - 6{{C}_{3}} - \frac{6}{\alpha }\left( {{{B}_{0}} + 2{{C}_{1}}} \right) + \frac{6}{{{{\alpha }^{2}}}}{{C}_{0}}} \right) + \frac{{G{{m}_{1}}{{e}_{1}}}}{{16\sqrt {{{a}_{2}}{{\kappa }_{2}}} }}\left( { - 6\alpha \left( {{{C}_{0}} + {{C}_{1}} + {{C}_{2}} - {{C}_{3}}} \right)\mathop + \limits_{_{{_{{_{{}}}}}}} } \right. \\ \left. { + \;6{{B}_{0}} + 10{{B}_{1}} - 2{{B}_{2}} + 6{{B}_{3}} - 3{{C}_{0}} + 27{{C}_{1}} - 3{{C}_{3}} + 3{{C}_{4}} - \frac{6}{\alpha }\left( {3{{C}_{0}} - {{C}_{2}}} \right)} \right). \\ \end{gathered} $

ЗАКЛЮЧЕНИЕ

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

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

Omarov T.B. (Ed.) Non-Stationary Dynamical Problems in Astronomy. N.Y.: Nova Sci. Publ., 2002.
Bekov A.A., Omarov T.B. The theory of orbits in non-stationary stellar systems // Astron. Astrophys. Transact. 2013. V. 22. № 2. P. 145–153.
Черепащук А.М. Тесные двойные звезды. Ч. II. М.: Физматлит, 2013. 572 с.
Eggleton P. Evolutionary processes in binary and multiple stars. Cambridge Univ. Press, 2006. 332 p.
Luk’yanov L.G. Dynamical evolution of stellar orbits in close binary systems with conservative mass transfer // Astron. Rep. 2008. V. 52. № 8. P. 680–693.
Минглибаев М.Дж. Динамика гравитирующих тел с переменными массами и размерами. LAMBERT Acad. Publ., 2012. 229 с.
Прокопеня А.Н., Минглибаев М.Дж., Маемерова Г.М. Символьные вычисления в исследованиях проблемы трех тел с переменными массами // Программирование. 2014. Т. 40. № 2. С. 51–59.
Minglibayev M.Zh., Mayemerova G.M. Evolution of the orbital-plane orientations in the two-protoplanet three-body problem with variable masses // Astron. Rep. 2014. V. 58. № 9. P. 667–677.
Minglibayev M.Zh., Prokopenya A.N., Mayemerova G.M., Imanova Zh.U. Three-body problem with variable masses that change anisotropically at different rates // Math. Comp. Sci. 2017. V. 11. № 3–4. P. 383–391.
Прокопеня А.Н., Минглибаев М.Дж., Шомшекова С.А. Применение компьютерной алгебры в исследованиях двухпланетной задачи трех тел с переменными массами // Программирование. 2019. Т. 45. № 2. С. 58–65.
Minglibayev M., Prokopenya A., Shomshekova S. Computing perturbations in the two-planetary three-body problem with masses varying non-isotropically at different rates // Math. Comp. Sci. 2020. V. 14. № 2. P. 241–251.
Wolfram S. An Elementary Introduction to the Wolfram Language. Champaign, IL: Wolfram Media, 2015. 324 p.
Прокопеня А.Н. Решение физических задач с использованием системы Mathematica. Брест: БГТУ, 2005. 260 с.
Minglibayev M.Zh., Omarov Ch.T., Ibraimova A.T. New forms of the perturbed motion equation // Rep. Nation. Acad. Sci. Republ. Kazakhstan. 2020. V. 2(330). P. 5–13.
Мещерский И.В. Работы по механике тел переменной массы. М.: Гос. изд-во тех.-теор. лит-ры, 1952. 281 с.
Дубошин Г.Н. Небесная механика. Основные задачи и методы. М.: Наука, 1975. 799 с.
Рой А.Э. Движение по орбитам. М.: Мир, 1981. 544 с.
Себехей В. Теория орбит: ограниченная задача трех тел. М.: Наука, 1982. 656 с.
Brouwer D., Clemence G.M. Methods of Celestial Mechanics. N.Y.: Acad. Press, 1961. 601 p.
Шарлье К. Небесная механика. М.: Наука, 1966. 628 с.
Murray C.D., Dermott S.F. Solar system dynamics. Cambridge University Press, New York, 1999. 592 p.

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

Инструменты

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

Журнал вычислительной математики и математической физики

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