Визуализация решений логистического уравнения, множеств Мандельброта и Жюлиа в Excel и Visual Basic Текст научной статьи по специальности «Математика»

Миллер Дейв Бриджес Ричард

ВИЗУАЛИЗАЦИЯ РЕШЕНИЙ ЛОГИСТИЧЕСКОГО УРАВНЕНИЯ, МНОЖЕСТВ МАНДЕЛЬБРОТА И ЖЮЛИА В EXCEL И VISUAL BASIC

В настоящей статье рассматриваются способы изображения множеств с помощью различных программным средств. Сначала мы1 изучаем логистическое уравнение, используя современные средства управления рабочим листом (Excel worksheet) - «движки» и макросы1, управляемые кнопками. Затем те же средства применяются к множествам Мандельброта и в заключение рассматриваются образы1 множеств Мандельброта и Жюлиа, полученные с выеоким разрешением благодаря использованию менее известных но гораздо более эффективным алгоритмов.

ЛОГИСТИЧЕСКОЕ УРАВНЕНИЕ

В 1976 году Р. Мэй [1] настоятельно рекомендовал, чтобы знакомство с логистическим уравнением происходило на ранней стадии математического образования. Правда тогда только калькуляторы были доступны широкой публике.

Теперь же, имея изощренные программные средства, мы можем исследовать это уравнение. Логистическое уравнение (1), названное так УегИиМ'ом в 40-х годах 19-го века [3] представляет собой модель изменения численности популяции, где рп - относительная величина популяции р = Р /И, а Р - абсолютная величина

' п п 7 п

Iterations 0 to 500

популяции к моменту п и N - максимальная популяция. Р. Мэй использовал эквивалентную модель (2), где А - константа.

Рп+1 = Рп + Г • Рп(1 - Рп) (1)

= A ■ Xn(1 - xn)

(2)

Диаграмма 1. Хаотическое Диаграмма 2. Цикл с пери-поведение при А = 3.701. одом 7 при А = 3.702.

Существует обширная литература, посвященная поведению этой модели (Peitgen et al, 1992). Мы предлагаем здесь три метода, использующих рабочий лист Excel'а. Заинтересовавшиеся могут использовать эти методы для самостоятельного изучения модели.

МОДЕЛИРОВАНИЕ В EXCEL

При открытии рабочего листа Excels с именем logisticterms (Miller, 2001) подсчитывают...................................... ся 500 итераций уравнения

...................................... (2) и демонстрируется гра-

..................................... фик зависимости между xn

и n. Значение параметра А связано с движком и мо...................................... жет меняться от 0 до 4 с

шагом 0.001. Другой движок связан с начальным ...................................... значением х0, которое меняется от 0 до 1 с шагом 0.001. Кнопка, связанная с макросом, позволяет пользователю увидеть больше итераций, сгруппированных по 500, и тем са-

мым увидеть, стабилизируются ли итерации при больших значениях п. Таким образом можно изучать поведение модели при разных значениях А и х0. Диаграммы 1 и 2 показывают эффект изменения А от 3.701 до 3.702 (х0 = 0.5)

Рабочий лист Excel'а с именем logisticcobweb (Miller, 2001) представляет те же результаты в другой форме - в форме паутинообразной диаграммы. Здесь строятся графики y = Ax(1-x) и y = x, а также отмечаются и объединяются точки (xn, xn+1), (xn+1, x 1), n = 0,1,2,... Это придает диаграмме вид паутины.

Рабочий лист Excel'а с именем logisticconsecutive (Miller, 2001) демонстрирует третий способ представления результатов. Те же результаты показаны на диаграммах 5 и 6.

МНОЖЕСТВО МАНДЕЛЬБРОТА

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

Множество Мандельброта определяется с помощью функции f(z,c) = z2 + c, где z и с - комплексные числа.

Если z = zr+izi и с = с+'с, то

f (z, c) = (zГ - zi + cr) + i(2zrzi + c) (3)

Последовательность (орбита) значений z получается посредством итераций

z1 = f(z0, c), z2 = f(zp

c), ... ,

zn = f(zn-1, c), ...

Очевидно, орбита зависит от начального значения z0 и от значения константы с. При этом обнаруживаются два явления. Либо zn стремится к бесконечности (обычно с большой скоростью), либо стремится к фиксированно-

Диаграмма 4. Цикл с периодом 7 при А = 3.702.

му числу или к предельному циклу, состоящему из конечного числа точек.

Множество Мандельброта - это множество тех значений с, для которых при г0 = 0 последовательность гп не стремится к бесконечности. Используя возможности Ехсе1'а, можно изобразить любую часть множества Мандельброта. На диаграмме 7 показан рабочий лист «квадрат», состоящий из 441 квадратной ячейки, где каждый квадрат представляет точку комплексной плоскости. Исходная область включает множество Мандельброта, где гг изменяется от -1.5 до 0.5, а ^ - от -1 до 1. Желтый (светло-серый) квадрат представляет число -0.1+0.9/, черный - 0+0/, а красный (темно-серый) - 0.5+/.

Применим итерации, положив с = 0.5+/. Таблица 1 показывает, что модуль г превосходит 5000 после 6 итераций.

На рабочем листе подсчитываются 50 итераций для каждой точки из заданного диапазона и определяется число ите-

1-1 Iterations 0 to 500

Диаграмма 6

Диаграмма 3. Хаотическое поведение при А = 3.701.

Диаграмма 5

№ итерации Zr çz ç

0 0 0 0

1 0.5 1 1.12

2 -0.25 2 2.02

3 -3.44 0 3.44

4 12.32 1 12.36

5 151.19 25.63 153.35

6 22203.04 7752.05 23517.43

Таблица 1

раций, при котором модуль z превзойдет заданное число (здесь 5000). Предполагается, что если это не произойдет за 50 итераций, то не произойдет и при большем числе итераций. Числа 50 и 5000 выбраны произвольно. Результаты расчетов с этими параметрами приведены на диаграмме 8.

Для раскраски диаграммы 8 использованы 4 цвета (в статье - оттенки серого). Точки множества Мандельброта окрашены черным цветом. Точки, для которых расходимость (I znl >5000) обнаруживается при n<17, выкрашены красным (темно-серым); если же это происходит при 17<n<33 или 34<n<49, то - желтым (светло-серым) или белым, соответственно. Выбор 4-х цветов связан со свойством рабочего листа (Microsoft Excel 97), которое позволяет связывать цвет отдельной ячейки с выполнением определенных условий (таких условий может быть только 3, что дает 4 цвета).

ИЗОБРАЖЕНИЯ

С ВЫСОКИМ РАЗРЕШЕНИЕМ

Наиболее распространенные «стандартные» алгоритмы построения образов множества Мандель-брота плохо работают при больших увеличениях, необходимых для демонстрации самых интригующих явлений, связанных с такими множествами. Здесь Диаграмма 7

мы представляем программу, реализующую менее известный, но гораздо более эффективный алгоритм.

Множество Мандельброта - это множество значений с, для которых последовательность итераций zn не расходится при z0 = 0. Множество Жюлиа для заданного параметра с - это множество тех z0, при которых последовательность итераций zn не расходится. Таким образом существует только одно множество Мандельброта и бесконечно много множеств Жюлиа. С этими простыми определениями связано большое количество интереснейших явлений, которых мы лишь слегка касаемся в данной статье.

Чтобы создать изображение одного из этих множеств, нужно установить соответствие между массивом пикселей на экране компьютера и множеством комплексных чисел z. Условимся, что х-коор-дината на экране представляет вещественную часть числа z, а y-координата - мнимую часть. Соответствующий программный код для среды программирования Visual Basic 5 приведен в листинге 1. В других средах могут быть отличия в деталях, но принципы сохраняются.

Первая секция листинга содержит объявления переменных. Те, что участвуют в итерациях, объявлены как Double для сохранения точности, иначе изображения могут потерять детали при большом увеличении. Подпрограмма Form_Load игра-

I Яй ■

■ ■ ■!■ I.в.

■

_ЙГ.

Диаграмма 8

ет роль головной программы. Подпрограмма SetScreen инициализирует экран Visual Basic'a с началом координат в левом нижнем углу и координатами пикселей в диапазоне от 0 до srange. Подпрограмма SetFractal инициализирует переменные, связанные с создаваемым изображением: zrmid, zimid и zrange описывают прямоугольник на комплексной плоскости, соответствующий экрану, а palno определяет используемую раскраску. Остальные переменные будут использованы позже. Подпрограмма SetColours сохраняет в массиве palcol() цвета, которые будут использованы для окраски разных частей образа.

Остальные короткие функции и подпрограммы:

- NotPlottedYet проверяет, был ли изменен цвет пикселя по сравнению с начальным; такая проверка упрощает дальнейшую работу с массивом пикселей;

- PlotPixel окрашивает пиксель;

- CheckForStop необходима в многозадачной среде Windows, чтобы иметь возможность прервать работу программы, когда приходится достаточно долго ждать завершения работы; подпрограмма StopButton реализует этот стоп-сигнал.

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

Математическая часть задачи выполняется в ComputeColour. Здесь экран представляет часть комплексной плоскости параметра с. Числа zn = f(zn-1, c), n = 1, 2, 3, ..., z0=0, полученные из уравнения (3) хранятся в массивах zrorb() и ziorb() (для стандартного алгоритма это не явля-

ется необходимым). Цикл Do loop работает, пока либо число итераций, либо модуль zn не превысят заданных границ. В первом случае с принадлежит множеству Мандельброта, и точка окрашивается одним цветом (обычно черным), а во втором - не принадлежит и окрашивается другим цветом. В простейшем варианте (palno =1) все точки, не принадлежащие множеству Мандельброта, окрашены одним цветом. Однако, как будет показано далее, целесообразней сделать цвет зависящим от числа итераций, проделанных до того, как достигнуто предельное значение модуля zn (maxmodzsq). Это приводит к появлению множества цветных контуров вокруг множества Мандельброта. Простейший выбор - чередование цветов контуров (рисунок 1).

Выбор предельного числа итераций (maxiter) зависит от степени увеличения. Если maxiter мало, то при большом увеличении граница изображаемого множества будет искажаться. Значение maxiter = 256 годится только для небольшого увеличения, а maxiter= 1024 подходит для любого увеличения, хотя и замедляет работу программы.

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

Рисунок 1

Пусть ваш код правильно изображает множество Мандельброта и окружающие его контуры. Теперь измените код, чтобы был реализован Case 2 (двухцветная схема -palno =1). В этом случае все еще используется стандартный алгоритм, и можно увидеть, насколько он неэффективен (рисунок 2) по сравнению с более мощным алгоритмом, реализованным в Case 3. В Case 2 теряется тонкая структура границы изображаемого множества (рисунок 3).

В чем же заключается трудность и как ее преодолевает этот более мощный алгоритм? Дело в том, что множество Мандельброта М состоит из чрезвычайно тонких нитевидных фрагментов. Такие фрагменты исчезают в любом двухцветном изображении, если разрешение не очень велико. Решение состоит в том, чтобы оценить расстояние от любой точки с, не принадлежащей М до ближайшей точки этого множества. Если это расстояние меньше одного пикселя, то точку с экрана следует окрашивать цветом границы М. Метод оценки расстояния был опубликован в [5]. Он основан на глубоком анализе последовательности итераций в [2]. Они доказали, что расстояние от с до М меньше, чем 2zn|l°g|zn|/|zn|, где n - номер первой итерации, для которой |zn| превосходит некоторое заданное значение (в программе это maxmodzsq должно быть весьма большим), а zn - производная zn по с. Дифференцируя уравнение zk+x = z\ + с по с, получаем рекуррентное соотношение для z

zk+1 = 2zkzk +1, k = 0,1,2,..., z0 = 0 (4)

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

Возвращаясь к листингу 2, заметим, что для использования соотношения (4) необходимо хранить в памяти все значения zk, 0<k<n; они сохраняются в массивах zrorb() и ziorb(). Если оценка расстояния предусмотрена (переменная DEMflag), то нужно вычислять zn, помещая вещественную и мнимую части в zdr и zdi. Так как величина z' может быть

n

очень большой вблизи множества М, нужно быть уверенным в отсутствии переполнения (в этом причина сравнения с maxmodzsq). Величина delta (ширина пикселя) подсчитывается в конце подпрограммы SetFractal, так как она зависит от степени увеличения. Дополнительная переменная dfactor (поправочный множитель) обеспечивает тонкую настройку, так как результат очень чувствителен к величине delta. Если Mdist (оценка расстояния между с и М) не превышает эту величину, то переменной colour присваивается значение 1 (черный), что обозначает точку, удаленную от М не более, чем на ширину пикселя, но не принадлежащую М. Если

А А

.....

Рисунок 2 Рисунок 3

Рисунок 4

colour = 0 соответствует серому цвету, а colour =1 - черному (как в SetColours), то при достаточно большом увеличении в правильно выбранном месте видны малые копии М, окрашенные серым цветом, сквозь которые просвечивают черные нитевидные структуры (Cases 2-15 в SetFractal и рисунок 4). Пиксели, не принадлежащие М или окрестности границы М, можно окрашивать стандартным алгоритмом.

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

де оценки расстояния теперь нужно для получения рекуррентного соотношения дифференцировать по z0, а не по с.

zk+1 = 2zkzk = 0,1,2,..., zо = 1 (5)

Case 18 в SetFractal дает с большим увеличением изображение части множества Жюлиа вокруг точки z = с, где с взято из малой копии множества М, построенной в Case 15. Изменяя увеличение, можно увидеть, что множество Жюлиа для с, взятого из окрестности границы множества М, само выглядит как эта часть М.

В заключение нужно предупредить о некоторых ограничениях, присущих данной программе. Наиболее серьезное из них связано с невозможностью изменять степень увеличения интересующего нас участка простым щелчком по экрану. Написать соответствующий код несложно, однако он будет существенно зависеть от используемой среды. Менее серьезным является неэффективность вычисления множества Жюлиа - значения z быстро сходятся к предельному значению или предельному циклу, однако вычисления продолжаются пока не проделано maxiter итераций. Это может быть легко исправлено. И, наконец, отсутствуют какие-либо средства сохранения результатов, но если вы работаете в Visual Basic, то можете связать управляющую кнопку с подпрограммой, помещенной в конце листинга 3. Нажатие кнопки посылает изображение в Clipboard в виде bitmap. Имеется также большой простор для экспериментирования с цветовыми схемами; другая возможность появляется при установке значения palno = 3.

Литература.

1. May R.(1976) Simple mathematical models with very complicated dynamics, Nature, 261, pp.459-467.

2. Milnor J., Thurston W.(1989) Selfsimilarity and hairiness in the Mandelbrot set, Lecture Notes, Pure and Appl. Math.,114.

3. Peitgen H-O., Jurgens H. & Saupe D.(1992) Fractals for the Classroom: Part One, Introduction to Fractals and Chaos, (Spring er-Verlag).

4. Peitgen H-O., Richter P.H.(1986) The Beauty of Fractals, (Springer-Verlag).

5. Peitgen H-O., Saupe D. (1988) The Science of Fractal Images, (Springer-Verlag).

Приложение.

Listing 1

Dim ssize, sxmin, symin, srange, scol, sbackcol As Double Dim sx, sy, s, stop flag

Dim palno, mandelcol, nearmandelcol, notmandelcol, palcol(1024) Dim iter, maxiter, maxmodsq, maxdmodsq, delta, pi, DEMflag, j, fractalno Dim zrmin, zimin, zrange, cr As Double, ci As Double, zr As Double Dim zi As Double, t As Double

Dim zrorb(1024) As Double, ziorb(1024) As Double, zdr As Double Dim zdi As Double

Private Sub Form_Load() SetScreen SetFractal 1 SetColours

s = 256: stopflag = False

While s > smin And Not CheckForStop

s = s / 2

For sy = symin To symin + srange Step s

For sx = sxmin To sxmin + srange Step s If NotPlottedyet(sx, sy) Then

zr = zrmin + (sx - sxmin) * zrange / srange zi = zimin + (sy - symin) * zrange / srange scol = ComputeColour(zr, zi) PlotPixel sx, sy, scol End If Next sx

If CheckForStop Then Exit For Next sy

Wend End Sub

Public Sub SetScreen()

ssize = 3900: sxmin = 0: symin = 0: srange = 256: smin = 1 screen.Height = ssize: screen.ScaleHeight = -srange screen.ScaleTop = srange

screen.Width = ssize: screen.ScaleWidth = srange: MainForm.Show End Sub

Public Sub SetFractal(fractalno)

maxmodzsq = 1E+21: maxdmodzsq = 1E+60: pi = 4 * Atn(1) Select Case fractalno

Case 1: j = 0: zrmid = -0.75: zimid = 0: zrange = 3: maxiter = 64

palno = 0: DEMflag = False: dfactor = 1 End Select

delta = dfactor * 0.71 * zrange / srange: zrmin = zrmid - zrange / 2 zimin = zimid - zrange / 2 End Sub

Public Sub SetColours()

sbackcol = vbCyan: screen.BackColor = sbackcol mandelcol = RGB(128, 128, 128): nearmandelcol = vbBlack notmandelcol = vbWhite

palcol(0) = mandelcol: palcol(1) = nearmandelcol palcol(2) = notmandelcol Select Case palno Case 0

For c = 3 To maxiter

palcol(c) = QBColor(1 + (c Mod 15)) Next c Case 1

For c = 2 To maxiter

palcol(c) = notmandelcol Next c Case 2

m = 12

For c = 2 To m

s = Log(c) / Log(m)

palcol(c) = RGB(256 * s, 256 * s, 256 * s) Next c

For c = m + 1 To maxiter

palcol(c) = notmandelcol Next c End Select

End Sub

Public Function NotPlottedyet(sx, sy)

If screen.Point(sx, sy) = sbackcol Then NotPlottedyet = True _ Else NotPlottedyet = False

End Function

Public Sub PlotPixel(sx, sy, scol)

screen.PSet (sx, sy) , palcol(scol)

End Sub

Public Function CheckForStop() DoEvents

If stopflag Then CheckForStop = True Else CheckForStop = False

End Function

Private Sub StopButton_Click() stopflag = True

End Sub

Listing 2

Public Function ComputeColour(zr As Double, zi As Double) If j = 0 Then cr = zr: ci = zi: zr = 0: zi = 0 iter = 0: zrorb(0) = zr: ziorb(0) = zi Do

iter = iter + 1

t = zr * zr - zi * zi + cr: zi = 2 * zr * zi + ci: zr = t zrorb(iter) = zr: ziorb(iter) = zi 'Comment: not needed for simple algorithm Loop Until iter >= maxiter Or zr * zr + zi * zi > maxmodzsq If iter >= maxiter Then ComputeColour = 0 If zr * zr + zi * zi > maxmodzsq Then ComputeColour = iter 'Comment: standard simple algorithm stops here If DEMflag Then

zdr = 0: zdi = 0: i = 0: If j = 1 Then zdr = 1 While zdr * zdr + zdi * zdi < maxdmodzsq And i <= iter - 1 t = 2 * (zdr * ziorb(i) + zdi * zrorb(i)) zdr = 2 * (zdr * zrorb(i) - zdi * ziorb(i)) If j = 0 Then zdr = zdr + 1 zdi = t: i = i + 1

Wend

If zdr * zdr + zdi * zdi > maxdmodzsq Then ComputeColour = 1

Else

modz = Sqr(zr * zr + zi * zi) modzd = Sqr(zdr * zdr + zdi * zdi) Mdist = 2 * modz * Log(modz) / modzd If Mdist < delta Then ComputeColour = 1 End If

End If End If End Function

Listing 3

Case Case Case Case Case Case Case Case Case Case Case Case Case Case Case Case

Case

2: j = palno = 3: j = palno = 4: j = palno = 5: j = palno = 6: j = palno = 7: j = palno = 8: j = maxiter 9: j = maxiter 10: j maxiter 11: j maxiter 12: j maxiter 13: j maxiter 14: j = maxiter 15: j = maxiter 16: j = maxiter 17: j = ci = -0 dfactor 18: j = cr = 0. DEMflag

zrmid = -0.75: zimid = 0 DEMflag = False: dfactor zrmid = -0.75: zimid = 0

DEMflag = True: dfactor = 1

zrmid = -0.75: zimid = 0

DEMflag = True: dfactor = 1 zrmid = -0.75: zimid = 0

zrange = 3 = 1 zrange = 3

zrange = 3

= 3

maxiter = 64

maxiter = 64

maxiter = 64

maxiter = 64

zrange

DEMflag = True: dfactor = 1

zrmid = 0: zimid = -1: zrange = 0.1: maxiter = 64 DEMflag = True: dfactor = 1

zrmid = -0.05: zimid = -1: zrange = 0.05: maxiter = 256 DEMflag = True: dfactor = 1

zrmid = 0.2849758: zimid = -0.0112716: zrange = 1 = 64: palno = 1: DEMflag = True: dfactor = 0.5 0: zrmid = 0.2849758: zimid = -0.0112716: zrange = 0.1 = 128: palno = 1: DEMflag = True: dfactor = 0.5

= 0

= 256: palno = 1: DEMflag = True: dfactor

= 0

= 0

= 0

zrmid = 0.2849758: zimid = -0.0112716

zrmid = 0.2849758: zimid = -0.0112716

= 1024: palno = 1: DEMflag = True: dfactor = 0.5

zrmid = 0.2849758: zimid = -0.0112716

= 1024: palno = 1: DEMflag = True: dfactor = 0.5

zrange = 0.01 = 0.5 zrange = 0.001

zrange = 0.0001

zrange = 0.00001

zrmid = 0.2849758: zimid = -0.0112716: = 1024: DEMflag = True: palno = 1: dfactor = 0.5 0: zrmid = 0.2849758: zimid = -0.0112716: zrange = 0.000001 = 1024: DEMflag = True: palno = 1: dfactor = 0.5 0: zrmid = 0.2849757: zimid = -0.0112718: zrange = 0.0000005 = 1024: DEMflag = True: palno = 1: dfactor = 0.5 1 : zrmid = 0 : zimid = 0 : zrange = 3 : cr = 0 : ci = 1 = 256: DEMflag = True: palno = 1: dfactor = 0.5 : 1: zrmid = 0: zimid = 0: zrange = 3: cr = 0.2849757 .0112718: maxiter = 1024: palno = 1: DEMflag = True = 0.5

= 1: zrmid = 0.2849757: zimid = -0.0112718: zrange = 0.0001 2849757: ci = -0.0112718: maxiter = 1024: palno = 1 = True : dfactor = 0.5

Private Sub copybutton_Click() ClipBoard.Clear

Clipboard.SetData screen.Image End Sub

Richard Bridges,

King Edwards School, Birmingham, UK. Dave Miller,

Department of Education, Keele University, UK.

© Наши авторы, 2001. Our authors, 2001