Перед нами возникает необходимость решить несколько проблем: приближение хореографии, уточнение, продолжение при изменении степени $a$, а также вычисление достаточно хороших приближений отображения Пуанкаре вблизи периодического решения. Фактически, последняя задача была решена только для восьмерки при $N=3$. Инструменты продолжения являются стандартными. Описание анализа бифуркаций и его приложений я рекомендую прочитать в [9] (см. также статью в настоящем сборнике К. Симо, Т. Стучи «Центральные устойчивые/неустойчивые многообразия и разрушение КАМ-торов в плоской задаче Хилла»).

7.1. Реализация вариационного метода

Функцию $q$ можно аппроксимировать с помощью функции $\hat{q}$, являющейся тригонометрическим полиномом со значениями на множестве равноотстоящих точек либо с помощью некоторого другого метода. Обозначим $\mathcal{P} \in R^{M}$ конечное множество параметров, необходимых для приближения. Тогда действие $A$ как в (3) аппроксимируется с помощью дискретного отображения $\hat{A}: R^{M} \rightarrow R$, где интеграл вычисляется с помощью численной квадратурной формулы. В соответствии с действием $Z / N Z$ (которое должно сохраняться при дискретизации), достаточно взять интегралы от 0 до $2 \pi / N$. Учитывая, что свободные от столкновений решения являются аналитическими, весьма удобна формула трапеции.

Затем переходим к минимизации $\hat{A}$, используя комбинацию вариантов градиентного метода. Этот метод имеет несколько
– Преимуществ. Он совершенно устойчив, и можно начать с очень приблизительных аппроксимаций типа нескольких гармоник или начерченной от руки кривой. Его легко программировать, использовать любой потенциал и т.д. Более того, градиент можно получить на основе множества $\mathcal{P}$ без необходимости численного дифференцирования. Он позволяет проводить проверку инвариантности энергии и значения остаточного ускорения: разницы между ускорением масс, вычисляемым на основе $\hat{q}$, и использованием (1).
– Неудобств. За исключением совершенно простых случаев, $M$ следует выбирать большим, особенно если прохождения осуществляются не слишком далеко от столкновения. Типичные значения $M$ находятся в области $\left[10^{3}, 10^{4}\right]$. Кроме того, $\hat{A}$ – очень плоская функция, со многими экстремумами. Число итераций для достижения достаточно хорошего приближения находится также в области $\left[10^{3}, 10^{4}\right]$, что замедляет процесс.

Рассмотренный метод ясно обнаруживает прохождения, близкие к седловой точке $\hat{A}$. Он может быть использован в будущем при попытке локализовать эти решения. Еще один возможный подход – это численная реализация леммы горного прохода, следуя дуге между двумя локальными минимумами градиентного потока.

Вариационный метод не предоставляет непосредственной информации об устойчивости периодических решений.
7.2. Уточнение решений

Имея некоторую начальную точку, определенную или на основе вариационного метода, или же как решение для близкого значения показателя $a$, мы можем попытаться использовать метод Ньютона. Рассмотрим значения
$z_{j}, j=1, \ldots, N$. Временной поток $2 \pi / N$ перемещает тело $j$ к телу $j-1$, где индексы следует рассматривать в $Z / N Z$. Таким образом, задачу можно преобразовать в поиск решения уравнения $G(Z):=\Phi(2 \pi / N, Z)-S(Z)=0$, где $Z$ обозначает все $z_{j}$ и $\dot{z}_{j}, \Phi$ – поток (1), а $S$ – смещение индексов.

Для решения $G(Z)=0$ методом Ньютона необходимо найти $D_{Z} G(Z)$. Его получают с помощью одновременного интегрирования обыкновенных дифференциальных уравнений (1) и уравнений первых вариаций. Этот процесс может оказаться утомительным, если $N$ велико. Например, $N=100$ приводит по существу к системе размерностью 160000 . Использование параллельных вычислений позволяет значительно уменьшить время вычисления, так как уравнения первых вариаций можно объединить в группы.

Одна из проблем состоит в том, что многие решения совершенно неустойчивы при больших $N$ и страдают от прохождений вблизи столкновений. Сходимость метода может быть гакже очень мала из-за существования очень большого количества близких решений. Эту проблему можно решить, ведя поиск значения $Z$ не при времени $t=0$, а также при некоторых промежуточных значениях $0=t_{0}<t_{1}<\ldots<t_{k-1}<t_{k}=2 \pi / N$. Пусть $Z_{m}$ – значение при $t_{m}$, тогда вместо одного уравнения $G=0$ получим систему уравнений $\Phi\left(t_{m+1} ; t_{m}, Z_{m}\right)-Z_{m+1}=0$. Это хорошо известный метод параллельной пристрелки (см. например, [10]). Система увеличилась в размере, но уравнения уже не требуют дополнительных усилий, и лучше наложены условия.

В качестве последнего комментария приведем следующие замечания: есть некоторая свобода в выборе исходного времени, а система также инварианта относительно вращения. Как здесь показано, изложенный метод недостаточен, потому что любое решение имеет связанное с ним шестимерное семейство решений, поэтому $D_{Z} G(Z)$ окажется сингулярным. Это затруднение легко преодолеть, выбрав начальное условие так, чтобы при $t=0$ некоторые из тел находились на оси $x$ с $\dot{x}=0$, а также используя приведение центра масс.

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

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

Выполнение этой работы стало возможным благодаря информации, предоставленной мне А. Шенсине, Р.Монтгомери и Дж. Джервером, а также многочисленным обсуждениям нами этой темы. Я также в долгу перед М. Хеноном, Р. Мартинесом и Р. Ортего за многие полезные замечания. Большая часть работы была выполнена во время академического отпуска, проведенного мной в Институте небесной механики (Парижская обсерватория), благодаря поддержке CNRS. Выражаю всем своим коллегам в этом институте свою признательность за гостеприимство и интерес к этой работе. Широко использовались также параллельные вычислительные возможности Группы динамических систем UB. Я благодарен также за поддержку со стороны грантов DGICYT PB 94-0215 (Испания) и CIRIT 1998SGR-00042.