Дискретная математика:
логика, группы, графы, фракталы

Акимов О.Е.

4.6. Динамические фракталы

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

Возьмем комплексную переменную z = x + iy и некоторую комплексную постоянную c = a + ib. Составим из них квадратичную зависимость рекурсивного типа:

z_{n + 1} = z_n² + c,  x_{n + 1} = x_n² – y_n² + a,  y_{n + 1} = 2x_n y_n + b.

Конструктивный подход всегда предполагает конкретность; пускай 

x₀ = 0, y₀ = 0, a = –0,1194, b = 0,6289.

Теперь с помощью компьютера нужно организовать цикл, обеспечивающий рекурсию, а результаты расчета — вывести на экран дисплея. В системе координат (x, y) после 1000 рекурсивных циклов мы увидим картинку, которая изображена на рис. 4.23а, а числовые значения первых шагов представлены в табл. 4.5

Таблица 4.5

Рис.	4.23	1	2	3	4	5	6	7
а	xn	–0,1194	–0,5007	–0,0979	–0,1322	–0,4615	–0,1277	–0,1410
	yn	0,6289	0,4787	0,1496	0,5996	0,4704	0,1948	0,5792
б	xn	–0,1244	–0,6805	0,0161	–0,1244	–0,6796	0,0149	–0,1244
	yn	0,7560	0,5679	–0,0169	0,7555	0,5680	–0,0161	0,7555
в	xn	–0,2370	–0,7433	0,1599	–0,2382	–0,8240	0,3066	–0,1637
	yn	0,7500	0,3945	0,1635	0,8023	0,3678	0,1438	0,8382
г	xn	0,2500	0,0724	–0,2850	–0,0245	0,2281	0,0691	–0,2496
	yn	0,4900	0,7350	0,5964	0,1501	0,4826	0,7102	0,5881
д	xn	0,2500	0,1100	–0,1935	–0,0708	0,2073	0,1173	–0,1253
	yn	0,4500	0,6750	0,5985	0,2184	0,4191	0,6238	0,5964
е	xn	–0,7000	–0,2164	–0,6542	–0,2808	–0,6230	–0,3227	–0,5983
	yn	0,0800	–0,0320	0,0938	–0,0428	0,1040	–0,0496	0,1120
ж	xn	0,3700	0,4454	0,3822	0,1161	–0,1514	0,2183	0,4029
	yn	–0,2480	–0,4315	–0,6324	–0,7314	–0,4179	–0,1215	–0,3010
з	xn	–0,4161	–0,6536	–0,1455	–0,9577	0,8342	0,3919	–0,6929
	yn	0,9093	–0,7568	0,9894	–0,2879	0,5514	0,9200	0,7210

а     б

в     г

д     е

ж     з

Рис. 4.23

Это — динамический фрактал, организованный на рекурсивной функции с комплексной переменной, весьма простого вида. Ломаная кривая называется орбитой фрактала. Стартовал он из точки с координатами (a, b) и пришел к точке притяжения с координатами x* = – 0,2381, y* = 0,4242. Фигуру, изображенную на рис. 4.23б, фракталом уже не назовешь, хотя стартовые условия для нее ненамного отличаются от первого случая. Данная орбита строго периодическая, и период ее равен трем. На следующем рис. 4.23в ломаная орбита уходит в бесконечность. На 341 цикле ордината принимает значение, равное 1,087, а y₃₄₅ = –19. Эти три конкретные орбиты представляют собой три класса орбит: фрактальные, стягивающиеся к центру притяжения (они показаны на рис. 4.23г – ж), строго периодические и хаотические, уходящие в бесконечность. Наибольший интерес вызывают, конечно, фрактальные орбиты.

Функция гарантированно уйдет в бесконечность, если на каком-то шаге |z| превысит единицу. Естественно, существует еще пограничная область, которая возникает при c = 0, |z| = 1, например, когда x₀ = cos (1,0), y₀ = sin (1,0). Точки сначала ложатся на окружность единичного радиуса (табл. 4.5), хотя и хаотическим образом, но затем, после 50-го цикла, они все больше отклоняются от окружности, что можно проследить как по числовым значениям

z₅₂ = (– 0,3008, 1,1229),  z₅₃ = (– 1,1704, – 0,6755),  z₅₄ = (0,9135, 1,5812), …,

так и по рис. 4.23е, где хорошо виден уход орбиты от линии окружности. Считается, что на границе единичного радиуса рекурсивные функции могут проявить себя по-разному: либо уйти на бесконечность, либо притянуться какой-нибудь внутренней точкой, предсказать что-либо определенное здесь невозможно.

Но вернемся назад, к первым трем рисункам (рис. 4.23 а, б и в). Давайте повнимательней присмотримся к ним: что скрывается за этими тремя видами орбит, откуда взялся центр притяжения и почему окружность единичного радиуса поделила все орбиты на два класса и как понимать периодичность. С этой целью на время перейдем от комплексных чисел z к вещественным x и исследуем ту же саму квадратичную зависимость на простейших числовых примерах

f (x):  x_{n + 1} = x_n² + a, где a — const.

Пусть сначала a = 0, тогда возможны два варианта, показанные на рис. 4.24а. Если начальная точка взята в интервале |x₀| < 1, то через несколько шагов рекурсивная функция f (x) достигает нуля, т.е. ноль является для нее неподвижной и притягивающей точкой, а единица — неподвижной и отталкивающей. Если |x₀| > 1, то орбита функции f (x) уходит на бесконечность. Для |x₀| = 1 (штрих-пунктирная линия) имеем неустойчивую ситуацию, когда малейшая флуктуация при вычислениях, приводит либо к уходу функции f (x) на бесконечность, либо она устремляется к нулю. Эти три варианта соответствуют случаям, которые и показаны на рис. 4.23 (а), (в) и (е). Теперь пускай a = 1, а процесс рекурсии стартует из точки 0 < x₀ < 1. После некоторого числа итераций f (x) выйдет на два состояния и начнет бесконечно долго колебаться, принимая значения равное то 1, то 0. Для f (x) эти точки уже не будут неподвижными (они будут неподвижны лишь для второй моды f ⁽²⁾(x); об этом ниже). Этот вариант соответствует ситуации, которая показана на рис. 4.24б.

а     б

Рис. 4.24

Теперь снова перейдем к комплексным числам. На рис. 4.25а показан тот же динамический фрактал, что и на рис. 4.23а, но только на фоне параболических кривых (обозначены пунктиром) w_{n + 1} = w_n² + c, где w = u + iv, при меняющихся значениях переменных 0,1 < u < 0,5 и –0,7 < v < 0,2. Орбита треугольной спирали начинается из точки c и целиком умещается на нижних ветвях семейства парабол. Орбита четырехугольной спирали, изображенная на рис. 4.25б, опирается уже на обе ветви семейства парабол, занимающих область 0,1 < u < 0,29 и –0,7 < v < 0,7. На рис. 4.25в закон рекурсии с простого квадратичного изменен на более сложный, дробно-квадратичный; при c = –0,11 + 0,52i имеем:

А на рис. 4.25г орбита строится в соответствии с дробно-линейным законом, где используются две комплексные константы:

c = –1,13 + 1,33i, d = 1,44 – 1,0i

В этом случае область функции w (пунктирные кривые) простирается в пределах 0,9 < u < –3,1 и –10 < v < 0,4; левую границу по v можно брать и –100; это зависит от того, насколько близко вы хотите подойти к точке (0,1) комплексной плоскости.

Эти симпатичные кривые очень заинтересовали математиков. Чтобы исследовать влияние константы c = a + ib на поведение функции z_{n + 1} = z²_n + c, поступают следующим образом. Для фиксированного c определяют множество точек z, модуль которых превысил единицу. Поскольку сама единица, как было сказано, еще ничего не гарантирует, то ориентируются на числа больше единицы. Если, скажем, |z| > 2, то далее орбита функции наверняка устремится к бесконечности. Множество значений z, для которых |z| < 2, называется заполняющим множеством Жюлиа.

Гастон Жюлиа (1893—1978) вместе с Пьером Фату (1878 — 1929) вскоре после Первой мировой войны исследовали проблемы, связанные с рекурсивной функцией от комплексной переменной. Известный математик Джон Милнор эти две фамилии, наверное, поставил бы в обратном порядке, поскольку считал Фату первопроходчиком и человеком, который внес, как он написал, «наиболее весомый и внушительный вклад» в развитие теории. Разъясняя, почему в 1918 г. «Гран-При математических наук» Парижской Академии наук получил все же 25-летний Жюлиа, а не кто-то другой, Милнор в своей книге «Голоморфная динамика» я звительно заметил, что молодой человек «имел некоторые преимущества, связанные с его статусом раненого героя войны». Он также не преминул сообщить о промахе Жюлиа в решении задачи локальной линеаризации: «В 1919 году Жюлиа объявил о полном разрешении этого вопроса для рациональных отображений степени два и выше, показав, что такая линеаризация никогда невозможна, однако его доказательство было ошибочным».

Вместо словосочетания «рекурсивная динамическая функция» Милнор пользуется словосочетанием «итерационное голоморфное отображение». Понятия итерация и рекурсия, отображение и функция можно считать синонимичными. Видимо, не будет большой ошибкой соотнести его понятие голоморфной функции с понятием гладкой или дифференцируемой функции. Милнор пишет: «Если V включено в C — открытое множество в комплексной плоскости, то функция f: V → C называется голоморфной, или комплексно аналитической, если ее первая производная z → f '(z), равная

определена и непрерывна, как функция из V в C, или, что эквивалентно, если в окрестности любой точки z₀ из V функция f разлагается в ряд по степеням z – z₀, который сходится к f (z) в некоторой окрестности z₀». Это определение понадобилось нам для того, что бы иметь представление, с какого рода функциями (неважно, как их называть и записывать) мы будем иметь дело в дальнейшем.

Сейчас мы вновь воспользуемся книгой Кроновера, чтобы привести из нее простой алгоритм для построения множества Жюлиа. Правда, в связи с этим Милнор предупреждает о некоторых трудностях, которые могут подкарауливать слишком самонадеянного программиста.

Вход:

c = a + ib — фиксированная точка
u, v — центр окна
s — сторона окна
p — число пикселов, умещающихся на стороне s

Выход: Графическое изображение множества Жюлиа

Инициализация:

(u – s/2, u + s/2) × (v – s/2, v + s/2) — графический экран

Шаги:

for m = 1 to p
x₀ = u – s/2 + ms/p
 for n = 1 to p
 y₀ = v – s/2 + ns/p
x = x₀, y = y₀, z = 0, iter = 1
while iter < 20 iter = iter + 1
x₁ = x² – y² + a, y₁ = 2xy + b, x = x₁, y = y₁, z = x² + y²
if z > 4, выход из цикла, end if
end while if z < 4, построить точку (x₀, y⁰), end if
end for
end for

Кроновер представил и результаты расчета по этому алгоритму — рис. 4.26.

а   б   в

г   д   е

Рис. 4.26

Множества Жюлиа симметричны относительно оси абсцисс и имеют фрактальную структуру, поскольку содержат элементы самоподобия. Рис. 4.26а соответствует рис. 4.23а для динамического фрактала с параметрами a = – 0,1194 и b = 0,6289; рис. 4.24б тоже соответствует рис. 4.23б; а вот рис. 4.26в уже не отвечает рис. 4.23в, так как при таких значениях a и b орбита уходит на бесконечность и никакого множества Жюлиа получиться не может; фрактал, изображенный на рис. 4.26в, рассчитан для a = – 0,20 и b = 0,75. Часто говорят не о множестве Жюлиа, а о границе Жюлиа; именно граница множеств Жюлиа представляет интерес, так как она определяет связанность множества Жюлиа. Следующие три фрактала показаны только своими границами: для рис. 4.26г a = 0,25 и b = 0,52; для рис. 4.26д a = 0,377 и b = – 0,248, что близко к параметрам фрактала, показанного на рис. 4.23ж, и для рис. 4.26е имеем a = – 0,7382 и b = 0,0827, что близко к параметрам фрактала, показанного на рис. 4.23е. Каждая из трех границ представляет собой один замкнутый контур, т.е. для последних трех фракталов множества Жюлиа будут односвязанными.

Немало превосходных рисунков приведено и в упомянутой книге Милнора. Автор детально обсуждает не только особенности множеств Жюлиа для квадратичной зависимости, но поясняет структуру более сложных голоморфных отображений. Этому посвящен больший объем его «Голоморфной динамики». Мы же ограничимся тем, что воспроизведем несколько характерных структур (рис. 4.27 – 4.29) и укажем функциональную зависимость, по которой они были построены (табл. 4.6).

а     б     в

г     д

Рис. 4.27

а     б     в

г     д

Рис. 4.28

а     б     в

г     д

е     ж

Рис. 4.29

Таблица 4.6

Рис.	f (z)
4.27a	z2 + (0,99 +0,14i)z
4.27б	z + exp(z) – 1
4.27в	z2 + i
4.27г	z2 – 1
4.27д	z2 – 3/2
4.28a	z2 – 0,00001/z3
4.28б	– 0,138 (z + 1/z) – 0,303
4.28в	(2z/(1 + z2))2
4.28г	z3 – 0,48z + 0,7063 + 0,5029i
4.28д	z2 – 0,7443 + 0,1212i
4.29a	z2 – 0, 4245 + 0,2075i
4.29б	sin(z)
4.29в	z5 + (0,8 + 0,8i)z4 + z
4.29г	z2 + zexp(6πi/7)
4.29д	z2 + zexp(1,26πi)
4.29е	z2exp(1,23πi)(z – 4)/(1 – 4z)
4.29ж	z2 – 0,765 – 0,12i

Если вернуться к первым трем рисункам, где множества Жюлиа затушеваны в черный цвет, то можно заметить следующее: первый фрактал (рис. 4.26а) явно односвязанный, второй (4.26б) находится на границе связанности, когда отдельные «листья кактуса» едва соединяются друг с другом, а третий фрактал (рис. 4.26в) распался на отдельные области, т.е. множество Жюлиа распределилось по огромному числу несвязанных друг с другом областей. Этот несвязанный фрактал особенно интересен, так как он построен для пограничных параметров a и b, когда небольшое их изменение приводит к уходу орбиты на бесконечность (рис. 4.23в). Таким образом, вся совокупность несвязанных фракталов-кактусов образуют новую границу, за пределами которой никакие множества Жюлиа существовать не могут. Эта граница называется границей Мандельброта. Внутри этой границы находится множество Мандельброта, каждым элементом которого является связанным множеством Жюлиа.

Для графического построения границы Мандельброта сгодится алгоритм для построения границы Жюлиа. Единственное отличие состоит в том, что теперь параметр c = a + ib становится варьируемой величиной, т.е. область нахождения множества Мандельброта лежит в координатной сетке (a, b), а не (x, y). Это приводит к замене равенств

x₀ = u – s/2 + ms/p,  y₀ = v – s/2 + ns/p,  x = x₀, y = y₀,

на равенства

a = u – s/2 + ms/p,  b = v – s/2 + ns/p,  x = 0, y = 0.

Результат работы такого модернизированного алгоритма представлен на рис. 4.30. На рисунке (а) дана укрупненная картина, на которой видно, что множество Мандельброта, как и множества Жюлиа, симметрично относительно оси абсцисс. На рисунке (а) имеется прямоугольная вырезка, которая на рисунке (b) изображена уже в увеличенном виде. Но и на этой увеличенной картинке имеется своя вырезка, развернутая на рисунке (c). Наконец, на рисунке (d) показана еще одна увеличенная вырезка предыдущего рисунка. На всех картинках видна фрактальная структура; самоподобие затрагивает не только общую конфигурацию, но и ее отдельные части. Сколько бы мы ни детализировали границу Мандельброта, мы всякий раз находили бы те или иные элементы подобия. По масштабной шкале можно углубляться до бесконечности, здесь все зависит от разрядности процессора вашего компьютера.

Рис. 4.30

На рис. 4.31а показана граница Мандельброта, размещенная в координатной системе (a, b), на которой видно, что самая большая часть множества Мандельброта (1) расположена в пределах – ¾ < a < ¼ и имеет форму известной кривой кардиоиды. Вторая по величине область (2) имеет форму круга, который занимает на оси абсцисс область –0,75 < a < –1,25. Слева от нее расположен крохотный кружок с числом 4, занимающий область –1,25 < a < –1,368. Еще левее видна едва заметная область с числом 8. Таким образом, появляется последовательность, которая имеет предел:

a₀ = 0,25, a₁ = – 0,75, a₂ = – 1,25, a₃ = – 1,368, 

a₄ = – 1,394, … a_∞ = – 1,401155...

Если из членов ряда составить отношение:

и найти его собственный предел, то получим константу, называющуюся постоянной Фейгенбаума F = 4, 6692016091029909…

Откуда взялся ряд чисел: 2, 4, 8, …? Его объяснение вытекает из рис. 4.31б, на котором изображены три моды квадратичной функции при c = –2:

f (z) = z² + c, f ⁽²⁾(z) = f (f (z)) = (z² + c)² + c,

f ⁽⁴⁾(z) = f ⁽²⁾(f ⁽²⁾(z)) = f (f (f (f (z)))).

а     б

Рис. 4.31

Наклонная прямая пересекает f ⁽ⁿ⁾(z) 2ⁿ раз. Таким образом, точки a_n отвечают высшим модам функции f ⁽ⁿ⁾(z), когда орбита периода 2ⁿ изменяется на орбиту периода 2^2^n (ниже об этом говорится более подробно). Первый нечетный период 3 возникает при a = –1,7549, где находится еще одна клякса в миниатюре; затем появляются следующие нечетные периоды 5, 7, 9 и т.д. Кроновер сообщает, что А.Н. Шарковский еще в 1964 г. доказал теорему: если функция f (z) имеет точки периода 3, 5, 7, ..., то будут существовать точки с периодом 3·2, 3·2², 3·2³, …, 5·2, 5·2², 5·2³, ..., 7·2, 7·2², 7·2³, … В силу фрактальной природы динамического процесса области с четными и нечетными периодами, в частности 3, 4 и 5, появляются вблизи определенных точек кардиоиды, т.е. вокруг области с единственной точкой притяжения.

Теперь разберемся, откуда здесь появилась кардиоида? Напомним, что кардиоида была впервые описана еще в конце XVII в. голландским математиком Коерсма, а в 1741 г. итальянский ученый Кастильон из-за сходства с сердцем окрестил ее принятым сейчас названием (буквально «сердцевидная»). Кардиоида, являясь частным случаем улитки Паскаля, относится к кривым четвертого порядка и имеет уравнение (x² + y² – px)² = q²(x² + y²) при p = q.

Рис. 4.32

Кардиоида прочерчивается точкой окружности радиуса ¼, катящейся по окружности с таким же радиусом, центр которой находится в начале координат. Слева от кардиоиды (рис. 4.32) как раз находится такая окружность. Площадь кардиоиды в точности равна площади шести кругов, а длина кривой — 16 радиусам, т.е. четырем.

Чтобы понять, как получается кардиоида, нужно помнить, как строится функция z_{n + 1} = z_n² + c. Она получается в результате рекурсии или итерационной процедуры, которая опирается на диагональ квадрата z_{n + 1} = z_n (рис. 4.24). Аналогичная итерационная процедура строится и для высших мод рекурсивной функции; прямая z = z_{n + 1} = z_n, пересекающая моды под углом 45° (рис. 4.31б), является здесь основополагающим элементом динамической модели. Отсюда для первой моды получаем условие z² + c = z. Чтобы определить, является ли та или иная точка притягивающей или отталкивающей, необходимо найти производную от функции и определиться: если |f'(z)| < 1, то z — притягивающая точка; если |f '(z)| > 1, то z — отталкивающая. Найдем производную функции и приравняем ее единице: 

|f '(z)| = |2z| = 1; так как z = exp (iφ), то z = [exp (iφ)]/2.

Подстановка этого значения z в условие z² + c = z дает уравнение кардиоиды, сдвинутой по оси x на величину ¼, или

     

Точки, лежащие внутри кардиоиды, являются центрами притяжения для всех орбит, так как для них выполняется условие |f '(z)| < 1. Внешние же точки являются отталкивающими центрами, так как для них |f '(z)| > 1, и все орбиты с этими точками уходят на бесконечность.

Для второй моды f ⁽²⁾(z) производная равна 4z(z² + c), а решение определяется равенством (z² + c)² + c = z, которое удобно представить в виде двух скобок:

(z² + c)² + c – z = (z² + c – z)(z² + z + c + 1) = 0.

Первая скобка, отвечающая уравнению z² + c – z = 0, дает две неподвижных точки z₁z₂, которые являются и решениями уравнения z² + z + c + 1 = 0. По теореме Виета произведение корней равно свободному члену z₁z₂ = c + 1, кроме того,

z₁² + c – z₁,  z₂² + c – z₂,  z₁² + c – z₂, z₂² + c – z₁.

Граница притяжения определяется равенством |4z(z² + c)| = |4z₁z₂| = |4(c + 1)| = 1, следовательно, круг с центрами притяжения периода 2 по оси абсцисс расположен в пределах –0,75 < a < –1,25. Как было сказано, этот круг является образующим для кардиоиды. Круг периода 4 на оси абсцисс занимает отрезок –1,25 < a < –1,368. Круг же периода 8 на оси абсцисс будет еще меньше, так как его диаметр определяется областью –1,368 < a < –1,394. Области периода 3 находятся из условия:

c³ + 2c² + (1 + z)c + (1 + z)² = 0.

Оно дает два круга радиусом 1/8, расположенных в точках (–0,125, ±0,775), и еще одну крохотную кардиоиду, расположенную по оси абсцисс в точке –1,755. Указанные моды (за исключением последней), изображенные на рис. 4.32, определяют общий вид «береговой линии» множества Мандельброта.

Рис. 4.33

Милнор привел границы множества Мандельброта для функции f (z) = z² + λ z (рис. 4.33); центр левого большого круга расположен в начале координат, центр правого большого круга расположен по оси абсцисс на расстоянии +2.