CC BY
40
Компьютерная алгебра
в научных вычислениях
УДК 519.61+517.925
Поиск семейств периодических решений обыкновенных дифференциальных уравнений с помощью метода нормальной формы. Часть 1
В. Ф. Еднерал*, О. Д. Тимофеевская^
* НИИ ядерной физики имени Д. В. Скобельцина
Московский государственный университет им. М. В. Ломоносова (НИИЯФ МГУ) Ленинские горы, д. 1, стр. 2, ГСП-1, Москва, Россия, 119991 Кафедра прикладной информатики и теории вероятностей Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198 ^ Кафедра физики высоких энергий и квантовой теории поля Физический факультет Московского государственного университета им.
М. В. Ломоносова Ленинские горы, д. 1, стр. 2, Москва, Россия, 119991
В настоящей работе кратко обсуждается применение метода резонансной нормальной формы к поиску семейств периодических решений автономных систем обыкновенных дифференциальных уравнений, разрешённых относительно производных и с полиномиальными нелинейностями в правых частях. При использовании сформулированного проф. А. Д. Брюно достаточном условии сходимости нормализующего преобразования, находятся локальные семейства периодических решений систем таких ОДУ в окрестностях стационарных точек. При этом в едином подходе исследуются как гамильтоновы, так и не гамильтоновы системы.
По соображениям объёма статья разбита на две части. В первой части описан алгоритм реализации метода нормальных форм. Отдельно кратко описаны созданные авторами программные пакеты. На языке RLISP разработан пакет для работы в системе REDUCE, а для работы с системой MATHEMATICA написан пакет на внешнем языке этой системы. Пакеты позволяют, в частности, получать формулы, описывающие локальные (содержащие неподвижную точку) семейства периодических решений. Результаты вычислений представляются в виде отрезков рядов Фурье заданной длины с частотой и коэффициентами, вычисленными в виде отрезков степенных рядов по параметру. Такое представление соответствует частному случаю отрезков рядов Пуассона. Важно, что при помощи единого алгоритма возможно изучать как двумерные, так и системы высоких порядков. Вторая часть статьи посвящена системам четвёртого порядка.
Сравнение табуляции полученных формул с численными решениями соответствующих уравнений показывает хорошее количественное согласие. Описываемый подход может быть использован при моделировании физических и биологических систем.
Ключевые слова: резонансная нормальная форма, динамические системы, локальные периодические семейства решений, компьютерная алгебра.
1. Введение
Метод нормальной формы основывается на преобразовании системы обыкновенных дифференциальных уравнений к более простой системе, называемой нормальной формой. На важность этого метода в случае исследования обыкновенных дифференциальных уравнений в окрестности неподвижной точки обратили внимание достаточно давно, например, [1] или [2].
Определение нормальной формы и нормализующего преобразования в литературе приводится в различных формах. Так, весьма развиты приложения для
Статья поступила в редакцию 26 июня 2014 г.
Первый автор поддержан грантом Президента РФ НШ-3042.2014.2.
гамильтоновых систем, например, [3-6], [7, (части 1,2)]. Резонансные нормальные формы и нормальные формы Белицкого исследуются в работах [8], [9, часть 5, §20], [10]. Существует немало алгоритмов (и их реализаций) для построения нормальных форм и соответствующих преобразований. Например для гамиль-тонова случая это улучшенный алгоритм Депри-Хори [11] и его компьютерно-алгебраическая реализация в системе REDUCE [12]. Сущность метода численного построения нормальных форм для Гамильтонианов описывается в работе [13]. Вопросы сходимости нормализующего преобразования обсуждаются в работах [5,6,14,15]. Что же касается построения нормальной формы в общем случае, мы упомянем здесь (в дополнение к книге А. Д. Брюно [16]) также статьи [17-19].
В настоящей работе мы будем использовать алгоритм, основанный на подходе, который был развит в работах А. Д. Брюно [5,6,8,9,16] для резонансной нормальной формы. Преимущество этого подхода состоит в возможности исследовать широкий класс автономных систем в рамках единой схемы, которая легко поддаётся алгоритмизации. В частности, этот подход обеспечивает конструктивный метод получения приближений для локальных семейств периодических и условно периодических решений в форме отрезков рядов Пуассона. Уделено особое внимание проблеме сходимости применяемых преобразований, что позволяет надеяться на то, что приближения для частот и семейств периодических решений около стационарных точек, полученные с помощью конечных формул, могут быть вычислены с нужной нам точностью. Кроме самих решений, мы можем также получать приближения для начальных условий, которые инициируют такие периодические решения. То есть, мы можем осуществлять элементы фазового анализа.
Другое достоинство используемого подхода состоит в алгоритмической простоте построения нормальной формы и соответствующего преобразования. Мы имеем прямую рекуррентную формулу для этой процедуры, вычисление которой не требует хранения большого числа промежуточных результатов, как этого требуют некоторые другие алгоритмы. Рассматриваемый подход свободен также от необходимости решать промежуточные системы уравнений и от любых ограничений на низко-резонансные случаи.
С помощью предлагаемого метода также возможно получать приближения для непериодических семейств. Результаты близки к результатам метода линеаризации Карлемана [20]. Для периодических и условно периодических случаев метод представляет собой обобщение метода Пуанкаре-Линдстеда [21] на многомерный случай.
Ниже мы опишем построение нормальных форм и их применение для построения приближений семейств периодических решений на примерах хорошо известных уравнений второго порядка. Во второй части статьи мы обсудим системы более высоких порядков.
2. Формулировка задачи
Рассмотрим систему автономных обыкновенных дифференциальных уравнений
Х = Ф(х), (1)
где x = (xi,... ,хп) — векторная функция времени, а X = dx/di — её полная производная по времени, Ф = (Ф1,..., Ф„) - вектор, который является функцией от x и, возможно, каких-либо параметров.
Данный тип уравнений возникает во множестве научных и инженерных задач, где имеют место осцилляции, вибрации или волновые процессы. Обычно при исследовании таких систем выполняются следующие действия:
1) бифуркационный анализ, то есть исследование картины поведения решений системы в зависимости от параметров. Особую важность имеет тот случай, когда это поведение резко изменяется при некоторых значениях параметров;
2) изучение фазового портрета, то есть исследование поведения решений системы
в зависимости от начальных условий;
3) поиск точных решений системы.
Если мы имеем решения системы в аналитической форме, мы имеем полную картину поведения системы, но получить такие решения удаётся крайне редко. Как правило, мы имеем лишь численные решения, а численные исследования пунктов 1 или 2 подчас являются очень сложной задачей. Например, совсем непросто получить численные решения в нестабильном случае. Поэтому есть интерес в некотором промежуточном подходе, который лежал бы между аналитическими и численными методами.
Метод нормальной формы широко используется для бифуркационного анализа. Об этих исследованиях можно прочитать, например, в работах [2,22-24]. Как продемонстрировано в этих книгах, численный бифуркационный анализ также основан на методе нормальной формы.
3. Предварительное преобразование системы
Изучение системы типа (1) в окрестности неподвижной точки х0, где Ф(х0) = 0, обычно включает в себя три предварительных шага. Во-первых, х сдвигается на —х0 так, что Ф(0) = 0, то есть, 0 — неподвижная точка. Каждая такая точка исследуется отдельно.
Второй шаг состоит в редукции системы к модельной форме, где вектор Ф(х) аппроксимируется вектором полиномов. Если в некоторой окрестности неподвижной точки Ф есть аналитическая функция х, тогда отрезок её степенного ряда может быть использован для получения аналитической аппроксимации. За счёт выполненного выше сдвига переменных правые стороны модельной системы будут состоять тогда из полиномов без свободных членов.
Третий шаг состоит в преобразовании матрицы линейной части системы к жордановой форме с помощью комплексного линейного преобразования х-пере-менных в у переменные.
После этих шагов система (1) принимает форму
Уг = А гУг + ^¿-1 + Ф ¿(у), г = 1,...,П, (2)
где Л = (А1,...,АП) — вектор собственных значений линейной части матрицы системы, аг — наддиагональные элементы жордановой формы и Ф = (Ф1,..., Фп) — вектор полиномов конечных степеней без констант и линейных членов.
В данной работе мы будем предполагать, что система (2) удовлетворяет следующим условиям:
— система автономна и имеет полиномиальные нелинейности;
— 0 — неподвижная точка и система будет изучаться в её окрестности;
— линейная часть правой части диагональна, то есть все а» = 0, причём не все собственные значения равны нулю, то есть Л = 0.
Это ограничение не является ограничением метода, оно сделано лишь для упрощения выкладок. В целом метод применим при любом виде жордановой матрицы [16], включая и вырожденный случай [25].
Не предполагается заранее, что система гамильтонова или сохраняет фазовый объем, или имеет какие-либо симметрии.
4. Метод нормальной формы 4.1. Нормализующее преобразование
При перечисленных выше ограничениях уравнения (2) могут быть записаны в форме
Уг = ^ + У1 ^ ЯчУ^ 1 = 1,...,n, (3)
где мы используем обозначение для мульти-индекса
п
Уч = П ^,
3 = 1
с векторной степенью разложения q == ..., дп), принадлежащей объединению множеств
Щ = {q е Ъп : щ > -1 и > 0, И] = г, ] = 1,...,п}.
Отрицательная степень может возникать в векторных показателях степеней за счёт выноса за скобки множителя yi из суммы в г-м уравнении системы (3). Нормализация осуществляется квазитождественным преобразованием
Уi = Zi + ^ ^ г = 1,...,п, (4)
чет
после которого получаем систему (3) в нормальной форме
^ = ф^г) = XiZi + г = 1,...,п. (5)
(ч,А)=о чет
Отметим, что в виду квазитождественности, преобразование (4) не меняет линейную часть системы.
Важное отличие (3) и (5) состоит в ограничении слагаемых суммы в правой части нормализованной системы лишь членами, удовлетворяющими уравнению
п
А) = £ ЬX, =0. (6)
3 = 1
Коэффициенты Ы и д в (4) и (5) определяются по рекуррентной формуле:
п
9^4 + ^ Л) • Ы,Ч = ^ + 6а) • Ы Р • Ззг + Фi,q, (7)
] = ! Р+г=ч
р,геи1м1 чет
где второе суммирование в правой части производится по всем векторам р, г е \JiNi, удовлетворяющим условию р + г = q е и Фi,ч — коэффициент множителя х^4 в полиноме в (2), аргументы которого были преобразованы при
помощи (4). Здесь ||р|| и ||г|| < где ||q|| = +... + дп, поэтому (7) является рекуррентной формулой.
Неоднозначность в (7) обычно устраняется наложением условий
Ы,ч = 0, если Л) = 0,
Яч = 0, если (q, Л) = 0. (8)
Такое нормализующее преобразование называется «отмеченным», а слагаемые, для которых выполняется условие (6) — резонансными [16].
Теорема 1 (см. [5]). Всегда существует формальное преобразование (4), приводящее систему (3) к её нормальной форме (5).
4.2. Сходимость нормализующего преобразования
Заметим, что суммы в (4) и (5) обычно включают в себя бесконечное число членов, хотя сумма в (3) может быть конечной. Свойства сходимости бесконечных сумм исследовались в работах [5,6,16].
Условие А. В нормальной форме (5) коэффициенты удовлетворяют условию
9з = Ч2) + ^зЪ(ъ), 3 = 1,...,n,
где а(ъ) и Ь(ъ) — некоторые степенные ряды, а черта означает комплексное сопряжение. Пусть также
п
шк = шш |(я, Л)| по q е М, (я, Л) = 0, ^ (ц < 2fc, к = 1, 2,....
3 = 1
ж
Условие ш (условие на малые знаменатели). Ряд ^ 2-klogШk > —то, сходит-
к= 1
ся.
Это условие выполняется почти для всех векторов Л. Во всяком случае оно выполняется для всех алгебраических собственных чисел.
Теорема 2 (см. [5]). Если вектор Л удовлетворяет условию ш и нормальная форма (5) удовлетворяет условию А, тогда нормализующее преобразование (4) сходится.
4.3. Локальные семейства периодических решений
В настоящей работе нас будут интересовать локальные семейства периодических решений, то есть периодические семейства, которые при некоторых начальных условиях включают в себя неподвижную точку (стягиваются в эту точку), причём предполагается также, что решения из этих семейств ограничены в фазовом пространстве, иными словами движение финитно. Пусть
(9)
представляют собой степенные ряды ъ без постоянных членов. Если они сходятся в некоторой окрестности начала координат ъ = 0, тогда решения системы
(ъ) = 0, 3 = 1,...,8 (10)
образуют локально аналитическую систему решений.
В работе [16] было показано, что локально аналитические системы периодических решений системы уравнений (3) могут быть найдены с помощью её нормальной формы (5). А именно, для нормальной формы (5) определим её формальную систему уравнений Л как
Л = {ъ : фг = АíZiШ, если И,еА^ = 0;
Хг = 0, если И,еА^ = 0; % = 1,...,п},
где ш — произвольный степенной ряд по переменным системы, который не зависит от номера г. Степенной ряд ^ тот же, что и в (5). Существование системы Л гарантирует выполнение условия А. Если все мнимые части собственных значений системы (2) соотносятся как рациональные числа (т.е. в резонансном случае с одной основной частотой) и Л = 0, то локальная формальная система Л является аналитической и содержит только периодические решения. Так как
идеал этой системы имеет сходящийся базис, соответствующая система (11) имеет смысл системы уравнений в сходящихся степенных рядах. С этой точки зрения в упоминавшемся выше резонансном случае нормальная форма адекватна первоначальной системе, по крайней мере на решениях системы Л. Эти решения содержат все периодические локальные семейства, имеющиеся вблизи рассматриваемой неподвижной точки, причём точность аппроксимации может быть задана произвольно выбором длины отрезков рядов.
Если не все собственные значения рационально совместимы, нам следует разбить систему Л на такие подсистемы, чтобы в каждой из них все прочие координаты, собственные числа которых рационально несовместимы с числами выбранной подсистемы, обращались бы в нуль. Каждая из этих подсистем — аналитическая система решений. Поэтому в фазовом пространстве система Л может иметь несколько компонент (ветвей). Каждая компонента может иметь свою собственную частоту ш.
Общий случай собственных значений, а также определение аналитической системы решений, которая содержит локальные семейства условно периодических решений, были рассмотрены в книге [16]. Возможно также применение подхода без строгого соблюдения условия единой резонансной частоты. Получающееся при этом приближение следует рассматривать не как аналитическое, но как гладкое [10]. Результаты вычислений, например, для случая двойного маятника [26], показывают неплохое количественное соответствие с численными расчётами.
5. Основной алгоритм
Алгоритм вычисления д и Ы в (4), (5) основан на формулах (7) и (8). Удобно выбрать компьютерное представление отрезков рядов gi,q и Ы,ч таким образом, чтобы их отдельные элементы были сгруппированы в однородные по степеням переменных подгруппы, когда в каждой подгруппе суммарная степень переменных (без учёта степеней параметров) равнялась бы к = 1,... ,т, для каждого % по отдельности. Далее можно вычислять д и Ы порядка п + 1, используя только ряды д и Ы порядка п, то есть используя (7) как рекуррентную формулу.
Алгоритм:
Пусть п размерность системы уравнений. Чтобы осуществить нормализацию до порядка т, следует осуществить следующие шаги (1). для г = 1, 2,... ,п :
Вычислить все квадраты в у элементах в правых частях нелинейностей Ф^у) в (2), то есть вычислить подгруппу первого порядка = 1) элементов
ряда fi,q в (3) и рассортировать их на два набора в зависимости от величины скалярного произведения (6). Первый набор, для которого это произведение равно нулю, будет подгруппой первого порядка элементов gi, а второй набор после деления на величину соответствующего скалярного произведения будет подгруппой первого порядка элементов Ы. (п). для к = 2, 3,... ,т :
(a) для г = 1,2,... ,п :
Вычислить ряд до порядка нелинейных членов из ФДу) в (2), для которых подстановка у выполняется при помощи (4) до порядка к — 1 и определить коэффициенты при мономах в качестве Дч.
(b) для г = 1,2,... ,п :
Вычислить отрезки рядов gi и hi до порядка к, подразделяя ряд Дч на два набора. После этого мы можем пополнить ряд gi до порядка к и ряд hi до порядка к без вклада от суммы в правой части (7).
(c) для г = 1,2,... ,п : для ] = 1, 2,...,п :
Дополнить ряд hi порядка к теми произведениями всех элементов рядов Ы,р и таких, что их полный порядок, т.е. ||р + г|| = к. Не все
эти произведения нужно в действительности вычислять, так как коэффициент (pj + 6i,j) для некоторых значений индекса j может быть равен нулю. Отметим, что сумма в правой части (7) никогда не даёт вклада в ряд д.
Алгоритмическая сложность вышеприведённого алгоритма невысока по сравнению со сложностью вычисления правой части нелинейной системы и поэтому нами не оценивалась. При данных обстоятельствах очень важно подсчитывать правые части очень экономно, используя, насколько это возможно, тот факт, что нам нужно вычислять на каждом шаге (ii) только члены Фi порядка к, а все члены более низких порядков при последующих операциях не изменяются. Проблема оптимизации этих вычислений представляет собой одно из наиболее важных ограничений при разработке автоматической генерации кодов вычисления правой части.
6. Компьютерная реализация символьных вычислений в методе нормальной формы
Вычисление коэффициентов нормальной формы (5) и соответствующее преобразование (4) при помощи формул (7) и (8) были реализованы как пакет программ NORT [27,28], написанных на языке RLISP, версии языка Standard LISP. Более ранние попытки вычисления высоких порядков нормальной формы с использованием верхнего уровня языка REDUCE [29] к успеху не привели. Пакет NORT содержит в настоящее время около 2000 операторов. Это пакет процедур для работы с отрезками многомерных степенных рядов, не содержащих свободных членов. В дополнение к процедурам для арифметических операций с рядами имеются специальные процедуры для построения нормальных форм, а также процедуры для подстановок, для вычисления корней (когда это возможно), для дифференцирования, для печати, для обращения многомерных степенных рядов и т.д. Он содержит также специальные процедуры для вычисления величин Ляпунова.
Комплексные численные коэффициенты рядов в пакете NORT могут обрабатываться по выбору четырьмя различными арифметиками: рациональной, модулярной, с плавающей точкой и приближённо-рациональной. Есть также несколько вариантов для вывода результатов на языке системы REDUCE. Заметим, что время на сборку мусора в LISPе при работе программы оценивается менее чем в 3% общего времени вычислений, что характеризует NORT как программу с достаточно хорошей внутренней организацией.
Самым длинным результатом, полученным при помощи пакета NORT, было вычисление нормальной формы 19-го порядка в системе с двумя малыми переменными и пятью параметрами, что заняло около 6.5 часов работы 3 ГГц Pentium IV процессора при объёме оперативной памяти в 2 Гбайт. Результирующая нормальная форма содержала 1174 члена, а нормализующее преобразование 226145 слагаемых [25].
К сожалению, к настоящему времени пакет NORT не имеет графического интерфейса, поэтому был создан пакет [30-32], написанный на входном языке системы MATHEMATICA. Этот пакет также содержит набор программ для работы с отрезками многомерных степенных рядов без свободных членов. Выражения в обоих пакетах могут содержать параметры, не являющиеся малыми величинами. Сравнение пакета MATHEMATICA с пакетом NORT показывает, что вычисления в системе MATHEMATICA более гибкие и удобные, но и значительно более медленные, чем в NORT.
7. Схема исследования нелинейных систем методом
нормальной формы
Предлагаемая здесь схема исследования системы обыкновенных дифференциальных уравнений методом нормальной формы выглядит так:
1) приведение системы к модельной форме (с полиномиальной правой частью без свободных членов). Очевидно, что таких форм несколько - вблизи каждой из неподвижных точек. Исследование нужно проводить в окрестности каждой;
2) линейная нормализация системы, то есть редукция линейной части в правых частях уравнений системы к жордановой форме, и исследование соответствующих линейных частей, в том числе поиск «резонансных» значений собственных чисел, то есть таких значений, при которых возникают группы, с рационально совместимыми чисто мнимыми собственными значениями. Систему следует исследовать для каждой такой группы. В случаях, когда пара собственных чисел близка к резонансным значениям, следует ввести дополнительно новую, постоянную по времени переменную так, чтобы получилась система в резонансе, как это сделано в примере исследования уравнения Ван дер Поля ниже;
3) нелинейная нормализация системы, то есть построение нормализующего преобразования и соответствующей нормальной формы;
4) бифуркационный анализ системы по параметрам с помощью анализа исчезающих и появляющихся при изменении параметров низших порядков соответствующих нормальных форм;
5) вычисление в виде формул, приближённо описывающих семейства периодических и условно периодических решений, содержащих в себе неподвижные точки, т.е. локальных решений;
6) нормализованная система имеет порядок меньше, чем исходная, поэтому возможно повторное исследование полученной системы по этой же схеме (вторичная нормализация) с дальнейшим понижением порядка.
8. Примеры исследования систем второго порядка
Ниже в этой секции мы будем обсуждать только двумерный случай. Вопросы сходимости и интегрируемости нормальных форм уравнений второго порядка типа (2) в окрестности неподвижной точки были детально исследованы А.Д. Брю-но [5,6,16] для случая, когда стационарная точка есть элементарная особая точка, то есть когда оба собственных значения системы не равны нулю одновременно. Результат исследования сходимости можно сформулировать следующим образом.
Пусть Л2 = 0 и Л = Л 1/ Л2. Если Im( Л) = 0 или если Л > 0, тогда преобразование сходится, и мы можем произвести аппроксимацию решения первоначальной системы из известных интегралов нормальной формы с помощью преобразования (4) с любой желаемой точностью. Другими словами, случаи «узла» и «фокуса» можно исследовать без каких-либо дополнительных условий. При отрицательном иррациональном Л сходимость будет иметь место, если для всех ненулевых векторов q с целыми элементами существуют положительные величины е, и, такие, что | (q, А) | > е(|q1l + |q2\)-w. Данное условие может быть проверено до построения нормальной формы. Этот случай есть особый случай седловой точки. Наконец, для действительных неположительных рациональных Л = —m/n ^ 0 мы будем иметь сходимость лишь при некоторых дополнительных требованиях на нормальную форму. Этот интересный случай включает так называемые случаи центра и предельного цикла. Рассмотрим пару соответствующих примеров.
Это уравнение второго порядка, которое берет своё начало из задачи о математическом маятнике:
где р — угол отклонения маятника. Простейшей модельной системой с полиномиальной правой частью для него будет кубическая система. Единицы измерения
8.1. Уравнение Дуффинга
(12)
могут быть выбраны так, что Шо = 1. Заменяя далее переменную р = \/Ъх, получим уравнение Дуффинга. Отметим, что если обсуждать наличие в задаче малого параметра, то его роль здесь играют малые начальные значения
^ = -ж + ж3. (13)
Это гамильтоново уравнение с полной энергией
Н = 2 ^ +2- 4(14)
Оно имеет три неподвижные точки х = —1, 0,1. Рассмотрим окрестность точки ж = 0.
Вводя комплексные переменные
Ж = У1 + У2, = г(У1 — У2), (15)
получаем запись (13) в диагонализованной форме
^ = т — 2 Ы + У2 )3, ¿Ц = —гУ2 + 2 (У1 + У2)3. (16)
Заметим, что в этой паре уравнений слагаемые правых частей могут быть получены друг из друга перестановкой у\ о у2 и комплексным сопряжением коэффициентов. Такая ситуация будет иметь место всегда, когда исходное уравнение имеет действительные коэффициенты.
Вектор собственных значений линейной части системы: Л = {I, —г}. В соответствии с определением нормальной формы (5), мы будем иметь в суммах правых частей этой формы только слагаемые, в которых (Л, р) = г (р1 — р2) = 0, то есть только члены со степенями переменных, для которых р\ = р2:
dzг , , 2 2 \
— = ггг + 21(31,1,1 ■ гг ■ + д1,2,2 ■ ^ ■ + ...),
¿Г (17)
^¿2 . , 2 2л
— = —1X2 + 22(32,1,1 ■ Х1 ■ %2 + д2,2,2 ■ г2 ■ г2 + ...).
Условие (11) для уравнения второго порядка с собственными значениями Х1 = —Х2 имеет вид
91,к,к ■ ■ г2)к = — 92,к,к ■ ■ г2)к.
к=1,... к=1,...
Можно показать, что для любого уравнения второго порядка с действительными коэффициентами имеет место равенство
91,г,г = 92АА, Ь,1,г,3 = ^ 2,3,г, (18)
поэтому условие А здесь имеет вид
^ Ве(д1,к,к) ■ (^1 ■ Х2)к = 0. (19)
к=1,...
После вычисления на компьютере нормальной формы для уравнения Дуффинга можно увидеть, что (19) выполняется автоматически, так как все д1¿,г и
§2,г,г для этого уравнения оказываются чисто мнимыми. Это означает, что для уравнения Дуффинга мы имеем ситуацию, называемую «центром», при которой периодическое решение существует на некотором многообразии начальных условий. В нашем случае — при достаточно малых начальных условиях.
Действительно, умножая первое уравнение (17) на и второе на Х\, и затем складывая их, получим ё(гхг^/йЪ = 0. Таким образом, система (17) имеет семейство решений
^(¿) = сге+1ш(с^с2)1, %(*) = С2 е-ш(с1 ^, (20)
где ш(г1 ■ г2) == 1 + д1,1,1/г + д1,2,2/г + ... — действительные постоянные, а с1, с2 — константы интегрирования.
Теперь мы можем получить приближение к решениям исходного уравнения (13), подставляя найденные выражения для в у^ (4), а после этого в х. Если мы выберем комплексно сопряжённые величины для С1 = С2, то получим аппроксимацию действительного семейства периодических решений в форме отрезка ряда Фурье с частотой ш(с 1 ■ С2) и коэффициентами, в виде отрезков степенных рядов от действительного параметра 1 ■ 2. Такие ряды являются частным случаем рядов Пуассона. Как несложно видеть из (20), фаза комплексного параметра С1 есть половина сдвига по времени, который является одной из констант интегрирования в решениях любой автономной системы ОДУ (трансляционная инвариантность решений автономных систем).
Мы получили эти ряды как ряды по переменной С1 ■ С2. Это неудобно, поэтому для конечного представления выразим Н в виде ряда по С1 ■ С2. Подставим затем найденные х и ёх/ё£ в (14). Обращая эти выражения согласно теореме о неявной функции и подставляя найденное выражение 1 ■ 2 как разложение по Н в выражения для ш их, получим окончательный результат как отрезки рядов по Н. Для экономии места мы приведём результат только до членов пятого порядка:
ш = 1 - 3н - 69н 2 - 633Н 3 - Н0421Н 4 - 5 +...,
4 64 256 16384 65536 + ,
х = /рШ х
" . / 9 271 т2 10779 т3 243613 4 2963587 Л
собШ ) 1 +--Н +--Н 2 +--Н3 +--Н4 +--Н5
^ ; V 16 256 4096 32768 131072 )
,п 1 3 тт 1209 т2 127233 , 6907221 Л
- с»^ )Н (16 + Т6Н + ШН 2 + б5736"Н 3 + 1т5гёН 4) +
+ еоз(5ш,)Н2 + 21н + + 255^3) -
+ ( ) \256 + 512 + 32768 + 65536 )
- )Н3 (Ш + ^ + 5ИИ +
+ cos(9wi )Н 4 ( + 21 Н ) - cos(11wi )Н5 (—1- ) +
v ; 1 65536 131072 J v ; V1°48576/
(21)
Мы опустили здесь константу сдвига по времени.
Результаты этих вычислений были проверены двумя способами. Во-первых, прямой подстановкой разложения (21) в первоначальное уравнение (13). После такой подстановки мы нашли в исходном уравнении лишь члены пренебрежимо малых порядков по Н. Второй способ проверки состоял в сравнении численных решений уравнения Дуффинга, полученных методом Рунге-Кутта 4-го порядка (процедура d02af пакета NAG), со значениями разложений, табулированных при различных значениях величины Н.
Принимая во внимание (14), имеем Н = Н(х = р/л/6, dж/dí = Выберем для проверки значения Н, не превышающие 0,163. Физически такие значения не малы, поскольку даже в случае, когда максимум угловой амплитуды маятника р достигается в верхней точке, то есть в точке р = ж/2 и d^/dt = 0, значение Нтах ~ 0,163, что нетрудно вычислить по формуле (14).
Введём функцию максимальной относительной ошибки в течение одного периода колебаний
г / (^series ^num)2 + (¿Жseries/di d$num/dt)2 /оо\
J err = sup W-2-ГТл-7"Т7л2-• (22)
ХПит + (¿Жпит/di)2
Эта функция даёт значения максимума относительного различия в фазовом пространстве между численными значениями приближения (21) (£series) с одной стороны и численными решениями уравнения (13) (жпит) с другой, при различных значениях полной энергии Н. Имеем
ferr(H = 0, 1) ~ 1, 8 X 10-8,
ferr(H = 0,125) ~ 7,4 X 10-7, ferr (Hmax = 0, 163) ~ 7, 4 X 10-5.
Видно, что достигнутая точность может быть использована для практических целей в физической области энергий. Заметим, что при вычислениях в реальном времени использование заранее полученных для нужных решений приближенных формул иногда может быть предпочтительнее численного интегрирования дифференциальных уравнений.
8.2. Уравнение Ван дер Поля
Это негамильтоново уравнение возникает в задачах о колебательных процессах в электрических цепях, при обсчёте ряда биологических моделей и при описании некоторых химических реакций
Л = + - )|. (23)
После линейной комплексной замены переменных (15) оно может быть приведено к диагональному виду. Заметим, что при этом величина £ входит в линейную часть диагонализованной системы. Собственными значениями её линейной части будут
Лх = (е2 - л/е4 - 4)/2, Л2 = (е2 + л/е4 - 4)/2.
Нам, однако, желательно получить систему в форме, подобной представлению уравнения Дуффинга выше, то есть в резонансной форме. Для этого следует представить £ как новую переменную (повысить размерность системы). Этот приём часто позволяет освободить собственные значения от зависимости от параметров и привести их к резонансному случаю. Таким образом, имеем
^ = гух + ^ - У2)[£2 - (Ух + У2)2},
dy2 . 1, чг 2 , de
= -гу2 + 2(У2 - У1П£ - (Ух + У2) }, ^ =0.
Изучим ниже окрестность нулевой неподвижной точки.
Как и для уравнения Дуффинга, эти уравнения имеют комплексно сопряжённые коэффициенты в членах, получающихся при перестановке у1 О у2. Также сопряжены при такой замене и пары ^ и Н (18).
Как и в уравнении Дуффинга, сумма в правой стороне нормальной формы будет включать только члены, где р1 = Р2:
1 =г^ + ¿1 (д 1,1,1,2 к ■ (¿1 ■ ¿2)£2к + 91,2,2,2к ■ (¿1 ■ ¿2)2£2к + ...),
к=0,1,. (24)
— = -гХ2 + *2 (92,1,1,2к ■ (21 ■ ¿2)е2 к + д2,2,2,2к ■ (¿1 ■ ^)2£ 2 к + . . .).
Третье «дополнительное» уравнение меняться не будет: ёе/— = 0. Условие (19) принимает теперь форму
йе ^ д1,з,з,2к£ 2к( ■ Z2У = 0, 31,0,0,0 = 0. (25)
3,к=0,1,...
В отличие от уравнения Дуффинга, это условие не удовлетворяется автоматически, но благодаря теореме о неявных функциях может быть разрешено в форме
2 к
1 ■ 2 = к 2 к. к=1,2,...
Можно заметить, что если это равенство удовлетворяется, то произведение 1 ■ 2 постоянно во времени. Поэтому мы можем продолжить вычисления таким же образом, как и в случае уравнения Дуффинга, но теперь свободна лишь константа интегрирования, отвечающая за сдвиг по времени. Это случай «предельного цикла». Условие (25) определяет его предельную периодическую траекторию, к которой решение стремится вне зависимости от начальных условий
Х1 (¿) = С1е+ш( с1^с2)4,
г2(Ь) = С2 е-ш(с1 , 2 к
С1 ■ С2 = 1к£ .
к=1,2,...
Подставляя далее найденные выше z^ в (4), получаем у^, а затем х. Если мы выберем комплексно сопряжённые величины для С1 = С2, мы получим так же, как для уравнения Дуффинга аппроксимации действительных решений в форме отрезков рядов Пуассона. Имеем
Ш = 1 - 1 £4 + Л_ £8 + £12 - 678899 £16 + ...
16 3072 884736 5096079360
, 1 4 23 8 51619 12 948555443 16\ С08(Ш) ' 2 + 64е - 49152£ - 169869312£ + 19568944742400£ +
2
,п . 4 I 3 101 4 24061 8 279818087 12\
+ еов(3ш^) £ 4---1--е 4 +--е8--е12 +
4 7 1 32 12288 28311552 815372697600 )
^ 5 1865 4 328835 8 111998015 12\
+ С08(5Ш^)£ 1 -96 + 110592е - 254803968е - 293534171136£ / +
х
8/ 1379 10923199 4 21049213549 Л + cos( ш^е ^ 110592 - 3185049600£ + 183458856960000£ ) +
8 f 61 1769369 4 161113663733 Л + cos(9cot)e ^20480 - 589824000£ + 237817036800000£ ) +
ч 12 ( 409871 1359229760383 Л
+ cos(11wt)e12---1--еМ +
v ; V 331776000 1872809164800000 )
12 f 715247 2076538440769 Л
+ COS(13wi)£ V- 3715891200 + 5243865661440000 £) +
, 16 ( 526426361 \ . 16 ( 392636471 \
+ COs(15wi)£ { 4661213921280/ + COs(1^)e [ 2996494663680^ +
2 f 1 15 4 779 8 4538017 12
+ кп(3ш)£ ( -4 + в - Ц79648£ - 6794772480£ ' + ( 85 8095 4 1252495
+ sin(5wt)e6---£4--е8 I +
v ; \2304 1327104 6115295232 1
7 99967 / 415949513
+ sin(7^)e%576 - т327Г040£4 + 382205952000£) +
+ sin(9wt)в 10 (-+ 117258703 в 4^ + + ( ) V 2457600 + 70778880000 ) +
• ч 10 ( 5533 1657839733 Л
+ Sin(11wi)£ {- 7372800 + 1486356480000£ ) +
. .1П . 14 ( 21731177 \
+ sin(13wi )£14 - +
v ; V57802752000)
. ,лг , 14 ( 138697 \
+ sin(15wi)£ \2774532096/ +
(26)
Здесь мы также положили произвольный сдвиг по времени равным нулю.
Вычисления в пакете NORT до 32 порядка по е заняло на PentiumPro-200 компьютере (200МГц) 1,5 минуты. Мы получили при этом 145 членов в каждой сумме нормальной формы (24) и 1773 члена для нормализующего преобразования от Zi к yi. Вычисленные выражения для частот содержат 9 слагаемых. Заметим, что степенной ряд для частоты уравнения Ван дер Поля уже вычислялся до 164 порядка по £ в работе [33].
Сравнение нашего результата с численным в терминах (22) даёт
ferr(e2 =0,5) ~ 8 X 10-10, ferr(e2 = 0, 75) ~ 4 X 10-8, ferr(e2 = 1,0) ~ 1 X 10-5.
Кроме решения уравнения мы можем также получить выражения для точек, которые могут быть выбраны в качестве начальных условий, лежащих непосредственно на траектории предельного цикла в уравнении (23). Это можно сделать, вычисляя степенной ряд по е, обратный к ряду (4):
ж = 0,
d;r _ / 17 4 _ 1577 8 _ 102956839 12 48722480822161 16 ~dt = £'\ + 96£ - 552960£ - 55738368000£ + 157315969843200000£ +
•••)
2 + 96 F' * -
1033 552960
е8 +
1019689 55738368000
е12 +
9835512276689 157315969843200000
е16 +
-х = 0. -
Обратим внимание на то, что траектория полученного предельного цикла при устремлении параметра е к нулю сжимается в точку. Говорят также, что возникновение предельного цикла с ростом этого параметра из нуля есть «бифуркация рождения предельного цикла» или бифуркация Хопфа [24].
и
х
9. Заключение
Таким образом, мы можем заключить, что построение высших порядков нормальных форм позволяет строить конечные формулы для количественных аппроксимаций периодических решений систем автономных нелинейных обыкновенных дифференциальных уравнений.
Мы также использовали этот метод для исследования периодических решений системы двойного маятника [26] и для оценок цикличности в планарной кубической системе [34] в связи с шестнадцатой проблемой Гильберта.
Литература
1. Arnold V. I., Anosov D. V. Dynamical Systems I (Encyclopaedia of Mathematical Sciences). — New York: Springer-Verlag, 1987.
2. Guckenheimer J., Holmes P. Nonlinear Oscillations. Dynamical Systems and Bifurcations of Vector Fields. — New York: Springer-Verlag, 1986.
3. Deprit A. Canonical Transformation Depending on a Small Parameter // Celestial Mechanics. — 1969. — Vol. 1, No 1. — Pp. 12-30.
4. Hori G. I. Theory of General Perturbations with Unspecified Canonical Variables // Publications of the Astronomical Society of Japan. — 1966. — Vol. 18, No 4. — Pp. 287-296.
5. Брюно А. Д. Аналитическая форма дифференциальных уравнений. I // Труды московского математического общества. — 1971. — Т. 25. — С. 119-262.
6. Брюно А. Д. Аналитическая форма дифференциальных уравнений. II // Труды московского математического общества. — 1972. — Т. 26. — С. 199-239.
7. Брюно А. Д. Ограниченная задача трех тел. Плоские периодические орбиты. — Москва: Наука, 1990.
8. Bruno A. D. Normal Forms // Mathematics and Computers in Simulation. — 1998. — Vol. 45. — Pp. 413-427.
9. Брюно А. Степенная геометрия в алгебраических и дифференциальных уравнениях. — Москва: Наука-Физматлит, 1998.
10. Bibikov Y. N. Local Theory of Nonlinear Analytic Ordinary Differential Equations. — New York: Springer-Verlag, Lect. Note Math., 1979. — Vol. 702.
11. Mersman W. A. A New Algorithm for Lie Transformation // Celestial Mechanics. — 1970. — Vol. 3. — Pp. 81-89.
12. Shevchenko I. I., Sokolsky A. G. Algorithms for Normalization of Hamiltonian Systems by Means of Computer Algebra // Computer Physics Communications. — 1993. — Vol. 77. — Pp. 11-18.
13. Godziewski K., Maciejewski A. J. System for Normalization of a Hamiltonian Function Based on Lie Series // Celestial Mechanics and Dynamical Astronomy. — 1990. — Vol. 49. — Pp. 1-10.
14. Ito H. Convergence of Birkhoff Normal Forms for Integrable Systems // Commen-tarii Mathematici Helvetici. — 1989. — Vol. 64. — Pp. 412-461.
15. Ito H. Integrability of Hamiltonian Systems and Birkhoff Normal Forms in the Simple Resonance Case // Mathematische Annalen. — 1992. — Vol. 292. — Pp. 411— 444.
16. Брюно А. Д. Локальный метод нелинейного анализа дифференциальных уравнений. — Москва: Наука, 1979.
17. Walcher S. On Differential Equations in Normal Form // Mathematische Annalen. — 1991. — Vol. 291. — Pp. 293-314.
18. Walcher S. On Transformations into Normal Form // Journal of Mathematical Analysis and Applications. — 1993. — Vol. 180. — Pp. 617-632.
19. Vallier L. An Algorithm for the Computation of Normal Forms and Invariant Manifolds // Proceedings of the International Symposium on Symbolic and Algebraic Computation, July 1993, Kiev, Ukraine. — New York: ACM Press, 1993. — Pp. 225-233.
20. Martin C. F., Zhou Y. Carleman Linearization of Linear Systems with Polynomial Output: Techrep 9 / Institut Mittag-Leffler. — Sweden, 2002-2003.
21. Verhulst F. Nonlinear Differential Equations and Dynamical Systems. — Berlin Heidelberg New York: Springer-Verlag, 1989. — 277 p.
22. Bruno A. D. Bifurcation of the Periodic Solutions in the Case of a Multiple Pair of Imaginary Eigenvalues // Selecta Mathematica Formerly Sovietica. — 1993. — Vol. 12. — Pp. 1-12.
23. Marsden J. E., McCracken M. The Hopf Bifurcation and its Applications. — New York: Springer Applied Mathematics Series, 1976. — Vol. 19, 408 p.
24. Hassard B., Kazarinoff N. D., Wan Y. H. Theory and Applications of Hopf Bifurcation. — Cambridge: Cambridge University Press, 1981. — 280 p.
25. Bruno A. D., Edneral V. F. On Possibility of Additional Solutions of the Degenerate System Near Double Degeneration at the Special Value of the Parameter // Proceedings of the 15th International Workshop Computer Algebra in Scientific Computing. September 9-13, 2013. Berlin, Germany. — Vol. 8136. — Cham Heidelberg: Springer-Verlag, series: LNCS, 2013. — Pp. 75-87.
26. Edneral V. F., Khanin R. Investigation of the Double Pendulum System by the Normal Form Method in MATHEMATICA // Programming and Computer Software. — 2004. — Vol. 30, No 2. — Pp. 115-117.
27. Edneral V. F., Khrustalev O. A. Propagation of Electromagnetic Waves in Thin-Film Structures with Smoothly Irregular Sections // Proceedings of International Conference on Computer Algebra and its Application in Theoretical Physics. September 1985. Dubna, USSR. — Dubna: JINR, 1985. — Pp. 219-224.
28. Edneral V. F., Khrustalev O. A. Program for Recasting ODE Systems in Normal Form // Sov. J. Programmirovanie. — 1992. — No 5. — Pp. 73-80.
29. Hearn A. C. REDUCE. User's Manual. — Berkeley: Rand Publication, CP87, 1987.
30. Edneral V. F., Khanin R. Multivariate Power Series and Normal Form Calculation in Mathematica // Proceedings of the Fifth Workshop on Computer Algebra in Scientific Computing. September, 2002, Big Yalta, Ukraine. — Munich: Tech. Univ. München, 2002. — Pp. 63-70.
31. Edneral V. F., Khanin R. Application of the Resonant Normal Form to High Order Nonlinear ODEs using MATHEMATICA // Nuclear Instr. and Methods in Physics Research, A. — 2003. — Vol. 502, No 2-3. — Pp. 643-645.
32. Edneral V. F. On Algorithm of the Normal Form Building // Proceedings of the 10th International Workshop on Computer Algebra in Scientific Computing. September 16-20, 2007. Bonn, Germany. — Vol. 4770. — Munich: SpringerVerlag, series LNCS, 2007. — Pp. 134-142.
33. Andersen G. M., Geer J. F. Power Series Expansions for the Frequency and Period of the Limit Cycle of the Van der Pol Equation // SIAM Journal on Applied Mathematics. — 1983. — Vol. 42. — Pp. 678-693.
34. Edneral V. F. Computer Evaluation of Cyclicity in Planar Cubic System // Proceedings of the ISSAC'97. July, 1997. Hawaii, USA. — New York: ACM Press, 1997. — Pp. 305-309.
EflHepan B.O., THMoi^eeBCKaa O.^. noacK ceMeäcTB nepHOflHHecKHX .
43
UDC 519.61+517.925
Looking for Families of Periodic Solutions of Ordinary Differential Equations Systems by Normal Form Method.
Part I
V. F. Edneral*, O. D. Timofeevskaya^
* Skobeltsyn Institute of Nuclear Physics Lomonosov Moscow State University (SINP MSU) 1(2), Leninskie gory, GSP-1, Moscow, Russian Federation, 119991 Department of Applied Informatics and Probability Theory Peoples' Friendship University of Russia 6, Miklukho-Maklaya str., Moscow, Russian Federation, 117198 ^ Physical Department of Lomonosov Moscow State University 1(2), Leninskie gory, Moscow, Russian Federation, 119991
In this paper, we discuss the application of resonant normal form method to the search of periodic solutions families of autonomous systems of explicit ordinary differential equations with polynomial nonlinearities in the right parts. Further, using formulated by Prof. A.D. Bruno sufficient convergence condition for the normalizing transformation, we find local families of periodic solutions of systems of such ODE in the vicinity of stationary points. In this unified approach both Hamiltonian and not Hamiltonian systems are investigated.
For reasons of volume the article is divided into two parts. In the first part we describe an algorithm of implementing the method of normal forms. Software packages created by the authors are briefly described separately. We have developed a RLISP language package for working in REDUCE system, and for MATHEMATICA system a package on the external language of this system. This packages allow us, in particular, to obtain formulas describing local (containing a fixed point) families of periodic solutions. The results of calculations are presented in the form of Fourier series segments of a given length with frequency and coefficients themselves calculated as parameter series segments. This representation corresponds to the special case of segments of Poisson series. It is important that using a single algorithm, one can study both two-dimensional and higher-order systems. The second part is devoted to fourth-order systems.
The comparison of tabulation of formulas obtained with numerical solutions of the corresponding equations shows good quantitative agreement. The approach described can be used for modeling of physical and biological systems.
Key words and phrases: resonant normal form, dynamical systems, local periodic families of solutions, computer algebra.
References
1. V. I. Arnold, D. V. Anosov, Dynamical Systems I (Encyclopaedia of Mathematical Sciences), Springer-Verlag, New York, 1987.
2. J. Guckenheimer, H. P., Nonlinear Oscillations. Dynamical Systems and Bifurcations of Vector Fields, Springer-Verlag, New York, 1986.
3. A. Deprit, Canonical Transformation Depending on a Small Parameter, Celestial Mechanics 1 (1) (1969) 12-30.
4. G. I. Hori, Theory of General Perturbations with Unspecified Canonical Variables, Publications of the Astronomical Society of Japan 18 (4) (1966) 287-296.
5. A. D. Bruno, Analytical Form of Differential Equations. I., Transactions of the Moscow Mathematical Society 25 (1971) 131-288.
6. A. D. Bruno, Analytical Form of Differential Equations. II, Transactions of the Moscow Mathematical Society 26 (1972) 199-239.
7. A. D. Bruno, The Restricted 3-Body Problem. Plane Periodic Orbits, Vol. 17, Berlin, New York, 1994.
8. A. D. Bruno, Normal Forms, Mathematics and Computers in Simulation 45 (1998) 413-427.
9. A. D. Bruno, The Power Geometry in Algebraic and Differential Equations, Vol. 17, Amsterdam, 2000.
10. Y. N. Bibikov, Local Theory of Nonlinear Analytic Ordinary Differential Equations, Vol. 702, Springer-Verlag, Lect. Note Math., New York, 1979.
11. W. A. Mersman, A New Algorithm for Lie Transformation, Celestial Mechanics 3 (1970) 81-89.
12. I. I. Shevchenko, A. G. Sokolsky, Algorithms for Normalization of Hamiltonian Systems by Means of Computer Algebra, Computer Physics Communications 77 (1993) 11—18.
13. K. Godziewski, A. J. Maciejewski, System for Normalization of a Hamiltonian Function Based on Lie Series, Celestial Mechanics and Dynamical Astronomy 49 (1990) 1-10.
14. H. Ito, Convergence of Birkhoff Normal Forms for Integrable Systems, Commentarii Mathematici Helvetici 64 (1989) 412-461.
15. H. Ito, Integrability of Hamiltonian Systems and Birkhoff Normal Forms in the Simple Resonance Case, Mathematische Annalen 292 (1992) 411-444.
16. A. D. Bruno, Local Method in Nonlinear Differential Equations, Springer-Verlag, Berlin, 1989.
17. S. Walcher, On Differential Equations in Normal Form, Mathematische Annalen 291 (1991) 293-314.
18. S. Walcher, On Transformations into Normal Form, Journal of Mathematical Analysis and Applications 180 (1993) 617-632.
19. L. Vallier, An Algorithm for the Computation of Normal Forms and Invariant Manifolds, in: Proceedings of the International Symposium on Symbolic and Algebraic Computation, July 1993, Kiev, Ukraine, ACM Press, New York, p. 225.
20. C. F. Martin, Y. Zhou, Carleman Linearization of Linear Systems with Polynomial Output, Tech. Rep. 9, Institut Mittag-Leffler, Sweden (2002-2003).
21. F. Verhulst, Nonlinear Differential Equations and Dynamical Systems, SpringerVerlag, Berlin Heidelberg, New York, 1989.
22. A. D. Bruno, Bifurcation of the Periodic Solutions in the Case of a Multiple Pair of Imaginary Eigenvalues, Selecta Mathematica Formerly Sovietica 12 (1993) 1-12.
23. J. E. Marsden, M. McCracken, The Hopf Bifurcation and its Applications, Vol. 19, Springer Applied Math. Series, New York, 1976.
24. B. D. Hassard, N. D. Kazarinoff, Y. H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981.
25. A. D. Bruno, V. F. Edneral, On Possibility of Additional Solutions of the Degenerate System Near Double Degeneration at the Special Value of the Parameter, in: Proceedings of the 15th International Workshop Computer Algebra in Scientific Computing. September 9-13, 2013. Berlin, Germany, Springer-Verlag, series: LNCS, Cham Heidelberg, p. 75.
26. V. F. Edneral, R. Khanin, Investigation of the Double Pendulum System by the Normal Form Method in MATHEMATICA, Programming and Computer Software 30 (2) (2004) 115-117.
27. V. F. Edneral, O. A. Khrustalev, Propagation of Electromagnetic Waves in Thin-Film Structures with Smoothly Irregular Sections, in: Proceedings of International Conference on Computer Algebra and its Application in Theoretical Physics. September 1985. Dubna, USSR, JINR publ., Dubna, p. 219.
28. V. F. Edneral, O. A. Khrustalev, Program for Recasting ODE Systems in Normal Form, Sov. J. Programmirovanie (5) (1992) 73-80.
29. A. C. Hearn, REDUCE. User's Manual, Rand Publication, CP87, Berkeley, 1987.
30. V. F. Edneral, R. Khanin, Multivariate Power Series and Normal Form Calculation in Mathematica, in: Proceedings of the Fifth Workshop on Computer Algebra in Scientific Computing. September, 2002, Big Yalta, Ukraine, Tech. Univ. München, Munich, p. 63.
31. V. F. Edneral, R. Khanin, Application of the Resonant Normal Form to High Order Nonlinear ODEs using MATHEMATICA, Nuclear Instr. and Methods in Physics Research, A 502 (2-3) (2003) 643-645.
32. V. F. Edneral, On Algorithm of the Normal Form Building, in: Proceedings of the 10th International Workshop on Computer Algebra in Scientific Computing.
EflHepan B.O., THMoi^eeBCKaa O.^. noacK ceMencTB nepwogHHecKHX .
45
September 16-20, 2007. Bonn, Germany, Vol. 4770, Springer-Verlag, series LNCS, Munich, 2007, pp. 134-142.
33. G. M. Andersen, J. F. Geer, Power Series Expansions for the Frequency and Period of the Limit Cycle of the Van der Pol Equation, SIAM Journal on Applied Mathematics 42 (1983) 678-693.
34. V. F. Edneral, Computer Evaluation of Cyclicity in Planar Cubic System, in: Proceedings of the ISSAC'97. July, 1997. Hawaii, USA, ACM Press, New York, 1997, pp. 305-309.
