В этом приложении мы описываем методы, используемые в численных экспериментах. Основная их цель – достичь достаточно хорошей точности за малое время вычислений. Несомненно для получения надежных численных результатов требуется точность. Возьмем за единицу времени вычислений время выполнения одной итерации (5) или время оценки $\lambda$ и $\rho$ вне начального переходного процесса. Тогда общее время вычислений проведенных исследований близко к $10^{13}$. В статье приведена только очень небольшая выборка полученных результатов. Для получения же всех результатов требуется высокая быстрота вычислений. Большая часть рассмотренной задачи была выполнена на кластере Беовульф, состоящем из 24 процессоров Pentium 266 и 2 процессоров Pentium 200. Каждый процессор способен выполнять примерно $2 \times 10^{5}$ итераций в секунду.
7.1. Вычисление $M(\theta)$

Нам необходимо проинтегрировать уравнение вида $\ddot{x}=q(t) x$. В рассматриваемом случае $q(t)=-\left(a^{2}+b(\cos (t)+\cos (\gamma t+\theta))\right.$, где $\gamma=\frac{1}{2}(1+\sqrt{5})$ и для большинства вычислений $a$ и $b$ были ограничены интервалом $[0,1]$.

Введем $x^{(k)}=\frac{1}{k !} \frac{d^{k} x}{d t^{k}}$. Аналогичным образом введем $q^{(k)}$. Тогда из правила Лейбница получаем соотношение
\[
x^{(k+2)}=\frac{1}{(k+1)(k+2)} \sum_{j=0}^{k} x^{(k-j)} q^{(k)},
\]

с помощью которого начиная с $x^{(0)}$ и $y^{(0)}=x^{(1)}$ можно рекуррентно вычислить все $x^{(k)}$ вплоть до любого порядка $N$ с вычислительными затратами $O\left(N^{2}\right)$. Затраты на тригонометрические вычисления заключаются только в нахождении $\cos$ и $\sin$ на каждом шагу, при этом степени $\gamma$ и обратная величина факториалов вычисляются только один раз.

Ряды для $x(t+h)$, зависящие от $h$ и начинающиеся с $x(t), y(t)$, мажорируются рядом функции $z(h)$, задаваемым уравнением $\ddot{z}(h)=r(h) z(h)$, с $r(h)=a^{2}+|b|(\exp (h)+\exp (\gamma h)), z(0)=|x(t)|$ и $\dot{z}(0)=|y(t)|$. Это показывает, что ряд сходится. Кроме того, если мы ограничимся указанным выше интервалом изменения $a$ и $b$ и возьмем $h<1$, тогда $r(h)<9$. Следовательно, $z(h)<A \exp (3 h)$, где $A$ зависит от $z(0)$ и $\dot{z}(0)$. Это показывает, что при выборе, скажем, $N=30$, последний член более чем в $10^{18}$ раз меньше первых двух.

Фактически мы интегрируем уравнение, начиная с векторов $(1,0)$ и $(0,1)$. При расчетах мы выбирали постоянный временной шаг $(2 \pi / 8$ или $2 \pi / 10$ обычно дают оптимальное время вычислений), а $N$ выбираем так, чтобы максимум двух последних членов, появляющихся в разложениях 4 элементов фундаментальной матрицы был более, чем в $10^{18}$ раз меньше максимума абсолютных значений элементов матрицы.

Вычисления были выполнены для большого числа $K$ значений $\theta$, равно распределенных в интервале $[0,2 \pi]$. При точных вычислениях, когда нас интересовало большое количество итераций, подходящее значение $K$ равнялось 5000. Большое время задания начальных условий (обычно 3 секунды) позволяет сократить время вычислений при итерациях.
7.2. Процесс интерполяции

Нам известны значения $M\left(\theta_{k}\right)$, где $\theta_{k}=2 \pi k / K, k=0, \ldots, K-1$. По этим значениям необходимо вычислить $M(\theta)$. Это делается с помощью интерполяции. Масштабируя и смещая начало координат, мы всегда можем предположить, что нам следует оценить функцию $f(s)$ в некоторой точке
$s \in[0,1)$ на основании значений $f_{j}=f(j), j=1-j_{0}, \ldots, j_{0}$. Интерполяционная формула имеет вид $\sum_{j=1-j_{0}}^{j_{0}} w_{j}(s) f_{j}$, где $w_{j}(s)$ – лагранжевы веса. Для еще большего ускорения процесса вычислений были предварительно вычислены веса $w_{j, i}$ для $s=i / 10000, i=0, \ldots, 10000$. Затем на основании предварительно вычисленных весов с помощью линейной интерполяции были получены соответствующие веса $w_{j}(s)$. Проверено, что это не приводит к значительным различиям в вычислениях $M(\theta)$. Интерполяцию следует применять к каждому из элементов матрицы. Подходящее значение для $j_{0}$ равно трем. Поэтому для интерполяции используется только 6 узлов. Численные оценки шестой производной элементов $M(\theta)$ в предварительных интегрированиях показывают, что ошибки в интерполяционной формуле имеют порядок $\epsilon$.
7.3. Оценки $\lambda$ и $\rho$

Формулы (6) и (7) определяют $\lambda$ и $\rho$ при условии существования пределов. Однако вследствие «неправильного» поведения итераций, чтобы получить достаточно хорошие оценки, возможно, число итераций $n$ должно быть взято очень большим. На рис. 25 приведен соответствующий пример, показывающий поведение частного из (6) после большого переходного процесса, который будет объяснен ниже. Типичное значение количества итераций в переходном процессе, в области, близкой к линии коллапса, равно $10^{6}$. Максимальное количество итераций было взято $10^{7}$. Результаты, приведенные в п. 5 , и поведение наблюдаемых на многих графиках оценок доказали пригодность такого эвристический подход для оценки $\lambda$. Этот подход также был использован и для $\rho$, несмотря на тот факт, что поведение частного из (7) более регулярно.

После переходного процесса, который в некоторых случаях достигает $2 \times 10^{6}$ итераций, частные в (6) вычисляются после каждой итерации, но записывается только самое большое частное в блоках, например, из 1000 последовательных итераций. Записывается также номер итерации, в которой имеет место максимум. Тогда после вычисления ряда блоков (скажем, 100 блоков; в последующих оценках используются 200,300 и т.д.) выполняется подбор значений $\alpha$ и $\beta$ для выражения $\lambda_{n_{i}}=\alpha+\beta / n_{i}$. Здесь $n_{i}-$ номера итераций, в которых встречаются максимумы, а $\lambda_{n_{i}}$ – соответствующие частные. Значение $\alpha$ равно улучшенной оценке показателя Ляпунова.

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

Рис. 25. Оценка $\lambda$ в зависимости от количества итераций для указанных значений $(a, b)$. Во всех случаях приведены максимумы и минимумы в блоках по $10^{3}$ итераций. В а) мы приводим только максимумы, при этом минимумы находятся ниже $2.1 \times 10^{-16}$. Отметим, что в случае с) для параметров выше линии коллапса и за пределами резонанса поведение наиболее нерегулярное, несмотря на переходный процесс из $2 \times 10^{6}$ итераций и вычисление дополнительных $7 \times 10^{6}$ итераций. Однако представляется, что даже в этом случае lim inf и limsup оценок $\lambda$ совпадают

бор значений $\alpha$ и $\beta$, отбрасываем все точки ниже полученной кривой, с оставшимися точками снова производим подбор, опять отбрасываем точки и т. д. до тех пор, пока число используемых точек в последнем подборе не окажется меньше некоторого заданного числа. На практике в качестве этого числа берется 50. Значение $\alpha$ последнего подбора используется как оценка $\lambda$. Наблюдалось, что полученные таким образом оценки становятся намного более согласованными при возрастании количества блоков. Когда три последовательные оценки отличаются меньше, чем на некоторое допустимое отклонение (для наиболее точных вычислений оно принималось равным $10^{-8}$ ), последнее из них используется как конечная оценка $\lambda$.

Для значений $(a, b)$, которым, очевидно, соответствует нулевое значение $\lambda$, результаты полученных описанным выше способом вычислений были всегда меньше $10^{-8}$, а в большинстве случаев даже меньше $10^{-10}$. Кроме того, в резонансных зонах поведение полученных значений гладкое с точностью до $10^{-8}$ (это можно определить, например, изучая конечные разности для равнораспределенных значений $a^{2}$ ). В областях за пределами резонансных полуостровов, выше линии коллапса, поведение также гладкое в пределах этого допустимого отклонения. Только вблизи пересечения линии коллапса, по-видимому, некоторые оценки могут иметь ошибки вплоть до $10^{-6}$.
7.4. Процесс сканирования

До сих пор мы описывали, как получить $\lambda$ и $\rho$ для заданных значений $(a, b)$. Задав значение $b$, мы определяли значения $\lambda$ и $\rho$ при изменении $a$ в заданном интервале. В этом «процессе сканирования» использовались следующие правила. Обозначим выбранные значения $a^{2}$, при которых выполнялись вычисления, как $a_{j}^{2}, a_{j}^{2}<a_{j+1}^{2}$. Выберем некоторые допустимые отклонения $\delta_{\lambda}$ и $\delta_{\rho}$ такие, что, $\left|\lambda\left(a_{j+1}^{2}\right)-\lambda\left(a_{j}^{2}\right)\right|<\delta_{\lambda}$ и $\mid \rho\left(a_{j+1}^{2}\right)-$ $-\rho\left(a_{j}^{2}\right) \mid<\delta_{\rho}$ при $a_{j+1}^{2}-a_{j}^{2}<\delta a_{\min }^{2}$. Типичные значения отклонений равны $\delta_{\lambda}=10^{-4}, \delta_{\rho}=10^{-5}, \delta a_{\min }^{2}=10^{-11}$. Эти значения были выбраны так, чтобы не пропустить большинство деталей графиков $\lambda$ и $\rho$ и в то же время эффективно осуществить сканирование.
7.5. Фурье-анализ

Как указывалось в п. 5, для выполнения фурье-анализа результатов в фазовом пространстве мы сначала должны исключить возможное экспоненциальное увеличение итераций. Это делается с помощью вычислений, после начального переходного процесса, матриц $\hat{P}_{n}=\exp (-n \lambda) P_{n}$, где $P_{n}$ – фундаментальные матрицы через время $2 \pi n$. (См. формулу (5).)

Мы рассматриваем только один из векторов матриц $P_{n}$. Точнее, мы использовали первый вектор. Если уравнение приводимо к уравнению с постоянными коэффициентами (форме Флоке), то фурье-анализ этого вектора должен содержать гармоники, содержащиеся в замене переменных (частоты которых обязательно должны быть целыми линейными комбинациями фундаментальных частот), а также гармоники частот, соответствующих приведенной матрице.

Обычно мы использовали первые $2^{18}$ итераций после переходного процесса, чтобы выполнить фурье-анализ. В этом случае с помощью алгоритма быстрого преобразования Фурье можно обнаружить главные гармоники. Некоторые полученные графики приведены в п. 5. Для определения основных частот существует несколько более точных методов (см., например, $[20],[21],[14])$. Однако можно использовать даже более простой метод, потому что известны основные частоты, и частота, соответствующая приведенному уравнению, приближенно равная $\rho$. Определение частот основных гармоник в случаях, которые выглядят приводимыми, выполняется с ошибкой меньше $10^{-7}$.

Благодарности. Авторы благодарят Ганса Яауслина за предложение изучить эту проблему им обоим, а также помощь в обмене идеями. Они также благодарны Расселу Джонсону, Анжелу Джорба и Йинфей Йи за полезное обсуждение проблемы. Исследования второго автора были поддержаны грантом DGICYT PB 94-0215 (Испания). Также выражается благодарность за частичную поддержку грантом Европейского Сообщества ERBCHRXCT940460 и грантом Каталонии CIRIT 1996S0GR-00105.