The existence of chaotic regimes of the fractional analogue of the Duffing-type oscillator

Full Text

Введение. Исследование случайных, хаотических, периодических или квазипериодических режимов различных колебательных систем (осцилляторов) является важной задачей в теории колебательных систем [1-6]. Существует достаточно много критериев и тестов для определения хаотических режимов. Одни методы основываются на исследовании амплитудночастотого спектра гармонических колебаний с помощью Фурье-анализа (statistical 0-1 test) [7]. Если в результате исследований спектр является дискретным, то речь идет о периодических или квазипериодических колебаниНаучная статьяях, а если спектр непрерывный, то колебания могут быть хаотическими или случайными. Другие методы основаны на отображении Пуанкаре для исследуемой колебательной системы, которое заключается в построении сечения фазовой траектории секущей поверхностью [4]. В случае, если отображение Пуанкаре представляет собой форму облака, то это соответствует случайному колебательному процессу, если же отображение Пуанкаре имеет форму некоторой линии, то возможны квазипериодические и хаотические режимы. Известно, что хаотические колебания обладают высокой чувствительностью к малым изменениям начальных условий. Поэтому оценка скорости разбегания фазовых траекторий с помощью максимальных показателей Ляпунова является одним из надежных методов определения хаоса в рассматриваемой динамической системе. Наличие положительного максимального показателя Ляпунова является критерием определения хаотического режима для рассматриваемой системы. Настоящая работа является продолжением работы [8], где в качестве объекта исследований мы выбрали дробный аналог осциллятора типа Дуффинга - динамическую систему, которая описывает нелинейные колебания с эффектом степенной памяти. Математическая модель такой системы представляет собой задачу Коши с дифференциальными уравнениями в производных дробных порядков. Дробный аналог осциллятора типа Дуффинга, как показали исследования [9], может обладать периодическими или квазипериодическими режимами, а также бистабильным поведением [10]. Поэтому возникает вопрос: возможны ли хаотические режимы для дробного аналога осциллятора типа Дуффинга Чтобы ответить на этот вопрос, мы в настоящей работе проведем исследование дробного осциллятора Дуффинга с помощью спектров максимальных показателей Ляпунова по аналогии с работой автора [11], построенных в зависимости от значений дробных порядков производной и коэффициента трения, бифуркационных диаграмм и фазовых траекторий. 1. Основные понятия и определения. Введем некоторые понятия и определения. Определение 1. Осциллятором Дуффинга называется нелинейный осцил- лятор вида () + () () + 3 () = cos(), (0) = 0 , (0) = 0 , (1) где () 2 (0, ) - функция смещения; [0, ] - время рассматриваемого процесса; > 0 - время моделирования; > 0 - коэффициент трения; и - амплитуда и частота внешнего воздействия; 0 и 0 - заданные константы, определяющие начальные условия. Определение 2. Дробным аналогом осциллятора типа Дуффинга (1) мы будем называть нелинейный осциллятор, который имеет математическое описание в виде задачи Коши: 0 () + 0 () () + 3 () = cos(), (0) = 0 , (0) = 0 , (2) где операторы 0 () 1 = (2 ) () ( )1 1 () 0 () = , (1 ) 0 ( ) , 1 < < 2, 0 0 < < 1, понимаются в смысле Герасимова-Капуто [12, 13]. Замечание 1. Свойства дробных производных Герасимова-Капуто можно найти, например, в работах [14, 15]. Замечание 2. Заметим, что в предельном случае, когда выполнены равенства = 2 и = 1, задача Коши (2) переходит в задачу Коши для классического осциллятора Дуффинга (1). Поэтому исследуемая система (2) является обобщением классического осциллятора Дуффинга и содержит не только новые, но и ранее известные свойства. Замечание 3. Отметим, что вопросы существования и единственности решения задачи Коши (2) были рассмотрены в работах [16], а в работе [17] предложены и исследованы конечно-разностные схемы для ее численного решения. Дробный аналог осциллятора типа Дуффинга (2) удобно записать в виде неавтономной системы: 1 () = 2 () , 1 = , 0 1 2 0 2 () = cos () 2 () + 1 () 31 (), 2 = , (3) 1 (0) = 10 , 2 (0) = 20 . Замечание 4. Отметим, что в системе (3) возможны случаи, когда 2 6 1 и 2 > 1, что необходимо учитывать при построении ее численного решения. Определение 3. Динамическую систему (3) будем называть соизмеримой, если 1 = 2 = , и несоизмеримой, если 1 = 2 . Рассмотрим более общую неавтономную колебательную систему: 1 () = 1 (1 () , 2 ()) , 0 1 2 0 2 () = 2 (1 () , 2 (), ) , (4) 1 (0) = 10 , 2 (0) = 20 , где функции 1 (1 (), 2 ()) и 2 (1 (), 2 (), ) обладают нужными свойствами гладкости для существования первых частных производных. Определение 4. Системой уравнений в вариациях 1 и 2 для динамической системы (4) будем называть систему вида 1 (1 , 2 ) 1 (1 , 2 ) 1 0 1 () = 1 + 2 , 1 2 (5) 2 2 () = 1 2 (1 , 2 , ) + 2 2 (1 , 2 , ) , 0 1 или в векторной форме: 0 = () , 2 1 2 } , 1 , 0 0 = {0 = {1 , 2 } , где () - якобиан системы (4), = {1 , 2 }. Система уравнений в вариациях (5) для исходной динамической системы (3) является ключевой при построении максимальных показателей Ляпунова по алгоритму Вольфа [18]. Эти уравнения описывают эволюцию бесконечно малого возмущения фазовой траектории исходной динамической системы (4). 2. Алгоритм Вольфа. Рассмотрим алгоритм Вольфа с ортогонализацией Грама-Шмидта для построения максимальных показателей Ляпунова динамической системы (3). Алгоритм состоит из следующих этапов. 1. Выбираем начальную точку - начальный вектор 0 и вместе с ней будем отслеживать возмущенных траекторий. В нашем случае = 2. 2. Решаем численно исходное уравнение совместно с двумя комплектами возмущенных уравнений или уравнений в вариациях (5). Причем общее количество решаемых уравнений можно определить числом ( + 1), где - количество уравнений в исходной системе. В качестве начальных векторов для уравнений в вариациях (5) необходимо выбрать набор векторов 00 , 00 , которые являются ортогональными и нормированными на единицу. 3. Через время траектория перейдет в некоторую точку 1 , векторы возмущения 1 , 1 перенормируем с помощью метода Грама-Шмидта по формулам 01 = 1 , 1 ( ) 0 1 = 1 1 , 01 1 , 10 = 1 , 1 где (, ) - скалярное произведение векторов. 4. Далее продолжаем счет от точки 1 и векторов возмущений 01 , 10 . Через очередной интервал времени получаем новый набор векторов возмущений 2 , 2 , который подвергается ортогонализации и перенормировке. 5. Этапы 2-4 повторяются раз и по ходу вычислений подсчитываются суммы 1 = ln ( ), 2 = ln ( ), =1 =1 в которых фигурируют векторы возмущений до перенормировки, но после ортогонализации. 6. Оценка максимальных показателей Ляпунова вычисляется по формулам max = , = 1, 2, . . . , . Замечание 5. Процедура ортогонализации Грама-Шмидта необходима для того, чтобы исключить доминирование составляющей максимального показателя Ляпунова на больших временах при вычислении векторов вдоль фазовой траектории. В противном случае задача будет являться плохо обусловленной. Замечание 6. Отметим, что для существования аттрактора в рассматриваемой динамической системе (3) необходимо, чтобы сумма всех показателей Ляпунова была отрицательной, т.е. max < 0. =1 Замечание 7. Необходимо отметить, что алгоритм Вольфа был выше приведен для случая 2 6 1, т.е. = 2. В случае, когда 2 > 1, размерность системы (3) увеличивается на единицу, поэтому количество возмущенных траекторий будет = 3, что приводит к увеличению количества уравнений в системах (3) и (4). Замечание 8. Большое значение имеет исследование спектра максимальных показателей Ляпунова, который строится в зависимости от значений интересующего нас управляющего параметра. В качестве управляющих параметров для системы (2) могут использоваться параметры , , . Поэтому мы будем исследовать следующие зависимости максимального показателя Ляпунова: max (), max () и max (). 3. Результаты исследований. Рассмотрим дробный аналог осциллятора типа Дуффинга (3) и будем считать, что в этой динамической системе 2 6 1. Значения параметров выберем следующими: = 0.3, = 1, [0, 100], (0) = 0.2, (0) = 0.3. Рассмотрим некоторый этапы алгоритма Вольфа для построения максимальных показателей Ляпунова. Пункт 1. Так как система (3) имеет два уравнения, количество уравнений в вариациях также будет равно двум. Для них мы выберем начальные (1) (1) (2) (2) условия 1 (0) = 1, 2 (0) = 0 и 1 (0) = 0, 2 (0) = 1. Легко проверить, что эти начальные условия ортогональны и нормированы на единицу. Пункт 2. Для построения уравнений в вариациях определим якобиан для системы (3): [ ] 0 1 () = , 1 32 тогда уравнения в вариациях и совместно с системой (3) имеют вид 1 0 1 () = 2 (), 2 0 2 () = 0.3 cos() 2 () + 1 () 31 (), (1) (1) (2) (2) 1 0 1 () = 2 () , ( ) (1) (1) 2 (1) 0 2 () = 1 () 1 321 () 2 () , 1 0 1 () = 2 () , ( ) (2) (2) 2 (2) 0 2 () = 1 () 1 321 () 2 (), 1 (0) = 0.2, 2 (0) = 0.3, (1) (1) 1 (0) = 1, 2 (0) = 0, (2) (2) 1 (0) = 0, 2 (0) = 1. Вся сложность алгоритма Вольфа заключается в решении системы (6), которая точного решения не имеет. Поэтому систему (6) необходимо решать численно, например, с помощью аппроксимации операторов дробных производных и перехода к алгебраической системе уравнений. После некоторых преобразований система (6) может иметь дискретный аналог, который можно записать следующим образом: ( ) 1 1 = ( ) + 1,+1 2, 1,1 1, 1, , =1 ( ) 1 (2,1 2, ) 0.3 cos ( ) + 1, 31, , 2,+1 = 1 ( ) 2, =1 ) ( 1 (1) (1) (1) (1) (1) 1 1,+1 = 2, (2,1 2, ) + 1, , =1 ( ) 1 (1) (1) (1) (1) (1) 1 (2,1 2, )+1, (1 321, ) , 2,+1 = ( ) 2, (7) =1 ( ) 1 (2) (2) (2) (2) (2) 1,+1 = 1 2, (1,1 1, ) + 1, , =1 ( ) 1 (2) (2) (2) (2) (2) 1 2 ( ) ( )+ (1 3 ) ; = 1, 2, 1, 2,+1 2,1 2, =1 10 = 0.2, 20 = 0.3, (1) (1) 10 = 1, 20 = 0, (2) (2) 10 = 0, 20 = 1, = 1, 2, . . . , . где = 1 , (2 1 ) = ( + 1)11 11 , = 2 ; (2 2 ) = ( + 1)22 22 - весовые коэффициенты, - количество расчетных узлов; = / - шаг сетки; () - гамма-функция Эйлера. Замечание 9. Вопросы устойчивости и сходимости численного алгоритма (7) представляют самостоятельный интерес и рассматриваться в рамках этой работы не будут. Заметим, что значения параметров в системе (7) будут выбраны, чтобы исключить «жесткость» системы, т.е. мы будем рассматривать «мягкую» задачу. Решение системы (7) при выборе параметров , [ (1) значений ] [управляющих ] (2) 1,+1 1,+1 , , дает пару векторов 1 = , 2 = , которые подверга(1) (2) 2,+1 2,+1 ются процедуре ортогонализации Грама-Шмидта и нормализации согласно пунктам 3 и 4 алгоритма. Потом вычисляются суммы и сами показатели Ляпунова согласно пунктам 5 и 6. Для нашей задачи Коши (1) мы будем получать два показателя Ляпунова: 1max и 2max , так как мы имеем два уравнения ( = 2) в исходной системе (3). Сумма показателей 1max и 2max является отрицательной, что говорит о том, что система (3) диссипативна и для такой системы существует аттрактор. Расчеты по алгоритму Вольфа показали, что один из показателей Ляпунова всегда отрицательный, а другой может менять свой знак, и поэтому мы его выберем для построения спектра показателя Ляпунова в зависимости от значений управляющих параметров , , . На рис. 1 приведена бифуркационная диаграмма для спектра максимального показателя Ляпунова в зависимости от значений управляющего параметра для соизмеримой системы (3) в случае 1 = 2 = 0.9. На рис. 1, a приведен спектр максимального показателя Ляпунова, построенный в зависимости от значений параметра [0, 1]. Видно, что есть области, в которых максимальный показатель Ляпунова положительный, например, [0.0074, 0.0689] [0.1247, 0.1516] и еще одна область в окрестности точки 0.1. Это свидетельствует о хаотических режимах, а остальные случаи отрицательности максимального показателя Ляпунова соответствуют аттрактору. На рис. 1, b приведен спектр максимального показателя Ляпунова, который был построен в зависимости от параметра [1, 2]. Видно, что на всем этом промежутке максимальный показатель Ляпунова отрицательный, что соответствует аттрактору - предельному циклу. Вышесказанное подтверждается построенной для этого случая фазовой траекторией (рис. 2). Следующий спектр максимальных показателей Ляпунова был построен в зависимости от значений дробного показателя [1, 2] при значениях параметров = 1 и = 0.15 (рис.3). На рис. 3 видно, что существуют области в окрестности точки 1.8 и [1.84, 2], в которой максимальный показатель Ляпунова положителен. Поэтому в этих областях должен существовать хаотический аттрактор. Другая область [1, 1.8], в которой показатель Ляпунова принимает отрицательные значения, соответствует аттрактору - предельному циклу. Фазовые траектории приведены для этих случаев на рис. 4. На рис. 4, а и 4, b фазовые траектории соответствуют хаотическому аттрактору при значении параметров = 1.81 и = 1.9, а на рис. 4, с при значении = 1.3 мы видим предельный цикл. Другой спектр максимального показателя Ляпунова мы можем получить, если его построить в зависимости от значений дробного параметра [0, 1] при фиксированных значениях параметров = 1.1 и = 0.15. На рис. 5 приведены результаты построения спектра максимального показателя Ляпунова. Мы видим, что максимальный показатель Ляпунова во всей области изменения параметра принимает отрицательные значения. Это дает основание считать, что в этом случае существует аттрактор, вид которого представлен на рис. 6. Из него видно, что этот аттрактор является предельным циклом. Рассмотрим теперь случай, когда в исходной системе (3) выполняется неравенство 2 = > 1. В этом случае систему (3) можно записать следующим образом: 1 0 1 () = 2 () , 1 = , 2 () = () , = 1, 3 2 0 2 3 () = cos () 2 () + 1 () 31 () , 3 = 1; 0 3 1 (0) = 10 , 2 (0) = 20 , 3 (0) = 0. Заметим, что в системе (8) выполняется условие 1 + 2 + 3 = < 2. Рис. 1. Бифуркационные диаграммы спектра максимальных показателей Ляпунова для дробного аналога осциллятора типа Дуффинга в зависимости от значений параметра : кривая (a) построена для параметра [0, 1] шагом = 0.01; кривая (b) построена для параметра [1, 2] с шагом = 0.01, = 2000 [Figure 1. Bifurcation diagrams of the spectrum of Lyapunov maximum exponents for the Duffing-type fractal oscillator depending on the values of the parameter : curve (a) constructed for the parameter [0, 1] step = 0.01; curve (b) constructed for the parameter [1, 2] with the step = 0.01, = 2000] Рис. 2. Фазовые траектории, построенные согласно бифуркационным диаграммам (рис. 1) в зависимости от значений параметра : фазовая траектория (a) построена для параметра = 0.01; фазовая траектория (b) построена для параметра = 0.9; фазовая траектория (c) построенная для параметра = 1.5 [Figure 2. Phase trajectories constructed according to the bifurcation diagrams (Fig. 1) depending on the values of the parameter : phase trajectory (a) constructed for the parameter = 0.01; phase trajectory (b) constructed for the parameter = 0.9; phase trajectory (c) constructed for the parameter = 1.5] Рис. 3. Бифуркационная диаграмма для спектра максимальных показателей Ляпунова для дробного аналога осциллятора типа Дуффинга в зависимости от значений дробного параметра [1, 2] с шагом = 0.01 при = 1 и = 0.15 [Figure 3. Bifurcation diagram of the spectrum of Lyapunov maximum exponents for the Duffingtype fractal oscillator depending on the values of the fractional parameter [1, 2] with the step = 0.01 when = 1 and = 0.15] Рис. 4. Хаотические аттракторы дробного аналога осциллятора типа Дуффинга, построенные согласно бифуркационной диаграмме (рис. 3) при значении параметра = 1 и = 0.15: кривая (a) построена при = 1.81; кривая (b) построена при = 1.9; кривая (c) построена при = 1.3 (предельный цикл) [Figure 4. Chaotic attractors of the Duffing-type fractional oscillator constructed according to bifurcation diagram (Fig. 3) depending on the values of the parameters = 1 and = 0.15: curve (a) constructed for the parameter = 1.81; curve (b) constructed for the parameter = 1.9; curve (c) constructed for the parameter = 1.3 (the limit cycle)] Рис. 5. Бифуркационная диаграмма спектра максимальных показателей Ляпунова для дробного аналога осциллятора типа Дуффинга в зависимости от значений дробного параметра [0, 1] с шагом = 0.01 при = 1.1 и = 0.15 [Figure 5. Bifurcation diagram of the spectrum of maximum Lyapunov exponents for the Duffingtype fractional oscillator depending on the values of the fractional parameter [0, 1] with the step = 0.01 when = 1.1 and = 0.15] Рис. 6. Предельные циклы для дробного аналога осциллятора типа Дуффинга, построенные согласно бифуркационной диаграмме (рис. 5), при значении параметра = 1.1 и = 0.15: кривая (a) построена при = 0.9; кривая (b) построена при = 0.5; кривая (c) построена при = 0.2 [Figure 6. Limit cycles of the Duffing-type fractional oscillator constructed according to bifurcation diagram (Fig. 5) depending on the values of the parameters = 1.1 and = 0.15: curve (a) constructed for the parameter = 0.9; curve (b) constructed for the parameter = 0.5; curve (c) constructed for the parameter = 0.2] Система (8) имеет три уравнения, поэтому необходимо составить 9 уравнений в вариациях, при этом якобиан для нее имеет вид 0 1 0 0 0 1 . () = (9) 2 1 3 0 Согласно (9), а также с учетом известных значений параметров, мы приходим к следующей системе: 1 0 1 () = 2 (), 2 0 2 () = 3 (), 03 3 () = 0.3 cos() 2 () + 1 () 31 (), (1) 1 (1) 0 1 () = 2 () , (1) 3 (1) 0 2 () = 3 () , ( ) (1) (1) 3 (1) 0 3 () = 1 () 1 321 () 1 () , (2) 1 (2) 1 () = 2 () , 0 (2) 2 (2) 0 2 () = 3 () , ( ) (10) (2) (2) 3 (2) 0 3 () = 1 () 1 32 () 3 () , (3) 1 (3) 0 1 () = 2 () , (3) 2 (3) 0 2 () = 3 () , ( ) 3 (3) () = (3) () 1 32 () (3) () ; 0 3 3 3 (0) = 0.2, (0) = 0.3, (0) = 0, 1 2 3 (1) (1) (1) 1 (0) = 1, 2 (0) = 0, 3 (0) = 0, (2) (2) (2) 1 (0) = 1, 2 (0) = 0, 3 (0) = 0, (3) (3) (3) 1 (0) = 0, 2 (0) = 0, 3 (0) = 1. Система (10) решается численно, так же как и система (6), но естественно, что для этого понадобится больше вычислительных ресурсов. Далее, следуя алгоритму Вольфа, после процедуры ортогонализации Грама-Шмидта, мы получим три максимальных показателя Ляпунова: 1max , 2max , 3max . В качестве примера приведем результаты расчетов максимальных показателей Ляпунова при различных значениях дробного параметра , которые приведены в таблице. Из таблицы мы видим эволюцию наследственной динамической системы Дуффинга. При значении = 0.1 все три максимальных показателя Ляпунова положительны и их знаки (сигнатура) имеют вид < +, +, + >. Это означает, что в этом случае отсутствует аттрактор. На рис. 7 приведена фазовая траектория, которая соответствует этому случаю. Видно, что эта траектория не является ни предельным циклом, ни предельной точкой и не соответствует хаотическому режиму. Интерес представляет следующий случай, когда = 0.2. В силу того что один из максимальных показателей Ляпунова является положительным, Таблица 1 Спектр максимальных показателей Ляпунова со знаками (сигнатурами), вычисленных при фиксированных значениях параметров = 1.8 и = 0.15 и различных значениях параметра [Spectrum of Lyapunov maximum exponents and its signature calculated by fixed values of the parameters = 1.8 and = 0.15, and various values of the parameter ] no. 1max 2max 3max Signature of the spectrum 1 2 3 4 5 6 7 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.22766 0.00269 0.01652 0.00566 0.00274 0.00503 0.02708 0.1152 0.3206 0.1472 0.0827 0.0400 0.0279 0.0358 0.059 0.914 1.379 1.684 1.805 1.707 1.429 < +, +, + > < +, , > < , , > < , , > < , , > < +, , > < +, , > Рис. 7. Фазовые траектории, построенные при ванных значениях параметров = 1.8 и = 0.15 ных значениях параметра [Figure 7. Phase trajectories constructed for fixed the parameters = 1.8 and = 0.15, and various the parameter ] фиксирои различvalues of values of а сумма всех показателей отрицательна, можно сделать вывод о существовании хаотического аттрактора (сигнатура < +, , > ). Далее мы видим из таблицы, что в диапазоне изменения параметра от 0.3 до 0.5 все показатели Ляпунова отрицательны (сигнатура < , , >), что соответствует предельному циклу (рис. 7). Следующая пара значений (0.6 и 0.7) для параметра приводит к сигнатуре < +, , > для максимальных показателей Ляпунова, что соответствует хаотическому аттрактору (рис. 7). Замечание 10. Стоит отметить, что в сигнатуре спектра максимальных показателей отсутствуют нулевые значения, которые характерны для автономных динамических систем, но динамическая система (3) является неавтономной, поэтому сигнатура спектра максимальных показателей не содержит нулевых значений. Заключение. В работе с помощью расчета максимальных показателей Ляпунова выполнено численное исследование дробного аналога осциллятора типа Дуффинга с целью идентификации хаотических колебательных режимов. Расчет максимальных показателей Ляпунова проведен по алгоритму Вольфа с ортогонализацией Грама-Шмидта, построены их спектры в зависимости от значений параметров [0, 2], [1, 2], [0, 1] и фазовые траектории для исследуемых колебательных режимов. Установлено, что все максимальные показатели Ляпунова могут принимать отрицательные значения, что соответствует на фазовой плоскости предельному циклу. Также возможен случай, когда первый максимальный показатель Ляпунова положительный, а остальные отрицательные, что указывает на существование хаотического аттрактора. Поэтому для дробного аналога осциллятора типа Дуффинга возможна бифуркация: хаотический аттрактор - предельный цикл. Этот результат был подтвержден с помощью исследования соответствующих фазовых траекторий. Отметим, что построение карт динамических режимов позволяет провести более детальный анализ наследственной динамической системы (2), однако на это требуются значительные вычислительные ресурсы и такой анализ не являлся целью данной работы.

About the authors

Roman Ivanovich Parovik

Institute of Cosmophysical Researches and Radio Wave Propagation, Far East Division, Russian Academy of Sciences; Kamchatka State University named after Vitus Bering

Email: romano84@mail.ru
Candidate of physico-mathematical sciences, Associate professor 7, Mirnaya st., Paratunka, Kamchatkiy kray, 684034, Russian Federation; 4, Pogranichnaya st., Petropavlovsk-Kamchatskiy, 683032, Russian Federation

References

Supplementary files