Локализация скрытых аттракторов обобщенной системы Чуа на основе метода гармонического баланса Текст научной статьи по специальности «Математика»

ЛОКАЛИЗАЦИЯ СКРЫТЫХ АТТРАКТОРОВ ОБОБЩЕННОЙ СИСТЕМЫ ЧУА НА ОСНОВЕ МЕТОДА ГАРМОНИЧЕСКОГО БАЛАНСА

В. И. Вагайцев1, Н. В. Кузнецов2, Г. А. Леонов3

1. С.-Петербургский государственный университет, аспирант, vladimir.vagaitsev@gmail.com

2. С.-Петербургский государственный университет, канд. физ.-мат. наук, доцент, nkuznetsov239@mail.ru

3. С.-Петербургский государственный университет, д-р физ.-мат. наук, профессор, leonov@math.spbu.ru

1. Введение. Аттракторы Чуа (double scroll attractors) были введены для дифференциальных уравнений, описывающих нелинейные электрические цепи [1—5]. Эти аттракторы были получены в результате возбуждения хаотических колебаний: состояние равновесия оказывалось неустойчивым, и в результате переходного процесса траектории стремились к аттрактору Чуа. Этот вычислительный процесс легко повторяем в стандартном компьютерном эксперименте [6]. В настоящей работе описан вычислительный процесс для определения скрытых аттракторов Чуа. Эти аттракторы сосуществуют с устойчивым состоянием равновесия и не могут быть вычислены упомянутым выше способом. Приведенный здесь алгоритм вычисления скрытых аттракторов Чуа является развитием идей, предложенных в работах [7-9].

2. Алгоритм локализации скрытых колебаний. В [7, 8] были предложены методы поиска периодических решений многомерных нелинейных динамических систем со скалярной нелинейностью. Напомним основные положения этого метода.

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

dx

— = Px + qV>(r*x), (1)

где x € R", P —постоянная n x n-матрица, q и r — n-мерные векторы, * — операция транспонирования, ф(х) —непрерывная скалярная функция и ф(0) = 0. В данной работе показано, что модификация этих алгоритмов позволяет находить не только периодические решения в системах вида (1), но и производить локализацию аттракторов хаотических систем [1-6].

Для поиска периодического решения, близкого к гармоническому колебанию, введем коэффициент гармонической линеаризации к так, чтобы матрица Po = P + kqr* линейной системы

dz , ,

~dt = °Z ^

имела пару чисто мнимых собственных значений ±шо(^о > 0), а остальные ее собственные значения имели отрицательные вещественные части. Тогда систему (1) можно

переписать в виде

dx

— = P0x + q^(r*x), (3)

где ^(а) = ф(а) — ка.

© Г. А. Леонов, Н. В. Кузнецов, В. И. Вагайцев, 2010

Поскольку нас интересуют периодические решения системы (3), представляется естественным ввести конечную последовательность непрерывных функций ф0(а), ф1(а), . фт(а) так, чтобы графики соседних функций ф3(а) и ф3+1 (а) в некотором смысле мало отличались друг от друга, функция ф0(а) была мала и фт(а) = ф(а).

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

Их

— = Pox + q^/(r*x), (4)

определив близкое к гармоническому устойчивое периодическое решение х0 (£). Все точки этого устойчивого периодического решения либо расположены в области притяжения устойчивого периодического решения х1 (I) системы

Их

— = Pox + q^(r*x) (5)

с Л = 1, либо при переходе от системы (4) к системе (5) с Л = 1 наблюдается бифуркация потери устойчивости и исчезновения периодического решения. В первом случае можно численно определить х1^), выпуская траекторию системы (5) с Л = 1 из начальной точки х0(0). Стартуя из точки х0 (0), вычислительная процедура после переходного процесса выходит на периодическое решение х1^) и вычисляет его. Для этого промежуток [0, Т], на котором происходит вычисление, должен быть достаточно большим.

После вычисления х1^) можно перейти к следующей системе (5) с Л = 2 и организовать процедуру вычисления периодического решения х2(£), выпуская из начальной точки х2(0) = х1(Т) траекторию, которая при возрастании Ь приближается к периодической траектории х2(£).

Продолжая эту процедуру далее и вычисляя х3 (£), используя траектории системы

(5) с начальными данными х3 (0) = х3-1 (Т), либо приходим к вычислению периодического решения системы (5) с Л = т (т. е. исходной системы (3)), либо на некотором шаге наблюдаем бифуркацию потери устойчивости и исчезновения периодического решения.

Наиболее простым и естественным классом функций ф3 в описанной выше процедуре являются функции ф0(а) = £ф(а), ф1 (а) = £1ф(а),..., фт-1 (а) =

£т-1ф(а), фт(а) = ф(а), где £ — «классический» малый положительный параметр и, например, £3 = Л/т, Л = 1,...,т.

При прямом применении метода гармонического баланса без введения малого параметра и проведения многошаговой процедуры решение с начальными данными, полученными с помощью метода гармонической линеаризации, может стремиться к устойчивому нулевому состоянию равновесия (см. пример 1, рис. 2). Таким образом, мы не всегда можем обнаружить «скрытые» колебания, опираясь лишь на классический метод гармонического баланса. Применение же метода описывающих функций совместно с вышеописанной многошаговой процедурой позволяет произвести численную локализацию «скрытых» колебаний системы. Стартовое значение параметра £ должно быть достаточно малым, чтобы вычислительная процедура смогла выйти на близкое к гармоническому «опорное» периодическое решение. Далее представляется возможным, применяя вышеописанную многошаговую процедуру, постепенно перейти от «опорного» периодического решения к искомому «скрытому» колебанию системы (например, к скрытому аттрактору в обобщенной системе Чуа [9]). При этом выбранный шаг изменения параметра £ должен быть достаточно малым.

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

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

Неособым линейным преобразованием х = 8у приведем систему (4) к виду

Здесь у1, у2 € К, уз € К"-2, А — постоянная (п — 2) х (п — 2)-матрица, все собственные значения которой имеют отрицательные вещественные части, Ь и с — (п — 2)-мерные векторы, &1 и 62 — некоторые вещественные числа.

Не умаляя общности, здесь можно принять, что для матрицы А существует положительное число а, такое что

Здесь В,а1,а2 —некоторые положительные числа. Из липшицевости функции ф(а) и системы вида (6) получим соотношения для решений с начальными данными из 0:

уі = -ш0у2 + Ь\Єф(уі + е*уз), У2 = ^0У1 + Ь-2Єф(уі + е*уз), Уз = Ауз + Ьєф(уі + е*уз).

(6)

Уз(А + А*)уз <-2а|уз|2, V уз € К"2.

Введем следующее множество в фазовом пространстве системы (6):

0 = {|уз|< 0£, у2 =0, у1 € [а1 ,а2}}.

(7)

уі(г) = сов(иог)уі (0) + о(є), У2(Ь) = вт(иоі)уі(0) + О (є), уз(і) = ехр(Аі)уз (0) + Оп-2(є),

(8)

где Оп-2(£) — (п — 2)-мерный вектор с компонентами 0(£). Введем описывающую функцию

2ъ/ ш о

ф(асов(ш01)) сов(^0Ь)д£.

0

Запишем передаточную функцию системы (4)

(9)

и передаточную функцию системы (6)

где I — единичная матрица, ц, в — некоторые вещественные числа, Q(p) —устойчивый полином степени п — 2, Н(р) —полином, степень которого меньше п — 2. Будем полагать, что полиномы Н(р) и Q(p) не имеют общих корней. В силу эквивалентности систем (4) и (6), передаточные функции этих систем совпадают. Отсюда получим соотношения

г1 = -Ъъ в = Ъ2си о, ^4 = с*(А-р1)-1Ь. (И)

Q(P)

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

Теорема 1. Если выполнены условия

ЗФ(а)

Ф(а0) = 0, г}-

За,

> 0,

то при достаточно малых £ > 0 система (4) с передаточной функцией (9) имеет Т-периодическое решение, такое что

2п

г*х(£) = аосо$(и)оЬ) + О(е), Т =--------Ь 0{е).

и0

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

Теорема 2. Если выполнены условия

ЗФ(а)

Ф(а0) = 0, г)-

За

< 0,

то при достаточно малых £ > 0 система (4) с передаточной функцией (9) имеет решение вида

2п

г*х(£) = ао сов(и>о^ + 0(е), 0,—

и0_

и поведение траекторий в окрестности этого решения имеет гиперболический характер.

Теоремы 1 и 2 совпадают с процедурой поиска устойчивых и неустойчивых периодических решений стандартным методом гармонической линеаризации [10].

Из решений вида (8) системы (6) при достаточно малых £ получим значения начальных данных для численной процедуры поиска колебаний:

а0

у(0) = = | 0 I , (12)

где амплитуда а,0 — это корень уравнения Ф(а) = 0.

4. Алгоритм локализации аттракторов для системы Чуа. Системы дифференциальных уравнений, описывающие поведение цепей Чуа [1-4, 6], являются трехмерными динамическими системами с непрерывной скалярной нелинейностью. Здесь мы

а=ао

а=ао

рассмотрим математическую модель, описывающую поведение цепи Чуа с пятью линейными элементами [З, 6], заданную в безразмерных координатах. Эта модель является модификацией широкоизвестной системы Чуа с четырьмя линейными элементами

[5, 11]:

x = а(у — x) — af (x),

i=(x — у + z), (1З)

Z = —(ві + Yz),

где x,y, z Є R, a, e,Y — параметры линейной части системы Чуа,

f(x) = mix + (то — mi)sat(x) = mix + — (то — ті)(|ж + 1| — \х — 1|) (14)

характеризует нелинейный элемент системы (1З), называемый также «диодом Чуа». Заметим, что классическая система Чуа может рассматриваться как система с нелинейностью типа sat(x) («насыщение») [12]. Рассмотрим обобщенную систему Чуа [9], т. е. систему (1З) с нелинейностью вида

ф(х) = mix + ~(то - mi)(|x + 1| - \х - 1|) + ~(s - т0)(|ж + (50| - \х - <5о|)• (15)

Здесь a,e,Y,mo,mi —параметры классической системы Чуа, а So и s —параметры, отвечающие за устойчивость нулевого положения равновесия. Параметр s выбирается так, чтобы линейная часть обобщенной системы Чуа была устойчивой, параметр So определяет область устойчивости нулевого положения равновесия.

На рисунке 1 приведены графики нелинейностей (14) и (15), заштрихованная область — это сектор устойчивости.

x

Рис. 1. Нелинейности стандартной и обобщенной систем Чуа и сектор устойчивости.

Описанная выше процедура позволяет провести численную локализацию аттракторов обобщенной системы Чуа. Для этого введем в систему (13) с нелинейностью (15)

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

Є1 = 0.1 до Є10 = 1.

Теорема 1 позволяет вычислить начальные данные для системы вида (6). Чтобы получить начальные данные для системы вида (4), нужно вычислить матрицу 8, приводящую систему (4) к виду (6).

Запишем систему Чуа в виде системы Лурье

Их

— = Рх + qV'(r*x),

(16)

где

/ —а а 0

х Є М3, Р = I 1 —1 1

\ 0 —в —1

, q

—а

0

0

Далее приведем систему (16) к виду (4):

Их

— = Р0Х + qє^(r*x), оЬ

(17)

где

Ро = Р + kqr*

—а — ка а 0 1 —1 1 0 —в —1

Ар°2 = ±г^о, = А = —3 < 0, ф(а) = Ф(а) — ка.

Далее, систему (17) при помощи неособого линейного преобразования х = 8у приведем к виду (6):

где

^ = Ну + Ьє^-(и*у),

і 0 —ио 0

Н = ио 0 0

! { 0 0 о

Ь

Ьі

1

(18)

1

0

—Н

Выпишем передаточные функции Шр0 (р) системы (17) и Шн(р) системы: (18)

Жр0 (р)= т*(Ро — р1 )-1д =

р2 + (1 + і)р +1 + в

р3 + (а + к +1+ 7 )р2 + (а7 + к7 + к + 7 + в )р + (ав + к7 + кв)’ (19)

Н (—Ьі + Н)р2 + (Ь2^о — ЬіО)р + Ни2 + Ь2^оО

ит < \ ~ьіР + ь2^о ,

№Н(р) = --------2 2 +

р2 + ^2

р + о

р3 + р20 + и2р + и“^0

Приравняв коэффициенты при соответствующих членах числителя и знаменателя в (19), получим значения коэффициентов передаточной функции приведенной системы

г

и

(6), выраженные через ио и параметры линейной части системы Чуа а, в и 7:

к

—аг( + и2 — 7 — в

а(1+ ч)

_ а + ш2 —/3+1+7 + 72 1 + 7 ’

_ а(7 + /3 - (1 + 7)й + И2)

а(7 + /3-^ - (1+7)й) =---------7^~р---------=

о

а

((1 + 7 — °)ио + (7 + в)о)

Ь2 - --------------------г~2

ио(^0 + С0)

Поскольку система (17) приводится к виду (18) неособым линейным преобразованием х = 8у, для матрицы 8 верны следующие соотношения:

Н = 8-1Р08, Ь = 8~^, и* = г*8,

равносильные соотношениям

8Н = Р08, 8Ь = q, и* = г* 8.

Решив эти матричные уравнения, найдем матрицу преобразования

где

«11

I «11 «12 «13

8 = I «21 «22 «23

\ у «31 «32 «33

. = 1, «12 = -0, «13

«22 и>о ' , «23 =

Н(0 — ка — а)

«21 = і + И, «22 =--, «23 = --------,

аа

ка — и2 (ав + кав + ка7 — ^и2) Н(02 — (1 + а + ка)0 + ка)

взі =-------, «32 =--------------------, взз =-----------------------•

а аио а

Используя соотношения (12), запишем соотношения между начальными данными систем (4) и (6):

ао \ / ао«и

х(0) = 8у(0) = 8 ( 0 І = I ао«2і

0 ) \ ао«зі

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

х(0) = а0, 2/(0) = а0(1 + к), г(0) = а0———. (20)

а

=5 0 5 =5 0 5

К X

Рис. 2. Пример 1. Обоснование многошаговой процедуры. Решения обобщенной системы Чуа при £\ =0.1 и £\ = 1 с одинаковыми начальными данными.

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

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

Пример 1. Рассмотрим обобщенную систему Чуа с параметрами в = 0.3, бо = 0.2, т0 = -0.7056, ті = -1.1468, а = 8.4562, в = 12.0732, 7 = 0.0052.

Вычислим значение «стартовой» частоты шо = 2.0392 и коэффициент гармонической линеаризации к = -0.9369. Используя соотношения (20), получаем начальные данные х(0) = 1.2282, у(0) = 0.0774, г(0) = -1.7547 для первого шага многошаговой процедуры построения решений. При є і = 0.1 после переходного процесса вычислительная процедура выходит на близкое к гармоническому периодическое решение. Далее, с увеличением параметра є, полученное близкое к гармоническому периодическое решение трансформируется в аттрактор Чуа [6].

На рисунке 2 показаны траектории обобщенной системы Чуа при є і = 0.1 и є і = 1 с одними и теми же начальными данными, полученными с помощью метода гармонического баланса описанным выше способом. На рисунке 3 изображены графики нелинейностей (14) и (15). Заштрихованная область — это секторы устойчивости [кі, к2] и [к3,к4), где кі = -0.9995, к2 = -0.9370, к3 = -0.1870 и к4 = +го. На рисунках 4, 5, 6 изображены проекции решений на плоскость {х,у} для значений є і = 0.1, єз = 0.3, є5 = 0.5, є7 = 0.7, є9 = 0.9, єі0 = 1.

Пример 2. Рассмотрим обобщенную систему Чуа с параметрами в = 0.2, бо = 0.1, т0 = -0.7056, ті = -1.1468, а = 8.4562, в = 12.62, 7 = 0.0052.

Вычислим значение «стартовой» частоты шо = 2.1892 и коэффициент гармонической линеаризации к = -0.9266. Используя соотношения (20), получим начальные данные х(0) = 1.3251, у(0) = 0.0973, г(0) = -1.9788 для первого шага многошаговой процедуры построения решений. При є і = 0.1 после переходного процесса вычислительная

■5 0 5

X

Рис. 4. Пример 1, £\

с=0.5

'5 0 5

X

Рис. 3. Пример 1. Нелинейности f (ж) и ф(х). Секторы устойчивости [^1,^2] и [кз, +го).

£=0.3

=5 0 5

X

процедура выходит на близкое к гармоническому периодическое решение. Далее, с увеличением параметра є полученное близкое к гармоническому периодическое решение трансформируется в аттрактор Чуа [6].

На рисунке 7 изображены графики нелинейностей (14) и (15). Заштрихованная область— это секторы устойчивости [кі, к2] и [кз,к4), где кі = -0.9995, к2 = -0.9267,

кз = -0.1973 и к4 = +то. На рисунках 8, 9, 10 изображены проекции решений на плоскость {х, у} для значений є і = 0.1, єз = 0.3, є5 = 0.5, є7 = 0.7, єд = 0.9, єіо = 1.

Пример 3. Рассмотрим обобщенную систему Чуа с параметрами в = -0.31, бо = 0.2, то = 0.1691, ті = -0.4768, а = -1.398, в = -0.0136, 7 = -0.0297.

Вычислим значение «стартовой» частоты шо = 0.6436 и коэффициент гармонической линеаризации к = -0.3067. Используя соотношения (20), получим начальные данные х(0) = -1.1061, у(0) = -0.7669, г(0) = 0.0115 для первого шага многошаговой процедуры построения решений. При є і = 0.1 после переходного процесса вычислительная

Рис. 11. Пример 3. Нелинейности /(х) и ф(ж). Сектор устойчивости [к1,к2].

Рис. 12. Пример 3, £1 = 0.1, £3 = 0.3.

Рис. 14. Пример 3, £9 = 0.9, £10 = 1.

Рис. 15. Пример 4. Нелинейности /(ж) и 'ф(х'). Сектор устойчивости [к1,к2].

£=0.3

процедура выходит на близкое к гармоническому периодическое решение. Далее, с увеличением параметра є полученное близкое к гармоническому периодическое решение трансформируется в аттрактор Чуа [6].

На рисунке 11 изображены графики нелинейностей (14) и (15). Заштрихованная область— это сектор устойчивости [к\, к2], где к і = —0.3142 и к2 = —0.3068. На рисунках 12, 13, 14 изображены проекции решений на плоскость {х,у} для значений єі = 0.1, є3 = 0.3, є5 = 0.5, є7 = 0.7, є9 = 0.9, є10 = 1.

Пример 4. Рассмотрим обобщенную систему Чуа с параметрами в = —0.29, бо =

0.1, т0 = 0.44, т1 = —0.4768, а = —1.3018, в = —0.0136, 7 = —0.0297.

Вычислим значение «стартовой» частоты шо = 0.5708 и коэффициент гармонической линеаризации к = —0.2616. Используя соотношения (20), получим начальные данные х(0) = —1.4439, у(0) = —1.0661, г(0) = 0.0164 для первого шага многошаговой про-

цедуры построения решений. При ei =0.1 после переходного процесса вычислительная процедура выходит на близкое к гармоническому периодическое решение. Далее, с увеличением параметра е полученное близкое к гармоническому периодическое решение трансформируется в аттрактор Чуа [6].

На рисунке 15 изображены графики нелинейностей (14) и (15). Заштрихованная область— это сектор устойчивости [ki, k], где ki = —0.3142 и k2 = —0.2617. На рисунках 16, 17, 18 изображены проекции решений на плоскость {x,y} для значений ei = 0.1, е3 = 0.3, е5 = 0.5, е7 = 0.7, е9 = 0.9, е10 = 1.

Литература

1. Chua L. O., Lin G.-N. Canonical realization of Chua’s circuit family // IEEE Transactions on Circuits and Systems. Vol. 37. N 4. July 1990. P. 885-902.

2. Chua L. O. The Genesis of Chua’s circuit // Archiv fur Elektronik und Ubertragung-stechnik. Vol. 46. 1992. P. 250-257.

3. Chua L. O. A Zoo of strange attractors from the canonical Chua’s circuits // Proceedings of the 35th Midwest Symposium on Circuits and Systems (Cat. No. 92CH3099-9), IEEE. Vol. 2. 1992. P. 916-926.

4. Chen G., Ueta T. Chaos in circuits and systems // World Scientific Series on Nonlinear Science. Series B. Vol. 11. 2002.

5. Barboza R., Chua L. O. The Four-element Chua’s circuit // International Journal of Bifurcation and Chaos. Vol. 18. N4. 2008. P. 943-955.

6. Bilotta E., Pantano P. A gallery of Chua attractors // World Scientific Series on Nonlinear Science. Series A. Vol. 61. 2008.

7. Леонов Г. А. О методе гармонической линеаризации // Доклады Академии наук. Сер. Механика. 424(4). 2009. С. 462-464.

8. Леонов Г. А. Эффективные методы поиска периодических колебаний в динамических системах // Прикладная математика и механика. Т. 74. Вып. 1. 2010. С. 37-73.

9. Леонов Г. А., Вагайцев В. И., Кузнецов Н. В. Алгоритм локализации аттракторов Чуа на основе метода гармонической линеаризации // Доклады Академии наук. Сер. Теория управления. 433(3). 2010.

10. Khalil H. K. Nonlinear Systems. New Jersey: Prentice Hall, 2002.

11. Matsumoto T. A Chaotic Attractor from Chua’s Circuit // IEEE Transactions on Circuits and Systems. Vol. CAS-31. N12. December 1984. P. 1055-1058.

12. Llibre J., Ponce E., Ros J. Algebraic Determination of Limit Cycles in a Family of 3-Dimensional Piecewise Linear Differential Systems.

Статья поступила в редакцию 2010 г.