Журнал вычислительной математики и математической физики, 2020, T. 60, № 12, стр. 2055-2072

1. ВВЕДЕНИЕ

Уравнения маятникового типа, которыми мы ограничимся в данной работе, понимаются как дифференциальные уравнения второго порядка с периодическими коэффициентами. Наиболее типичный (но не наиболее общий) пример – это уравнение Хилла. Математический маятник (с различными модификациями), уравнение Ван дер Поля, уравнение Дюфинга, уравнение Белецкого и т.д. – привести полный список было бы невозможно.

Такого рода уравнения либо являются содержательными задачами сами по себе, либо появляются в результате редукции более сложных физических задач.

Большинство из упомянутых задач считаются решенными, но продолжают использоваться уже в качестве модельных задач для демонстрации возможностей новых (или хорошо забытых старых) методов.

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

Например, в работах [1], [2] уравнение Белецкого исследовалось методом усреднения при малых значениях инерциального параметра $\mu $. Однако в то время как само уравнение Белецкого зависит от параметра $\mu $ аналитическим образом, для метода усреднения потребовалось вводить малый параметр $\sqrt {{\text{|}}\mu {\kern 1pt} {\text{|}}} $, что нарушает аналитичность уравнения и повышает коразмерность вырождения. Иными словами, полученные ряды по малому параметру состоят наполовину из нулей.

Другой подход к этой же задаче предложен в [3], [4], где никаких замен параметров не требуется.

Пример с уравнением Белецкого показывает, что вполне понятное желание как-то упростить систему и свести ее к уже изученному стандарту может иметь обратный эффект даже для ручного счета.

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

Мы используем эти соображения для конструирования некоторого универсального алгоритма, позволяющего вычислять периодические решения уравнений маятникового типа в виде рядов по малому параметру.

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

Таким образом, согласно регулярной теории возмущений Пуанкаре, решение задачи обычно ищут в виде ряда

где $h$ – это выбранный малый параметр, а ${{y}_{0}}(t)$ – это решение невозмущенной задачи, которое считается известным. Подстановка решения (1) в исходное уравнение дает треугольную систему ОДУ, которая решается последовательно в виде квадратур. При этом также возникают условия разрешимости (коэффициенты при резонансных членах), которые должны выполняться, чтобы полученное решение было периодическим.

Все эти манипуляции относительно легко осуществимы в полуавтоматическом (т.е. интерактивном) режиме, когда требуется небольшое количество членов разложения (1).

Однако этот подход заведомо неоптимален, когда планируется получить достаточно большое число членов разложения (1).

Во-первых, исходное уравнение (вообще говоря) нелинейно. Поэтому для того, чтобы вычислить решение с точностью $O({{h}^{n}})$, может потребоваться подстановка (1) с некоторым запасом. В любом случае подстановка (1) сразу же дает весьма громоздкое выражение.

Во-вторых, решение полученной треугольной системы ОДУ требует манипуляций с тригонометрическими полиномами и анализа квадратур на предмет выявления секулярных членов.

Мы предлагаем подход, который, по сути, реализует метод малого параметра, схематично изложенный выше, но не требует ни разложения уравнений в ряды, ни манипуляций с тригонометрическими полиномами.

Рассмотрим разложение

которое формально используется вместо (1) при подстановке в исходное уравнение.

В случае если ищется малое решение (а невозмущенное решение ${{y}_{0}}(t)$ не обязано быть малым), правую часть (2) необходимо умножить на $h$, т.е. предварительно раздуть решение.

Однако никакой фактической подстановки (2) в исходное уравнение не требуется, если использовать формальные уравнения в вариациях (см. [3]).

Исходное уравнение зависит от малого параметра $h$ явным образом (вообще говоря). Решение $y(t)$ также зависит от $h$, но неявно. Предлагается следующий формализм для вычисления треугольной системы ОДУ для формальных вариаций решения.

Обозначим решение $y(t) = {{y}_{0}}(t)$ как вариацию нулевого порядка и подставим в исходное уравнение. Затем формально дифференцируем исходное уравнение по параметру $h$ по правилу

нужное число раз.

Здесь, разумеется, возникает задача: как научить CAS правильно дифференцировать произвольные выражения по переменной, которая не входит явным образом в некоторые функции. Мы вернемся к этой задаче позже.

После формального дифференцирования (нужное число раз) получим треугольную систему ОДУ для формальных вариаций решения при $h \ne 0$. Теперь подставляем $h = 0$ в эту систему и получаем треугольную систему ОДУ для вариаций на невозмущенном решении ${{y}_{0}}(t)$, т.е. по сути уравнения для функций, входящих в разложение (2).

Отметим некоторые общие моменты, которые появляются при решении любой из задач указанного типа.

Поскольку период решения (вообще говоря) зависит от малого параметра, то вычисление вариаций следует делать не в исходном уравнении, а в уравнении с предварительной заменой времени $t = P(h)s$ (или $t = s{\text{/}}\omega (h)$) так, что по времени $s$ все функции уже будут $2\pi $-периодическими.

Это единственное (и очевидное) преобразование исходного уравнения, которое применяется только в случаях переменного периода.

Таким образом, помимо вариаций ${{y}_{k}}(t)$, в уравнения будет входить неизвестная функция $P(h)$, которая обрабатывается точно таким же образом, как и ${{y}_{k}}(t)$ в (3).

Соответственно, наряду с разложением (2) получим также разложение (нормированного) периода (или частоты)

где величины ${{P}_{k}}$ – это однозначно определяемые вещественные константы, причем ${{P}_{0}} = 1$, если невозмущенная задача имеет $2\pi $-периодические решения.

В результате формального дифференцирования (нужное число раз) и подстановки $h = 0$ получим треугольную систему ОДУ для вариаций ${{y}_{k}}(s)$ и величин ${{P}_{k}}$ на невозмущенном решении.

Заметим, что первая вариация в этой системе автоматически окажется линеаризацией исходного уравнения на невозмущенном решении.

Таким образом, общая ситуация выглядит так.

Нулевая вариация на невозмущенном решении – это линейное (и часто однородное) дифференциальное уравнение второго порядка, имеющее $2\pi $-периодические решения.

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

Последовательное решение неоднородных ОДУ требует выполнения условий разрешимости, которые в механике обычно связывают с резонансами.

Иными словами, на каждом шаге решение выражается в виде квадратуры, но оно будет периодической функцией только при условии обнуления секулярных (т.е. полиномиальных по времени) членов. Коэффициенты при резонансных членах, которые породили секулярные, дают алгебраические уравнения, необходимые для определения коэффициентов решения.

При этом необходимо учесть, что выполнение условий разрешимости на шаге $n$ позволяет определить коэффициенты решения (вообще говоря) только на шаге $n,n - 1,n - 2, \ldots $, в зависимости от вырожденности рассматриваемой задачи.

Иными словами, задача сводится к решению бесконечной линейной системы алгебраических уравнений для $M$-диагональной матрицы, где $M = 1,2, \ldots $ .

В невырожденном случае все нужные на шаге $n$ величины определяются на том же шаге $n$, т.е. этот случай соответствует диагональной линейной системе.

Как мы отмечали, работа с квадратурами в системах CAS – не самый оптимальный способ решать треугольные системы ОДУ данного типа. Поэтому предлагается следующий подход.

Пусть уравнение для нулевой вариации на невозмущенном решении имеет вид

где $L$ – это линейный дифференциальный оператор второго порядка (хотя данный подход обобщается и на уравнения высших порядков). Тогда уравнение для $n$-й вариации имеет вид где ${{R}_{n}}(s)$ – это $2\pi $-периодическая функция, зависящая от предыдущих вариаций.

Функции ${{R}_{n}}(s)$ будем записывать в виде полиномов Лорана от переменной $u = exp(is)$.

Таким образом, все уравнения $L(y) = {{u}^{k}}$, $k \in \mathbb{Z}$ будут однозначно разрешимы в виде полинома Лорана $y(u)$, если $k$ не является резонансным значением.

Какие $k$ являются резонансными – это определяется один раз для каждой задачи. Обычно $k = \pm 1$.

После этого остается только разложить ${{R}_{n}}(u)$ на мономы по $u$ и разделить их на резонансные и нерезонансные.

Нерезонансные члены дадут компоненты решения ${{y}_{n}}(s)$ после применения оператора ${{L}^{{ - 1}}}$, а коэффициенты при резонансных членах дадут уравнения для определения неизвестных коэффициентов решения на шаге $n,n - 1,n - 2, \ldots $, в зависимости от вырожденности рассматриваемой задачи.

Ну и, разумеется, на каждом шаге к полученному частному решению прибавляется решение однородного уравнения с неизвестными коэффициентами, которые и определятся из резонансных членов на последующих шагах.

Далее на ряде примеров мы покажем, как предложенный алгоритм работает в конкретных задачах.

2. УРАВНЕНИЕ ВАН ДЕР ПОЛЯ

Уравнение Ван дер Поля

где $h$ – это малый параметр, является одной из наиболее популярных модельных задач маятникового типа.

Ясно, что при $h = 0$ уравнение (4) имеет двупараметрическое семейство $2\pi $-периодических решений.

Можно легко показать, что уравнение (4) является диссипативным при $h > 0$ и имеет предельный цикл, период которого зависит от $h$.

Однако мы не будем пользоваться известной информацией об уравнении Ван дер Поля либо использовать какие-либо преобразования, кроме описанных в разд. 1.

Сделаем в уравнении (4) замену времени $t = P(h)s$ так, что по времени $s$ функция $x(s)$ будет $2\pi $-периодической:

Очевидно, $P(0) = 1$, поэтому дифференциальный оператор $L$, т.е. левая часть всех уравнений в вариациях, имеет вид

Обозначим решение $x(s)$ как вариацию нулевого порядка, т.е. $x(s) = {{y}_{0}}(h)$, где зависимость решения от $s$ нас пока что не интересует, а неявную зависимость решения от $h$ мы сделали явной. Аналогично, $P(h) = {{P}_{0}}(h)$. Тогда получим выражение, которое требуется формально дифференцировать по $h$

где ${{z}_{0}}(h)$ обозначает производную ${{y}_{0}}(h)$ по $s$.

Обратим внимание, что функции, входящие в выражение (5), явно зависят от $h$, поэтому алгоритмы символьного дифференцирования в системах CAS их распознают как функции от $h$. Однако вместо обычного символьного дифференцирования применяется правило (3). Например, первая производная (5) имеет вид

Ясно, что формулы быстро становятся необозримыми. Поэтому, как мы отмечали, данные вычисления не предназначены для ручного счета. В то время как для CAS такие вычисления не представляют трудностей.

Нетрудно догадаться, что, вычислив несколько вариаций и подставив в них $h = 0$, мы автоматически получим уравнения для коэффициентов тейлоровских разложений всех нужных нам функций. Это делается следующим образом.

Нулевая вариация на невозмущенном решении – это просто линейный осциллятор, поэтому

где ${{a}_{0}}$ и ${{b}_{0}}$ – это некоторые пока неизвестные вещественные коэффициенты. Однако исходное уравнение автономно, т.е. инвариантно относительно сдвига по времени, поэтому сразу можно взять ${{b}_{0}} = 0$ (и далее аналогично).

Поскольку нам нужны также производные функций по независимой переменной $s$, в то время как вычисления ведутся с полиномами Лорана от $u = exp(is)$, то производная по $s$ запишется как

по правилу дифференцирования сложной функции.

Таким образом, получим первую вариацию

Нетрудно проверить, что уравнение

имеет решение т.е. для нерезонансных значений $k$ интегрирование ОДУ сводится к линейной операции, в то время как резонансные значения $k = \pm 1$ дают секулярные члены, которые требуется обнулить.

Поэтому для первой вариации (и всех последующих) надо разложить правую часть в сумму мономов по $u$ и обработать их соответственно их степени.

Поскольку задача вещественна, то два комплексных коэффициента при $k = 1$ и $k = - 1$ комплексно сопряжены. Они дают два уравнения относительно вещественных неизвестных ${{a}_{0}}$ и ${{P}_{1}}$.

На самом деле мы получаем уравнение разветвления, так как полученная система имеет решения

Первое из этих решений даст просто тождественный нуль, а два других эквивалентны, поэтому выберем ${{a}_{0}} = 1$.

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

Определившись с нужной ветвью и отбросив уже использованные резонансные члены, получаем общее решение для первой вариации

и ее производной ${{z}_{1}}$ по правилу (6).

Далее все вычисления однотипны, но достаточно громоздки, чтобы их приводить в явном виде. Они легко программируются в любой системе CAS, так как требуются только простые операции с полиномами.

Из второй вариации получаем резонансный член при степени $u$

откуда однозначно определяем ${{a}_{1}} = 0$ и ${{P}_{2}} = 1/8$, и т.д.

Таким образом, получим ${{a}_{{2n - 1}}} = 0$, $n \in \mathbb{N}$, и

А также, ${{P}_{{2n - 1}}} = 0$, $n \in \mathbb{N}$, и

Напомним, что мы использовали не обычные производящие функции, т.е. степенные ряды, а экспоненциальные производящие функции. Поэтому, например, нормированный период

Хотя полученный ряд для $P$ является четной функцией $h$, сэкономить на этом не получится, так как исходное уравнение зависит аналитически именно от $h$. Возвращаясь к исходным переменным, получим

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

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

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

Что касается уравнения Ван дер Поля, то никаких закономерностей в коэффициентах решения, данных выше, не наблюдается. Что неудивительно, так как римановы поверхности решений этого уравнения устроены фрактально сложно (см. [5]).

Мы предлагаем читателю сравнить наш подход с традиционным методом усреднения, применявшимся к этой задаче в [6, с. 63].

3. УРАВНЕНИЕ ДЮФФИНГА

Уравнение Дюффинга

также весьма популярно как модельная задача для демонстрации работы алгоритмов вычисления периодических решений (см., например, [6, с. 71], [7]).

Кроме того, это уравнение может рассматриваться как упрощенный вариант уравнения математического маятника, которое можно обработать точно таким же образом, как показано ниже.

Нас интересуют малые решения уравнения (7) и как они эволюционируют, когда становятся не малы.

Однако уравнение (7) консервативно (и обратимо), поэтому все достаточно малые решения будут периодическими.

Это означает, что в качестве малого параметра следует взять некоторый диаметр самого решения, т.е. первый интеграл. Положим

Раздуем решение так, что нулевое приближение будет уже не мало, т.е. сделаем замену $x(t) = Hw(t)$ в уравнении (7) и его первом интеграле (8). Получим уравнение

и его первый интеграл что являет собой регулярно возмущенный линейный осциллятор.

Теперь видно, что первоначально введенный малый параметр неудобен, так как уравнения (9) и (10) аналитически зависят от параметра $h = {{H}^{2}}$. Его мы и будем далее использовать.

Сделаем в уравнениях (9) и (10) замену времени $t = P(h)s$ так, что по времени $s$ функция $w(s)$ будет $2\pi $-периодической. Получим уравнения

и

Теперь также видно, что можно немного сэкономить, взяв замену времени $t = \sqrt {P(h)} s$, и получить более простые формулы для вариаций. Но мы не будем этого делать.

Поскольку $P(0) = 1$, то дифференциальный оператор $L$, т.е. левая часть всех уравнений в вариациях, такой же, как и для уравнения Ван дер Поля. Обозначим решение $w(s)$ как вариацию нулевого порядка, т.е. $w(s) = {{y}_{0}}(h)$, $dw(s){\text{/}}ds = {{z}_{0}}(h)$, как мы делали ранее. Аналогично, $P(h) = {{P}_{0}}(h)$. Тогда получим выражения, которые требуется формально дифференцировать по $h$

Нулевая вариация на невозмущенном решении – это просто линейный осциллятор, поэтому

где мы положим ${{b}_{0}} = 0$, пользуясь автономностью исходного уравнения (и далее аналогично).

Нулевая вариация первого интеграла $I(h)$ дает

где по правилу дифференцирования (6).

Подставив ${{y}_{0}}$ и ${{z}_{0}}$ в нулевую вариацию $I(h)$, получим уравнение разветвления

т.е. ${{a}_{0}} = \pm 1{\text{/}}2$. Поскольку эти два решения эквивалентны в силу обратимости системы, то выберем ${{a}_{0}} = 1{\text{/}}2$.

В отличие от уравнения Ван дер Поля, здесь все нужные величины определяются на каждом шаге, поэтому далее просто выполняется цикл до нужного порядка. А именно, на шаге $n$:

${{P}_{n}}$ определяется из резонансного члена; далее

${{y}_{n}} = {{a}_{n}}u + {{a}_{n}}{\text{/}}u + {{L}^{{ - 1}}}$ (от нерезонансных членов);

${{z}_{n}} = d{{y}_{n}}{\text{/}}ds$ по правилу (6);

${{a}_{n}}$ находится из $n$-й вариации $I(h)$.

Таким образом, получим для $n \in {{\mathbb{N}}_{0}}$

и

Напомним, что мы использовали экспоненциальные производящие функции. Поэтому (нормированный) период решения для исходного времени $t$

Возвращаясь к исходным переменным, получим решение $x(t)$, аналогичное данному в разд. 2 для уравнения Ван дер Поля.

Однако здесь нас интересует вопрос, есть ли какие-либо закономерности в коэффициентах решения. Мы не нашли на него ответа в литературе, поэтому применили программу “Guess”, реализованную на языке Maple, и которая находится в открытом доступе.

Можно также применить встроенный пакет “gfun”, который когда-то находился в открытой библиотеке Maple-программ “Share”.

В результате получим

где ${\mathbf{K}}$ – это эллиптический интеграл I рода.

Из этой формулы непосредственно видно, что радиус сходимости степенного ряда для $P$ равен $h = 1{\text{/}}2$, что, впрочем, можно получить прямо из анализа первого интеграла, так так при $h = 1{\text{/}}2$ фазовый портрет имеет сепаратрису.

Обратим также внимание, что при $h > 0$ аргумент эллиптического интеграла является чисто мнимым, так что $P \to \infty $ при $h \to 1{\text{/}}2$ и имеет очень сложную степенно-логарифмическую асимптотику.

Отметим также, что в уравнении математического маятника

которое выглядит существенно более “нелинейным”, чем уравнение Дюффинга, период малых колебаний устроен значительно более просто, и имеет вид хотя наша программа в CAS не смогла распознать этот ряд.

4. СИСТЕМА ХЕНОНА–ХЕЙЛЕСА

Система Хенона–Хейлеса (СХХ) (см. [8]) интенсивно изучалась на протяжении десятилетий, и поэтому весьма удобна как модельная задача для демонстрации возможностей различных алгоритмов и методов. Например, в работе [9] (и см. ссылки там) вычисляются малые периодические решения этой задачи с помощью приведения к нормальной форме.

Приведение систем такого рода к нормальной форме является весьма нетривиальной задачей само по себе и требует владения (в нескольких смыслах) специальным программным обеспечением.

Нормальная форма удобна тем, что в ней видны некоторые скрытые симметрии задачи, которые используются для разбиения задачи на случаи и подслучаи, которые изучаются отдельно.

Это особенно относится к гамильтоновой нормальной форме СХХ, где используются переменные действие–угол (см. [10]).

Также определенные свойства нормальной формы используются для доказательства сходимости полученных рядов (см. [9]).

Однако платой за эти удобства (помимо обладания нужными навыками и владения пакетом программ) является значительное усложнение задачи.

Так, первым же шагом нормализации системы является ее комплексификация, т.е. удвоение размерности системы. Уже это обстоятельство может привести к существенному усложнению решений. Например, в уравнении Ван дер Поля самое сложное поведение решений – это предельный цикл. В то время как комплексификация этого уравнения демонстрирует хаотическое поведение решений (см. [5]).

В СХХ комплексификация приводит впоследствии к необходимости отлавливать комплексные решения, которых нет в исходной задаче (см. [9]).

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

Причем здесь требуется как сходимость рядов по малому параметру, так и сходимость нормализующего преобразования.

Заметим, что сходимость самих рядов по малому параметру часто довольно очевидна, так как они получаются в рамках регулярной теории возмущений. А вот сходимость нормализующего преобразования проблематична и требует дополнительного анализа.

Поскольку сходимость рядов по малому параметру влечет аналитичность решения, которую мы интенсивно эксплуатируем, то остановимся на этом вопросе подробнее.

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

Поэтому докажем следующее

Предложение. Пусть имеется аналитическая система ОДУ, разрешенная относительно производной и аналитически зависящая от малого параметра $h \in \mathbb{C}$ при достаточно малых ${\text{|}}x{\kern 1pt} {\text{|}}$ и ${\text{|}}t{\kern 1pt} {\text{|}}$:

Тогда решение $x(t,h)$, $x(0) = {{x}_{0}}$, аналитически зависит от $h$ при достаточно малых ${\text{|}}{{x}_{0}}{\kern 1pt} {\text{|}}$, ${\text{|}}t{\kern 1pt} {\text{|}}$ и ${\text{|}}h{\kern 1pt} {\text{|}}$.

Доказательство. Согласно теореме Коши решение $x(t,h)$ является голоморфной функцией времени $t$, т.е. разлагается в абсолютно сходящийся степенной ряд с коэффициентами, которые являются голоморфными фукциями $h$ по условию. Поэтому можно переставить порядок суммирования и представить $x(t,h)$ как голоморфную функцию $h$ с коэффициентами, которые являются голоморфными фукциями $t$. Q.E.D.

Разумеется, аналитическая продолжаемость решения по (вещественному) времени $t$ является свойством конкретной задачи. Но в случае существования малых периодических решений, т.е. продолжаемости решения при $h = 0$ на конечный интервал, возмущенные малые решения также продолжаются на некоторый интервал, аналитически зависящий от $h$, что следует из леммы Бореля о конечном покрытии и теоремы о монодромии (аналитической продолжаемости).

Совершенно аналогичные утверждения доказаны в несколько более частном случае, но зато гораздо подробнее в [11, с. 50, с. 55].

Возвращаясь к основной теме данной статьи, мы применим наш минималистский подход к СХХ и покажем, что все малые периодические решения находятся непосредственно из самих уравнений и без каких-либо преобразований, кроме описанных во введении.

СХХ является гамильтоновой системой четвертого порядка и, на первый взгляд, не относится к уравнениям маятникового типа, изучаемым в этой статье. Однако это не так.

Запишем эту задачу в комплексной форме:

где ${{y}_{1}}(t)$ и ${{y}_{2}}(t)$ – это обычные вещественные координаты в классическом гамильтониане (см. [9]) т.е.

Тогда первый интеграл системы (он же гамильтониан) принимает вид

где малый параметр $h$ вводится обычным образом, а параметр $H$ введен для сравнения с работой [9].

В такой записи СХХ становится весьма похожей на уравнение Дюффинга.

Заметим, что эта компактная форма записи системы является “бонусом”, но не является необходимой для нашего исследования.

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

Отметим также, что инвариантные подпространства, найденные в СХХ в [9] с помощью нормальной формы, легко находятся непосредственно из уравнений (12) и (13).

Так, взяв $(a)$: ${{x}_{1}}(t) = {{y}_{1}}(t) \equiv 0$ или $(b)$: ${{x}_{1}}(t) \equiv \pm \sqrt 3 {{x}_{2}}(t)$, ${{y}_{1}}(t) \equiv \pm \sqrt 3 {{y}_{2}}(t)$ в (12), получим вещественное уравнение

Однако мы не будем пользоваться этой информацией и найдем все семейства решений непосредственно из уравнений (12) и (13).

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

Как было отмечено, вычисления ведутся в интерактивном режиме вплоть до уравнения разветвления, после чего вычисления организуются в виде цикла.

В зависимости от вырожденности задачи, т.е. от некоторых внутренних “резонансов”, связанных отчасти с соотношениями между собственными значениями линейной части системы, уравнения разветвления могут появляться на разных шагах алгоритма.

Причем это свойство вырожденности, будучи по сути асимптотическим, совершенно не зависит от приведения системы к нормальной форме или других регулярных (или формальных) преобразований.

Так, из результатов работы [9] следует (хотя там это явно не отмечено), что в СХХ имеются два уравнения разветвления. Одно на 2-м, а другое на 4-м шаге алгоритма. Только в [9] эти шаги делаются для нормальной формы СХХ, а у нас эти вырождения проявятся непосредственно в СХХ.

Следуя описанному алгоритму, раздуем решение так, что нулевое приближение будет уже не мало, т.е. сделаем замену $z(t) = hw(t)$ в уравнении (12) и его первом интеграле (13).

Затем в полученных уравнениях сделаем замену времени $t = \sqrt {P(h)} s$ так, что по времени $s$ функция $w(s)$ будет $2\pi $-периодической. Получим уравнение

и его первый интеграл

Дадим некоторые пояснения по поводу проделанных замен переменных.

Замена времени выбрана в такой форме для (незначительного) упрощения последующих вычислений.

Отказ от использования функции ${\text{Re}}(\,)$ связан с тем, что мы работаем с полиномами Лорана от переменной $u = exp(is)$. Это, в свою очередь, влечет отказ от встроенных в CAS функций, работающих с комплексными величинами. Вместо этого мы используем одну линейную функцию ${\text{C}}(\,)$ (комплексное сопряжение) со следующими свойствами:

Однако при этом мы по-прежнему пользуемся встроенными в CAS средствами работы с комплексными переменными.

Таким образом, обозначив $w(s) = {{y}_{0}}(h)$, $dw(s){\text{/}}ds = {{z}_{0}}(h)$, $P(h) = {{P}_{0}}(h)$, как мы делали ранее, получим выражения, которые требуется формально дифференцировать по $h$

где $L$ – это тот же оператор (т.е. линейный осциллятор), что и для всех предыдущих случаев.

Нулевая вариация на невозмущенном решении – это линейный (комплексный) осциллятор, а последующие вариации будут уже задаваться неоднородными уравнениями, поэтому

где все величины ${{a}_{n}},{{b}_{n}},{{c}_{n}},{{d}_{n}}$, $n \in {{\mathbb{N}}_{0}}$, являются вещественными.

На самом деле, в силу автономности системы часть величин можно сразу же положить равными нулю, что следует из следующей леммы.

Лемма 1. Пусть периодическое решение автономного ОДУ раскладывается в абсолютно сходящийся ряд по малому параметру $h$ в виде

где внутренняя сумма по $k$ на самом деле конечна. Тогда сдвигом по времени $s$ можно добиться, чтобы все величины ${{A}_{{n,1}}}$ или ${{B}_{{n,1}}}$ стали вещественными (или чисто мнимыми).

Доказательство. Поменяем порядок суммирования в решении, что допустимо по условию. Тогда

Сдвиг по времени $s \to s + {{s}_{0}}(h)$ умножит сумму

на величину $exp(i{{s}_{0}}(h))$. Поэтому эту сумму можно сделать вещественной. Но это аналитическая функция вещественной переменной $h$, поэтому все коэффициенты ${{A}_{{n,1}}}$ следует выбрать вещественными. Q.E.D.

Эта лемма имеет непосредственное отношение к структуре множества периодических решений автономной системы ОДУ. Однопараметрическое по $h$ семейство таких решений лежит на двумерном множестве в фазовом пространстве в силу возможности сдвига вдоль траектории. Поэтому при продолжении семейства периодических решений по $h$ всегда существует некоторый произвол, который реализуется в виде появления произвольных постоянных, от которых решение на самом деле не зависит.

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

Нулевая вариация интеграла $I(h)$ после подстановки ${{y}_{0}}$ и ${{z}_{0}}$ дает

и пока что это выражение никак не может быть использовано.

Здесь и далее следует воздержаться от назначения коэффициентам некоторых значений для решения получаемых уравнений.

Условия разрешимости для первой вариации $L(h)$ (коэффициенты при $u$ и $1{\text{/}}u$) дают

поэтому необходимо ${{P}_{1}} = 0$ (напомним, что ${{P}_{0}} = 1$).

Первая вариация интеграла $I(h)$ дает уравнение

и это уравнение также остается пока в запасе.

Первое уравнение разветвления возникает на втором шаге (который соответствует членам при ${{h}^{2}}$).

Коэффициенты при резонансных членах (которые зависят только от ${{a}_{0}}$, ${{b}_{0}}$, ${{c}_{0}}$, ${{d}_{0}}$ и ${{P}_{2}}$) и их комплексные сопряжения, плюс уравнение (18) дают пять уравнений от пяти неизвестных. Эту систему проще всего решить, вычислив базис Гребнера. В результате получаем, что либо ${{P}_{2}} = 5{\text{/}}3$ (случай 1), либо ${{P}_{2}} = - 2{\text{/}}3$ (случай 2).

После выбора ${{P}_{2}}$ используем вычисленный базис Гребнера повторно, но уже без переменной ${{P}_{2}}$. В результате получаем

Таким образом, ветвление решений не является (вообще говоря) бинарным, хотя случаи 2а и 2б оказываются эквивалентными в силу обратимости СХХ.

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

Поскольку эта задача уже подробно исследована, то мы не станем подробно описывать все различные ветви периодических решений (многие из которых эквивалентны в силу симметрий СХХ и ее обратимости), а ограничимся только вычислением семейств с различными периодами.

Разберем сперва случай 2(а, б) как более простой.

Вычисления необходимо проводить в интерактивном режиме до шага 4 (т.е. членов при ${{h}^{4}}$ в разложении решения). Этот факт априори ниоткуда не следует, но только начиная с пятого шага становится ясной структура решения, и вычисления программируются в виде цикла.

Пусть ${{c}_{0}} = {{d}_{0}} = 0$ (случай 2а). Тогда оказывается, что ${{c}_{k}} = {{d}_{k}} = 0$, $k \in \mathbb{N}$. Кроме того, условие $a_{0}^{2} + b_{0}^{2} = 1{\text{/}}2$ никак далее не упрощается, что соответствует цилиндричности множества в фазовом пространстве, соответствующего семейству периодических решений в гамильтоновой системе с двумя степенями свободы (см. [13, теорема 6.5.1.]). Иными словами, коэффициенты ${{a}_{0}}$ и ${{b}_{0}}$ можно выбрать произвольно, но так, чтобы выполнилось условие $a_{0}^{2} + b_{0}^{2} = 1{\text{/}}2$. А также коэффициенты ${{a}_{n}}$ или ${{b}_{n}}$ далее можно выбирать произвольно, т.е. ${{a}_{n}} = {{a}_{n}}({{b}_{n}})$ или наоборот, что определяется из $n$-й вариации первого интеграла (17).

Поэтому можно выбрать ${{b}_{n}} = 0$, $n \in {{\mathbb{N}}_{0}}$. Тогда, вернувшись в начало алгоритма, начиная с шага $n = 3$ вычисления программируются в виде цикла. Получаем ${{P}_{n}}$ из резонансного члена при $u$ или ${{u}^{{ - 1}}}$, в то время как ${{a}_{n}}$ находятся из $n$-й вариации первого интеграла (17). Таким образом, получим ${{a}_{{2n + 1}}} = {{P}_{{2n + 1}}} \equiv 0$, $n \in {{\mathbb{N}}_{0}}$, и

Вычислим разложение (нормированного) периода

а затем подставим в это разложение $h = \sqrt {2E} $ или $h = \sqrt {2H} $ для соответствия малого параметра работам [12] или [9] соответственно. Тогда получим коэффициенты разложения в табл. 3 в [12], т.е. коэффициенты так называемого “circular” периода. Мы посчитали эти коэффициенты до степени ${{h}^{{40}}}$ (или ${{H}^{{20}}}$).

Если теперь разложить $1{\text{/}}\tilde {P}(h) = \omega (H)$ в ряд по $H$, то получим случай ${{\omega }_{3}}$ (см. [9, формула (17)]).

К сожалению, никаких закономерностей у “circular” семейства для коэффициентов периода или амплитуды не обнаружилось.

Рассмотрим теперь случай 1.

Здесь встретилась ситуация, не описанная в литературе, в том числе и в литературе по СХХ. А именно, в СХХ имеется второе уравнение разветвления, которое появляется на 4-м шаге алгоритма. В работах [10], [12], [9], где были найдены все семейства малых периодических решений, этот вопрос обойден вниманием, а различные семейства описываются просто как случаи и подслучаи.

В случае 1 уместно применить лемму 1 и сразу положить ${{d}_{n}} \equiv 0$, $n \in {{\mathbb{N}}_{0}}$. Тогда ${{c}_{0}} = \pm 1{\text{/}}2$. Поскольку оба эти случая эквивалентны (что априори ниоткуда не следует и требует проверки, которую мы опускаем), то положим ${{c}_{0}} = - 1{\text{/}}2$.

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

На шаге $n = 1$ и $n = 3$ получим ${{P}_{1}} = 0$ и ${{P}_{3}} = 0$. На шаге $n = 4$ из базиса Гребнера получаем ${{P}_{4}} = 33{\text{/}}4$ или ${{P}_{4}} = 145{\text{/}}4$.

Поскольку нас интересует, главным образом, период, выберем наиболее просто устроенные семейства для каждого из периодов.

Один из случаев – это “curvilinear” период в табл. 2 в [12]. В наших обозначениях получаем

Здесь, как и для “circular” периода, коэффициенты ${{b}_{n}}$, $n \geqslant 1$, можно выбирать произвольно, что соответствует различным путям (различной параметризации) на цилиндричной поверхности в фазовом пространстве, соответствующей семейству. Этот эффект не был ранее описан в литературе по СХХ.

Выберем ${{b}_{n}} \equiv 0$, $n > 1$, что соответствует некоторому сечению семейства решений. Тогда вычисления проводятся в виде цикла, начиная с шага $n = 5$ следующим образом.

На шаге $n$, ${{P}_{n}}$ и ${{a}_{{n - 2}}}$ находятся из системы 2 уравнений, которые дают коэффициенты при резонансных членах. Далее, ${{c}_{n}}$ находятся из вариации первого интеграла. Таким образом, получим ${{b}_{n}} = {{d}_{n}} = {{P}_{{2n + 1}}} \equiv 0$, ${{a}_{{2n}}} = - {{c}_{{2n}}}$, $n \in {{\mathbb{N}}_{0}}$, и

Проделав преобразования для периода, как в (19), и разложив $1{\text{/}}\tilde {P}(h) = \omega (H)$ в ряд по $H$, получим случай ${{\omega }_{{(1.4)}}}$ [9, формула (17)].

Продолжим рассмотрение cлучая 1, но возьмем “rectilinear” семейство (см. [12]). Так же как и ранее, ограничимся рассмотрением самого простого пути на цилиндричной поверхности, соответствующей этому семейству решений.

В наших обозначениях “rectilinear” семейство начинается как

и вычисляется далее в виде цикла так же, как и предыдущее “curvilinear” семейство.

Таким образом, получим ${{b}_{n}} = {{d}_{n}} = {{P}_{{2n + 1}}} \equiv 0$, ${{a}_{{2n}}} = {{c}_{{2n}}}$, $n \in {{\mathbb{N}}_{0}}$, и

Проделав преобразования для периода, как в (19), и разложив $1{\text{/}}\tilde {P}(h) = \omega (H)$ в ряд по $H$, получим случай ${{\omega }_{{(2.5)}}}$ [9, формула (17)]. В данном случае для периода существует выражение в замкнутой форме

где $F(\,)$ – это классическая гипергеометрическая функция (см. также [12]).

Это единственный случай в СХХ, когда период семейства найден в явном виде.

Доказательство этого факта можно получить из приведения СХХ к виду (14). Оба случая в (14) имеют один и тот же период (20), так как значение гамильтониана для них вычисляется по-разному.

Заметим, что наиболее подробное исследование семейств периодических решений ССХ было проделано в [12], где использовался метод Линштедта–Пуанкаре. Наш метод, по сути, можно рассматривать как модернизацию и упрощение этого классического метода, адаптированную для вычислений в CAS.

5. СИСТЕМЫ, ЗАВИСЯЩИЕ ОТ ВРЕМЕНИ

Если в уравнения явно входят периодические по времени функции, то периодическое решение имеет смысл искать в виде ряда Фурье, т.е. как функцию того же периода. Подстановка такого ряда с неизвестными коэффициентами в исходное уравнение и поиск этих коэффициентов носит название “метода гармонического баланса” и широко используется в инженерных задачах.

Однако и в этом случае наш метод вычисления периодических решений имеет явные вычислительные преимущества. Мы проиллюстрируем этот тезис на некоторых известных примерах.

Рассмотрим уравнение

Спрашивается, при каких вещественных $a$ и $b$ существуют малые по $h$ периодические решения $z(t)$ уравнения (21)?

Эта задача при $a = b = 1$ рассматривалась в [14] как пример гамильтоновой системы, не имеющей малых периодических решений.

Уравнение

имеет решение поэтому ясно, что в данном случае резонансным значением будет $k = 0$.

Действуя по нашей схеме, на нулевом шаге получим $z(t) = {{a}_{0}}$, где ${{a}_{0}}$ пока не определено.

На первом шаге коэффициент при резонансном члене дает уравнение разветвления

откуда находим, что если $a < 0$, то $z(t) \equiv \sqrt { - 1{\text{/}}a} $ будет единственным периодическим решением уравнения (21) при $b \ne - 2$.

Если же $b = - 2$, то на шаге 2 найдем, что либо ${{a}_{0}} = \sqrt { - 1{\text{/}}a} $ (и т.д., тривиальное решение), либо ${{a}_{0}} = 0$. В этом случае все произвольные постоянные далее определяются с запаздыванием на 2 шага как нули. Таким образом, получим

Рассмотрим еще один пример из статьи [14]: частный случай хорошо известной задачи – маятник с вибрирующей точкой подвеса при отсутствии силы тяжести

Как пример применения теоремы 3 в [14] автор нашел малые по $\mu $ периодические решения этой задачи, которые ответвляются от положений равновесия $z = 0$ и $z = \pi $, и исследовал их устойчивость.

Однако наши вычисления показывают, что таких решений не существует. Это следует также из формулы (5.5) в [14], где приведены вращательные периодические решения уравнения (22), т.е. $z(t) = nt + \theta (t)$, где $\theta (t)$ – $2\pi $-периодическая функция. Формула (5.5) в [14] формально применима при числе вращения $n = 0$, т.е. для просто периодических решений, и в этом случае дает $z \equiv 0$. Заметим также, что при $n = 1$ решением будет функция

Отметим, что линеаризованное в нуле уравнение (22) – это уравнение Матье (после замены времени $t = 2\tau $) с параметрами $a = 0$ и $q = - 2\mu $, которое не имеет периодических решений, так как для их существования необходимо ненулевое собственное значение $a \ne 0$ (см. [15, с. 951]).

Действуя по нашей схеме, на нулевом шаге получим $z(t) = {{a}_{0}}$, где ${{a}_{0}}$ пока не определено. Наши вычисления дают уравнение разветвления на втором шаге

Нетрудно проверить, что ${{a}_{0}} = 0$ или ${{a}_{0}} = \pi $ дадут тривиальные решения. Другая ветвь ${{a}_{0}} = \pi {\text{/}}2$ дает решение Уместно сравнить наш подход с тем, как эта задача решалась в [16].

Данная задача для вращательных периодических решений оказалась значительно более нетривиальной. Следуя [14], положим

и будем искать $2\pi $-периодические решения $\theta (t)$, малые вместе с $\mu $.

Как мы показали, при $n = 0$ единственным периодическим решением (23) является тождественный нуль. При $n > 1$, действуя по нашей схеме, получим

что совпадает с формулой (5.5) в [14], за исключением произвольной константы интегрирования ${{a}_{1}}$.

Далее формулы получаются весьма громоздкими, и мы их опускаем.

Мы не можем априори положить константы интегрирования равными нулю, так как они могут понадобиться на следующих шагах для устранения резонансных членов. Однако оказалось, что все эти константы появляются в виде множителя при резонансных членах, т.е. необходимо обнуляются, но на более поздних шагах.

Если принять этот экспериментальный факт как теорему, то это означает, что уравнение (23) имеет единственное и притом нечетное $2\pi $-периодическое решение $\theta (t)$ при любом $n \in \mathbb{N}$.

Далее оказалось, что четные и нечетные числа вращения $n$ дают совершенно разное поведение решений.

При $n$ нечетном все константы интегрирования оказываются необходимо равными нулю, но с некоторым запаздыванием. А именно, константу ${{a}_{k}}$ необходимо положить равной нулю только на $(n + k)$-м шаге алгоритма (т.е. при вычислении члена разложения при ${{\mu }^{{n + k}}}$). Таким образом, задача тем более вырождена, чем больше (нечетное) число вращения.

При $n$ четном все константы интегрирования необходимо положить равными нулю, но со значительно большим запаздыванием. А именно, константа ${{a}_{k}}$ должна стать нулем только на $(2n + 2 + k)$-м шаге алгоритма (т.е. при вычислении члена разложения при ${{\mu }^{{2n + 2 + k}}}$). Таким образом, при четном $n$ задача оказывается значительно более вырожденной, чем при нечетном числе вращения.

Приведем для ясности таблицу, дающую шаг алгоритма $N$, где константу ${{a}_{1}}$ необходимо положить равной нулю в зависимости от числа вращения $n$:

Столь необычно большая вырожденность задачи (далее мы уточним этот термин) приводит к появлению одного эффекта, не описанного в литературе (насколько нам известно).

Возьмем число вращения $n = 4$ и посчитаем разложение решения $\theta (t)$ до членов 10-го порядка включительно (т.е. с точностью до членов порядка $o({{\mu }^{{10}}})$). Это намного превышает обычную точность вычисления подобных разложений.

Полученная функция $\theta (t)$ зависит от 10 произвольных констант, которые входят в четную часть этой функции. Положим ${{a}_{k}} = k$ и вычислим начальное приближение функции $\theta (t)$ при $\mu = 0.1$ и $t = 0$. Получим (все знаки верны)

Численным интегрированием (с большой точностью) убеждаемся, что ошибка за период $2\pi $ не превосходит $4 \times {{10}^{{ - 12}}}$, что также превосходит обычную точность численных расчетов в обычной арифметике.

Итак, получено решение, которое неотличимо от периодического как численно, так и в символьном виде (до 10-го порядка) и не является нечетной функцией.

Однако резонансный член на 11-м шаге даст секулярную добавку к решению в виде

что объясняет, почему эту погрешность почти невозможно отследить численно.

Таким образом, при большом четном числе вращения существует множество решений задачи (23), которые практически неотличимы от периодических, но таковыми не являются. Причем численно этот эффект также практически не наблюдаем.

6. ЗАКЛЮЧЕНИЕ

Нами предложен (как мы полагаем) универсальный метод вычисления периодических решений в системах ОДУ с малым параметром. Причем для обоснования существования разложений и их сходимости используются только две фундаментальные теоремы: теорема Коши об аналитической зависимости решения от начальных данных и теорема Пуанкаре об аналитической зависимости решения от параметра.

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

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

Решение всех ОДУ сводится к линейным операциям с полиномами после того, как разрешилось уравнение разветвления (обычно нелинейное), т.е. после того, как ветвь решений выбрана однозначно.

При этом, как мы видели, уравнение разветвления может быть неединственным и появляться на разных шагах алгоритма.

Чем позже появляется уравнение разветвления, тем более сложным является решение задачи с алгоритмической точки зрения. Поэтому уместно вырожденность задачи измерять номером шага, на котором появляется уравнение разветвления. Примем это в качестве определения меры вырожденности задачи. Ясно, что данное определение вырожденности не сводится к резонансам для собственных значений линейной части системы, а является существенно нелинейным.

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

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

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

Данное определение вырожденности задачи является по сути асимптотическим и потому не зависит от регулярных или даже формальных замен переменных (например, приведение системы к нормальной форме).

Как показано в разд. 5, простая механическая система может обладать богатым набором “периодических” решений, которые, на самом деле, таковыми не являются. В природе такие явления наблюдаются довольно часто. Например, две нейтронные звезды, вращающиеся вокруг их общего центра масс, постепенно сближаются, а Луна, обращаясь вокруг Земли, постепенно от нее удаляется.

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