1.1. Постановка задачи

В этой статье изучается уравнение Хилла с квазипериодической вынуждающей силой
\[
\ddot{x}+\left(a^{2}+b p(t)\right) x=0,
\]
${ }^{*}$ Broer H., Simó C. Hill’s equation with quasi-periodic forcing: resonance tongues, instability pockets and global phenomena, Bul. Soc. Bras. Mat. 29, 1998, 253-293. Перевод Богатыревой Е. В., Килина А. А.

где $a$ и $b$ – параметры. Параметрическая вынуждающая сила $p=p(t)$ является квазипериодической функцией $p(t)=P\left(t \omega_{1}, t \omega_{2}\right)$ с рационально независимыми частотами $\omega_{1}, \omega_{2}$ и дважды периодической вещественной аналитической функцией $P: T^{2} \rightarrow \mathbb{R}$. Точнее говоря, мы главным образом рассматриваем частный случай
\[
p(t)=\cos \left(\omega_{1} t\right)+\cos \left(\omega_{2} t\right) \text {, или } P\left(\theta_{1}, \theta_{2}\right)=\cos \left(\theta_{1}\right)+\cos \left(\theta_{2}\right) \text {. }
\]

Выберем $\omega_{1}=1$ и $\omega_{2}=\gamma$, где $\gamma=\frac{1}{2}(1+\sqrt{5})$ – золотое сечение. Кроме того, будут рассмотрены некоторые возмущения (2). При конечных значениях параметра $a$ исследуются два случая. Случай малых $|b|$, когда некоторые резонансные явления изучаются с помощью методов усреднения. И случай больших $|b|$, для которого приводятся численные результаты и их эвристические объяснения.

Уравнение (1) можно также записать как $-\partial^{2} u / \partial x^{2}-b V(x) u=a^{2} u$, просто заменив $x$ на $u$, а $t$ на $x$. Это классическое уравнение Шредингера (см, например, [10]), где потенциал $V$ заменяет квазипериодическую функцию $p(t)$, а $a^{2}$ играет роль энергии.
Уравнение второго порядка (1) можно переписать в виде системы
\[
\left\{\begin{array}{l}
\dot{\theta}=\omega, \\
\dot{x}=y, \\
\dot{y}=-\left(a^{2}+b P(\theta)\right) x,
\end{array}\right.
\]

где $(\theta,(x, y)) \in \mathbb{T}^{2} \times \mathbb{R}^{2}$ и $\omega=\left(\omega_{1}, \omega_{2}\right)$. Заметим, что двумерный тор $\mathbb{T}^{2} \times\{(0,0)\}$, заданный равенством $x=y=0$, инвариантен относительно квазипериодического потока с частотным вектором $\omega$. Первая цель исследования заключается в изучении нормального линейного поведения этого инвариантного тора, в частности, его устойчивости и его приводимости к форме Флоке. Вторая цель – исследовать случай большого $b$, чтобы увидеть, какого рода интересные явления происходят при его увеличении.

Примечание. В случае необходимости мы можем расширить систему (3) до гамильтоновой системы с фазовым пространством $\mathbb{T}^{2} \times \mathbb{R}^{2} \times \mathbb{R}^{2}=$ $=\{\theta, J,(x, y)\}$, где $J$ канонически сопряжен $\theta$. Несомненно, если мы возьмем в качестве гамильтониана $H(x, y, \theta, J)=\frac{1}{2} y^{2}+\omega J+\frac{1}{2}\left(a^{2}+b P(\theta)\right) x^{2}$, то система принимает вид
\[
\dot{\theta}=\omega, \quad J=-\frac{b}{2} \frac{\partial P}{\partial \theta} x^{2}, \quad \dot{x}=y, \quad \dot{y}=-\left(a^{2}+b P(\theta)\right) x .
\]

1.2. Задачи и результаты
1.2.1. Резонансные полуострова

Сначала рассмотрим несколько существенных моментов классической теории уравнений Хилла и Матье с периодической вынуждающей силой. Уравнение Хилла (1) с периодической по времени вынуждающей силой $p$ широко изучалось в литературе, см., например, [5,6], а также указанные там источники. Соответствующие работы с нелинейной параметрической вынуждающей силой см., например, в [2,7]. В этом случае на плоскости $(a, b)$ образуется бесконечно много резонансных полуострова, начинающихся в резонансных точках $(a, b)=\left(\frac{k}{2}, 0\right), k \in N$, в которых эти полуострова имеют более или менее острые «вершины». С внешней стороны каждого полуострова тривиальное $2 \pi$-периодическое решение $T^{1} \times\{(0,0)\}$, заданное равенством $x=y=0$, устойчиво. Классический периодический случай уравнения Матье получается из (1) при $p(t)=\cos (t)$. В этом случае $k$-й резонансный полуостров имеет верхушку со степенью заостренности порядка $(k-1)$. Заметим, что в периодическом случае, множество резонансов дискретно. На рис. 1 показаны резонансные полуострова для классического периодического случая.

В рассматриваемом квазипериодитеском случае (1) с $p$, задашшы урав нением (2), можно ожидать чего-нибудь подобного в отношении инвариантного двумерного тора $\mathbb{T}^{2} \times\{0\}$ векторного поля (3). В этом случае множество резонансов
\[
(a, b)=\left(\frac{1}{2}\left(k_{1}+k_{2} \gamma\right), 0\right),\left(k_{1}, k_{2}\right) \in \mathbb{Z}^{2}
\]

плотно заполняет ось $a$. Можно поинтересоваться, взаимодействуют ли «близкие» резонансные полуострова, при малых значениях $|b|$, или же они отделены вплоть до средних значений $b$. На рис. 2 показаны резонансные полуострова в случае, когда $\rho$ задано выражением (2) с соотношением частот, равным золотому сечению.
1.2.2. Очаги неустойчивости

Еще одно интересное явление касается очагов неустойчивости. В периодическом случае они встречаются в возмущениях классических случаев Матье типа $p(t)=\cos (t)+c \cos (2 t)$, где $c$ – малый параметр. В действительности границы второго резонансного полуострова, возникающего при $(a, b)=(1,0)$, пересекаются на расстоянии порядка $|c|$ от оси $a$, вследствие

Рис. 1. Полуострова классического уравнения Матье на плоскости $\left(a^{2}, b\right)$

чего образуется очаг неустойчивости. В [5,6] дано качественное объяснение этого явления на основе теории особенностей. Здесь мы рассмотрим аналог такого возмущения с помощью следующей модификации уравнения (2):
\[
p(t)=\cos (t)+\cos (\gamma t)+c \cos \left(\left(k_{1}+k_{2} \gamma\right) t\right)
\]

при фиксированных целых чисел $k_{1}, k_{2}$ и малом параметре $c$.
Приведем некоторый численный аргумент в пользу появления очагов неустойчивости. Резонансный полуостров с верхушкой в $(a, b)=\left(\frac{1}{2}\left(k_{1}+\right.\right.$ $\left.\left.+k_{2} \gamma\right), 0\right)$ имеет две границы, которые можно достаточно точно приблизить с помощью процедуры усреднения или нормализации, при которой величина $a-\frac{1}{2}\left(k_{1}+k_{2} \gamma\right)$ раскладывается в степенной ряд по $b$. В нормальной форме каждая граница устойчивости определяется двумя типами вкладов, один порожден членом $b(\cos (t)+\cos (\gamma t))$, а другой $-b c \cos \left(\left(k_{1}+k_{2} \gamma\right) t\right)$. Кроме того, члены, являющиеся комбинациями обоих возмущений, могут

Рис. 2. Основные полуострова уравнения Хилла с квазипериодической вынуждающей силой $p(t)=\cos t+\cos (\gamma t)$ на плоскости $\left(a^{2}, b\right)$

играть важную роль. Все эти типы вклада могут иметь главные члены, зависящие от $c$ и различных степеней $b$. Следовательно, они «сокращаются» при уменьшении значений $b$, создавая, таким образом, очаги неустойчивости. Этот процесс подобен тому, что случается в периодическом случае, см. [5,6]. На рис. 3 показан очаг неустойчивости для $p$, заданного уравнением (4), где $k_{1}=1, k_{2}=1$, а $c=0.3$. Позже мы вернемся к анализу этих полуостровов, основанному на нормальных формах.
1.2.3. Дальнейшие задачи, численное исследование

Линейное уравнение (1) задает нормализованную линейную часть инвариантного двумерного тора $\mathbb{T}^{2} \times\{0\}$. Будем говорить, что (1) или, что то же самое, (3), приводимо (к форме Флоке), если с помощью линейной замены переменных $x$ и $y$, квазипериодически зависящей от времени, его можно привести к уравнениям с постоянными коэффициентами. В периодическом

Рис. 3. Очаг неустойчивости уравнения Хилла с квазипериодической вынуждающей силой $p(t)=\cos t+\cos (\gamma t)+0.3 \cos (1+\gamma) t)$ на плоскости $\left(a^{2}, b\right)$

случае эту задачу разрешает классическая теория Флоке. В квазипериодическом случае приводимость не всегда возможна. В $[22,10,11],[15],[17,18]$, $[19]$, а также $[3,4,8,9]$ можно найти несколько результатов, относящихся к неприводимости в подобном или более широком смысле.

При больших значениях $|b|$ метод возмущений менее эффективен, и мы прибегаем к численному исследованию. Максимальный показатель Ляпунова и число вращения являются хорошими инструментами изучения уравнения (1) или же системы (3) при любом значении параметров $(a, b)$. С целью такого изучения мы вводим алгоритмический подход, который оказывается очень полезным на практике.

Сначала мы переформулируем задачу следующим образом. В периодическом случае вся информация относительно свойств (1) содержится в матрице монодромии, т.е. фундаментальной матрице, вычисленной при $t=T$ ( $T$ является периодом $p$ ), и которая при $t=0$ равна единичной матрице. Также можно использовать матрицу Флоке, т. е. логарифм матрицы монодромии, деленный на $T$. Сравните, например, [5,6].

В рассматриваемом случае (2) двум углам $\theta_{1}, \theta_{2}$ соответствуют частоты $\omega_{1}=1$ и $\omega_{2}=\gamma$, зафиксируем один угол (например $\theta_{1}$ с периодом $2 \pi$ ). Тогда угол $\theta_{2}$ будет принимать последовательные значения $\{2 n \pi \gamma \bmod 2 \pi\}, n \in N$. Поэтому матрица «монодромии» также будет зависеть от значения $\theta_{2}$ в начале каждого временного интервала. Пусть $M(\theta)-$ фундаментальная матрица (1) с $p(t)=\cos (t)+\cos (\theta+\gamma t)$, начинающаяся с единичной матрицы при $t=0$ и вычисленная при $t=2 \pi$. Далее мы определяем последовательность матриц $\left\{P_{n}\right\}$, используя преобразование косого произведения:
\[
(P, \theta) \mapsto(M(\theta) P, \theta+2 \pi \gamma \bmod 2 \pi)
\]

начиная с $(P, \theta)=(I, 0)$. Таким образом, $P_{n}(\theta)=M(\theta) P_{n-1}(\theta), n \in N$, и все матрицы $P_{n}(\theta)$ симплектические. Отметим, что в периодическом случае для линейного отображения Пуанкаре $P$ выполняется $P_{n}=P^{n}$.

Нам интересны «средние» свойства этой последовательности. Мы изучаем асимптотическое поведение (его существование) как свойства «расширения», так и «вращения». Пусть $\mu_{n}$ – главное собственное значение $P_{n}$. Тогда максимальный показатель Ляпунова определяется выражением
\[
\lambda=\lim _{n \rightarrow \infty} \frac{1}{n} \log \left|\mu_{n}\right|,
\]

которое оказывается независимым от исходного значения $\theta$. Аналогичным образом мы рассматриваем вращение, соответствующее $P_{n}$. При этом следует рассматривать вращение, не учитывая $\bmod 2 \pi$. Для проведения эффективных вычислений достаточно начать с любого нормированного вектора (например, $(1,0)$ ) и затем рассмотреть аргументы последовательности векторов, полученных из начального последовательным применением матрицы $M(2 n \pi \gamma \bmod 2 \pi)$. Пусть $\alpha_{n}$ – полученные значения аргументов (для их вычисления необходимо знать свойства последовательного действия матриц $M$ на вектора). Тогда число вращения можно определить как
\[
\rho=\lim _{n \rightarrow \infty} \frac{1}{n} \alpha_{n},
\]

и оно также не будет зависеть от исходного значения $\theta$.
Примечания.
1) Мы приведем некоторые полезные советы по эффективным вычислениям, просто описав методы, использованные при приближенном исследовании. Конкретные детали и дальнейшие уточнения см. в приложении.
Один из возможных методов состоит в интегрировании (1) на длительном промежутке времени (скажем, вплоть до $t=2 \pi \times 10^{5}$ или $t=2 \pi \times 10^{6}$ ). Однако этот процесс занимает много времени и, кроме того, может привести к сильному распространению ошибок. Поэтому для него требуются большая эффективность и точность.
2) Поскольку необходимо знать множество матриц $M(\theta)$, то для уменьшения вычислений используется интерполяционный метод. При этом вычисляются $M(\theta)$ для равноудаленных значений $\theta$ в интервале $[0,2 \pi)$, а для всех остальных значений $M(\theta)$ определяется с помощью интерполяции. Обычно число точек берется в пределах от 200 до 1000, а порядок полиномиальной интерполяции от 7 до 15.
3) Вычисление $M(\theta)$ для каждого из выбранных значений $\theta$ выполняется точно и быстро с помощью метода рядов Тейлора. Этот метод требует вычисления производных высокого порядка. Зная $x$ и $\dot{x}$, значение $\ddot{x}$ и производные более высокого порядка можно получить используя уравнение (1) и правило Лейбница. Для эффективности вычислений необходимо определить оптимальные порядок разложения и длину шага. Типичный оптимальный порядок колеблется от 20 до 30.
4) Действуя описанным выше способом, можно легко получить до $10^{6}$ итераций. Для проверки точности найденных оценок $\lambda$ и $\rho$ применяется следующий метод. Через каждые 10000 итераций оценивается значение $\lambda$ и $\rho$. Полученная оценка считается хорошей, если три последовательных значения (например, после 130000, 140000 и 150000 итераций) равны друг другу с заданной точностью. Типичная величина точности колеблется в пределах от $10^{-5}$ от $10^{-7}$. Более того, вычисление $\lambda$ и $\rho$ в зависимости от параметров проводится так, чтобы от шага к шагу их значения мало изменялись (это необходимо, чтобы не пропустить интересные явления). Во время проведения эксперимента шаг по $a$ при необходимости уменьшался до $10^{-9}$.
Первая теоретическая задача заключается в определении существования этих двух чисел. Ясно, что если (1) приводимо, то соответствующая матрица Флоке и приводящее преобразсвание непосредственно обеспечивают нам их существование. Для малого $b$ множество значений $a$ с неустойчивым двумерным тором имеет число вращения из множества
\[
\mathcal{R}=\left\{\left.\frac{1}{2}\left(k_{1}+k_{2} \gamma\right) \in \mathbb{R} \right\rvert\, k_{1}, k_{2} \in \mathbb{Z}\right\}
\]

и по определению точно соответствует резонансному полуострову. Теперь зададим малое значение $b$ и рассмотрим изменение $a$. Несмотря на тот факт, что множество пересечений полуостровов плотно, мера его дополнения все еще относительно большая, что доказано, например, в [17]. Более конкретно, множество значений $a$, для которого приводимость невозможна, экспоненциально мало по $b$, при общих условиях. В действительности, результат даже шире, поскольку в [17] рассмотрена только приводимость к устойчивому или «эллиптическому» случаю, и, следовательно, резонансные зоны отброшены. Однако хорошо известно, что в гиперболическом случае приводимость всегда возможна. Эти результаты также подтверждаются в рамках КАМ-теории.

При малом $b$ большая часть значений $a$ должна соответствовать точкам за пределами резонансных полуостровов. Действительно, расчетные значения $\lambda$ незначительно отличаются от нуля, что можно увидеть на рис. 4 слева, где для $b=0.2$ построен график показателя Ляпунова в зависимости от $a^{2}$ в интервале $[0.25,1]$. На рисунке справа показаны соответствующие знатения числа вращения в том же интервале. Существенную роль, как видно из рисунков, играют только резонансы $\left(k_{1}, k_{2}\right)=(1,0)$ и $(0,1)$ (максимумы $\lambda$ слева и в центре рисунка). Очень малый рост $\lambda$, связанный с резонансами $(-2,2)$ и $(3,-1)$ обнаруживается для $a^{2}$, близких к 0.382 и 0.477 соответственно. Ясно видно, что показатель Ляпунова положителен только в резонансных полуостровах, когда значение $a$ принадлежит $\mathcal{R}$. Кроме «больших» ступенек, видимых на графике числа вращения, существует счетное множество шагов едва видимых ступенек (соответствующих каждому значению $\rho \in \mathbb{Q}$ ). На дополнении $\mathcal{R}$ значение $\lambda$ равно нулю.

Заметим также, что из результатов работы [16] следует, что в некотором интервале по $a$ возможна некоторая «почти» приводимость, равномерно по a. То есть уравнение (1) можно привести к уравнениям с постоянными коэффициентами, с экспоненциально малым по $b$ остатком.

С другой стороны, для больших значений $b$ (в нашем случае это начинает происходить при $b$ около 0.5 ) и для $a \in[0.5,1]$ мы также все еще можем говорить о «резонансных полуостровах». Здесь, как и выше, число вращения принадлежит $\mathcal{R}$, но утверждение о том, что точки вне этих полуостровов имеют $\lambda=0$ уже не является истинным. На рис. 5 показаны те

Рис. 4. Графики показателя Ляпунова и числа вращения для $b=0.2$. Слева: максимальный показатель Ляпунова $\lambda$. Справа: число вращения $\rho$.

же данные, что и на рис. 4 , но для $b=0.6$. Здесь видно намного больше ступеней, и к тому же, за их пределами показатель Ляпунова отличается от нуля. На рис. 6 и 7 изображены увеличения малых областей изменения $a$ для близких значений $b$, а дополнительные сведения будут приведены в разделе 5.

Очевидно, что для значений $(a, b)$ за пределами полуостровов с $\lambda>0$ приведение к форме Флоке невозможно, поскольку в матрицу $2 \times 2$ невозможно включить $\pm \lambda, \pm \rho$.

Интересное глобальное явление касается линии $\mathcal{C}$ в плоскости $(a, b)$, см. рис. 9 и 10. Здесь, по-видимому, происходит «коллапс резонансов». Более того, при фиксированном $b$ слева от этой кривой точки пространства параметров, лежащие вне полуостровов, по-видимому, имеют положительное значение $\lambda$. Численные исследования происхождения коллапса см. в разделе 5. Краткое изложение см. также в выводах.