Текст
                    ЭЛЕМЕНТЫ

НЕЛИНЕЙНОЙ

ДИНАМИКИ:

ОТ ПОРЯДКА К ХАОСУ

В. В. Васин, Л. Б. Ряшко ЭЛЕМЕНТЫ НЕЛИНЕЙНОЙ ДИНАМИКИ: ОТ ПОРЯДКА К ХАОСУ Рекомендовано Учебно-методическим советом по математике и механике УМО по классическому университетскому образованию в качестве учебного пособия для студентов физико-математических и технических специальностей /?&с Москва ♦ Ижевск 2006
УДК 517.938 ББК В161.618 В195 Рецензенты: кафедра прикладной математики Уральского государствен- ного технического университета зав. кафедрой докт. физ.- мат. наук, проф. А. Н. Сесекин; докт. физ.-мат. наук, проф. А. Л. Агеев В195 Васин В. В., Ряшко Л. Б. Элементы нелинейной динамики: от порядка к хаосу. — М.-Ижевск: НИЦ «Регулярная и хаотическая динами- ка»; Институт компьютерных исследований, 2006. — 164 с. В пособии излагаются элементы теории хаотического по- ведения (детерминированного хаоса) простейших дискрет- ных и непрерывных динамических систем и обсуждаются основные понятия фрактальной геометрии. ISBN 5-93972-469-8 ББК В161.618 © В. В. Васин, Л. Б. Ряшко, 2006 © НИЦ «Регулярная и хаотическая динамика», 2006 http://rcd.ru http://ics.org.ru
Оглавление Введение ........................................... 7 Глава 1. Дискретные динамические системы 14 § 1. Простая модель — сложная динамика........... 14 §2. Анализ одномерной системы................... 17 2.1. Основные понятия. (17). 2.2. Отыскание точек покоя и циклов. (19). 2.3. Устойчивость инвариант- ных подмножеств. (19). Упражнения (23). § 3. Модели динамики популяции................... 24 3.1. Линейная модель. (24). 3.2. Нелинейная мо- дель. (27). Упражнения (28). §4. Поведение нелинейной системы xt+i = /ixt(l — Xt) . 29 4.1. Случай 0<^^1. (29). 4.2. Случай 1<^3. (30). 4.3. Случай //>3. Бифуркация рождения цикла (30). 4.4. Каскады бифуркаций, удвоение периода и пере- ход к хаосу (32). Упражнения (37). § 5. Хаотическое поведение системы............... 38 5.1. Поведение критических значений параметра. (38). 5.2. Показатель Ляпунова. (41). Упражнения (43). § 6. Универсальность в поведении систем...........43 6.1. Определение константы а. (43). 6.2. Универсаль- ность констант. (44). 6.3. Класс функций, порождаю- щий бифуркацию удвоения. (47). Упражнения (49). § 7. Уравнение удвоения.......................... 50 § 8. Численное нахождение универсальной константы ар 55 §9. Численное нахождение универсальной константы <5Р 58 9.1. Общий метод. (58). 9.2. Прямой алгоритм. (60). Упражнения (62).
Глава 2. Элементы фрактальной геометрии . . 64 § 10. Комплексные динамические системы.......... 64 10.1. Множества Жюлиа, Фату и Мандельброта. (64). 10.2. Классификация множеств Жюлиа. (69). 10.3. N-фуркации системы. (76). Упражнения (78). §11. Отображение подобия и самоподобие множеств ... 79 § 12. Топологическая и фрактальная размерности .... 82 Упражнения (84). § 13. Галерея классических фракталов............ 85 13.1. Канторово множество. (85). 13.2. Кривая Кох (85). 13.3. Фрактал Мандельброта-Гивена (86). 13.4. Решето Серпинского (клиновидная кривая) (86). 13.5. Ковер Серпинского (двумерное множество). (87). 13.6. Фрактал Давида (двумерное множество).(88). 13.7. Пятиугольник Дюрера (двумерное множест- во) (88). 13.8. Губка Серпинского (трехмерное мно- жество) (89). 13.9. Кривая Гильберта. (89). 13.10. Структура аттракторов. (90). Упражнения (91). § 14. Функциональное уравнение для фракталов.....92 14.1. Уравнение для канторова множества. (92). 14.2. Уравнение для кривой Кох. (93). 14.3. Уравнение для решета Серпинского. (94). § 15. Итерационная аппроксимация фракталов.......95 15.1. Свойство сжимаемости оператора Хатчинсона в метрике Хаусдорфа. (95). 15.2. Метод последова- тельных приближений. (98). § 16. Проблема сжатия информации ...............101 16.1. Общий формализм. (101). 16.2. Случай мно- жеств. (103). 16.3. Случай функции. (103). 16.4. Чис- ленный пример. (106). Глава 3. Непрерывные динамические системы 109 § 17. Модели динамических процессов ............109 17.1. Одномерная модель динамики популяции. (109). 17.2. Модель «хищник-жертва». (ПО). 17.3. Линей-
ный осциллятор. (Ш). 17.4. Электронный осцилля- тор. Уравнение Ван-дер-Поля. (ИЗ). 17.5. Химиче- ский осциллятор (брюсселятор). (ИЗ). 17.6. Хао- тический осциллятор. Модель Лоренца. (115). 17.7. Модель Ресслера. (115). § 18. Фазовый портрет системы дифференциальных урав- нений и его свойства ............................116 18.1. Основные понятия. (116). 18.2. Фазовые пор- треты линейных систем. (120). 18.3. Численные ме- тоды решения дифференциальных уравнений. (121). Упражнения (123). §19. Анализ нелинейной системы в окрестности точки покоя............................................126 19.1. Система первого приближения. (126). 19.2. Ус- тойчивость точки покоя. (128). 19.3. Примеры (129). Упражнения (132). §20. Анализ системы в окрестности цикла...........133 20.1. Основные понятия. Система первого прибли- жения. (133). 20.2. Линейные системы с периоди- ческими коэффициентами. Элементы теории Фло- ке. (135). 20.3. Экспоненциальная устойчивость цик- ла.(137). Упражнения (139). §21. Бифуркации...................................140 21.1. Структурная устойчивость и бифуркации. (140). 21.2. Бифуркация рождения цикла. (142). 21.3. По- рядок и хаос в модели Лоренца. (151). Упражне- ния (156). Список литературы...........................158 Предметный указатель .......................162
Я приветствую тех, кто любит хаос. Ведь без хаоса нельзя родить танцу- ющую звезду. Ф. Ницше Введение В бесконечном многообразии явлений окружающего ми- ра человек всегда выделял простое и сложное, упорядочен- ное и непредсказуемое, однозначно определенное и случай- ное. Явления первого типа составляют мир порядка, а все остальное относилось к зоне хаоса. Развитие наук, выяв- ляя неизвестные ранее закономерности явлений, относящих- ся к зоне хаоса, позволяло переводить эти закономерности в зону порядка.
Классическими математическими моделями для про- цессов, наблюдаемых в природе, стали дифференциальные и разностные уравнения. Закон, выраженный таким урав- нением, позволял по известному начальному состоянию ис- следуемой системы однозначно определять ее состояние в любой последующий момент времени. Успехи механики и электродинамики, опирающиеся на уравнения Ньютона - Максвелла, в решении задач описания движения материаль- ных тел и электромагнитных процессов позволяли надеяться сделать предсказуемыми явления в этих важных областях знания. Казалось, дело за немногим: выявить закономерности, лежащие в основе этих, еще мало изученных процессов, за- писать найденные законы в виде дифференциальных или разностных уравнений и «... будущее предстанет перед нами с полной определенностью». Эти слова Лапласа, основопо- ложника идеи детерминизма, рисовали счастливое будущее, когда наконец воцарится желанный порядок, а хаос останет- ся лишь как воспоминание о далеком прошлом, когда циви- лизация делала свои первые шаги. Лапласовские мечты о детерминизме столкнулись с но- вой реальностью. В 60-х годах прошлого века были об- наружены весьма простые по форме записи динамиче- ские модели, имеющие чрезвычайно сложное поведение. Ха- ос появлялся там, где его никак не ждали. Оказалось, что дифференциальные (модель Лоренца) и разностные (модель Ферхюльста) уравнения, задаваемые простейши- ми квадратичными функциями, при сколь угодно малых изменениях их параметров, могут резко менять харак- тер своих решений и от упорядоченных регулярных дви- жений (положение равновесия или периодические колеба- ния) переходить к хаотическим с непредсказуемой дина- микой.
Таким образом, было установлено, что существуют ди- намические системы, поведение которых фактически нель- зя проследить на достаточно широком промежутке времени или, как говорят, система имеет малый горизонт прогнозиро- вания или пределы предсказуемости (например,прогноз по- годы) ввиду чрезвычайно сильной неустойчивости основных характеристик системы. По остроумному замечанию одного из авторов популярной статьи, нелинейная динамика лиши- ла иллюзии глобальной предсказуемости. К концу прошлого века возникла теория динамическо- го хаоса, которая позволяет в некоторых ситуациях описать универсальные сценарии перехода от упорядоченного пове- дения системы к хаосу и наоборот. Было установлено, что хаотическое поведение присуще многим развивающимся си- стемам с нелинейной динамикой и что хаос — достаточно глубокая характеристика природных явлений. Характерным примером проявления этого феномена является любая эко- номическая система, в которой период устойчивого разви- тия, характеризующийся стабильностью базовых макроэко- номических показателей, обязательно сменяется периодом нестабильности, коллапса и хаоса. Это сопровождается, как правило, существенной структурной перестройкой самой си- стемы, после чего система выходит на качественно новый уровень развития. Поэтому можно утверждать, что переход через хаос есть необходимое условие развития любой эконо- мической системы. Под динамической системой (ДС) условим- ся понимать всякую записанную в математических тер- минах совокупность соотношений, которая однозначно оп- ределяет некоторую функцию (в общем случае вектор-фун- кцию), зависящую от дискретного или непрерывного аргу- мента, играющего роль времени, и параметра (скалярного или векторного).
Конечно, это слишком общее определение. Обычно этот термин (т. е. ДС) используют, когда речь идет об объектах или явлениях, описываемых системами уравнений обыкно- венных или в частных производных эволюционного типа. Таким образом, динамическую систему мы определяем в математическом смысле, как правило, отвлекаясь от явле- ний (процессов), которые она описывает. Однако в инженер- ной трактовке (для физиков) обычно различают собственно динамическую систему и ее математическое описание. Приведем для сравнения соответствующее определение из учебного пособия В. С. Анищенко: «Под динамической си- стемой понимают любой объект или процесс, для которого однозначно определено понятие состояния как совокупности некоторых величин в данный момент времени, и задан закон, который описывает изменение (эволюцию) начального состо- яния с течением времени. Динамические системы — это ме- ханические, физические, химические и биологические объек- ты, вычислительные процессы преобразования информации, совершаемые в соответствии с конкретным алгоритмом» [2]. Значительный всплеск интереса исследователей к хао- тическим ДС возник после выхода работ М. Фейгенбаума в конце 70-х годов прошлого столетия, который, исследуя простейший итерационный процесс с квадратичной нелиней- ностью, обнаружил хаотическое поведение итераций и нашел две универсальные константы, характеризующие поведение целого класса ДС, включая систему Лоренца [19]. Примерно в этот же период произошло еще одно со- бытие, вызвавшее большой резонанс, это открытие Б. Ма- ндельброта, который, изучая с помощью компьютерного моделирования итерации с комплексной квадратичной фу- нкцией (т. е. двумерный итерационный процесс), обнаружил, что область изменения комплексного параметра (множество Мандельброта), при котором траектория критической точки
системы удерживается в ограниченной области, имеет уди- вительно причудливую форму, состоящую из самоподобных фрагментов (см. об этом [29]). Самоподобную и нерегулярную структуру образуют также множества Жюлиа (первые компьютерные изображе- ния этих множеств получены Дж. X. Хаббардом и Б. Ман- дельбротом), которые порождаются ограниченными траек- ториями комплексной квадратичной функции с фиксирован- ным параметром из множества Мандельброта. Множества, обладающие свойством подобия и некото- рой «геометрической хаотичностью», были названы Б. Ма- ндельбротом фракталами. Таким образом, возникла новая дисциплина — фрактальная геометрия. Было обнаружено, что аттракторы (предельные точки квадратичного итерационного процесса) в бесконечной це- пи бифуркаций (удвоение периода) образуют самоподобные, т. е. фрактальные, структуры с величиной скейлинга, равной 2-й универсальной константе Фейгенбаума. Также выяснилось, что еще в 20-х годах, в работах французских математиков Г. Жюлиа и П. Фату были пред- сказаны и изучены многие свойства комплексных ДС. Но эти работы не были оценены, даже были забыты на многие деся- тилетия. Отчасти это произошло, может быть, потому, что из-за отсутствия компьютеров в то время не было возможно- сти проиллюстрировать их результаты наглядно в графиче- ском виде и показать завораживающую красоту фракталь- ных форм. Причем, что интересно отметить, фракталы, рожденные как математические объекты в процессе исследования очень простых динамических систем, удивительным образом напо- минают то причудливую форму облака, структуру кристал- ла или фрагмент Галактики, то сильно изрезанную берего- вую линию, трещину в металле или ювелирное украшение.
В пособии излагаются элементы теории хаотического по- ведения на примере простейших динамических систем, в ко- торых переход к хаосу сопровождается бифуркацией (п-фур- кацией) периода системы. Исследуются два класса динами- ческих систем: дискретные ДС, задаваемые одномерным или двумерным итерационным процессом, и непрерывные ДС, описываемые нелинейными дифференциальными уравнени- ями специального вида и обладающие сильной чувствитель- ностью к изменениям входящих параметров. Учебное пособие состоит из трех глав, разделенных на параграфы, и организовано следующим образом. Глава 1 целиком посвящена дискретным ДС, главным образом, итерационному процессу с действительной квад- ратичной функцией шага. Рассматривается переход к хао- су через цепочку бифуркаций удвоения периода. Вводятся базовые понятия, определяются универсальные константы, устанавливается связь этих констант с типом нелинейности итерируемой функции, предпринимается попытка описания класса отображений, порождающих явления удвоения пери- ода. Глава 2 — это введение во фрактальную геометрию. Выясняется поведение итераций с комплексной квадратич- ной функцией (двумерный итерационный процесс), описыва- ются основные классы множеств Жюлиа и их связь с мно- жеством Мандельброта. Далее, определяется фрактальная размерность, позволяющая охарактеризовать класс фрак- тальных множеств, описываются классические фракталы и выводятся уравнения, которым они удовлетворяют, из- лагаются методы итерационной аппроксимации фракталов и подходы к проблеме сжатия информации. Глава 3 посвящена непрерывным ДС. Рассматривают- ся классические модели, задаваемые дифференциальными уравнениями: модель хищник - жертва, механическая, элек-
тронная (уравнение Ван-дер-Поля) и химическая (брюсселя- тор) колебательные системы, хаотические осцилляторы (мо- дель Лоренца и Ресслера). Дается способ описания динамики систем при помощи фазовых портретов. Излагаются методы анализа ДС вблизи точки покоя и цикла. Детально рассмат- ривается бифуркация рождения цикла. На примере циклов модели Лоренца представлен каскад бифуркаций удвоения периода, приводящий к хаосу. Содержание пособия основано на материале, который читался авторами в течение ряда лет в Уральском государ- ственном университете. Оно может быть рекомендовано сту- дентам физико-математических и технических специально- стей для первоначального знакомства с хаотической дина- микой и фрактальной геометрией.
Глава 1 Дискретные динамические системы В последние годы стало ясно, что высо- кая чувствительность к начальным услови- ям, приводящая к хаотическому поведению во времени, никоим образом не исключение, — это типичное свойство многих систем. Г. Шустер You can’t know how happy I am that we met I’m strangely attracted to you. C. Porter § 1. Простая модель — сложная динамика Многие считают, что сложное поведение динамической системы определяется большим количеством переменных, связанных громоздкими формулами. Однако существуют примеры очень простых систем, демонстрирующих доста- точно сложную динамику. Одним из таких примеров явля- ется рассматриваемая ниже система, в поведении которой присутствуют основные черты таких важных явлений нели- нейной динамики, как циклы различной кратности, переме- шивание, хаос.
Рассмотрим динамику системы, заданной итерацион- ным процессом xt+1 = {10жД, £ = 0,1,... (1.1) Функция шага <^(ж) = {Юж}, где {у} означает дробную часть числа у, по заданному начальному значению Жо однозначно определяет последовательность Ж1 = у?(жо),Ж2 = <£>(Ж1),... . Отображение <£>(ж) переводит полуинтервал [0,1) в себя. Это означает, что, начиная с произвольной точки жо € [0,1), все последующие элементы Ж1,ж2,... будут также лежать в полуинтервале [0,1). Примеры последовательностей ж0 = 0.27, Ж1 = 0.7, Ж2 = 0, ... хп = 0, ... Жо = 0.51(3), Ж1 = 0.1(3), ж2 = 0.(3), ... жп = 0.(3), ... показывают, что для некоторых начальных значений состоя- ние системы уже через несколько шагов перестает изменять- ся, т. е. последовательность {ж4} становится стационарной. При этом £0 = 0,= 0.(1),...,^8 = 0.(8) являются непо- движными точками отображения </?(ж) или точками покоя системы (1.1). Взяв другие начальные значения, можно получить по- следовательности ж0 = 0.27(35), ж3 = 0.(53), ж0 = 0.8(307), ж3 = 0.(730), Ж1 = 0.7(35), ж4 = 0.(35), Ж! = 0.(307), Ж4 = 0.(307), ж2 = 0.(35), ж5 = 0.(53), ж2 = 0.(073), ж5 = 0.(073), (1-2) (1-3) элементы которых через некоторое число шагов начинают повторяться, образуя цикл. Для последовательности (1.2)
цикл состоит из двух элементов (2-цикл) £i = 0.(35), £2 = = 0.(53). Для последовательности (1.3) наблюдается 3-цикл £1 = 0.(307),£2 = 0.(073),£3 = 0.(730). С какого бы рацио- нального числа Xq = O.oi • • • an(6i • • • Ьт) мы ни стартовали, итерационный процесс (1.1) через конечное число шагов ста- новится периодическим: Х1 = О.О2 . . . ап(61 . . . 6т)> • • • , Хп = 0.(61 • • • ^т)> п+1 = О.(&2 • • • ЬтЬ\), ... , Хп+т = 0.(61 . . . 6т), .... При этом длина переходного процесса равна п, а длина цик- ла определяется количеством знаков т в повторяющейся ча- сти десятичного позиционного представления начального со- стояния xq. Как видим, поведение системы является вполне упорядоченным и предсказуемым. В этом случае говорят, что в системе наблюдается порядок. Последовательность, стартующая с иррационального числа, ведет себя гораздо сложнее. Например, взяв Xq = = {л} = 0.14159265..., получим последовательность, кото- рая бесконечное число раз попадает в сколь угодно малую окрестность каждой точки полуинтервала [0,1). Поведение последовательности выглядит как случайное или хаотиче- ское. Тогда говорят, что в системе наблюдается хаос. Рассмотрим задачу приближенного описания динамики системы (1.1), выходящей из точки xq — {л}. Десятичная дробь, представляющая а?о, содержит бесконечную последо- вательность цифр после запятой. В качестве приближения для Xq возьмем х*0 = 0.14159265. Тогда найденные по форму- ле = y{x*t} последующие элементы х^,х%,... являются соответствующими приближениями для Xi, Х2, •. • • При этом ошибка Д( = Xf-x* меняется следующим обра- зом: 10-7 < До < ю~8, 10“6 < Д1 < 10"7,..., 10"1 < Д8 < 1. Десятикратный рост ошибки на каждом шаге приводит к то- му, что уже через восемь шагов мы ничего не можем сказать
о поведении полученного решения. Действительно, элемен- ты приближения х*8 = 0 = х$ = х}0... уже не несут никакой информации об истинных значениях Х8,Хд,Хю, • • • • Отображение промежутка [0,1) в себя, задаваемое фун- кцией </?(х) = {Юх}, можно разделить на этапы. Разобьем весь промежуток [0,1) на 10 частей: Д = [0,0.1), I2 = = [0.1,0.2),...,/ю = [0.9,1.0). Сначала функция ipi(х) = = Юх растягивает каждый из полуинтервалов 1п в десять раз: /[ = </?i[/i] = [0,1.0),/2 = уд[/2] = [1,2), l'lo = V’lf/io] = = [9,10). Затем функция </?2(х) = {х} переводит каждый из этих растянутых полуинтервалов в один исходный по- луинтервал [0,1), отождествляя (склеивая) соответствую- щие точки. В результате функция <^(х) = ^2[9?i (х)] есть суперпозиция двух функций — растяжения <pi (х) и склеива- ния <^2(х). Последовательные итерации растяжения и скле- ивания, удаляя друг от друга близкие точки и сближая да- лекие, хорошо перемешивают точки полуинтервала [0,1). Замечание. Если считать хо случайной величиной, распре- деленной на интервале [0,1) с плотностью ро(х), то щ(х) — плот- ность распределения xi = <р(хо) — получается из Ро(х) усредне- нием 10 Pi (ж) = у^52Ро(Ю(х - ^)). г=1 § 2. Анализ одномерной системы 2.1. Основные понятия. Рассмотрим одномерную систему, определяемую итерационным процессом xt+i = 9?(xt), (2.1) где х — скалярная переменная, X — область определения, a Y = <р(Х) — область значений функции <р(х). Предпола-
гается, что Y Q X. Тогда для любого xq € X процесс (2.1) задает последовательность xt (t = 0,1,...), которая называ- ется орбитой ТОЧКИ Xq. Определение 2.1. Множество М С X называется инвариантом системы ( 2 . 1), если <р(М) С М. Если xq € М, то и все последующие элементы xt G М. Простейшим примером инвариантного множества является точка покоя, т. е. неподвижная точка отображения <р. Определение 2.2. Точка С € X называется точкой покоя системы (2.1), если £ = 92(C). ЕСЛИ Xq = С, ТО Xq = Хг = Xq = . . . = Xt = Xt+i = . . ., все последующие элементы не меняются. Важным примером инвариантного множества является цикл. Определение 2.3. Множество М = {Ci, С2, • • •, Cfc} на- зывается к - циклом системы (2.1), если между его эле- ментами имеется следующая связь: €2 = <XCi), Сз = 9?(Сг), • • • , С* = y’(Cfe-i), €1 = <£(&)• Определение 2.4. Точка С называется Апериоди- ческой, если С = 9?fc(C) = ¥’(9?(- - • (^(С)))- Здесь функция <р применяется к раз. Пусть у системы (2.1) имеется 2-цикл М = {Ci, €2}- То- гда ДЛЯ Xq = Cl получим ПОСЛеДОВатеЛЬНОСТЬ Xi = С2,^2 = = Сь • • •, x2t = Ci^2t+1 = С2, - • -, а для Xq = С2 — последова- тельность Xi = Cl, х2 = &,---ix2t = C2,^21+l = Cl,---- Оба эти решения являются 2-периодическими: xt+2 = xt. В общем случае каждая последовательность, стартую- щая из произвольного элемента Ацикла х0 = является Апериодической: xt+k = xt.
2.2. Отыскание точек покоя и циклов. Точка по- коя £ системы (2.1) удовлетворяет равенству £ — р(£). Таким образом, для отыскания всех точек покоя системы (2.1) тре- буется найти все решения уравнения х = </>(т). Поиск цикла также сводится к решению некоторого уравнения. Действительно, рассмотрим fc-кратную суперпо- зицию функций </?(х) 9?(т) = ^(...<^(х))) и связанную с <рк(х) динамическую систему »+1 = Лй)- (2-2). При одинаковых начальных значениях у$ — Xq элементы ор- бит для систем (2.1) и (2.2) связаны соотношением yt = Xkt- Как видим, {yt} является подпоследовательностью для по- следовательности xt. При этом элементы yt получаются вы- делением (сечением) в последовательности {xt} элементов, отстоящих друг от друга на к шагов. Каждый элемент & А:-цикла {£i, удовлетво- ряя соотношению £, = </?fe(^), является точкой покоя систе- мы (2.2). При этом все элементы £i,...,& являются корнями уравнения х = <рк(х). Среди корней этого уравнения могут быть не только элементы /с-цикла. Действительно, очевид- ным корнем этого уравнения является точка покоя £ систе- мы (2.1). Если система (2.1) имеет более короткий /-цикл, где I кратно к и I < к, то все элементы этого Z-цикла также являются корнями уравнения х — <рк(х). 2.3. Устойчивость инвариантных подмножеств. Пусть р(хо, М) = inf |то — х\ — расстояние от точки Хо до хем подмножества М.
Определение 2.5. Инвариантное подмножество М в системе (2.1) называется устойчивым по Ляпуно- ву, если в некоторой его окрестности U (М С U) справед- ливо следующее: Ve > 0 35 > 0: Ух0 е U V* р(х0, М) < 5 => p(xt, М) < е. В противном случае подмножество М называется не- устойчивым. Определение 2.6. Устойчивое по Ляпунову подмно- жество М называется асимптотически устойчи- вы м, если Vzq G U lim p(xt,M) = 0, где U — некоторая t—>оо окрестность множества М. Для простейшего инвариантного множества, состоящего из одной точки покоя С, достаточное условие асимптотиче- ской устойчивости дается следующей теоремой. Теорема 2.1. Пусть в некоторой окрестности Ue = = (€ — точки покоя £ системы (2.1) существу- ет непрерывная производная !р'(х) и выполняется неравен- ство |^(€)| < 1- Тогда точка покоя £ является асимпто- тически устойчивой. Доказательство. Благодаря непрерывности р'(х), в некоторой окрестно- сти Ug = (£—5, £+5) выполняется неравенство |<p'(z)| < 1- При этом, используя формулу Лагранжа, получим |</?(z) — — £| = W)||z - е| — £|- Данное неравенство означает, что при любом Хо G Us справедлива оценка |zt — £| Q11жо — — С| < 5, гарантирующая асимптотическую устойчивость точки покоя £. Теорема 2.2. Пусть |<^'(£)| > 1- Тогда С — неустойчи- вая точка покоя.
Отметим, что величина q = является количе- ственной характеристикой степени устойчивости точки по- коя £. Действительно, вблизи £ справедливо приближенное равенство l^t+i -€1« W)lkt-€l- Величина q показывает, во сколько раз решение систе- мы (2.1) приближается к точке покоя £ за один шаг. При q = 1 (критический случай) точка может быть как устойчивой, так и неустойчивой. Исследование устойчивости цикла систе- мы (2.1), как следует из приведенных выше рассуждений, сводится к анализу устойчивости его элементов £i,...,£fc как точек покоя системы (2.2). В соответствии с теорема- ми 2.1, 2.2 все решает производная (</?fc(x)) в этих точках. Теорема 2.3. Пусть Xq, х^, ..., Xk-i — последователь- ные итерационные точки в процессе (2.1). Тогда справедли- ва формула k—1 (у‘Ы)' = ПЛ)' (2.3) t=0 Доказательство. Применим индукцию по к. Для к = 0 формула очевидна. Предположим, что она справедлива для /с—1, и докажем, что она имеет место для к. Действительно, имеем (</(х0))' = [^(<pfc-1(x0))]' = (/((/_1(хо)) • (</-1(яо))' = к-2 к-1 = 9?'(arfc_i) • JJyCzt) = J]V(zt). t=0 t=0 Следствие 2.1. Пусть {£i,. ..,£&} — к-цикл. Тогда производная (<^fc(x)) во всех точках & (г = 1,2,..., к) име-
ет одно и то же значение k «= (Ле.))'= Пу «-) t=0 Теорема 2.4. Пусть в некоторой окрестности к-цик- ла системы (2.1) существует непрерывная производ- ная <р'(х) и выполняется неравенство k пу&) i=l (2-4) Тогда к-цикл является асимптотически устойчивым. Доказательство. Из выражений (2.3), (2.4) вытекают неравенства |(</’/с(&))/| < 1 (г = 1,..., к), из которых по теореме 2.1 сле- дует асимптотическая устойчивость всех точек покоя & для системы (2.2). Рассмотрим точку Ее асимптотиче- ская устойчивость означает существование окрестности U = = (€i ~ <5)€i + <^)> для которой соответствующие решения yt системы (2.2) обладают следующими свойствами: Ууо Vi уо € U => yt G U, lim yt = t—>oo Свяжем с последовательностью yt последователь- ность Xj, порождаемую системой (2.1) с начальным усло- вием х0 = у0. Отметим, что yt = Xkt. Теперь для каждого i G {1,..., к} в последовательности Xj выделим подпоследо- вательность yl — Xkt+i-i- Отметим, что последовательность ylt есть решение системы (2.2) с начальным условием уг0 = х^, при этом = (2.5)
Из непрерывности ip'(х) следует существование констан- ты К >0, для которой справедлива оценка ИД») (2.6) Из (2.5),(2.6) следует неравенство 1?4 -61 K\yt -б|. Теперь сходимость yt к £i влечет для каждого i сходимость у\ к £г, при t —* оо. Таким образом, последовательность Xj схо- дится к циклу {6> •••>€*;}• Асимптотическая устойчивость цикла доказана. Упражнения 2.1. Для системы zt+1 = zf + 0.1 а) доказать существование трех точек покоя 6 <6<6; б) найти непересекающиеся интервалы, содержащие эти точки покоя; в) доказать, что из этих трех точек только 6 является устойчивой; г) найти максимальный инвариантный интервал, содер- жащий д)доказать, что, начиная с любой точки этого интерва- ла, последовательность xt стремится к £2- 2.2. Доказать теорему 2.2. 2.3. Для критического случая, когда в точке покоя £ вы- полняется равенство |<^'(^)| = 1, привести примеры устойчи- вости и неустойчивости. 2.4. Указать условия, при которых последовательность xt системы (2.1) стремится к устойчивой точке покоя £ моно- тонно.
2.5. Для системы т<+1 = p,xt + 7 указать значения пара- метров, при которых система имеет 2-цикл. Будет ли этот цикл устойчивым (асимптотически устойчивым)? § 3. Модели динамики популяции Сообщество животных одного вида, населяющих опре- деленную территорию, называют популяцией. Изучать жизнь отдельной популяции можно с самых различных то- чек зрения — разнообразие природы бесконечно. Нас будет интересовать лишь численность популяции Nt в различные моменты времени t = 1,2,... и ее динамика. 3.1. Линейная модель. Для построения модели тре- буется учесть основные факторы, влияющие на изменение численности. Таковыми являются рождаемость и смерт- ность. Предполагается, что за время, прошедшее между со- седними моментами наблюдений t и t + 1, появилось на свет aNt и умерло (3Nt особей. Данное предположение кажется достаточно естественным: количество как родившихся, так и умерших особей должно быть пропорционально общему числу Nt особей популяции. Соответствующие пропорции за- даются параметрами: а — коэффициент рождаемости и /3 — коэффициент смертности. В результате учета этих факторов получаем линейную модель динамики популяции Nt+i = Nt + cxNt — 0Nt = (iNt. (3.1) Эта модель в действительности зависит только от одного па- раметра /1 = 1 + а — (3 — коэффициента естественного при- роста. Используя уравнение (3.1), легко предсказывать значе- ния численности популяции. Действительно, уравнение (3.1)
означает, что последовательность есть геометри- ческая прогрессия со знаменателем д, каждый элемент ко- торой Nt может быть выражен через начальный No соотно- шением Nt = у1 Nq- Что же ожидает популяцию в будущем? В рамках модели (3.1) в зависимости от величины д для данной популяции возможны три варианта динамики: 1) при 0 < д < 1 численность популяции монотонно убывает к нулю. В условиях превышения смертности над рождаемостью (а < /3) популяция вымирает; 2) при д = 1 численность популяции не изменяется: Nt = No. Популяция, благодаря балансу между рождаемо- стью и смертностью (а = /3), находится в состоянии равно- весия; 3) при д > 1 численность популяции монотонно возрас- тает к бесконечности. Превышение рождаемости над смерт- ностью (а > /3) ведет к демографическому взрыву. Наглядное сравнение динамики численности популяции для этих трех случаев дано на рис. 3.1. Для первого и третьего случаев на рис. 3.2 и рис. 3.3 с помощью графиков функции у = <д(х) = дх иллюстриру- ется геометрический метод получения по начальному значе- нию No последующих элементов М, N2, -... Для этого следу- ет из точки No на оси ОХ сначала провести вертикаль до гра- фика ц>(х), а затем — горизонталь до графика у = х. Абсцис- са найденной точки и есть N\. Далее вертикаль проводим уже из точки N\ и т. д. В результате проделанных построе- ний последовательность Nt изображается в виде «лестницы Ламерея». Направленное перемещение по этой «лестнице» дает возможность наглядно судить о характере поведения последовательности Nt. Как видим, величина Nq численности популяции в на- чальный момент времени влияет лишь на количественную
0123456789 t Рис. 3.2. Динамика при 0< д< 1 Рис. 3.1. Линейная модель. Три варианта динамики сторону дела. Качественная картина динамики определя- ется исключительно параметром /л. Значение параметра д = д* называют критическим, или бифуркационным (от лат. bifurcus — раздвоенный), если при переходе д че-
рез д* характер динамики системы качественно изменяется. В данной модели критическим значением, отделяющим один вариант от другого, является д* = 1. 3.2. Нелинейная модель. Неограниченный рост численности популяции, следуемый из модели (3.1) при д >1 (см. вариант 3), в реальности не наблюдается. Ограничен- ность жизненного пространства и, прежде всего, недоста- ток продуктов питания приведут к неизбежному замедлению роста численности. Территория, на которой обитает данная популяция, в состоянии прокормить лишь определенное ко- личество особей. По мере заполнения экологической ниши внутри популяции нарастает напряженность: жизненных ре- сурсов на всех не хватает. На динамику численности все сильнее начинает воздействовать новый фактор — голод. Для учета этого фактора П. Ф. Ферхюльст еще в 1845 го- ду предложил добавить в уравнение (3.1) нелинейный член. Новая модель, называемая в биологии логистическим урав- нением, имеет вид М+1 = дМ - ~(N2t. (3.2) Здесь 7 — коэффициент смертности, связанной с ограничен- ностью ресурса. Выбор квадратичной зависимости в (3.2) можно пояснить следующим образом. Недостаток ресурсов порождает внутривидовую борь- бу, интенсивность которой пропорциональна количеству воз- можных контактов между отдельными особями. В популя- ции из N особей количество возможных парных (другие учи- тывать не будем) контактов пропорционально N2. Не все контакты заканчиваются летальным исходом. Соответству- ющий процент и задается коэффициентом 7. Займемся анализом этой нелинейной модели. Сначала упростим ее заменой xt = -^Nt, перейдя от (3.2) к системе,
зависящей уже только от одного параметра: Tt+i = <^(irt), = рх(1 - х). (3.3) Как видим, именно д — коэффициент естественного приро- ста — является тем параметром, который определяет каче- ственную картину динамики данной модели. При этом пара- метр 7 играет роль масштабирующего множителя. Модель (3.3) сохраняет биологический смысл (числен- ность популяции не может быть отрицательным числом) лишь при ц > 0 и xt G [0,1]. Причем условие xt G [0,1] гарантируется лишь при д 4, поскольку максимум функ- ции <р(х) на этом отрезке равен д/4. В строгой математиче- ской формулировке это выглядит так: при любом р € [0,4] функция <д(ж) задает отображение отрезка [0,1] в себя. Та- ким образом, при 0 р С 4, стартуя из любой точки Хо, лежащей на отрезке [0,1], все последующие элементы Xi, Х2,..., xt,..., формируемые системой (3.3), будут принад- лежать этому же отрезку. Отрезок [0,1] есть инвариантное множество отображения <д(х). Упражнения 3.1. Найти точки покоя и исследовать их устойчивость для линейной модели динамики популяции с внешним воз- действием и: М+1 = Р-Nt 4- щ а) при 0<д<1им4 = 7>0, 7 — величина постоянного внешнего притока особей; б) при д > 1 и и* = —7 < 0, 7 — величина постоянного оттока особей; в) при р > 1 и щ = —7 — k(Nt — N), 7 > 0. Здесь внешнее воздействие строится по принципу обратной связи
по отклонению Nt — N состояния системы Nt от желаемого уровня численности N, к — коэффициент обратной связи. 3.2. Для модели Риккера Nt+1 = /j,Ntexp(-Nt) найти точки покоя и исследовать их устойчивость в зависи- мости от параметра ц > 0. § 4. Поведение нелинейной системы xt+i = Mxt(l - xt) Рассмотрим динамику системы (3.3) при различных зна- чениях параметра ц, двигаясь по отрезку [0,4] от левого кон- ца к правому. Найдем сначала точки покоя. Корнями уравнения х = = <Дт) являются = 0, £2 = 1 — д- Отметим, что в интере- сующий нас отрезок [0,1] точка £2 попадает лишь при ц 1. 4.1. Случай 0 < у, 1. В данном диапазоне изме- нения параметра ц система имеет единственную точку по- коя & = 0. Анализ устойчивости £i связан (см. теорему 1) с величиной q = |<Д(£1)| = ц. При 0 < /1 1 точка покоя яв- ляется асимптотически устойчивой. Для любого То из [0,1] последовательность xt, монотонно убывая, сходится к £i как геометрическая прогрессия со знаменателем, близким к д. Для /д близких к нулю, скорость сходимости высокая. При приближении ц к единице скорость сходимости падает. Слу- чай /л = 1 не охватывается теоремой 1 и является критиче- ским. Однако и здесь последовательность xt монотонно схо- дится к точке £1 = 0 (см. рис. 4.1). При ц > 1 точка покоя £i = 0 становится (см. теорему 2) неустойчивой. Для исследования устойчивости точки покоя
Рис. 4.1. Динамика системы при ц = 1 £2 = 1 — д рассмотрим производную у/(£2) = 2 — д. Неравен- ство q = |2 — д| < 1 имеет решение 1 < /л < 3, что приводит к следующему случаю. 4.2. Случай 1 < р, 3. В этом случае у систе- мы (3.3) имеется две точки покоя: £1 — неустойчивая, £2- устойчивая. Для любого Xq из (0,1) последовательность xt сходится к £2. Скорость сходимости определяется величи- ной q = |2 — д| < 1. В окрестности точки ^2, при 1 < д 2 последователь- ность сходится монотонно (рис. 4.2), а при 2 < д 3 моно- тонность нарушается (рис. 4.3). При переходе д через значение /л^ — 3 точка покоя £2 становится неустойчивой. Что же происходит в системе при когда обе точки покоя являются неустойчивыми? 4.3. Случай (л > 3. Бифуркация рождения цик- ла. Рассмотрим отображение xt —* xt+2 за два шага, зада- ваемое функцией <^2(ж) — <^(<^(ж)) = ц2ж(1 — х)(1 — щг(1 — ж)), х1+г = (4.1) Анализ системы (4.1) начнем, как обычно, с отыска- ния точек покоя. Уравнение х = <р2(х) — алгебраическое
Рис. 4.2. Динамика системы при 1 < ц 2 Рис. 4.3. Динамика системы при 2 < д 3 уравнение четвертой степени — имеет своими корнями ра- нее найденные £i и Это и понятно: точки покоя £1,^2 си- стемы (3.3) остаются точками покоя также и системы (4.1). Понижая степень на две единицы, получим квадратное урав- нение /лх2 - (д + l)z + (1 + д) = О с корнями 1Л + 1 ± ^/(д+1)(д-3) &’4 = 2д В интересующем нас здесь случае д > 3 корни ^3,^4 — вещественны и различны. Каков же смысл у этой пары с точки зрения исходной модели (3.3)? Пусть хо = £3. Най- дем далее ад = у>(Сз)- Из равенств (/?2(ад) = </?(<^(xi)) = — <^(<д(<д(£з))) = ^(€з) = Ti следует, что ад является точкой покоя системы (4.1). Поскольку Xi не совпадает ни с одной из точек £1,£2>£з, то Xi — £4. Далее находим xz = = = Т?2(€з) = €з,жз = £4. Формируемая здесь последователь-
ность получается периодической: x^t = Сз> x2t+i = €4- Ее эле- менты составляют цикл S2 = {£3, £4} (рис. 4.4). Таким образом, при переходе параметра ц через кри- тическое значение /zi = 3 в системе (3.3) происходит би- фуркация: одновременно с потерей устойчивости точкой покоя £2 рождается цикл S2 периода 2. Эта бифуркация со- провождается рождением двух новых точек покоя £3,^4 си- стемы (4.1). При этом устойчивость этих точек покоя в си- стеме (4.1) эквивалентна устойчивости цикла S2 в систе- ме (3.3). 4.4. Каскады бифуркаций, удвоение периода и переход к хаосу. При дальнейшем увеличении д точ- ки покоя Сз и £4 сохраняют устойчивость лишь до некото- рого следующего бифуркационного значения д2. Переход д через fj,2 ведет к одновременной потере устойчивости у то- чек £з и £4 . При этом у отображения <д4(ж) = <д2(<д2(х)) = = <£>(<£>(<£>(<д(я)))) — многочлена 8-й степени — рождаются че- тыре новые устойчивые точки покоя ^5,^6, ^7,^8-
Данная бифуркация с точки зрения исходной системы означает рождение устойчивого цикла S4 — {£5, £7, £8} : £б = <£(&),& = ^(£б),& = <£>(&),€5 = ^(&)- Численность популяции повторяется через четыре шага (см. рис. 4.5). Дальнейшее увеличение д обнаруживает аналогичные бифуркационные значения дз, д4, д5..связанные с рожде- нием циклов S8, S™, S32 .... При этом каждый раз, прохо- дя очередное бифуркационное значение, соответствующий цикл S2k теряет устойчивость, происходит бифуркация уд- воения периода, и рождается устойчивый цикл S4k. Как видим, появление достаточно разнообразных цик- лов (ритмов) [8] в жизни популяции можно объяснить сугубо внутренними факторами, не изобретая в качестве первопри- чины каких-то периодических внешних воздействий. Общая картина усложнения циклов, происходящих в ре- зультате бифуркаций удвоения периода, представлена на рис. 4.6. Здесь при каждом д указано устойчивое ин- вариантное предельное множество (множество предельных точек) — аттрактор (attract — притягивать) систе- мы (3.3). По мере прохождения бифуркационных значений д0, Д1,... картина усложняется. На интервале [0, до] аттракто- ром является устойчивая точка покоя £Дд) = 0. На интер- вале (до,А<1]5 где £1 теряет устойчивость и далее не изобра- жается, расположен график £2(д). На следующем интерва- ле (Д1,Д2], где £2 теряет устойчивость, появляются графи- ки £з(д) и £4(д) — точки аттрактора S2. Далее, на (д3,д4] представлены графики £«(д) для i = 5,6,7,8 — точки ат- трактора S4 и т. д. Отметим своеобразие нелинейных систем. В линейном случае появление неустойчивой точки всегда сопровожда- ется уходом решений в бесконечность (система разруша- ется).
Рис. 4.6. Бифуркационная диаграмма В нелинейном — потеря устойчивости ведет к появлению новых качественных особенностей в ее поведении, рождению новых аттракторов. Последовательность цп при п —> оо имеет предел (обо- значим его /Zoo). Система (3.3) при /л, = /л^ — 3.5699456... формирует сложную непериодическую последовательность — хаос. Соответствующие предельные множества получили на- звание странных аттракторов. Некоторое пред-
ставление о характере поведения системы (3.3) в зоне (Доо,4] можно получить из рис. 4.7. Его левая часть повторяет би- фуркационную схему рис. 4.6 и отражает зону параметров, для которых система ведет себя периодически ((0, Доо) — зо- на «порядка»). Здесь аттракторы состоят из конечного на- Рис. 4.7. Самоподобие
бора точек. При переходе параметра р, через р1ГХ. «порядок» сменяется «хаосом». Аттракторы начинают выглядеть как множества, состоящие из сплошных (зачерненных) интерва- лов. В этой зоне выделяются и просветы — окна, в которых снова виден «порядок». В этих окнах расположены аттрак- торы, отвечающие периодическим циклам вида Sp 2 , Р = = 3, 5, 7,.... Наиболее отчетливо в самом широком окне ви- ден цикл S3. В зоне (моо,4] содержится бесконечное количество би- фуркаций «хаос» —> «порядок» и «порядок» —> «хаос». Сле- дует подчеркнуть, что в режиме «хаоса» при больших п практически невозможно предсказать значение хп. Как же так, может спросить читатель, ведь состояние xt однозначно определяется из (3.3) по t и Xq? Дело в том, что неизбежные, пусть даже очень малые, ошибки в определении начального значения Xq приведут к тому, что вместо последовательно- сти Xq, Xi, х2, ..., xt,... вы будете получать другую последо- вательность — Xq, х\, Х%, • • • Xf ,.... В режиме «хаоса», стартуя с практически неотличимых значений xq и Xq, уже через несколько итераций элементы х^ потеряют всякую связь с истинными значениями xt. В «хаотическом» режиме элементы последовательно- сти хп начинают вести себя как случайные величины. Для описания таких последовательностей используют при- емы, принятые в теории вероятности и математической ста- тистике. В бифуркационной картине на рис. 4.7 а, б, в вы- делена последовательность кадров а, б, в, вложенных один в другой. При этом каждый последующий кадр представля- ет собой увеличенный фрагмент предыдущего. Как видим, все эти кадры удивительным образом воспроизводят прак- тически одну и ту же картину. Здесь повторяются не толь- ко бифуркации удвоения периода, но и цепочки «хаос» —* —> «порядок», «порядок» —> «хаос». Отмеченная регуляр-
§4. Поведение нелинейной системы xt+i=/za?t(l—аД 37 ность и повторяемость этих структур (самоподобие) позволяет надеяться на получение простого описания форм хаотического поведения, обнаруженных при исследовании!! различных нелинейных моделей. Подводя итог приведенного анализа разнообразных яв- лений, порожденных нелинейной моделью динамики популя- ции, можно сказать следующее. Одна группа явлений связа- на с существованием регулярных и весьма упорядоченных процессов типа предельных точек покоя или циклов. В та- ких системах царит порядок, позволяющий по данным о про- шлом и настоящем предсказывать будущее. Другую группу составляют хаотические процессы, возможности предсказа- ния которых являются весьма ограниченными. Как правило, в анализе таких систем используется статистический подход (см., например, [14]), позволяющий получить лишь некото- рые усредненные характеристики. При этом между хаосом и порядком существует глубо- кая внутренняя связь. Хаотическое поведение возникает как предел усложняющейся последовательности периодических движений. Рассмотренная нами простейшая одномерная мо- дель с квадратичной нелинейностью наглядно это демон- стрирует. Дополнительный материал по обсуждаемым здесь во- просам, включая другие сценарии перехода к хаосу, можно найти, например, в работах [6], [12], [13], [15], [21], [23]. Упражнения 4.1. Из условия |(</’2(^з,4))/| = 1 найти бифуркационное значение /л2, при котором теряет устойчивость 2-цикл. 4.2. Из условия |(^’2(^з,4))/| = 0 найти значение ц, соот- ветствующее наиболее устойчивому 2-циклу.
4.3. Для системы xt+i = pg(xt), д(х) = < 2х, 2-2z, а) указать интервал значений р, при которых элементы последовательности xt не выходят из отрезка [0,1]; б) найти точки покоя и исследовать их устойчивость в зависимости от д; в) построить бифуркационную диаграмму. 4.4. Для системы <р(х) = рх(1 — х2) найти точки покоя и исследовать их устойчивость. § 5. Хаотическое поведение системы 5.1. Поведение критических значений парамет- ра. В предыдущем параграфе был подробно описан сцена- рий поведения динамической системы, порождаемый итера- ционным процессом xt+1 = д xt(l - Xt). (5.1) Этот сценарий таков: при определенных значениях па- раметра д, Д1,д2,...,Дп происходит бифуркация системы, сопровождаемая удвоением числа предельных точек-аттрак- торов итерационного процесса (5.1), а именно: при д = д„ точки цикла порядка 2П-1 становятся неустойчивыми и ро- ждается цикл порядка 2П. Приведем приближенные числен- ные значения для нескольких критических значений пара-
метра д : Д1 < С д < с 3.449499. •• = Д2 — цикл порядка 2; Д2 < с д < С 3.544090. •• = Дз — цикл порядка 22; РЗ < с д < С 3.564407. = Ш — цикл порядка 23; Д4 < с д < С 3.568759. •• = Дб — цикл порядка 24; Д5 < с д < ; 3.569692. • • = Де — цикл порядка 25; Дб * с д < : 3.569891. • • = Д7 — цикл порядка 26; д7 < с д < ; 3.569934. = № — цикл порядка 27. Американский физик М. Фейгенбаум, экспериментируя на калькуляторе с процессом (5.1), обнаружил, что последо- вательность д„ критических значений параметра сходится к некоторому предельному значению Доо с геометрической скоростью, которая характеризуется величиной lim Р“П Цп—1 п->оо Дп+1 ’ Р-П = /, (5.2) и вычислил приближенно ее значение, которое составило 6 = = 4.6692016... [19]. Величина 8 = 4.6692016... получила название первой универсальной константы Фейгенбаума. Можно приближенно вычислить предельное значение Доо последовательности критических значений параметра д, ис- пользуя формулу (5.2) и найденные значения дп (п — = 1,2,...). Для этого перепишем (5.2) в виде приближенного соотношения Р"П Р'П— 1 ~ ^(/^п+1 Р'п)
и положим в левой части рп ~ р^,, а в правой цп+1 ~ р^. Подставляя найденные значения рп = ps, pn_i = /./7, имеем mUoo др8 -нт _ 4.669202 • 3.569934 - 3.569891 <5 — 1 ~ 3.669202 ~ 3.569946. Дальнейшие эксперименты показали, что при р > р^ поведение итераций (5.1) становится абсолютно неупорядо- ченным, или, как принято говорить, хаотическим. Чтобы пояснить хаотический характер итераций при р > примем в процессе (5.1) р = 4 и возьмем в каче- стве начальной точки xq = sin2 тг/З, где/? — иррациональное число. Тогда, подставляя xq в формулу (5.1), последовательно находим а?1 = 4 sin2 7г/3(1 — sin2 лД) = sin2 2яД, Хэ = sin2 22лД, ..., хп = sin2 2пяД. Полученные соотношения показывают, что хп ведет себя как псевдослучайная величина, которая хаотично заполняет от- резок [0,1]. Эту ситуацию хорошо иллюстрирует рис. 5.1, на котором представлены итерации (5.1) при р = 3.9 после до- статочно большого числа шагов (N = 1000). Видно, что хп достаточно плотно заполнили отрезок [xmjn, zmax], где ^min $max = <^(z) = Zf 2=3.9 4 = 0.975, /1=3.9 Ж—X max = 0.095. /1=3.9 м2 Л _ /А 4 \ 4/
Рис. 5.1 5.2. Показатель Ляпунова. Введем некоторую ха- рактеристику — функцию параметра д, по знаку которой можно судить, в какой области значений параметра находит- ся система: область регулярного поведения, точки бифурка- ции или область хаоса? Строгому определению предпошлем некоторые наводя- щие соображения. На отрезке [0,1] возьмем две близкие точ- ки xq, х0 + £, где £ — малый параметр. После п итераций процессом (5.1) точки займут положение </?”(то), ^(^o+s) со- ответственно. Введем величину А, для которой приближенно выполнено соотношение г K(zo+£) -<^(ЯО)| exp [nA] ~.
Поскольку правая часть соотношения при s —> О стремит- ся к производной то можно определить функ- цию А*’,а:о(ц), не зависящую от п: х^°М= пт 11п!ады'|. п—>оо Введенная характеристика называется показате- лем Ляпунова. Этот показатель будет интересовать нас прежде всего как функция параметра ц. Однако ясно, что он зависит также от начальной точки хо и функции 99, хо- тя можно предположить, что в результате осреднения по п зависимость от Xq будет слабой. Теперь, с учетом формулы (2.3), показатель Ляпунова принимает вид п— 1 А^»(м)= lim (5.3) п—>ОО г=0 Утверждение 5.1. Если A(/z) < 0, то процесс (5.1) осуществляет периодический режим с циклом некоторого порядка. Если А(ц) = 0, то это соответствует режиму бифуркации, т. е. удвоению периода. Если A(ju) > 0, то по- ведение процесса хаотическое. Рис. 5.2. иллюстрирует факт, составляющий содержание утверждения 5.1. Как будет показано в следующем параграфе, явле- ния бифуркации, удвоения периода и перехода к хаотиче- скому поведению системы не являются уникальным свой- ством процесса (5.1), а имеют место для целого класса функций 9?м(ж). В частности, аналогичным свойством об- ладает итерационный процесс с треугольным отображени- ем <р^(х) = ^(1 - |1 - 2т|).
У пражнения 5.1. Вычислить показатель Ляпунова для треугольного отображения Ыж) = ^(1 - |1 - 2т|). 5.2. Вычислить приближенно константу 6 для треуголь- ного отображения. Сравнить полученную величину с первой универсальной константой Фейгенбаума 6 = 4.669201. § 6. Универсальность в поведении систем 6.1. Определение константы а. М. Фейгенбаум [19] открыл еще одну константу. Она определяется следующим образом. При переходе от цикла порядка 2n (zz = 1,2,...) к циклу порядка 2n+1 в процессе (5.1) при некотором значе-
нии параметра /j, = Мп точка х* — 1/2 становится аттрак- тором (см. рис. 7.1, 7.2). Обозначим через dn алгебраиче- ское расстояние (т. е. с учетом знака) от аттрактора х* = 1/2 до ближайшего к нему аттрактора (см. рис. 4.6). Оказыва- ется, что существует предел lira = —а, а > 0, (6.1) п^оо ап+1 который определяет вторую универсальную кон- стан т у. Ее численное значение с точностью до семи знаков после запятой есть а = 2.5029078. Позднее было установлено (см., например, [25]), что эту константу также можно опре- делить соотношением d! lim —— = —а, (6.1а) п^°° 4+1 где d'n — расстояние между ближайшими к 1/2 аттракторами при значении параметра д = дп+1, когда цикл порядка 2П разрушается и рождается цикл порядка 2n+1. 6.2. Универсальность констант. Почему же кон- станты <5, а назвали универсальными? Изложим свою версию происхождения этого термина. В своей работе [19] М. Фейгенбаум пишет, что после его чис- ленных экспериментов с процессом (5.1) и обнаружением за- кономерностей, характеризующихся константами (5.2), (6.1), коллега по Лос-Аламосской лаборатории П. Стайн обратил его внимание, что итерационный процесс xt+1 = д sin7r;rt (t = 0,1,2,...) с функцией шага <дм(т) = д sin хх также обладает свойством удвоения периода (числа предельных точек) при изменении
параметра д и, что самое удивительное, с теми же констан- тами 6, а ! Дальнейшие исследования показали, что явление би- фуркации (резкого изменения поведения системы ) обнару- живается для двумерного процесса Хеннона: zs+i = 1 - ax2t + yt yt+1 = bxt, причем последовательность an критических значений пара- метра а (для некоторых фиксированных значений Ь), при которых происходит качественный скачок в процессе (6.2), подчиняется соотношению (5.2). Конечно, в данном случае можно говорить о некотором сходстве систем (5.1) и (6.2), например, о наличии квадра- тичной нелинейности. Однако та же закономерность при изменении одного из параметров была обнаружена в знаменитой системе Лоренца х = —ах + ау, у = гх — у — xz, z = ху — bz. А именно: при а = 10, 6 = 8/3 отношение Гп ~ Гп-1 ап — гп+1 - гп критических значений параметра гп очень хорошо аппрок- симирует константу 6, например 65 ~ 4.670 (см. [7] ). Рис. 17.6 иллюстрирует сложное поведение решения си- стемы Лоренца. На этом рисунке изображена траектория, порождаемая системой Лоренца при г = 28, а = 10,5 = 8/3. Оказывает- ся, что а) она притягивается к ограниченной области в фа- зовом пространстве (аттрактору Лоренца); б) движение ее
блуждающее, т. е. траектория делает один виток направо, затем несколько витков налево, затем снова направо и т. д.; в) траектория очень чувствительна к малым изменениям на- чальных условий. Факт геометрической скорости сходимости (с констан- той 6) критических значений параметра был подтвержден и во многих других системах. Это обстоятельство, по-ви- димому, породило у некоторых исследователей уверенность в том, что переход ДС к хаосу происходит по одному и то- му же сценарию и может быть описан одной и той же константой 8. Что касается константы а, которая отвечает за структуру аттракторов в одномерном итерационном про- цессе, то были найдены другие функции <£>Дж), отличные от 9?м(ж) = цж(1 — х),<Рц(х) = /zsin-Trz, которые генериру- ют итерационные последовательности (орбиты) со свойством удвоения числа предельных точек, расположение которых на отрезке подчинено соотношению (6.1). По-видимому, М. Фейгенбаумом впервые была высказа- на гипотеза, которая до сих пор строго не доказана, но и не опровергнута, что установленные закономерности, характе- ризуемые константами 8 = 4.6692016..., а = 2.5029078..., справедливы по крайней мере для всех унимодальных функ- ций с квадратичной нелинейностью в точке максимума, т. е. функций </?(т), удовлетворяющих условию 9?(.т) — <£>(ж) = О(х — ж)2 в окрестности единственной точки максимума ж на отрезке. Таким образом, можно сказать, что сценарий удвоения ' числа предельных точек процесса (5.1) и переход к хаосу не зависят от конкретной функции р(х) в пределах клас- са с квадратичной нелинейностью и, следовательно, носят универсальный характер.
Таким образом, в математике, наряду с числами е = =• 2.718281..., тг = 3.141592..., найдены еще две замеча- тельные иррациональные константы а и 6, которые характе- ризуют сценарий бифуркации в одномерном итерационном процессе с функцией перехода <рр(х), обладающей квадра- тичной нелинейностью. Вскоре выяснилось, однако, что численные значения констант о, <5, определяемых соотношениями (5.2), (6.1), за- висят от характера нелинейности переходной функции рр(х). Проще всего это можно пояснить, если от системы (5.1) перейти к процессу vt+1 = 1 - <z|v£|2, (6.3) где функция перехода 'фа(х) = 1 — ах2 задана на отрез- ке [—1,1]. Заменой vt = —^-x(xf — к), ° = -- один JJL Л A ~L процесс сводится к другому, но процесс (6.3) более удобен, так как в нем ясно виден квадратичный характер нелиней- ности и его легко можно варьировать, меняя степень. Оказалось, что если теперь рассмотреть одномерный итерационный процесс t>t+i = 1 - a\vt|₽ = -ipa(vt) (6.4) при p > 2 (например, p = 3 ), то он, как и (6.3), испыты- вает бифуркации (удвоение периода), существуют пределы (5.2), (6.1), но численные значения констант 6, а уже дру- гие. Итак, константы 6, а зависят от степени р, и следует писать 8Р, ар, отмечая зависимость констант от показателя степени р. 6.3. Класс функций, порождающий бифуркацию удвоения. Естественно возникает вопрос, можно ли опи-
сать совокупность всех функций ^(ж), для которых ите- рационный процесс претерпевает бифуркации, причем сце- нарий удвоения периода, структура аттракторов и переход к хаотическому поведению системы характеризуется некото- рыми константами д, а, которые определяются соотношени- ем (5.2), (6.1). Попытаемся это сделать из существующего к настоящему времени теоретического анализа и вычисли- тельного опыта. Обозначим этот класс через К и будем считать для опре- деленности, что речь идет о функциях на отрезке [0,1]. Не лишено достоверности (это лишь гипотеза!) следую- щее утверждение (для простоты в записи функции мы опус- каем параметр ц ). Утверждение 6.1. Для того чтобы функция <р(х), где <р: [0,1] —> [0,1], трижды дифференцируемая почти всю- ду, принадлежала классу К, достаточно выполнение следу- ющих условий: 1) существует единственная точка максимума х* G (0,1) этой функции; 2) Д(х) >0 \/х > 0 € [0, х*), <Д{х) < 0 Vx(t*, 1]; 3) производная Шварца 5И(х) = <(*) </'(*) '(*) 3 2 почти для всех х Е [0,1]. Строгое доказательство этого факта (как, впрочем, и пример, опровергающий его) отсутствует, но некоторую аргументацию в пользу этого утверждения, в частности важ- ности условия 3, можно найти в работе Г. Шустера ([23], см. приложение 3).
Однако справедлива следующая теорема, в которой для 6 С3[0,1] с условиями 1 — 3 утверждается несколь- ко более слабое свойство, чем принадлежность классу К. Теорема 6.1. Пусть <р € О'3[0,1] и выполнены усло- вия 1-3 из утверждения 6.1. Тогда существует не более одного устойчивого цикла. В случае его существования эле- менты цикла являются предельными точками последова- тельности ип = <рп(х*), где х* — точка из условия 2 (<р'(х*) = 0 ). Доказательство теоремы приведено в работе Д. Зингера. Там же построен пример, показывающий, что без условия 3 утверждение, вообще говоря, не имеет места (см. [31]). Замечание 6.1. Условие 3 в утверждении 1 нельзя, вообще говоря, заменить более слабым — Sl[</?](«) < 0. Например, для линейной функции 5[т?](х) = 0, но эта функция, очевидно, не принадлежит классу К. С другой стороны, для кусочно-линейной функции tp(x) = р(1 —11 —2т|) также 5[</?](т) = 0 (за исключением точки х = 1/2, где она является недифференцируемой), тем не менее эта функция из класса К, в чем непосредственно можно убедиться, исследуя итерационный процесс с такой переходной функцией. Следовательно, условие 3 не является необходимым для е. к. Упражнения 6.1. Доказать следующие свойства производной Швар- ца: 1) если f — линейная функция, аveC3,roS[/MJW = = S[</?](x) для всех х; 2) если — дробно-линейная функция, то S^]^) = 0;
3) если £[99] (z) < 0 (5'[9?](ж) > 0) и 5[/](ж) < О (ЭД W ? 0), то S(№)]W < 0(ЭД(9)]Ы > 0); 4) если для функции 99: [0,1] —> [0,1] ^[^(ж) < О, то для любого целого п выполняется неравенство [9/1] (ж) < 0. Здесь 9?п(ж) = 9?(9р(. .. у?(ж)))у п раз 6.2 А. Проверить аналитически, что следующие функции удовлетворяют условиям утверждения 6.1: 1)Ы‘Т) = W1 - ж)> х е [о, 1], 0 ц 4; 2)^(^) = ехр-^2, х G [-1,1], 0 1; З^Дх) = ц sin ж, х G [0,7г], 0 ц 1; 4)9?м(а;) = 1 — ц|ж|р, х G [—1,1], 0 < д < 1, 1 < р < сю. Б. Численно найти показатель Ляпунова на некоторой сетке по р, для функций из пункта А. Построить бифурка- ционную диаграмму. В. Построить другие функции из класса К и протести- ровать их аналогично пунктам А, Б. Г. Построить функцию <р(ж), удовлетворяющую услови- ям 1, 2 утверждения 6.1 и не принадлежащую классу К. Д. Построить функцию <р(т) из класса К, удовлетво- ряющую условиям 1 — 2 и не удовлетворяющую условию 3 утверждения 6.1. § 7. Уравнение удвоения Преобразуем формулу (6.1) для второй универсальной константы, выразив ее непосредственно через итерируемую функцию 9?M(x). Эта константа была определена как рассто- яние между Xq = и ближайшей к ней точкой хкогда обе
являются элементами аттрактора. Данная ситуация харак- теризуется понятием суперцикла. Определение 7.1. Суперциклом называется цикл, в котором k=a>o= О’ где Xq — некоторая предельная точка-аттрактор (элемент цикла), Мп — значение параметра, при котором производ- ная 2п-й итерации функции обращается в нуль. На основании формулы (2.3) имеем 2n—1 1®=*о= П = О' г=0 Так как х = 1/2 — единственная точка, в которой (х) = = 0, то х = 1/2 является элементом цикла. Эту ситуацию иллюстрирует рис. 7.1 (а, б). Заметим, что аналогично определяется суперцикл для функции "фа/х) = 1 — ах2, где роль точки х = 1/2 игра- ет х = 0. Принимая во внимание факты, изложенные в § 4, мож- но написать следующие соотношения: = — ^,^2 = ^2(5) — ip • • • ,<4 = <Рмп (^) — а с учетом формулы (6.1) — цепочку приближенных ра- венств , _ du ____ dn—\ __ d^ n+1 ~ (^) ~ (^Zr ~ (=^T что позволяет заключить о существовании предела = Jim <г„+1(-о)" = Ип+1(|) - |](-о)". (7.1)
Рис. 7.1 Сделаем замену переменных х = Тогда точка Z Z z перейдет в точку 0, а функция в /M(i) = ^(1 — i2). Сохранив для новой переменной старое обозначение, пере- пишем формулу (7.1): </1= 1ЙП J0) - 0]. п—>се Численный анализ показывает, что для любого х суще- ствует также предел д^х) = lira (-а)7мя+1(тА- п—>оо п+± (—а)
Введем семейство функций gi — lim (—а)п/м ..................). Определим оператор удвоения T(j) = -о9(9(-§)) и покажем, что справедливо соотношение 5г-1(ж) = Тд^х). Действительно, преобразуя выражение для получаем представление = lim (—а)п/м . (, 37 а V ' n-»ooV ’ JMn+^-l к(-а\п п—>оо п^г 1 ~ (а)П Г2”-1 (________X______ (—а)'г“1 ,/м’г+-1 ^(_a)(_a)n-iy — ('лЛт из которого, в предположении справедливости предельного перехода, получаем = (-а)^(^(-тг^)) = T(9i). В условиях существования предела при i —> оо приходим куравнению удвоения 9(х) = -о9(9(-§)) = Г(9), (7.2)
которое называют уравнением Цвитановича- Фейгенбаума (см. [25]). Вывод уравнения (7.2), при- веденный выше, принадлежит, по-видимому, Фейгенбауму (см. [19]). Нетрудно проверить, что если д(х) — решение уравне- ния (7.2), то /(х) = 7^(т/7), где 7 — числовой параметр, также является решением. Чтобы избавиться от неоднознач- ности, нормируем решение, потребовав выполнения усло- вия <?(0) = 1. Подставляя условие нормировки в уравне- ние (7.2), приходим к соотношению а = -1/р(1), (7.3) которое позволяет переписать уравнение в виде ^(1)£?(т) = р(р(р(1)ж)), (7.4) т. е. избавиться от 2-й универсальной константы. Уравнение (7.4) можно решать приближенными метода- ми (см. § 8). Найдя решение, можно вычислить универсаль- ную константу а по формуле (7.3). Как обстоит дело с существованием решения? Этому во- просу посвящены многие статьи. Мы ограничимся формули- ровкой одного утверждения из работы О. Ландфорда, дока- зательство которого было получено с привлечением компью- терных вычислений (см. [28]). Теорема 7.1. Существует функция д, аналитическая и четная в круге {z € С : |z| < \/8}, сужение которой на [—1,4-1] есть решение уравнения (7.4), причем производ- ная Шварца функции д отрицательна на отрезке [—1,4-1]. При дальнейшем рассмотрении было установлено [24], что решение уравнения (7.4) не единственно, каждое реше-
ние представимо в виде степенного ряда ОО д(х) = 1 + ^ШР (7-5) i=l при некотором 1 < р < оо, а найденная по этой функции уни- версальная константа ар = —(1/^(1)) отвечает за скейлинг аттракторов в динамической системе, описываемой процес- сом (6.4) со степенной нелинейностью р. § 8. Численное нахождение универсальной константы ар Ограничимся константой а2, которая соответствует квадратичной функции 1ра(х) = 1 — ах2, х G [—1,4-1]. На основании представления (7.5) из предыдущего парагра- фа решение уравнения удвоения (7.2) следует искать в виде степенного ряда по четным степеням. При грубой аппрок- симации можно ограничиться двумя членами, что, с учетом нормировки <?(0) = 1, дает представление д(х) = 1 + Ьх2 с неизвестным коэффициентом 6. Для нахождения b подста- вим выражение для д(х) в уравнение (7.2) и приравняем ко- эффициенты при степенях х, отбрасывая члены четвертого порядка в правой части. Тогда имеем 1 + Ьх2 = —а(1 + 6(1 + 6()j)2)2), откуда 1 = —а(1 + 6), а — —2Ь. Исключая а, приходим к уравнению 2Ь2 + 26 — 1 = 0, поэтому 61э2 = (—2 ± л/12)/4. Из некоторых соображений следует, что нужно взять от- рицательный корень, что дает 6 ~ —1.36, следователь- но, «2 ~ 2.72.
Для получения более точного приближения, естествен- но, нужно взять отрезок степенного ряда с большим числом членов и решать систему нелинейных уравнений, получен- ную с помощью схемы коллокации. Опишем подробнее эту технологию в случае динамической системы (6.4) с нелиней- ностью степени р. На основании представления (7.5) ищем решение в виде отрезка ряда п ди = 1 + £е,|тГ г=1 (8-1) Зададим некоторую сетку {тД,г = 1,2, ...,п на отрез- ке [— 1,4-1]. Подставим выражение (8.1) для д(х) в уравне- ние (7.4) и потребуем выполнения равенства для каждой точ- ки х — Xi сетки. Окончательно приходим к системе нелиней- ных уравнений п п (1 + £&)(1 + £ед,П- i=l г=1 п п п -1 - £&ii + Е^-К1 + ЕьИ’Т = о. з = тд (8.2) г=1 г=1 г=1 относительно коэффициентов {£Д разложения (8.1). Для приближенного нахождения искомого вектора £ = = (€ii€25 • • •, €п) можно применять подходящий метод ре- шения систем п нелинейных уравнений с п неизвестными. К. Бриггс использовал для этой цели классический метод Ньютона и, привлекая специальные программные средства, позволяющие проводить вычисления с 200 знаками, вычис- лил константы ар, 6Р для р = 2,3,..., 12 со 100 знаками после запятой (см. [24]). Полученные значения для ар,ёр с тремя знаками после запятой приведены в таблице 1.
Таблица 1 р 2 3 4 5 6 7 2.502 1.927 1.690 1.555 1.467 1.405 4.669 5.967 7.284 8.349 9.296 10.222 Р 8 9 10 И 12 1.358 1.321 1.291 1.267 1.246 8Р 10.948 11.768 12.341 13.076 13.535 При больших значениях р метод Ньютона расходится (неустойчивость!). А. Б. Смирнова [18], используя итеративно регуляризо- ванные методы Ньютона и Гаусса-Ньютона (см. [32]), повто- рила результаты К. Бриггса с точностью до пяти знаков по- сле запятой и вычислила константу ар для р = 13,14,..., 24 с той же точностью. Полученные значения с тремя знаками после запятой приведены в таблице 2. Таблица 2 Р 13 14 15 16 17 18 Qp 1.229 1.213 1.201 1.189 1.178 1.169 р . 19 20 21 22 23 24 Qp 1.161 1.153 1.146 1.140 1.134 1.129 Приведем выражение для функции д(х), полученной в результате решения системы (8.2) при р = 2 (квадратич- ный случай) и п = 7 с девятью знаками после запятой д(х) = 1 - 1.527632997т2 + 0.104815194т4 + 0.026705673т6- - 0.003527413т8 + О.ОООО81581т10 + 0.000025368т12-
Вычисленное по формуле а = —1/^(1) оказывается равной а = 2.502907875.... Заметим, что коэффициенты степенного ряда довольно быстро убывают, поэтому высокую точность для решения д{х) (а следовательно, для а) можно получить при сравнительно небольшой размерности системы (8.2). § 9. Численное нахождение универсальной константы 5Р 9.1. Общий метод. Перейдем теперь к вычислению универсальной константы которая отвечает за скейлинг по параметру д (см. формулу (5.2) из § 5). М.Фейгенбаумом было доказано следующее утверждение. Утверждение 9.1. Универсальная константа 5Р сов- падает с наибольшим по модулю собственным числом опе- ратора производной L{g)=T'{g), где Т — оператор удвоения {см. (7.2)), ад— решение уравнения универсальности (7.4) для случая степенной нелинейности с показателем р. Доказательство приведено в книге Г. Шустера [23] (см. §32). Вычислим производную Гато оператора удвое- ния Т(д) = -ад(д(-%У). L(t/)/i==lim(— {g+^(g(-^+Xh{-^-g{g{-^} ---------------------- > = lim(—а)<------------;-------------1- А—>0ч А А/г(9(-§) + АЛ(-§))] = -Л<№(-а)) ))} (9-П
Согласно утверждению 9.1, универсальная константа 8Р является решением задачи на собственные значения для ли- нейного оператора £(<?), т. е. L(g)h = 6h. (9.2) Как и в случае константы сначала найдем грубое приближение для 52- Для этой цели воспользуемся следу- ющим приемом. При решении задачи (9.2) будем считать функцию h константой. Тогда, подставляя в (9.2) найденное выражение для L(g) и полагая х = 0, приходим к соотноше- нию —a(gf(l)h + h) = 6h, которое после сокращения на h переходит в следующее: -а(У(1) + 1) = <5. (9.3) Для нахождения </(1) продифференцируем два раза уравнение универсальности д(х) = -ад(д(-£у), что дает s'W = я'(9(-г)) • 9'(-§), = s"(9(-5» (э'(-§))2(-з) + э'(9(-2)) 9"(§)(-з). Подставим в последнее соотношение х = 0, воспользуемся условием д'(0) — 0, которое следует из представления реше- ния (7.5), и сократим на д"(0). Окончательно получаем 1=9'(1)(-1М Вместе с (9.3) это дает выражение для 6: о? — а = 6.
Воспользовавшись найденным в предыдущем параграфе значением ~ 2.72, получаем ~ 4.67. Для сравнения численное значение этой константы с девятью знаками есть 62 = 4.669016091... 9.2. Прямой алгоритм. Можно предложить прямой метод (см. [25]) вычисления константы 62, в котором нет необходимости решать уравнение удвоения и находить мак- симальное собственное значение производной оператора Т. Он основан на гипотезе К. Бриггса [25], что константу д2 можно вычислять также по формуле 62 = (9.4) где ai — не критические значения параметра, в котором про- исходит бифуркация удвоения периода в процессе с функци- ей V’a(^) = 1 — ах<2 (ср. с (5.2)), а значение параметра, при котором реализуется суперцикл порядка 2г. Рассмотрим последовательность полиномов опреде- ляемых рекуррентно: 60(а) = 0, Ьк(а) = а - [5fe_i(a)]2, /с =1,2,.... Утверждение 9.2. Пусть к — 2г. Процесс с функци- ей фа(х) имеет суперцикл порядка к = 2г тогда и только тогда, когда bk(a) = 0. Доказательство вытекает из определения суперцикла для функции и того факта, что х = 0 является корнем полинома фа(х) = 0 (см. Рис- 7.1).
Установленное свойство позволяет находить значения параметра а, входящие в формулу (9.4), как корень уравне- ния bk(a) = 0. Это позволяет выписать следующий алгоритм приближенного вычисления <52 О . ®г—1 2 'по /п tr\ ®г—1 4 j $ 2,3,..., (9'5) j+1 j b2i(aJ{) а------------ j = 0,1,2,..., ^к(а) ~ 1 26^_1(а)Ьд._1(а), к — 1,2,3, (ц = lim a3i} ci —1 —2 г -i. 6 =--------------------, o2 = lim о . &i — CLi—l i—>oo (9.6) (9.7) (9.8) (9-9) Поясним эти формулы. Соотношение (9.5) предназначе- но для вычисления начального приближения при вычисле- нии очередного значения параметра а^. Формулы (9.6)-(9.8) реализуют метод Ньютона для нахождения корня уравне- ния bfc(a) = 0. Наконец, формула (9.9) есть просто переза- пись формулы (9.4). Оказывается, что описанный алгоритм для нахожде- ния 62 можно использовать для вычисления универсальной константы а2. В работе К. Бриггса [25] была предложена формула ,. bi+1(ai+1) •lim "и( \ ^2 а2’ (9.10) где Ь((йг) определяются рекуррентно соотношением (9.7). Таким образом, зная приближенное значение предела в соотношении (9.10) и одной из универсальных констант, можно вычислить другую.
Существует ли какая-либо аналитическая функциона- льная зависимость между константами «2 и <$2? Являются ли они алгебраическими числами, т. е. корнями многочлена с целыми числами? До сих пор эти вопросы остаются откры- тыми. Однако на основе численного моделирования установ- лен следующий интересный факт (см. об этом [25]). Утверждение 9.3. Если бы константа а2 была кор- нем полинома степени 20 или меньше с целыми коэффици- ентами, то по крайней мере один из коэффициентов пре- восходил бы 2 • 1015. Если бы константа была корнем по- линома степени 20 или меньше с целыми коэффициентами, тогда по крайней мере один из коэффициентов превосходил бы 5 • 1015. Иными словами, универсальные константы (*2,62 не являются корнями полиномов с целыми коэффициента- ми до 20-й степени, коэффициенты которых не превосхо- дят 2 • 1015 или 5 • 1015 соответственно. У пражнения 9.1. Вычислить константу 5г, проведя предварительно дискретизацию задачи (9.2) и используя итерационный про- цесс хп+х = Ахп, lim ^—^ = 82, П->ОО (xn.Xn) где А — оператор, полученный в результате дискретиза- ции оператора L(g), g — решение уравнения удвоения из §8 (см. формулу (8.7)). 9.2. Используя алгоритм (9.5)-(9.9) и начальные аппрок- симации oq = 0,й1 = 1,£>о = 0,51 = 3.2, вычислить Ui (г = = 2,3,..., 10) и 510 ~ 52.
9.3. Используя формулу (5.10) и найденное выше ^2~<510, вычислить приближенное значение а2- 9.4. Модифицировать процесс (9.5)-(9.9) и выписать ал- горитм для вычисления <53.
Глава 2 Элементы фрактальной геометрии § 10. Комплексные динамические системы Почему геометрию называют холодной и сухой? Одна из причин заключается в ее неспособности описать форму облака, горы, дерева или берега моря. Облака — это не сферы, горы — это не ко- нусы, линии берега — это не окружности, а кора не является гладкой, а молния не распространя- ется по прямой... . Природа демонстрирует нам не просто более высокую степень, а совсем другой уровень сложности. Число различных масштабов длин в структурах всегда бесконечно. Существова- ние этих структур бросает нам вызов в виде труд- ной задачи изучения тех форм, которые Эвклид отбросил как бесформенные, — задачи исследова- ния морфологии аморфного. Б. Мандельброт A fractal can be beatiful as a sunset. A. Phillips 10.1. Множества Жюлиа, Фату и Мандельбро- та. В предыдущих параграфах были рассмотрены динами-
ческие системы, порождаемые итерационными процессами с функциями действительного переменного (одномерный случай). В частности, были рассмотрены функции шага = дт(1 - х), х G [0,1]; = 1 - av2, х 6 [-1,1], которые сводятся одна к другой линейной заменой перемен- ных при некоторой связи параметров (см. § 6). Если в ите- рационном процессе vt+i = 1 - av2 сделать замену vt = zt/(—a) и положить (—а) = с, то прихо- дим к динамической системе вида zt+1 = z2 + c. (10.1) Если zt принимает действительные значения, то наблю- дается аналогичная картина бифуркации удвоения периода, которая характеризуется теми же константами Пусть теперь zt 6 С, с 6 С, т. е. принимают комплекс- ные значения. Если воспользоваться представлением zt = — xt + iyt, с = р + iq, то после подстановки в (10.1) и вы- деления действительных и мнимых частей приходим к дву- мерной системе Zt+1 = х2 - у2 + р Vt+i = 2xtyt + q, эквивалентной системе (10.1). Переход к комплексным системам привел к появле- нию множеств очень необычной, можно сказать, фантастиче- ски причудливой формы. Речь идет, конечно, о множествах Жюлиа, Фату и Мандельброта.
Комплексные динамические системы в начале прошло- го века интенсивно исследовали французские математики Г. Жюлиа и П. Фату. Их исследования не получили даль- нейшего развития и фактически были забыты на несколько десятилетий, «поскольку в отсутствие современной компью- терной графики было почти невозможно передать их тонкие идеи» [16]. Всплеск интереса к динамическим системам, порожда- емым итерационными процессами с рациональными функ- циями перехода, произошел в конце 1970-х годов, когда Б. Мандельброт и Дж. X. Хаббард получили компьютерное изображение некоторых множеств Жюлиа, т. е. множеств, имеющих ограниченную траекторию. Кроме того, Б. Ман- дельброт продемонстрировал множество, носящее теперь его имя, значений параметра с, для которых траектория нуля (критической точки функции z2+c ) остается ограниченной. Формы таких множеств настолько причудливы и чре- звычайно изящны, что никого не оставляют равнодушным и производят сильное эстетическое впечатление. Отличительной чертой этих множеств является самопо- добие: например, на границе основного множества Мандель- брота (см. рис. 10.1) можно увидеть многочисленные его ко- пии сколь угодно малого размера, все множества Жюлиа (см. рис. 10.3-10.9) состоят из самоподобных фрагментов. Можно также обратить внимание, что границы этих мно- жеств сильно изрезаны, что, по-видимому, дало основание Б.Мандельброту назвать такого типа множества фрак- талами (от fractus — изломанный). Этот термин также хорошо согласуется с английским словом fractional (дроб- ный), что подчеркивает еще одно характерное свойство та- ких множеств, связанное с дробностью размерности Хаус- дорфа (строгие определения подобия,, фрактала и размерно- сти будут даны в следующих двух параграфах).
Перейдем теперь к строгим определениям множеств Жюлиа и Мандельброта, приняв обозначение Рс(г) = z2 + с, р? = ре(ре(...ре)). п раз Определение 10.1. Наполненным множе- ством Жюлиа для Pc(z) называется множество К(РС) = {г € С : Vn |Рсп(г)| < R < оо}, где R — некоторая положительная константа; т. е. это мно- жество точек комплексной плоскости, орбиты которых огра- ничены. Граница множества К называется просто множе- ством Жюлиаи обозначается J(PC). Определение 10.2. Дополнение C\J(Pc) называется множеством Фату для отображения Рс. Определение 10.3. Множеством Мандель- брота для Pc(z) называется множество М = {с е С : Vn |Р”(0)| < R < оо}, где R — некоторая фиксированная положительная констан- та, т. е. М — множество параметров с, для которых орбиты нуля ограничены (итерации zn+\ = Рс(гп), п = 0,1,..со- держатся в круге фиксированного радиуса Р > 0 ). Множество Мандельброта представлено на рис. 10.1. Оно расположено в круге {z: \z\ С г}. Предоставим читате- лю судить, является ли это множество прекрасным или без- образным с художественной точки зрения, но определенно можно сказать, что оно имеет необычно причудливую фор- му с точки зрения классической геометрии. Можно дать эквивалентные определения множеств J(PC) и М на основе понятия чувствительности отображения.
Определение 10.4. Пусть X — метрическое простран- ство, а А — его подмножество. Говорят, что отображе- ние F: X X чувствительно к начальным данным на мно- жестве А, если 3/?>0:VyeA V Ое(у) 3peOs(y) Зп > О p(Fn(y), Fn(p)) > (3, Fn = F(F(.^..F)), п раз т. е. сколь угодно близкие точки на некотором шаге итерации разойдутся на расстояние, большее чем /3. Определение 10.1'. Множество Жюлиа J(PC) — это множество всех точек из С, на котором отображение Pc(z) чувствительно к начальным данным. — 1.Х, —2.4 -2.0 -1.6 -1.2 -0.6 -0.2 0 0.2 0.4 0.6 0.8
Определение 10.3'. Множество Мандельброта — это множество значений параметра с, для которого наполненное множество Жюлиа К(РС) связно. Доказательство эквивалентности определений и неко- торых других свойств множеств К(РС\ J(PC), М можно най- ти в монографии [26]. Замечание 10.1. Не следует думать, что множества Жюлиа и Мандельброта определяются только для отображе- ния Pc(z) = z2 + с. Просто исторически изучение комплекс- ных систем начиналось с отображения Рс(г), и все определе- ния были сначала введены для этого отображения. Аналогичным образом вводится определение наполненного множе- ства Жюлиа для произвольного отображения /: С —> С, а именно K(f) - {z : г € С, |/п(г)| < R < оо}. (10.2) Чтобы обобщить определение множества Мандельброта, напо- мним сначала определение критической точки г* отображения f(z) как корня уравнения — 0, т. е. f'(z*) = 0. Пусть те- перь fc(z) = f(z) + с, f'c(z*) = 0. Тогда множество Ман- дельброта для отображения fc и критической точки г* определяется как множество M(fc, Z) = {с : с € С, |/сп(г*)| R < оо}. (10.3) Таким образом, сколько критических точек (корней уравне- ния fc(z) = 0 ), столько и множеств Мандельброта. Поскольку нуль — критическая точка отображения Рс(^) = z2 + с, то данное определение согласуется с введенным ранее. Наиболее полно изу- чены множества K(f\ M(fc, г*) для случаев, когда f — полином или рациональная функция (см., например, [16], [26]). 10.2. Классификация множеств Жюлиа. В соот- ветствии с определением 10.3' при изменении параметра с
в пределах множества Мандельборта наполненное множе- ство Жюлиа остается связным. Более того, оказывается, что эти множества имеют некоторые характерные формы в за- висимости от того, из какой части множества М взят па- раметр с. Иными словами, положение параметра с во мно- жестве Мандельброта определяет свой специфический класс множеств Жюлиа. Это позволяет провести классификацию основных типов множеств Жюлиа. Проследим, как меняется множество Жюлиа, когда па- раметр с последовательно принимает действительные отри- цательные значения. Начнем с параметра с = 0.0 + 0i. Очевидно, что мно- жество Жюлиа в этом случае есть круг единичного радиуса с аттрактором z = 0.0 + 0г. При с — —0.25 + 0г это множе- ство уже представляет собой несколько деформированную окружность с аттрактором (устойчивой неподвижной точ- кой) z2 = 1 — \/2- Действительно, решая квадратное уравнение Рс(г) = z2 + с = z, находим два корня zii2 = (1 ± — 4с)/2 и проверяем усло- вие |.Pc'(z)| = \2z\ < 1. При с = —0.25 + 0г этому условию удовлетворяет лишь корень г2 = 1 — л/2. На рис. 10.2 изображены множества Жюлиа при раз- личных значениях параметра с. Обратим особое внимание на значение с = —0.75 + 0г. Напомним, что при переходе от функции = цх(1 — х) к Рс = z2 + с была установлена связь параметров д(ц - 2) = -4с, поскольку д(р. — 2)/4 = а (см. § 6) и а = —с (см. п. 10.1). По- этому с = —0.75 соответствует единственное положительное
значение параметра д = 3. А это, как было установлено в § 4, есть критическое значение параметра д, при котором проис- ходит бифуркация, т. е. удвоение аттракторов (предельных точек). Таким образом, с = — 0.75+ 0г также критическое значе- ние параметра с, и, следовательно, при дальнейшем умень- шении с, что соответствует увеличению параметра д, уже по- являются два аттрактора, которые условно помечены жир- ными точками на трех последних фигурах рис. 10.2. С другой стороны, первые четыре фигуры на этом рисунке соответствуют значениям параметра с, которые расположены внутри основной области Мандельброта (ти- па кардиоиды). Границы этих фигур представляют со- бой деформированную окружность, охватыва- ющую собой единственную притягивающую неподвижную точку. Оказывается, что подобная картина наблюдается для всех значений параметра с, не выходящих из кардиоиды. с = 0.0 + 0 г с = —0.25 + 0г с = —0.5 + 0г с = -0.75 + 0? с = -0.8 + 0г с = -0.9 + 0г
Например, при с = —0.12375 + 0.56508г множество Жюлиа, изображенное на рис. 10.3, также имеет вид деформирован- ной окружности. Возьмем теперь значение с = —0.12 + 0.74г, которое рас- положено в центре самой большой почки («луковки») сверху от основной части М. Множество Жюлиа при этом значении параметра с представлено на рис. 10.4. Цикл периода 3 появ- ляется здесь в результате трифуркации неподвижной точки, когда параметр с переходит из основной части в соответству- ющую почку. Значение с = —0.481762 — 0.531657г соответствует ме- сту прорастания средней нижней почки, дающей устойчивые циклы периода 5, когда с переходит внутрь почки. Для этого значения параметра реализуется так называемый параболический случай. Соответствующее множе- ство Жюлиа изображено на рис. 10.5: здесь неподвижная точка является точкой ветвления. Значение с = —0.39054 — 0.58679г является граничной точкой основного множества Мандельброта, которое не сов-
Рис. 10.4 Рис. 10.5 падает с точкой прорастания почек. Внутри области, огра- ниченной множеством Жюлиа, процесс итераций протекает следующим образом: сначала итерационные точки переска- кивают из меньших, периферийных, точек в большие до тех пор, пока не попадут внутрь диска, содержащего неподвиж- ную точку. Этот диск назван диском Зигеля (в честь немецкого математика К. Л. Зигеля). После того как точки попадают в диск, они начинают вращаться по своим инва- риантным окружностям, не покидая их (см. рис. 10.6). Четыре примера, рассмотренные выше, охватывают фактически все типичные случаи, когда множество Жюлиа содержит внутренние точки. Итак, - - если с лежит внутри основного множества Мандельбро- та, то множество Жюлиа представляет собой дефор- мированную окружность, охватывающую единственную притягивающую неподвижную точку-аттрактор (рис. 10.2, 10.3); — если с лежит внутри одной из почек, то множество Жюлиа состоит из бесконечного числа (фрактально) де- формированных окружностей, охватывающих устойчи- вый цикл некоторого порядка (рис. 10.4);
— если с является точкой прорастания почки, то имеет ме- сто параболический случай, граница имеет усики, дотя- гивающиеся до устойчивого аттрактора (рис. 10.5); — если с является любой другой точкой границы или поч- ки, то реализуется случай с диском Зигеля (рис. 10.6). Рис. 10.6 Оказывается, возможен еще один класс множеств Жю- лиа — так называемые кольца Эрмана, которые не ре- ализуются в случае отображения Pc(z), но могут возникать при других отображениях (см. разд. 3 книги [16] X. Пайтгена и П. Рихтера, где изложены результаты Д. Салливана о пол- ной классификации множеств Фату и Жюлиа). Рассмотренные случаи относятся к множествам Жю- лиа с внутренними точками. Однако существуют множе- ства Жюлиа, не содержащие внутренних точек. Множе- ство Мандельброта окружено иглоподобными разветвленны- ми антеннами. Если поместить параметр с на самый конец антенны, то получится множество Жюлиа подобной формы.
На рис. 10.7 показано множество для с — i; оно носит назва- ние дендрит. На рис. 10.8 представлено множество при выборе с из другой части антенны (вторичного множества Мандельброта). Рис. 10.7 Какова структура множеств Жюлиа, если взять значе- ние параметра с вне множества Мандельброта М? В этом случае множество Жюлиа перестает быть связным и рас- падается в облако точек, называемое пылью Фату. Два примера таких множеств показано на рис. 10.9. Рис. 10.9 Простейший способ приближенного построения напол- ненного множества Жюлиа основан непосредственно на определении 10.1 и состоит из следующих этапов:
1. Выбирается квадрат с центром в начале координат, и строится достаточно мелкая сетка. 2. Для каждого узла z0 сетки вычисляется P^zq), О k N, для достаточно большого N. 3. Если | Pcfc(2o)| > В при некотором к, то z0 Кс. 4. Если |Pcfc(zo)| В, 0 п N, то zq € Кс. При этом в качестве константы В можно взять В — = тах{|с|,2}. Действительно, если \z\ > В, то |Рс(-г)| И |г2 + с| И И и, следовательно, P?(z) —► оо. Содержание этого пункта в основном составляет мате- риал из книги X. Пайтгена и П. Рихтера (см. [16]). 10.3. N-фуркации системы. В отличие от действи- тельного случая с квадратичной функцией итерации (10.1), в комплексной плоскости при изменении комплексного пара- метра с порождают не только бифуркации, но п-фуркации (п — любое целое), когда А:-цикл становится неустойчивым и рождается /cn-цикл. Если z0 — точка /с-цикла, то условием n-фуркации является выполнение соотношения = ехр(2тггт/п), dzo что влечет равенство \dzk/dzo\ = 1. Числа т/п, где т/п — взаимно простые и т — = 1,2,..., п — 1, называются поворотными ч и с л а - м и n-фуркации системы. Если обозначить через Cj после- довательные значения параметра с, при которых происхо- дит n-фуркация системы (10.1) с поворотным числом $т/п,
то оказывается, что, как и в действительном случае, суще- ствует предел (комплексное число) Cj — Cj-i iim ------— = от/п, J->OO Cj+1 — Cj т. e. 6m/n — комплексный вариант универсальной константы Фейгенбаума S?. Аналогом уравнения удвоения является уравнение п-фуркации g(z) = ag(n\z/a), где д^ — п раз взятая суперпозиция функции д. Если нор- мировать функцию условием g(0) = 1, то для параметра а получаем формулу Подобным же образом определяются константы 6т/п, ат/п для других функций, например для полиномов про- извольного порядка и рациональных функций. Значения этих величин для некоторых полиномов вида P(z) = с + + zd (2 < d 5) приведены в работе К. Бриггса [25], а мно- жества Жюлиа для рациональных функций, возникающих при реализации метода Ньютона для кубических полиномов, представлены в монографии X. Пайтгена и П. Рихтера [16]. Известно, что для действительных полиномов Р,-, име- ющих только действительные корни, метод Ньютона P(zk\ г‘+1 =г - (10-4) p'(z*) сходится к какому-либо корню почти при всех начальных значениях R. Оказывается, что подобным свойством обла- дает метод Ньютона и для полинома P(z) — z3 — 1 (с дву- мя комплексными корнями), т. е., за исключением началь- ных точек из множества нулевой плоской меры (множества
Жюлиа), наблюдается сходимость к одному из корней поли- нома P(z) = z3 — 1 = 0. Интересно отметить, что граница области притяжения (т. е. множества начальных точек, для которых метод (10.4) сходится к одному из корней полинома) корней имеет фра- ктальную структуру. Это и не удивительно, поскольку ме- тод Ньютона (10.4) для полинома P(z) представляет собой не что иное, как дискретную комплексную динамическую си- стему с рациональной функцией перехода N(z). Фракталь- ные множества, которые возникают в нелинейных динами- ческих системах (в частности в дискретных ДС), называют динамическими фракталами, в отличие от кон- структивных (классических) фракталов, о которых пойдет речь в § 13. Упражнения 10.1. Найти численно область притяжения для каждого из корней Zi(i = 1,2,3) полинома P(z) = z3 — 1, т. е. найти множества A(zt) = {z е С : lim Nk(z) = Zi}, г = 1,2,3. k—>oo Изобразить эти множества графически. 10.2. Построить множества A(zi) для корней полино- ма P(z) = %3 + (А — 1)х — Л = 0 при различных значени- ях А. 10.3. Доказать теоретически и проверить численно, что параметры аи св процессах (6.3), (10.1) определяют в дей- ствительном случае одну и ту же универсальную констан- ту 6, которая введена формулой (5.2).
§11. Отображение подобия и самоподобие множеств Одним из наиболее отличительных свойств множеств Жюлиа является их самоподобие. Если, например, рассмот- реть под микроскопом некоторый фрагмент на границе мно- жества Кс, то нам предстанет картина, которая, во-первых, практически не зависит от того, в каком месте этот фрагмент расположен, а во-вторых, существенно не отличается от той, которую мы видели без микроскопа. Множество Мандель- брота также обладает (с некоторыми оговорками) этим уди- вительным свойством самоподобия, на его границу нанизано бесконечное число копий самого себя. Попытаемся перенести наши интуитивные представле- ния об этом свойстве на язык математических терминов. Определение 11.1. Пусть X — метрическое простран- ство. Отображение S: X —* X называется отображени- ем подобия, если для любых х, у G X p(S(x), S(y}) = = гр(х, у) для некоторого г > 0. Определение 11.2. Два множества метрического про- странства называются подобными, если взаимно-одно- значное соответствие между ними можно установить с по- мощью некоторого преобразования подобия S. ПРИМЕР 1. Пусть X — R2. Преобразования а) растяжения (гомотетии) рг: (х,у) —> (ж', у'), х' = гж, У' = гу, б) сдвига ть : х' = х - h, у' = у - Ь2,Ь = (Ьг, Ь2)г; в) вращения (поворота на угол в ) О: х' — х cos 0 — у sin 0 у' = xsmO + у cos
г) симметрии относительно оси OY: х' = — х, у' = +у, а также их суперпозиция (в любом порядке) являются преобразованием подобия. Справедлива следующая теорема о представлении отображения подобия [27]. Теорема 11.1. Пусть X — сепарабельное гильберто- во пространство. Для того чтобы S было преобразовани- ем подобия, необходимо и достаточно, чтобы имело ме- сто представление S = рггьО, где уг. Тъ — отображение растяжения и сдвига, а О — ортогональное отображение, т. е. О* = О-1. Доказательство . Необходимость. Пусть S — преобразование подо- бия. Обозначим g(x) = (S(x) — S(0))/r. Тогда, согласно опре- делению 11.1, имеем .. . .... нэд — эдп .. .. ПяИ-зМН = "г = Ik - у\\. в частности, д(0) = 0, ||р(ж)|| = ||ж||. Покажем, что отображение д сохраняет скалярное про- изведение. Действительно, {д(х),д(у)) = |(||ЭДН2 + 11ЭД112- ||ЭД -ЭДН2) = = |(1И12 + IMI2- Н*ЭД12) = Поэтому, если {еД — ортонормированный базис, то таковым будет и {ЭДД. Из соотношений ОО ОО 9(х) = 52<ЭД,ЭД))ЭД) = вг)ЭД) i=l i=l
следует, что д(х) — линейное отображение. Умножая послед- нее соотношение скалярно на у, находим ОО Х/> = 22<^еО^(ег),г/) = г=1 ОО ОО = = ^{х,а){еь д-\у)), г=1 г=1 откуда следует, что д* = д~\ Для д(х) введем новое обозначение О(х). Тогда о(х) = (ад - s(o))/r, S(x) = r(O(x) + ±S(0)) = дДтДОД)), где b = —S(0)/r. Достаточность. С учетом того факта, что ортого- нальное отображение сохраняет норму, требуемое свойство проверяется непосредственно. Определение 11.3. Множество Q метрического про- странства X называется самоподобным, если его мож- но разбить на произвольно малые части так, что каждая часть будет подобна целому множеству. Таким образом, для таких множеств при некотором раз- биении каждая часть подобна исходному множеству, сле- довательно, существует отображение подобия, их связы- вающее. Поэтому открывается возможность представления самоподобного множества К как конечного объединения его образов для некоторой совокупности функций подо- бия <jJi, i = 1, 2,..., т, т. е. т K = \Jut(K). (11.1) г=1
При дополнительном условии, что каждое из Wj являет- ся сжимающим отображением, оператор Т = U™ называ- ется оператором Хатчинсона. § 12. Топологическая и фрактальная размерности Самоподобие не является характеристическим свойст- вом множеств Жюлиа, Фату и Мандельброта. Отрезок и прямоугольник, например, очевидно являются самоподоб- ными фигурами. Однако следует отметить, что множества, подобные множеству Жюлиа, обладают некоторой нерегу- лярностью, или, как иногда говорят, геометрической хаотич- ностью. Чтобы охарактеризовать оба этих свойства — и са- моподобие и нерегулярность, а также выделить класс таких множеств с геометрической хаотичностью, нам понадобит- ся понятие фрактальной размерности множе- ства. Напомним сначала определение топологической размерности. Используемое здесь пространство X яв- ляется метрическим пространством со счетной базой. Из все- го многообразия подходов к определению размерности мы остановимся на таком, в котором пространство имеет раз- мерность п, если произвольно малые куски пространства, окружающие каждую точку, могут быть ограничены под- множествами размерности < п — 1. Этот метод определения является индуктивным, причем исходным пунктом индук- ции является принятие пустого множества в качестве (-За- мерного пространства. Определение 12.1. Пустое множество и только пустое множество имеет размерность —1.
Пространство X имеет размерность п (п > 0) в точ- ке р, если р обладает произвольно малыми окрестностями, границы которых имеют размерность С п — 1. X имеет раз- мерность п (обозначаем dim X п ), если X имеет раз- мерность < п в каждой своей точке. X имеет размерность п в точке р, если верно, что X имеет размерность п в р и неверно, что X имеет размер- ность п — 1 в р. X имеет размерность п, если X С п верно, a dim X < п — 1 неверно. X имеет размерность оо, если dim X п неверно для каждого п. Перейдем к определению размерности Хаусдорфа. Пусть А — подмножество метрического пространства, diam(a) = = sup{p(rr, у): х, у € А} и {U} — счетное открытое покрытие множества А, т. е. А С Пусть she — положительные числа. Определим функцию /i®(A) = inf{^2 diam(f4)s : А С [J U, diam(t4) < е}, 2=0 2=1 где нижняя грань берется по всевозможным покрытиям. По- скольку /i®(A) как функция s не убывает, то существует пре- дел, конечный или бесконечный: h8(A) = lim^( А). Хаусдорфом было установлено, что для любого множе- ства А (А С /?п) существует такое число £>н(А), что hUA} = / 00 для s < ' ' [0 для 8 > Dh(A). Определение 12.2. Число Ря(А) = inf{s: h,s(A) = = 0} = sup{s: hs(A) = оо} называется размерностью Хаусдорфа, или фрактальной размерностью, множества А.
Множества с нецелой размерностью Хаусдорфа приня- то называть фракталами. Чтобы охватить некоторые экзотические случаи (кривая Гильберта), когда множество имеет целую размерность Хаусдорфа, но сохраняет все при- знаки фрактальности, дадим более общее определение, сле- дуя Б. Мандельброту. Определение 12.3. Фракталом называется мно- жество, размерность Хаусдорфа которого строго больше его топологической размерности. Заметим, что для отрезка, прямоугольника и паралле- лепипеда эти размерности, очевидно, совпадают. Необходимо отметить, что из-за трудности вычисле- ния величины Dh обычно на практике используется другое, более удобное для вычисления определение фрактальной размерности (box-counting dimension) как величины Рь(А), определяемой соотношением А(А) = Ит (12-1) где Ne(A) — наименьшее число шаров (открытых или за- мкнутых) диаметром е, покрывающих множество А. Как по- казывают примеры, величины Djj(A') и Db(A) не всегда сов- падают. Например, если А — множество рациональных то- чек на отрезке [0.1], то Db(A) = 1, a Dh(A) = 0. Однако для многих классических фракталов эти величины совпада- ют, что делает использование Рь(А) в качестве фрактальной размерности оправданным. Упражнения 12.1. Вычислить топологическую и фрактальную раз- мерности множества точек с рациональными координатами в n-мерном кубе.
12.2. Вычислить Dh(A) и Db(A) для множества А = = {0,1/2,1/3,..., 1/п,...}. 12.3. Вычислить фрактальную размерность Dh(A) для канторова множества и кривой Кох из § 13. §13. Галерея классических фракталов Рассмотрим примеры так называемых конструк- тивных фракталов, которые строятся с помощью некоторой основы и фрагмента, повторяющегося при каж- дом уменьшении масштаба. Вычислим их фрактальную раз- мерность, используя формулу (12.1). 13.1. Канторово множе- ство. Оно строится следующим образом. Отрезок [0,1] делится на три равные части и удаляется средний интервал (1/3,2/3). К ос- тавшимся двум отрезкам при- меняется снова та же процеду- ра. Схематично это показано на рис. 13.1. Для покрытия n-й ста- дии построения канторова множе- ства требуется Ne(A) = 2п шаров (отрезков) диаметром в = 3~п. Поэтому £)В(Л) = lim = |. п—>ОС 1П о о о 1 Ч I---------1 1 2 1 3 3 Рис. 13.1 о 13.2. Кривая Кох. Здесь также отрезок [0,1] разби- вается на три равные части. Над средним интервалом стро- ится равносторонний треугольник, а основание его убирает- ся. Затем к каждому звену полученной ломаной линии при- меняется аналогичное построение и так далее.
Рис. 13.2 Четыре этапа (начиная с ну- левого) построения кривой Кох изображены на рис. 13.2. Ясно, что для покрытия кривой, по- лученной на п-м этапе, потре- буется 4П шаров (отрезков) диа- метром £ — (1/3)”, что дает Цв(Л) = lim п—>ОО 1п4га 1пЗп In 4 1пЗ’ Заме- тим, что имеются многочислен- ные другие варианты кривой Кох (см. [30] стр. 90-93). Рис. 13.3 13.3. Фрактал Мандель- брота-Гивена. Снова отре- зок [0, 1] делится на три части, и над средним отрезком выстра- ивается квадрат, две стороны ко- торого продолжаются вниз на ве- личину длины стороны квадра- та, как показано на рис. 13.3. Теперь к каждому отрезку дли- ны 1/3 полученной фигуры сно- ва проводятся те же построения. Два этапа описанной процедуры представлены на рис. 13.3. Ясно, что Db(A) = lim = -jaj. 13.4. Решето Серпинско- го (клиновидная кривая). В правильный треугольник встра- ивается равносторонний тре- угольник с вершинами, лежащи- ми на серединах сторон исход-
ного треугольника. Процедура повторяется для получен- ных треугольников, за исключением треугольника, лежа- щего в центре. Два этапа таких построений показаны на рис. 13.4. Объектом нашего анализа является кривая, со- Рис. 13.4 ставленная из сторон всех тре- угольников. Имеем Db(A) = = Ит = иг!’ ЧаСТ° рассматривают вариант реше- та, когда в качестве исходного объекта берется равнобедрен- ный прямоугольный треуголь- ник. 13.5. Ковер Серпинско- го (двумерное множество). В общем случае квадрат делит- ся на п2 равных частей (квад- ратов), и затем изымается к2 квадратов. На рис. 13.5 пред- ставлен случай п = 5, к = 3. К оставшимся квадратам снова применяется описанная процедура, т. е. каждый из них де- лится на п2 частей и изымаются к2 ячеек по выбранному правилу (через один). В этом случае Рь(А) = Ит \п= п—>ос т5" = in 16 In 5
Рис. 13.7 13.6. Фрактал Давида (дву- мерное множество). Исходным множеством является правильный ше- стиугольник. Каждая его сторона де- лится на три равные части, и про- водятся построения, как показано на рис. 13.6. В результате получается 6 ше- стиугольников меньшего размера. Для каждого из них процесс снова повторя- ется (заштрихованные части не учиты- ваются). Фрактальная размерность полученного в пределе множества со- ставляет величину А(Л) = Ит п-*<ю m3 m3 13.7. Пятиугольник Дюрера (двумерное множество). Сторо- на правильного пятиугольника делит- ся по правилу золотого сечения. По- лученный в результате этой процеду- ры меньший отрезок используется как сторона нового пятиугольника, примыкающего к вершине исходного пятиугольника (см. рис. 13.7). При этом еще один пятиугольник (шестой) выстраива- ется в центре. Для каждого из 6 полученных пятиугольников схема повторяется. Диаметр описанного (покрывающего) круга на каждом очередном этапе уменьшается в 2/(3 — -\/5) раз, следовательно, Д>(А) = lim — п—>оо In In 6 З + ч/б" 2
13.8. Губка Серпинского (трехмерное множе- ство). Единичный куб разбивается на кубические ячейки со стороной, равной 1/3, и изымается 7 кубиков (на рис. 13.8 они заштрихованы). Процесс повторяется для каждого из оставшихся кубиков. Поэтому для фрактальной размерно- сти предельного множества получаем значение Db(A) = lim ln20n 1пЗп In 20 1пЗ ’ Рис. 13.8 Список классических фракталов не ограничивается мно- жествами, приведенными выше. Упомянем некоторые из них: дерево Пифагора, кривая дракона, 3/2 -кривая, кривая Гильберта и др. Остановимся лишь на последней, посколь- ку она обладает одной особенностью, а именно имеет целую фрактальную размерность. 13.9. Кривая Гильберта. Квадрат делится на четы- ре квадрата, и середины полученных квадратов соединяются отрезками, как показано на первой фигуре рис. 13.9. Затем
такие же построения проделываются для каждого из мень- ших четырех квадратов. Причем на нижних двух квадра- тах фигура разворачивается на 90°, как показано на втором квадрате. Четыре фигуры соединяются между собой, как показано на 3-м квадрате. Снова те же построения проде- лываются для образованной фигуры, уменьшенной в 2 раза, и т. д. Простой подсчет показывает, что П,,А 1п22” — 1 о Р’(Л) = М 1„ 2" =2- Таким образом, кривая Гильберта имеет целую размер- ность Db = 2, хотя налицо все атрибуты фрактальности. Вот почему предпочтительней определение фрактала как множе- ства, для которого топологическая и фрактальная размерно- сти различны (см. определение 12.3). 13.10. Структура аттракторов. На рис. 13.10 изоб- ражены аттракторы (предельные точки) для критических значений параметра щ (г = 1,2, ...,5), т. е. значений, при которых происходит бифуркация в системе (5.1) с квадра-
тинной функцией. Можно видеть, что кластер аттракторов, отмеченный на самой нижней линии, повторяет в миниатю- ре совокупность аттракторов, отмеченную на линии выше, если ее отобразить симметрично относительно крайней пра- вой точки. При этом кластер на каждой последующей линии упакован примерно в а раз плотнее (см. формулу (6.1 о)) кластера на предыдущей линии, т. е. наблюдается само- подобие. Вычислена приближенно фрактальная размерность DB(A) для 2п-цикла в пределе при п—*оо (см. [23]), а именно РВ(Л) = -1п2/1п[|(± + ^)] ~ 0.543. Упражнения 13.1. Вывести формулу для длины n-й стадии построе- ния кривой Кох и клиновидной кривой (Рис. 13.4). 13.2. Построить фракталы, отличные от представлен- ных выше. Вычислить их фрактальные размерности. 0.0 0.25 0.5 0.75 1.0 Mi ।-----1----1—|—।-----1 М2 1---------1-------1—1 Мз I-------1---1-----1Д-Н
§ 14. Функциональное уравнение для фракталов Как было отмечено в § 11 для самоподобных множеств, каковыми являются фракталы, открывается принципиаль- ная возможность конструирования функционального ура- внения второго рода K = W(K\ (14.1) где оператор Хатчинсона W, действующий на множестве всех подмножеств некоторого исходного метрического про- странства X, имеет представление W — U”=1Wj, a Wi — сжи- мающее отображение подобия. Оказывается, что для каждого фрактала можно постро- ить отображения Wj таким образом, что данный фрактал бу- дет единственным решением уравнения (14.1), и это решение можно получить методом последовательных приближений кп+1 = W(Kn), начиная с произвольного начального приближения К° О (Х° — любое замкнутое множество). Появляется возмож- ность восстанавливать фракталы с любой точностью. 14.1. Уравнение для канторова множества. Обо- значим это множество через С. Введем отображения ?щ(т) = |, w2(x) = | + |, W(A) = W1(A) и w2(A). О ОО Покажем, что W(С) = С. Как известно, канторово мно- жество можно представить как совокупность точек, троич- ное разложение которых состоит только из 0 и 2, т. е. С — {-X : х = O.di<i2 • • •, € {0, 2}}.
Пусть х = O.aia,2 ... = O13-1 + а23~2 + • • • € С. Тогда Wi(x) = ^х = G13 2 + о23 $ Т • • • = 0.0aia2 ... С, О = 7?{O.6I16I2 • • •} + S — 0.0<Z-i<Z-2 • • • + 0.20... = о о = 0.2ai а2 ... С. Итак, РГ(С') С С. С другой стороны, если у = 0.Oia2 ... и у G С, то существует х €Е С и точно одна из двух функций иу (г = 1,2), что Wi(x) = у. Именно, положим х — О.а2оз .... Если у элемента у oi = 0, то выбираем wi, что дает ^{«2^ + + •••} — 0.0а2аз... = у. В противном случае, если ai = 2, то выбираем w2 и имеем ^{«2^ + +•••} + | = 0.2о2оз ... — у 14.2. Уравнение для кривой Кох. Прежде чем пе- рейти к построению образующих функций для кривой Кох, обратимся к рис. 13.2, где представлены четыре стадии этой кривой. Каждая стадия кривой не является, конечно, само- подобной фигурой в строгом смысле. Но можно заметить, что уже для четвертой стадии уменьшенная в 3 раза целая фигура мало отличается от каждой из образующих ее четы- рех частей. Ясно, что в пределе, т. е. для истинной кривой Кох, можно говорить об их совпадении. Поэтому левая часть может быть получена как умень- шенная в три раза копия целой кривой. Правая часть — сдвиг полученной копии по оси х на величину 2/3. Левая боковая часть строится поворотом копии на 60° и сдвигом
по оси х на 1/3. Правая боковая часть получается поворо- том на —60°, сдвигом по оси ж на 1/2 и по оси у на величи- ну \/3/2. Выражения для функций Wi(x, у), входящие в оператор Хатчинсона W из уравнения (14.1), принимают вид ^i(*^5?/) 3^)’ ^2(«£> ?/) (3х + З'З^’ ™з(я, У) = - ^У + + |у), 1 V 3 1 w4(x,y) = (qX + -g-!/+ 2’ 14.3. Уравнение для решета Серпинского. В ка- честве исходного объекта возьмем прямоугольный треуголь- ник с катетами единичной длины. На рис. 14.1 представлены (0,0) (1,0) Рис. 14.1 две стадии построения клиновид- ной кривой Р. Определим функ- ции Wi(x,y) (г = 1,2,3): wi(x,y) = (±х, ±у), w2(x, у) = W!(x, у) + (|, 0), ш3(х,у) = wi(x,y) -I- (0,|). По построению кривая Р может быть представлена в виде Р={(ж,у): ж=(0.аха2 •••),?/= = (0.М2...)}, где в двоичном разложении координат.bk не могут одно- временно быть единицы. Очевидно, что множество Р явля-
ется объединением трех множеств Pi (i = 1,2,3) : Pi = {(ж, у) : х = (0.0а2а3 (0.06263 • • •)}; Р2 = {(яг, у) : х = (0.1а2а3 ...),?/ = (0.06263 Рз = {(ж, у) : X = (0.0а2а3 ...),?/ = (0.16263...)}. Покажем, что шДР) С Pit Действительно, если z = = (0.0102 ..., 0.6ib2 ...) G Р, то Wi(z) = ^z = (0.0oio2..., 0.06i62 )€. Pi. и Обратно, пусть и = (0.0а2о3 ..., 0.06263...) G Pi, тогда, взяв z = (О.о2а3..., 0.6263,...), имеем wi(z) = и. Таким образом, установлено,что гиДР) = Pi. Аналогич- но рассматриваются случаи г = 2,3. Тем самым установлено, что для решета Серпинского оператор Хатчинсона имеет вид IV(A) = wi(A) U w2(A) U w3(A). § 15. Итерационная аппроксимация фракталов 15.1. Свойство сжимаемости оператора Хатчин- сона в метрике Хаусдорфа. Напомним определение ме- трики Хаусдорфа между двумя множествами некоторого ис- ходного метрического пространства X. Определим сначала полуотклонение между мно- жествами А п В как величину /3(А, В) = supр(х, В) = sup inf р(х, у). хеА хеА vG-B
Как видно из рис. 15.1, полуотклонение не обладает свойством симметричности. Полуотклонение может быть за- писано также в более наглядной эквивалентной форме /?(А, В) = inf{г: А С В£}, В£ = {z G X: р(х, В) < г}. Рис. 15.1 Определение 15.1. Пусть А и В — множества метри- ческого пространства X с метрикой р между его элементами. Величина h(A, В) = тах{/?(А, В), /3(В, А)} называется метрикой Хаусдорфа. Определение 15.2. Пусть F — некоторое отображение из X в X. Если для любых х,у 6 X p(F(x), F(y)) < ар(х, у), то говорят, что F удовлетворяет условию Липшица (обозна- чим этот класс через Lip а ). Лемма 15.1. Если F G Lipa, то h(F(A), F(B')') < < ah(A, В).
Доказательство. Действительно, исходя из определения метрики Хаус- дорфа, имеем /г(Г(А), F(B)) = max{/3(F(A), F(B)), /?(F(B), F(A))} = = max{sup inf p(F(x), F(?/)), sup inf p(F(x), F(?/))} < хеА У^в уев xeA max { о sup inf p(x, y), tv sup inf p(x, z/)} = ah(A, B). хеА У^в y^B xeA Лемма 15.2. Пусть U”=1B, — ограниченные множества. Справедливо соотношение h^Ai^Bi) < sup А(Л, Bi). (15.1) Доказательство . Пусть для определенности Mu"=i A, u?=1Bi) = /?(и?=1 А, и”=1 в€). Тогда имеем /i(U?=1 А, и”=1В,) = sup inf р(х, у) = = sup inf p(x, y) sup inf p(x, y) = xeAi0 y^Bjo xeAio yeBiQ = /3(Aio,Bio) MAo,^io) < sup MA, Bi). Замечание 15.1. Соотношение (15.1) справедливо и для счетной совокупности множеств.
Из установленных лемм вытекает Теорема 15.1. Если для оператора Хатчинсона W -- = U"=1Wj, каждое отображение Wi е Lip г^, то h(W(A), W(ВУ) rh(A, В), г = max г* 1<г<п т. е. W 6 Lip г относительно хаусдорфовой метрики. В частности, если w, — сжимающее отображение, т.е. Vi < 1, то оператор W — сжимающий в метрике Хаусдорфа. 15.2. Метод последовательных приближений. Сначала сформулируем хорошо известную теорему Банаха для сжимающих отображений. Теорема 15.2 (Принцип сжимающих отоб- ражений). Пусть X — полное метрическое простран- ство, S: X X — сжимающее отображение, т. е. для некоторого г < 1 p(S(x), Sty)) < гр(х, у) Ух, уеХ. Тогда: 1) существует единственное решение х* уравнения х = = S(i); 2) для любого х° € X метод последовательных прибли- жений xk+i — S(xk) сходится к решению х* при k —> оо; 3) справедлива оценка к ^ptAS^T). (15-2)
Доказательство . Для любых целых к, р имеем неравенства р(хк, хк+р) р(хк, xk+1) Ч- p(rrfe+1, хк+2') Ч-F +р(хк+р~\ хк+р) С гкр(х°, х1) Ч- тк+1р(х°, х1) Ч-I- 4-rfe+p-1/?(2?0, ж1) < rfc(l Ч- г Ч-1- гк+р~1)р(х°, х1) S^- Так как правая часть стремится к нулю при к —> оо, то последовательность фундаментальна. Поэтому хк схо- дится к некоторому элементу х* G X. Переходя к пределу при р —► оо, получаем оценку (15.2). Теорема 15.3. Пусть X — полное метрическое про- странство с метрикой р и W = U”=1Wj — оператор Хат- чинсона, где Wi : 2х —> 2х — сжимающие отображения с константой Г{ < 1. Тогда для любого ограниченного замкнутого множе- ства итерационный процесс Ak+1 = W(Ak), к = 0,1,..., (15.3) сходится к единственному компактному множеству К, которое является решением уравнения А = ИДА), причем справедлива оценка (15.2) с метрикой Хаусдорфа. Доказательство . Известно,что множество всех замкнутых подмножеств полного метрического пространства с метрикой Хаусдорфа образует полное метрическое пространство, в котором мно-
жество всех компактных подмножеств является замкнутым подмножеством. По теореме 15.1 W — сжимающее отображение в мет- рике Хаусдорфа. По теореме 15.2 для любого замкнутого множества А имеет место сходимость к некоторому множе- ству А*, т. е. lim h(Ak,A*) = 0, к—>оо причем справедлива оценка к h(Ak, Л*) < h(A°, W(A°Y), г = max n. 1 — Г 1<г<п Так как в качестве Л° можно взять компактное множество, а непрерывное отображение переводит компактное множе- ство в компактное, то Л* — компактно. Замечание 15.2. В предыдущем параграфе были получе- ны уравнения с оператором Хатчинсона для некоторых фракта- лов. Очевидно, что все входящие в них отображения wi — сжи- мающие, поэтому для аппроксимации фракталов (с любой точ- ностью) применим итерационный процесс (15.3). Обратим также внимание на то, что все гщ в упомянутых примерах являются афинными преобразованиями, которые в общем случае предста- вимы в форме Wi(x) = AiX + bi (г = 1,2,..., n), где Ai : X —+ X — линейный ограниченный оператор (в примерах это матрица), X — линейное нормированное пространство, bi — фиксированный элемент из X, а коэффициент сжатия совпа- дает с нормой оператора Л$, поскольку sunlb(A~My)ll f г — ьир х^у IF У\\ = НАН-
Таким образом, на трех примерах (множества Канто- ра, Кох, Серпинского) показано, что конструктивные фрак- талы допускают не только простое геометрическое построе- ние, но и численное восстановление с помощью рекуррентной формулы (15.3). При этом в качестве начального приближе- ния А0 может быть взято любое подмножество исходного пространства. § 16. Проблема сжатия информации 16.1. Общий формализм. Предположим, что име- ется некоторый объект х°, для которого требуется большой объем машинной памяти для его хранения. Из результатов предыдущего параграфа следует, что если известен сжимаю- щий оператор S, для которого х° является неподвижной точ- кой и для задания S нужен сравнительно небольшой объем памяти, то целесообразно хранить этот оператор, а объект х° восстанавливать в процессе счета методом последовательных приближений. Однако построение такого оператора для заданного гг° весьма непростая, а в некоторых случаях неразрешимая за- дача. В этой ситуации естественно изменить постановку за- дачи и попытаться построить близкий к упомянутому опе- ратор, неподвижная точка которого не совпадает с х°, но ап- проксимирует гг° с нужной точностью. Теоретическая возможность построения такого операто- ра вытекает из следующей теоремы. Теорема 16.1. Пусть X— метрическое пространство и S X —> X — сжимающий оператор с константой г < 1, который удовлетворяет условию р(х°, S(x°Y) < г(1 — г).
Тогда для любого у € X найдется такой номер N, что при k N p(Sk(y),x°) < г. Доказательство. I Используя неравенство треугольника и условия теоре- I мы, имеем Я Р<х°,3\уУ) < р(х°, 5\г“)) + р(5‘(А.5‘(»)) « I « р(х°, S(z“)) + p(S(x°), S2(i")) + • • • + ptS*-1^0-), S‘(x°))+ I + rp(Sk-\x°-), S^fs)) « P(X°, S(z“)(l + r + • • + r‘-‘)+ I + ^pfx0,y) Й + г*р(х°,у) < s + rtp{x°,y), откуда следует требуемая оценка. j Обратимся теперь к оценке (15.2). Если положить в ней к = 0, то приходим к неравенству Барнсли (16.1) которое показывает, что аппроксимация интересующего нас элемента х° неподвижной точкой х* оператора S, вообще говоря, тем лучше, чем меньше величина р(х°, S(x0)). | Таким образом, чтобы восстановить интересующий нас 1 элемент х° с точностью £ последовательностью итера- | ций Sk(y), необходимо: | 1) построить сжимающий оператор S, удовлетворяющий I неравенству р(х°, S(a:0)) < г(1 — г); т 2) для заданного начального приближения у найти но- 1 мер к, для которого выполнено неравенство |
16.2. Случай множеств. Из результатов предыду- щего пункта следует, что точность восстановления элемен- та х° фактически зависит от величины р(т°, <$,(т0)), поэтому желательно строить оператор S, который минимизирует эту величину для фиксированного элемента х° : min{p(x0,S'(T0)) : S € F}. (16.2) Если речь идет об аппроксимации множеств (фракта- лов), то подлежит минимизации функционал W где W(B°) = U”=1Wj(B°), В°— фиксированное множество. Обычно ограничиваются рассмотрением оператора W из класса афинных отображений, т. е. Wj(-) = АД-) + bi. Поскольку нахождение точного решения задачи мини- мизации весьма затруднительно, естественно ограничиться нахождением поправки ДИ7 = U”=1 [ДАД-) + Аб,], которая задает направление убывания функционала /(А), т. е. f(W + AAVE) < f(W\ (16.3) Достаточные условия, которым должна удовлетворять по- правка AVE оператора W, чтобы выполнялось неравен- ство (16.3), описаны, например, в книге В. И. Бердышева и Л. В.Петрак [5] (см. §3.2 главы 2). 16.3. Случай функции. Когда обсуждают пробле- му сжатия информации, то обычно имеют в виду циф- ровую информацию, касающуюся некоторого изображения. Математической моделью черно-белого изображения может служить функция z = f(x,y), заданная на прямоугольни- ке П € R2, значение которой выражает глубину (оттенок)
серого в точке (х, у). Функцию можно считать масштаби- рованной так, что 0 < f(x,y) С 1- Предположим также, что f(x, у)— элемент некоторого нормированного простран- ства: например, X = £2(П). Пусть {£>;}— некоторая совокупность подобластей, на- зываемых доменами, из исходного множества П. Предположим, что, наряду с выделенным набором об- ластей Dj; С П, прямоугольник П разбит на регионы Ri (например прямоугольники) таким образом, что П = U”=1/?i, A Rj) = 0 (i 7^ J), где ц — мера Лебега. Пусть для каждого i е {1,2, ...,п} указан номер j(г) и взаимно-однозначное отображение щ: Dj^ —» Ri (г = = 1,2,..., ri), которое задается формулой Vi-.t^q q= (gi,g2) € Ri, t = (^,i2) € D^, qk=ai1t1 + al2t2 + /3ik, к = 1,2, (16.4) где верхний индекс i = 1,2,..., п означает номер региона. Теперь определим отображение S на каждом из регио- нов Ri по формуле /-ад, ад|Я,(<?) = а</«1(9)) + (>, (16.5) (? Ri, % 1) 2, ... , Л.), где flj, Ь~ найденные из некоторых соображений коэффици- енты. Формулу (16.4) можно также интерпретировать следу- ющим образом. Оператор S имеет представление S(f) (<?) = = U”=1w,(/)(q), где Wi — отображения, действующие на функции /, которые определены на множестве Dj^ форму- лой (16.5).
Поскольку качество аппроксимации зависит от близости элементов f и S(f'), то естественно выбирать коэффициен- ты из условия = inf ||Ш - - /ШН12(Лц- (16-6) Из следующих соотношений IIW) - s(/)||,.2(n)2 = у |S(/)W - 5(/)(3)|= dq = п п г = Еа* / ~ = i=1 i = Е«2 / 1Ж-Ж12л^ /т-ж12^Е1«г|2л, i=1 Dj(i) П i=1 где Ц = агпа22 — аг12а21, получаем достаточное условие для сжимаемости оператора S: »=1 Если вместо £2 (И) привлечь пространство С(П) с нор- мой 11/11= sup{|/(i)|: ten}, то, поскольку имеем оценку IW) - -$(7)НС(П) = тгах N su₽ х(?)) ~ Ч?))! QE:Ri <тах|(ц|||/-7НС(П)’
условие сжимаемости базового оператора принимает вид max{|aj| : i = 1, 2,..., n} < 1. Наилучшие значения параметров ai,bi при заданных функциях определяются из решения задачи (16.6). За- дачу нахождения подходящего индекса j(z) для отображе- ния 1^(г)(<1,<г) : Dj —> Ri, задаваемого формулой (16.4), можно осуществить в два этапа: сначала для каждой па- ры (j, г), j = 1,2,..., I, i = 1,2,..., п, находится величина = inf „ in.f П[Ж - - Ь]|/гЛ2, Vji‘.Dj—+Ri a,b J а затем находится номер j(i), для которого реализуется min{£y : j = 1,2, Итак, для реализации описанной процедуры необходи- мо: 1) выбрать набор домен Dj(j = 1,2,...,/) в заданной области П; 2) разбить область П на конечное число регионов Ri (г = 1,2, ...,п); 3) для каждого региона Ri выбрать подходящий до- мен ; 4) определить отображение —> Ri вида (16.4); 5) построить отображение S : X —> X так, чтобы S было сжимающим в пространстве^ и неподвижная точка S была близка к исходной функции (изображению) f(x,y'). 16.4. Численный пример. Фрактальное сжатие чер- но-белого изображения «Лена» выступает часто как тесто- вый пример, на котором проверяется качество предлага- емой методики. Это изображение (см. рис. 16.1) состоит из 256 х 256 графических элементов (пикселей) с глубиной
Рис. 16.1 цвета 8 бит на точку (256 оттенков серого, поскольку 255 = = 1111111 = 1 + 2 + 22 + • • • + 27 ). Таким образом, исходное изображение требует 8 • 256- •256/8 = 65536 байт памяти (1 байт =8 бит). В качестве регионов были выбраны 1024 непересекаю- щихся квадратов 7?i, R2,..., Я1024 размером 8x8 точек. Кол- лекция доменов D — {РД состояла из всевозможных квад- ратов размером 16 х 16 точек. Их общее количество состав- ляет 241 х 241 = 58 081. В соответствии с описанной в предыдущем пункте процедурой для каждого региона находится домен Dj^) и определяются отображения (формула (16.4)) и парамет- ры ai,bi. Тем самым строится базовый оператор = <</«'(</)) + А
С помощью оператора S генерируется итерационный процесс fk+i = S(fk), где в качестве начального при- ближения /° выбиралось равномерно серое изображение. На рис. 16.2 представлены начальное приближение, 1, 2 и 10-я итерации [30]. Рис. 16.2 Поскольку отображение р, задается шестью параметра- ми агк1, агк2113к (к ~ 1,2), то нужно хранить восемь парамет- ров для каждого номера г {г = 1,2,..., 1024). Задание восьми параметров требует 31 бит памяти (см. [30]), поэтому всего нужно 1024 • 31/8 = 3968 байт. Коэффициент сжатия состав- ляет к = 65536/3968 = 16.5. Этот параграф написан на основе материала, почерп- нутого из монографии В. И. Бердышева и Л. В. Петрак [5] и книги X. Пейтгена, X. Юргена, Д. Соупа [30].
Глава 3 Непрерывные динамические системы С точки зрения математики в нелинейных динамических системах с числом степеней свободы больше 2 (особенно во многих биологических, метеорологических и экономических моделях) можно обнаружить хаос, и, следовательно, на достаточно больших временах их поведение становится непредсказуемым. Г. Шустер § 17. Модели динамических процессов Стандартными моделями динамических процессов с не- прерывным временем являются дифференциальные уравне- ния. Впервые появившиеся в задачах механики, они в на- стоящее время широко используются и в других областях при моделировании различных процессов. 17.1. Одномерная модель динамики популяции. Простейшей моделью динамики популяции с непрерывно ме- няющимся временем является дифференциальное уравнение х = ах — /Зх — 7т2. (17-1)
Здесь x(t) — численность популяции в момент времени t, а > 0 — коэффициент рождаемости, /3 > 0 — коэффициент смертности, вызванной старением организма, 7 > 0 — коэф- фициент смертности, связанный с ограниченностью ресурса. Для заданного гс(О) = xq (значения численности популяции в начальный момент времени t = 0) решение этого уравне- ния можно найти аналитически. При а 0 aeat при а = 0 x(i) = —Ц- 7* + *(<) = Здесь динамика численности определяется коэффициентом естественного прироста а = а — /3. При а > 0 численность популяции x(t), независимо от начального значения х$, при t —> 00 стремится к стационар- ному значению (см. рис.17.1) lim xit} — х t—oo % > 0. В случае а С 0 численность популяции стремится к ну- лю, популяция вымирает. 17.2. Модель «хищник —жертва». Классической моделью динамики двух взаимодействующих популяций яв- ляется модель Лотке - Вольтерра, задаваемая системой двух дифференциальных уравнений х = ах — (Зху, У = —1У + 6ху, а > 0, /3 0, 7 > б, 6 0. (17.2)
Рис. 17.1. Динамика решений уравнения (17.1) при а > в Здесь x(f) — численность жертв, y(f) — численность хищ- ников. В условиях изоляции (/? = 8 = 0) численность жертв x(f) = е°^Жо экспоненциально возрастает, а численность хищ- ников y(i) = экспоненциально убывает. При взаи- модействии (/? > 0, 6 > 0) в системе может наблюдать- ся равновесное состояние, отвечающее стационарному реше- нию x(f) = х = y(t) = у = Отклонение от равновесия ведет к колебательному решению (см. рис. 17.2). 17.3. Линейный осциллятор. Движение бруска, связанного с пружиной при наличии трения, задается диф- ференциальным уравнением второго порядка тх + кх + 1х = 0. Здесь x{t) — отклонение от положения равновесия, т — мас- са бруска, к — коэффициент трения, I — коэффициент жест- кости пружины. Будем далее считать, что т = 1.
Рис. 17.2. Периодические решения в модели «хищник - жертва» От уравнения перейдем к системе х = у у = —1х — ку. (17.3) В отсутствие трения (к = 0) брусок совершает гармониче- ские колебания по закону x(i) = Tq coswi + sinwi, w = VI, где До = д(0), Уо = д(0) — начальные положение и скорость, aw- частота колебаний (рис. 17.3, а). Слабое трение (0 < к < 2ш) ведет к затуханию коле- баний (рис. 17.3,6). Действительно, в этом случае решение имеет вид x(t) = I x0cosa?it + 2г/0 + хок 2u?i sinwit --t е 2 ,<Д1 = и, следовательно, lim^oo x(t) = 0. При сильном трении (к > 2ш) система стремится к по- ложению равновесия без колебаний (рис. 17.3, е).
Рис. 17.3. Линейный осциллятор: а — отсутствие трения {к = 0); б — слабое трение (0 < к < 2о>); в — сильное трение (к > 2о>) 17.4. Электронный осциллятор. Уравнение Ван- дер-Поля. Классической моделью электронного генерато- ра является предложенное Ван-дер-Полем дифференциаль- ное уравнение х + х = <5(1 — <г2)х. Перепишем его в виде системы х = у у = -х + <5(1 - х2)у. (17-4) При <5 > 0 траектории системы с ростом t стремятся к неко- торой замкнутой кривой — предельному циклу (рис. 17.4). В результате в системе поддерживаются колебания фик- сированной частоты и амплитуды (автоколебания). 17.5. Химический осциллятор (брюсселятор). Изучая неравновесные химические процессы, Тьюринг (1952), Пригожин и Лефевр (1976) предложили и исследо- вали систему х — 1 — (5 + 1)ж -Т ах2у, а > 0, Ь > 0 у - Ьх — ах2у.
Рис. 17.4. Предельный цикл уравнения Ван-дер-Поля при <5 = 1 Данная модель, получившая название брюсселятора, подоб но уравнению Ван-дер-Поля, имеет устойчивые периодиче ские режимы (рис. 17.5). Рис. 17.5. Предельный цикл брюсселятора (а = 0.2, b = 1.065)
17.6. Хаотический осциллятор. Модель Лоренца. В 1962 году Лоренц, специалист по физике атмосферы, пред- ложил простую модель тепловой конвекции х = а(у — х) у = гх — у — xz z = ху — bz. (17-6) Прямое компьютерное моделирование этой системы при а = о = 10, b = 2, г = 28 показало сложное нерегулярное поведе- ние ее решений (рис. 17.6). Рис. 17.6. Хаотический аттрактор системы Лоренца в проекции на плоскость xOz Система Лоренца стала классической моделью хаотиче- ской динамики. 17.7. Модель Ресслера. Следующая система, пред- ложенная Ресслером (1976), является хорошей моделью, де-
монстрирующей периодическое и хаотическое поведение х = ~(y + z) у = х + ay z = а + z(x — /л). (17.7) §18. Фазовый портрет системы дифференциальных уравнений и его свойства 18.1. Основные понятия. В предыдущем парагра- фе представлены модели, задаваемые одно-, двух- и трех- мерными системами дифференциальных уравнений. В об- щем п-мерном случае автономная система дифференциаль- ных уравнений задается следующим образом: Х1 = f1(x1,X2,...,Xn) •^2 == /2(^1, 2'2, • • • , Хп) Хп = /П(Х1,Х2,...,ХП). Ее векторная запись имеет вид х = /(ж), (18-1) где Х1 Х2 л*) = /1(Ж1, ... ,хп) /2(^1, ...,хп)
Пусть X с Rn — область определения функции f(x). Пред- полагается, что при любых х^ из X дифференциальное уравнение (18.1) имеет решение х = <£>(£, х^), определен- ное для всех t 0, с начальным условием х(0) = Подробное геометрическое изображение решения x(t) с по- мощью графика — интегральной кривой — требует (п + 1)-мерного пространства переменных t, Xi,..., хп. Если при п = 1 интегральные кривые располагаются на плоскости t,xi и их изображение не вызывает особых затруднений, то уже при п = 2 соответствующие интегральные кривые лежат в трехмерном пространстве t,Xi,X2, что резко усложняет их наглядное представление. Для сокращения размерности можно пожертвовать переменной t, оставив только так назы- ваемые фазовые переменные Xi,X2,..., хп, составляющие п- мерное фазовое пространство. При п = 2 фазовое пространство двумерно и называется фазовой плос- костью. Проекция интегральной кривой x(t) на фазовое пространство называется фазовой кривой, или фа- зовой траекторией. Множество фазовых кривых, отвечающих различным начальным данным, называется фазовым портретом системы. Во многих случаях фазовый портрет позволяет получить достаточно наглядное представление о динамике системы. В каждой точке х фазового пространства системы (18.1) вектор /(ж) есть вектор скорости движения системы вдоль фазовой кривой, проходящей через эту точку. Век- тор /(ж) указывает направление касательной к соответству- ющей фазовой кривой. Множество точек фазового простран- ства с указанными в них направлениями составляют поле направлений системы (18.1). Поле направлений позво- ляет построить, хотя бы приближенно, фазовый портрет ис- следуемой системы. Для этого линии, изображающие фазо- вые кривые, следует провести так, чтобы в каждой своей
точке они имели касательную, совпадающую с полем направ- лений. Определение 18.1. Решение x(t) системы (18.1) назы- вается устойчивым по Ляпунову, если Ve > 0 35 > 0 : VQ О Vz(0) ||х(0) — т^|| < 6 => ||x(t) — <p(t, ж^)|| < е. В противном случае решение x(f) называется неустой- чивым. Определение 18.2. Решение z(t) системы (18.1) на- зывается асимптотически устойчивым, если оно устойчиво по Ляпунову и 35 > 0 : ||ж(0) — < 5 => lim ||ж(<)~а/°))|| = 0. t—>+оо Определение 18.3. Решение x(t) системы (18.1) назы- вается экспоненциально устойчивым, если За > 0 ЗК > 0 35 > 0 : VO 0 Wo) ||rr(O) — < 5 => ||x(t) — ¥>(i,ж^)|| Ke а‘||ж(£) — т^||. Определение 18.4. Множество М С X называется инвариантом системы (18.1), если V .r(0) ЕМ Vt > 0 9?(i, х(0)) G М. I Если & М, то и во все последующие моменты време- 1 ни <p(t, х^) Е М. Простейшим примером инвариантного | множества является точка покоя. 1
Определение 18.5. Точка х Е X называется точ- кой п о к о я системы (18.1), если Vi 0 <p(t, х) = х. Если х — точка покоя, то f(x) — 0. Все точки покоя системы (18.1) находятся из решения системы /(*) = 0. (18.2) Другим примером инвариантного множества является цикл. Определение 18.6. Пусть £(£) является Т-периоди- ческим решением системы (18.1): £(£ + Т) — £(£). Множе- ство Г = {£(£)|0 ^ £ < Т} называется циклом. В фазовом пространстве цикл изображается в виде за- мкнутой кривой. Возьмем в качестве начальной произволь- ную точку цикла. Можно показать, что фазовая кривая со- ответствующего решения совпадает с циклом. Введем функцию p(x,Y) = inf ||х — г/||, задающую рас- y&Y стояние от фиксированной точки х до множества Y. Определение 18.7. Компактное инвариантное мно- жество М системы (18.1) называется устойчивым по Ляпунову, если справедливо следующее: Ve > 0 3<5 > 0 : V£ >0 Уж(0) р(х(0>, М) < 8 => /?(</?(£, ж^),М) < е. Определение 18.8. Компактное инвариантное множе- ство М системы (18.1) называется асимптотически устойчивым, если оно устойчиво по Ляпунову и Э<5 > 0 : р(х^\М}<8=> lim p(<p(t, М) = 0.
При этом множество U = {х G Х\ lim p(<p(t, х), М) = 0} на- зывается областью (бассейном) притяжения инвариантного множества М. Определение 18.9. Компактное инвариантное множе- ство М системы (18.1) называется экспоненциально устойчивым, если За > 0 ЭК > 0 36 > 0 : Vt 0 Vx(0) р(х^, М) < 6 => p(<p(t, т(0)), М) Ke~atp(x^, М). 18.2. Фазовые портреты линейных систем. Рас- смотрим двумерную (п = 2) линейную систему х = Ах, <2ц П12 <121 а22 Пусть Ai, Л2 — собственные числа, a hi, h% — линейно неза- висимые собственные векторы матрицы А. По этим данным общее решение системы записывается аналитически х = С1вЛ1(/г2 + С1вЛ2<Л2. Здесь возможны следующие случаи [17]: a) Ai, А2 — вещественные одного знака. Фазовый пор- трет — узел (рис. 18.1,а); б) Ai, А2 — вещественные разных знаков. Фазовый пор- трет — седло (рис. 18.1, б); в) А12 = ос ± г(3 — комплексно сопряженные (а^О). Фазовый портрет — фокус (рис. 18.1, в); г) Ai,2 = ±г/? — чисто мнимые. Фазовый портрет — центр (рис. 18.1, г). При ReAi^ < 0 движение вдоль фазовых траекторий идет в направлении точки покоя х = 0 й lim x(t) =0. Точка t—>oo
покоя х — 0 — асимптотически устойчива. Тогда говорят, что узел (фокус) является устойчивым. При ReAi 2 > 0 движение вдоль фазовых траекторий идет по направлению от точки покоя в бесконечность. Точка покоя х = 0 неустойчива. В этом случае говорят, что узел (фокус) является неустойчивым. В случае центра 7?еА1д = 0 движение происходит по за- мкнутым фазовым траекториям вокруг точки покоя. Точка покоя х = 0 устойчива (но не асимптотически). В случае седла всегда имеется направление, движение по которому идет от точки покоя в бесконечность. Поэтому здесь точка покоя х = 0 является неустойчивой. Приведенная здесь детальная классификация фазовых портретов получена на основе аналитического представле- ния для общего решения рассматриваемой линейной си- стемы. Построение точных аналитических решений для нелинейных систем возможно лишь в каких-то частных случаях. В общем случае при исследовании отдельных траекторий и построении фазовых портретов нелинейных дифференциальных уравнений используют численные ме- тоды. 18.3. Численные методы решения дифференци- альных уравнений. Разобьем временной отрезок [to, to+T] на N частей узлами t0 < ti < ... < ty = to+T с шагом h = = —: tm+i = tm + h. Пусть x(t) — решение задачи Коши системы дифференциальных уравнений с начальным усло- вием X = /(t,x), x(t0) = Х0. Обозначим через хт приближенное значение для неизвест- ного точного решения x(tm) в момент tm. Для расчета хт используют различные методы.
Рис. 18.1. Фазовые портреты линейной системы: а — узел, б — седло, в — фокус, г — центр Метод Эйлера. Расчет приближенных значений ведется по формуле Получаемые приближения имеют погрешности первого по- рядка: ||x(tm) - жто|| = O(h).
Метод Рунге-Кутта. Расчет приближенных зна- чений ведется по формулам: Хт+1 = хт ± ^(-Kl ± 2К2 ± 2Кз ± К4); Ki = hf(tm, хт), К2 = hf(tm + %,хт + К3 = hf(tm ± 4, хт 4- -у-), К4 = hf(tm + h,xm + К$). Получаемые по методу Рунге-Кутта приближения имеют по- грешности четвертого порядка ||x(im) — xm|| = О(Л4). Вывод формул, анализ их погрешностей, описание других числен- ных методов можно найти в книгах [4], [20]. У пражнения 18.1. Построить поле направлений и фазовые траекто- рии следующих систем: а) модель Лотке - Вольтерра (17.2) для а = /3 = 7 = 5 = = 1; а = 6 = 1, /3 = 7 = 4; б) линейный осциллятор (17.3) для I = 1 при к = 0; к = ±0.5; к = ±1; к = ±2; в) осциллятор Ван-дер-Поля (17.5) при 5 = 0, 5 = ±0.5, 6 = ±1, 8 = ±2; г) х ± sin х = 0 д) х ± х — ах3 = 0, а = —1, ±1. 18.2. Построить фазовые портреты линейных систем х = Зх у = 2х ± у ’ ( х = х \ У = 2х — у ' х = х + 3у у = — 6х ± 5у ’ х = — 2х — 5у у = 2х + 2у
18.3. Пусть s = ап + a22 (след матрицы Л), А = = ацагг — Oi2«2i (определитель А). В плоскости парамет- ров (s, А) изобразить зоны, соответствующие узлу (устой- чивому и неустойчивому), фокусу (устойчивому и неустой- чивому), седлу и центру. 18.4. Сравнить фазовые кривые системы (17.2), полу- ченные методами Эйлера и Рунге-Кутта при различных ша- гах h = 0.5; 0.1; 0.01; 0.001. 18.5. Используя методы Эйлера и Рунге-Кутта, постро- ить фазовые портреты систем из задачи 18.1. 18.6. Построить фазовые портреты следующих систем при различных значениях параметров: а) уравнение Хайрера х + х = ех — (х)3; б) гликолитический осциллятор ' х = 1 — ху < • ( у = РУ [х- 1 + ; q + у J в) осциллятор Ван-дер-Поля (мягкий режим) х + х = 5(1 — гг2)х; г) осциллятор Ван-дер-Поля (жесткий режим) х + х = 5(1 + ах2 — 6х4)±, а = 5, b = 0.5; д) уравнение Дуффинга х + х + (Зх3 = 0; е) модель «хищник-жертва» с ограниченностью ресур- са ' х = ж(1 — ау) — ух2 у = у(/Зх-1)
ж) модель «хищник-жертвах насыщением хищника х = т(1 - 7^) - а-^—у X I- JU у = -У + Ртг^У -L JU з) маятник с трением х + 2бх 4- sin х = 0; и) брюсселятор ' х — 1 — (Ъ + 1)т + ах2у, а > 0, Ь > 0 у = Ьх — ах2у ’ к к) модель Хопфа (мягкий режим) ' ±1 = у,Х\ — Х2 — Zl(x2 + Т2) Х2 = X! + ух2 - х2(х% + Х$) ’ к л) модель Хопфа (жесткий режим) ' it = Х\(у + 2х2 + 2®2 — (х2 + Т2)2) — х2 х2 = Х2(]Л + 2x1 + 2^2 _ (zi + xl)2) + Х1' 18.7. Написать программу динамики множества точек, равномерно распределенных в начальный момент времени в заданном прямоугольнике. Положение каждой точки в по- следующие моменты времени определяется численно по за- данной системе дифференциальных уравнений. При помо- щи этой программы исследовать динамику различных си- стем в прямом и обратном времени. К каким предельным множествам стягиваются эти точки? Как меняется характер сходимости и форма предельных множеств при изменении параметров динамической системы?
§19. Анализ нелинейной системы в окрестности точки покоя Если для построения фазового портрета нелинейной си- стемы в основном используются лишь численные методы, то при исследовании характерных особенностей вблизи точ- ки покоя возможен общий аналитический подход. Этот под- ход состоит в отыскании для исследуемой нелинейной систе- мы некоторой близкой линейной, с тем чтобы по результатам анализа последней можно было судить об основных чертах динамики исходной нелинейной системы. 19.1. Система первого приближения. Рассмот- рим способ построения линейной системы первого прибли- жения в окрестности точки покоя исходной нелинейной си- стемы. Для наглядности ограничимся сначала случаем дву- мерной (п = 2) системы Х1 = /1(^1, х2) х2 = /2(^1, ж2)‘ (19-1) Пусть х — (zi,t2) — точка покоя системы (19.1). Это озна- чает, что Л(ж1,ж2) = О /2(^1,ж2) = О' (19.2) Разложим в окрестности точки покоя х правые части систе- мы (19.1) — функции /1 и /2 — в ряды Тейлора: ОТ О т Х1 = /101, Ж2)+^(^1, Ж2)(Х1-Х1)+^|(Ж1, х2)(х2-х2)+.. . , от о т • Х2 = f2(xi, ж2) + (Ж1, ж2) (а?! - Xi) 4- (Ж1, х2) (х2 - х2) +...
Первый член в каждом из этих разложений, благодаря (19.2), равен нулю. Далее идут линейные члены, за ними — слага- емые более высоких порядков, которые вблизи точки покоя существенно меньше линейных. Отбрасывая эти малые слагаемые и делая замену Z1 = Xi - Ж1, Z2 = Х2 - Х2, dh(- dfa, «и = «12 = ^(^1,^2), dh(- - ч df2,_ _ . 021 = ^(Х1’Ж2), «22 = ^(Х1,Х2), получим линейную систему ( Z\ = «1121 + (Z12Z2 Z2 = 02121 + (l22Z2 Данная система получила название системы перво- го приближения для исходной нелинейной системы в окрестности точки покоя. Системы первого приближения играют весьма важную роль в исследовании нелинейных си- стем. Как правило (если не рассматривать особые вырожден- ные случаи), общий характер фазового портрета нелинейной системы вблизи точки покоя совпадает с фазовым портре- том соответствующей системы первого приближения. Тип фазового портрета системы первого приближения, благода- ря ее линейности, определяется достаточно просто — анали- тически (см. §18.2). Примеры анализа некоторых двумерных нелинейных систем по системам первого приближения дают- ся в §19.3. В случае общей n-мерной системы (18.1) в окрестности точки покоя х замена х = х + z (Ц2Ц—мало) и разложение
Тейлора приводят к равенствам х = x + z = f(x + z) = f(x)+ ^(x)z +... . Отбрасывая нелинейные члены с учетом равенств — = f(x) = 0, получим для малых отклонений z = х—х состоя- ния х от положения равновесия х линейную систему первого приближения 4 = Az, А = (19.3) С/Ju 19.2. Устойчивость точки покоя. Исследование ус- тойчивости начнем с линейной системы (19.3). Для этой си- стемы точкой покоя является вектор z — 0. Анализ общего решения системы (19.3) позволяет получить (см. об этом [9]) следующие критерии. Теорема 19.1. Для того, чтобы точка покоя z = 0 системы (19.3) была устойчивой по Ляпунову, необходимо и достаточно, чтобы все собственные значения Ai (г = = 1,...,п) матрицы А имели неположительные веще- ственные части: ReAj 0. При этом собственным зна- чениям, лежащим на мнимой оси (ReAj = 0), должны со- ответствовать клетки Жордана размерности единица. Теорема 19.2. Для того, чтобы точка покоя z = 0 системы (19.3) была экспоненциально устойчивой, необхо- димо и достаточно, чтобы все собственные значения A^i = = 1,...,п) матрицы А имели отрицательные веществен- ные части: Re Aj < 0. Отметим, что в случае линейных систем понятия асим- птотической и экспоненциальной устойчивости эквивалент- ны.
Если система (19.3) является экспоненциально устойчи- вой, то она остается таковой и при малых изменениях ее па- раметров. В случае, если система (19.3) просто устойчива по Ляпунову, отмеченное свойство уже не выполняется. Сколь угодно малые изменения параметров могут перевести соб- ственные значения матрицы системы, лежащие на мнимой оси, в правую часть комплексной полуплоскости, что сдела- ет систему уже неустойчивой. Сформулируем теперь критерий устойчивости ния х(1) = х, являющегося точкой покоя нелинейной мы (18.1). Теорема 19.3. Для того, чтобы точка покоя реше- систе- X си- стемы (18.1) была экспоненциально устойчивой, необходи- мо и достаточно, чтобы у системы первого приближения л л t-\ z — Az, A = д-(х), ox была экспоненциально устойчивой точка покоя z = 0. Как видим, исследование экспоненциальной устойчиво- сти точек покоя нелинейной системы сводится к выяснению знаков вещественных частей собственных значений матриц соответствующих систем первого приближения. 19.3. Примеры Пример 1. Рассмотрим модель «хищник-жертва». Для отыскания точек покоя системы (17.2) составим систему (ах — (Зху = 0 ( — 'УУ + $ХУ — 0 Получаем две точки покоя: Xi — 0, у\ — 0 — полное от- сутствие какой-либо жизни, Х2 = уъ — ~~ равновесное состояние.
Матрицы соответствующих систем первого приближе- ния имеют вид Ai а 0 0 —у ^2 0 да L - 0у ' 6 0 У матрицы Xj собственные значения Ai = а, А2 = —у — ве- щественны и имеют разные знаки. Фазовый портрет первой точки покоя — седло. У матрицы А2 собственные значения Л12 = ±у/ау i — чисто мнимые. Фазовый портрет второй точки покоя — центр. Как видим (см. рис. 17.2), фазовые портреты систем первого приближения в окрестности точек покоя соответствуют фазовому портрету исходной нелиней- ной системы. Пример 2. Рассмотрим уравнение Ван-дер-Поля. Для отыскания точек покоя составим систему У = 0 —х + 5(1 — х2\у = 0 Здесь единственной точкой покоя является х = 0, у = 0. Матрица соответствующей системы первого приближения равна и имеет характеристическое уравнение Л2 — 6Х +1 = 0. В за- висимости от значений параметра 6 возможны следующие случаи: 1) При о < — 2 собственные значения ----- вещественные и отрицательные. Фазовый портрет — устой- чивый узел.
2) При — 2 < 6 < 0 собственные значения становятся комплексно сопряженными с отрицательной вещественной . 5 ± i\/4 — 52 _ „ частью: Ai^ =-----5-----• Фазовый портрет — устойчивый фокус. 3) При 5 = 0 собственные значения А12 = ±г — чисто мнимые. Фазовый портрет — центр. 4) При 0 < 5 < 2 собственные значения становятся ком- плексно сопряженными с положительной вещественной ча- . 5 ± г\/4 — 52 _ „ стыо: Ai;2 =----2-----• Фазовыи портрет — неустойчивый фокус. г\ гг х о \ ± V52 — 4 5) При о > 2 собственные значения Ai^ = ----%----- вещественные, положительные. Фазовый портрет — неустой- чивый узел. Здесь по системе первого приближения при различ- ных значениях параметра удается достаточно точно предста- вить фазовый портрет исходной нелинейной системы вбли- зи точки покоя (рис. 21.9). Однако в области, удаленной от точки покоя, нелинейная система имеет предельный цикл (рис. 17.4) — важнейшую особенность, которую система пер- вого приближения уже никак не отражает. Исследованию предельных циклов посвящен §20 нашей книги. Пример 3. Проведем анализ устойчивости точки покоя брюсселятора. Из системы {1 — (Ь + 1)я? + ах2у = 0 Ьх — ах2у = 0 найдем единственную точку покоя х = 1, у = Матрица соответствующей системы первого приближения равна А = Г 1 ° ~ ~b —а
и имеет характеристическое уравнение А2 — (6—а— 1)А+аЬ = = 0. Для того, чтобы при положительных коэффициентах а и b выполнялось условие ДеАх^ < 0, необходимо и достаточ- но, чтобы выполнялось неравенство b < а + 1. Полученное неравенство задает в пространстве значений па- раметров а > 0, b > 0 область, внутри которой точка покоя брюсселятора сохраняет экспоненциальную устойчи- вость. Граница области устойчивости задается уравнени- ем Ь* = а + 1. При переходе параметров системы через эту границу точка покоя теряет устойчивость. Потеря устойчи- вости сопровождается появлением у системы качественно но- вого типа решения — предельного цикла (см. далее в § 21, рис. 21.8). Упражнения 19.1. Доказать, что для уравнения (17.1) а) решение х = при а > 0 является экспоненциально устойчивым; б) решение х = 0 при а = 0 асимптотически устойчиво, но не экспоненциально; в) решение х = 0 при а > 0 — неустойчиво. 19.2. Доказать, что необходимым и достаточным усло- вием асимптотической устойчивости точки покоя (0,0) ли- нейной двумерной системы являются неравенства s = Оц + + 022 < 0, △ = ®11<J22 ~ а12а21 > 0. 19.3. Провести анализ устойчивости точек покоя нели- нейных систем из задачи 18.6 (в пространстве их парамет- ров) по системам первого приближения.
19.4. Для следующих систем У - х3 . -X - у3 ’ х У У 2-т’ ( X = у + X3 Y у = х + у3 исследовать устойчивость точки покоя х = 0, у = 0. 19.5. Доказать критерии устойчивости и асимптотиче- ской устойчивости теорем 19.1 и 19.2 для случая п = 2. § 20. Анализ системы в окрестности цикла 20.1. Основные понятия. Система первого при- ближения. Рассмотрим систему i = /(.г), (20.1) имеющую Т-периодическое решение x = £(t) (£(t+T) = £(£)). Графиком такого решения будет замкнутая фазовая кри- вая — цикл Г. Точка х0 = £(0) отмечает на Г начальное поло- жение этого решения. Всякое решение, стартующее с любой другой точки цикла, будет двигаться по этой же замкнутой кривой. Цикл Г — инвариантное множество системы (20.1). Если начальную точку взять в окрестности цикла, то траек- тория соответствующего решения может вести себя различ- ным образом. Здесь возможны следующие варианты: а) решение приближается к циклу так, что отклонение от цикла стремится к нулю; другими словами, фазовая тра- ектория наматывается на цикл; б) решение движется вдоль цикла и формирует замкну- тую фазовую кривую — новый цикл, расположенный рядом с исходным; в) решение удаляется от цикла; фазовая кривая разма- тывается по спирали.
Рис. 20.1. Циклы: а — устойчивый, б— неустойчивый, в — полу- устойчивый Среди возможных сочетаний динамики снаружи и вну- три цикла обычно выделяют следующие: устойчивый пре- дельный цикл (рис. 20.1, а); неустойчивый цикл (рис. 20.1, б); полуустойчивый цикл (рис. 20.1, в). Нас интересуют условия, при которых предельный цикл имеет сильное экспоненциальное притяжение. Предполагается, что в U — некоторой окрестности кри- вой Г - определена функция = argmin^p ||j/ — ж||, за- дающая для каждого х из окрестности U ближайшую к ней точку 7(а?) с цикла Г. Тогда функция д(а;) = х—'у(х) задает отклонение точки х от цикла Г. Определение 20.1. Решение £(£) будем называть экспоненциально орбитально устойчи- в ы м, если найдутся такие а > 0, К > 0, что справедли- во неравенство IW))II для любого решения x(t) системы (20.1) с начальным усло- вием ж(0) = Xq Е U. Динамика малых отклонений, как и в случае точки по- коя, определяется системой первого приближения.
Рассмотрим для Т-периодического решения £(£) систе- мы (20.1) новую переменную z — х — £(£). Подставив х = = £(£) + г в систему (20.1) и разложив ее правую часть в ряд Тейлора, получим € + i = /« W + г) = /К(г)) + g(C(t))2 +... . Отбрасывая нелинейные члены с учетом тождества £ = = /(£(£)), получим линейную систему первого приближения z = F(t)z, F(t) = g(€(t)). (20-2) Матрица F(t) этой системы вслед за функцией £(£) являет- ся Т-периодической. 20.2. Линейные системы с периодическими ко- эффициентами. Элементы теории Флоке. Рассмот- рим для системы z = A(t)z, (20.3) где A(t) — произвольная Т-периодическая (п х п)-матрица, фундаментальную матрицу Z(t) = [^i(i) z2(£) ... zra(£)], со- ставленную из линейно независимых решений системы (20.3) с начальными условиями zi(0) = (1,0,... ,0)т, г2(0) = = (0,1,..., 0)т,..., гп(0) = (0,0,..., 1)т. Лемма 20.1. Для фундаментальной матрицы спра- ведливо тождество Z(t + T) = Z(t)- Z(T). (20.4) Доказательство. Рассмотрим вектор-функции ii(t) = zi(t + Т),..., zn(t) = = zn(t + Т) — столбцы матрицы Z(t + Т). Благодаря Т-пе-
риодичности матрицы A(t), справедливы соотношения Zi(t) = = A(t + T)zi(t+T) = A(t)zi(t+T) = A(t')zi(t'), означающие, что функции Zi(f) также являются решениями системы (20.3). Последнее позволяет связать их при помо- щи фундаментальной матрицы Z(t) со своими начальными значениями Zj(t) = Z(t)z(0). Переписывая эти соотношения в исходных обозначениях Zi(t + T) — Z(t)zi(T), получим тре- буемое тождество (20.4). Матрица В = Z(T), задающая отображение за период Т системы (20.3), называется матрицей монодромии. Любое решение z(t) системы (20.3) в моменты времени, крат- ные периоду, благодаря (20.4) легко выражается с помощью матрицы монодромии через начальные данные: z(kT) = Z(kT)z(0) = Z((k - l)T)Z(T)z(0) = ... = Bkz(0). Рассмотрим постоянную матрицу A = ^1пВ (В = егл) и матричную функцию Ф(£) = Z(t)e~tA. Лемма 20.2. Матричная функция Ф(£) является Т-пе- риодической. Доказательство следует из цепочки равенств: Ф(* + Т) = Z(t + T)e-(t+T)A = Z(t)Z(T)e-TXe-tK = = {Z(T)e~TA = 1} = Z(t)e"tA = Ф(£). Приводимость. В системе (20.3) сделаем замену пере- менных z — Ф(1)у, где у — новая переменная. Теорема 20.1. Система (20.3) в новых переменных у имеет вид У = Лу. (20.5)
Доказательство. Равенство (20.5) следует из (20.3) и соотношений z = Ф(б)у + Ф(£)£ = Z(t)e~tXy + Z(f)e-tA(-A)?/ + Ф(Д)у = = A(t)Z(t)e~tAy - Z(t)e~txAy + Ф^)у A(t}z = А(б)Ф(б)у = A(f)Z[t)e~tKy. Как видим, линейная система с периодическими коэффици- ентами (20.3) подходящей заменой переменных приводится к системе (20.5) с постоянной матрицей. Собственные значения pi (г = 1,...,п) матрицы моно- дромии В = Z(T) называют мультипликаторами системы (20.3). Собственные числа А; матрицы А — харак- теристические показатели — связаны с мульти- пликаторами соотношениями Aj Ln Pi, Pi e Критерием асимптотической устойчивости решения y=Q си- стемы (20.5) является условие Re Аг < 0 (см. теорему 19.2). В силу равенства |рг| = е^еХгТ это эквивалентно усло- вию |р,| < 1. Полученное неравенство и теорема 20.1 дают следующий результат. Теорема 20.2. Для асимптотической устойчивости решения z = 0 системы (20.3) необходимо и достаточ- но, чтобы все мультипликаторы удовлетворяли неравен- ствам <1, (г = 1,... ,п). 20.3. Экспоненциальная устойчивость цикла. Си- стема первого приближения (20.2) в общем классе линей- ных систем с периодическими коэффициентами (20.3) имеет
важную особенность. Действительно, вектор-функция r(t) = = в силу равенств r(£) = |^(£(£))£(f) и f(t) = Ж*)) является частным решением системы первого приближе- ния (20.2). Отсюда, в частности, следует равенство г(Т) — — Вг(0), или, с учетом Т-периодичности функции r(t), ра- венство г(0) = Вг(0). Как видим, вектор г(0) является собственным вектором матрицы В с соответствующим соб- ственным значением, равным единице. Таким образом, в слу- чае цикла матрица монодромии В обязательно имеет муль- типликатор pi = 1. Полученное означает (см.теорему 19.2), I что точка покоя z = 0 системы первого приближения (20.2) никогда не может быть асимптотически устойчивой. Следует • подчеркнуть, что для экспоненциальной устойчивости цик- ла этого и не требуется. Вопрос об экспоненциальной устой- чивости цикла Г решается в зависимости от расположения остальных мультипликаторов р2,..., рп. ; Теорема 20.3 (Андронова —Витта). Для экспонен- циальной орбитальной устойчивости решения £(£) систе- < мы (20.1) необходимо и достаточно, чтобы мультиплика- торы Р2,---,Рп удовлетворяли неравенствам |рД < 1, i = - — 2,..., п. Доказательство см. в монографии Б. П. Демидовича [9]. Из теоремы Виета и формулы Лиувилля следуют равен- ства Pi • Р2 • • • • • рп = detB = e fotrF^dt. Неравенство [ trF(t)dt < 0 (20.6) Jo по теореме Андронова — Витта является в общем п-мерном случае необходимым условием экспоненциальной устойчиво- сти цикла Г. f I 1
В двумерном случае справедливо равенство р2 = det В, что делает неравенство (20.6) не только необходимым, но одновременно и достаточным условием экспоненциальной устойчивости цикла (критерий Пуанкаре). В случае цикла на плоскости мультипликатор р2 имеет простой геометрический смысл, показывая при малых откло- нениях, во сколько раз траектория приближается к циклу за один оборот. Величина р2 может служить мерой устой- чивости предельного цикла к начальным возмущениям. Ма- лость р2 означает высокую степень устойчивости. При зна- чениях р2 < 1, но близких к единице, цикл устойчив слабо. При р2 > 1 цикл неустойчив. Случай р2 = 1 — критиче- ский. Здесь возможны различные варианты: цикл устойчив, но не экспоненциально; цикл полуустойчив; цикл находится в окружении других близких циклов и т. д. Упражнения 20.1. Написать программу отыскания устойчивого пре- дельного цикла с заданной точностью. При помощи этой программы найти предельные циклы и соответствующие им мультипликаторы для следующих систем: а) уравнение Хайрера, б) гликолитический осциллятор, в) осциллятор Ван-дер-Поля (мягкий режим), г) осциллятор Ван-дер-Поля (жесткий режим), д) брюсселятор. 20.2. Исследовать изменения циклов из задачи 20.1 и их мультипликаторов при изменении параметров порождаю- щих их систем. Получить графики зависимости мультипли- каторов от параметров.
§21. Бифуркации 21.1. Структурная устойчивость и бифуркации. В исследовании зависимости фазового портрета системы от изменения входящих в нее параметров выделяют два слу- чая. К первому относят системы, изменение параметров ко- торых в некоторой области сопровождается сохранением ка- чественной картины их фазовых портретов. Два фазовых портрета называют качественно одинаковыми, если суще- ствует взаимно-однозначное и взаимно-непрерывное отобра- жение, переводящее один фазовый портрет в другой. Об- ласть параметров, внутри которой сохраняется качественная картина фазовых портретов системы, называется обла- стью структурной устойчивости. Второй слу- чай составляют системы, в которых при прохождении пара- метра у через некоторое значение у* происходит качествен- ное изменение фазового портрета (см.[1]). При этом говорят, что в системе при у = у* происходит бифуркация, а са- мо значение у* называют точкой бифуркации. Пример. Фазовые портреты линейной системы ( хг = yx-i | х2 = —х2 для различных у представлены на рис. 21.1. Как видим, система имеет две области структурной устойчивости. При —оо < у < 0 фазовый портрет — устой- чивый узел (рис. 21.1, а; 21.1,6). При 0 < у < оо фазовый портрет — седло (рис. 21.1, г). Эти области разделены един- ственной точкой бифуркации д* = 0 (рис. 21.1, е). Переход параметра у через бифуркационное значе- ние у* = 0 слева направо сопровождается потерей устой- чивости ТОЧКИ ПОКОЯ Х1 =0, х2 = 0.
Рис. 21.1. Фазовые портреты: а — при р = —0.5; б— при р = —0.1; в — при р = 0; г — при р = 1 Бифуркация потери устойчивости точки покоя в случае линейных систем обязательно приводит к тому, что фазовые траектории уходят в бесконечность. В нелинейном случае та- кая бифуркация может сопровождаться появлением у систе- мы новых инвариантных множеств.
21.2. Бифуркация рождения цикла. Рассмотрим двумерную систему — модель Хопфа Xi = рХ\ — х2 — £l(Xi + х?) Х2 = Хг + рх2 - Т2(Х1 + х%) (21-1) Перейдя от декартовых координат хг, х2 к полярным г, tp по формулам Xi = г cos </>, х2 = г sin </?, получим совсем простую систему . Q Г = fJLr — Г ф = 1 (21-2) состоящую из двух независимых уравнений. Решением вто- рого уравнения будет функция <р — po+t, где </?о — начальное положение угла Решения первого уравнения в зависимости от парамет- ра р имеют следующий вид: а) при // < О го = 0 — единственная точка покоя, к которой экспоненциально стремятся все другие решения при t —> -Ьоо (рис. 21.2,а); б) при р = 0 общая картина предыдущего случая со- храняется, однако скорость стремления решений к го = О перестает быть экспоненциальной: r(t) = г(0) ^/2tr2(0) + 1 (как видим, для всех р 0 точка покоя го = 0 является асимптотически устойчивой); в) при р > 0 у первого уравнения системы (21.2) наряду с го = 0 появляется еще одна точка покоя и = ^/р, к кото- рой стремятся экспоненциально другие решения. Располо- жение интегральных кривых для этого случая представлено на рис. 21.3, а.
Рис. 21.2. Динамика систем при у, = —1: а — система (21.2); б — система (21.1) Как видно, при переходе параметра у, через бифуркаци- онное значение ц* = 0 в область ц > О точка покоя г0 = О теряет устойчивость. При этом появляется новая устойчивая точка покоя п. Проведенный анализ позволяет проследить изменение фазовых портретов системы (21.1): а) при ц < О точка покоя (0,0) системы экспоненциально устойчива. Фазовый портрет — фокус (рис. 21.2, б); б) при ц = 0 точка покоя (0,0) — асимптотически устой- чива; в) при ц > 0 точка покоя (0,0) — неустойчива. В систе- ме появляется устойчивый предельный цикл — окружность радиуса ri = y/Jl (рис. 21.3, б). Таким образом, переход параметра ц через бифуркаци- онное значение ц* — 0 сопровождается качественным изме- нением фазового портрета системы (21.1). При этом потеря устойчивости точки покоя (0,0) сопровождается рождени- ем устойчивого предельного цикла — окружности х% + х% = = ц. Следует отметить, что размер цикла (радиус окруж- ности ri = уД) непрерывно меняется (возрастает) по ме-
Рис. 21.3. Динамика систем при у = 0.5: а — система (21.2); б — система (21.1) ре удаления параметра от своего бифуркационного значе- ния д* = 0. Такой вариант бифуркации называется мяг- ким рождением цикла, или бифуркацией Андроно- ва-Хопфа. Другой вариант бифуркации рождения цикла рассмот- рим на примере системы ±i = х^у, + 2х2 + 2^2 — (т2 + Х2)2) — х2 ±2 = Х2(у + 2х% + 2х% - (х% + Х%)2) + Хг ' В полярных координатах данная система распадается на два независимых уравнения ( г = г(д + 2г2 — г4) 1^ = 1 (21-4) Решения первого уравнения в зависимости от параметра у имеют следующий вид:
а) При ц, < — 1 единственной точкой покоя является Го — = 0. Все другие решения стремятся к го = 0 экспоненциаль- но. Точка покоя г0 = 0 устойчива (рис. 21.4, а). Рис. 21.4. Динамика систем при д = —1.2: а — система (21.4); б — система (21.3) б) При /л = — 1 наряду с Го = 0 появляется еще од- на точка покоя и = 1. Расположение интегральных кривых становится другим (см. рис. 21.5,а). Решения с начальными точками, лежащими выше ri = 1, монотонно убывая, стре- мятся к и = 1. Решения с начальными точками, лежащими ниже ri = 1, уже стремятся к го = 0. Точка покоя г0 — 0 устойчива, а точка покоя и = 1 полуустойчива. в) При — 1 < ц < 0 имеем уже три точки покоя: то = = 0, и = \/1 — \/1 + г2 = \/1 + д/1 + Д- Точка покоя гг при переходе д в интервал —1 < ц < 0 расщепляется в две. Соот- ветственно меняется и расположение интегральных кривых (см. рис. 21.6, а). Решения с начальными точками, лежащими в интерва- ле (ri, оо), стремятся к гг. Решения с начальными точками,
Рис. 21.5. Динамика систем при /z = —1: а — система (21.4); б — система (21.3) Рис. 21.6. Динамика систем при fj, = —0.5: а — система (21.4); б— система (21.3) лежащими в интервале (0, гД, стремятся к г\. Точки покоя г0 и г2 — устойчивы, точка покоя ri — неустойчива. г) При д = 0 имеем две точки покоя: tq = 0, Гг = Вид интегральных кривых показан на рис. 21.7,а. Точка по- коя г0 — неустойчива, точка покоя и — устойчива.
Рис. 21.7. Динамика систем при р, = 0 : а — система (21.4); б — система (21.3) д) При ц, > О имеем две точки покоя г0 = 0, г2 = = -(/1 + Общая картина подобна случаю г. Проследим, как представленные изменения в поведении интегральных кривых r(t) скажутся на изменении фазовых портретов системы (21.3). а) При д < —1 у системы (21.3) — единственная точка покоя (0,0). Фазовый портрет — устойчивый фокус (рис. 21.4, б). Все траектории стремятся к точке покоя. б) При ц. = — 1 у системы (21.3) появляется периоди- ческое решение (цикл), фазовая траектория которого есть окружность + $2 = 1 единичного радиуса. Фазовые траек- тории с началом вне единичного круга асимптотически стре- мятся к этому циклу. Все остальные с началом внутри круга стремятся к точке покоя (0,0) (рис. 21.5,6). В системе (21.3) при ц, = — 1 произошла жест- кая бифуркация рождения цикла (полуустой- чивого) фиксированного радиуса.
в) При — 1 < jtz < 0 на фазовом портрете мы видим результат расщепления полуустойчивого цикла из б на два новых цикла. Один из них (внешний, радиусом гд) — устой- чивый, другой (внутренний, радиусом ri) — неустойчивый (рис. 21.6, б). Точка покоя (0,0) по-прежнему остается устойчивой. При приближении ц к 0 внутренний цикл, уменьшая свой радиус, стремится к точке покоя (0,0). г) При у = 0 происходит слияние внутреннего цикла с точкой покоя (0,0). В результате данной бифуркации точ- ка покоя теряет устойчивость и у системы на фазовой плос- кости остается единственное притягивающее множество — предельный цикл (рис. 21.7,6). д) При // > 0 и дальнейшем увеличении параметра ц фазовый портрет системы (21.3) качественно не изменяет- ся. В нелинейных системах потеря устойчивости точки по- коя при изменении параметров часто сопровождается ро- ждением устойчивого предельного цикла. Отмеченная связь (потеря устойчивости точки покоя — появление предельно- го цикла) позволяет искать бифуркации рождения цикла по бифуркациям потери устойчивости точек покоя, используя для отыскания последних анализ линейных систем первого приближения. Возможности такого подхода иллюстрируются на при- мерах брюсселятора и уравнения Ван-дер-Поля. У брюсселя- тора (см. анализ системы первого приближения в примере 3 из §19) бифуркация потери устойчивости точки покоя х = Г = 1, у = - происходит при пересечении прямой b — а + 1. Например, для а = 1 бифуркационным значением параметра будет Ь* = 2. При b > 2 у брюсселятора вокруг уже неустой- чивой точки покоя появляется (рис. 21.8) предельный цикл.
Рис. 21.8. Фазовые портреты брюсселятора при а = 1: а — для Ь = 1.5; б — для b = 2; в — для b = 2.05; г — для b = 2.3 Размеры цикла по мере увеличения параметра Ь непре- рывно возрастают. Рис. 21.8 демонстрирует для брюсселято- ра бифуркацию мягкого рождения цикла. У модели Ван-дер-Поля (см. анализ системы перво- го приближения в примере 2 из § 19) потеря устойчиво- сти точки покоя х = 0, у = 0 (см. рис. 21.9) происходит при <5, = 0.
Рис. 21.9. Фазовые портреты модели Ван-дер-Поля: а — для 6 = = —0.2; б — для 6 = 0; в — для <5 = 0.2; г — для <5=1 Для <5 = 0 уравнение Ван-дер-Поля вырождается в ли- нейное с фазовым портретом типа центр. Фазовыми траек- ториями здесь являются окружности различных радиусов. Последнее означает возможность в данной системе гармони- ческих колебаний любой амплитуды в зависимости от выбо- ра начальных данных. При переходе в область S > 0 на фазовом портрете оста- ется единственная замкнутая фазовая кривая (цикл), близ-
кая при малых 6 к окружности радиуса 2. Остальные фа- зовые кривые асимптотически, при t —» +оо, приближают- ся к этому циклу. Фазовые портреты уравнения Ван-дер- Поля при 8 = —0.2, 5 = 0, 8 = 0.2, 5=1 представлены на рис. 21.9. Наличие у системы устойчивого периодического реше- ния с фиксированной частотой и амплитудой, не завися- щей от выбора начальных данных, позволяет использовать это электронное устройство в качестве генератора электри- ческих колебаний. Следует отметить, что подобные режи- мы, получившие название автоколебаний, возможны лишь в нелинейных системах. 21.3. Порядок и хаос в модели Лоренца. Отме- тим, что, наряду с решением x(t), y(t),z(t), для модели Ло- ренца (17.6) всегда будет решением и набор функций — x(t), —y(t'),z(t'). Поэтому фазовый портрет симметричен относи- тельно оси OZ. Рассмотрим для функции v = х2 + у2 + (z — ст — г)2 производную в силу системы Лоренца / \ 2 =-2w + |(<т-Ь г)2, w = сгх2 + у2 + b ( z — + r I . При достаточно больших К ^К2 шах сфе- ра v = /С2(сг+г)2 целиком содержит эллипсоид w = ^(tr+r)2. При этом во всех точках сферы выполняется неравен- ство i> < 0, означающее, что все фазовые трактории пере- секают сферу снаружи вовнутрь и далее из нее не выходят. £ Для любого решения след матрицы F = — имеет по- ох стоянное отрицательное значение trF = —(ст + b + 1) < 0,
что означает равномерное сжатие фазового объема (см. [5]). Таким образом, всякое притягивающее множество системы Лоренца имеет нулевой объем. Проследим изменение фазового портрета в зависимости от параметра г. При всех значениях параметров система Лоренца име- ет точку покоя х = 0, у = 0, z = 0. Характеристическое уравнение соответствующей системы первого приближения (Л + 6)(А2 + (<т + 1)А + сг( 1 — г)) = 0 имеет вещественные корни - (а + 1) ± у/(сг- 1)2 + 4<7Г А1 = — О, Л2.3 — ---------2-------------• При 0 < г < 1 все корни отрицательны и точка покоя х = — у = z асимптотически устойчива. При г > 1 эта точка становится неустойчивой. В этом случае у системы появляются еще две точки покоя с коор- динатами х = у = ±-\/6(г — 1) и z = г — 1. В силу симметрии они имеют одинаковый тип. Соответствующее характеристи- ческое уравнение Л3 + (а + b + 1)Л2 + Ь(а + г)Л 4- 2ba(r — 1) = 0 в силу критерия Рауса-Гурвица приводит к условию устой- чивости сг(а + Ь + 3) г < г — -----------. a — b— 1 При г > г* все три точки покоя становятся неустойчивы- О ми. Для рассматриваемых здесь параметров и = 10, b = - бифуркационное значение г* = 24.737.
При г = 28 Лоренц обнаружил притягивающее множе- ство — странный аттрактор (см. рис. 17.6). В модели Лоренца при дальнейшем увеличении пара- метра можно увидеть большое разнообразие как хаотиче- ских, так и периодических режимов. Так, например, известно окно периодичности — интер- вал 99.524 < г < 100.795. Этот интервал разбивается на подынтервалы Ц = (99.98,100.795), /2 = (99.629,99.98), Ц = — (99.547,99.629), ..., /2п, ... с предельными циклами Г1, Г2, Г4,..., Г2«,... Здесь — цикл, наблюдаемый на подынтервале Д. Переход параметра г из одного интервала в другой со- провождается бифуркациями удвоения периода. Так, напри- мер, при переходе из Ц в 12 цикл Г1 расщепляется и образует- ся 2-цикл Г2, при переходе из 12 в Ц 2-цикл Г2 расщепляется в 4-цикл Г4 и т. д. (рис. 21.10). Наблюдаемые здесь бифуркации удвоения периода сопровождаются характерными изменениями мультиплика- торов. Один из мультипликаторов, отражающий динамику возмущений вдоль цикла, всегда равен единице. Степень притяжения возмущенных траекторий к циклу характери- зуется двумя оставшимися: р\ и р2 (|pi| |р2|). Неравенство |pi)2| < 1 является необходимым и достаточным условием экспоненциальной устойчивости цикла. В стандартном сце- нарии бифуркации удвоения периода старший мультиплика- тор pi меняется следующим образом. При переходе парамет- ра г через точку бифуркации (границу между соседними ин- тервалами Д и I2k) значение pi стремится к —1, достигает — 1 и после появления цикла с удвоенным периодом становится равным 1. При этом модуль мультипликатора р2 остается строго меньше единицы. На интервалах Л, 12,..., /2«,... выделим значения и, г2,..., г2п,... параметра г, соответствующие минимуму фун-
Рис. 21.10. Бифуркации удвоения периода: а — 1-цикл для г = — 100.563; б — 2-цикл для г = 99.803; в — 4-цикл для г = 99.5866 кции |pi(r)|: rk = argmin |pi (г) |. reik Предельный цикл для значения г = гк является наибо- лее устойчивым в классе /с-циклов на интервале 1к- Такой цикл будем называть суперциклом. Для модели Ло- ренца ri = 100.563, г-2 = 99.803, г4 = 99.5866. Суперциклы, отвечающие этим параметрам, представлены на рис. 21.10. На рис. 21.11 дан график функции рДг). Здесь хорошо виден бифуркационный механизм. На каждом из интерва- лов 1г, I2, Ц, при убывании параметра г, функция pi(r) монотонно убывает от +1 до —1, имея скачки в точках би- фуркации. Что касается мультипликатора /?2(т), то он прак- тически не изменяется, сохраняя в данном диапазоне значе- ние |pi(r)| = 1.75 • 10-2.
В поведении мультипликаторов необходимо отметить одну важную деталь. Функция pi(r), непрерывно изменяя свои значения на интервале Д от +1 до —1, не может при- нимать значение, равное нулю (матрица монодромии явля- ется всегда невырожденной). Обход нуля связан с выходом в комплексную плоскость. Рис. 21.11. График мультипликатора pi(r) в окне периодичности Функция Р1(г) (совместно с р2(^) = /h(r)) становится комплекснозначной только в малой окрестности точки В остальной части интервала Д она имеет вещественные зна- чения. Символически эта деталь в поведении рДт) изображе- на на рис. 21.11 маленькими окружностями. После каскада бифуркаций удвоения периода в системе наблюдается хаос. Поведение решения в проекции на плоскость xOz представ- лено на рис. 21.12.
Рис. 21.12. Хаотический аттрактор системы Лоренца при <т = = 10, Ъ = |, г = 99 О У пражнения 21.1. Найти интервалы структурной устойчивости и ис- следовать бифуркации в системах ( х = У [ у = —х + ау' = х — ау = ах — у ‘ 21.2. Для системы «хищник-жертва»(с насыщением) в плоскости параметров а и b изобразить зоны структур- ной устойчивости. Выделить зону существования предель- ного цикла. Проиллюстрировать бифуркации рождения пре- дельного цикла при пересечении границ этой зоны. 21.3. Для модели Лоренца (<т = 10, b = 8/3) при раз- личных значениях г G (0, 28) найти решения O(t), ?/(£), Д*))
и выходящие из близких точек z(0) = 1, у(0) = О, ДО) = 0 и ДО) = 1 + 10-6, z/(0) = О, ДО) = 0. Изобразить x(t) и x(t) на графике. Найти момент времени Т, при котором ||ДТ)—ДТ)|| > 10. Исследовать зависимость Т от г. 21.4. Проследить бифуркации модели Лоренца для а = = 10, b = 8/3 при изменении г на интервалах а) (0,28); б)(99.5,100.7). 21.5. Исследовать зависимость мультипликаторов моде- ли Лоренца для ст = 10, b = 8/3 от параметра г на интервале (99.5,100.7). 21.6. Для модели Ресслера (17.7) для а = 0.2 проследить бифуркации при д € (2, 5).
Список литературы [1] Андронов А. А. Бифуркации динамических систем. — М.: Наука, 1962. [2] Анищенко В. С. Знакомство с нелинейной динами- кой: Учеб, пособие. — Саратов: Изд-во Гос.УНЦ «Колледж». 2000. [3] А р н о л ь д В. И. Теория катастроф. М.: Наука, 1990. [4] Бахвалов Н. С., Жидков Н. П., Кобель- ков Г. Менные методы. — М.: Наука, 1987. [5] Б е р д ы ш е в В. И., П е т р а к Л. В. Аппроксимация функций, сжатие численной информации, приложения. — Екатеринбург: УрО РАН, 1999. [6] Берже П., П о м о И., Видаль К. Порядок в хао- се. О детерминистическом подходе к турбулентности. — М.: Мир, 1991. [7] В у л Е. В., С и н а й Я. Г., X а н и н К. М. Универ- сальность Фейгенбаума и термодинамический формализм // Успехи физических наук. — 1984. — Т. 39, № 9. — С. 3-37. [8] Г л а с с Л., М э к и М. От часов к хаосу. Ритмы жизни. — М.: Мир, 1991. [9] Д е м и д о в и ч Б. П. Лекции по математической теории устойчивости. — М.: Наука, 1967.
[10] Жигулев В. М. Динамика неустойчивостей (динансти- ка). — М.: Изд-во МФТИ, 1996. [11] Кузнецов А. П. Колебания, катастрофы, бифуркации, хаос. — Саратов: Изд-во Гос. УНЦ «Колледж», 2000. [12] Кузнецов С. П. Динамический хаос. — М.: Наука. Гл. ред. физ.-мат. лит., 2001. (Сер. Современная теория колеба- ний и волн). [13] Лихтенберг А., Либерман М. Регулярная и хаотическая динамика. — М.: Мир, 1984. [14] Малинецкий Г. Г., Потапов А. Б. Современные проблемы нелинейной динамики. — М.: Эдиториал УРСС, 2000. [15] Н е й м а р к Ю. И., Ланда П. С. Стохастические и хаотические колебания. — М.: Наука, 1987. [16] Пайтген X. О., Рихтер П. X. Красота фракталов. Об- разы комплексных динамических систем. — М.: Мир, 1993. [17] Понтрягин Л.С. Обыкновенные дифференциальные уравнения. М.: Наука, 1974. [18] Смирнова А. Б. Итеративные методы решения нели- нейных операторных уравнений 1 рода и их приложения: Дис....канд. физ.-мат. наук. — Екатеринбург, 1995. [19] Фейгенбаум М. Универсальность в поведении нели- нейных систем // Успехи физических наук. — 1983. — Т. 141, № 2. - С.343-374. [20] Хайрер Э., Нерсетт С., Ваннер Г. Реше- ние обыкновенных дифференциальных уравнений: Нежест- кие задачи. — М.: Мир, 1990.
[21] Чуличков А. И. Математические методы нелинейной динамики. М.: Наука. Гл. ред. физ.-мат. лит., 2000. [22] Шредер М. Фракталы, хаос, степенные законы. Мини- атюры из бесконечного рая / НИЦ «Регул, и хаотич. дина- мика». — Ижевск, 2001. [23] Шустер Г. Детерминированный хаос. Введение. — М.: Мир, 1988. [24] Briggs К. A precise calculation of the Feigenbaum constants // Math. Comput. — 1991. — Vol. 57, nr. 195. — P.435-439. [25] Briggs K. Feigenbaum Scaling in Discrete Dynamical Systems. Dissertation of the Degree Doctor. — Melburne: University of Melburne, 1997. [26] Complex Dynamic System: Mathematics Behiend Mandelbrot and Julia Sets. Proc. Sympos. Appl. Math. (Held in Cincinatti, OHIO Jan 10-11. 1994) 1994. Vol. 49. [27] Hutchinson J. E. Fractals and selfsimilarity //Indian Univ. Math. J. - 1981. - Vol. 30, nr. 5. - P. 713-747. [28] Lanford О. E. A computer — assisted prof of the Feigenbaum conjectures // Bull. Amer. Math. Soc. — 1982. — Vol. 6, nr. 3. P. 427-434. [29] Mandelbrot B.B. Fractal aspects of the iteration of z —> Az(l — z) for complex A and z. In. Nonlinear Dynamics, Helleman R.H(ed). Annals New-York Acad. Sciences. — 1980. — Vol. 357. - P. 249-259. [30] P e i t g e n H. O., Jurgen H., Saupe D. Chaos and Fractals. New Frontiers of Science. N.Y. etc.: Springer-Verlag, 1992.
[31] S i n g e r D. Stable orbits and bifurcation of maps on the interval // SIAM J. Appl. Math. Comput. — 1978. — Vol. 85, — nr. 2. - P. 260-267. [32] Smirnova A. B., Vas in V. V. Iterative approctimation of solutions of non-linear unstable problems in a Hilbert space // Russ. J. Numer. Anal. Math. Model. — 1993. — Vol. 8, nr. 2. P. 127-145.
Предметный указатель К-цикл 18 TV-фуркации системы 76 Автоколебания 151 Аттрактор 33 — странный 34 Бифуркационное значение па- раметра 26 Бифуркация 140 — рождения цикла 142 ---жесткая 147 ---мягкая 144 Вторая универсальная констан- та 44 Губка Серпинского 89 Дендрит 75 Деформационная окружность 71 Динамическая система (ДС) 9, 10 Динамические фракталы 78 Диск Зигеля 73 Домен 104 Инвариант системы 18, 118 Инвариантное множество, асимптотически устойчивое 20 — устойчивое по Ляпунову 20 Интегральная кривая 117 Канторово множество 85 Ковер Серпинского 87 Кольца Эрмана 74 Кривая Гильберта 89 Кривая Кох 85 Матрица монодромии 136 Метод Рунге-Кутта 123 — - Эйлера 122 Метрика Хаусдорфа 96 Множества подобные 79 Множество Жюлиа 67 — Мандельброта 67 — - Фату 67 — самоподобное 81 Модель «хищник - жертва» 110 — Лоренца 115 — Ресслера 115 Мультипликатор 137 Оператор Хатчинсона 82 Орбита точки 18
Осциллятор линейный 111 — хаотический 115 — химический 113 — электронный 113 Отображение подобия 79 Параболический случай 72 Первая универсальная констан- та 39 Поворотные числа 76 Показатель Ляпунова 42 Поле направлений 117 Полуотклонение между множе- ствами 95 Пыль Фату 75 Пятиугольник Дюрера 88 Размерность Хаусдорфа 83 — топологическая 82 — фрактальная 83 Регион 104 Решето Серпинского 86 Самоподобие 35, 37, 81, 91 Седло 120 Система первого приближения 127 Странный аттрактор 34 Суперцикл 51 Теорема Банаха 98 Точка периодическая 18 — покоя 18, 119 Точка бифуркации 140 Узел 120 Уравнение Ван-дер-Поля 113 — Цвитановича-Фейгенбаума 54 — удвоения 53 Условие Липшица 96 Устойчивость асимптотическая 20, 118, 119 — по Ляпунову 20, 118, 119 — экспоненциальная 118, 120 — экспоненциально-орбиталь- ная 134 Фазовая плоскость 117 — траектория 117 Фазовое пространство 117 Фазовый портрет 117 Фокус 120 Фрактал 84 — Давида 88 — Мандельброта-Гивена 86 — динамический 78 Фрактальная геометрия 11, 64 Хаос 36, 41 Характеристический показа- тель 137 Центр 120 Цикл 18, 119
Васин Владимир Васильевич Ряшко Лев Борисович Элементы нелинейной динамики: ОТ ПОРЯДКА К ХАОСУ Учебное пособие Дизайнер М. В. Ботя Технический редактор А. В. Широбоков Компьютерная верстка Д. П. Вакуленко Корректор Г. Г. Тетерина Подписано в печать 29.03.2006. Формат 60 х 84У16. Печать офсетная. Усл. печ.л. 9,53. Уч. изд. л. 8,61. Гарнитура Таймс. Бумага офсетная №1. Заказ №113. Научно-издательский центр «Регулярная и хаотическая динамика» 426034, г. Ижевск, ул. Университетская, 1. http://rcd.ru E-mail: mail@rcd.ru Тел./факс: (+73412) 500-295
ISBM 5-93972-469-8