*Прим. В. Шипкова: отдельная благодарность Венедикту Ли за помощь в оформлении главы.
23-й час
Набор Мандельброта
В предыдущей главе мы завершили рассмотрение функций рисования и менеджеров геометрии объектов библиотеки Tk. В этой главе мы применим на практике знания, полученные в предыдущих шести главах, и создадим завершённое приложение для расчёта и вывода наборов Мандельброта (Mandelbrot sets), а также связанных с ними наборов Джулия(Julia sets).
Это было в 1985 г., только вышел августовский номер Scientific American. На следующее утро большинство инженеров компании, в которой я работал, принесли с собой написанные на скорую руку простенькие программы Мандельброта, а к концу недели у каждого уже была своя такая программа. В следующий понедельник тем, кто пытался выполнить какую-либо полезную работу на любом из семи компьютеров компании, оставалось только чертыхаться, поскольку компьютеры зависали из-за того, что на них несколько программ одновременно вычисляли наборы Мандельброта.
Я думаю, что та же самая картина повторилась во всех компаниях страны, где работали программисты, хотя в нашей компании концентрация инженеров-компьютерщиков на единицу простого населения была, пожалуй, самой высокой. В то время наиболее мощной машиной с графическим интерфейсом, к которой я имел доступ, была Intelecolor, способная выгружать информацию со скоростью 9 600 бод (1 бод = 1 бит/с). И даже на этом чуде техники вывод набора Менделброта занимал не менее 5 минут. А такому давно вымершему динозавру, как Gould NP1, требовалось несколько часов для выполнения той же непростой работы. Компьютер буквально дымился, и это при том, что палитра состояла только из восьми выводимых цветов. Я пытался уговорить шефа подключить сетевую карту на мой терминал, но он сказал: "Обойдешься, я знаю, зачем тебе это нужно!" И он был прав.
Я полагаю, что в течение нескольких месяцев после выхода этого номера Scientific American 90—95% национальных компьютерных ресурсов было потрачено на выведение наборов Мандельброта. Изучая материал данной главы, Вы сможете испытать то же удовольствие, которое испытывали мы.
Если вычисление наборов Мендельброта реализовать на языке С, то программа будет выполняться относительно быстро. Конечно, следует учесть, что в 1985 г. самый мощный Gould NP1 имел быстродействие, как мы недавно подсчитали с друзьями, аж 33 МГц. И это был один из мощнейших четырёхпроцессорных компьютеров. Поэтому программы вычисления наборов Мандельброта можно было писать только на С, иначе программы выполнялись бы сутками. Но в наше время, когда быстродействие большинства домашних персональных компьютеров перевалило за 400 МГц, мы можем позволить себе реализовать такую программу полностью на языке Python. Безусловно, что в таком виде программа будет выполняться медленнее, чем если бы она была написана на С. Для компьютеров большинства читателей время выполнения нашей программы, показанной далее в этой главе, составит где-то около 5 минут. Разобравшись в принципах решения задачи в Python, Вы сможете разыскать в Internet более быстрый вариант программы. Для создания программ, подобных той, что приведена в качестве примера в этой главе, могут оказаться полезными дополнительный модуль NumPy, который можно загрузить с сервера http://numpy.sourceforge.net/, и модуль вывода окон в Windows wxPython (http://alldunn.com/wxPython/). Рассмотрение этих модулей выходит за рамки данной книги.
Я уверен, что Вы уже видели в своей жизни графику, построенную на основе наборов Мандельброта. Если нет, то один пример показан на рис. 23.1.
Это изображение было создано лично Беноитом Мандельбротом (Benoit Mandelbrot), описавшим в 1980 г. математический подход вычисления наборов его имени, и подарено им фантасту Артуру С. Кларку (Arthir С. Clarke). Др. Кларк воскликнул, что наборы Мандельброта вмещают в себе больше, чем Вселенная, или, по крайней мере, являются её отражением. Этой идеей он поделился со Стефаном Хавкином (Stephen Hawking). В поддержку др. Хавкин привел пример минимальной длины Планка. (Planck Length — 1.616x10 33 см, меньше которой не может быть ни одна элементарная частица. Меньше кварка - это одна из фундаментальных констант Вселенной.) Тем не менее в математической модели Мандельброта могут существовать частицы меньшего размера, по крайней мере виртуально. Единственный, но безусловный предел — это минимально допустимое значение числа с плавающей запятой на Вашем компьютере.
*Прим. Венедикта Ли: Здесь, после сверки с оригиналом, заменены рис. 23.1, 23.3 и добавлены отсутствовавшие фрагменты текста и рисунки 23.4-23.9.
Рис. 23.1. Графическое отображение набора Мандельброта
Математическая сторона вопроса
Формула вычисления набора Мандельброта чрезвычайно проста. В общем виде она выглядит так:
z=z2+с
Некоторые соглашения остались невидимыми в уравнении. Значения z и с — это комплексные числа. Важным моментом также является то, что данное уравнение применяется для вычисления цвета каждого пикселя изображения. Для наборов Мандельброта принято, что изображение формируется в системе координат размером -2,0x2,0 по осям Х и У с точкой (-2,0, -2,0) в верхнем левом углу. Напомним, что комплексные числа состоят из двух компонент — реального и мнимого. Графически принято отображать комплексное число как точку в системе координат, где по оси X откладывается реальная составляющая, а по оси Y— мнимая. Любая точка в этом пространстве — комплексное число, такое как с. В Python создать комплексное число очень просто: c=complex(x,y). Для наглядности система координат комплексных чисел показана на рис. 23.2.
Как в любой другой системе координат, оси X и Y в системе комплексных чисел бесконечны. Но нам для графического отображения набора Мандельброта нет необходимости вычислять все множество комплексных чисел, а достаточно ограничиться только набором точек, выводимых на экран. В программе, представленной ниже в этой главе, для показа графики я выбрал холст размером 400x400 пикселей, т.е. в нашей системе координат по осям X и Y отложено по 400 точек. Число точек на видимой части оси координат называется разрешением. Таким образом, в нашей системе разрешение по оси Х и по оси Y равняется 400, а в целом на экран выводится 160000 точек.
Рис. 23.2. Система комплексных чисел
Увеличивать разрешение без видимой нужды не рекомендую, ведь уравнение z = z2 + с многократно выполняется для всех выводимых пикселей. Но мы до сих пор не знаем,, что такое z. Это значение можно представить как третье измерение — расстояние от точки до плоскости, в которой лежит набор Мандельброта. Первое, что мы делаем, — это определяем для каждого пикселя, принадлежит ли он набору. Если да, то пиксель выводится чёрным цветом. Если нет, то цвет пикселя определяется по тому, насколько он отдалён от плоскости набора Мандельброта. Но как определить, принадлежит ли пиксель (или точка) этому набору? Для каждой точки мы решаем указанную выше формулу и следим, как увеличивается расстояние от плоскости (0, 0j). Затем применяем следующий алгоритм.
Точки, которые никогда далеко не удаляются от плоскости набора, принадлежат этому набору и окрашиваются в чёрный цвет.
Точки, которые стабильно удаляются от плоскости набора в бесконечность, не принадлежат набору, а их цвет определяется скоростью удаления.
Как Вы знаете, это недостойное занятие, издеваться над компьютером, оперируя такими понятиями, как "никогда" и "бесконечность". Поэтому мы задаём некоторые произвольные пределы, в которых определяется цвет пикселей. Первый "верстовой столб" на пути от плоскости набора установим по длине выводимой части оси координат. Так как значения на оси X изменяются в пределах от -2,0 до 2,0, то первое граничное значение z можно сделать равным 4,0. Но как же не хочется иметь дело с отрицательными значениями. Проще будет учитывать абсолютное расстояние от точки до плоскости, которое в этом случае будет изменяться от 0 до 2,0. Для этого после вычисления расстояния по формуле z=z2+с используется функция abs(). Обратите внимание, что функция abs() в качестве параметра принимает числовые значения любого типа, независимо от того, был импортирован модуль math или cmath.
Для обычных чисел функция abs() просто удаляет значение знака. Так, abs(-4) возвращает 4, a abs(3.14159) — 3.14159. Вычисление абсолютного значения комплексного числа выполняется сложнее. Сначала мнимая и реальная составляющие возводятся в квадрат, затем суммируются, а от результата берется корень. Например, при вычислении абсолютного значения комплексного числа (-0.5, 1j) выполняются следующие действия:
sqrt((-0.5*-0.5)+(1*1))
sqrt((0.25)+(1))
sqrt(1.25)=1.11803398875
Каждую из показанных выше строк можно выполнить в окне интерпретатора по отдельности, в любом случае результат будет один и тот же. Только помните, что в данном случае мы имеем дело не с комплексным числом, а с его имитацией. А теперь попробуем вычислить абсолютное значение действительно комплексного числа. Введите следующие строки в окне интерпретатора:
с=0.5+lj
abs(c)
В результате получим то же самое число: 1.11803398875 (Если у Вас получилось что-то другое, то либо Вы работаете с весьма древней операционной системой с устаревшим способом представления чисел с плавающей запятой, либо установка Python на Вашем компьютере прошла катастрофически неудачно.) В языке С, где нет встроенных комплексных чисел, для вычисления абсолютного значения комплексного числа потребуется самостоятельно написать код, что усложнит всю программу. К счастью, в Python приведённое выше выражение можно использовать для возвращения искомого значения таким, как оно есть, поэтому я избавлен от необходимости рассказывать Вам о том, как много проблем и побочных явлений ожидает того, кто пытается работать с комплексными числами в других языках программирования. (Раз уж я взялся за написание этой книги, то чувствую своим долгом убедить Вас, что Python — самый лучший язык программирования!)
*Прим. В. Шипкова: и я, и я, и я того же мнения. (с) Ослик. "Винни-Пух и все все все".
Ну а теперь рассмотрим псевдокод нашей программы, чтобы понять её структуру и принципы работы.
for всех пикселей по оси X:
for всех пикселей по оси X:
с=complex(X,у)
for всех пикселей по оси Z:
if abs(z)>2.0:
break
z=(z**2)+с
присвоить пикселю z чёрный цвет, если цикл не был прерван инструкцией break в противном случае вычислить цвет пикселя по числу итераций
Точки по краям выводимого поля выходят за пределы набора после первой итерации. Например, abs(-2.0+-2.0j) сразу возвращает 2.82842712475. Небольшая программа, показанная в листинге 23.1, позволит Вам проследить динамику изменения значения z для каждой точки в пределах выводимого поля комплексных чисел.
Листинг 23.1. Программа ctest.py
import sys
import string
c=-0.5+1j
z=0j
if len(sys.argv)>2:
r=string.atof(sys.argv[1])
i=string.atof(sys.argv[2])
c=complex(r,i)
print "initial value of c", abs(c)
for i in range(64):
print "step",i,"value of z", z, \"distance (abs()) of z",abs(z)
if abs(z)>2.0:
break
z =(z**2)+c
print "color z", iВ строках 6, 7 по умолчанию задаются те значения переменных с и z, которые мы уже использовали в предыдущих примерах. Вывод программы показан в листинге 23.1.
Листинг 23.2. Вывод программы ctest.py
initial value of с 1.11803398875
step 0 value of z Oj
distance (abs)) of z 0.0
step 1 value of z (-0.5+lj)
distance (abs()) of z 1.11803398875
step 2 value of z (-1.25+0j)
distance (abs()) of z 1.25
step 3 value of z (1.0625+lj)
distance (abs()) of z 1.45907719124
step 4 value of z (-0.37109375+3.125j)
distance (abs()) of z 3.14695655694
color z 4
Очевидно, что точки, сосредоточенные вокруг центра (0,0j), будут окрашены в чёрный цвет, а точки по краю изображения будут окрашены в тот цвет, который мы определим как "самый далекий". Можно вывести вполне очевидный общий принцип, что по мере удаления от центра точки с большей скоростью улетают вдаль от плоскости набора Мандельброта, а точки, приближенные к центру, являются членами этого набора. Трудно сказать, действительно ли точки, не входящие в набор, удаляются в бесконечность или стремятся к какому-либо своему пределу. Тем, кому это интересно, следует заняться теорией хаоса и воспользоваться некоторыми её концепциями, в особенности концепцией атракторов (центров притяжения, к которым стремятся элементы стохастической системы, предоставленные сами себе). Но нас сейчас эти вопросы не интересуют.
Выбор цвета для пикселей основывается на числе итераций, необходимых для того, чтобы точка вышла за пределы набора Мандельброта. В листинге 23.2 пикселю был присвоен цвет под номером 4. Мы будем использовать те же цветовые схемы, задаваемые функциями модуля colormap, которые уже использовали при создании анимации в предыдущей главе. Каким будет 4-й цвет, зависит от выбранной цветовой схемы. Для функции, выбираемой по умолчанию, это будет цвет "#F2620C". В листинге 23.3 показана измененная программа ctest2.py, которая присваивает полученный цвет атрибуту фона холста.
Листинг 23.3. Программа ctest2.py
import sys
import string
from colormap import *
from Tkinter import *
cmap = SetupColormap0(64)
c=-0.5+1j
z=0j
if len(sys.argv)>2:
r=string.atof(sys.argv[1])
i=string.atof(sys.argv[2])
c=complex(r,i)
print "initial value of c", abs(c)
for i in range(64):
print "step",i,"value of z", z, \"distance (abs()) of z",abs(z)
if abs(z)>2.0:
break
z = (z**2) + c
print "color z", i, cmap[i]
root=Tk()
root["background"]=cmap[i]
root.mainloop()Программа вычисления набора Мандельброта
Изображение на рис. 23.1 было получено с помощью программы, код которой показан в листинге 23.4.
Листинг 23.4. Программа mandelbrot.py
*Прим. В. Шипкова: код программы - 610 строк. Поэтому текст не привожу. Код, который был приведён здесь содержал более 100 ошибок(!!!). Длина его была 709 строк. Скачать файл можно здесь. Здесь и далее ссылки на скачку одного и того же архива zip.
Код, безусловно, слишком длинный, чтобы обсуждать его в деталях. Мы только в общих чертах рассмотрим наиболее важные части программы. На рис. 23.3 показано окно программы в том виде, в каком оно открывается сразу после запуска. Можно установить альтернативный код кнопок, как показано на рис. 23.4, если в меню Файл выбрать команду Значки на кнопках.
Класс Msize, опредёленный в строках 10—55, позволяет выбрать для показа часть изображения и устанавливать разрешение осей X и Y. В результате мы можем управлять масштабом вывода изображения. Чтобы выбрать часть изображения для просмотра под увеличением, следует воспользоваться полосами прокрутки.
Рис. 23.3. Окно программы mandelbrot.py, которое открывается по умолчанию
Рис. 23.4. Альтернативный вид отображения кнопок со значками
Рис. 23.5. Подготовка к увеличению в 5.55 раз
Щелкните по кнопке Считать, чтобы посмотреть выбранный участок в увеличенном виде, как показано на рис. 23.6
Рис. 23.6. Завершенное вычисление
Все приложение реализовано как один класс Mandel, что облегчает возможность его модификации с помощью наследования, например для получения наборов Джулия, которые мы рассмотрим вкратце в следующем разделе. Список im метода __init__() в строке 81 содержит имена графических файлов, используемых для быстрого вывода исходных наборов Мандельброта при запуске приложения. Если бы программа вынуждена была вычислять этот набор, то её запуск растянулся бы на несколько минут. Как Вы видели, с помощью этого изображения можно выбрать для вычисления только некоторую часть набора, чтобы сократить время и не вычислять значения для всех пикселей.
Метод colorpoll(), опредёленный в строках 363-379, обслуживает прокручиваемые списки количества цветов и цветовой схемы. Выбрав значение в списке, Вам не нужно вводить его каким-либо способом в программу. Каждые 250 мс метод самоактивизируется и считывает текущие установленные значения, которые затем передаются в переменные-члены self.ncolor и self.сmар.
Методы xliner(), yliner() и zliner() (строки 259-346) используются для выделения части общего плана для увеличенного просмотра. Методы xliner() и yliner() прорисовывают вертикальную и горизонтальную линии в соответствии с положением бегунков на верхней и левой полосах прокрутки, отмечая в точке пересечения центр квадрата выделения. Метод zliner() устанавливает размер квадрата выделения по позиции бегунка на нижней полосе прокрутки.
Основная нагрузка при выполнении программы ложится на метод calc(), опредёленный в строках 399-473. При исходном разрешении 400 точек на ось этому методу приходится вычислять цвет для 160000 пикселей экрана. Для определения цвета одних пикселей достаточно одной итерации, тогда как для всех пикселей, входящих в набор, вычисление повторяется как минимум 64 раза или даже больше, в зависимости от того, какое значение присвоено переменной dwell.
Большинство остальных методов класса Mandel достаточно просты и понятны. Хотя код в строках 517-549 требует дополнительного объяснения.
Гастон Маурис Джулия (Gaston Maurice Julia), служивший солдатом во время Первой мировой войны, где-то незадолго до 1918 г. описал математическую матрицу, которая теперь носит его имя. Хотя Джулия был известен среди математиков своего времени, его труды оставались забытыми до тех пор, пока Беноит Мендельброт не описал свои наборы. По крайней мере частично, как признает и сам Мандельброт, он базировался на расчетах Джулия. Вас может удивить, что наборы Джулия были открыты раньше. Но ещё больше Вас должен удивить тот факт, что этот ученый вообще сумел открыть фрактальные объекты, не только не имея современного мощного компьютера, но даже не представляя, что это такое.
В чём же состоит разница между набором Мандельброта и набором Джулия. В случае с наборами Мандельброта мы последовательно выбираем все точки на плоскости, преобразовываем их координаты х и у в комплексное число и затем определяем, стремится ли эта точка улететь в бесконечность при выполнении ряда итераций базового уравнения. В случае с набором Джулия мы делаем почти то же самое, только в качестве значения с выбираем одну точку, принадлежащую набору Мандельброта. Мы также последовательно переходим от точки к точке комплексного пространства и определяем степень стремления z к бесконечности, но значение с при этрм остаётся константным. Координаты х и у, заданные точкой с, определяют исходное значение переменной z, тогда как в наборе Мендельброта исходное значение всегда (0, 0j). Таким образом, проанализировав все вышесказанное, нетрудно понять, что ту же программу, которую мы использовали для вычисления наборов Мандельброта, после небольшой модификации можно настроить на вычисление наборов Джулия. Тем более, что для выбора исходной точки с нам в любом случае нужен набор Мандельброта. На рис. 23.7 показан набор Джулия, вычисленный для исходной точки из набора Мандельброта с координатами (-0.79,-0.15j).
Как и при работе с наборами Мандельброта, в окне программы мы можем выбрать часть изображения и увеличить её. Пример выбора части изображения для увеличения показан на рис. 23.8, а сам увеличенный участок на рис. 23.9.
Рис. 23.7. Набор Джулия, вычисленный по точке (-0.79, -0.15j)
Рис. 23.8. Подготовка к увеличению
Рис. 23.9. Результат увеличения
Поскольку предыдущую программу mandelbrot.py мы выполнили практически полностью как один класс Mandel, то теперь можем унаследовать от него новый класс Julia. В листинге 23.5 показан код программы Julia.ру.
Листинг 23.5. Программа julia.py
*Прим. В. Шипкова: длина приведённого листинга - 269 строк. Столько же строк в оригинальном файле. Скачать файл можно здесь.
*Прим. Венедикта Ли: В этом архиве нет программы julia.py, а есть julia2.py. Её выполнение даёт те же результаты, что приведены в оригинале и показаны выше на этой странице.
Немного прокомментируем код программы julia.py. Главное, что нам нужно сделать в этой программе, — это изменить функциональность исходного класса Mandel. Это достигается замещением некоторых методов родительского класса в дочернем. Изменения коснулись методов calc(), clicker() и popJulia(). Эти методы теперь выполняются иначе или не выполняются вообще. Прежде всего мы изменили надпись и назначение кнопки Выход. На месте этой кнопки в окне приложения Julia Cruiser теперь появилась кнопка Закрепить с, благодаря которой закрепляются координаты точки с, после чего появляется возможность с помощью верхней и боковой полос прокрутки перемещать центр квадрата выделения, оставляя неизменными координаты исходной точки. Если не закрепить точку с, то, перемещая бегунки полос прокрутки, мы сместим исходную точку в неизвестном направлении, в результате чего не увеличим масштаб отображения, а получим новый непредсказуемый набор. Для закрепления исходной точки используется метод lock(), опредёленный в строках 32-43. Вычисление набора Джулия (строка 517 и далее в листинге 23.4) начинается с установки центральной точки и её закрепления. Собственно вычисление набора происходит в строках 529-532. (Обратите внимание, что эти строки не переписываются в программе julia.ру, а импортируются из класса Mandel.)
В данной главе Вы познакомились с наборами Мандельброта и Джулия и принципами их вычисления. Кроме того, на примерах окон приложений мы на практике применили почти все графические объекты библиотеки Tkinter для создания интерфейса GUI. В следующей главе мы поговорим об интерфейсе CGI, используемом при программировании для Web, рассмотрим несколько примеров программных кодов и обсудим ряд философских проблем.
Помимо загрузки процессора на максимальную мощность, чем ещё полезны наборы Мандельброта?
Если Вы бесчувственны к красоте или не так чувствительны, как моя жена ("Боже, что это? Как прекрасно"), то, возможно, Вас приведёт в восторг тот факт, что фрактальная геометрия и теория хаоса используются довольно широко, начиная от прогнозирования погоды и заканчивая сжатием графических файлов.
Вы не упоминали до сих пор о фрактальной геометрии. Что это такое?
Этот термин и направление в математике были предложены Беноитом Мандельбротом. С его точки зрения фрактальная геометрия — это геометрия живой природы. Крупные объекты являются подобием мелких объектов, из которых они созданы. Представьте себе береговую линию материка, как она видна из космоса. Теперь будем увеличивать изображение, как мы это делали в окне приложения Mandelbrot Cruiser. Чем ближе мы приближаемся, тем меньшая часть береговой линии нам видна, но выглядит она все так же, как и из космоса. Тонну материалов по фрактальной геометрии и её применении Вы сможете найти в Internet. Некоторые полезные ссылки Вы найдёте ниже, в разделе "Примеры и задания".
- С помощью какого уравнения вычисляются наборы Мандельброта?
а) Е = mс2
б) а = r2
в) z = z2 + с
г) F = G(m1 m2)/d2
Приблизительно сколько вычислений чисел с плавающей запятой нужно выполнить для определения набора Мандельброта в поле размером 400x400 пикселей с 32 установленными цветами?
а) 160000
б) 1000000000
в) 5120000
г) 2560000
Как много различных наборов Джулия содержит в себе один набор Мандельброта?
а) Один.
б) 160000.
в) По одному для каждой точки.
г) Бесконечное множество.
- в. Наборы Мандельброта вычисляются по формуле z = z2 + с.
г. Последнее значение — 2560000, пожалуй, ближе всего к правильному. Это не может быть число 160000, поскольку для каждой точки плоскости выполняется несколько итераций вычислений, неодинаковое для разных точек, так различно их отдаление от плоскости набора. Если это число равно 5120000, то все точки поля должны быть черными. Мы бы не стали вычислять такой набор, поскольку это скучно.
в. Для каждой точки набора Мандельброта можно рассчитать один набор Джулия. Такой набор Мандельброта называют каталогом наборов Джулия. Но справедливым будет также последний ответ — бесконечное множество, так как любой набор Мандельброта, по сути, содержит бесконечное число точек.
Интересную информацию, где упоминается как константа Планка, так и число гуголплекс, Вы найдёте на Web-странице по адресу http://www.informatic.unifrankfurt.de/~fp/Tools/GetAGoogol.html.
Дополнительную информацию о наборах Мандельброта Вы найдёте в следующих источниках.
- А. К. Dewdney. Computer Recreations. Scientific American 253, № 2 (August 1985), p. 16-24.
- Benoit B. Mandelbrot. The Fractal Geometry of Nature. — San Francisco: W. H. Freeman and Company, 1977, 1982, ISBN 0-7167-1186-9.
- htpp://linas.org/art-gallery/escape/escape.html.
Общепризнанной "библией" по компьютерной графике считается книга Computer Graphics: Principles and Practice, Second Edition in С by James D. Foley, Andries van Dam, Steven K. Feiner and John F. Hughes. — New York: Addison-Wesley, 1990.
