Сложно сказать, насколько форма напоминает корабль, верхняя часть под увеличением скорее напоминает береговую линию, видимую с высоты:
Кстати, фрактальные алгоритмы - неотъемлемая часть современных методов рисования ландшафтов, таких как горы, деревья, облака и пр, эти объекты действительно имеют фрактальную природу.
В youtube можно найти интересные трехмерные визуализации фракталов, например, набрав в поиске “mandelbrot 3d”:
Для желающих поэкспериментировать самостоятельно, исходные коды программы генерации фракталов, приведены в приложении к книге. В приложении также есть готовые изображения фракталов, сгенерированные с высоким разрешением более 10000х10000 пикселов.
22. Горн Гавриила
В математике есть некоторые парадоксы, связанные с бесконечностью.
Так называемый “Горн Гавриила” (Gabriel's Horn) - это фигура, образованная функцией y = 1/x, если создать из нее 3х-мерную фигуру:
Интерес данного объекта в том, что “горн” имеет конечный объем - но при этом бесконечную площадь поверхности. Если бы такая фигура существовала в реальности, ее можно было бы заполнить краской, а вот покрасить целиком не получилось бы.
Действительно, формулы площади (V) и объема (A) “горна” легко найти в Википедии:
Из них видно, что при a = ∞, V становится равным ℼ, а A действительно стремится к бесконечности.
23. Построение графиков функций
Очень многие закономерности гораздо легче исследовать графически. Очень просто и эффективно строить графики с помощью языка 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 без каких-либо изменений.
24. Точность компьютерных вычислений
Запустим интерпретатор 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, пользователь не задумывается о битах. Но как в поговорке про суслика, мы его не видим, а он есть.
Введем команды в языке Python:
a=0.8
print a
Получим 0.8 - команда print достаточно “умна” чтобы округлить результат. Но сделаем вывод с большим числом знаком после запятой:
print("{0:.20f}".format(a))
Получим 0.8000000000000004441, что как говорится, и требовалось доказать.
Что с этим делать? Ну в общем-то ничего. Реальная ошибка в точности достаточно мала, чтобы это было проблемой. Но при желании можно использовать типы данных большей точности, правда ценой увеличения объема памяти и замедления скорости вычислений. Например, в С++ можно выбрать между типами float и double, второй имеет большую точность, но занимает больше памяти.
В целом, выполняя вычисления на компьютере, следует помнить что:
1) Если есть возможность использовать целочисленные вычисления - лучше использовать их, это помимо точности, эффективнее как по быстродействию, так и по занимаемой памяти.
2) Вещественные числа нельзя сравнивать оператором ==, как было показано выше, результат будет неверен. К примеру, следующий код не будет работать:
a = 0.5
b = 0.4
if a-b == 0.1:
print "0.1!"
Вместо этого достаточно убедиться, что результат мал:
if abs(a-b-0.1) < 1e-8:
print "0.1!"
3) Для денежных расчетов, где ошибки округления весьма нежелательны и баланс должен сойтись, существуют специальные библиотеки и типы данных, например библиотека python-money для Python.
4) Если все-таки желательно повысить точность, в Python можно использовать тип Decimal. Он хранит, по возможности, все значимые биты числа.
from decimal import Decimal
a = Decimal("0.5")
b = Decimal("0.4")
print a-b == Decimal("0.1")
Данная программа действительно вернет True. Но достигается это за счет значительного (до 10 раз или даже выше) замедления скорости расчета, да и не для всех вычислений это работает, например для уже упомянутого числа 1/3 точное значение сохранить невозможно. Но для таких случаев можно хранить числа в виде дробей, сохраняя и числитель и знаменатель. В общем, единого решения тут нет, каждый случай нужно рассматривать отдельно.
Но опять же, стоит повторить, что в большинстве случаев можно не беспокоиться, точности стандартных типов данных вполне достаточно. К примеру, если сохранить в переменной расстояние от Москвы до Владивостока, равное 9136.9км, то ошибка округления составит 10-11м, что лишь немногим больше размера одного атома. Скорее всего, такой ошибкой можно пренебречь.
25. Математика и календарь
Еще одна большая тема, где соединяются математика, астрономия и реальная жизнь - это наш обычный календарь. Например, часто можно услышать что високосный год считается “тяжелым”. Может ли такое быть?
Для начала, вспомним что Земля вращается по эллиптической орбите вокруг Солнца (что обеспечивает смену времен года), и в то же время, вокруг своей оси (что обеспечивает чередование дней). Второй период всем известен, и составляет 24 часа, а вот первый период известен уже не всем, и составляет 365.242187 суток. А в году, как известно, 365 дней.
На картинке это выглядит примерно так:
Первую неточность в календаре знали уже в древнем Риме, еще во времена Юлия Цезаря. Каждые 4 года добавлялся один день, чтобы скомпенсировать те самые недостающие 1/4 суток. В какой именно год добавлять - несложно догадаться, что принципиальной разницы тут нет, главное чтобы число лет было кратно 4. Можно даже добавлять по 2 дня каждые 8 лет, принципиально это ни на что не повлияет. Может ли високосный год быть “тяжелее” другого, соседнего? Очевидно что нет, т.к. все равно, какой именно год считать високосным - это лишь общепринятая условность, и на орбиту Земли выбор года никак не влияет.
Но если с орбитой все более-менее понятно, то с календарем не все оказалось так просто. Внимательный читатель наверное заметил, что разница в 0.242187 суток вовсе не кратна ¼. И совершенно справедливо. В юлианском календаре год оказывался длиннее астрономического на 11.2 минуты. Казалось бы немного, но за 128 лет накапливался целый день разницы. Впрочем, в древние времена ритм жизни был не столь быстр как сейчас, и ошибку заметили … только спустя 1.5 тысячи лет, когда “набежало” уже более 10 дней разницы. Это стало заметно не только астрономам - такой сдвиг календаря влиял и на даты посевов, и на празднование дней равноденствия. Тогда, в 1582 году папа Григорий XIII ввел новый календарь, названный теперь “григорианским” - в нем убрано часть високосных годов. Современное правило звучит так: год високосный, если он делится на четыре без остатка, но если он делится на 100 без остатка, это не високосный год, но если он делится без остатка на 400, это високосный год. Таким образом, 1996, 2000, 2004 - високосные годы, но 1900й год високосным не был. Правило выглядит запутанным, но зная, что день разницы накапливается за 128 лет, нетрудно понять почему именно так. Каждый 100й год 1 день “убирается” чтобы скомпенсировать разницу, набегающую за 100 лет, но каждый 400й год делается еще одно исключение, чтобы скомпенсировать оставшиеся 28 лет.