Что такое бассейны ньютона

Зона кода

Продолжим заниматься построением простейших анимаций с использованием языка C99, библиотеки pgraph и видеоредактора VirtualDub. Напомню, что ранее с помощью перечисленных инструментов были построены анимированные фигуры Лиссажу. А в данной статье будет рассмотрено создание анимированных бассейнов Ньютона.

Бассейны Ньютона представляют собой множества точек комплексной плоскости, построенные по определённому закону (о нём мы поговорим далее). Другое название бассейнов Ньютона — фракталы Ньютона. Фрактал — это множество, обладающее свойством самоподобия, заключающимся в возможности его разбиения на части, имеющие такую же структуру, что и само множество. В свою очередь, каждая из частей может быть разбита на собственные части, обладающие таким же свойством, и т. д.

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

В нашей статье мы будем заниматься визуализацией одного частного случая бассейнов Ньютона. Надеюсь, что нам удастся сгенерировать какие-нибудь более или менее забавные картинки и видеоролики.

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

Алгоритм построения бассейнов Ньютона

В основе алгоритма лежит один из методов приближённых итераций нахождения корня действительного уравнения вида f(x) = 0, где f(x) — это функция, имеющая не обращающуюся в 0 производную на некотором промежутке X, которому принадлежит, также, и искомый корень уравнения s. Метод заключается в построении последовательности <x_k> по следующему правилу:

x 0 ∈ X , x k = x k — 1 — f x k — 1 f ‘ x k — 1 , k = 1 , 2 , . . . .

Начальное значение x₀ стараются выбрать как можно более близким к s. Известно, что при выполнении некоторых условий, на которых мы сейчас не будем останавливаться, все члены последовательности также принадлежат X, причём

Другими словами, рассматриваемая последовательность при выполнении упомянутых условий сходится к искомому корню уравнения.

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

Описанный метод называется методом Ньютона или методом касательных. Последнее название обусловлено следующим геометрическим смыслом метода: x_k при k ≥ 1 представляет собой абсциссу точки пересечения касательной к графику функции y = f(x) в точке (x_k-1, f(x_k-1)) с осью абсцисс.

Этот метод естественным образом можно переформулировать и на случай решения комплексного уравнения вида f(z) = 0, где f(z) — это комплексная функция, дифференцируемая на некотором комплексном множестве W , производная которой не обращается на этом множестве в 0. Соответствующая последовательность также задаётся рекуррентно:

z 0 ∈ W , z k = z k — 1 — f z k — 1 f ‘ z k — 1 , k = 1 , 2 , . . . .

Будем говорить, в дальнейшем, что точка z₀ «порождает» последовательность <z_k>.

Если комплексный корень s ∈ W уравнения f(z) = 0 существует, то, как и в предыдущем случае, при выполнении определённых условий он может быть найден через предел:

Поскольку уравнения вида f(z) = 0 могут иметь несколько корней, то и данный предел может иметь разные значения в зависимости от выбора начального значения z₀. Множество всех точек плоскости, таких, что последовательности, порождаемые ими, сходятся к корню s_n уравнения f(z) = 0, называется областью притяжения данного корня. Эти области притяжения и называются бассейнами Ньютона.

Визуализация бассейнов Ньютона заключается в том, что области притяжения разных корней мы «раскрашиваем» разными цветами. В результате получаются довольно забавные картинки. Именно таким раскрашиванием мы и будем заниматься.

Давайте теперь перейдём к конкретному уравнению, области притяжения корней которого мы будем строить. Вот оно:

Несложно найти корни этого уравнения. Для этого достаточно извлечь комплексный корень пятой степени из единицы. В результате извлечения получаем 5 значений корня: s_n = exp(2iπn/5), где n = 0, 1, 2, 3, 4.

Наше уравнение можно представить в виде f(z) = 0, где f(z) = z 5 — 1. Имеем:

z — f ( z ) f ‘ ( z ) = z — z 5 — 1 z 5 — 1 ‘ = z — z 5 — 1 5 z 4 = z — z 5 + 1 5 z 4 = 0 , 8 z + 0 , 2 / z 4 ,

откуда получаем следующую рекуррентную формулу общего члена последовательности, сходящейся к корню исходного уравнения:

z k = 0 , 8 z k — 1 + 0 , 2 / z k — 1 4 , k = 1 , 2 , . . . .

Известно, что данная последовательность существует и сходится к одному из корней уравнения при любых комплексных начальных значениях z₀, за исключением точки 0.

Таким образом, можно уже сформулировать алгоритм построения бассейнов Ньютона.

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

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

Критерий окончания вычисления членов последовательности будет таким же, как и в действительном случае. Оно прервётся, как только перестанет выполняться условие

Здесь ε — это некоторое заранее заданное положительное число, подбираемое эмпирически.

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

Проведём из центра этой окружности, совпадающего с началом координат, 5 лучей, проходящих через середины сторон пятиугольника. Эти лучи будут иметь уравнения вида arg z = φ, где φ принимает значения ± π / 10, ± 3 π / 10, π. Лучи разбивают комплексную плоскость на 5 секторов. В каждом секторе находится по одному корню, равноудалённому от лучей, образующий данный сектор. Нам нужно лишь выяснить, в каком из 5 секторов лежит точка z_m.

Пронумеруем секторы числами от 0 до 4 таким образом, чтобы номер сектора совпадал с лежащим в нём корнем. Для того, чтобы определить номер сектора, содержащего z_m, нужно выполнить целочисленное деление аргумента числа z_m на π / 5. Если аргумент измеряется от —π до π, то результат, равный 0, соответствует 0-му сектору, результат, равный 1 или 2 — 1-му сектору, равный 3 или 4 — 2-му, равный -1 или -2 — 3-му, а равный -3 или -4 — 4-му.

Таким образом, номер корня, к области притяжения которого принадлежит число z₀, находится как номер сектора, в котором лежит член последовательности, порождённой числом z₀, вычисленный последним.

Теперь уже можно переходить к тексту программы, рисующей бассейны Ньютона.

Статичное изображение бассейнов Ньютона

Ниже размещён текст программы, выполняющей построение бассейнов Ньютона

Программа состоит из двух функций. Все действия по построению фрактала Ньютона выполняются функцией fractal1() , которая вызывается из функции main() . Такой подход использован для того, чтобы для других построений фракталов к программе можно было бы добавлять функции fractal2() , fractal3() и т. д., а внутри main() лишь заменять вызов одной функции вызовом другой.

В строке 1 подключаем нашу графическую библиотеку pgraph. Очень хорошо, что язык C99 поддерживает работу с комплексными числами на уровне синтаксиса. Благодаря этому, нам не потребуется реализовывать эту поддержку самостоятельно. Для использования данной возможности, подключаем заголовочный файл (см. стр. 2), содержащий объявления функций cpow() , carg() и cabs() и определения макросов complex и I .

В строках 4-10 определены следующие макросы-конcтанты:

• W и H — ширина и высота изображения соответственно. В качестве размера изображения я выбрал 854×480, поскольку данный размер является одним из стандартных размеров кадра, поддерживаемых сервисом Ютьюб, на который будет загружена анимация (но пока мы рисуем лишь статичную картинку, построение анимации будет рассмотрено позже).
• X0 и Y0 — координаты точки O (начала координат) в системе координат, поддерживаемой библиотекой pgraph. Точка O будет располагаться почти в центре нашей прямоугольной области.
• PXL — длина одного пиксела. Другими словами, величина, обратная PXL , равна количеству пикселей, содержащихся в отрезке единичной длины. В нашем случае значение PXL равно 0,1, т. е. отрезок длины 1 состоит из 10 пикселей. Чем больше PXL , тем меньше площадь отображаемой прямоугольной области комплексной плоскости.
• EPS — число ε , участвующее в условии продолжения вычисления членов последовательности. Чем меньше EPS , тем большее количество членов последовательности будет найдено и тем больше времени на это будет потрачено. А слишком большое значение EPS может привести к неверному определению того, к области притяжения какого корня принадлежит исследуемая точка. Я взял значение 0,1, но, скорее всего, его можно немного увеличить.
• PI5 — число π / 5 с точностью 15 знаков после запятой.

Переходим к описанию функции fractal1() . В строке 15 создаём изображение. Далее посредством двух циклов for , вложенных один в другого (стр. 16, 17) перебираем все пиксели изображения.

Для каждого пикселя вычисляем координаты x и y точки, которая ему соответствует (стр. 19, 20), а саму точку комплексной плоскости помещаем в переменную z (стр. 21). Далее, в случае, если точка не совпадает с началом координат (см. условие в стр. 22) выясняем, области притяжения какого корня она принадлежит.

Для этого в цикле do-while (стр. 25-30) вычисляем члены последовательности, порождённой нашей комплексной точкой, до тех пор, пока модуль разности двух членов последовательности, вычисленных последними, не станет меньше, чем EPS .

На всякий случай поясним, что функция cpow( ) возводит комплексное число, переданное ей в качестве первого аргумента, в комплексную степень, равную второму аргументу, а функция cabs() возвращает модуль комплексного числа.

Объявляем переменную col для хранения цвета, которым будет закрашен текущий пиксель (стр. 31). Далее с помощью функции carg() находим аргумент последнего вычисленного члена последовательности и делим полученный результат на PI5 , после чего, в зависимости от сектора, в котором оказывается лежащим данный член, присваиваем col один из пяти цветов (стр. 32-53).

Пиксели, соответствующие точкам, попадающим в область притяжения корней с номерами 0, 1, 2, 3 или 4, мы закрашиваем красным, зелёным, синем, жёлтым или фиолетовым цветом соответственно. Само раскрашивание осуществляется в строке 54.

Заметим, что точка с координатами (0, 0) (начало координат) не подлежит окраске и имеет цвет по умолчанию (белый).

Нам остаётся лишь сохранить изображение в файле (стр. 57), после чего удалить его (стр. 58).

Из функции main() (стр. 61) вызывается функция fractal1() (стр. 63) и возвращается нулевое значение (стр. 64).

Выполнение программы приводит к появлению в корневом каталоге исполняемого файла графического файла result.bmp, содержащего изображение, приведённое ниже (его можно увеличить, щёлкнув по нему):

Мы видим 5 лучей, выходящих из начала координат (это, кстати, те самые лучи, проходящие через середины сторон воображаемого пятиугольника) с «насаженными» на них своеобразными «сердечками». Если мы будем «скользить» взглядом по лучу по направлению к центру, то заметим, что эти «сердечки» имеют тенденцию уменьшаться по мере приближения взгляда к началу координат. Однако, если присмотреться, то можно заметить, что уменьшение прекращается, как только мы достигаем точки, удалённой от начала координат на единицу. За этой точкой на луче имеется одно-единственное «сердечко», превышающее предыдущее по размерам.

Но давайте, всё же, «приблизим» картинку, для того, чтобы детально рассмотреть эти 5 «сердечек», являющихся последними на каждом луче. Для изменения масштаба необходимо переопределить макрос PXL . Если мы хотим увеличить масштаб, например, в 20 раз, то значение PXL должны уменьшить в 20 раз, т. е. принять равным 0,005. Переопределим строку 8 следующим образом:

Теперь, после выполнения программы получаем следующее изображение (как и прежде, его можно увеличить с помощью щелчка):

Увеличенные бассейны Ньютона

Падение в центр бассейнов Ньютона

А теперь давайте выполним плавное приближение (увеличение масштаба) и сохраним его в виде анимации.

Теперь уже будем строить не одно изображение, а несколько, являющихся кадрами анимации. Количество кадров будет задаваться макросом FRAMES . Длина пикселя будет уже не фиксированной, а изменяющейся. После построения очередного кадра будет вычисляться новое значение длины посредством умножения старого на постоянный множитель, задаваемый макросом FACTOR . Значение множителя подбирается эмпирически; оно должно обеспечить плавность анимации.

Определим эти два макроса следующим образом:

Итак, анимация будет состоять из 300 кадров, а длина пикселя каждый раз будет уменьшаться на 2%, что соответствует увеличению масштаба. Построение анимации реализуем в функции fractel2() , немного переделав функцию fractal1() . Вот код новой функции:

Как мы видим, появилась переменная pxl (стр. 6), предназначенная для хранения изменяющейся длины пикселя. Начальное значение pxl определяется макросом PXL (считается, что оно равно 0,1). Также появился новый цикл по номерам кадров (стр. 7). В конце каждой его итерации изображение записывается в графический файл, имеющий имя вида xxx.bmp (стр. 50, 51), где xxx — трёхзначный номер кадра (если номер состоит из одной или двух цифр, то он дополняется до трёхзначного нулями слева). Все файлы сохраняются в директории result, которая к моменту выполнения программы должна находиться в корневом каталоге исполняемого файла.

Ну а самое последнее действие, выполняемое на текущей итерации, — это изменение длины пикселя (стр. 52).

Вот и все отличия fractal2() от fractal1() .

Нам осталось лишь в теле функции main() заменить вызов одной функции вызовом другой:

После выполнения программы в директории result появляются 300 графических файлов. Наша задача — создать из них видеофайл. О том, как это сделать с помощью видеоредактора VirtualDub, мы уже рассказывали ранее. Но в этот раз, в отличие от предыдущего, мы изменим частоту видеокадров с 10 кадров/сек. на 25 кадров/сек. Для этого перед сохранением видеофайла выберем в меню «Видео» пункт «Частота кадров. «, после чего в открывшемся диалоговом окне активизируем текстовое поле, расположенное после надписи «Изменить на», введём в него значение 25 и нажмём на кнопку «OK» (см. рисунок).

Изменение частоты кадров

К сожалению, использовать для сжатия без потерь кодек Huffyuv не удалось: итоговый видеофайл оказался испорченным. Причина этого кроется, видимо, в том, что размер кадра по вертикали не кратен 4. Поэтому ролик был сохранён без сжатия. Размер его получился равным 352МБ, а его продолжительность составила 12с. Видеофайл был загружен на Ютьюб. Вот какой получился «кислотный» мультфильм:

К сожалению, качество ролика после его загрузки на Ютьюб оказалось отвратительным из-за сильного сжатия, использующегося данным сервисом. Представление об анимации получить можно, но читателю стоит самостоятельно создать ролик и просмотреть его у себя на компьютере.

Падение в центр бассейнов Ньютона с вращением

Напоследок — ещё один вариант анимации для разнообразия. На этот раз изменение масштаба будет сопровождаться вращением бассейнов Ньютона вокруг начала координат по часовой стрелке. Реализовано это будет очень просто: каждый кадр будет повёрнут относительно другого на небольшой постоянный угол, определяемый макросом STEP :

Значение этого угла, опять же, подбирается эмпирически, из соображений плавности анимации. Функцию, в которой изменение масштаба с вращением будет реализована, назовём fractal3() . Она будет получена из функции fractal2() с помощью небольших изменений. Вот её код:

Рассмотрим отличия этой функции от fractal2() . В строке 7 создаётся переменная alpha , которой присваивается начальное значение 0. Она предназначена для хранения угла поворота бассейнов Ньютона, содержащихся в текущем кадре, относительно начального положения, соотвествующего значению 0. Таким образом, в ходе цикла по кадрам значение alpha будет увеличиваться с шагом STEP .

На каждой итерации цикла вычисляются значения косинуса и синуса угла поворота, которые сохраняются в перменных c и s соотвественно (стр. 10). Эти значения используются для нахождения текущей комплексной точки (см. стр. 16). Теперь эта точка, помещённая в переменную z , уже не совпадает с точкой, координаты которой содержатся в перменных x и y , а «смещена» от неё на угол поворота, хранящийся в alpha , против часовой стрелки (благодаря этому, сама «картинка» оказывается повёрнутой по часовой стрелке).

Последнее отличие от fractal2() — это строка 55, в которой изменяется угол поворота. Данное действие являеся последним в текущей итерации цикла по кадрам.

Разумеется, нельзя забывать и об изменени, вносимом в функцию main() :

Новый набор кадров точно так же, как и ранее, был «превращён» в видеоролик, который был загружен на Ютьюб:

И снова отвратительное качество! Ну не любит Ютьюб мелкие детали, тем более, находящиеся в сравнительно быстром движении. А у нас этого «добра» — предостаточно! Опять призываю читателя «не верить» Ютьюбу, а сгенерировать ролик самостоятельно и просмотреть его с помощью плеера на компьютере.

Заключение

Ну вот, пожалуй, и всё, о чём я хотел рассказать. Это была история о том, как с помощью достаточно простых программ и с привлечением несложного математического аппарата можно рисовать довольно забавные картинки.

Источник