Часть методологии, представленной в этом и в следующем параграфах, была получена в сотрудничестве с Г.Гомесом, А.Джорбой и Дж. Масдемонтом по контракту с Европейским космическим агентством. Она применялась к точке $L_{1}$ ОЗТТ в случае системы Солнце-Земля и к точкам $L_{4}$ и $L_{5}$ системы Земля – Луна. Позднее она применялась к некоторым улучшенным моделям физической задачи и, наконец, к наилучшей из имеющихся моделей, которая была выведена из численных .IPI. эфемерил.

Рассмотрим пространственную ограниченную задачу трех тел (ОЗТТ), т. е. движение безмассовой частицы (космического корабля) под действием гравитации двух массивных тел (Земли и Луны). Предполагаем, что массивные тела движутся по круговым орбитам вокруг их центра масс. При использовании вращающейся системы координат с периодом одной из круговых орбит и с началом координат в центре масс, можно считать, что два массивных тела покоятся. Масштабируя массы и расстояния, можно принять сумму масс Земли и Луны и расстояние между ними равными единице. Тогда масса Луны, являющаяся параметром, будет равна $\mu_{\text {moon }}=$ $=0.012150582$.

Во вращающейся (синодической) системе координат существует пять точек равновесия (или точек либрации в классической терминологии). Нас интересует точка либрации $L_{2}$, расположенная за Луной на продолжении линии, соединяющей Землю и Луну. Точку равновесия между Землей и Луной обозначают $L_{1}$.

В пространственной ОЗТТ точки либрации $L_{1}$ и $L_{2}$ в линейном приближении ведут себя по типу седло $\times$ центр $\times$ центр. В общем случае вблизи
${ }^{1}$ ОЗТТ – ограниченная задача трех тел.

них существуют двумерные торы. Из-за проблемы малых знаменателей на этих торах частоты должны удовлетворять некоторому диофантову условию. Когда нас интересуют короткие временные промежутки, проблемы малых знаменателей не возникают, если они не связаны с резонансами малого порядка. Пусть $\omega_{1}, \omega_{2}$ – две основные частоты в точках $L_{i}, i=1,2$. Для малых значений $\mu$ частоты достаточно близки друг другу. Для предельной задачи Хила ( $\mu=0$ ) они равны $\left(28^{1 / 2}-1\right)^{1 / 2}$ и 2 соответственно. Для $\mu=\mu_{\text {moon }}$ в точке $L_{2}$ значения $\omega_{1}$ и $\omega_{2}$ близки к 1,86264588 и 1,78617616 , и после резонанса, возникающего при отношении один к одному, следующий сильный резонанс возникает лишь в 47 -ом порядке. Очень сильный резонанс возникает в 8398-ом порядке.

При изменении амплитуд орбит частоты также меняются. Известно, что плоские периодические орбиты Ляпунова, выходящие из точки либрации, бифурцируют в гало-орбиты с амплитудой $x$ ( $x$ направлена к планетам), приближенно равной 0,1845 , если за единицу расстояния принято расстояние от меньшей массы до точки либрации. В представленном случае она равна значению $\gamma=0,16783273$, умноженному на расстояние Земля-Луна. Таким образом это расстояние соответствует резонансу один к одному, приводящему к гало-орбитам. Если мы будем проводить наши исследования в области малых амплитуд, эффект этого резонанса будет не слишком сильным. Резонансы высокого порядка имеют малую амплитуду и большой период. Таким образом, мы можем попробовать формально получить инвариантные торы до некоторого порядка (меньшего 47). В конкретном примере, представленном ниже, порядок равен 45 . Обрывая выражение на этом порядке, мы получим слоение инвариантными торами в некоторой области вокруг точки равновесия. Для некоторых амплитуд торы точно не существуют, но малая ширина соответствующих стохастических областей делает этот факт несущественным для приложений на коротких временных промежутках, в нашем случае полет космического корабля происходит в течение всего нескольких лет.
2.1. Уравнения движения

Рассмотрим следующую систему координат: начало координат помещено в точку $L_{2}$, ось $x$ направлена в сторону меньшей планеты, а ось $y$ расположена в плоскости движения планет и повернута относительно оси $x$ на угол $\frac{\pi}{2}$ против часовой стрелки, и, наконец, ось $z$ дополняет их до правой системы координат. Примем расстояние $\gamma$ между малой планетой $m_{2}$ и $L_{2}$ за единицу. Единицу времени выберем так, чтобы период обращения малой планеты был равен $2 \pi$. Через $\mu$ обозначим массу малой планеты (в нашем случае – Луны), а через 1 – $\mu$ – массу бо́льшей планеты (в нашем случае Земли). Значения $\gamma$ могут быть получены как решения уравнения Эйлера пятой степени
\[
\gamma^{5}+(3-\mu) \gamma^{4}+(3-2 \mu) \gamma^{3}-\mu \gamma^{2}-2 \mu \gamma-\mu=0 .
\]

Его можно разрешить методом Ньютона, выбирая в качестве первоначальной итерации значение $(\mu / 3)^{1 / 3}$.

Пусть оси $X, Y, Z$ образуют обычную синодическую систему координат с началом в центре масс. Уравнения движения в этой системе имеют вид
\[
\ddot{X}-2 \dot{Y}=\Omega_{X}, \quad \ddot{Y}+2 \dot{X}=\Omega_{Y}, \quad \ddot{Z}=\Omega_{Z},
\]

где
\[
\Omega=\frac{1}{2}\left(X^{2}+Y^{2}\right)+(1-\mu) r_{1}^{-1}+\mu r_{2}^{-1},
\]

причем
\[
r_{1}^{2}=(X-\mu)^{2}+Y^{2}+Z^{2}, \quad r_{2}^{2}=(X-\mu+1)^{2}+Y^{2}+Z^{2} .
\]

Выполним следующее преобразование координат
\[
x=-\frac{1}{\gamma}(X-\mu+1+\gamma), \quad y=-\frac{1}{\gamma} Y, \quad z=\frac{1}{\gamma} Z .
\]

Тогда получим
\[
\begin{aligned}
r_{1}^{-1} & =(1-\gamma)^{-1}\left(1+\frac{2 \gamma}{1-\gamma} x+\frac{\gamma^{2}}{(1-\gamma)^{2}} \rho^{2}\right)^{-1 / 2}= \\
& =(1-\gamma)^{-1} \sum_{n \geqslant 0}\left(\frac{\gamma \rho}{1-\gamma}\right)^{n}(-1)^{n} P_{n}\left(\frac{x}{\rho}\right), \\
r_{2}^{-1} & =\gamma^{-1}\left(1-2 x+\rho^{2}\right)^{-1 / 2}=\gamma^{-1} \sum_{n \geqslant 0} \rho^{n} P_{n}\left(\frac{x}{\rho}\right),
\end{aligned}
\]

где $\rho^{2}=x^{2}+y^{2}+z^{2}, P_{n}$ – полиномы Лежандра, а
\[
c_{n}=\gamma^{-3}\left(\mu+(-1)^{n}(1-\mu)\left(\frac{\gamma}{1-\gamma}\right)^{n+1}\right)
\]

– вспомогательные коэффициенты, зависящие только от $\mu$. Из формулы (3) легко видеть, что в случае задачи Хилла, когда $\mu$ стремится к нулю, $c_{2}=4$ и $c_{k}=3 \times(1)^{k}$ для $k>2$. Таким образом, уравнения движения можно записать в виде
\[
\begin{aligned}
\ddot{x}-2 \dot{y}-\left(1+2 c_{2}\right) x & =\frac{\partial}{\partial x} \sum_{n \geqslant 3} c_{n} \rho^{n} P_{n}\left(\frac{x}{\rho}\right), \\
\ddot{y}+2 \dot{x}+\left(c_{2}-1\right) y & =\frac{\partial}{\partial y} \sum_{n \geqslant 3} c_{n} \rho^{n} P_{n}\left(\frac{x}{\rho}\right), \\
\ddot{z}+2 \dot{x}+\left(c_{2}-1\right) y & =\frac{\partial}{\partial z} \sum_{n \geqslant 3} c_{n} \rho^{n} P_{n}\left(\frac{x}{\rho}\right) .
\end{aligned}
\]

Напомним также, что уравнения (1) можно записать в гамильтоновой форме с гамильтонианом
\[
\begin{array}{c}
H=\frac{1}{2}\left(P_{X}^{2}+P_{Y}^{2}+P_{Z}^{2}\right)+Y P_{X}-X P_{Y}-(1-\mu) r_{1}^{-1}-\mu r_{2}^{-1}= \\
=\frac{1}{2}\left(\dot{X}^{2}+\dot{Y}^{2}+\dot{Z}^{2}\right)-\Omega,
\end{array}
\]

где введены переменные
\[
P_{X}=\dot{X}-Y, \quad P_{Y}=\dot{Y}+X, \quad P_{Z}=\dot{Z} .
\]

Уравнения (4) также можно явно записать в гамильтоновой форме. Пусть $p_{x}=\dot{x}-y, p_{y}=\dot{y}+x, p_{z}=\dot{z}$. Соответствующий гамильтониан имеет вид
\[
K=\frac{1}{2}\left(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}\right)+y p_{x}-x p_{y}-\sum_{n \geqslant 2} c_{n} \rho^{n} P_{n}\left(\frac{x}{\rho}\right) .
\]

Гамильтонианы $H$ и $K$ связаны соотношением
\[
H=K \gamma^{2}-\frac{1}{2}(1-\gamma-\mu)^{2}-\frac{\mu}{\gamma}-\frac{1-\mu}{1-\gamma} .
\]
2.2. Решение в формальных рядах

Если отбросить из правой части системы (4) члены выше второй степени (то есть исключить экспоненциально убывающие и растущие члены), то решения системы, ограничивающие центральные многообразия, имеют вид
\[
x=\alpha \cos \left(\omega, \tau+\phi_{1}\right), \quad y=-\bar{k} \alpha \sin \left(\omega_{1} \tau+\phi_{1}\right), \quad z=\beta \cos \left(\omega_{2} \tau+\phi_{2}\right), \text { (7) }
\]

где
\[
\bar{k}=\frac{\left(\omega_{1}^{2}+2 c_{2}+1\right)}{2 \omega_{1}} .
\]

Параметры $\alpha, \beta$ и $\phi_{1}, \phi_{2}$ представляют собой произвольные амплитуды и фазы соответственно. Частоты в (7) выражаются следующим образом: $\omega_{1}=$ $=\left(2-c_{2}+\left(9 c_{2}^{2}-8 c_{2}\right)^{1 / 2} / 2\right)^{1 / 2}, \omega_{2}=c_{2}^{1 / 2}$. Так как система (4) автономна, то можно выбрать начальное время таким, что $\phi_{1}=0$. Кроме того, обозначим $\omega_{1} \tau$ и $\omega_{2} \tau+\phi_{2}$ соответственно через $\theta_{1}$ и $\theta_{2}$. При рассмотрении правой части (4), согласно методу Линдштедта-Пуанкаре, необходимо исключить секулярные члены для различных частот, зависящих от амплитуд. Обозначим $\theta_{1}$ через $d \tau$, а $\theta_{2}$ – через $f \tau+\phi_{2}$, где $d=\omega_{1}+O(\alpha, \beta), f=$ $=\omega_{2}+O(\alpha, \beta)$.

Далее будем искать решения в комплексной экспоненциальной форме вида
\[
\begin{array}{c}
x=\sum x_{i, j, k, m} \alpha^{i} \beta^{j} \gamma^{k} \gamma_{2}^{m}, \quad y=\sqrt{-1} \sum y_{i, j, k, m} \alpha^{i} \beta^{j} \gamma_{1}^{k} \gamma_{2}^{m}, \\
z=\sum z_{i, j, k, m} \alpha^{i} \beta^{j} \gamma_{1}^{k} \gamma_{2}^{m},
\end{array}
\]

где $\gamma_{s}=\exp \left(i \theta_{s}\right), s=1,2$. Вследствие симметрии задачи коэффициенты в (8) удовлетворяют соотношениям
1. $i, j \geqslant 0, i+j \geqslant 1,|k| \leqslant i,|m| \leqslant j, i-k \equiv 0 \bmod 2, j-m \equiv 0 \bmod 2$.
2. $x_{i, j,-k,-m}=x_{i, j, k, m}, y_{i, j,-k,-m}=-y_{i, j, k, m}, z_{i, j,-k,-m}=z_{i, j, k, m}$. Следовательно, для вычислений достаточно оставить лишь члены с $k>0$ или, если $k=0$, с $m \geqslant 0$.
3. Для $x$ и $y$ необходимо предполагать $j \equiv 0 \bmod 2$, а для $z-j \equiv 1$ $\bmod 2$.
4. Метод Линдштедта-Пуанкаре можно использовать так, что при $(k, m)=(1,0)$ единственный член, возникающий в выражении для $x$, соответствует $(i, j)=(1,0)$. Аналогичным образом, если $(k, m)=$ $=(0,1)$ то единственный член в $z$ соответствует $(i, j)=(0,1)$.

5. В выражениях для $d$ и $f$ содержатся только члены вида $\alpha^{i} \beta^{j}$, где $i$ и $j$ являются четными.
Вычисления членов в разложении (8) производятся по возрастании общего порядка $n=i+j$. Для $n=1$ в (7) получим значения
\[
x_{1,0,1,0}=\frac{1}{2}, \quad y_{1,0,100}=\frac{\bar{k}}{2}, \quad z_{0,1,0,1}=\frac{1}{2} .
\]

Используя свойства 1) и 5), можно получить $x_{n}, y_{n}, z_{n}$ (т. е. члены порядка $n$ по $\alpha$ и $\beta$ ), а если $n$ нечетно, то еще и $d_{n-1}$ и $f_{n-1}$. Для вычисления $x_{n}, y_{n}, z_{n}$ необходимо положить $x_{i, j, 1,0}=0$, если $(i, j)
eq(1,0)$, и $z_{i, j, 0,1}=0$, если $(i, j)
eq(0,1)$ (см. свойство 4), что и позволяет получить $d_{n-1}$ и $f_{n-1}$ для нечетного $n$. Подробнее изложение рекуррентной процедуры и особенности ее практической реализации, связанные с хранением коэффициентов, можно найти в [5].

Наконец, сделаем несколько априорных замечаний относительно объема памяти и времени вычислений. Программы были структурированы таким образом, чтобы мы могли изменять несколько параметров вычислений с помощью значений
$N G=\left(j^{4}+7 j^{3}+20 j^{2}+26 j+6\right) / 6, \quad N G 1=2 j+1, \quad N G 2=\left(j^{2}+3 j+2\right) / 2$, если максимальный порядок при вычислении программы равен $n=2 j+$ +1 . Всего для вычислений требуется $20+2 n$ векторов размерности $N G$ и несколько векторов размерностей $N G 1, N G 2$. Таким образом, в целом потребуется объем памяти (при использовании арифметики с удвоенной точностью), равный
\[
\frac{\left(4 j^{5}+56 j^{4}+236 j^{3}+691 j^{2}+881 j\right)}{6}+\text { const. }
\]

При $n=35$ это выражение приближенно равно $15.5 \mathrm{Mb}$, а при $n=45$ $48 \mathrm{Mb}$. Для простоты можно считать, что при больших $n$ объем памяти, необходимый для вычислений до порядка $n$, пропорционален $n^{4}(n+23)$. Анализ программы показывает, что время вычислений пропорционально примерно $n^{8}$.
2.3. Результаты и тесты

Программа, описанная в предыдущем разделе, выполнялась для $L_{1}$ и $L_{2}$ при различных значениях $n$ и $\mu$. Результаты, представленные ниже при

Рис. 1. Проекции $2 D$-тора, соответствующего $\alpha=0,15$ и $\beta=0,40$, на плоскости $(X, Y)$ и $(X, Z)$, полученные с помощью метода Линдштедта-Пуанкаре

$\mu=\mu_{\mathrm{moon}}$, были вычислены для порядка, равного 45 . Время вычисления на компьютере НР 9000/735 равно $35^{m}$.

На рисунках 1 и 2 показаны проекции $2 D$-тора, соответствующего значениям $\alpha=0,15, \beta=0,40$, на плоскости $(X, Y),(X, Z),(Y, Z)$. Разложения, приведенные выше, могут быть использованы при получении соответствующих ляпуновских семейств «горизонтальных» и «вертикальных» орбит, когда либо $\beta=0$, либо $\alpha=0$. На рис. 3 показаны некоторые из них для нескольких значений $\alpha$ (проекции $(X, Y)$ ) и для нескольких значений $\beta$ (проекции $(X, Y)$ и $(Y, Z)$ ). На этом рисунке мы воспользовались заменой (2), чтобы представить результаты в синодических координатах. Отметим, что «вертикальные» орбиты проектирутся на плоскость $(X, Y)$ как кривые, проходящиеся дважды. Проекция на $(X, Z)$, не представленная здесь, также выглядит как дуга кривой, проходящейся дважды. Орбиты расположены приблизительно на вертикальном цилиндре.
Рис. 2. Проекция на плоскость $(Y, Z)$ предыдущего рисунка

Для полноты также покажем некоторые гало-орбиты. Эти орбиты возникают из плоского ляпуновского семейства, когда они теряют устойчивость внутри центрального многообразия точки $L_{2}$. То есть кроме начальной гиперболической пары переменных, общей для всех орбит вблизи $L_{2}$, другая пара также становится неустойчивой. Это происходит, когда вертикальные колебания вокруг плоских .яяпуновских орбит входят в резонанс один к одному с периодом орбиты. Эти орбиты можно получить с помощью процедуры Линдштедта-Пуанкаре и разложения типа (8), где теперь оба угла $\theta_{s}, s=1,2$ равны, а амплитуды $\alpha, \beta$ удовлетворяют некоторому соотношению $g(\alpha, \beta)=0$ (см. [3]). На рис. 3 изображены проекции этих орбит на плоскости $(X, Y),(X, Z),(Y, Z)$ для нескольких значений $\beta$. Соответствующие значения $\alpha$ получаются из соотношения $g(\alpha, \beta)=0$, в показанном на рисунках интервале $\beta$ амплитуда $\alpha$ меняется лишь между 0.1853 и 0.2318 .

Аналитические решения (1) проверялись численно следующим образом: для заданных $\rho>0$ и $\theta \in\left[0, \frac{\pi}{2}\right] \alpha$ и $\beta$ вычислялись с помощью формул $\alpha=\rho \cos \theta, \beta=\rho \sin \theta$. Затем мы брали фазы $\phi_{1}$ и $\phi_{2}$ кратными $\frac{\pi}{4}$. Вследствие симметрии достаточно выбрать $\phi_{2}<\pi$. Из этих данных мы получаем начальные условия для $x, y, z$ и $\dot{x}, \dot{y}, \dot{z}$ при $t=0$ и, используя (2), находим соответствуююие синодические координаты. Затем начинается численное интегрирование и через определенные интервалы времени вплоть до $\pi$ вычисляются расстояния (в фазовом пространстве) между численным и аналитическим решениями. Это время приближенно равно периоду горизонтальных и вертикальных колебаний. Затем при заданном $\theta$ получаем значения $\rho$ и $\rho_{\max }\left(\rho<\rho_{\operatorname{m\varepsilon x}}\right.$ ) такие, чтобы максимальная ошибка во все рассматриваемые моменты времени для всех начальных фаз была меньше, чем $10^{-6}$. Принимая во внимание, что неустойчивая компонента на этом временном интервале возрастает более, чем в 1000 раз, начальные ошибки физически имеют порядок лишь несколько метров. Значения $\rho_{\max }$ в зависимости от $\theta$ на плоскости $(\alpha, \beta)$ показаны на рис. 4 (внутренняя кривая).

На том же самом рисунке показаны еще две кривые, ограничивающие области на плоскости $(\alpha, \beta)$. При заданном $\theta$ в выражения для $d$ и $f$ можно подставить $\alpha$ и $\beta$ как функции $\rho$, при этом в разложении возникают только четные члены. «Практический радиус сходимости» (несмотря на то, что ряды не сходятся) для разложений вплоть до порядка $n=45$ можно оценить различными способами. Один из методов заключается в изучении поведения формы $A \times B^{m}$ коэффициента при степени $m$ в выражении для $\rho$,

Рис. 3. Некоторые проекции периодических орбит, полученные с помощью разложения (8). (а) Проекции плоских орбит на ( $X, Y$ ) при $\alpha=$ $=0.05,0.1,0.15,0.2,0.25,0.3 . \mathrm{Ha}$ (b) и (c) изображены проекции на плоскости $(X, Y)$ и $(Y, Z)$ орбит, соответствующих значениям $\beta=$ $=0.1,0.2,0.3,0.4,0.5,0.6,0.7 . \mathrm{Ha}$ (d), (e) и (f) показаны проекции на плоскости $(X, Y),(X, Z)$ и $(Y, Z)$ орбит из гало семейства для $\beta=$ $=0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4$

последующем приведении его к логарифмической шкале, нахождении значений $B$ для $d$ и $f$ и использовании обратной величины максимума этих значений в качестве радиуса сходимости. Результатом применения этого метода является внешняя кривая на рис. 4. Другая возможность – это выбор $\left|t_{44}\right|<10^{-3}$, где $t_{44}$ обозначает наибольший из последних членов в $d$ и $f$ для заданного значения $\rho$. Так возникает промежуточная кривая, расположенная между двумя предыдущими.

Рис. 4. Три кривые, проходящие от оси $\alpha$ к оси $\beta$ определяют несколько оценок области «практической сходимости». Внутренняя кривая дает очень хорошие результаты при прямой численной проверке. Два семейства кривых, проходящих из правого нижнего угла в левый верхний соответствуют постоянным значениям разности (непрерывной линии) и частного (пунктирной линии) частот

На том же рисунке также изображены два семейства кривых. Они соответствуют постоянным значениям разности и частного $d$ и $f$. Разность меняется от 0.005 до 0.110 с шагом 0.005 . При увеличении разности кривые переходят из правого нижнего угла в левый верхний угол. Частное меняется в диапазоне от 1.0025 до 1.0625 с шагом 0.0025. Дальнейшее изучение и возможность применения всех этих семейств орбит к полетам космических кораблей можно найти в [3, 4, 5] и [6].
2.4. Редукция на центральное многообразие вблизи $L_{2}$

Мы хотим понять, как организованы орбиты, вычисленные в предыдущих параграфах, чтобы получить общее описание динамики в окрестности $L_{2}$.

Сначала уберем неустойчивые члены гамильтониана. Это достигается редукцией на центральное многообразие $W^{c}$ размерности 4. Опишем, как происходит вычисление. После линейной замены переменных, переходя от $x, y, p_{x}, p_{y}$ к $x_{1}, x_{2}, y_{1}, y_{2}$ и полагая $x_{3}=z$ и $y_{3}=p_{z}$, получим, что гамильтониан $K$, заданный формулой (6), можно записать в виде
\[
M=\lambda x_{1} y_{1}+\frac{1}{2} \omega_{1}\left(x_{2}^{2}+y_{2}^{2}\right)+\frac{1}{2} \omega_{2}\left(x_{3}^{2}+y_{3}^{2}\right)+\sum_{j \geqslant 3} M_{j}\left(x_{1}, x_{2}, x_{3}, y_{1}, y_{2}, y_{3}\right),
\]

где $M_{j}$ – однородный полином степени $j$.
Затем мы переходим к промежуточной нормальной форме, пытаясь уничтожить в $M$ все члены, у которых общая степень по переменным $\left(x_{1}, y_{1}\right)$ равна единице. Это достигается использованием канонического преобразования с подходящей порождающей функцией (которую предстоит определить). Данная процедура не вызывает затруднений, поскольку здесь нет проблемы малых знаменателей. Знаменатели, возникающие при определении порождающей функции (при использовании некоторых комплексификаций для уменьшения объема вычислений), имеют вид:
\[
\left(k_{1}-l_{1}\right) \lambda+\sqrt{-1}\left(k_{2}-l_{2}\right) \omega_{1}+\left(k_{2}-l_{3}\right) \omega_{2}
\]

с ограниченным снизу модулем $|\lambda|$. Пусть $\widetilde{M}$ – гамильтониан этой промежуточной нормальной формы. С помощью перестановки слагаемых его можно записать в виде
\[
\widetilde{M}=M^{0}\left(x_{2}, x_{3}, y_{2}, y_{3}\right)+\sum_{j_{1}+j_{2}>1}\left(x_{1}^{j_{1}} y_{1}^{j_{2}}\right) M^{j_{1} j_{2}}\left(x_{2}, x_{3}, y_{2}, y_{3}\right) .
\]

Следует иметь в виду, что новые переменные отличаются от предыдущих, несмотря на то, что мы сохранили прежние обозначения. Преобразование является тем более близким к тождественному, чем ближе траектория к $L_{2}$.

Рис. 5. Сечение Пуанкаре гамильтонова потока, редуцированного на центральном многообразии точки $L_{2}$ в случае Земля-Луна на заданном уровне энергии

То, что проблема малых знаменателей не возникла, ни в коей мере не означает, что процедура является сходящейся. В общем случае получается только $C^{\infty}$-многообразие. Норма коэффициентов порядка $n$ возрастает как $n !$.

Из формы нового гамильтониана видно, что $I_{1}=x_{1} y_{1}$ – первый интеграл. Проведем редукцию на центральное многообразии, положив $I_{1}=0$. Таким образом остается гамильтониан с двумя степенями свободы
\[
M^{0}=\frac{1}{2} \omega_{1}\left(x_{2}^{2}+y_{2}^{2}\right)+\frac{1}{2} \omega_{2}\left(x_{3}^{2}+y_{3}^{2}\right)+\sum_{j \geqslant 3} M_{j}^{0}\left(x_{2}, x_{3}, y_{2}, y_{3}\right) .
\]

Отметим, что если в (9) положить $x_{3}=y_{3}=0$, мы получим инвариантное множество относительно потока, связанного с $M^{0}$. То есть мы приходим к гамильтониану с одной степенью свободы, орбиты которого являются просто ляпуновскими периодическими орбитами вблизи $L_{1}$.

Фиксируя $M^{0}=h$, мы можем использовать $x_{3}=0$ при $y_{3}>0$ в качестве поверхности сечения. Это сечение, $\Sigma$, не является глобальным, поскольку оно полностью содержит соответствующую ляпуновскую орбиту. Однако можно считать, что все точки рассматриваемой траектории (в пределах сечения) определены. Таким образом мы получим сферу $S^{2}$ (классическое слоение Хопфа), которую можно рассматривать как пространство «оскулирующих траекторий». Для некоторого значения $h$, например $h_{H}$, уровень $M^{0}=h$ содержит ляпуновскую орбиту, бифурцирующую в гало-орбиты. Начиная с этого значения энергии, сечение $\sigma$ содержит две неподвижные точки, связанные с двумя (симметричными) гало-орбитами. Известно, что по крайней мере для средних значений $h$ эти орбиты эллиптического типа. Для значения $h$, не слишком далеких от $h_{N}$, имеем на $S^{2}$ четыре неподвижные точки. Две из них соответствуют ляпуновской орбите (гиперболической) и почти вертикальной периодической траектории (эллиптической) соответственно. Другие две являются гало-орбитами (эллиптическими).

Программа, выполняющая эту редукцию, была разработана, реализована и выполнена вплоть до общего порядка, равного 16. Пусть $\xi$ и $\eta$ обозначают новые переменные после опускания гиперболических слагаемых в промежуточной нормальной форме. После получения $M^{0}(\xi, \eta)$ мы можем провести моделирование. На рис. 5 показаны результаты этого моделирования для такого значения энергии, при котором явно видны гало-орбиты. Рисунок соответствует сечению Пуанкаре $\xi_{2}=0$ (в исходных переменных сечение лежит вблизи $z=0$ ). Почти в центре рисунка видна неподвижная точка, соответствующая «вертикальной» периодической орбите. Она окружена инвариантными торами, которые принадлежат семейству, найденному в начале этого параграфа, несмотря на то, что в данном случае амплитуды слишком велики и не могут быть хорошо аппроксимированы с помощью результатов метода Линштедта-Пуанкаре. Внешняя кривая является ляпуновской плоской траекторией, расположенной на этом уровне энергии. Две другие неподвижные точки соответствуют двум гало-орбитам, которые симметричны друг другу относительно $z=0$. Они, в свою очередь, также окружены инвариантными $2 D$-торами. Между $2 D$-торами вокруг вертикальной орбиты и $2 D$-торами вокруг гало-орбит расположены следы устойчивого и неустойчивого многообразий плоской ляпуновской орбиты. Однако, несмотря на то, что значение энергии отличается от энергии в точке $L_{2}$, даже стохастические области, связанные с этими многообразиями, сложно увидеть.

Отметим, что рис. 5 изображен в масштабе, отличающемся от масштаба, использованного в рис. 3. Переменные соответствуют промежуточной нормальной форме. Луна расположена вне рисунка на отрицательной вертикальной оси рис. 5 .