Кроме коллинеарных точек либрации, в ОЗТТ также существуют лагранжевы или треугольные точки либрации, расположенные вне линии, проходящей через массивные тела. В обозначениях предыдущего параграфа в синодической системе они покоятся в точках $\left(\frac{\mu-1}{2}, \frac{ \pm \sqrt{3}}{2}, 0\right)$. Все эти точки либрации (как коллинеарные, так и треугольные) возникают также в общей задаче трех тел. Одна из классических задач небесной механики связана с устойчивостью треугольных точек. Приведем здесь некоторые результаты
– Треугольные точки линейно устойчивы для массового параметра $\mu \in$ $\left[0, \mu_{1}\right]$, где
\[
\mu_{s}=\frac{\omega_{\text {short }}}{\omega_{\text {long }}}=s, \quad s \in \mathbb{N},
\]
a $\omega_{\text {short }}$ и $\omega_{\text {long }}$ – частоты в $L_{4,5}$ в плоском случае [15]. Их значения равны
\[
\left[\frac{1}{2}\left(1 \pm(1-27 \mu(1-\mu))^{1 / 2}\right)\right]^{1 / 2} .
\]

Следовательно, $\mu=\left(1-(23 / 27)^{1 / 2}\right) / 2 \approx 0.038521$.
– В плоском случае [11] треугольные точки устойчивы по Ляпунову при $\mu \in\left[0, \mu_{1}\right] \backslash\left\{\mu_{2}, \mu_{3}\right\}$.
– Относительно устойчивости по Ляпунову в пространственной задаче ничего неизвестно. Возможно в этом случае происходит диффузия Арнольда.
– Существует область «практической устойчивости», которая не слишком мала. Приведем здесь кратко то, как был получен этот результат [2], $[12]$.

Пусть $H=H_{2}+H_{3}+\ldots-$ разложение гамильтониана по степеням вблизи $L_{4,5}$, где $H_{k}$ содержит однородные члены степени $k$ по координатам и импульсам. С помощью линейного преобразования $H_{2}$ можно свести к виду
\[
-\frac{1}{2} \omega_{s}\left(q_{1}^{2}+p_{1}^{2}\right)+\frac{1}{2} \omega_{l}\left(q_{2}^{2}+p_{2}^{2}\right)+\frac{1}{2} \omega_{z}\left(q_{3}^{2}+p_{3}^{2}\right),
\]

где $\omega_{s}, \omega_{l}$ обозначают $\omega_{\text {short }}$, $\omega_{\text {long }}$, а $\omega_{z}$ – частота в вертикальном направлении ( $\omega_{z}=1$ ).
Каноническое преобразование $(q, p) \rightarrow(Q, P)$, полученное с помощью порождающей функции $G=G_{3}+G_{4}+\ldots+G_{n}$, применяется для того, чтобы привести $H$ к нормальной форме вплоть до порядка $n$. С его помощью получаем новый гамильтониан
\[
\mathscr{H}=N_{2}+N_{3}+\ldots+N_{n}+\mathscr{R}_{n+1},
\]

где $N_{k}$ – члены степени $k$ в нормальной форме и $N_{2}=H_{2}$. Если взять только $N_{2}+N_{3}+\ldots+N_{n}$, система будет интегрируема. Пусть
\[
I_{j}=\frac{1}{2}\left(Q_{j}^{2}+P_{j}^{2}\right),
\]

новые импульсы. Тогда диффузия по импульсам происходит из-за вклада
\[
\dot{I}_{j}=\left\{I_{j}, R_{n+1}\right\},
\]

где $\{\cdot, \cdot\}$ обозначают скобку Пуассона. Так как остаток $\mathscr{R}_{n+1}$ имеет вид $R_{n+1}+R_{n+2}+\ldots$, введем некоторую норму $\left\|R_{k}\right\|$. Тогда можно получить последовательность ограничений $\left\|H_{k}\right\|, k \geqslant 2,\left\|G_{k}\right\|, 2<$ $k \leqslant n$ (возникает лишь конечное число малых знаменателей) и $\left\|R_{k}\right\|$, $k>n$. Ограничения $\left\|R_{k}\right\|$, задаются в виде рекурсии, которая зависит лишь от нормы однородных частей исходного гамильтониана $\left\|H_{k}\right\|$ и от текущих малых знаменателей, которые появляются до $n$-го порядка. В шаре радиуса $\rho$ в переменных $(Q, P)$ имеем
\[
\left|\mathscr{R}_{n+1}\right|<\sum_{k>n} \rho^{k}\left\|R_{k}\right\|,
\]

где $|\cdot|$ обозначает супремум нормы. Аналогично мы можем ограничить скорость диффузии $\left|\dot{I}_{j}\right|$.

Для заданных $T, \delta$ существует начальный радиус $\rho_{0}$ такой, что если $(Q, P)_{t=0} \in B_{\rho_{0}}$, тогда $(Q, P)_{t} \in B_{\rho_{0}(1+\delta)}$ для всех $|t|<T$, где $B_{\rho}$ обозначает шар радиуса $\rho$ с центром в начале координат. Отметим, что можно получить лучшие результаты, если $H_{k}, N_{k}, G_{k}$ явно вычислены аналитически вплоть до некоторого порядка для выбранного значения $\mu$. В общем случае в качестве результата получается оценка типа Нехорошева (т. е. для фиксированного $\delta$ имеем $T \approx \exp \left(c / \rho^{d}\right)$, где $c$ и $d$ – некоторые положительные константы).
– Численное моделирование ([5], [6]) показывает «область устойчивости» даже бо́льшую, чем описанную в плоском случае [10].
Кроме вышеизложенного, также известны некоторые частные результаты для эллиптической ограниченной задачи трех тел (т. е. массивные тела движутся по эллиптическим траекториям), связанные с локальным анализом «практическая устойчивость» ([9]) и с численным определением «области устойчивости» ([1], [6]). Для реальных задач ситуация может быть более сложной, например, в случае треугольных точек системы Земля – Луна наиболее важные отклонения от модели ОЗТТ связаны с влиянием Солнца. Было установлено, что эти точки неустойчивы, но на некотором расстоянии от них существуют области практической устойчивости (по крайней мере для временных интервалов, намного превышающих продолжительность полета космического корабля), расположенные на расстоянии от $1 / 4$ до $3 / 4$ расстояния Земля – Луна ([6], [13]).

В этом параграфе мы исследуем модель ОЗТТ при $\mu=0.0002$. Существует несколько причин для выбора такого значения. С одной стороны оно достаточно мало, чтобы считать, что теория возмущений может дать хорошую аппроксимацию. С другой стороны, оно достаточно велико, чтобы не требовалось очень большого количества времени для обнаружения некоторых уходов. И наконец, она также близка к действительному соотношению масс в системе Сатурн- Титан.
3.1. Границы области практической устойчивости

Вблизи треугольных точек систему можно рассматривать как совокупность трех гармонических осцилляторов с частотами, равными $\omega_{\text {short }}$, $\omega_{\text {long }}$ и $\omega_{\text {vert }}$, последняя из которых равна единице. При малых $\mu$ значения $\omega_{\text {short }}, \omega_{\text {long }}$ соответственно равны $1-27 \mu / 8+O\left(\mu^{2}\right)$ и $(27 \mu / 4)^{1 / 2}+O\left(\mu^{3 / 2}\right)$. Следовательно, в качестве базиса в окрестности $L_{5}$ можно рассматривать частоты $\omega_{\text {vert }}, \omega_{\text {long }}$ и $\omega_{\text {vert }}-\omega_{\text {short }}$, значение которых имеет порядок $1, \mu^{1 / 2}$ и $\mu$ соответственно.

В случае, если частоты удовлетворяют некоторым условиям (не слишком близки к сильному резонансу), существует область «практической устойчивости». Предыдущее обсуждение порядков величин показывает, что если $\mu$ мало, никаких сильных резонансов не возникает, за исключением резонанса один к одному между коротким периодом и вертикальными колебаниями. Следовательно, потребуется чрезвычайно большой промежуток времени, чтобы увеличить действие на значительную величину. Нас интересует форма этих областей, поскольку результаты работ [2] и [12] дают достаточно пессимистичные оценки.

Рис. 6. Проекция некоторых сечений в области глобальной устойчивости вокруг $L_{5}$ для $\mu=0.0002$

Приведем теперь описание следующего эксперимент. Для заданных значений $\rho, \alpha$ и $Z$ рассмотрим начальную точку с синодическими координатами
\[
X=\mu+(1+\rho) \cos (2 \pi \alpha), \quad Y=(1+\rho) \sin (2 \pi \alpha), \quad Z,
\]

и с нулевой начальной скоростью $\dot{X}=0, \dot{Y}=0, \dot{Z}=0$. Точка $L_{5}$ соответствует $\rho=0, \alpha=1 / 3, Z=0$. Затем начнем численное интегрирование вплоть до заданного конечного времени. Точки вблизи $L_{5}$ просто продолжают двигаться вокруг нее (обычно на $3 D$-торе), но если начальные условия достаточно далеки от точки либрации, проекция этой траектории на плоскость $(X, Y)$ пересекает ось $X$ и затем ее можно рассматривать, как убежавшую из «большой» окрестности $L_{5}$. Ввиду формы некоторых больших торов, которые, похоже, существуют даже для больших значений $Z$, удобно принять в качестве критерия убегания из окрестности $L_{5}$ условия $Y(t)<Y^{*}$ для некоторого отрицательного значения $Y^{*}$. В качестве практического значения было взято $Y^{*}=-0.5$. Для рассматриваемого значения $\mu$ мы нашли траектории, лежащие на торах «вокруг» $L_{5}$ такие, что $Y(t)$ достигает значения вплоть до -0.233 . Промежуток времени, использованный для вычислений равен 10000 периодам обращения массивных тел. В случае системы Сатурн – Титан это означает более четырехсот лет. Такой промежуток может показаться коротким, однако посмотрим комментарии ниже.

Рис. 6 отображает часть результатов проведенного эксперимента. На нем приведена проекция на плоскость $(X, Y)$ тех точек, которые не убегают согласно принятому критерию, стартуя с плоскостей $Z=$ $=0.0,0.3,0.4,0.5,0.6,0.7$ и 0.8 . Орбиты с начальным $Z=0.1$ и 0.2 накладываются на проекции случая $Z=0.0$. Размеры шагов по $\rho$ и $\alpha$ были приняты равными $10^{-5}$ и $10^{-3}$ соответственно.

Следует сделать несколько замечаний, относительно полученных результатов.
– Для малых значений $\mu$, как в рассмотренном случае, переход из «устойчивой» области к неустойчивости достаточно резок. Это означает, что, например, увеличение от 10000 итераций до 100000 приводит к очень небольшим изменениям. С другой стороны, бо́льшая часть «убегающих» точек, уходит из окрестности $L_{5}$ за несколько итераций. Бо́льшая часть точек, которые убегают после более длительного времени, уходят за время порядка $1 / \mu$, что связано с наименьшей из основных частот $\omega_{\text {vert }}-\omega_{\text {short }}$.
– Для плоского случая неубегающие точки заполняют область, которую мы хотим определить. В этом случае возьмем кривую нулевой скорости (к.н.с), проходящую через коллинеарную точку $L_{3}$, расположенную за бо́льшим телом на расстоянии, довольно близком к расстоянию между массивными телами. Она соответствует точкам на плоскости $(X, Y)$, в которых частица, лишенная массы, покоится при данном значении энергии. В нашем случае это та же энергия, что и в точке $L_{3}$. Верхняя часть этой кривой (при $Y>0$ ) проходит по угловой переменной от $\alpha=0$ до значения, которое стремится к $\alpha_{\lim }=\cos ^{-1}\left(-\left(2^{1 / 2}-\right.\right.$ $-0.5) / \pi) \approx 156.09$ градусов, при $\mu \rightarrow 0$. Рассмотрим кривую $\rho=0$ (т.е. на единичном расстоянии от бо́льшего тела). К.н.с. расположена на расстоянии $O\left(\mu^{1 / 2}\right)$ от нее. Если говорить более точно, существует функция $K(\alpha)$, имеющая двойной ноль при $\alpha=0$ и простой ноль при $\alpha=\alpha_{\text {lim }}$, такая, что к.н.с описывается с помощью равенства $\rho= \pm(K(\alpha) \mu)^{1 / 2}+O(\mu)$. Граница области устойчивости на плоскости $Z=0$ имеет вид $\rho= \pm \frac{1}{2}(K(\alpha) \mu)^{1 / 2}+O(\mu)$.
– Существует численное свидетельство того, что в плоском случае граница области «устойчивости» связана с устойчивым и неустойчивым многообразием центрального устойчивого многообразия в $L_{3}$ (см. [6]). Это многообразие размерности 1, в частности, оно содержит устойчивое и неустойчивое многообразия семейства плоских ляпуновских орбит вокруг $L_{3}$.
– Полное множество точек в области устойчивости с нулевой синодической скоростью в пространственной ОЗТТ, похоже, является тонкой оболочкой. Форма оболочки, круговая в горизонтальном направлении и параболическая в вертикальном. Разрезая зону «устойчивости» значением $\alpha=1 / 3$, центральная точка, как функция от $Z$, близка к $\rho=-$ $-0.245 Z^{2}$.
– Как мы покажем в следующем параграфе, существуют семейства неустойчивых $2 D$-торов, которые также играют роль границы области устойчивости. Если начальное значение $Z$ мало, эти торы достигают окрестности точки $L_{3}$. Их сечения плоскостью $Z=0$ имеют положительные значения по $X$. Для больших начальных значений $Z$ сечение плоскостью $Z=02 D$-торов, которые, вероятно, находятся на границе области устойчивости, имеют отрицательные значения $X$.
Подводя итоги численных экспериментов, в плоском случае получаем, что бо́льшая часть устойчивых движений происходит на $2 D$-торах и они «ограничены» плоскими периодическими орбитами ляпуновского семейства вокруг $L_{3}$. В пространственном случае бо́льшая часть устойчивых движений происходит на $3 D$-торах и они «ограничены» $2 D$-торами, связанными с резонансом один к одному между коротким периодом и вертикальными частотами. Эти $2 D$-торы расположены по-разному, в основном в зависимости от начального значения $Z$. Разумеется, резонансы высокого порядка играют важную роль и приводят к появлению цепей гетероклинических связей, уводящих точки, вероятно, из любой окрестности любой устойчивой точки, но это требует очень большого количества времени и этот механизм можно увидеть лишь для достаточно больших значений $\mu$ (скажем, порядка $10^{-2}$ ).
3.2. Неустойчивые $2 D$-торы: обнаружение и численные эксперименты

На рис. 7 показаны первые сто тысяч итераций сечения Пуанкаре плоскостью $Z=0$ орбит, начинающихся в достаточно близких точках. Первая была вычислена при $\rho=-0.1506066340, \alpha=1 / 3, z=0.8$ и с нулевой синодической скоростью. Вторая отличается лишь значением $\rho$, которое равно -0.1506066339 . Рисунок показывает не только чувствительность траекторий к начальным условиям, но также существование $2 D$-тора, который разделяет оба вида движения.

Чтобы сделать последнее утверждение более ясным, на рис. 8 на одном том же графике показаны первые итерации отображения Пуанкаре, начинающиеся при обоих начальных условиях. После прохождения вблизи гиперболической инвариантной кривой (пересечение $2 D$-тора с $Z=0$ ) орбиты расходятся экспоненциальным образом. Стоит отметить, что движение точки, не лежащей на $3 D$-торе приближается последовательно к нескольким другим $2 D$-торам. При рассмотрении орбиты в подходящие временные интервалы кажется, что она достаточно близка к некоторому $3 D$-тору, однако время от времени она меняется (обычно очень быстро). Это явление свидетельствует о существовании гомоклинического/гетероклинического пересечения в пространственной ОЗТТ.

Вычислим теперь неустойчивые $2 D$-торы. Зафиксируем энергию и используем сечение Пуанкаре плоскостью $Z=0$. Таким образом, задача сводится к вычислению инвариантной гиперболической кривой четырехмерного симплектического преобразования. Т. к. в общем случае существует канторово однопараметрическое семейство таких объектов, мы также должны зафиксировать некоторое число вращений. Чтобы выбрать значение энергии и число вращений, будем действовать следующим образом.

Энергия выбирается (несколько произвольным образом) как среднее энергий двух рассмотренных выше орбит. Для выбора числа вращений сначала изучим итерации двух предыдущих орбит, которые остаются вблизи искомой кривой. Рассмотрим последовательные максимумы и минимумы по любой из текущих переменных ( $X, Y, \dot{X}, \dot{Y}$, значение $\dot{Z}$ выводится из соотношения для энергии (5)), затем с помощью использования интерполяции различных данных в зависимости от числа итераций получим число вращения $
u \approx 0.02189307$, которое мы будем использовать впоследствии. Также важно и обратное к этому числу $N=1 /
u \approx 45.67656$.

Вычисления периодической орбиты (устойчивой или нет) для отображения являются, в теории, простой задачей. Предположим, что мы обозначили отображение как $T$ и ищем периодическую точку с минимальным периодом $k$. Для точки орбиты $p$ можно записать уравнение $T^{k} p-p=$ $=0$ и попытаться его решить, используя, например, метод Ньютона. Для успеха нам нужно наложить некоторые условия на свойства этой орбиты, а начальные данные должны быть достаточно близки к $p$. Далее некоторые технические сложности могут возникнуть, например, если какие-либо собственные значения $D\left(T^{k}\right)(p)$ очень велики. В этом случае общая и очень эффективная процедура разрешения этого затруднения основана на методе параллельной пристрелки (см., например, [14], а также [3]).

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

Поскольку мы не имеем обратного отображения как в случае периодических орбит, его можно синтезировать. Будем искать неподвижную точку $\mathscr{R}(p)$ обратного отображения $\mathscr{R}$, где $\mathscr{R}=T^{N}, N=45.67656 . \mathrm{Pa}-$ зумеется, мы не можем вычислить степень 45.67656 отображения Пуанкаре, но можно вычислить степени в промежутке $\left[N_{1}, N_{2}\right]$, содержащем искомое значение $N$, вычислить $\mathscr{R}$, а также его дифференциал с помощью интерполяции. Однако эта численная задачи не имеет единственного решения. Действительно, все точки на траектории будут неподвижными точками синтезированного обратного отображения. Для определения единственного решения мы можем зафиксировать некоторые координаты точки $p$.

Рис. 7. Первые 100000 итераций отображения Пуанкаре $Z=0$, начиная с нулевой синодической скорости при $\alpha=1 / 3, Z=0.8$. Сверху начальное $\rho=$ $=-0.1506066340$, внизу начальное $\rho=-0.1506066339$

Рис. 8. Первые 5000 итераций отображения Пуанкаре на $3 D$-торе и первые 4165 итераций убегающей точки. См. начальные данные в подписи к рис. 7. Обе траектории остаются в течение почти 2000 итераций вблизи инвариантной кривой отображения Пуанкаре

В рассматриваемой задаче мы зафиксировали значение $Y=Y^{*}=$ $=0.3897977$. Выбор был произвольным, поскольку можно выбирать любое значение, не слишком близкое к экстремальным значениям $Y$ на кривой. Значение, которое было взято, соответствует одной из итераций отображения Пуанкаре для точки, расположенной на $3 D$-торе. Эта итерация будет взята в качестве начальной точки метода Ньютона. Лишь три координаты $(X, \dot{X}, \dot{Y}$ ) остались свободными, и они успешно вычисляются. Заметим, что координата $Y$ обратного отображения не обязана совпадать с $Y^{*}$. Совпадение должно произойти автоматически, и это было использовано в качестве проверки вычислений. Метод Ньютона останавливается, когда ошибка возвращения меньше $10^{-13}$. Такая точность получена с помощью четырех итераций процедуры. Проверка начального и конечного значения $Y$ также приводит к разности, меньшей $10^{-13}$. Результаты вычислений представлены на рис. 9.

Дифференциал обратного отображения $D \mathscr{R}$ также дает оценку гиперболичности инвариантной кривой, главное собственное значение $\lambda_{\max } \approx$ 1.24658 , откуда следует, что главное собственное значение отображения Пуанкаре вдоль кривой $\lambda_{\operatorname{map}} \approx 1.004821$. Это очень малая неустойчивость, что делает численные вычисления намного проще.

Следует отметить, что можно использовать процедуру продолжения для вычисления других инвариантных кривых либо с помощью изменения энергии, числа вращения либо их обоих. Когда число вращения принимает рациональное значение (а в компьютерной аппроксимации это происходит всегда!), инвариантной кривой обычно не существует и она, вероятно, заменяется парой периодических орбит эллиптическо-гиперболического или гиперболическо-гиперболического типа. Следовательно, не существует причины, по которой какая-либо точка этих периодических орбит будет находиться на гиперплоскости $Y=Y^{*}$. Однако численный метод по-прежнему может привести к некоторому результату без каких либо проблем. Я пробовал провести такие вычисления для чисел вращения $
u=1 / 45$ и $
u=1 / 46$ и не встретил каких-либо затруднений. Причина, очевидно, заключается в том, что обратное отображение (в данном случае $T^{45}$ или $T^{46}$ ) имеет настолько медленную динамику, что она ниже допустимого значения для метода ( $\left.10^{-13}\right)$. Чтобы увидеть, что метод неверен, в этом случае достаточно провести все вычисления с более высокой степенью точности.
3.3. Обобщение: численное получение инвариантных торов.

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

Сначала переформулируем вышеприведенную процедуру. Пусть $\mathbb{T}=$ $=\mathbb{R} / \mathbb{Z}$ – стандартный $1 D$-тор, параметризованный значениями, меняющимися в диапазоне $[-0.5,0.5]$. Отображение на инвариантной кривой сопряжено отображению $x \rightarrow x+
u$, и, следовательно, после $k$ итераций текущее значение $x$ равно $x+k
u(\bmod 1)$. Выбирая несколько значений $k$ (в конкретном примере, приведенном выше, мы выбрали $k$ от 38 до 54), получим несколько значений $x$ вблизи нуля. А затем с помощью интерполяции можно вычислить обратное отображение $\mathscr{R}$.
В общем случае мы ищем $m$-мерный тор, на котором динамика сопря-

Рис. 9. Инвариантная гиперболическая кривая отображения Пуанкаре $Z=0$ ОЗТТ, изображенная на рис. 8. (a) – проекция кривой на $X Y$. (b), (c), (d) – проекции точек соответствующего $2 D$-тора на плоскости $(X, Y),(X, Z)$ и $(Y, Z)$

жена отображению $x \rightarrow x+
u$ стандартного тора $\mathbb{T}^{m}=\mathbb{R}^{m} / \mathbb{Z}^{m}$, а $
u$ из $\mathbb{R}^{n}$ удовлетворяет условию нерезонансности: $(k,
u) \in \mathbb{Z}$ для любого $k \in \mathbb{Z}^{m}$ влечет за собой $k=0$, здесь $(\cdot, \cdot)$ обозначает скалярное произведение в $\mathbb{R}^{n}$. Мы предполагаем, что все ограничения, накладываемые первыми интегралами, были убраны, и мы сразу рассматриваем редуцированную задачу. Поскольку нас интересует $m$-мерный тор, зафиксируем $m$ координат искомой точки $p$. Далее предположим, что общая размерность $n$ по меньшей мере $2 m$, и что изменение остальных $n-m$ координат приводит к невырожденному изменению частот $
u$ при условии, что они существуют (вероятно, они существуют только на каком-то канторовом множестве). Невырожденность понимается в том смысле, что дифференциал частот относительно оставшихся переменных имеет максимальный ранг. Также предполагается, что оставшиеся $n-2 m$ координат, если такие имеются, связаны с собственными числами, отличными от единицы. Таким образом, мы неявно предполагаем, что динамика вдоль тора сводится к динамике с постоянными коэффициентами (см. [7] и [8]).

Затем рассматриваются различные итерации $T$ с числами итераций $k_{1}, k_{2}, \ldots, k_{q}$ такими, что на стандартном торе соответствующие значения $x\left(k_{1}
u, \ldots, k_{q}
u(\bmod 1)\right)$ попадают в заданную окрестность нуля. После этого можно синтезировать обратное отображение $\mathscr{R}$ как значение в нуле некоторой функции, интерполированной по значениям $T^{k_{j}}(p)$ в точках $x_{j}^{
u}, j=1, \ldots, q$. Условие, накладываемое на $p$, заключается в том, что $n-m$ свободных координат должны удовлетворять $\mathscr{R}(p)=p$, а оставшиеся $m$ координат используются как дополните.тьная проверка. Уравнение может быть разрешено методом Ньютона, необходимый дифференциал также получается с помощью интерполяции. И наконец, из дифференциала обратного отображения можно получить собственные значения $T$ по оставшимся $n-$ $-2 m$ переменным.

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

Благодарности

Эта работа частично поддерживалась грантом РB94-0215 (Испания). Вычислительные средства были представлены по грантам ESP91-0403 (Испания) CIRIT GRQ93-1135 (Каталония). Также была оказана поддержка Европейским Космическим Агентством European Network Grant ERBCHRXCT 940460 и European INTAS grant 93-0339. Кроме того, я благодарю Г. Гомеса, А. Джорбу и Дж. Масдемонта за многочисленные полезные обсуждения.