Графический процессор и одержимость фракталами

Вы одержимы фракталами и не знаете, как можете их использовать по максимуму вместе с MATLAB? Тогда эта статья для вас. В ней я расскажу, как можно с помощью графического процессора улучшить производительность вычисления фракталов и многое другое.

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

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

Рисунок 1. Фрактальная структура в дереве

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

Это множество точек C на комплексной плоскости, для которых рекуррентное соотношение

Z_n+1 = Z_n² + C при Z₀ = 0 задает ограниченную последовательность. Иными словами, на каждой итерации число возводится в квадрат и прибавляется константа. Если за определенное число шагов последовательность осталась в заданных пределах, то значение константы принадлежит множеству. Выглядит это так:

Рисунок 2. Множество Мандельброта в комплексной плоскости

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

В данной статье мы затронем несколько фракталов - "Пылающий корабль", "Башня из возведения в степень" и глобальную сходимость метода Ньютона. 
Также задействуем некоторое оборудование - графический процессор NVIDIA^® Titan V, расположенный в специальной док-станции, подключенной к ноутбуку через интерфейс Thunderbolt 3, что позволит уменьшить время вычислений фракталов почти в 300 раз по сравнению с вычислительным процессором (ЦП).

Визуализация множества Мандельброта

Первое опубликованное изображение множества было представлено на строчном принтере в 1978 году Робертом Бруксом и Питером Мательским:

Рисунок 3. Изображение множества Мандельброта на строчном принтере

Первое же цветное изображение впервые появилось на обложке журнала Scientific American в 1985 году. Это было во время, когда графические дисплеи компьютеров становились общедоступными. 

Рисунок 4. Цветное изображение множества Мандельброта.

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

Внутри множества Мандельброта

Что случается в Мандельброте, остается в Мандельброте.

Мандельброт охватывает простые итерации с комплексными числами, начиная из начальной точки z₀. Множество Мандельброта - область в комплексной плоскости, состоящая из чисел z₀, для которых траектории определены как

                                                            z_k+1 = z²_k + z₀, k = 1, 2, ...

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

Цветное изображение, представленное выше, - общее представление множества Мандельброта.

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

Изображение ниже - эскиз геометрии множества Мандельброта. Наибольший элемент - кардиода в виде сердца, окруженная параметрической кривой вида

                                                            z = e^it/2 − e^2it/4, −? ≤ t ≤ ?

Рисунок 5. Эскиз геометрии множества Мандельброта

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

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

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

Вычисление Мандельброта

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

Ниже представлена функция, которая лежит в основе вычислений.

function kz = mandelbrot(z0,depth)
    z = z0;
    kz = 0;
    while (z*conj(z) <= 4) & (kz <= depth)
        kz = kz + 1;
        z = z*z + z0;
    end
end

Входными данными является комплексная скалярная начальная точка z₀ и действительный скалярный параметр, назовем его depth. На выходе - счетчик kz.

Если итерация прекращается, потому что z находится за пределами круга радиуса два, тогда z было суждено стать бесконечным, и счетчик kz можно использовать в качестве индекса в цветовой карте. С другой стороны, если z выходит за пределы итерации глубины, то объявляется, что он находится в наборе Мандельброта.

Давайте создадим сетку точек на комплексной плоскости, покрывающую квадрат шириной w с центром в комплексной точке zc.

grid = 512;
s = w*(-1:1/grid:1);
[u,v] = meshgrid(s+real(zc),s+imag(zc));

Есть единственная разница между программой, которая работает на ЦП, и той, которая работает на GPU. Чтобы сгенерировать сетку на ЦП:

z0 = u + i*v;

Для графического процессора:

z0 = gpuArray(u + i*v);

Тогда для любого другого процессора выражение

kz = arrayfun(@mandelbrot, z0);

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

Теперь давайте посмотрим на два необычных изображения, заимствованных из внешней границы Мандельброта.

На внешней границе

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

Рисунок 6. Долина морских коньков

Поскольку множество Мандельброта является самоподобным, оно содержит бесконечное количество миниатюрных Мандельбротов, каждый из которых имеет ту же форму, что и большой. Тот, что показан на рисунке ниже, имеет коэффициент увеличения 10¹⁰.

Рисунок 7. Миниатюрный Мандельброт

Скрытая ?

Замечательный результат был обнаружен в 1991 году Дэйвом Боллом, тогда аспирантом Государственного университета Колорадо. Болл исследовал поведение итерации Мандельброта в Долине морских коньков. Долины сужаются по мере приближения к оси, которую они встречаются точке (-3/4,0), красной точке на Рисунке 5.

Мы можем повторить вычисление Болла на маленькой сетке с центром рядом с осью в точке

-3/4 + yi для крошечного значения y и мнимой единицы i.

y = 1.0e-07
zc = -3/4 + y*i

Сделайте сетку достаточно большой, чтобы она касалась оси:

width = 2*y
grid = 4

Выберите параметр depth, обратно пропорциональный y, что сделает его огромным:

depth = 4/y

С этими параметрами запустите приведенный выше код. Тогда вы получите:

kz =

40000000 40000000 40000000 40000000 40000000
40000000 40000000 40000000 40000000 40000000
40000000 40000000 31415926 40000000 40000000
40000000 40000000 20943951 40000000 40000000
40000000 40000000 15707963 40000000 40000000

Посмотрите на количество итераций в центре сетки. Увидели знакомое значение? (число ?)

Это не случайность. В статье 2001 года Аарона Клебаноффа «? в множестве Мандельброта» анализируются аналогичные вычисления в точке перед кардиоидой.

Рассмотрим еще любопытный вариант итерации Мандельброта.

Пылающий корабль

Пылающий корабль произошел из странной итерации:

z_k+1 = F(zk) + z₀, k = 1, 2, ...

где

F(z) = (|Re(z)| + i|Im(z)|) 2

Данная итерация странная, потому что функция F(z) не является аналитической. Меня заинтересовала эта итерация из-за сверхъестественного сходства на следующих изображениях.

Начальная область, показанная на рисунке 8, представляет собой квадрат шириной 3,5 с центром в -0.5-0.5i. Я добавил стрелку, как указатель на интересующую нас область следом за кораблем.

Рисунок 8. Горящий корабль, начальная область.

Увеличьте след корабля в 500 раз до точки -1.861-.002i. Примените bone палитру, чтобы рисунок казался леденеющим, а не пылающим. На рисунке 9 показан полученный фрактал рядом с фотографией 1915 года корабля «Эндьяранс», исследователя Антарктиды Эрнеста Шеклтона, захваченного льдами в море Уэдделла.

Рисунок 9. Слева: увеличенный вид следа от Пылающего Корабля. Справа: фотография корабля Шеклтона, вмерзшего в лёд в Антарктиде, сделанная Херли в 1915 году.

Башня из возведения в степень

Начните с любого комплексного числа z и многократно возводите его в степень.

z, z^z, z^{z^z},...

Мы можем выразить это как итерацию. Начнем с y₀ = 1 и далее

y_k+1=z^y_k

Если z слишком велико, эта итерация уйдет в бесконечность. Для некоторых же z она будет сходиться к конечному пределу. Например, если z = sqrt(2), то y_k сходится к 2. Самый интересный случай, когда y_k зацикливается. Например, если z взять близким к 2.5+4i, то цикл будет иметь длину 7.

2.4684 + 4.0754i
-0.6216 + 0.3634i
0.2603 - 0.0184i
1.4868 + 0.3613i
-3.4877 + 6.1054i
7.7632e-06 - 2.6617e-06i
1.0000 + 0.0000i
2.4684 + 4.0755i

Эта длина цикла является основой фрактала башни возведения в степень. На рисунке 10 показан общий фрактал:

Рисунок 10. Фрактал башни возведения в степень

Увеличение в 10⁵ раз и изменение цветовой карты дает изображение, показанное на рисунке 11.

Рисунок 11. Подробный вид башни возведения в степень

Метод Ньютона

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

Итерация схожа. Выберите функцию f(z) с производной f′( z). Начните с z₀ и выполняйте итерации

z_к+1 = z_k- f(z_k) /f′(z_k)

В конечном итоге это приведет к нулю функции f. Подсчитайте количество итераций, необходимых, чтобы сойтись.

Есть много функций (и цветовых карт) на выбор. Мой любимый - кубический многочлен:

f(z) = z³ – 2z– 5

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

Рисунок 11. Итерация Ньютона по z³ – 2z– 5

На рисунке 12 показано глобальное поведение метода Ньютона, ищущего нули функции

f(z) = tan(sin(z)) − sin(tan(z))

Такая функция имеет бесконечно много нулей и неограниченную первую производную.

Рисунок 12. Подробный вид из итерации Ньютона по tan(sin(z)) – sin(tan(z))

Функция на рисунке 13: f(z) = z*sin(1/z)

Наиболее заметные синие области окружают нули на уровне ± 1/?.

Рисунок 13. Подробный вид из итерации Ньютона по z*sin(1/z)

Многие из описанных здесь изображений были для меня новыми - я никогда их раньше не видел. Это стало возможным благодаря интерактивным экспериментам с графическим процессором и использовании среды MATLAB для расчета фракталов.