Рассмотрим плоское движение материальной точки в гравитационном поле Земли. Если ее движение происходит вне атмосферы, то основной силой, дей-
*) Н. Н. Моисеев, Methods non-linear mechanics in the problems of the dynamics of setellites, Trans of XIV-th Astronautical congress, 1962.

ствующей на точку, будет сила притяжения Земли. Поэтому движение рассматриваемой точки будет описываться уравнением
\[
\frac{d^{2} r}{d t^{2}}=-\frac{\mu}{r^{3}} r+\varepsilon f
\]

где $\boldsymbol{r}$-радиус-вектор точки, $\mu$-гравитационная постоянная, $\boldsymbol{\varepsilon} \boldsymbol{f}$ — вектор возмущающих сил, $\varepsilon$-малый параметр (обозначения даны на рис. 32). Первое слагаемое, стоящее в правой части уравнения (8.1), — это напряженность гравитационного поля одного притягивающего центра. В действительности, силы притяжения Земли несколько отличаются от этой величины, поскольку Земля не представляет собой точечную массу. Разность вектора напряженности реального поля Земли и величины $\left(\mu / r^{3}\right) \boldsymbol{r}$ является одной из составляющих вектора $\varepsilon f$. Кроме того, к числу составляющих этого вектора мы отнесем притяжение Луны и планет, сопротивление атмосферы и т. д. Может оказаться также, что аппарат (который мы схематизируем как материальную точку) находится под дей.

Рис. 32. ствием тяги двигателя «малой тяги». В этом случае ускорение, которое создается этим двигателем, мы также отнесем к числу возмущающих сил.

Перепишем уравнение (8.1) в проекциях на оси полярной системы координат
\[
\left.\begin{array}{c}
\frac{d v_{r}}{d t}=\frac{v_{\varphi}^{2}}{r}-\frac{\mu}{r^{2}}+\varepsilon f_{r}, \\
\frac{d v_{\varphi}}{d t}=-\frac{v_{r} v_{\varphi}}{r}+\varepsilon f_{\varphi^{\cdot}}
\end{array}\right\}
\]

Здесь $v_{r}$ и $v_{\varphi}$ — проекции вектора скорости
\[
\boldsymbol{V}=\frac{d \boldsymbol{r}}{d t}=v_{r} \boldsymbol{r}^{0}+v_{\varphi} \boldsymbol{p}^{0}
\]

на направление радиуса-вектора и трансверсаль
\[
v_{r}=\frac{d r}{d t} ; \quad v_{\varphi}=r \frac{d \varphi}{d t} .
\]

Величина $\boldsymbol{h}$, определенная равенством
\[
\boldsymbol{h}=\boldsymbol{V} \times \boldsymbol{r} \text {, }
\]

носит название секториальной скорости. Используя равенства (8.3) и (8.4), легко найдем, что
\[
\boldsymbol{h}=\boldsymbol{r}^{2} \frac{d \varphi}{d t}\left(\boldsymbol{p}^{0} \times \boldsymbol{r}^{0}\right)
\]

Вводя обозначение
\[
u=\frac{1}{r},
\]

можно выражение угловсй скорости $d \varphi / d t$ записать так:
\[
\frac{d \varphi}{d t}=u^{2} h .
\]

Рассмотрим первое из уравнений (8.2). Заметим прежде всего, что, используя (8.5), его можно представить в виде
\[
\frac{d}{d t}\left(\frac{d r}{d t}\right)=u^{3} h^{2}-\mu u^{2}+\varepsilon f_{r} .
\]

Составим теперь уравнение, которому удовлетворяет функция $u(\varphi)$. Для этого проделаем следующие вычисления:
\[
\begin{array}{c}
\frac{d r}{d t}=\frac{d r}{d \varphi} u^{2} h=u^{2} h \frac{d}{d \varphi}\left(\frac{1}{u}\right)=-h \frac{d u}{d \varphi}, \\
\frac{d}{d t}\left(\frac{d r}{d t}\right)=-\frac{d}{d \varphi}\left(h \frac{d u}{d \varphi}\right) u^{2} h=-u^{2} h\left\{\frac{d h}{d \varphi} \frac{d u}{d \varphi}+h \frac{d^{2} u}{d \varphi^{2}}\right\} .
\end{array}
\]

Используя эти вспомогательные выражения и уравнение (8.6), получим
\[
\frac{d^{2} u}{d \varphi^{2}}+\frac{1}{h} \frac{d u}{d \varphi} \frac{d h}{d \varphi}+u=\frac{\mu}{h^{2}}-\frac{\varepsilon f_{r}}{u^{2} h^{2}} .
\]

Составим далее уравнение, которому удовлетворяет величина $h$; так как
\[
v_{\varphi}=r \frac{d \varphi}{d t}=u h, \quad v_{r}=\frac{d r}{d t}=-h \frac{d u}{d \varphi},
\]

то второе из уравнений системы (8.2) мы можем переписать в виде
\[
\frac{d}{d t}(u h)=u^{2} h^{2} \frac{d u}{d \varphi}+\varepsilon f_{\varphi} .
\]

Переходя к переменному $\varphi$, найдем
\[
\frac{d}{d t}(u h) \equiv\left[\frac{d u}{d \varphi} h+\frac{d h}{d \varphi} u\right] u^{2} h .
\]

Из этих двух выражений получим искомое уравнение
\[
\frac{d h}{d \varphi}=\frac{\varepsilon f_{\varphi}}{u^{3} h} .
\]

Используя в свою очередь уравнение (8.8), мы можем исключить из уравнения (8.7) производную $d h / d \varphi$
\[
\frac{d^{2} u}{d \varphi^{2}}+u=\frac{\mu}{h^{2}}-\varepsilon\left\{\frac{f_{r}}{u^{2} h^{2}}+\frac{d u}{d \varphi} \frac{f_{\varphi}}{u^{3} h^{2}}\right\} .
\]

Уравнения (8.8) и (8.9) образуют замкнутую систему уравнений относительно $u$ и $h$. В результате ее интегрирования мы получаем элементы орбиты как функции полярного угла ч. Для того чтобы перейти к аргументу $t$, необходимо еще вычислить квадратуру
\[
t=\int_{0}^{\varphi} \frac{d \varphi}{u^{2} h} .
\]

При $\varepsilon=0$ система уравнений (8.8), (8.9) имеет точное решение
\[
\begin{array}{c}
u(\varphi)=\frac{\mu}{h^{2}}+x \cos \varphi, \\
h=\mathrm{const}=h_{0},
\end{array}
\]

где $x$ и $h_{0}$ — некоторые постоянные.
Следуя общей схеме метода усреднений, решение системы (8.8), (8.9) будем искать в виде
\[
u=\frac{\mu}{h^{2}}+x \cos y, \quad \frac{d u}{d \varphi}=-x \sin y,
\]

где $x$ и $y$ — неизвестные функции угла $\varphi$. Из (8.12) находим
\[
\begin{array}{c}
x=\sqrt{\left(u-\frac{\mu}{h^{2}}\right)^{2}+\left(\frac{d u}{d \varphi}\right)^{2}}, \\
y=-\operatorname{arctg} \frac{d u}{d \varphi} \frac{1}{u-\frac{\mu}{h^{2}}} .
\end{array}
\]

Дифференцируя (8.13), в силу уравнений (8.8) и (8.9) получим $\frac{d x}{d \varphi}=\frac{1}{x}\left[-x \cos y \frac{d}{d \varphi}\left(\frac{\mu}{h^{2}}\right)-\frac{d u}{d \varphi} \varepsilon\left(\frac{f_{r}}{u^{2} h^{2}}+\frac{d u}{d \varphi} \frac{f_{\varphi}}{u^{3} h^{2}}\right)\right]=$
\[
=\varepsilon\left\{\frac{2 \mu f_{\varphi} \cos y}{\left(\frac{\mu}{h^{2}}+x \cos y\right)^{3} h^{4}}+\sin y\left[\frac{f_{r}}{\left(\frac{\mu}{h^{2}}+x \cos y\right)^{2} h^{2}}-\frac{x \sin y f_{\varphi}}{\left(\frac{\mu}{h^{2}}+x \cos y\right)^{3} h^{3}}\right]\right\}
\]

и окончательно
\[
\begin{aligned}
\frac{d x}{d \varphi}=\varepsilon\left\{\frac{2 \mu f_{\Phi} h^{2} \cos y}{\left(\mu+x h^{2} \cos y\right)^{3}}+\frac{f_{\tau} h^{2} \sin y}{\left(\mu+x h^{2} \cos y\right)^{2}}-\right. & \left.\frac{x f_{\varphi} h^{4} \sin ^{2} \varphi}{\left(\mu+x h^{2} \cos y\right)^{3}}\right\} \equiv \\
& \equiv \varepsilon F_{1}(x, h, y, \varphi) .
\end{aligned}
\]

Продифференцируем теперь равенство (8.14)
\[
\frac{d y}{d \varphi}=-\frac{\frac{d^{2} u}{d \varphi^{2}}\left(u-\frac{\mu}{h^{2}}\right)-\frac{d u}{d \varphi}\left(\frac{d u}{d \varphi}-\frac{d}{d \varphi}\left(\frac{\mu}{h^{2}}\right)\right)}{\left[1+\left(\frac{d u}{d \varphi} \frac{1}{u-\frac{\mu}{h^{2}}}\right)^{2}\right]\left(u-\frac{\mu}{h^{2}}\right)^{2}},
\]

откуда после замены (9.12) получим
\[
\begin{aligned}
\frac{d y}{d \varphi} & =-\frac{1}{x^{2}}\left\{x \operatorname { c o s } y \left(-x \cos y-\varepsilon\left[\frac{f_{r}}{\left(\frac{\mu}{h^{2}}+x \cos y\right)^{2} h^{2}}-\right.\right.\right. \\
& \left.\left.\left.-\frac{x \sin y f_{\varphi}}{\left(\frac{\mu}{h^{2}}+x \cos y\right)^{3} h^{2}}\right]\right)+x \sin y\left(-x \sin y-\frac{2 \mu \varepsilon f_{\varphi}}{h^{4}\left(\frac{\mu}{h^{2}}+x \cos y\right)^{3}}\right)\right\}
\end{aligned}
\]

или окончательно
\[
\frac{d y}{d \varphi}=1+\varepsilon\left\{\frac{f_{r} h^{2} \cos y}{x\left(\mu+x h^{2} \cos y\right)^{2}}-\frac{i_{\varphi} h^{4} \sin 2 \varphi}{2\left(\mu+x h^{2} \cos y\right)^{3}}+\frac{2 \mu f_{\varphi} h^{2} \sin y}{x\left(\mu+x h^{2} \cos y\right)^{3}}\right\} \equiv
\]
\[
=1+\varepsilon F_{3}(x, h, y, \varphi) \text {. }
\]

Система уравнений (8.8), (8.15) и (8.16) может быть представлена в виде
\[
\left.\begin{array}{l}
\frac{d x}{d \varphi}=\varepsilon F_{1}(x, h, y, \varphi), \\
\frac{d h}{d \varphi}=\varepsilon F_{2}(x, h, y, \varphi), \\
\frac{d y}{d \varphi}=1+\varepsilon F_{3}(x, h, y, \varphi),
\end{array}\right\}
\]

где
\[
F_{3}=\frac{f_{\varphi} h^{5}}{\left(\mu+x h^{2} \cos y\right)^{2}} .
\]

Если функции $f_{\varphi}, f_{r}$ являются периодическими функциями полярного угла $\varphi$, го система уравнений (8.17) относится к типу уравнений, изучаемому в этой главе. Она содержит две медленные переменные $x$ и $h$ и две быстрые $y$ и $\varphi$, поэтому для ее исследования могут быть использованы изложенные выше методы анализа. При изучении этой системы может встретиться вся гамма случаев, изученных в этой главе. В системе (8.17) могут возникать резонансные ситуации.

В том случае, когда возмущающие силы не зависят от полярного угла $\varphi$, система (8.17) — это система с одной вращающейся фазой. Укороченные уравнения Ван-дер-Поля можно записать для нее в следующем виде:
\[
\left.\begin{array}{l}
\frac{d x}{d \varphi}=\varepsilon \bar{F}_{1}(x, h), \\
\frac{d y}{d \varphi}=1+\varepsilon \bar{F}_{3}(x, h), \\
\frac{d h}{d \varphi}=\varepsilon \bar{F}_{2}(x, h) .
\end{array}\right\}
\]

Переход от системы (8.17) к системе (8.18) связан с операцией усреднений по быстрой переменной $y$. Заметим, что $\bar{F}_{i}$ в общем случае, когда $f_{i}$ зависит от $x, h$ и $y$, — это некоторые довольно сложные квадратуры. Таким образом, переход от системы (8.17) к (8.18) приводит к отделению медленных переменных, что позволяет при той же точности проводить интегрирование с большим шагом и, следовательно, уменьшить число шагов численного интегрирования. В то же время на каждый шаг численного интегрирования в системе (8.18) затрачивается значительно большее число машинных операций, поскольку он требует на каждом шаге вычисления квадратур. Отсюда следует, что вопрос о целесообразности перехода от одной системы к другой в каждом отдельном случае должен быть предметом специального исследования.

Примечание. Метод численного интегрирования, по идее близкий к изложенному, но не использующий непосредственно приемов теории усреднений движения, был разработан еще в 1960 г. Г. П. Таратыновой *). Метод осреднений в его стандартном виде в траекторных задачах теории искусственных спутников Земли, по-видимому, впервые был применен Лассом и Лореллом **).

Рассмотрим снова систему (8.17) и предположим, что требования точности определяют величину шага по фазовым переменным. Пусть шаг $\Delta x \sim \delta, \Delta y \sim \delta, \Delta h \sim \delta$. Тогда шаг по времени $\Delta t$, если принять простейшую разностную схему, определится из равенства
\[
\begin{aligned}
\Delta t_{1}=\min \left\{\frac{\Delta x}{\varepsilon \max \left|F_{1}(x, y, h)\right|} ; \quad \frac{\Delta h}{\left.\varepsilon \max \mid F_{2}(x, y, h)\right\}} ;\right. \\
\left.\frac{\Delta y}{1+\varepsilon \max \left|F_{3}(x, y, h)\right|}\right\}=\frac{\delta}{1+\varepsilon \max \left|F_{3}\right|} \sim \delta,
\end{aligned}
\]
*) Г. П. Таратынов а, Методы численного решения уравнения в конечных разностях и их применение к расчетам орбит искусственных спутников Земли, сб. «Искусственные спутники Земли», №4, 1960.
$\left.{ }^{* *}\right)$ H. Lass, J. Lore11, Low acceleration takcoff from a satellite orbit, J. Amer Rocket. Soc., vol. 31, № 1 (1961).

которое справедливо для достаточно малых значений $\varepsilon$. Таким образом, число шагов $N_{1}$ определится так:
\[
N_{1}=\frac{T}{\delta},
\]

где $T$ — промежуток интегрирования.
Обозначим через $k_{1}$ число машинных операций, необходимых для реализации одного шага численного интегрирования системы (8.17). Тогда $R_{1}$ — общее число операций, необходимых для вычисления $x$ и $h$, будег
\[
R_{1}=\frac{T}{\delta} k_{1} .
\]

Рассмотрим теперь систему (8.18) и поставим себе ту же задачу: вычислить $x(T)$ и $h(T)$.
Теперь величина $\Delta t_{2}$ шага численного интегрирования будет
\[
\Delta t_{2}=\frac{\delta}{\boldsymbol{E} F},
\]

где
\[
F=\max \left\{\bar{F}_{1}(x, h), \bar{F}_{2}(x, h)\right\} .
\]

Обозначая через $k_{2}$ число машинных операций, необходимых для реализации одного шага численного интегрирования системы (8.18), мы получаем такую оценку $R_{2}$ — числа операций:
\[
R_{2}=\frac{T}{\delta} \mathrm{E} k_{2} .
\]

Вопрос о целесообразности замены в процессе численного интегрирования системы (8.17) системой (8.18) может быть в дан ном случае решен сравнением чисел $R_{1}$ и $R_{2}$.

Примечания. 1. Для изучения задачи, рассмотренной в данном параграфе, мы составили уравнения для величин $x$ и $h$, которые в случае $\varepsilon=0$ являются интегралами изучаемой системы уравнений. Такие величины в небесной механике носят название оскулирующих элементов. Выбранные оскулирующие элементы были очень удобны для демонстрации процедуры применения методов осреднения. Однако уравнения для этих элементов оказались довольно сложными: при усреднении возникала проблема вычисления сложных вспомогательных квадратур. Поэтому в прикладных расчетах используют обычно другие оскулирующие элементы.
2. Заметим, что если величины возмущающих сил постоянны (или не зависят от $y$ ), то квадратуры в уравнениях (8.18) вычисляются в явном виде.