Легко видеть, что ряд расходящийся, числа увеличиваются, значит (0.5, 0.5) не попадает в множество.
При с = (0.25, 0.25) = 0.25 + 0.25j, имеем следующий ряд:
Z0 = 0
Z1 = Z0 + c = 0.25+0.25j
Z2 = Z1 + c = 0.25+0.375j
Z3 = Z2 + c = 0.171875+0.4375j
Z4 = Z3 + c = 0.088134+0.400390625j
Z5 = Z4 + c = 0.09745508+0.32057666j
Z6 = Z5 + c = 0.15672809+0.31248365j
Z7 = Z6 + c = 0.17691766+0.34794993j
Z8 = Z7 + c = 0.160230+0.37311697j
Z9 = Z8 + c = 0.1364575+0.36956959j
Как можно видеть, ряд не расходится, т.е. точка (0.25, 0.25) попадает в множество.
Разумеется, вручную делать такие проверки было бы крайне неудобно, собственно это одна из основных причин, по которой фракталы не были известны до появления компьютеров. Чтобы заметить какие-либо закономерности фрактальной природы, необходима возможность их автоматического рисования. Кстати, сам Мандельброт работал в IBM и имел доступ к весьма мощным компьютерам своего времени.
На языке Python проверку попадания точки в множество легко записать в виде функции. Для нас удобно то, что Python умеет работать с комплексными числами, что делает запись кода более короткой.
Сама функция имеет следующий вид:
def countMandelbrotIterationsForPt(pt):
z = 0
c = pt
threshold = 64
for iteration in xrange(threshold):
z = z*z + c
if abs(z) > 4:
break
return iteration
Здесь вычисляется 64 итерации, которых вполне достаточно чтобы определить, является ли ряд расходящимся или нет.
Чтобы сохранить весь фрактал, нужно перебрать все точки в диапазоне [-2, 2], программа, сохраняющая фрактал в файл, приведена ниже:
from PIL import Image
points = 1000
img = Image.new('RGB', (2*points,2*points), "black")
pixels = img.load()
# Range.X: -2..2,
# Range.Y: -2..2
for ix in xrange(-points, points, 1):
for iy in xrange(-points, points, 1):
pt = complex(2.0*ix/points, 2.0*iy/points)
i = countMandelbrotIterationsForPt(pt)
if i > 16:
colR, colG, colB = 4*i, 4*i, 0
if colR >= 255: colR = 255
if colB >= 255: colB = 255
if colB >= 255: colB = 255
img_x = points + ix
img_y = points + iy
pixels[img_x, img_y] = (colR, colG, colB)
if ix % 10 is 0:
print "Done: {}%".format(100.0*(ix + points)/(2*points))
img.save("fractal.png")
Как можно видеть, программа содержит 2 вложенных цикла ix, iy, затем значение преобразуются в комплексное число, для которого и выполняется описанная выше функция проверки.
Результат выполнения программы для 1000 точек, показан на рисунке:
Если в несколько раз увеличить количество точек, можно увидеть структуру верхней части более детально:
Немного изменив формулу, можно получить другие виды фракталов. Например, множество Жюлиа, описывается также, но Z0 = pt, и c = const. Казалось бы, небольшое отличие, приводит к совершенно другому рисунку фрактала:
Кстати, этот рисунок напоминает сорт цветной капусты “Романеско”:
Наконец, если в формуле фрактала Мандельброта Zn+1 = Zn2 + c, вычислять абсолютное значение Z*Z в виде (|a| + j|b|)*(|a| + j|b|), мы получим фрактал с названием “горящий корабль”:
Сложно сказать, насколько форма напоминает корабль, верхняя часть под увеличением скорее напоминает береговую линию, видимую с высоты:
Кстати, фрактальные алгоритмы - неотъемлемая часть современных методов рисования ландшафтов, таких как горы, деревья, облака и пр, эти объекты действительно имеют фрактальную природу.
В youtube можно найти интересные трехмерные визуализации фракталов, например, набрав в поиске “mandelbrot 3d”:
Для желающих поэкспериментировать самостоятельно, исходные коды программы генерации фракталов, приведены в приложении к книге. В приложении также есть готовые изображения фракталов, сгенерированные с высоким разрешением более 10000х10000 пикселов.
23. Горн Гавриила
В математике есть некоторые парадоксы, связанные с бесконечностью.
Так называемый “Горн Гавриила” (Gabriel's Horn) - это фигура, образованная функцией y = 1/x, если создать из нее 3х-мерную фигуру:
Интерес данного объекта в том, что “горн” имеет конечный объем - но при этом бесконечную площадь поверхности. Если бы такая фигура существовала в реальности, ее можно было бы заполнить краской, а вот покрасить целиком не получилось бы.
Действительно, формулы площади (V) и объема (A) “горна” легко найти в Википедии:
Из них видно, что при a = ∞, V становится равным ℼ, а A действительно стремится к бесконечности.
24. Построение графиков функций
Очень многие закономерности гораздо легче исследовать графически. Очень просто и эффективно строить графики с помощью языка Python и библиотеки matplotlib.
Выведем графики функций y = x2 и у = 0.5x2.
import matplotlib.pyplot as plt
plt.title('Graph')
plt.xlabel('X')
plt.ylabel('Y')
x_values = range(0,10)
y_values1 = map(lambda x: 0.5*x*x, x_values)
y_values2 = map(lambda x: x*x, x_values)
plt.axis([0, 10, 0, 20]) # xmin, xmax, ymin, ymax
plt.plot(x_values, y_values1, 'ro', label='0.5*x*x')
plt.plot(x_values, y_values2, marker='o', linestyle='--', label='x*x')
plt.legend()
plt.show()
Результат работы программы выглядит так:
Как можно видеть, мы имеем массив входных данных x_values = range(0,20), и 2 массива выходных данных y_values1 и y_values2, вычисляемых по соответствующим формулам. Все остальное за нас делает библиотека matplotlib. Команда plt.plot выводит график на экран, причем можно задать разные варианты оформления. Команда выводит подсказку (“легенду”), и является опциональной. Команды plt.title, plt.xlabel и plt.ylabel задают названия графика и осей. Команда plt.show выводит окно с графиком на экран.
Рассмотрим более сложный пример.
Выведем график так называемой “функции распределения простых чисел”, показывающий количество простых чисел меньше либо равных числу x.
Программа вывода графика выглядит так:
import matplotlib.pyplot as plt
import math
def is_prime(n):
if n % 2 == 0 and n > 2:
return False
for i in xrange(3, int(math.sqrt(n)) + 1, 2):
if n % i == 0:
return False
return True
def get_primes(n):
cnt = 0
for i in xrange(1, n):
if is_prime(i):
cnt += 1
return cnt
N = 400
x_values = range(2, N)
y_values = map(lambda n: get_primes(n), x_values)
plt.axis([0, N, 0, N/4])
plt.plot(x_values, y_values, 'ro', markersize=1)
plt.show()
Кстати, согласно Википедии, еще в 18м веке Гауссом и Лежандром было высказано предположение что функция распределения простых чисел выглядит так:
Это легко проверить графически, добавив строчку вывода второго графика:
plt.plot(x_values,map(lambda n:n/math.log(n),x_values), marker='o', linestyle='--')
Окончательный график выглядит так:
С помощью библиотеки numpy строить графики еще проще.
График функции sin(t) в диапазоне 0..2𝜋:
import matplotlib.pyplot as plt
import numpy as np
t = np.arange(0.0, 2.0, 0.01)
s = np.sin(2*np.pi*t)
plt.plot(t, s)
plt.show()
Важно заметить, что t здесь - не одно число, а целый массив - в данном примере используется возможность numpy работать непосредственно с массивами.
И наконец, с помощью matplotlib можно строить даже 3х-мерные графики функции двух переменных, для этого используется модуль mplot3d:
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import math
fig = plt.figure()
ax = fig.gca(projection='3d')
# Axis
ax.set_xlabel('X')
ax.set_xlim(0, 20)
ax.set_ylabel('Y')
ax.set_ylim(0, 20)
ax.set_zlabel('Z')
ax.set_zlim(0, 20)
# Graph
for x in range(0,20):
for y in range(0,20):
z = 10*math.sin(0.1*x + 0.1*y)
ax.scatter(x, y, z, c='r', s=1)
plt.show()
Результат показан на рисунке:
Интересно отметить кросплатформенность кода на Python - приведенная выше программа работает и на Windows и на OSX без каких-либо изменений.
25. Точность компьютерных вычислений
Запустим интерпретатор Python и введем простую программу:
a1 = 0.5 - 0.4
a2 = 1.5 - 1.4
print a1
print a2
Получим результаты 0.1 и еще раз 0.1. Пока все логично.
Введем новую команду:
print a1 == a2
В результате получаем … False. 0.1 != 0.1? Как такое может быть?
Чтобы убедиться, что значения не равны, введем:
print a1 - a2
Получаем вовсе не 0, как хотелось бы, а -1.1102230-16. Почему так? Разберемся.
Полученный результат - вовсе не ошибка, а вполне известная особенность компьютерных вычислений с вещественными числами. Любое число в языке программирования имеет свой определенный тип. Обычно это или целое число (1,2,3…), или так называемое число с плавающей точкой (1.1, 2.5, 3.8,...).
С целыми числами все просто - в памяти компьютера они хранятся в двоичном виде как суммы степеней 2. Например для числа “5”:
5 = 0*128 + 0*64 + 0*32 + 0*16 + 0*8 + 1*4 + 0*2 + 1*1 = 00000101b.
Любое десятичное число можно перевести в двоичное “туда” и “обратно”, соответствие тут однозначно, никакой потери точности нет. Ограничения лишь в количестве используемых байт, например для 1 байта (8 бит) максимальное число будет 128+64+32+16+8+4+2+1 = 255. Для 2х бит максимально хранимое значение равно 65535, для 4х бит 2147483647, и так далее.
C вещественными числами все сложнее. Для их хранения используется так называемый стандарт IEEE-754. Согласно этому стандарту, значения хранятся в так называемом формате “чисел с плавающей точкой”, состоящем из трех величин - знака (+ или -), мантиссы и экспоненты.
value = (1 + b0/2 + b1/4 + b2/8 + b3/16 + … ) * 2e-127
(b - биты мантиссы, e - экспонента)
К примеру, значение 3.14 в двоичном виде будет храниться так: 01000000010010001111010111000011. Первый 0 - это знак (+), 1000000 = 128 - это экспонента, а 10010001111010111000011 - это мантисса.
Достаточно проверить первые несколько чисел:
(1 + 1/2 + 1/16 + ...) * 2128-127 = 3.14
Ключевая особенность таких чисел - это их ограниченная точность, ведь на мантиссу и экспоненту отводится вполне определенное число бит. Ошибка весьма невелика, но тем не менее она есть. Например, число 0.5 будет храниться как 00111111000000000000000000000000, что действительно дает 0.5. А вот 0.6 будет храниться как 00111111000110011001100110011010, что в результате дает 0.60000002384185791015625.
Разумеется, записывая в языке программирования строку a=0.6, пользователь не задумывается о битах. Но как в поговорке про суслика, мы его не видим, а он есть.