Космические исследования, 2019, T. 57, № 2, стр. 117-131

Кватернионные уравнения возмущенного движения искусственного спутника Земли

Ю. Н. Челноков *

Институт проблем точной механики и управления РАН
г. Саратов, Россия

* E-mail: chelnokovyuN@gmail.com

Поступила в редакцию 20.04.2017
После доработки 25.06.2018
Принята к публикации 20.09.2018

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

Аннотация

Получены кватернионные уравнения возмущенного движения искусственного спутника в гравитационном поле Земли с учетом его зональных, тессеральных и секториальных гармоник в четырехмерных переменных Кустаанхеймо–Штифеля и в модифицированных четырехмерных переменных, в которых уравнения движения спутника имеют более простую и симметричную структуру в сравнении с уравнениями движения в переменных Кустаанхеймо–Штифеля. Полученные уравнения линейны для невозмущенных кеплеровских движений, в общем случае имеют вид уравнений движения возмущенного четырехмерного осциллятора. Они, в отличие от классических уравнений, регулярны (не содержат особых точек типа сингулярности) для движения спутника в центральном гравитационном поле Земли под действием возмущающих сил, в описании которых не содержатся отрицательные степени расстояния спутника до центра Земли выше первой. В этих уравнениях основными переменными являются переменные Кустаанхеймо–Штифеля или модифицированные четырехмерные переменные, предложенные автором статьи, а также энергия движения спутника и время. Новая независимая переменная связана со временем дифференциальным соотношением, содержащим расстояние спутника до центра масс Земли. Уравнения удобны для применения методов нелинейной механики и высокоточных численных расчетов. В случае движения искусственного спутника в гравитационном поле Земли, в описании которого не учитываются тессеральные и секториальныее гармоники, но учитываются его зональные гармоники, найдены первые интегралы полученных уравнений движения спутника, предложены замены переменных и преобразования этих уравнений, позволившие получить для изучения движения спутника замкнутые системы дифференциальных уравнений меньшей размерности, в частности, системы уравнений четвертого и третьего порядков.

1. КВАТЕРНИОННЫЕ УРАВНЕНИЯ ВОЗМУЩЕННОГО ЦЕНТРАЛЬНОГО ДВИЖЕНИЯ МАТЕРИАЛЬНОЙ ТОЧКИ

Векторное дифференциальное уравнение возмущенного центрального движения материальной точки $M$ с массой $m,$ движущейся в центральном силовом поле с потенциалом $\Pi $ под действием возмущающей силы, равной геометрической сумме силы, имеюшей потенциал $\Pi *,$ и силы $m{\mathbf{p}},$ имеет вид

(1.1)
$\begin{gathered} \frac{{{{d}^{2}}{\mathbf{r}}}}{{d{{t}^{2}}}} = - \frac{1}{m}\left( {\frac{{d\Pi }}{{dr}}\frac{{\mathbf{r}}}{r} + \frac{{\partial \Pi {\text{*}}}}{{\partial {\mathbf{r}}}}} \right) + {\mathbf{p}}, \\ r = \left| {\mathbf{r}} \right|,\,\,\,\,\Pi = \Pi (r),\,\,\,\,\Pi {\text{*}} = \Pi {\text{*}}(t,{\mathbf{r}}), \\ {\mathbf{p}} = {\mathbf{p}}(t,{\mathbf{r}},{{d{\mathbf{r}}} \mathord{\left/ {\vphantom {{d{\mathbf{r}}} {dt}}} \right. \kern-0em} {dt}}). \\ \end{gathered} $
Здесь ${\mathbf{r}}$ – радиус-вектор, проведенный из центра $O$ притяжения, $t$ – время, ${\mathbf{p}}$ – возмущающее ускорение, обусловленное силой $m{\mathbf{p}}.$

Уравнение невозмущенного центрального движения материальной точки получается из уравнения (1.1), если в нем положить $\Pi * = 0,$ ${\mathbf{p}} = 0.$

Введем в рассмотрение следующие системы координат: $X(O{{X}_{1}}{{X}_{2}}{{X}_{3}})$ – система координат с началом в точке $O,$ движущаяся относительно инерциальной системы координат поступательно, $Y(M{{Y}_{1}}{{Y}_{2}}{{Y}_{3}})$ – система координат с началом в точке $M$ и осью $M{{Y}_{1}},$ направленной по радиусу-вектору ${\mathbf{r}}.$ Тогда

$\begin{gathered} \frac{{\partial \Pi {\text{*}}}}{{\partial {\mathbf{r}}}} = {\text{grad}}\Pi * = \sum\limits_{i = 1}^3 \,\frac{{\partial \Pi {\text{*}}}}{{\partial {{x}_{i}}}}{{{\mathbf{x}}}_{i}},\,\,\,\,{\mathbf{p}} = \sum\limits_{i = 1}^3 \,{{p}_{i}}{{{\mathbf{x}}}_{i}}, \\ {\mathbf{r}} = \sum\limits_{i = 1}^3 \,{{x}_{i}}{{{\mathbf{x}}}_{i}} = r{{{\mathbf{y}}}_{1}},\,\,\,\,\Pi * = \Pi {\text{*}}(t,{{x}_{1}},{{x}_{2}},{{x}_{3}}). \\ \end{gathered} $
Здесь ${{{\mathbf{x}}}_{i}}$ и ${{{\mathbf{y}}}_{i}}$ – орты осей ${{X}_{i}}$ и ${{Y}_{i}};$ ${{x}_{i}}$ и ${{p}_{i}}$ – проекции векторов ${\mathbf{r}}$ и ${\mathbf{p}}$ на ось ${{X}_{i}}$ (${{x}_{i}}$ – декартовые координаты точки $M$ в системе координат $X$).

В дальнейшем уравнения движения точки будем рассматривать, записывая их во вращающейся системе координат $Y.$ Ориентацию системы координат $Y$ в основной системе координат $X$ будем задавать с помощью параметров Эйлера (Родрига–Гамильтона) ${{\lambda }_{j}}(j = 0,1,2,3),$ связанных соотношением $\lambda _{0}^{2} + \lambda _{1}^{2} + \lambda _{2}^{2} + \lambda _{3}^{2} = 1.$ Декартовые координаты ${{x}_{i}}$ точки $M$ связаны с параметрами Эйлера и расстоянием $r$ от точки до центра Земли соотношениями

(1.2)
$\begin{gathered} {{x}_{1}} = r\left( {\lambda _{0}^{2} + \lambda _{1}^{2} - \lambda _{2}^{2} - \lambda _{3}^{2}} \right), \\ {{x}_{2}} = 2r({{\lambda }_{1}}{{\lambda }_{2}} + {{\lambda }_{0}}{{\lambda }_{3}}),\,\,\,\,{{x}_{3}} = 2r({{\lambda }_{1}}{{\lambda }_{3}} - {{\lambda }_{0}}{{\lambda }_{2}}). \\ \end{gathered} $
В кватернионной записи эти соотношения принимают вид
$\begin{gathered} {{{\mathbf{r}}}_{x}} = {{x}_{1}}{\mathbf{i}} + {{x}_{2}}{\mathbf{j}} + {{x}_{3}}{\mathbf{k}} = r{\mathbf{\lambda }} \circ {\mathbf{i}} \circ {\mathbf{\tilde {\lambda }}}, \\ {\mathbf{\lambda }} = {{\lambda }_{0}} + {{\lambda }_{1}}{\mathbf{i}} + {{\lambda }_{2}}{\mathbf{j}} + {{\lambda }_{3}}{\mathbf{k}}. \\ \end{gathered} $
Здесь и далее запись вида ${{{\mathbf{a}}}_{\xi }}$ означает отображение вектора ${\mathbf{a}}$ на базис $\xi (\xi = X,Y),$ т.е. ${{{\mathbf{a}}}_{\xi }}$ – кватернион с нулевой скалярной частью, составленный из проекций вектора ${\mathbf{a}}$ на базис $\xi ;$ ${\mathbf{\lambda }}$ – кватернион ориентации системы координат $Y$ в $X,$ норма которого равна единице; верхняя волна – символ кватернионного сопряжения: ${\mathbf{\tilde {\lambda }}} = {{\lambda }_{0}} - {{\lambda }_{1}}{\mathbf{i}} - {{\lambda }_{2}}{\mathbf{j}} - {{\lambda }_{3}}{\mathbf{k}},$ ${\mathbf{i}},{\mathbf{j}},{\mathbf{k}}$ – векторные мнимые единицы Гамильтона.

Доопределим вращательное движение системы координат $Y,$ полагая произвольно задаваемую проекцию ${{\omega }_{1}}$ вектора $\omega $ абсолютной угловой скорости вращения системы координат $Y$ на ось ${{Y}_{1}}$ (т.е. на направление радиуса-вектора ${\mathbf{r}}$ точки $M$) равной нулю:

(1.3)
${{\omega }_{1}} = 2( - {{\lambda }_{1}}\mathop {\dot {\lambda }}\nolimits_0 + {{\lambda }_{0}}\mathop {\dot {\lambda }}\nolimits_1 + {{\lambda }_{3}}\mathop {\dot {\lambda }}\nolimits_2 - {{\lambda }_{2}}\mathop {\dot {\lambda }}\nolimits_3 ) = 0.$

Здесь и далее верхняя точка – символ дифференцирования по времени $t.$

Введем переменную $h,$ являющуюся энергией материальной точки при $\Pi * = 0{\text{:}}$

$h = ({1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2})m{{\text{v}}^{2}} + \Pi (r).$

Переменная $h$ удовлетворяет уравнению

$\frac{{dh}}{{dt}} = \frac{{d{\mathbf{r}}}}{{dt}} \cdot \left( {m{\mathbf{p}} - \frac{{\partial \Pi {\text{*}}}}{{\partial {\mathbf{r}}}}} \right),\,\,\,\,\frac{{d{\mathbf{r}}}}{{dt}} = {\mathbf{\dot {r}}} = {\mathbf{v}}.$
Здесь и далее центральная точка – символ скалярного произведения векторов, ${\mathbf{v}}$ – вектор скорости точки в основной системе координат $X,$ $\text{v} = \left| {\mathbf{v}} \right|.$

Введем также полную энергию $h{\text{*}}$ материальной точки, определяемую соотношением

$h* = h + \Pi {\text{*}}(t,{\mathbf{r}})$
и меняющуюся в соответствии с дифференциальным уравнением

$\frac{{dh{\text{*}}}}{{dt}} = \frac{{\partial \Pi {\text{*}}}}{{\partial t}} + m{\mathbf{\dot {r}}} \cdot {\mathbf{p}}.$

Энергии $h$ и $h{\text{*}}$ будем использовать в качестве дополнительных переменных при рассмотрении движения материальной точки $M.$

Перейдем в уравнениях движения материальной точки $M$ от декартовых координат ${{x}_{i}}(i = 1,2,3)$ к новым переменным – переменным Кустаанхеймо–Штифеля ${{u}_{j}}(j = 0,1,2,3)$ ($KS$-переменным) в соответствии с формулами [1]

(1.4)
$\begin{gathered} {{x}_{1}} = u_{0}^{2} + u_{1}^{2} - u_{2}^{2} - u_{3}^{2}, \\ {{x}_{2}} = 2({{u}_{1}}{{u}_{2}} - {{u}_{0}}{{u}_{3}}),\,\,\,\,{{x}_{3}} = 2({{u}_{1}}{{u}_{3}} + {{u}_{0}}{{u}_{2}}). \\ \end{gathered} $

Из соотношений (1.2) и (1.4) следует, что переменные Кустаанхеймо–Штифеля ${{u}_{j}}$ и параметры Эйлера ${{\lambda }_{j}}$ связаны формулами

${{u}_{0}} = {{r}^{{1/2}}}{{\lambda }_{0}},\,\,\,\,{{u}_{i}} = - {{r}^{{1/2}}}{{\lambda }_{i}},\,\,\,\,i = 1,2,3.$

Таким образом, $KS$-переменные ${{u}_{j}}$ получаются из параметров Эйлера ${{\lambda }_{j}},$ характеризующих ориентацию введенной вращающейся системы координат $Y$ в основной системе координат $X,$ их умножением (со своими знаками) на множитель, равный квадратному корню из расстояния $r$ материальной точки $M$ до притягивающего центра, т.е. они являются нормированными определенным образом параметрами Эйлера. Поэтому переход от декартовых координат ${{x}_{i}}$ материальной точки к $KS$-переменным ${{u}_{j}}$ фактически означает запись исходных уравнений движения материальной точки (1.1) во вращающейся системе координат $Y.$

В кватернионной записи формулы (1.4) принимают следующий вид:

(1.5)
${{{\mathbf{r}}}_{x}} = {{x}_{1}}{\mathbf{i}} + {{x}_{2}}{\mathbf{j}} + {{x}_{3}}{\mathbf{k}} = {\mathbf{\tilde {u}}} \circ {\mathbf{i}} \circ {\mathbf{u}},$
где новая кватернионная переменная u = u0 + u1i + $ + \,{{u}_{2}}{\mathbf{j}} + {{u}_{3}}{\mathbf{k}}.$

Соотношение (1.3) для проекции угловой скорости ${{\omega }_{1}}$ в новых переменных принимает вид

${{\omega }_{1}} = 2{{r}^{{ - 1}}}({{u}_{1}}{{\dot {u}}_{0}} - {{u}_{0}}{{\dot {u}}_{1}} + {{u}_{3}}{{\dot {u}}_{2}} - {{u}_{2}}{{\dot {u}}_{3}}) = 0.$

Из этого соотношения при замене независимой переменной $t$ на новую независимую переменную $\tau $ по формуле $dt = rd\tau $ следует билинейное соотношение Е. Штифеля и Г. Шейфеле

(1.6)
${{{{u}_{1}}d{{u}_{0}}} \mathord{\left/ {\vphantom {{{{u}_{1}}d{{u}_{0}}} {d\tau }}} \right. \kern-0em} {d\tau }} - {{{{u}_{0}}d{{u}_{1}}} \mathord{\left/ {\vphantom {{{{u}_{0}}d{{u}_{1}}} {d\tau }}} \right. \kern-0em} {d\tau }} + {{{{u}_{3}}d{{u}_{2}}} \mathord{\left/ {\vphantom {{{{u}_{3}}d{{u}_{2}}} {d\tau }}} \right. \kern-0em} {d\tau }} - {{{{u}_{2}}d{{u}_{3}}} \mathord{\left/ {\vphantom {{{{u}_{2}}d{{u}_{3}}} {d\tau }}} \right. \kern-0em} {d\tau }} = 0,$
связывающее между собой переменные Кустаанхеймо–Штифеля и их первые производные и играющее, по словам Е. Штифеля и Г. Шейфеле, основную роль в их построении регулярной небесной механики [1, с. 29].

Таким образом, регуляризующему преобразованию декартовых координат Кустаанхеймо–Штифеля (1.4) и их билинейному соотношению (1.6) можно дать следующую кинематическую интерпретацию [2, 3]: регуляризующее преобразование координат Кустаанхеймо–Штифеля означает переход от декартовых координат материальной точки $M$ в основной системе координат к новым переменным, являющимся нормированными определенным образом компонентами сопряженного кватерниона поворота системы координат $Y,$ ось ${{Y}_{1}}$ которой направлена вдоль радиуса-вектора ${\mathbf{r}}$ точки $M$. Нормирующий множитель равен квадратному корню из расстояния $r$ от точки $M$ до центра притяжения. Билинейное соотношение Кустаанхеймо–Штифеля накладывает на движение трехгранника $Y$ дополнительное (неголономное) условие, заключающееся в равенстве нулю проекции ${{\omega }_{1}}$ вектора абсолютной угловой скорости трехгранника $Y$ на направление радиуса-вектора ${\mathbf{r}}$ (ось ${{Y}_{1}}$).

В результате преобразований векторного дифференциального уравнения (1.1) возмущенного центрального движения материальной точки, включающих переход от декартовых координат ${{x}_{i}}$ этой точки к переменным Кустаанхеймо–Штифеля ${{u}_{j}}$ в соответствии с формулами (1.4) или (1.5), введение в качестве дополнительных переменных энергии $h$ или $h{\text{*}}$, а также времени $t$ и использование вместо времени $t$ новой независимой переменной $\tau ,$ вводимой в соответствии с дифференциальным соотношением $dt = rd\tau ,$ получаем следующие системы кватернионных дифференциальных уравнений возмущенного центрального движения материальной точки в переменных Кустаанхеймо–Штифеля ${{u}_{j}}$ [48]:

1) Система уравнений для произвольного вида потенциала $\Pi (r)$ центрального силового поля, содержащая в качестве дополнительной переменной энергию $h = (1{\text{/}}2)m{{\text{v}}^{2}} + \Pi (r){\text{:}}$

$\begin{gathered} 2\frac{{{{d}^{2}}{\mathbf{u}}}}{{d{{\tau }^{2}}}} + \frac{1}{m}\left[ {\frac{{d(r\Pi )}}{{dr}} - h} \right]{\mathbf{u}} = r{\mathbf{Q}}, \\ \frac{{dh}}{{d\tau }} = 2m{\text{scal}}\left( {\frac{{d{\mathbf{\tilde {u}}}}}{{d\tau }} \circ {\mathbf{Q}}} \right),\,\,\,\,\frac{{dt}}{{d\tau }} = r. \\ \end{gathered} $

Дифференциальное уравнение для расстояния $r$ от точки до центра притяжения

$\frac{{{{d}^{2}}r}}{{d{{\tau }^{2}}}} + \frac{1}{m}\left[ {\frac{{d({{r}^{2}}\Pi )}}{{dr}} - 2hr} \right] = r{\text{scal}}({\mathbf{\tilde {u}}} \circ {\mathbf{Q}}).$

2) Система уравнений для произвольного вида потенциала $\Pi (r)$ центрального силового поля, содержащая в качестве дополнительной переменной полную энергию $h* = h + \Pi *{\text{:}}$

$\begin{gathered} 2\frac{{{{d}^{2}}{\mathbf{u}}}}{{d{{\tau }^{2}}}} + \frac{1}{m}\left[ {\frac{{d(r\Pi )}}{{dr}} - h{\text{*}}} \right]{\mathbf{u}} = r{\mathbf{q}} - \frac{1}{{2m}}\frac{{\partial (r\Pi *)}}{{\partial {\mathbf{u}}}}, \\ \frac{{dh{\text{*}}}}{{d\tau }} = r\frac{{\partial \Pi {\text{*}}}}{{\partial t}} + 2m{\text{scal}}\left( {\frac{{d{\mathbf{\tilde {u}}}}}{{d\tau }} \circ {\mathbf{q}}} \right),\,\,\,\,\frac{{dt}}{{d\tau }} = r. \\ \end{gathered} $

Дифференциальное уравнение для расстояния $r$

$\frac{{{{d}^{2}}r}}{{d{{\tau }^{2}}}} + \frac{1}{m}\left[ {\frac{{d({{r}^{2}}\Pi )}}{{dr}} - 2h{\text{*}}r} \right] = r\left[ { - \frac{2}{m}\Pi {\text{*}} + {\text{scal}}({\mathbf{\tilde {u}}} \circ {\mathbf{Q}})} \right].$

3) Система уравнений для возмущенного кеплеровского движения, когда потенциал центрального силового поля $\Pi (r) = - m\mu {{r}^{{ - 1}}},$ содержащая в качестве дополнительной переменной кеплеровскую энергию $h = (1{\text{/}}2)m{{\text{v}}^{2}} - m\mu {{r}^{{ - 1}}}{\text{:}}$

(1.7)
$2\frac{{{{d}^{2}}{\mathbf{u}}}}{{d{{\tau }^{2}}}} - \frac{h}{m}{\mathbf{u}} = r{\mathbf{Q}},$
(1.8)
$\frac{{dh}}{{d\tau }} = 2m{\text{scal}}\left( {\frac{{d{\mathbf{\tilde {u}}}}}{{d\tau }} \circ {\mathbf{Q}}} \right),\,\,\,\,\frac{{dt}}{{d\tau }} = r.$

Дифференциальное уравнение для расстояния $r$

$\frac{{{{d}^{2}}r}}{{d{{\tau }^{2}}}} - \frac{{2h}}{m}r - \mu = r{\text{scal}}({\mathbf{\tilde {u}}} \circ {\mathbf{Q}}).$

4) Система уравнений для возмущенного кеплеровского движения, когда потенциал центрального силового поля $\Pi (r) = - m\mu {{r}^{{ - 1}}},$ содержащая в качестве дополнительной переменной полную энергию $h* = h + \Pi *{\text{:}}$

(1.9)
$2\frac{{{{d}^{2}}{\mathbf{u}}}}{{d{{\tau }^{2}}}} - \frac{{h{\text{*}}}}{m}{\mathbf{u}} = r{\mathbf{q}} - \frac{1}{{2m}}\frac{{\partial (r\Pi *)}}{{\partial {\mathbf{u}}}},$
(1.10)
$\frac{{dh{\text{*}}}}{{d\tau }} = r\frac{{\partial \Pi {\text{*}}}}{{\partial t}} + 2m{\text{scal}}\left( {\frac{{d{\mathbf{\tilde {u}}}}}{{d\tau }} \circ {\mathbf{q}}} \right),\,\,\,\,\frac{{dt}}{{d\tau }} = r.$

Дифференциальное уравнение для расстояния $r$

(1.11)
$\frac{{{{d}^{2}}r}}{{d{{\tau }^{2}}}} - \frac{{2h{\text{*}}}}{m}r - \mu = r\left[ { - \frac{2}{m}\Pi {\text{*}} + {\text{scal}}({\mathbf{\tilde {u}}} \circ {\mathbf{Q}})} \right].$

В приведенных уравнениях $\mu $ – произведение гравитационной постоянной на массу притягивающего тела, ${\text{scal}}()$ – скалярная часть кватернионного произведения, стоящего в круглых скобках,

$\begin{gathered} {\mathbf{Q}} = {\mathbf{q}} - \frac{1}{{2m}}\frac{{\partial \Pi {\text{*}}}}{{\partial {\mathbf{u}}}},\,\,\,\,{\mathbf{q}} = - {\mathbf{i}} \circ {\mathbf{u}} \circ {{{\mathbf{p}}}_{x}}, \\ \frac{{\partial \Pi {\text{*}}}}{{\partial {\mathbf{u}}}} = \frac{{\partial \Pi {\text{*}}}}{{\partial {{u}_{0}}}} + \frac{{\partial \Pi {\text{*}}}}{{\partial {{u}_{1}}}}{\mathbf{i}} + \frac{{\partial \Pi {\text{*}}}}{{\partial {{u}_{2}}}}{\mathbf{j}} + \frac{{\partial \Pi {\text{*}}}}{{\partial {{u}_{3}}}}{\mathbf{k}}, \\ \frac{{\partial \Pi {\text{*}}}}{{\partial {{u}_{j}}}}{\text{ = }}\sum\limits_{i = 1}^3 {\frac{{\partial \Pi {\text{*}}}}{{\partial {{x}_{i}}}}\frac{{\partial x}}{{\partial {{u}_{j}}}}} ,\,\,\,\,j = 0,1,2,3, \\ \Pi * = \Pi {\text{*}}(t,{{{\mathbf{r}}}_{x}}),\,\,\,\,{{{\mathbf{p}}}_{x}} = {{{\mathbf{p}}}_{x}}(t,{{{\mathbf{r}}}_{x}},{{{\mathbf{v}}}_{x}}); \\ r = {\mathbf{u}} \circ {\mathbf{\tilde {u}}} = \sum\limits_{j = 0}^3 {u_{j}^{2}} ,\,\,\,\,{{{\mathbf{r}}}_{x}} = {\mathbf{\tilde {u}}} \circ {\mathbf{i}} \circ {\mathbf{u}}, \\ {{{\mathbf{v}}}_{x}} = {{{{\mathbf{\dot {r}}}}}_{x}} = {{{\dot {x}}}_{1}}{\mathbf{i}} + {{{\dot {x}}}_{2}}{\mathbf{j}} + {{{\dot {x}}}_{3}}{\mathbf{k}} = 2{\mathbf{\tilde {u}}} \circ {\mathbf{i}} \circ {\mathbf{\dot {u}}} = \\ = 2{{r}^{{ - 1}}}{\mathbf{\tilde {u}}} \circ {\mathbf{i}} \circ (d{\mathbf{u}}{\text{/}}d\tau ),\,\,\,\,h = 2mr\sum\limits_{j = 0}^3 {\mathop {\dot {u}}\nolimits_j^2 } + \Pi (r) = \\ = 2m{{r}^{{ - 1}}}\sum\limits_{j = 0}^3 {{{{\left( {\tfrac{{d{{{\mathbf{u}}}_{j}}}}{{d\tau }}} \right)}}^{2}}} + \Pi (r). \\ \end{gathered} $

Каждая из приведенных систем уравнений образует замкнутую систему нелинейных дифференциальных уравнений десятого порядка относительно переменных Кустаанхеймо–Штифеля ${{u}_{j}}(j = 0,1,2,3),$ их первых производных ${{d{{u}_{j}}} \mathord{\left/ {\vphantom {{d{{u}_{j}}} {d\tau }}} \right. \kern-0em} {d\tau }},$ энергии $h$ или $h{\text{*}}$ и времени $t$ (приведенное дифференциальное уравнение для расстояния $r$ не входит в состав этих систем и может рассматриваться дополнительно). В качестве новой независимой переменной выступает переменная $\tau ,$ связанная со временем $t$ дифференциальным соотношением $dt/d\tau = r.$

Отметим следующие основные достоинства уравнений движения материальной точки (1.7), (1.8) и (1.9), (1.10) в переменных Кустаанхеймо–Штифеля [116]:

– они, в отличие от традиционно используемых уравнений в декартовых и криволинейных координатах, регулярны (не содержат особых точек типа сингулярности) для возмущенного движения точки в ньютоновском центральном поле сил тяготения (при условии, что в описании действующих возмущающих сил не содержатся отрицательные степени расстояния спутника до центра Земли выше первой);

– линейны для невозмущенных кеплеровских движений и имеют в этом случае вид

${{{{d}^{2}}{{u}_{j}}} \mathord{\left/ {\vphantom {{{{d}^{2}}{{u}_{j}}} {d{{\tau }^{2}}}}} \right. \kern-0em} {d{{\tau }^{2}}}} - ({h \mathord{\left/ {\vphantom {h {(2m)}}} \right. \kern-0em} {(2m)}}){{u}_{j}} = 0,\,\,\,\,h = {\text{const}}\,\,(j = 0,1,2,3)$

(для эллиптического кеплеровского движения, когда кеплеровская энергия $h < 0,$ эти уравнения эквивалентны уравнениям движения четырехмерного одночастотного гармонического осциллятора, квадрат частоты которого равен $ - (h{\text{/}}(2m)) = {\text{const}}$);

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

– близки к линейным уравнениям для возмущенных кеплеровских движений;

– позволяют во многих случаях представить правые части дифференциальных уравнений движения небесных и космических тел в полиномиальной форме, удобной для их решения с помощью ЭВМ.

Эти свойства уравнений в переменных Кустаанхеймо–Штифеля позволили разработать эффективные методы нахождения решений в аналитической или численной форме таких трудных для классических методов задач как исследование движения вблизи притягивающих масс или движения по орбитам с большими эксцентриситетами. Так, Е. Штифелем, Г. Шейфеле, Т.В. Бордовицыной и др. [1, 11, 15] показано, что использование регулярных уравнений в переменных Кустаанхеймо–Штифеля позволяет повысить точность численного решения ряда задач небесной механики и астродинамики (например, задачи о движении ИСЗ по орбитам с большими эксцентриситетами) от трех до пяти порядков по сравнению с решениями, полученными при использовании классических (ньютоновских) уравнений. Кроме того, кватернионные уравнения в переменных Кустаанхеймо–Штифеля позволили построить (с использованием принципа максимума Понтрягина) эффективные решения ряда задач оптимального управления орбитальным движением космического аппарата (Я.Г. Сапунков, Ю.Н. Челноков, В.А. Юрко [7, 8, 1722]).

2. УРАВНЕНИЯ ВОЗМУЩЕННОГО ДВИЖЕНИЯ СПУТНИКА В ГРАВИТАЦИОННОМ ПОЛЕ ЗЕМЛИ В ПЕРЕМЕННЫХ КУСТААНХЕЙМО–ШТИФЕЛЯ

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

$\frac{{{{d}^{2}}{\mathbf{r}}}}{{d{{t}^{2}}}} = - \frac{{\partial {{\Pi }_{{\text{E}}}}}}{{\partial {\mathbf{r}}}} + {\mathbf{p}} = - \left( {\frac{{d\Pi }}{{dr}}\frac{{\mathbf{r}}}{r} + \frac{{\partial \Pi {\text{*}}}}{{\partial {\mathbf{r}}}}} \right) + {\mathbf{p}},$
где
$\begin{gathered} r = \left| {\mathbf{r}} \right|,\,\,\,\,{{\Pi }_{{\text{E}}}} = \Pi + \Pi *,\,\,\,\,\Pi = \Pi (r) = - f{{m}_{{\text{E}}}}{{r}^{{ - 1}}}, \\ \Pi * = \Pi {\text{*}}(t,{\mathbf{r}}),\,\,\,\,\,{\mathbf{p}} = {\mathbf{p}}(t,{\mathbf{r}},d{\mathbf{r}}{\text{/}}dt), \\ \end{gathered} $
${\mathbf{r}}$ – геоцентрический радиус-вектор спутника, ${{m}_{{\text{E}}}} = {{m}_{{{\text{earth}}}}}$ – масса Земли, ${{\Pi }_{{\text{E}}}} = {{\Pi }_{{{\text{earth}}}}} = \Pi + \Pi {\text{*}}$ – потенциал гравитационного поля Земли, $\Pi = \Pi (r)$ – его центральная составляющая, $\Pi * = \Pi _{z}^{ * }({\mathbf{r}}) + \Pi _{{ts}}^{ * }(t,{\mathbf{r}})$ – составляющая, обусловленная нецентральностью гравитационного поля Земли ($\Pi _{z}^{ * }({\mathbf{r}})$ – составляющая потенциала, содержащая зональные гармоники гравитационного поля Земли, $\Pi _{{ts}}^{ * }(t,{\mathbf{r}})$ – составляющая потенциала, содержащая тессеральные и секториальные гармоники гравитационного поля Земли [2325]), $f$ – постоянная тяготения, ${\mathbf{p}}$ – возмущающее ускорение центра масс спутника от действующих на спутник негравитационных сил.

Кватернионные дифференциальные уравнения возмущенного движения искусственного спутника в гравитационном поле Земли в переменных Кустаанхеймо–Штифеля получаются из уравнений (1.9), (1.10) и имеют следующий вид:

(2.1)
$2\frac{{{{d}^{2}}{\mathbf{u}}}}{{d{{\tau }^{2}}}} - h{\text{*}}{\mathbf{u}} = r{\mathbf{q}} - \frac{1}{2}\frac{{\partial \left( {r\Pi {\text{*}}} \right)}}{{\partial {\mathbf{u}}}},$
(2.2)
$\frac{{dh{\text{*}}}}{{d\tau }} = r\frac{{\partial \Pi _{{ts}}^{ * }}}{{\partial t}} + 2{\text{scal}}\left( {\frac{{d{\mathbf{\tilde {u}}}}}{{d\tau }} \circ {\mathbf{q}}} \right),\,\,\,\,\,\frac{{dt}}{{d\tau }} = r.$

Дифференциальное уравнение для расстояния $r$ спутника до центра масс Земли имеет в соответствии с (1.11) следующий вид:

(2.3)
$\frac{{{{d}^{2}}r}}{{d{{\tau }^{2}}}} - 2h{\text{*}}r - f{{m}_{{\text{E}}}} = r\left[ { - 2\Pi {\text{*}} + {\text{scal}}({\mathbf{\tilde {u}}} \circ {\mathbf{Q}})} \right].$

В уравнениях (2.1) и (2.2) расстояние r = $ = {\mathbf{u}} \circ {\mathbf{\tilde {u}}} = u_{0}^{2} + u_{1}^{2} + u_{2}^{2} + u_{3}^{2},$ полная энергия спутника $h{\text{*}}$ определяется соотношением

$\begin{gathered} h* = h + \Pi {\text{*}}({{{\mathbf{r}}}_{x}}), \\ h = 2r\sum\limits_{j = 0}^3 \dot {u}_{j}^{2} + \Pi (r) = 2{{r}^{{ - 1}}}\sum\limits_{j = 0}^3 \mathop {\left( {\frac{{d{{{\mathbf{u}}}_{j}}}}{{d\tau }}} \right)}\nolimits^2 + \Pi (r); \\ {\mathbf{Q}} = {\mathbf{q}} - \frac{1}{2}\frac{{\partial \Pi {\text{*}}}}{{\partial {\mathbf{u}}}},\,\,\,\,{\mathbf{q}} = - {\mathbf{i}} \circ {\mathbf{u}} \circ {{{\mathbf{p}}}_{x}},\,\,\,\,\Pi * = \Pi {\text{*}}(t,{{{\mathbf{r}}}_{x}}), \\ {{{\mathbf{p}}}_{x}} = {{{\mathbf{p}}}_{x}}(t,{{{\mathbf{r}}}_{x}},{{{\mathbf{v}}}_{x}}),\,\,\,\,\,{{{\mathbf{r}}}_{x}} = {\mathbf{\tilde {u}}} \circ {\mathbf{i}} \circ {\mathbf{u}}, \\ {{{\mathbf{v}}}_{x}} = {{{{\mathbf{\dot {r}}}}}_{x}} = {{{\dot {x}}}_{1}}{\mathbf{i}} + {{{\dot {x}}}_{2}}{\mathbf{j}} + {{{\dot {x}}}_{3}}{\mathbf{k}} = 2{\mathbf{\tilde {u}}} \circ {\mathbf{i}} \circ {\mathbf{\dot {u}}} = \\ = 2{{r}^{{ - 1}}}{\mathbf{\tilde {u}}} \circ {\mathbf{i}} \circ (d{\mathbf{u}}{\text{/}}d\tau ). \\ \end{gathered} $

2.1. Уравнения возмущенного движения спутника в гравитационном поле Земли без учета его тессеральных и секториальных гармоник в переменных Кустаанхеймо–Штифеля. Отбрасывая тессеральные и секториальные гармоники гравитационного поля Земли (учитывая лишь его центральную составляющую и зональные гармоники), имеем [2325]:

$\begin{gathered} {{\Pi }_{{\text{E}}}} = {{\Pi }_{{\text{E}}}}(r,\varphi ) = \\ = - \frac{{f{{m}_{{\text{E}}}}}}{r} + \frac{{f{{m}_{{\text{E}}}}}}{r}\sum\limits_{n = 2}^\infty {{J}_{n}}\mathop {\left( {\frac{R}{r}} \right)}\nolimits^n {{P}_{n}}(sin\varphi ),\,\,\,\,r = \left| {\mathbf{r}} \right|, \\ \end{gathered} $
где $R$ – средний экваториальный радиус Земли, ${{J}_{n}}$ – безразмерные постоянные, характеризующие фигуру Земли, ${{P}_{n}}$ – полином Лежандра n-го порядка, $\varphi $ – геоцентрическая широта.

Систему координат $X,$ в которой рассматривается движение спутника, введем следующим образом: ее начало $O$ поместим в центр Земли, ось $O{{X}_{3}}$ направим к северному полюсу Земли, а ось $O{{X}_{1}}$ – в точку весеннего равнодействия.

Потенциал ${{\Pi }_{{\text{E}}}}$ представим в виде

${{\Pi }_{{\text{E}}}} = \Pi (r) + \Pi _{z}^{ * }(r,\gamma ),\,\,\,\,\gamma = sin\varphi = \cos \vartheta = {{{{x}_{3}}} \mathord{\left/ {\vphantom {{{{x}_{3}}} r}} \right. \kern-0em} r},$
(2.4)
$\Pi (r) = - \frac{{f{{m}_{{\text{E}}}}}}{r},\,\,\,\,\Pi _{z}^{ * }(r,\gamma ) = \frac{{f{{m}_{{\text{E}}}}}}{r}\sum\limits_{n = 2}^\infty {{J}_{n}}\mathop {\left( {\frac{R}{r}} \right)}\nolimits^n {{P}_{n}}(\gamma ),$
где $\vartheta $ – угол между осью $O{{X}_{3}}$ и вектором ${\mathbf{r}}.$

В соответствии с (1.2) и (1.4)

(2.5)
$\begin{gathered} \gamma = cos\vartheta = {{r}^{{ - 1}}}{{x}_{3}} = \\ = 2({{\lambda }_{1}}{{\lambda }_{3}} - {{\lambda }_{0}}{{\lambda }_{2}}) = 2{{r}^{{ - 1}}}({{u}_{1}}{{u}_{3}} + {{u}_{0}}{{u}_{2}}). \\ \end{gathered} $

Обозначим

(2.6)
$\Pi _{z}^{ + }(r,\gamma ) = r\Pi _{z}^{ * }(r,\gamma ) = f{{m}_{{\text{E}}}}\sum\limits_{n = 2}^\infty {{J}_{n}}\mathop {\left( {\frac{R}{r}} \right)}\nolimits^n {{P}_{n}}(\gamma ).$

Используя (2.5), (2.6), находим

(2.7)
$\begin{gathered} \frac{\partial }{{\partial {\mathbf{u}}}}\left( {\Pi _{z}^{ + }(r,\gamma )} \right) = \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}\frac{{\partial r}}{{\partial {\mathbf{u}}}} + \frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }}\frac{{\partial \gamma }}{{\partial {\mathbf{u}}}} = \\ = 2[( - {{l}_{1}}{{u}_{0}} + {{l}_{2}}{{u}_{2}}) + ( - {{l}_{1}}{{u}_{1}} + {{l}_{2}}{{u}_{3}}){\mathbf{i}} + \\ + \,\,( - {{l}_{1}}{{u}_{2}} + {{l}_{2}}{{u}_{0}}){\mathbf{j}} + ( - {{l}_{1}}{{u}_{3}} + {{l}_{2}}{{u}_{1}}){\mathbf{k}}], \\ \end{gathered} $
(2.8)
${{l}_{1}} = \frac{\gamma }{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}},\,\,\,\,{{l}_{2}} = \frac{1}{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }}.$

Подставляя в уравнения (2.1), (2.2) соотношения (2.6)–(2.8) (с учетом того, что в рассматриваемом случае $\Pi {\text{*}} = \Pi _{z}^{ * }$), и переходя к скалярной записи, получаем следующие дифференциальные уравнения возмущенного движения спутника в гравитационном поле Земли без учета его тессеральных и секториальных гармоник в переменных Кустаанхеймо–Штифеля:

(2.9)
$\begin{gathered} \frac{{{{d}^{2}}{{u}_{0}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{0}} = \frac{1}{2}({{l}_{1}}{{u}_{0}} - {{l}_{2}}{{u}_{2}} + r{{q}_{0}}), \\ \frac{{{{d}^{2}}{{u}_{1}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{1}} = \frac{1}{2}({{l}_{1}}{{u}_{1}} - {{l}_{2}}{{u}_{3}} + r{{q}_{1}}), \\ \frac{{{{d}^{2}}{{u}_{2}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{2}} = \frac{1}{2}({{l}_{1}}{{u}_{2}} - {{l}_{2}}{{u}_{0}} + r{{q}_{2}}), \\ \frac{{{{d}^{2}}{{u}_{3}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{3}} = \frac{1}{2}({{l}_{1}}{{u}_{3}} - {{l}_{2}}{{u}_{1}} + r{{q}_{3}}); \\ \end{gathered} $
(2.10)
$\frac{{dh{\text{*}}}}{{d\tau }} = 2\sum\limits_{j = 0}^3 \,\left( {\frac{{d{{u}_{j}}}}{{d\tau }}{{q}_{j}}} \right),\,\,\,\,\frac{{dt}}{{d\tau }} = r,$
где величины ${{l}_{1}}$ и ${{l}_{2}}$ определяются соотношениями (2.8), $r = u_{0}^{2} + u_{1}^{2} + u_{2}^{2} + u_{3}^{2},$ $\gamma = 2{{r}^{{ - 1}}}({{u}_{1}}{{u}_{3}} + {{u}_{0}}{{u}_{2}}),$ а величины ${{q}_{j}}$ – компоненты кватерниона ${\mathbf{q}} = - {\mathbf{i}} \circ {\mathbf{u}} \circ {{{\mathbf{p}}}_{x}}{\text{:}}$

(2.11)
$\begin{gathered} {{q}_{0}} = {{u}_{0}}{{p}_{1}} - {{u}_{3}}{{p}_{2}} + {{u}_{2}}{{p}_{3}}, \\ {{q}_{1}} = {{u}_{1}}{{p}_{1}} + {{u}_{2}}{{p}_{2}} + {{u}_{3}}{{p}_{3}}, \\ {{q}_{2}}\, = \, - {{u}_{2}}{{p}_{1}} + {{u}_{1}}{{p}_{2}} + {{u}_{0}}{{p}_{3}},\,\,\,{{q}_{3}} = - {{u}_{3}}{{p}_{1}} - {{u}_{0}}{{p}_{2}} + {{u}_{1}}{{p}_{3}}. \\ \end{gathered} $

В соотношениях (2.11), по-прежнему, ${{p}_{i}}$ – проекции вектора возмущающего ускорения ${\mathbf{p}}$ на оси системы координат $X.$

Дифференциальные уравнения (2.9), (2.10) для переменных Кустаанхеймо–Штифеля можно также записать в следующей форме:

$\begin{gathered} \frac{{{{d}^{2}}{{u}_{j}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{j}} = \\ = \frac{1}{2}\left( {\left( {\frac{\gamma }{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right){{u}_{j}} - \left( {\frac{1}{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }}} \right)u_{j}^{ * } + r{{q}_{j}}} \right), \\ j = 0,1,2,3; \\ u_{0}^{{\text{*}}} = {{u}_{2}},\,\,\,\,u_{1}^{{\text{*}}} = {{u}_{3}},u_{2}^{{\text{*}}} = {{u}_{0}},\,\,\,\,u_{3}^{{\text{*}}} = {{u}_{1}}. \\ \end{gathered} $

Дифференциальное уравнение (2.3) для расстояния $r$ спутника до центра масс Земли с учетом равенства ${\text{scal}}({\mathbf{\tilde {u}}} \circ {\mathbf{Q}}) = {\text{scal}}({\mathbf{\tilde {u}}} \circ {\mathbf{q}}) - r\partial \Pi _{z}^{{\text{*}}}{\text{/}}\partial r$ принимает вид

$\begin{gathered} \frac{{{{d}^{2}}r}}{{d{{\tau }^{2}}}} - 2h{\text{*}}r = f{{m}_{{\text{E}}}} - \frac{\partial }{{\partial r}}(r\Pi _{z}^{ + }) + r{\text{scal}}({\mathbf{\tilde {u}}} \circ {\mathbf{q}}), \\ {\mathbf{q}} = - {\mathbf{i}} \circ {\mathbf{u}} \circ {{{\mathbf{p}}}_{x}}. \\ \end{gathered} $

Отметим, что в случае отсутствия возмущающего ускорения ${\mathbf{p}}$ центра масс спутника от действующих на спутник негравитационных сил величины ${{q}_{j}} = 0.$ Полная энергия $h{\text{*}}$ спутника становится постоянной величиной, поэтому дифференциальное уравнение для полной энергии выпадает из рассмотрения. Уравнения движения спутника в этом случае принимают вид

(2.12)
$\begin{gathered} \frac{{{{d}^{2}}{{u}_{j}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{j}} = \frac{1}{2}\left( {\left( {\frac{\gamma }{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right){{u}_{j}} - \left( {\frac{1}{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }}} \right)u_{j}^{{\text{*}}}} \right), \\ h* = {\text{const}};\,\,\,\,j = 0,1,2,3, \\ u_{0}^{{\text{*}}} = {{u}_{2}},\,\,\,u_{1}^{{\text{*}}} = {{u}_{3}},\,\,\,\,u_{2}^{{\text{*}}} = {{u}_{0}},\,\,\,\,u_{3}^{{\text{*}}} = {{u}_{1}}. \\ \end{gathered} $

Эти уравнения необходимо дополнить уравнением для времени $t{\text{:}}$ $dt{\text{/}}d\tau = r.$

2.2. Уравнения возмущенного движения спутника в гравитационном поле Земли с учетом его тессеральных и секториальных гармоник в переменных Кустаанхеймо–Штифеля. Потенциал гравитационного поля Земли c учетом его тессеральных и секториальных гармоник имеет вид

${{\Pi }_{{\text{E}}}} = \Pi + \Pi _{z}^{ * } + \Pi _{{ts}}^{ * },$
где составляющие потенциала $\Pi $ и $\Pi _{z}^{ * }$ имеют вышеприведенный вид (2.4), а составляющая потенциала $\Pi _{{ts}}^{ * }$ описывает тессеральные и секториальные гармоники гравитационного поля Земли и имеет следующий вид [23, 24]:
$\begin{gathered} \Pi _{{ts}}^{ * }(r,\gamma ,\lambda ) = \frac{{f{{m}_{{\text{E}}}}}}{r}\sum\limits_{n = 2}^\infty \sum\limits_{k = 1}^n \mathop {\left( {\frac{R}{r}} \right)}\nolimits^n {{P}_{{nk}}}(\gamma ) \times \\ \times \,\,\left( {{{C}_{{nk}}}cos(k\lambda ) + {{S}_{{nk}}}sin(k\lambda )} \right),\,\,\,\,r = \left| {\mathbf{r}} \right|, \\ \end{gathered} $
где ${{C}_{{nk}}},$ ${{S}_{{nk}}}$ – безразмерные постоянные, характеризующие фигуру Земли, ${{P}_{{nk}}}$ – присоединенные функции Лежандра, $\lambda $ – географическая долгота.

Географическая долгота определяется через декартовые координаты и переменные Кустаанхеймо–Штифеля посредством соотношений

(2.13)
$\begin{gathered} \lambda = {{\lambda }_{a}} - {{\Omega }_{{\text{E}}}}t,\,\,\,\,\,{{\lambda }_{a}} = {\text{arctg}}\frac{{{{x}_{2}}}}{{{{x}_{1}}}} = \\ = {\text{arctg}}\frac{{2({{u}_{1}}{{u}_{2}} - {{u}_{0}}{{u}_{3}})}}{{u_{0}^{2} + u_{1}^{2} - u_{2}^{2} - u_{3}^{2}}}, \\ \end{gathered} $
где ${{\lambda }_{a}}$ – абсолютная долгота, ${{\Omega }_{{\text{E}}}}$ – угловая скорость суточного вращения Земли.

Тригонометрические функции $cos(k\lambda )$ и $sin(k\lambda )$ находятся через функции $cos\lambda $ и $sin\lambda $ с помощью кватернионных соотношений, вытекающих из формулы Муавра:

$\begin{gathered} cos(2\lambda ) + {\mathbf{i}}sin(2\lambda ) = \\ = (cos\lambda + {\mathbf{i}}sin\lambda ) \circ (cos\lambda + {\mathbf{i}}sin\lambda ) = \\ = {{(cos\lambda + {\mathbf{i}}sin\lambda )}^{2}}, \\ \end{gathered} $
$\begin{gathered} cos(3\lambda ) + {\mathbf{i}}sin(3\lambda ) = \\ = \,(cos{\kern 1pt} \lambda \, + \,{\mathbf{i}}sin{\kern 1pt} \lambda )\, \circ \,(cos{\kern 1pt} \lambda \, + \,{\mathbf{i}}sin{\kern 1pt} \lambda )\, \circ \,(cos{\kern 1pt} \lambda \, + \,{\mathbf{i}}sin{\kern 1pt} \lambda )\, = \\ = {{(cos\lambda + {\mathbf{i}}sin\lambda )}^{3}},...,cos(k\lambda ) + {\mathbf{i}}sin(k\lambda ) = \\ = {{(cos\lambda + {\mathbf{i}}sin\lambda )}^{k}}. \\ \end{gathered} $
Обозначим

(2.14)
$\begin{gathered} \Pi _{{ts}}^{ + }(r,\gamma ,\lambda )\, = \,r\Pi _{{ts}}^{{\text{*}}}(r,\gamma ,\lambda ) = f{{m}_{{\text{E}}}}\sum\limits_{n = 2}^\infty \sum\limits_{k = 1}^n \mathop {\left( {\frac{R}{r}} \right)}\nolimits^n {{P}_{{nk}}}(\gamma ) \times \\ \times \,\,({{C}_{{nk}}}cos(k\lambda ) + {{S}_{{nk}}}sin(k\lambda )),\,\,\,\,r = \left| {\mathbf{r}} \right|. \\ \end{gathered} $

Используя (2.5), (2.13), (2.14) находим

(2.15)
$\begin{gathered} \frac{\partial }{{\partial {\mathbf{u}}}}\left( {\Pi _{{ts}}^{ + }(r,\gamma ,\lambda )} \right) = \frac{{\partial \Pi _{{ts}}^{ + }}}{{\partial r}}\frac{{\partial r}}{{\partial {\mathbf{u}}}} + \frac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \gamma }}\frac{{\partial \gamma }}{{\partial {\mathbf{u}}}} + \frac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \lambda }}\frac{{\partial \lambda }}{{\partial {\mathbf{u}}}} = \\ = 2[( - {{l}_{1}}{{u}_{0}}\, + \,{{l}_{2}}{{u}_{2}})\, + \,( - {{l}_{1}}{{u}_{1}}\, + \,{{l}_{2}}{{u}_{3}}){\mathbf{i}}\, + \,( - {{l}_{1}}{{u}_{2}}\, + \,{{l}_{2}}{{u}_{0}}){\mathbf{j}} + \\ + \,\,( - {{l}_{1}}{{u}_{3}}\, + \,{{l}_{2}}{{u}_{1}}){\mathbf{k}}]\, + \,\frac{1}{{(x_{1}^{2}\, + \,x_{2}^{2})}}\frac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \lambda }}\left( {{{x}_{1}}\frac{{\partial {{x}_{2}}}}{{\partial {\mathbf{u}}}}\, - \,{{x}_{2}}\frac{{\partial {{x}_{1}}}}{{\partial {\mathbf{u}}}}} \right), \\ {{l}_{1}} = \tfrac{\gamma }{r}\tfrac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \gamma }} - \tfrac{{\partial \Pi _{{ts}}^{ + }}}{{\partial r}},\,\,\,\,{{l}_{2}} = \tfrac{1}{r}\tfrac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \gamma }}. \\ \end{gathered} $

Подставляя в уравнения (2.1), (2.2) соотношение $\Pi * = \Pi _{z}^{{\text{*}}} + \Pi _{{ts}}^{{\text{*}}}$ и соотношения (2.6)–(2.8), (2.14), (2.15) и переходя к скалярной записи, получаем следующие уравнения возмущенного движения спутника в гравитационном поле Земли с учетом его тессеральных и секториальных гармоник в переменных Кустаанхеймо–Штифеля:

(2.16)
$\begin{gathered} \frac{{{{d}^{2}}{{u}_{j}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{j}}\, = \,\frac{1}{2}\left( {\left( {\frac{\gamma }{r}\frac{{\partial {{\Pi }^{ + }}}}{{\partial \gamma }}\, - \,\frac{{\partial {{\Pi }^{ + }}}}{{\partial r}}} \right){{u}_{j}}\, - \,\left( {\frac{1}{r}\frac{{\partial {{\Pi }^{ + }}}}{{\partial \gamma }}} \right)u_{j}^{{\text{*}}}} \right) - \\ - \,\,\frac{1}{{4\left( {x_{1}^{2} + x_{2}^{2}} \right)}}\frac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \lambda }}\left( {{{x}_{1}}\frac{{\partial {{x}_{2}}}}{{\partial {{u}_{j}}}} - {{x}_{2}}\frac{{\partial {{x}_{1}}}}{{\partial {{u}_{j}}}}} \right) + \frac{1}{2}r{{q}_{j}}, \\ j = 0,1,2,3; \\ \frac{{dh{\text{*}}}}{{d\tau }} = r\frac{{\partial \Pi _{{ts}}^{{\text{*}}}}}{{\partial t}} + 2\sum\limits_{j = 0}^3 \left( {\frac{{d{{u}_{j}}}}{{d\tau }}{{q}_{j}}} \right),\,\,\,\,\frac{{dt}}{{d\tau }} = r. \\ \end{gathered} $
Здесь

$\begin{gathered} u_{0}^{{\text{*}}} = {{u}_{2}},\,\,\,\,u_{1}^{{\text{*}}} = {{u}_{3}},\,\,\,\,u_{2}^{{\text{*}}} = {{u}_{0}},\,\,\,\,u_{3}^{{\text{*}}} = {{u}_{1}}, \\ {{\Pi }^{ + }} = \Pi _{z}^{ + } + \Pi _{{ts}}^{ + };\,\,\,\,{{x}_{1}} = u_{0}^{2} + u_{1}^{2} - u_{2}^{2} - u_{3}^{2}, \\ {{x}_{2}} = 2({{u}_{1}}{{u}_{2}} - {{u}_{0}}{{u}_{3}}), \\ \end{gathered} $

$\Pi _{z}^{ + }$ и $\Pi _{{ts}}^{ * }$ имеют вид (2.6) и (2.14).

3. УРАВНЕНИЯ ВОЗМУЩЕННОГО ДВИЖЕНИЯ СПУТНИКА В ГРАВИТАЦИОННОМ ПОЛЕ ЗЕМЛИ В МОДИФИЦИРОВАННЫХ ЧЕТЫРЕХМЕРНЫХ ПЕРЕМЕННЫХ

Получим другие уравнения движения спутника, которые, обладая всеми достоинствами выше приведенных уравнений в переменных Кустаанхеймо–Штифеля, имеют более простую и симметричную структуру. Для этого необходимо вместо переменных Кустаанхеймо–Штифеля использовать другие четырехмерные переменные, введенные автором в [26] (см. также [7, 8]).

В случае Кустаанхеймо–Штифеля ось $M{{Y}_{1}}$ ранее введенной вращающейся системы координат $Y$ была направлена по радиусу-вектору ${\mathbf{r}}$ центра масс спутника. Координаты ${{x}_{i}}$ спутника в системе координат $X$ в соответствии с (1.2) и (1.4) связаны с переменными Кустаанхеймо–Штифеля ${{u}_{j}}$ соотношениями

(3.1)
$\begin{gathered} {{x}_{1}} = r\left( {\lambda _{0}^{2} + \lambda _{1}^{2} - \lambda _{2}^{2} - \lambda _{3}^{2}} \right) = u_{0}^{2} + u_{1}^{2} - u_{2}^{2} - u_{3}^{2}, \\ {{x}_{2}} = 2r({{\lambda }_{1}}{{\lambda }_{2}} + {{\lambda }_{0}}{{\lambda }_{3}}) = 2({{u}_{1}}{{u}_{2}} - {{u}_{0}}{{u}_{3}}), \\ {{x}_{3}} = 2r({{\lambda }_{1}}{{\lambda }_{3}} - {{\lambda }_{0}}{{\lambda }_{2}}) = 2({{u}_{1}}{{u}_{3}} + {{u}_{0}}{{u}_{2}}); \\ {{u}_{0}} = {{r}^{{1/2}}}{{\lambda }_{0}},\,\,\,\,{{u}_{i}} = - {{r}^{{1/2}}}{{\lambda }_{i}},\,\,\,\,i = 1,2,3, \\ r = u_{0}^{2} + u_{1}^{2} + u_{2}^{2} + u_{3}^{2}, \\ \end{gathered} $
где ${{\lambda }_{j}}$ – по-прежнему, параметры Эйлера (Родрига–Гамильтона), характеризующие ориентацию системы координат $Y$ относительно $X$.

Направим по радиусу-вектору ${\mathbf{r}}$ не ось $M{{Y}_{1}}$ системы координат $Y$, а ось $M{{Y}_{3}}$. В этом случае все кватернионные уравнения разделов 1 и 2 сохраняют свой вид, лишь вместо орта ${\mathbf{i}}$ необходимо взять орт ${\mathbf{k}}$ (это, кстати, демонстрирует удобство использования кватернионных моделей астродинамики). Новые переменные ${{u}_{j}},$ определяемые через параметры Родрига–Гамильтона, как и в случае Кустаанхеймо–Штифеля, формулами (3.1), будут связаны с координатами ${{x}_{i}}$ соотношениями

(3.2)
$\begin{gathered} {{x}_{1}} = 2r({{\lambda }_{1}}{{\lambda }_{3}} + {{\lambda }_{0}}{{\lambda }_{2}}) = 2({{u}_{1}}{{u}_{3}} - {{u}_{0}}{{u}_{2}}), \\ {{x}_{2}} = 2r({{\lambda }_{2}}{{\lambda }_{3}} - {{\lambda }_{0}}{{\lambda }_{1}}) = 2({{u}_{2}}{{u}_{3}} + {{u}_{0}}{{u}_{1}}), \\ {{x}_{3}} = r\left( {\lambda _{0}^{2} - \lambda _{1}^{2} - \lambda _{2}^{2} + \lambda _{3}^{2}} \right) = u_{0}^{2} - u_{1}^{2} - u_{2}^{2} + u_{3}^{2}, \\ \end{gathered} $
которые в кватернионной записи принимают вид
(3.3)
${{{\mathbf{r}}}_{x}} = {{x}_{1}}{\mathbf{i}} + {{x}_{2}}{\mathbf{j}} + {{x}_{3}}{\mathbf{k}} = r{\mathbf{\lambda }} \circ {\mathbf{k}} \circ {\mathbf{\tilde {\lambda }}} = {\mathbf{\tilde {u}}} \circ {\mathbf{k}} \circ {\mathbf{u}},$
где новая кватернионная переменная u = u0 + $ + \,{{u}_{1}}{\mathbf{i}} + {{u}_{2}}{\mathbf{j}} + {{u}_{3}}{\mathbf{k}}$ имеет смысл, отличный от кватернионной переменной, использованной в разделах 1 и 2.

Расстояние $r,$ по-прежнему, находится через новые переменные ${{u}_{j}}$ по формуле $r = u_{0}^{2} + u_{1}^{2} + u_{2}^{2} + u_{3}^{2}.$

Переменная $\gamma ,$ фигурирующая в потенциале гравитационного поля Земли, выражается через переменные Кустаанхеймо–Штифеля с помощью соотношений (2.5). Для новых переменных имеем:

(3.4)
$\gamma = \lambda _{0}^{2} - \lambda _{1}^{2} - \lambda _{2}^{2} + \lambda _{3}^{2} = {{r}^{{ - 1}}}(u_{0}^{2} - u_{1}^{2} - u_{2}^{2} + u_{3}^{2}).$

Отсюда, учитывая (3.1) и равенство $\lambda _{0}^{2} + \lambda _{1}^{2}$ + $ + \,\lambda _{2}^{2} + \lambda _{3}^{2} = 1,$ получаем

(3.5)
$\begin{gathered} \gamma = 1 - 2(\lambda _{1}^{2} + \lambda _{2}^{2}) = 2(\lambda _{0}^{2} + \lambda _{3}^{2}) - 1 = \\ = 1 - 2{{r}^{{ - 1}}}(u_{1}^{2} + u_{2}^{2}) = 2{{r}^{{ - 1}}}(u_{0}^{2} + u_{3}^{2}) - 1, \\ r = u_{0}^{2} + u_{1}^{2} + u_{2}^{2} + u_{3}^{2}. \\ \end{gathered} $

Из сопоставления (2.5) с (3.5) видно, что выражения переменной $\gamma $, от которой зависит потенциал $\Pi _{z}^{ * },$ через новые переменные ${{u}_{j}}$ могут быть представлены в двух различных формах и имеют более простую и симметричную структуру, что и позволяет получить более простые и симметричные, чем в случае использования переменных Кустаанхеймо–Штифеля, уравнения движения спутника.

Потенциал гравитационного поля Земли

$\begin{gathered} {{\Pi }_{{\text{E}}}} = \Pi + \Pi *,\,\,\,\,\Pi = \Pi (r), \\ \Pi * = \Pi _{z}^{ * }(r,\gamma ) + \Pi _{{ts}}^{ * }(r,\gamma ,\lambda ). \\ \end{gathered} $

Обозначим

(3.6)
$\begin{gathered} {{\Pi }^{ + }}(r,\gamma ,\lambda ) = r\Pi * = r\left( {\Pi _{z}^{ * }(r,\gamma ) + \Pi _{{ts}}^{ * }(r,\gamma ,\lambda )} \right) = \\ = \Pi _{z}^{ + }(r,\gamma ) + \Pi _{{ts}}^{ + }(r,\gamma ,\lambda ), \\ \end{gathered} $
где $\Pi _{z}^{ * }(r,\gamma )$ и $\Pi _{{ts}}^{ * }(r,\gamma ,\lambda )$ определены формулами (2.4) и (2.14).

Используя (3.4), (3.6), находим

(3.7)
$\begin{gathered} \frac{\partial }{{\partial {\mathbf{u}}}}\left( {{{\Pi }^{ + }}(r,\gamma ,\lambda )} \right) = \frac{{\partial {{\Pi }^{ + }}}}{{\partial r}}\frac{{\partial r}}{{\partial {\mathbf{u}}}} + \frac{{\partial {{\Pi }^{ + }}}}{{\partial \gamma }}\frac{{\partial \gamma }}{{\partial {\mathbf{u}}}} + \frac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \lambda }}\frac{{\partial \lambda }}{{\partial {\mathbf{u}}}} = \\ = 2\frac{{\partial {{\Pi }^{ + }}}}{{\partial r}}{\mathbf{u}} - 2{{r}^{{ - 1}}}\frac{{\partial {{\Pi }^{ + }}}}{{\partial \gamma }}[(\gamma - 1){{u}_{0}} + (\gamma + 1){{u}_{1}}{\mathbf{i}} + \\ + \,\,(\gamma + 1){{u}_{2}}{\mathbf{j}} + (\gamma - 1){{u}_{3}}{\mathbf{k}}] + \frac{1}{{(x_{1}^{2} + x_{2}^{2})}}\frac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \lambda }} \times \\ \times \,\,\left( {{{x}_{1}}\frac{{\partial {{x}_{2}}}}{{\partial {\mathbf{u}}}} - {{x}_{2}}\frac{{\partial {{x}_{1}}}}{{\partial {\mathbf{u}}}}} \right). \\ \end{gathered} $

Общая запись кватернионных уравнений движения ИСЗ в переменных Кустаанхеймо–Штифеля и в модифицированных переменных имеет один и тот же вид (2.1), (2.2). Поэтому подставляя в уравнения (2.1), (2.2) соотношения (3.6), (3.7) и переходя к скалярной записи, получаем следующие уравнения возмущенного движения спутника в гравитационном поле Земли в новых (модифицированных) четырехмерных переменных ${{u}_{j}}{\text{:}}$

(3.8)
$\begin{gathered} \frac{{{{d}^{2}}{{u}_{k}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{k}} = \frac{1}{2}\left( {\frac{{\gamma - 1}}{r}\frac{{\partial {{\Pi }^{ + }}}}{{\partial \gamma }} - \frac{{\partial {{\Pi }^{ + }}}}{{\partial r}}} \right){{u}_{k}} - \\ - \,\,\frac{1}{4}\frac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \lambda }}\frac{{\partial \lambda }}{{\partial {{u}_{k}}}} + \frac{1}{2}r{{q}_{k}},\,\,\,\,k = 0,3, \\ \end{gathered} $
(3.9)
$\begin{gathered} \frac{{{{d}^{2}}{{u}_{s}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{s}} = \frac{1}{2}\left( {\frac{{\gamma + 1}}{r}\frac{{\partial {{\Pi }^{ + }}}}{{\partial \gamma }} - \frac{{\partial {{\Pi }^{ + }}}}{{\partial r}}} \right){{u}_{s}} - \\ - \,\,\frac{1}{4}\frac{{\partial \Pi _{{ts}}^{ + }}}{{\partial \lambda }}\frac{{\partial \lambda }}{{\partial {{u}_{s}}}} + \frac{1}{2}r{{q}_{s}},\,\,\,\,s = 1,2, \\ \end{gathered} $
(3.10)
$\frac{{dh{\text{*}}}}{{d\tau }} = - r{{\Omega }_{{\text{E}}}}\frac{{\partial \Pi _{{ts}}^{ * }}}{{\partial \lambda }} + 2{\text{scal}}\left( {\frac{{d{\mathbf{\tilde {u}}}}}{{d\tau }} \circ {\mathbf{q}}} \right),\,\,\,\,\frac{{dt}}{{d\tau }} = r.$
Здесь

$\begin{gathered} \frac{{\partial \lambda }}{{\partial {{u}_{j}}}} = \frac{1}{{x_{1}^{2} + x_{2}^{2}}}\left( {{{x}_{1}}\frac{{\partial {{x}_{2}}}}{{\partial {{u}_{j}}}} - {{x}_{2}}\frac{{\partial {{x}_{1}}}}{{\partial {{u}_{j}}}}} \right), \\ \frac{{\partial \lambda }}{{\partial {{u}_{0}}}} = \frac{2}{{x_{1}^{2} + x_{2}^{2}}}({{x}_{1}}{{u}_{1}} + {{x}_{2}}{{u}_{2}}),\,\,\,\,\frac{{\partial \lambda }}{{\partial {{u}_{1}}}} = \frac{2}{{x_{1}^{2} + x_{2}^{2}}} \times \\ \times \,\,({{x}_{1}}{{u}_{0}} - {{x}_{2}}{{u}_{3}}),\,\,\,\,\frac{{\partial \lambda }}{{\partial {{u}_{2}}}} = \frac{2}{{x_{1}^{2} + x_{2}^{2}}}({{x}_{1}}{{u}_{3}} + {{x}_{2}}{{u}_{0}}), \\ \frac{{\partial \lambda }}{{\partial {{u}_{3}}}} = \frac{2}{{x_{1}^{2} + x_{2}^{2}}}({{x}_{1}}{{u}_{2}} - {{x}_{2}}{{u}_{1}}),\,\,\,\,{{\Pi }^{ + }} = \Pi _{z}^{ + }(r,\gamma ) + \\ + \,\,\Pi _{{ts}}^{ + }(r,\gamma ,\lambda ),\,\,\,\,\Pi _{z}^{ + }(r,\gamma ) = r\Pi _{z}^{ * }(r,\gamma ),\,\,\,\,\Pi _{{ts}}^{ + }(r,\gamma ,\lambda ) = \\ = r\Pi _{{ts}}^{ * }(r,\gamma ,\lambda ),\,\,\,\,{{x}_{1}} = 2({{u}_{1}}{{u}_{3}} - {{u}_{0}}{{u}_{2}}),{{x}_{2}} = 2({{u}_{2}}{{u}_{3}} + {{u}_{0}}{{u}_{1}}). \\ \end{gathered} $

Компоненты кватерниона

${\mathbf{q}} = {{q}_{0}} + {{q}_{1}}{\mathbf{i}} + {{q}_{2}}{\mathbf{j}} + {{q}_{3}}{\mathbf{k}} = - {\mathbf{k}} \circ {\mathbf{u}} \circ {{{\mathbf{p}}}_{x}}$
имеют следующий вид:

$\begin{gathered} {{q}_{0}} = {{u}_{0}}{{p}_{3}} + {{u}_{1}}{{p}_{2}} - {{u}_{2}}{{p}_{1}},\,\,\,\,{{q}_{1}} = {{u}_{0}}{{p}_{2}} - {{u}_{1}}{{p}_{3}} + {{u}_{3}}{{p}_{1}}, \hfill \\ {{q}_{2}} = - {{u}_{0}}{{p}_{1}} - {{u}_{2}}{{p}_{3}} + {{u}_{3}}{{p}_{2}},\,\,\,\,{{q}_{3}} = {{u}_{1}}{{p}_{1}} + {{u}_{2}}{{p}_{2}} + {{u}_{3}}{{p}_{3}}. \hfill \\ \end{gathered} $

В правые части уравнений (3.8)–(3.10) входят частные производные от потенциалов ${{\Pi }^{ + }}$ и $\Pi _{{ts}}^{ * }$ по переменным $r,\gamma ,\lambda .$ После их вычисления в полученные выражения должны быть подставлены величины $r,\gamma ,\lambda ,$ выраженные через переменные ${{u}_{j}}$ в соответствии с соотношениями (3.5) и соотношениями

$\lambda = {{\lambda }_{a}} - {{\Omega }_{{\text{E}}}}t,\,\,\,\,{{\lambda }_{a}} = {\text{arctg}}\frac{{{{x}_{2}}}}{{{{x}_{1}}}} = {\text{arctg}}\frac{{{{u}_{1}}{{u}_{3}} - {{u}_{0}}{{u}_{2}}}}{{{{u}_{0}}{{u}_{1}} + {{u}_{2}}{{u}_{3}}}}.$

Полученные уравнения движения спутника в новых (модифицированных) переменных ${{u}_{j}}$ (3.8)–(3.10) имеют более простую и симметричную структуру в сравнении с уравнениями в переменных Кустаанхеймо–Штифеля (2.16) что позволяет построить более эффективные вычислительные алгоритмы численного интегрирования уравнений движения спутника. Уравнения образуют замкнутую систему нелинейных дифференциальных уравнений десятого порядка относительно новых переменных ${{u}_{j}}(j = 0,1,2,3),$ их первых производных $d{{u}_{j}}{\text{/}}d\tau ,$ полной энергии $h{\text{*}}$ и времени $t.$

Для нахождения координат ${{x}_{i}}$ и проекций скорости $\mathop {\dot {x}}\nolimits_i $ спутника через новые переменные ${{u}_{j}}$ необходимо воспользоваться формулами (3.2) или (3.3) и кватернионным соотношением

${{{\mathbf{v}}}_{x}} = 2{\mathbf{\tilde {u}}} \circ {\mathbf{k}} \circ {\mathbf{\dot {u}}} = 2{{r}^{{ - 1}}}{\mathbf{\tilde {u}}} \circ {\mathbf{k}} \circ (d{\mathbf{u}}{\text{/}}d\tau ).$

Дифференциальное уравнение для расстояния $r$ от спутника до центра масс Земли имеет следующий вид:

$\begin{gathered} \frac{{{{d}^{2}}r}}{{d{{\tau }^{2}}}} - 2h{\text{*}}r = f{{m}_{{\text{E}}}} - \frac{\partial }{{\partial r}}\left( {r{{\Pi }^{ + }}} \right) + r{\text{scal}}({\mathbf{\tilde {u}}} \circ {\mathbf{q}}), \\ {\mathbf{q}} = - {\mathbf{k}} \circ {\mathbf{u}} \circ {{{\mathbf{p}}}_{x}}. \\ \end{gathered} $

4. УРАВНЕНИЯ ДВИЖЕНИЯ СПУТНИКА В ГРАВИТАЦИОННОМ ПОЛЕ ЗЕМЛИ С УЧЕТОМ ЕГО ЗОНАЛЬНЫХ ГАРМОНИК В МОДИФИЦИРОВАННЫХ ЧЕТЫРЕХМЕРНЫХ ПЕРЕМЕННЫХ И ИХ АНАЛИЗ

Уравнения движения спутника в гравитационном поле Земли с учетом его зональных гармоник в модифицированных переменных ${{u}_{j}}$ получаются из уравнений (3.8)–(3.10) в случае, когда возмущающее ускорение ${\mathbf{p}}$ спутника отсутствует, т.е. когда величины ${{q}_{j}} = 0,$ а также когда потенциал ${{\Pi }^{ + }} = \Pi _{z}^{ + }$ ($\Pi _{{ts}}^{ * } = 0$). Эти уравнения имеют следующий вид:

(4.1)
$\frac{{{{d}^{2}}{{u}_{k}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{k}} = \frac{1}{2}\left( {\frac{{\gamma - 1}}{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right){{u}_{k}},\,\,\,\,k = 0,3,$
(4.2)
$\frac{{{{d}^{2}}{{u}_{s}}}}{{d{{\tau }^{2}}}} - \frac{1}{2}h{\text{*}}{{u}_{s}} = \frac{1}{2}\left( {\frac{{\gamma + 1}}{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right){{u}_{s}},\,\,\,s = 1,2,$
(4.3)
$\frac{{dt}}{{d\tau }} = r,\,\,\,\,h* = {\text{const}}.$
Здесь $\Pi _{z}^{ + }$ определяется формулой (2.6), а $r$ и $\gamma $формулами (3.5).

Дифференциальное уравнение для расстояния $r$ от спутника до центра масс Земли в этом случае принимает следующий вид:

(4.4)
$\frac{{{{d}^{2}}r}}{{d{{\tau }^{2}}}} - 2h{\text{*}}r = f{{m}_{{\text{E}}}} - \frac{\partial }{{\partial r}}(r\Pi _{z}^{ + }).$

Отметим, что уравнения (4.1)(4.3) можно также получить из уравнений, приведенных в [26] (см. также [7, 8]), если фигурирующий в них массовый параметр $m = {{m}_{0}}{\text{/}}({{m}_{0}} + {{m}_{1}})$ (где ${{m}_{0}}$ – масса Земли, ${{m}_{1}}$ – масса спутника) положить равным единице.

Уравнения (4.1)(4.3) движения спутника, дополненные соотношениями (2.6) и соотношениями

(4.5)
$\begin{gathered} r = u_{0}^{2} + u_{1}^{2} + u_{2}^{2} + u_{3}^{2}, \\ \gamma = 1 - 2{{r}^{{ - 1}}}(u_{1}^{2} + u_{2}^{2}) = 2{{r}^{{ - 1}}}(u_{0}^{2} + u_{3}^{2}) - 1, \\ \end{gathered} $
образуют замкнутую систему дифференциальных уравнений девятого порядка относительно модифицированных переменных ${{u}_{j}}(j = 0,1,2,3)$ и времени $t.$ Эти уравнения обладают всеми достоинствами уравнений движения спутника в переменных Кустаанхеймо–Штифеля (2.12), образующими в рассматриваемом случае систему дифференциальных уравнений того же порядка. Однако правые части уравнений (4.1) и (4.2) имеют в сравнении с уравнениями в переменных Кустаанхеймо–Штифеля более простую и симметричную структуру.

Отметим, что уравнения (4.1) и (4.2) могут рассматриваться независимо от уравнения (4.3) для времени $t$ и образуют замкнутую систему дифференциальных уравнений восьмого порядка относительно модифицированных переменных ${{u}_{j}}(j = 0,1,2,3).$ Отметим также, что в отличие от уравнений в переменных Кустаанхеймо–Штифеля (2.12), уравнения (4.1), (4.2), как это видно из (3.5), распадаются на две независимые симметричные подсистемы, если расстояние $r$ определять из уравнения (4.4). Поэтому из приведенных уравнений выделяется замкнутая система дифференциальных уравнений шестого порядка. Она образуется либо уравнениями (4.1), (4.4) и соотношением (2.6), в которых после вычисления частных производных следует положить

$\gamma = 2{{r}^{{ - 1}}}(u_{0}^{2} + u_{3}^{2}) - 1,\,\,\,\,r = u_{0}^{2} + u_{1}^{2} + u_{2}^{2} + u_{3}^{2},$

либо уравнениями (4.2), (4.4) и соотношением (2.6), в которых после вычисления частных производных следует положить

$\gamma = 1 - 2{{r}^{{ - 1}}}(u_{1}^{2} + u_{2}^{2}),\,\,\,\,r = u_{0}^{2} + u_{1}^{2} + u_{2}^{2} + u_{3}^{2}.$

Если предположить, что решения ${{u}_{0}}(\tau ),$ ${{u}_{3}}(\tau ),$ $r(\tau )$ (или решения ${{u}_{1}}(\tau ),$ ${{u}_{2}}(\tau ),$ $r(\tau )$) системы уравнений (4.1), (4.4) (или (4.2), (4.4)) найдены, то уравнения (4.2) (или (4.1)) будут являться независимыми линейными дифференциальными уравнениями второго порядка с переменными коэффициентами.

4.1. Первые интегралы уравнений движения спутника. В дальнейшем будем обозначать первую производную по независимой переменной $\tau $ одним верхним штрихом, а вторую производную по этой переменной – двумя верхними штрихами.

Умножим каждое из уравнений системы (4.1) и (4.2) на соответствующую этому уравнению первую производную $u_{j}^{'}.$ Затем сложим полученные уравнения и проинтегрируем результат сложения. В итоге получим следующий первый интеграл дифференциальных уравнений движения спутника (4.1) и (4.2):

(4.6)
$H\, = \,\frac{1}{2}\sum\limits_{j = 0}^3 u_{j}^{{'2}}\, - \,\frac{1}{4}h{\text{*}}\sum\limits_{j = 0}^3 \mathop {{{u}_{j}}}\nolimits^2 \, + \,\frac{1}{4}\Pi _{z}^{ + }(r,\gamma )\, = \,{{c}_{H}}\, = \,{\text{const}}.$
Здесь потенциал $\Pi _{z}^{ + }$ определен формулой (2.6), величины $r$ и $\gamma $ определены соотношениями (4.5), постоянная интегрирования ${{c}_{H}}\, = \,(1{\text{/}}4)\mu \, = \,(1{\text{/}}4)f{{m}_{{\text{E}}}},$ $h{\kern 1pt} *\, = \,{\text{const}}.$

Умножим первое уравнение подсистемы (4.1) на ${{u}_{3}},$ а второе уравнение этой подсистемы – на ${{u}_{0}}.$ Вычтем из первого полученного уравнения второе полученное уравнение и проинтегрируем результат вычитания. Аналогичное проделаем с уравнениями подсистемы (4.2). В итоге получим следующие первые интегралы дифференциальных уравнений движения спутника (4.1) и (4.2):

(4.7)
${{u}_{3}}u_{0}^{'}\, - \,{{u}_{0}}u_{3}^{'}\, = \,{{c}_{1}}\, = \,{\text{const}},\,\,\,{{u}_{2}}u_{1}^{'}\, - \,{{u}_{1}}u_{2}^{'}\, = \,{{c}_{2}}\, = \,{\text{const}}.$

Билинейному соотношению (1.6) в переменных Кустаанхемо–Штифеля будет соответствовать другое билинейное соотношение в модифицированных переменных ${{u}_{j}},$ имеющее вид

${{u}_{3}}u_{0}^{'} - {{u}_{0}}u_{3}^{'} + {{u}_{2}}u_{1}^{'} - {{u}_{1}}u_{2}^{'} = 0.$

Из сопоставления первых интегралов (4.7) с этим билинейным соотношением следует, что постоянные интегрирования ${{c}_{1}}$ и ${{c}_{2}}$ должны быть связаны условием ${{c}_{1}} = - {{c}_{2}}.$

Таким образом, уравнения движения спутника (4.1), (4.2) в модифицированных переменных ${{u}_{j}}$ имеют три первых интеграла (4.6) и (4.7).

4.2. Уравнения движения спутника в гамильтоновой форме. Уравнения движения спутника (4.1), (4.2) в переменных ${{u}_{j}}$ можно записать в следующей гамильтоновой форме:

$\frac{{d{{u}_{j}}}}{{d\tau }} = \frac{{\partial H}}{{\partial {{w}_{j}}}},\,\,\,\,\frac{{d{{w}_{j}}}}{{d\tau }} = - \frac{{\partial H}}{{\partial {{u}_{j}}}},\,\,\,\,j = 0,1,2,3.$
Здесь ${{u}_{j}}$ – обобщенные координаты, ${{w}_{j}} = u_{j}^{'}$ – обобщенные импульсы, $H$ – функция Гамильтона, определяемая формулой (4.6).

4.3. Преобразование уравнений движения спутника и построение замкнутых систем уравнений меньшей размерности. Из системы дифференциальных уравнений движения спутника (4.1), (4.2) в модифицированных переменных ${{u}_{j}}$ восьмого порядка можно получить системы уравнений движения спутника в других переменных меньшей размерности.

Для этого выполним следующие преобразования уравнений (4.1), (4.2). Умножим первое уравнение подсистемы (4.1) на $2{{u}_{0}},$ а второе уравнение – на $2{{u}_{3}}$ и сложим полученные уравнения. Аналогично, умножим первое уравнение подсистемы (4.2) на $2{{u}_{1}},$ а второе уравнение – на $2{{u}_{2}}$ и также сложим полученные уравнения. Затем умножим первое уравнение подсистемы (4.1) на $u_{0}^{'},$ а второе уравнение – на $u_{3}^{'}$ и сложим полученные уравнения. Далее умножим первое уравнение подсистемы (4.2) на $u_{1}^{'},$ а второе уравнение – на $u_{2}^{'}$ и также сложим полученные уравнения. Введем новые переменные

(4.8)
$\begin{gathered} \kappa = u_{0}^{2} + u_{3}^{2},\,\,\,\,\nu = u_{1}^{2} + u_{2}^{2}, \\ \alpha = u_{0}^{{'{\text{2}}}} + u_{3}^{{'{\text{2}}}},\,\,\,\beta = u_{1}^{{'{\text{2}}}} + u_{2}^{{'{\text{2}}}}. \\ \end{gathered} $

В итоге указанных преобразований получим следующую систему уравнений движения спутника в переменных $\kappa ,\alpha ,\nu ,\beta $ шестого порядка:

(4.9)
$\kappa {\text{''}} - 2\alpha - h{\text{*}}\kappa = \left( {\frac{{\gamma - 1}}{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\kappa ,$
(4.10)
$\alpha {\text{'}} - \frac{1}{2}h{\text{*}}\kappa {\text{'}} = \frac{1}{2}\left( {\frac{{\gamma - 1}}{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\kappa {\text{'}},$
(4.11)
$\nu {\text{''}} - 2\beta - h{\text{*}}\nu = \left( {\frac{{\gamma + 1}}{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\nu ,$
(4.12)
$\beta {\text{'}} - \frac{1}{2}h{\text{*}}\nu ' = \frac{1}{2}\left( {\frac{{\gamma + 1}}{r}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\nu {\text{'}},$
где

(4.13)
$\begin{gathered} r = \kappa + \nu ,\,\,\,\,\gamma = (\kappa - \nu ){\text{/}}(\kappa + \nu ),\,\,\,\,\gamma - 1 = - 2\nu {\text{/}}r,\,\,\,\, \\ \gamma + 1 = 2\kappa {\text{/}}r,\,\,\,\,h* = {\text{const}}. \\ \end{gathered} $

С учетом соотношений (4.13) уравнения (4.9)(4.12) принимают вид

(4.14)
$\kappa {\text{''}} - 2\alpha - h{\text{*}}\kappa = - \left( {\frac{{2\nu }}{{{{{(\kappa + \nu )}}^{2}}}}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} + \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\kappa ,$
(4.15)
$\alpha {\text{'}} - \frac{1}{2}h{\text{*}}\kappa {\text{'}} = - \frac{1}{2}\left( {\frac{{2\nu }}{{{{{(\kappa + \nu )}}^{2}}}}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} + \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\kappa {\text{'}},$
(4.16)
$\nu {\text{''}} - 2\beta - h{\text{*}}\nu = \left( {\frac{{2\kappa }}{{{{{(\kappa + \nu )}}^{2}}}}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\nu ,$
(4.17)
$\beta {\text{'}} - \frac{1}{2}h{\text{*}}\nu {\text{'}} = \frac{1}{2}\left( {\frac{{2\kappa }}{{{{{(\kappa + \nu )}}^{2}}}}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\nu {\text{'}}.$

Умножим уравнение (4.9) или (4.14) на $\kappa {\text{',}}$ а уравнение (4.10) или (4.15) – на $2\kappa $ и вычтем из первого полученного уравнения второе полученное уравнение. Умножим уравнение (4.11) или (4.16) на $\nu {\text{'}},$ а уравнение (4.12) или (4.17) – на $2\nu $ и вычтем из первого полученного уравнения второе полученное уравнение. Проинтегрируем уравнения, полученные в итоге вышеуказанных преобразований. В итоге получим следующие первые интегралы дифференциальных уравнений (4.9)–(4.12) или (4.14)–(4.17):

(4.18)
$\begin{gathered} {{\kappa }^{{'2}}} - 4\alpha \kappa = {{c}_{\kappa }} = {\text{const}}, \\ {{\nu }^{{'2}}} - 4\beta \nu = {{c}_{\nu }} = {\text{const}}. \\ \end{gathered} $

Перейдем в первых интегралах (4.18) от переменных $\kappa ,\alpha ,\nu ,\beta $ к переменным ${{u}_{j}}$ и $u_{j}^{'}$ в соответствии с формулами (4.8). Получим соотношения

$ - {{\kappa }^{{'2}}} + 4\alpha \kappa + {{c}_{\kappa }} = 4{{({{u}_{3}}u_{0}^{'} - {{u}_{0}}u_{3}^{'})}^{2}} + {{c}_{\kappa }} = 0,$
$ - {{\nu }^{{'2}}} + 4\beta \nu + {{c}_{\nu }} = 4{{({{u}_{2}}u_{1}^{'} - {{u}_{1}}u_{2}^{'})}^{2}} + {{c}_{\nu }} = 0.$

Из сравнения этих соотношений с первыми интегралами (4.7) уравнений движения спутника в переменных ${{u}_{j}}$ и $u_{j}^{'}$ следует, что постоянные интегрирования ${{c}_{\kappa }},$ ${{c}_{\nu }}$ и ${{c}_{1}},$ ${{c}_{2}},$ фигурирующие в первых интегралах уравнений движения спутника в переменных $\kappa ,\alpha ,\nu ,\beta $ и в переменных ${{u}_{j}}$ и $u_{j}^{'},$ связаны соотношениями

${{c}_{\kappa }} = - 4c_{1}^{2},\,\,\,\,{{c}_{\nu }} = - 4c_{2}^{2}.$

Поскольку ${{c}_{1}} = - {{c}_{2}},$ то в силу предыдущих равенств ${{c}_{\kappa }} = {{c}_{\nu }}.$

Сложим левые и правые части уравнений (4.15) и (4.17) и проинтегрируем полученное уравнение. Получим следующий первый интеграл дифференциальных уравнений движения спутника (4.14)–(4.17):

(4.19)
$\begin{gathered} \alpha + \beta - \frac{1}{2}h{\text{*}}(\kappa + \nu ) + \frac{1}{2}(\kappa + \nu )\Pi _{z}^{ * }(r,\gamma ) = \\ = 2H = 2{{c}_{H}} = {\text{const}}, \\ \end{gathered} $
где $r = \kappa + \nu ,\gamma = (\kappa - \nu ){\text{/}}(\kappa + \nu ).$

Полученные первые интегралы (4.18) и (4.19) уравнений движения спутника (4.14)–(4.17) или (4.9)–(4.12) в переменных $\kappa ,\alpha ,\nu ,\beta $ соответствуют первым интегралам (4.7) и (4.6) уравнений движения спутника (4.1), (4.2) в переменных ${{u}_{j}}.$

С помощью первых интегралов (4.18) уравнения (4.15) и (4.17) для переменных $\alpha $ и $\beta $ могут быть исключены из рассмотрения и вместо системы дифференциальных уравнений (4.14)–(4.17) шестого порядка будем иметь для изучения движения спутника следующую систему дифференциальных уравнений четвертого порядка:

(4.20)
$\begin{gathered} \kappa {\text{''}} - \frac{1}{{2\kappa }}\left( {{{\kappa }^{{'2}}} - {{c}_{\kappa }}} \right) - h{\text{*}}\kappa = \\ = - \left( {\frac{{2\nu }}{{{{{(\kappa + \nu )}}^{2}}}}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} + \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\kappa , \\ \end{gathered} $
(4.21)
$\nu {\text{''}} - \frac{1}{{2\nu }}({{\nu }^{{'2}}} - {{c}_{\nu }})\, - \,h{\text{*}}\nu \, = \,\left( {\frac{{2\kappa }}{{{{{(\kappa \, + \,\nu )}}^{2}}}}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }}\, - \,\frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\nu .$

Кроме системы дифференциальных уравнений (4.20) и (4.21) четвертого порядка для изучения движения спутника можно использовать следующую систему дифференциальных уравнений того же порядка в переменных $\kappa ,\alpha ,\nu ,\beta ,$ получающуюся в результате объединений первых интегралов (4.18) и уравнений (4.15), (4.17):

(4.22)
${{\kappa }^{{'2}}}\, - \,4\alpha \kappa \, = \,{{c}_{\kappa }}\, = \,{\text{const}},\,\,\,{{\nu }^{{'2}}}\, - \,4\beta \nu \, = \,{{c}_{\nu }}\, = \,{\text{const}},$
(4.23)
$\alpha {\text{'}} - \frac{1}{2}h{\text{*}}\kappa {\text{'}} = - \frac{1}{2}\left( {\frac{{2\nu }}{{{{{(\kappa + \nu )}}^{2}}}}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} + \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\kappa {\text{'}},$
(4.24)
$\beta {\text{'}} - \frac{1}{2}h{\text{*}}\nu {\text{'}} = \frac{1}{2}\left( {\frac{{2\kappa }}{{{{{(\kappa + \nu )}}^{2}}}}\frac{{\partial \Pi _{z}^{ + }}}{{\partial \gamma }} - \frac{{\partial \Pi _{z}^{ + }}}{{\partial r}}} \right)\nu {\text{'}}.$
Здесь в правые части двух последних уравнений после нахождения частных производных необходимо использовать равенства $r = \kappa + \nu ,$ γ = $ = (\kappa - \nu ){\text{/}}(\kappa + \nu ).$

Переменные $\alpha $ и $\beta $ связаны соотношением (4.19) (входят в первый интеграл уравнений движения спутника в виде их суммы). Поэтому уравнение для переменной $\alpha $ или $\beta $ может быть исключено из состава системы уравнений (4.22)–(4.24). В итоге будем иметь замкнутую систему трех дифференциальных уравнений первого порядка относительно переменных $\kappa ,\nu ,\alpha $ или $\kappa ,\nu ,\beta .$

В уравнения движения спутника в гравитационном поле Земли с учетом его зональных гармоник входит потенциал $\Pi _{z}^{ * }(r,\gamma ),$ являющийся функцией расстояния $r$ спутника до центра масс Земли и геоцентрической широты $\varphi $ (точнее, переменной $\gamma = sin\varphi $). Поэтому представляет интерес получение замкнутой системы дифференциальных уравнений относительно лишь переменных $r$ и $\gamma .$

Для этого умножим первое уравнение подсистемы (4.22) на $\nu ,$ а второе уравнение этой подсистемы – на $\kappa $ и сложим левые и правые части полученных уравнений. Преобразуем новое полученное уравнение с учетом равенства ${{c}_{\kappa }} = {{c}_{\nu }}.$ В итоге получим уравнение

(4.25)
${{r}^{2}}{{\gamma }^{{'2}}} + \left[ {{{r}^{{'2}}} - 4r(\alpha + \beta )} \right]\left( {1 - {{\gamma }^{2}}} \right) = 4{{c}_{\kappa }}.$
Здесь $\alpha + \beta = 2H + (1/2)r(h{\text{*}} - \Pi _{z}^{ * }(r,\gamma ));$ $H = fm{\text{/}}4,$ $h*,$ ${{c}_{\kappa }}$ – постоянные величины.

Уравнение (4.25) и уравнение

$\frac{{{{d}^{2}}r}}{{d{{\tau }^{2}}}} - 2h{\text{*}}r = f{{m}_{{\text{E}}}} - \frac{\partial }{{\partial r}}({{r}^{2}}\Pi _{z}^{ + }(r,\gamma ))$
образуют искомую замкнутую систему дифференциальных уравнений третьго порядка относительно расстояния $r$ и переменной $\gamma = sin\varphi ,$ которая может быть использована для изучения движения спутника в гравитационном поле Земли с учетом его зональных гармоник.

ЗАКЛЮЧЕНИЕ

В работе получены кватернионные уравнения возмущенного движения искусственного спутника в гравитационном поле Земли с учетом его зональных, тессеральных и секториальных гармоник в четырехмерных переменных Кустаанхеймо–Штифеля и в модифицированных четырехмерных переменных. В этих уравнениях основными переменными являются переменные Кустаанхеймо–Штифеля или модифицированные четырехмерные переменные, предложенные автором статьи, а также энергия движения спутника и время. Новая независимая переменная $\tau $ связана со временем $t$ дифференциальным соотношением $dt = rd\tau ,$ содержащим расстояние $r$ спутника до центра масс Земли. Эта замена времени эквивалентна аналитическому регулированию длины шага интегрирования дифференциальных уравнений движения спутника.

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

– они регулярны (не содержат особых точек типа сингулярности (деления на ноль)) для возмущенного движения спутника в ньютоновском (центральном) поле сил тяготения (при условии, что в описании действующих возмущающих сил не содержатся отрицательные степени расстояния спутника до центра Земли выше первой);

– линейны для невозмущенных кеплеровских движений, для эллиптического кеплеровского движения, когда кеплеровская энергия $h = {\text{const}} < 0,$ эти уравнения эквивалентны уравнениям движения четырехмерного одночастотного гармонического осциллятора, квадрат частоты которого равен $ - h{\text{/}}2;$

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

– близки к линейным уравнениям для возмущенных кеплеровских движений;

– позволяют во многих случаях представить правые части дифференциальных уравнений движения спутника в полиномиальной форме, удобной для их решения с помощью ЭВМ.

Уравнения удобны для применения методов нелинейной механики и высокоточных численных расчетов. Вместе с тем отметим, что при учете несферичности Земли в предлагаемых уравнениях движения спутника помимо слагаемых, описывающих ньютоновскую силу тяготения, появляются слагаемые, описывающие зональные, тессеральные и секториальные гармоники потенциала гравитационного поля Земли, нарушающие регулярность этих уравнений.

В отношении модифицированных четырехмерных переменных, использованных в статье наряду с переменными Кустаанхеймо–Штифеля, отметим следующее. Автором статьи были ранее даны [2, 3, 7] наглядные геометрическая и кинематическая интерпретаци регуляризующему преобразованию Кустаанхеймо–Штифеля (1.4) декартовых координат и билинейному соотношению (1.6), играющему, по словам Е. Штифеля и Г. Шейфеле, основную роль в их построении регулярной небесной механики [1, с. 29]. В соответствии с ними преобразование координат Кустаанхеймо–Штифеля означает переход от трехмерных декартовых координат спутника в основной системе координат к новым четырехмерным переменным, являющимся нормированными (с помощью расстояния $r$ спутника до притягивающего центра) параметрами Эйлера (Родрига–Гамильтона), характеризующими ориентацию введенной автором статьи вращающейся системы координат $Y,$ т.е. фактически означает запись исходных уравнений движения спутника во вращающейся системе координат с использованием для описания вращательного движения этой системы координат параметров Эйлера, являющихся компонентами кватерниона поворота Гамильтона.

Эти интерпретации позволили автору статьи ввести в рассмотрение модифицированные четырехмерные переменные [26] (см. также [7, 8]). Модифицированные переменные $u_{j}^{ * }$ связаны с переменными Кустаанхеймо–Штифеля ${{u}_{j}}$ соотношениями

$\begin{gathered} u_{0}^{{\text{*}}} = (1{\text{/}}2)({{u}_{0}} + {{u}_{1}} + {{u}_{2}} + {{u}_{3}}), \\ u_{1}^{{\text{*}}} = - (1{\text{/}}2)({{u}_{0}} - {{u}_{1}} - {{u}_{2}} + {{u}_{3}}), \\ u_{2}^{{\text{*}}} = - (1{\text{/}}2)({{u}_{0}} + {{u}_{1}} - {{u}_{2}} - {{u}_{3}}), \\ u_{0}^{{\text{*}}} = - (1{\text{/}}2)({{u}_{0}} - {{u}_{1}} + {{u}_{2}} - {{u}_{3}}) \\ \end{gathered} $
и являются их линейными композициями (ранее в статье верхняя звездочка у модифицированных переменных опускалась).

В кватернионной записи эти соотношения принимают вид ортогонального преобразования

${\mathbf{u}}* = \pi \circ {\mathbf{u}},\,\,\,\,\pi = \frac{1}{2}\left( {1 - {\mathbf{i}} - {\mathbf{j}} - {\mathbf{k}}} \right).$

С геометрической точки зрения эта замена переменных означает, что по радиусу-вектору ${\mathbf{r}}$ центра масс спутника направляется не ось $M{{Y}_{1}}$ вращающейся системы координат $Y,$ введенной автором статьи в случае использования переменных Кустаанхеймо–Штифеля, а ось $M{{Y}_{3}}$. В этом случае все кватернионные уравнения разделов 1 и 2 статьи в переменных Кустаанхеймо–Штифеля сохраняют свой вид, лишь вместо орта ${\mathbf{i}}$ берется орт ${\mathbf{k}}$ (это, как уже отмечалось, демонстрирует удобство использования кватернионных моделей астродинамики). При этом компонентами кватернионных переменных, присутствующих в новых уравнениях, являются уже не переменные Кустаанхеймо–Штифеля, а модифицированные четырехмерные переменные.

Отметим, что эффективность применения кватернионов для решения проблем регуляризации моделей небесной механики и астродинамики впервые была продемонстрирована автором статьи в работах [2, 3], а затем в работах [48, 27] (на приоритет автора статьи в области кватернионной регуляризации уравнений небесной механики указывается в работе [16]).

Полученные в статье кватернионные уравнения движения спутника в модифицированных четырехмерных переменных обладают всеми достоинствами уравнений движения спутника в переменных Кустаанхеймо–Штифеля (также полученных в статье), но имеют более простую и симметричную структуру. Это обусловлено тем, что выражения переменной $\gamma = sin\varphi ,$ от которой зависит потенциал гравитационного поля Земли, через новые переменные $u_{j}^{{\text{*}}}$ могут быть представлены в двух различных, более компактных симметричных формах

$\gamma = 1 - 2{{r}^{{ - 1}}}\left( {u_{1}^{{ * 2}} + u_{2}^{{ * 2}}} \right) = 2{{r}^{{ - 1}}}\left( {u_{0}^{{ * 2}} + u_{3}^{{ * 2}}} \right) - 1,$
в сравнении с ее представлением
$\gamma = 2{{r}^{{ - 1}}}({{u}_{1}}{{u}_{3}} + {{u}_{0}}{{u}_{2}})$
в переменных Кустаанхеймо–Штифеля (здесь, напомним, $\varphi $ – геоцентрическая широта), что и позволяет получить более простые и симметричные, чем в случае использования переменных Кустаанхеймо–Штифеля, уравнения движения спутника.

Более простые и симметричные структуры уравнений приводят к более эффективным вычислительным алгоритмам при численном интегрировании дифференциальных уравнений движения спутника.

Удобство и эффективность использования полученных уравнений движения спутника в модифицированных четырехмерных переменных для аналитического исследования движения спутника показано в статье на примере рассмотрения движения спутника в гравитационном поле Земли, в описании которого учитываются только зональные гармоники (помимо центральной (ньютоновской) составляющей потенциала поля тяготения). В последнем (4-ом) разделе работы, найдены первые интегралы этих уравнений движения спутника в рассматриваемом случае, предложены замены переменных и преобразования этих уравнений, позволившие получить для изучения движения спутника замкнутые системы дифференциальных уравнений меньшей размерности, в частности, системы уравнений четвертого и третьего порядков.

Полученные уравнения движения спутника в четырехмерных переменных Кустаанхеймо–Штифеля и в модифицированных четырехмерных переменных в отличие от классических уравнений, регулярны для движения спутника в ньютоновском гравитационном поле Земли под действием возмущающих сил, в описании которых не содержатся отрицательные степени расстояния спутника до центра Земли выше первой. При учете несферичности Земли в предлагаемых уравнениях помимо слагаемых, описывающих ньютоновскую силу тяготения, появляются, как уже отмечалось, слагаемые, описывающие зональные, тессеральные и секториальные гармоники потенциала гравитационного поля Земли, нарушающие регулярность этих уравнений. Проблема устранения особенностей, порождаемых этими гармониками потенциала гравитационного поля Земли, остается открытой и представляет несомненный интерес.

Отметим в этой связи, что в работах [57] автора статьи идеи кватернионной регуляризации уравнений возмущенной пространственной задачи двух тел были использованы для разработки теории кватернионной регуляризации векторного дифференциального уравнения возмущенного центрального движения материальной точки (1.1), в которых потенциал $\Pi = \Pi (r)$ – произвольная дифференцируемая функцияй расстояния $r$ от точки до центра силового поля. Были получены общие кватернионные дифференциальные уравнения возмущенного центрального движения материальной точки с регуляризующими функциями, установлены необходимые и достаточные условия их приводимости к удобному для аналитического и численного исследования осцилляторному виду (к виду уравнений движения четырехмерного возмущенного осциллятора, совершающего в случае невозмущенного центрального движения гармонические колебания с одинаковой частотой); получены различные (в том числе новые регулярные) системы кватернионных дифференциальных уравнений возмущенного центрального движения материальной точки в нормальной и осцилляторной формах, указаны их свойства и области использования.

Так, были получены кватернионные регулярные дифференциальные уравнения возмущенного центрального движения материальной точки осцилляторного вида, регулярные для потенциала

$\Pi (r) = - {{a}_{1}}{{r}^{{ - 1}}} - {{a}_{2}}{{r}^{{ - 2}}} - {{a}_{3}}{{r}^{{ - 3}}} - {{a}_{4}}{{r}^{{ - 4}}},\,\,\,\,{{a}_{i}} = {\text{const}}$
являющегося полиномом четвертой отрицательной степени расстояния $r$ до центра притяжения.

Построено с использованием кватернионов решение пространственной задачи невозмущенного центрального движения материальной точки для произвольного потенциала $\Pi = \Pi (r)$ в униформизированной форме, позволяющей избавиться от необходимости рассмотрения ветвления решений, возникающего при обходе критических точек типа полюса. Для одного из вариантов решения этой задачи было использовано униформизированное решение плоской задачи невозмущенного центрального движения материальной точки, полученное И.М. Беленьким в работе [28].

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

– запись уравнений движения во вращающейся системе координат с использованием в качестве параметров ориентации этой системы координат четырехмерных параметров Эйлера (Родрига–Гамильтона),

– использование вместо трехмерных декартовых координат новых четырехмерных координат типа переменных Кустаанхеймо–Штифеля,

– использование кватернионов для компактности и наглядности записи уравнений и проводимых преобразований,

– введение в качестве дополнительных переменных энергии движения, момента количества движения,

– введение регуляризующего преобразования времени (введение новой независимой переменной) с использованием расстояния и других характеристик движения спутника.

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

  1. Stiefel E.L., Scheifele G. Linear and Regular Celestial Mechanics. Berlin: Springer, 1971. (Штифель Е., Шейфеле Г. Линейная и регулярная небесная механика. M.: Наука, 1975.).

  2. Челноков Ю.Н. К регуляризации уравнений пространственной задачи двух тел // Изв. АН СССР. Механика твердого тела. 1981. № 6. С. 12–21.

  3. Челноков Ю.Н. О регулярных уравнениях пространственной задачи двух тел // Изв. АН СССР. Механика твердого тела. 1984. № 1. С. 151–158.

  4. Челноков Ю.Н. Применение кватернионов в теории орбитального движения искусственного спутника. Ч. 1 // Космич. исслед. 1992. Т. 30. Вып. 6. С. 759–770 (Cosmic Research. P. 612).

  5. Челноков Ю.Н. Кватернионная регуляризация и стабилизация возмущенного центрального движения. Ч. 1 // Изв. РАН. Механика твердого тела. 1993. № 1. С. 20–30.

  6. Челноков Ю.Н. Кватернионная регуляризация и стабилизация возмущенного центрального движения. Ч. 2 // Изв. РАН. Механика твердого тела. 1993. № 2. С. 3–11.

  7. Челноков Ю.Н. Кватернионные модели и методы динамики, навигации и управления движением. М.: ФИЗМАТЛИТ, 2011. 560 с.

  8. Челноков Ю.Н. Кватернионная регуляризация в небесной механике и астродинамике и управление траекторным движением. I // Космич. исслед. 2013. Т. 51. № 5. С. 389–401 (Cosmic Research. P. 350).

  9. Брумберг В.А. Аналитические алгоритмы небесной механики. М.: Наука, 1980. 208 с.

  10. Vivarelli M.D. The KS transformation in hypercomplex form // Celest. Mech. Dyn. Astron. 1983. № 29. P. 45–50.

  11. Бордовицына Т.В. Современные численные методы в задачах небесной механики. М.: Наука, 1984.

  12. Vrbik J. Celestial mechanics via quaternions // Can. J. Phys. 1994. № 72. P. 141–146.

  13. Vrbik J. Perturbed Kepler problem in quaternionic form // J. Phys. 1995. A 28. P. 193–198.

  14. Waldvogel J. Quaternions and the perturbed Kepler problem // Celest. Mech. Dyn. Astron. 2006. № 95. P. 201–212.

  15. Бордовицына Т.В., Авдюшев В.А. Теория движения искусственных спутников Земли. Аналитические и численные методы. Томск: Изд-во Том. ун-та, 2007.

  16. Waldvogel J. Quaternions for regularizing Celestial Mechanics: the right way // Celest. Mech. Dyn. Astron. 2008. № 102. P. 149–162.

  17. Челноков Ю.Н., Сапунков Я.Г. Построение оптимальных управлений и траекторий космического аппарата на основе регулярных кватернионных уравнений задачи двух тел // Космич. исслед. 1996. Т. 34. № 2. С. 150–158 (Cosmic Research. P. 137).

  18. Сапунков Я.Г. Применение KS-переменных к задаче оптимального управления космическим аппаратом // Космич. исслед. 1996. Т. 34. № 4. С. 428–433 (Cosmic Research. P. 395).

  19. Челноков Ю.Н., Юрко В.А. Кватернионное построение оптимальных управлений и траекторий движения космического аппарата в ньютоновском гравитационном поле // Изв. РАН. Механика твердого тела. 1996. № 6. С. 3–13.

  20. Сапунков Я.Г. Решение задач оптимального управления космическим аппаратом с ограниченной и импульсной тягой в KS-переменных // Мехатроника, автоматизация, управление. 2010. № 3. С. 73–78.

  21. Сапунков Я.Г., Челноков Ю.Н. Построение оптимальных управлений и траекторий движения центра масс космического аппарата, снабженного солнечным парусом и двигателем малой тяги, с использованием кватернионов и переменных Кустаанхеймо–Штифеля // Космич. исслед. 2014. Т. 52. № 6. С. 489–499 (Cosmic Research. P. 450).

  22. Челноков Ю.Н. Кватернионная регуляризация в небесной механике и астродинамике и управление траекторным движением. III // Космич. исслед. 2015. Т. 53. № 5. С. 430–446 (Cosmic Research. P. 394).

  23. Абалакин В.К., Аксенов Е.П., Гребеников Е.А., Демин В.Г., Рябов Ю.А. Справочное руководство по небесной механике и астродинамике. М.: Наука, 1976.

  24. Дубошин Г.Н. Небесная механика: Методы теории движения искусственных небесных тел. М.: Наука, 1983.

  25. Демин В.Г. Движение искусственного спутника в нецентральном поле тяготения. М.–Ижевск: НИЦ “Регулярная и хаотическая динамика”, Ижевский институт компьютерных исследований. 2010.

  26. Челноков Ю.Н. Применение кватернионов в теории орбитального движения искусственного спутника. Ч. 2 // Космич. исслед. 1993. Т. 31. Вып. 3. С. 3–15 (Cosmic Research. P. 409).

  27. Челноков Ю.Н. Кватернионная регуляризация уравнений возмущенной пространственной ограниченной задачи трех тел. I // Изв. РАН. Механика твердого тела. 2017. № 6. С. 24–54.

  28. Беленький И.М. Об одном методе униформизации решений в задачах центрального движения // Прикладная математика и механика. 1981. Т. 45. Вып. 1. С. 34–41.

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