CC BY
12
ISSN 1810-0198. Вестник Тамбовского университета. Серия: естественные и технические науки
Том 24, № 126 2019
( Пчелинцев А.Н., Полуновский А.А., Юханова И.Ю., 2019 DOI 10.20310/1810-0198-2019-24-126-187-203 УДК 519.622.2
Метод гармонического баланса для отыскания приближённых периодических решений системы Лоренца
Александр Николаевич ПЧЕЛИНЦЕВ1 , Андрей Андреевич ПОЛУНОВСКИЙ2 , Ирина Юрьевна ЮХАНОВА1
1 ФГБОУ ВО «Тамбовский государственный технический университет» 392000, Российская Федерация, г. Тамбов, ул. Советская, 106 ORCID: https://orcid.org/0000-0003-4136-1227, e-mail: pchelintsev.an@yandex.ru ORCID: https://orcid.org/0000-0002-8339-0459, e-mail: irina_yu_10@mail.ru 2 ФГБОУ ВО «Московский государственный технический университет им. Н.Э. Баумана (национальный исследовательский университет)» 105005, Российская Федерация, г. Москва, 2-я Бауманская ул., 5 ORCID: https://orcid.org/0000-0002-6557-3649, e-mail: apap2009@yandex.ru
The harmonic balance method for finding approximate periodic
solutions of the Lorenz system
Alexander N. PCHELINTSEV1 , Andrey A. POLUNOVSKIY2 , Irina Yu. YUKHANOVA1
1 Tambov State Technical University 106 Sovetskaya St., Tambov 392000, Russian Federation ORCID: https://orcid.org/0000-0003-4136-1227, e-mail: pchelintsev.an@yandex.ru ORCID: https://orcid.org/0000-0002-8339-0459, e-mail: irina_yu_10@mail.ru 2 Bauman Moscow State Technical University 5 Baumanskaya 2-ya, Moscow 105005, Russian Federation ORCID: https://orcid.org/0000-0002-6557-3649, e-mail: apap2009@yandex.ru
Аннотация. Рассматривается метод гармонического баланса для нахождения приближённых периодических решений динамической системы Лоренца. При разработке программного обеспечения, реализующего описываемый метод, был выбран математический пакет Maxima. Показаны недостатки символьных вычислений для получения системы нелинейных алгебраических уравнений относительно циклической частоты, постоянных членов и амплитуд гармоник, составляющих искомое решение. Для ускорения расчётов впервые эта система была получена в общем виде. Приведены результаты вычислительного эксперимента - коэффициенты тригонометрических полиномов, приближающих найденное периодическое решение, начальное условие и период цикла. Полученные результаты были проверены с помощью описанного ранее в работах авторов высокоточного метода интегрирования, основанного на аппроксимации степенными рядами.
Ключевые слова: система Лоренца; аттрактор; метод гармонического баланса; ряды Фурье
Для цитирования: Пчелинцев А.Н., Полуновский А.А., Юханова И.Ю. Метод гармонического баланса для отыскания приближённых периодических решений системы Лоренца // Вестник Тамбовского университета. Серия: естественные и технические науки. Тамбов, 2019. Т. 24. № 126. С. 187-203. DOI 10.20310/1810-0198-2019-24-126-187-203.
Abstract. We consider the harmonic balance method for finding approximate periodic solutions of the Lorenz system. When developing software that implements the described method, the math package Maxima was chosen. The drawbacks of symbolic calculations for obtaining a system of nonlinear algebraic equations with respect to the cyclic frequency, free terms and amplitudes of the harmonics, that make up the desired solution, are shown. To speed up the calculations, this system was obtained in a general form for the first time. The results of the computational experiment are given: the coefficients of trigonometric polynomials approximating the found periodic solution, the initial condition, and the cycle period. The results obtained were verified using a high-precision method of numerical integration based on the power series method and described earlier in the articles of the authors.
Keywords: Lorenz system; attractor; harmonic balance method; Fourier series For citation: Pchelintsev A.N., Polunovskiy A.A., Yukhanova I.Yu. Metod garmonichesko-go balansa dlya otyskaniya priblizhonnykh periodicheskikh resheniy sistemy Lorentsa [The harmonic balance method for finding approximate periodic solutions of the Lorenz system]. Vestnik Tambovskogo universiteta. Seriya: estestvennye i tekhnicheskie nauki - Tambov University Reports. Series: Natural and Technical Sciences, 2019, vol. 24, no. 126, pp. 187-203. DOI 10.20310/1810-0198-2019-24-126-187-203. (In Russian, Abstr. in Engl.)
Введение
Рассмотрим нелинейную систему дифференциальных уравнений, введённую Э. Ло-ренцом в работе [1],
x 1 = a(x2 — x\),
x 2 = rx1 — x2 — x1x3, (0.1)
x 3 = x1x2 — bx3,
где a =10 , r = 28 , b = 8/3 (эти значения параметров системы (0.1) теперь называют классическими, а саму систему — системой Лоренца).
В статье [1] для системы Лоренца с классическими параметрами доказано следующее утверждение: существует такое C > 0, что для любого решения
X(t) = [xi(t) x2(t) x3(t)]T
начиная с некоторого момента времени становится справедливым неравенство |X(t) | <C, и дивергения векторного поля скоростей системы (0.1) отрицательна всюду в R3. Тогда [1] существует предельное множество — аттрактор Лоренца, — к которому притягиваются все траектории динамической системы при t ^ то. Таким образом, аттрактор определяет поведение решений динамической системы на больших отрезках времени.
У. Такер в работе [2] доказал гиперболичность аттрактора в системе (0.1), т. е. аттрактор состоит из траекторий, всюду плотных на нём (континуум седловых циклов), вдоль которых близкие траектории экспоненциально разбегаются; это и создает их хаотическое поведение. Тогда (как отмечает Д.В. Аносов в послесловии к книге [3, с. 285]) в аттракторе системы (0.1) может существовать бесконечное число асимптотически устойчивых периодических траекторий, но их область притяжения может быть достаточно малой (трудно улавливаемой в численном эксперименте).
Как известно (см., например, [4,5]), символическую динамику используют для отслеживания циклов в системе Лоренца. Разбивают область в фазовом пространстве, содержащую аттрактор, на конечное число подобластей. Обозначая каждый элемент разбиения буквой, траектории на аттракторе, проходящие через соответствующие области, кодируются последовательностями таких символов. Если в последовательности имеется регулярность — повторяемость групп символов, — то соответствующая траектория считается циклом. Однако возвращаемость траектории в некоторую окрестность своей части не говорит о её замкнутости. Критику результатов подобных вычислительных экспериментов можно найти, например, в [6].
В 2004 г. Д. Вишванат опубликовал работу [7], в которой привел начальные условия и периоды для трёх циклов в аттракторе Лоренца с достаточно большой точностью. Алгоритм вычислений основан на методе Линдштедта-Пуанкаре (ЛП), на который (в отличие от методов численного интегрирования) не влияет устойчивость цикла, к которому строятся приближения. Полученные в [7] вычислительные данные можно проверить, решая задачу Коши высокоточными численными методами (см, например, [9]).
Анализ работ [7,8] Д. Вишваната показал, что автор приводит общее описание алгоритма без ссылок на программную реализацию (в МАТЬАБ, как указано в его статьях). При этом не ясно, как для ЛП-метода символьно решается получаемая неоднородная линейная система дифференциальных уравнений с периодическими коэффициентами (например, для уравнения Ван дер Поля это сделать можно без особых проблем). Таким образом, актуальной задачей остаётся разработка алгоритма поиска циклов системы (0.1), детальное описание его реализации, получение начальных значений и периода цикла с заданной точностью.
Целью данной работы является отыскание приближённых периодических решений в системе Лоренца на основе метода гармонического баланса, являющегося более простым в реализации, чем ЛП-метод, используемый в [7,8]. При этом будет получена в общем виде система нелинейных алгебраических уравнений относительно циклической частоты, постоянных членов и амплитуд гармоник, составляющих искомое решение.
1. Метод гармонического баланса
Попытки построить приближённые периодические решения системы (0.1) предпринимались и до Д. Вишваната (см., например, [10]) методом гармонического баланса, но с малой точностью представления вещественных чисел, при этом в статье [10] не указаны начальные условия и периоды найденных циклов (приведены только рисунки с цикла-
ми). Сейчас этот метод активно развивается в работах А. Луо [11-13] для отыскания периодических решений нелинейных систем дифференциальных уравнений.
Будем использовать метод гармонического баланса для получения приближений к периодическим решениям системы (0.1). Для этого сделаем аппроксимацию фазовых координат на периоде Т тригонометрическими полиномами в общем виде с неизвестной циклической частотой ш (поскольку мы не знаем значение Т; в общем случае оно может быть иррациональным числом):
н
~ XI(¿) = £1,0 + ^^ (см сов(гш£) + вт(гш£)),
г=1
н
£2(£) « £2^) = £2,0 + ^ (с2,г сов(гш£) + 32,г 8т(гш£)),
г=1
н
Хэ(£) « Хэ(^) = £3,0 + ^ (сз,г сов(гш£) + 3з,г вт(гш£)),
г=1
где к — заданное количество гармоник. Если г > к, то мы полагаем
С1,г = 3М = С2,г = 32,г = Сз,г = 3з,г = 0. (1.2)
В силу правой части системы (0.1) составим невязки
¿1СО = £1(0 - ^СО - £1(^)1,
$2(£) = £'2(£) - [г£1 (¿) - £2^) - £1 (¿)£з(^)1,
5з(1) = £'з(*) - [£1(*)£2(*) - 6£з(*)1,
где штрихом переобозначена производная функции по времени. Если производить вычисления в аналитическом виде, то для каждой невязки нужно следующее:
1. Продифференцировать по времени соответствующий тригонометрический полином.
2. Где имеются произведения фазовых координат, перемножить соответствующие тригонометрические полиномы, преобразовав при этом произведения тригонометрических функций в суммы.
3. Привести подобные слагаемые для каждой функции еов() и в1п() с соответствующим аргументом.
4. В силу равенств (1.2), отсечь от полученной невязки гармоники более высокого порядка.
5. Приравнять полученную невязку к нулю, т. е. коэффициенты при её гармониках.
Если собрать в единое целое найденные алгебраические уравнения для каждой невязки, то получим пока ещё незамкнутую систему нелинейных уравнений относительно неизвестных амплитуд См , , с2,^, , сз,^ и (г = 1,Н), постоянных
членов Х1,о , ж2,о и ж3,о и циклической частоты ш. Количество неизвестных в системе равно 3(1 + 2к) + 1 = 6к + 4, а уравнений — на единицу меньше.
Дополнительное уравнение можно получить исходя из следующих соображений. Известно (см. [5,7]), что искомые циклы пересекают плоскость, проходящую через положения равновесия системы (0.1)
О1 (-/б(г - 1), -/б(г - 1), г - 1) , О2 (/б(г - 1), /Ь(г - 1), г - 1) (1.3)
и параллельную плоскости х1 Ох2 (сечение Пуанкаре). Таким образом, третья координата в начальном условии для искомых циклов равна величине г - 1 , откуда Х3(0) = г - 1. Тогда дополнительное уравнение системы имеет вид:
н
Х3,0 + ^ Сз,г - 27 = 0.
г=1
Других дополнительных сведений о периодических решениях системы Лоренца авторы не встречали. Заметим, что для трёх циклов, найденных Д. Вишванатом, в начальном условии для третьей координаты было взято число 27. Далее приведем пример системы уравнений при к = 2:
ш51,1 - 10с2,1 + 10с1,1 = 0, - 10з2д + 10в1,1 - с1,1ш = 0, 2ш«1,2 - 10С2,2 + 10С1,2 = 0,
-10«2,2 + 10в 1,2 - 2С1,2Ш = 0, 10х1,0 - 10х2,0 = 0,
. . «1,1«3,2 «1,2«3,1 . . С1,1С3,2 С1,2С3,1 ос п
С1,1Хз,о + С3,1Х1,0 +--2--1--2--+ ш«2,1 +--2--1--2--+ С2,1 - 28с1,1 = 0,
. с1,1 «3,2 с1,2«3,1 , . с3,1«1,2 С3,2«1,1 «1,1X3,0 + «3,1X1,0 +--2---2--+ «2,1 +--2---2--28^1,1 - С2,1Ш = 0,
, «1,1«3,1 0 с1,1с3,1
С1,2Х3,0 + С3,2Х1,0--2--+ 2ш«2,2 +--2--+ С2,2 - 28С1,2 = 0,
«1,2X3,0 + «3,2X1,0 +--2--+ «2,2 - 28«1,2 +--^--2С2,2Ш = 0,
, «1,2«3,2 , «1,1«3,1 , С1,2С3,2 , С1,1С3,1
Х1,0Х3,0 + Х2,0 - 28X1,0 +--2--+--2--'--2--'--2— = 0,
. «1,1«2,2 «1,2«2,1 . 8С3,1 С1,1С2,2 С1,2С2,1 _ -см^о- c2,lxl,о + ш«3,1--2---2—' 3---2---2— =
8«3
,1 с1,1«2,2 , с1,2«2,1 с2,1«1,2 , с2,2«1,1 п -«1,1X2,0 - «2,1X1,0 +—з---2—1—2---2—'—2--с3,1ш = 0,
, о , «1,1«2,1 , 8С3,2 С1,1С2,1 _
£1,2X2,0 - ^2,2X1,0 + 2ш«3,2 +--2--+ _3---2— = 0,
8«3
,2 с1,1 «2,1 с2,1 «1,1 -«1,2X2,0 - «2,2X1,0 +--3---2---2--2С3,2Ш = 0,
8X3,0 «1,2«2,2 «1,1«2,1 С1,2С2,2 С1ДС2Д _
Т"- Xl,0X2,0 - ""2 2 2 = 0,
Xз,о + С3,1 + С3,2 - 27 = 0.
Отметим, что для любого h подобная система имеет решения
xi,0 = Х2,о = ±\/b(r - 1), Хз,0 = r - 1, cfc,i = 0, sfc,i = 0, ш — любое число, k = 1, 3, i = 1, h,
соответствующие указанным положениям равновесия (1.3).
Таким образом, полученная нелинейная система алгебраических уравнений имеет неединственное решение. Для отыскания её приближённых решений будем использовать численный метод Ньютона, сходимость которого к нужному решению (т. е. отличному от положения равновесия) зависит от выбора начального приближения.
2. Символьные вычисления для получения системы алгебраических
уравнений
Итак, для получения приближения к периодическому решению мы должны получить нелинейную систему относительно неизвестных коэффициентов разложения и частоты. Как показано в п. 1, даже для двух гармоник система имеет громоздкий вид. Поэтому рассмотрим алгоритм проведения символьных вычислений для её получения.
При разработке программного обеспечения был выбран математический пакет Maxima. Программа получения амплитуд и постоянных членов невязок при h = 2 представлена далее.
/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [wxMaxima: input start ] */
display2d:false$
x1:x10+c1c1*cos(1*omega*t)+s1c1*sin(1*omega*t)+
c1c2*cos(2*omega*t)+s1c2*sin(2*omega*t)$
x2:x20+c2c1*cos(1*omega*t)+s2c1*sin(1*omega*t)+
c2c2*cos(2*omega*t)+s2c2*sin(2*omega*t)$
x3:x30+c3c1*cos(1*omega*t)+s3c1*sin(1*omega*t)+
c3c2*cos(2*omega*t)+s3c2*sin(2*omega*t)$
assume(omega >0)$
delta1:trigreduce(diff(x1,t)-(10*(x2-x1)),t)$ delta2:trigreduce(diff(x2,t)-(28*x1-x2-x1*x3),t)$ delta3:trigreduce(diff(x3,t)-(x1*x2-8/3*x3),t)$ expand(diff(delta1,cos(1*omega*t))); expand(diff(delta1,sin(1*omega*t))); expand(diff(delta1,cos(2*omega*t))); expand(diff(delta1,sin(2*omega*t)));
expand(integrate(delta1,t,0,2*°/„pi/omega)*omega/(2*°/„pi)); expand(diff(delta2,cos(1*omega*t))); expand(diff(delta2,sin(1*omega*t))); expand(diff(delta2,cos(2*omega*t)));
expand(diff(delta2,sin(2*omega*t)));
expand(integrate(delta2,t,0,2*°/„pi/omega)*omega/(2*°/„pi)); expand(diff(delta3,cos(1*omega*t))); expand(diff(delta3,sin(1*omega*t))); expand(diff(delta3,cos(2*omega*t))); expand(diff(delta3,sin(2*omega*t)));
expand(integrate(delta3,t,0,2*°/„pi/omega)*omega/(2*°/„pi)); /* [wxMaxima: input end ] */
Выражение display2d:false$ выключает многострочное рисование дробей, степеней и т. п. Знак $ позволяет вычислить результат выражения, но не выводить на экран (вместо ;). Функция trigreduce(выражение,t) свёртывает все произведения тригонометрических функций относительно переменной t в комбинации сумм. Дифференцирование невязок по гармоническим функциям необходимо для получения соответствующих амплитуд. Функция expand(выражение) раскрывает скобки (выполняет умножение, возведение в степень, приводит подобные слагаемые).
Для нахождения постоянных членов невязок применяется интегрирование на периоде, т. е. постоянный член k -ой невязки равен
Чтобы при символьном интегрировании пакет не задавал вопрос о знаке частоты, дается команда assume(omega > 0)$.
Файл с командами пакета формируется аналогично для любого количества h гармоник. После выполнения данной программы пакет выведет в консоли символьные выражения для левой части системы алгебраических уравнений, которая будет решаться в нём же методом Ньютона.
Заметим, что самая затратная по времени операция здесь — символьное интегрирование. Например, для 120 гармоник время формирования системы — более 2-х суток. Здесь можно распараллелить вычислительный процесс на три компьютера, но это значительного эффекта не даст. Поэтому систему алгебраических уравнений нужно формировать сразу. Далее получим общий вид этой системы.
Отметим, что при решении системы нелинейных уравнений методом Ньютона матрица Якоби для левой части системы не обращается — в пакете Maxima используется LU-разложение для решения системы линейных уравнений на каждой итерации метода.
3. Общий вид системы алгебраических уравнений
Поскольку правая часть системы (0.1) содержит нелинейности в виде произведений фазовых координат, получим соотношения, выражающие коэффициенты тригонометрических полиномов, получаемых при умножении приближений Xi(t)Xs(t) и Xi(t)x2(t).
2п
Рассмотрим две функции f (¿) и ^(¿) , представимые рядами Фурье
те
f (¿) = а0 + ^^ (а, сов(гш£) + 6, в1п(гш£)),
,= 1 те
^(¿) = А0 + ^ (А, сов(гш^) + В, в1п(гшг)).
,= 1
Пусть
те
f (¿) = а0 + ^^ (а, сов(гш^) + в в1п(гш£)).
г=1
Следуя книге [14, с. 123-125], имеем следующие соотношения:
1те
а0 = а^Л + — ^ (атАт + 6тВт),
т=1
1те
а, = а^А + 1 ( ) + 6т(Вт+, + Вт-,)) , (3.4)
т=1
1те
= а0В, + — ^ ] (ат(Вт+, - Вт_,) - 6т(Ат+, - Ат-,)) . (3.5)
2 ' ^ т=1
Будем предполагать, что при г > к
а, = = А, = В, = 0.
Поскольку для нашей задачи мы ищем приближение до к -ой гармоники включительно, занулим все амплитуды в произведении при г > к , т. е.
а, = в = 0.
Таким образом, мы перейдем от произведения рядов к произведению тригонометрических полиномов.
Также в соотношениях (3.4) и (3.5) мы будем предполагать [14, с. 124], что
Ат_, то Вт_, В,—В0
Тогда получим:
1н
а0 = а^Л + — ^ (атАт + 6тВт),
т=1
1 ^ 1 ^ 1 ^ 1 ^
аг а0Аг + 2 ^ у атАт+г + 2 ^ ^ атАт—г + 2 ^ ^ Вт+г + 2 ^ у ЬтВт—г = т=1 т=1 т=1 т=1
1 Н—г 1 1 г—1 1 Н
= а0Аг + 2 ^ у атАт+г + ^агА0 + 2 ^ у атАг—т + 2 ^ у атАт—г + т=1 т=1 т=г+1
1 Н 1 1 г—1 1 Н
+ 2 ^ у ЬтВт+г + 2ЬгВ0 2 ^ у ЬтВг—т + 2 ^ у ЬтВт—г
т=1 т=1 т=г+1
1 Н—г 1 г—1
= а0Аг + агА0 + 2 ^ у (атАт+г + ЬтВт+г) + 2 ^ у (атАг—т - ЬтВг—т) + т=1 т=1
1 Н
+ 2 ^ у (атАт—г + ЬтВт—г)
т=г+1
1 ^ 1 ^ 1 ^ 1 ^
вг а0Вг + 2 ^ у атВт+г 2 ^ у атВт—г 2 ^ у ЬтАт+г + 2 ^ у ЬтАт—г
т=1 т=1 т=1 т=1
1 Н—г 1 г—1 1 Н
а0Вг + 2 ^ у атВт+г + 2 ^ у атВг—т 2 ^ у атВт—г
т=1 т=1 т=г+1
1 Н—г 1 г—1 1 Н
2 ^ у ЬтАт+г + ЬгА0 + 2 ^ у ЬтАг—т + 2 ^ у ЬтАт—г т=1 т=1 т=г+1
1 Н—г 1 г—1
аоВг + ЬгАо +
т=1 т=1
1 Н
+ 2 ^ у ( атВт—г + ЬтАт—г) .
т=г+1
Применяя полученные формулы для вычисления произведений тригонометрических полиномов к невязкам, мы можем записать уравнения для г -ых гармоник (г = 1, к — номер гармоники, к =1, 3 — номер невязки):
к = 1:
гш«1,г - 10с2,г + 10с1,г = 0, -гшс1,г - 10«2,г + 10«1,г = 0,
уравнение, соответствующее постоянному члену для первой невязки,
Xl,о - X2,0 = 0,
к = 2:
н,
гшв2,, - —8сМ + с2,, + £1,0с3,, + с1,,£з,0 + 2 ^ (с1,тсз,т+, + в1,т3з,т+,) +
2
т=1
1 ,— 1
+ — (с1
,т сз,, т - 31 ,т т )+
2
т=1
1н
+ — ^ у (с1,тсз,т_, + s1,ms3,m_г) 0,
т=,+ 1
1 н_ ,
-гшс2,, - 2831,, + 32,, + £1,0вз,, + «М£з,0 + 1 ^ (^^3,™+* - «1,тсз ,т+,) +
2
т=1
1 ,— 1
+ — (с1
,т т + 31 ,т сз,, )+
т=1
н
+ 1 ^ ^ ( c1,ms3,m_г + s1,mc3,m_г) 0,
т=,+ 1
уравнение, соответствующее постоянному члену для второй невязки, —
1н
-—8£1,0 + £2,0 + £1,0£з,0 + — ^ (с1,тсз,т + 51,т«3,т) = 0,
т=1
к = 3:
н_ ,
гш33,, - £1,0с2,, - с1 ,,£2,0 - 1 ^ ] (с1,тс2,т+, + s1,ms2,m+г) -
т=1
1 г_1
— ^ ^ (с1,тс2,,_т s1,ms2,г_m) т=1
— ^ у (с1,тс2,т_, + s1,ms2,m_г) + З3с3,, 0, т=,+1
1 н_ ,
-гшс3,, - £1,032,, - s1,г£2,0 - 1 ^ ] (с1,т32,т+, - 31,тс2,т+,) -
т=1 1 г_1
— ^ ^ (c1,ms2,г_m + s1,mc2,г_m) т=1
1 н 8
— / у ( c1,ms2,m_г + s1,mc2,m_г) + 333,, 0, т=,+1
уравнение, соответствующее постоянному члену для третьей невязки, —
-Х1,0Х2,0 — 2 ^^ (с1,тс2,т + + 88Х3,0 = 0,
т=1
дополнительное уравнение системы —
н
Х3,0 + ^ Са,г - 27 = 0.
г=1
4. Результаты вычислительного эксперимента
В результате многочисленных вычислительных экспериментов было подобрано начальное приближение для циклической частоты, постоянных членов и амплитуд при Л = Л = 5 :
Ш = 4, Х1,0 = Х2,0 = Хз,0 = 0, С1,г = -1, % = 1, 5, ви = 0, 2 = 1, 3, 4, 5, в1,2 = 1.
Данный результат замечателен тем, что метод Ньютона сходится к решению, отличному от положений равновесия. Поэтому для улучшения точности приближенного периодического решения мы рассматриваем систему алгебраических уравнений для значения Л , равного некоторому Л2 > Л1. Полученное численное решение системы при Л = Л1 берётся как начальное приближение для амплитуд с индексами % < Л1 у системы с Л = Л2 , а значения начального приближения для амплитуд с индексами % > Л1 полагаются равными нулю.
В таблицах 1-3 приведён результат решения системы при Л = 35, точность метода Ньютона - 10-8. Значение периода получается равным Т = 1.558652210, начальное условие для полученного приближённого периодического решения -
Х1(0) = -2.147367631, Х2(0) = 2.078048211, Х3(0) = 27. (4.6)
Начальные значения (4.6) были проверены на периоде в компьютерной программе, реализующей численное интегрирование системы (0.1) модифицированным методом степенных рядов [9] с точностью оценки общего члена ряда 10 -25 , 100 бит под мантиссу вещественного числа и машинным эпсилон 1.57772 • 10-30 . При таких параметрах метода приближённые значения фазовых координат, полученные с помощью численного интегрирования, были также проверены тем же численным методом, но в обратном времени. Значения в обратном времени совпадают с (4.6) до 9-го знака включительно после точки. Результирующие же значения х1(Т) , х2(Т) и х3(Т) совпадают с (4.6) до 8-го знака включительно.
Таблица 1
Амплитуды гармоник для £1(^) , £1,0 = 0
г С1,г
1 —5.780478259196228 8.56017654325353
2 0 0
3 3.160762628380509 2.239212141102876
4 0 0
5 0.6958870387616096 -0.7979388979225431
6 0 0
7 -0.1891992374027477 -0.1864921358925765
8 0 0
9 -0.04770429623010056 0.04554044367245914
10 0 0
11 0.01112322884679491 0.01209138588669679
12 0 0
13 0.003061207095371694 -0.002735092350544739
14 0 0
15 -6.744578887916229 • 10-4 -7.748319471034087 • 10-4
16 0 0
17 -1.960718247379475 • 10-4 1.665584161919807 • 10-4
18 0 0
19 4.116738805347028 • 10-5 4.960493476144467 • 10-5
20 0 0
21 1.254757391175977 • 10-5 -1.018054283421179 • 10-5
22 0 0
23 -2.518375902000733 • 10-6 -3.173486439630506 • 10-6
24 0 0
25 -8.025338211960923 • 10-7 6.230623750431923 • 10-7
26 0 0
27 1.541534734542893 • 10-7 2.0292802821633 • 10-7
28 0 0
29 5.130649139299358 • 10-8 -3.813725452268523 • 10-8
30 0 0
31 -9.43393531993558 • 10-9 -1.297038481588497 • 10-8
32 0 0
33 -3.278552746800046 • 10-9 2.333260259021725 • 10-9
34 0 0
35 5.76957885768651 • 10-10 8.28626640138045 • 10-10
Таблица 2
Амплитуды гармоник для Хг^) , х2,0 = 0
г С2,г «2,г
1 —2.32972926505593 10.89038310357172
2 0 0
3 5.86875317198698 —1.5832552129833
4 0 0
5 —0.9124249133801483 —2.200556873678218
6 0 0
7 —0.7154457265566421 0.3473932955614448
8 0 0
9 0.1175186702136983 0.2186139734768588
10 0 0
11 0.06473984670858603 —0.03723215039412078
12 0 0
13 —0.01127208646321726 —0.01877739524860192
14 0 0
15 —0.005359671824365359 0.003303445299126894
16 0 0
17 9.453499475830811 • 10-4 0.001510235036151227
18 0 0
19 4.211022386354685 • 10-4 —2.657049331814368 • 10-4
20 0 0
21 —7.363528144366622 • 10-5 — 1.164013765469982 • 10-4
22 0 0
23 —3.19419300699788 • 10-5 2.017609175377016 • 10-5
24 0 0
25 5.47663534401654 • 10-6 8.710929378319451 • 10-6
26 0 0
27 2.362852034076972 • 10-6 — 1.474901091428546 • 10-6
28 0 0
29 —3.94532524722541 • 10-7 —6.379296603810031 • 10-7
30 0 0
31 — 1.715198229248314 • 10-7 1.049218598356554 • 10-7
32 0 0
33 2.776045093375681 • 10-8 4.59473450493284 • 10-8
34 0 0
35 1.22681173575872 • 10-8 —7.31171826830086 • 10-9
Таблица 3
Амплитуды гармоник для я3(£) , х3,0 = 23.04210397942006
г С3,г «3,г
1 0 0
2 7.568410271550653 -9.50386584559212
3 0 0
4 -3.555327211552558 -1.844710563805469
5 0 0
6 -0.4741220131932616 1.279043179069961
7 0 0
8 0.4227292179138024 0.1274574086305204
9 0 0
10 0.03498415351761577 -0.1315337800809524
11 0 0
12 -0.03934013541135439 -0.009645786231708874
13 0 0
14 -0.002660052258813564 0.01145537653603837
15 0 0
16 0.003271688724557337 7.33752523103949 • 10-4
17 0 0
18 2.024982256871223 • 10-4 -9.206266886554897 • 10-4
19 0 0
20 -2.560063570343799 • 10-4 -5.58964460662525 • 10-5
21 0 0
22 -1.542436654918173 • 10-5 7.050327849098175 • 10-5
23 0 0
24 1.926014222030195 • 10-5 4.25261452471065•10-6
25 0 0
26 1.170939944189529 • 10-6 -5.225643926851625 • 10-6
27 0 0
28 -1.409525591131397 • 10-6 -3.21879984959824 • 10-7
29 0 0
30 -8.83134288999026 • 10-8 3.782652721710986 • 10-7
31 0 0
32 1.010610960272394 • 10-7 2.418021923473667 • 10-8
33 0 0
34 6.606163280924149 • 10-9 -2.689431432873997 • 10-8
35 0 0
Цикл, соответствующий (4.6), показан на рис. 1. Отметим, что найденный цикл совпадает с первым циклом Вишваната в [7], все знаки после точки для величины Т также совпадают с данными из работы [7].
х3
40 35 30 25 20 15 10 5 0
х1
20
30
Рис. 1. Цикл, полученный методом гармонического баланса
Список литературы
[1] E. N. Lorenz, "Deterministic nonperiodic flow", Journal of the Atmospheric Sciences, 20:2 (1963), 130-141.
[2] W. Tucker, "A rigorous ODE solver and Smale's 14th problem", Foundations of Computational Mathematics, 2:1 (2002), 53-117.
[3] Ж. Палис, В. Ди Мелу, Геометрическая теория динамических систем, Мир, М., 1986.
[4] М. И. Рабинович, "Стохастические автоколебания и турбулентность", Успехи физических наук, 125 (1978), 123-168.
[5] Z. Galias, W. Tucker, "Validated study of the existence of short cycles for chaotic systems using symbolic dynamics and interval tools", International Journal of Bifurcation and Chaos, 21:2 (2011), 551-563.
[6] R. Lozi, "Can we trust in numerical computations of chaotic solutions of dynamical systems?", Topology and Dynamics of Chaos. In Celebration of Robert Gilmore's 70th Birthday. - World Scientific Series in Nonlinear Science Series A, 84 (2013), 63-98.
[7] D. Viswanath, "The fractal property of the Lorenz attractor", Physica D: Nonlinear Phenomena, 190:1-2 (2004), 115-128.
[8] D. Viswanath, "The Lindstedt-Poincare technique as an algorithm for computing periodic orbits", SIAM Review, 43:3 (2001), 478-495.
[9] А.Н. Пчелинцев, "Численное и физическое моделирование динамики системы Лоренца", Сибирский журнал вычислительной математики, 17:2 (2014), 191-201.
[10] K. Neymeyr, F. Seelig, "Determination of unstable limit cycles in chaotic systems by method of unrestricted harmonic balance", Zeitschrift fur Naturforschung A, 46:6 (1991), 499-502.
[11] A. C. J. Luo, J. Huang, "Approximate solutions of periodic motions in nonlinear systems via a generalized harmonic balance", Journal of Vibration and Control, 18:11 (2011), 1661-1674.
[12] A.C.J. Luo, Toward Analytical Chaos in Nonlinear Systems, John Wiley & Sons, Chichester, 2014.
[13] A.C.J. Luo, S. Guo, "Analytical solutions of period-1 to period-2 motions in a periodically diffused brusselator", Journal of Computational and Nonlinear Dynamics, 13:9 (2018), 090912.
[14] G.P. Tolstov, Fourier Series, Dover Publications, New York, 1962.
References
[1] E.N. Lorenz, "Deterministic nonperiodic flow", Journal of the Atmospheric Sciences, 20:2 (1963), 130-141.
[2] W. Tucker, "A rigorous ODE solver and Smale's 14th problem", Foundations of Computational Mathematics, 2:1 (2002), 53-117.
[3] J. Palis, Jr. W. de Melo, Geometric Theory of Dynamical Systems An Introduction, SpringerVerlag, New York, 1982.
[4] M.I. Rabinovich, "Stochastic self-oscillations and turbulence", Soviet Physics Uspekhi, 21:5 (1978), 443-469.
[5] Z. Galias, W. Tucker, "Validated study of the existence of short cycles for chaotic systems using symbolic dynamics and interval tools", International Journal of Bifurcation and Chaos, 21:2 (2011), 551-563.
[6] R. Lozi, "Can we trust in numerical computations of chaotic solutions of dynamical systems?", Topology and Dynamics of Chaos. In Celebration of Robert Gilmore's 70th Birthday. - World Scientific Series in Nonlinear Science Series A, 84 (2013), 63-98.
[7] D. Viswanath, "The fractal property of the Lorenz attractor", Physica D: Nonlinear Phenomena, 190:1-2 (2004), 115-128.
[8] D. Viswanath, "The Lindstedt-Poincare technique as an algorithm for computing periodic orbits", SIAM Review, 43:3 (2001), 478-495.
[9] A. N. Pchelintsev, "Numerical and physical modeling of the dynamics of the Lorenz system", Numerical Analysis and Applications, 7:2 (2014), 159-167.
[10] K. Neymeyr, F. Seelig, "Determination of unstable limit cycles in chaotic systems by method of unrestricted harmonic balance", Zeitschrift fur Naturforschung A, 46:6 (1991), 499-502.
[11] A. C. J. Luo, J. Huang, "Approximate solutions of periodic motions in nonlinear systems via a generalized harmonic balance", Journal of Vibration and Control, 18:11 (2011), 1661-1674.
[12] A.C.J. Luo, Toward Analytical Chaos in Nonlinear Systems, John Wiley & Sons, Chichester, 2014.
[13] A. C. J. Luo, S. Guo, "Analytical solutions of period-1 to period-2 motions in a periodically diffused brusselator", Journal of Computational and Nonlinear Dynamics, 13:9 (2018), 090912.
[14] G.P. Tolstov, Fourier Series, Dover Publications, New York, 1962.
Информация об авторах
Пчелинцев Александр Николаевич, кандидат физико-математических наук, заведующий кафедрой «Высшая математика». Тамбовский государственный технический университет, г. Тамбов, Российская Федерация. E-mail: pchelintsev.an@yandex.ru ORCID: https://orcid.org/0000-0003-4136-1227
Полуновский Андрей Андреевич, студент. Московский государственный технический университет им. Н. Э. Баумана (национальный исследовательский университет), г. Москва, Российская Федерация. E-mail: apap2009@yandex.ru ORCID: https://orcid.org/0000-0002-6557-3649
Юханова Ирина Юрьевна, магистр. Тамбовский государственный технический университет, г. Тамбов, Российская Федерация. E-mail: irina_yu_10@mail.ru
ORCID: https://orcid.org/0000-0002-8339-0459
Конфликт интересов отсутствует.
Для контактов:
Пчелинцев Александр Николаевич
E-mail: pchelintsev.an@yandex.ru
Поступила в редакцию 21.01.2019 г.
Поступила после рецензирования 19.03.2019 г.
Принята к публикации 20.05.2019 г.
Information about the authors
Alexander N. Pchelintsev, Candidate of Physics and Mathematics, Head of the Higher Mathematics Department. Tambov State Technical University, Tambov, the Russian Federation. E-mail: pchelintsev.an@yandex.ru ORCID: https://orcid.org/0000-0003-4136-1227
Andrey A. Polunovskiy, Student. Bauman Moscow State Technical University, Moscow, the Russian Federation. E-mail: apap2009@yandex.ru ORCID: https://orcid.org/0000-0002-6557-3649
Irina Yu. Yukhanova, Master. Tambov State Technical University, Tambov, the Russian Federation. E-mail: irina_yu_10@mail.ru ORCID: https://orcid.org/0000-0002-8339-0459
There is no conflict of interests.
Corresponding author:
Alexander N. Pchelintsev E-mail: pchelintsev.an@yandex.ru
Received 21 January 2019 Reviewed 19 March 2019 Accepted for press 20 May 2019
