Реализация в Matlab алгоритмов построения фрактальных объектов

Задание

Тема курсового проекта: Реализация в MATLAB алгоритмов построения фрактальных объектов.

Исходные данные: Журнал Exponenta Pro. Математика в приложениях. № 3, 2003 г.

Основные вопросы, подлежащие исследованию:

1. Применение рекурсивного алгоритма.

2. Применение L-системы и терл — графики.

3. Применение системы итерированных функций.

К защите представить:

· пояснительную записку объемом не менее 40 листов:

· теоретическое обоснование;

· листинг программы;

· результаты расчетов.

· электронную версию всех материалов курсового проекта;

· слайды плакатов в Ms Off PowerPoint 2003:

· тема и цель курсового проекта;

· алгоритм решения задачи;

· результаты решения задачи в виде графиков и таблиц;

· экранный вид интерфейса задачи.

• Введение
• 1. Рекурсивный алгоритм
• 1.1 Ковёр Серпинского
• 1.1.1 Код программы «Serpinsky.m»
• 1.1.2 Изображения ковра Серпинского
• 1.2 Квадратный ковёр Серпинского
• 1.2.1 Код программы «Serpinsky2.m»
• 1.2.2 Изображения ковра Серпинского
• 1.3 Кривая Коха
• 1.3.1 Код программы «Koch.m»
• 1.3.2 Изображения кривой Коха
• 2. L — системы и терл — графика
• 2.1 Снежинка Коха
• 2.1.1 Код программы «RuleKoch.m» (возвращает функцию)
• 2.1.2 Код программы «Snowflake.m» (вывод изображения)
• 2.1.3 Изображение снежинки Коха
• 2.2 Дракон Хартера-Хайтвея
• 2.2.1 Код программы «Dracon.m»
• 2.2.2 Изображение дракона Хартера-Хайтвея
• 2.2.3 Изображение кривой Гильберта
• 2.2.4 Изображение кривой Госпера
• 2.2.5 Изображение кривой Серпинского
• 2.3 Ветвление
• 2.3.1 Код программы «Flower.m»
• 2.3.2 Изображение цветка
• 2.3.3 Изображение куста
• 2.3.4 Изображение Снежинки
• 3. Системы итерированных функций
• 3.1 Построение ковра Серпинского с помощью ДСИФ
• 3.1.1 Код программы «SerpDSIF.m»
• 3.1.2 Изображение ковра Серпинского построенного при помощи ДСИФ
• 3.2 Кристалл построенный по алгоритму РСИФ
• 3.2.1 Код программы «Cristal.m»
• 3.2.2 Изображение кристалла
• 3.2.3 Код программы «Maple.m»
• 3.2.4 Изображение Листа
• 3.2.5 Код программы «Paporotnic.m»
• 3.2.6 Изображение папоротника
• Заключение
• Список использованных источников
• Введение

Фракталы — математические объекты дробной размерности, название которых было введено в математику Б. Мандельбротом, являются в настоящее время, как предметом самостоятельных математических исследований, так и инструментарием, используемым в целом ряде прикладных задач нелинейной динамики, теории хаоса, обработки сигналов. Однако только относительно недавно появилось первое полноценное учебное пособие по новой быстро развивающейся математической дисциплине, основой которого стал учебный курс, преподававшийся автором книги в течение ряда лет в университете Миссури Колумбия. Так как при изучении фракталов и хаоса большую роль играет компьютерное моделирование, в курсе предусмотрено параллельное изучение теоретических вопросов и проведение компьютерных экспериментов. Это отличает его структуру от традиционной структуры большинства математических курсов: теорема-доказательство-пример-задача. В обобщенном виде подробно описаны известные алгоритмы построения фрактальных объектов (L-системы и терл-графика, аффинные преобразования, системы итерированных функций, случайные системы итерированных функций). Однако соответствующих программ, созданных на каком-либо языке программирования или в математическом пакете, не приводится. В то же время опыт практической реализации алгоритмов построения фрактальных объектов, описанных в каком-либо из современных математических пакетов (MATLAB, Mathcad, Maple, Matematica и т. д.), широко используемых в настоящее время в преподавании целого ряда физико-математических дисциплин, показывает, что существует необходимость внесения в них определенных корректировок, учитывающих особенности выбранного пакета (в первую очередь графические). В курсовой работе отражены алгоритмы построения классических фракталов и их программные реализации в пакете MATLAB.

1. Рекурсивный алгоритм

1.1 Ковёр Серпинского Рассматривается алгоритм, основанный на использовании рекурсивной функции. Приведён пример построения простого самоподобного фрактала — ковра Серпинского. В рассматриваемом в данном разделе алгоритме используется способ построения, основанный на последовательном удалении из начальной области внутренних подобластей в соответствии с заданными правилами.

Выберем в качестве начального множества S₀ равносторонний треугольник вместе с областью, которую он замыкает. Удалим внутренность центральной треугольной области и назовем оставшееся множество S₁. Затем повторим описанный процесс для каждого из трех оставшихся треугольников и получим приближение S₂. Продолжая, таким образом, получим последовательность вложенных множеств, S_n пересечение которых и образует ковер Серпинского S (рис. 1). Из построения видно, что ковер является объединением N=3 непересекающихся уменьшенных в два раза копий (коэффициенты подобия по горизонтали и вертикали в данном случае оказываются одинаковыми, r = ½). Фрактальная размерность ковра Серпинского d равна:

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

1. Задать порядок ковра N.

2. Задать координаты вершин исходного треугольника ABC: (X_A, Y_A), (X_B, Y_B), (X_C, Y_C).

3. Построить равносторонний треугольник ABC и залить его синим цветом.

Вычислить координаты середин сторон треугольника ABC:

5. Построить треугольник A? B?C? и залить его красным цветом.

6. Повторить N раз действия, описанные в пп. 4, 5, для треугольников AA? C?, A? BB?, C? B?C, соответственно.

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

Для реализации алгоритма в пакете MATLAB создана специальная функция, возвращающая изображение ковра Серпинского. Для этого использовался встроенный текстовый редактор MATLAB. Код программы приведён ниже. Изображение ковра Серпинского на рисунках 1 — 6.

1.1.1 Код программы «Serpinsky.m»

% Листинг файла Serpinsky. m

function z = Serpinsky (Lmax)

% функциЯ, возвращающая изображение ковра Серпинского

% Lmax — порядок ковра

% задание координат вершин равнобедренного треугольника

x1=0; y1=0; x2=1; y2=0; x3=0.5; y3=sin (pi/3);

h=figure (1); % инициализациЯ графического окна

hold on; % включение режима рисования фигур в одном графическом окне

fill ([x1 x2 x3],[y1 y2 y3],'b');

% прорисовка равностороннего треугольника

set (gca,'xtick',[],'ytick',[]); % отключение режима оцифровки осей

set (gca,'XColor','w','YColor','w'); % установка цвета рисования осей

Simplex (x1,y1,x2,y2,x3,y3,0,Lmax);

% обращение к функции, прорисовывающей равносторонние треугольники

% белого цвета

hold off % отключение режима рисования фигур в одном графическом окне

function z=Simplex (x1,y1,x2,y2,x3,y3,n, Lmax)

% рекурсивная функция, прорисовывающая равносторонние треугольники

% белого цвета

if n

% задание координат вершин текущего равностороннего треугольника

dx=(x2-x1)/2; dy=(y3-y1)/2; x1n=x1+dx; y1n=y1; x2n=x1+dx+dx/2;

y2n=y1+dy;x3n=x1+dx/2; y3n=y1+dy;

fill ([x1n x2n x3n],[y1n y2n y3n],'r');

n=n+1;

% рекурсия

Simplex (x1,y1,x1n, y1n, x3n, y3n, n, Lmax); Simplex (x1n, y1n, x2, y2,x2n, y2n, n, Lmax);

Simplex (x3n, y3n, x2n, y2n, x3, y3,n, Lmax); end

1.1.2 Изображения ковра Серпинского Рисунок 1 — Ковёр Серпинского 0 порядка

Рисунок 2 — Ковёр Серпинского 1 порядка

Рисунок 3 — Ковёр Серпинского 2 порядка

Рисунок 4 — Ковёр Серпинского 3 порядка

Рисунок 5 — Ковёр Серпинского 4 порядка

Рисунок 6 — Ковёр Серпинского 5 порядка

1.2 Квадратный ковёр Серпинского Алгоритм визуализации фрактального объекта, основные этапы построения которого представлены на рисунках 7 — 12, реализуется следующей последовательностью действий.

1. Задать порядок ковра N.

2. Задать координаты вершин исходного квадрата ABCD: (X_A, Y_A), (X_B, Y_B) (X_C, Y_C),(X_D, Y_D)_.

3. Построить квадрат ABCD и залить его синим цветом.

4. Вычислить координаты точек, делящих стороны квадрата ABCD на три равные части: dx=(X_B-X_A)/3, dy=(Y_B-Y_A)/3, X_A=X_A+dx, Y_A=Y_A+dy, X_B=X_A+dx+dx, Y_B=Y_A+dy, X_C=X_A+dx+dx, Y_C=Y_A+dy+dy, X_D=X_A+dx, Y_D=Y_A+dy+dy.

5. Построить квадрат A? B?C?D? и залить его белым цветом.

6. Повторить N раз действия, описанные в пп. 4, 5, для квадратов с вершинами, имеющими следующие координаты:

(X_A, Y_A), (X_A', Y_A), (X_A', Y_A'), (X_A, Y_A'),

(X_A', Y_A), (X_B', Y_A), (X_B', Y_B'), (X_A', Y_B'),

(X_B', Y_A), (X_B, Y_B), (X_B, Y_B'), (X_B', Y_B'),

(X_B', Y_B'), (X_B, Y_B'), (X_B, Y_C'), (X_C', Y_C'),

(X_C', Y_C'), (X_B, Y_C'), (X_C, Y_C), (X_C', Y_C'),

(X_D', Y_D'), (X_C', Y_C'), (X_C', Y_C'), (X_D', Y_C'),

(X_A, Y_D'), (X_D', Y_D'), (X_D', Y_C), (X_D, Y_D),

(X_A, Y_A'), (X_A', Y_D'), (X_D', Y_D'), (X_A, Y_D')

соответственно. Ниже приведен листинг файла Serpinsky2. m, который содержит описание функции, возвращающей изображение квадратного ковра Серпинского.

1.2.1 Код программы «Serpinsky2.m»

% Листинг файла Serpinsky2. m

function z=Serpinsky2(Lmax)

% функция, возвращающая изображение квадратного ковра Серпинского

% Lmax — порядок ковра

% Задание координат вершин исходного квадрата

x1=0; y1=0; x2=1; y2=0; x3=1; y3=1; x4=0; y4=1;

figure (1); hold on; fill ([x1 x2 x3 x4],[y1 y2 y3 y4],'b');

set (gca,'xtick',[],'ytick',[]); set (gca,'XColor','w','YColor','w');

Quadrate (x1,y1,x2,y2,x3,y3,x4,y4,0,Lmax);

hold off

function z=Quadrate (x1,y1,x2,y2,x3,y3,x4,y4,n, Lmax)

% рекурсивная функция, прорисовывающая квадраты белого цвета

if n

dx=(x2-x1)/3; dy=(y3-y1)/3;

x1n=x1+dx; y1n=y1+dy; x2n=x1+dx+dx; y2n=y1+dy;

x3n=x1+dx+dx; y3n=y1+dy+dy; x4n=x1+dx; y4n=y1+dy+dy;

fill ([x1n x2n x3n x4n],[y1n y2n y3n y4n],'r');

n=n+1;

Quadrate (x1,y1,x1+dx, y1, x1+dx, y1+dy, x1, y1+dy, n, Lmax);

Quadrate (x1+dx, y1, x1+2*dx, y1, x1+2*dx, y1+dy, x1+dx, y1+dy, n, Lmax);

Quadrate (x1+2*dx, y1, x2,y1,x2,y1+dy, x1+2*dx, y1+dy, n, Lmax);

Quadrate (x1+2*dx, y1+dy, x2, y1+dy, x2, y1+2*dy, x1+2*dx, y1+2*dy, n, Lmax);

Quadrate (x1+2*dx, y1+2*dy, x2, y1+2*dy, x2, y3,x1+2*dx, y3, n, Lmax);

Quadrate (x1+dx, y1+2*dy, x1+2*dx, y1+2*dy, x1+2*dx, y4, x1+dx, y4, n, Lmax);

Quadrate (x1,y1+2*dy, x1+dx, y1+2*dy, x1+dx, y4, x1,y4,n, Lmax);

Quadrate (x1,y1+dy, x1+dx, y1+dy, x1+dx, y1+2*dy, x1, y1+2*dy, n, Lmax);

End

1.2.2 Изображения ковра Серпинского Рисунок 7 — Квадратный ковёр Серпинского 0 порядка

Рисунок 8 — Квадратный ковёр Серпинского 1 порядка

Рисунок 9 — Квадратный ковёр Серпинского 2 порядка

Рисунок 10 — Квадратный ковёр Серпинского 3 порядка

Рисунок 11 — Квадратный ковёр Серпинского 4 порядка

Рисунок 12 — Квадратный ковёр Серпинского 5 порядка

1.3 Кривая Коха Еще одним примером фрактального объекта, для построения которого оказывается удобным использовать рекурсивный алгоритм, является кривая Коха. Построение данной кривой начинается с отрезка K₀ единичной длины.

Удалим из отрезка K₀ отрезок длины 1/3 и добавим два новых отрезка такой же длины. Назовем полученное множество K₁. На следующем шаге разделим каждый отрезок длины 1/3 на три части длины 1/9 и повторим описанную процедуру, заменяя на каждом шаге среднюю ветвь двумя новыми отрезками. Обозначим через K_n фигуру, получившуюся после n-го шага. Можно строго доказать, что последовательность кривых сходится к предельной кривой K бесконечной длины, фрактальная размерность которой равна Ниже приводится листинг рекурсивной функции, возвращающей изображение кривой Коха.

Изображение кривой Коха, возвращенное описанной выше функцией, представлено на рисунках 13−18.

1.3.1 Код программы «Koch.m»

% Листинг рекурсивной функции, возвращающей изображение кривой Коха

function z=Koch (N)

% функция, возвращающая изображение кривой Коха

x1=0; y1=0; % левая точка начального отрезка

x2=1; y2=0; % праваяя точка начального отрезка

figure (1); axis ([0 1 0 1]); set (gca,'XColor','w','YColor','w');

hold on;

Coord (x1,y1,x2,y2,N);

% вызов рекурсивной функции, прорисовывающей кривую Коха

function z=Coord (x1,y1,x2,y2,n)

% рекурсивная функция, прорисовывающая кривую Коха

if n>0

% вычисление координат концов отрезков на очередном шаге рекурсии

dx=(x2-x1)/3; dy=(y2-y1)/3;

x1n=x1+dx; y1n=y1+dy;

x2n=x1+2*dx; y2n=y1+2*dy;

xmid=dx/2-dy*sin (pi/3)+x1n; ymid=dy/2+dx*sin (pi/3)+y1n;

% рекурсиЯ

Coord (x1,y1,x1n, y1n, n-1);

Coord (x1n, y1n, xmid, ymid, n-1);

Coord (xmid, ymid, x2n, y2n, n-1);

Coord (x2n, y2n, x2, y2,n-1);

else

r1=[x1 y1]; r2=[x2 y2]; R=cat (1,r1,r2);

plot (R (, 1), R (, 2),'Color','b'); % построение кривой Коха

end;

1.3.2 Изображения кривой Коха Рисунок 13 — Кривая Коха 0 порядка

Рисунок 14 — Кривая Коха 1 порядка

Рисунок 15 — Кривая Коха 2 порядка

Рисунок 16 — Кривая Коха 3 порядка

Рисунок 17 — Кривая Коха 4 порядка

Рисунок 18 — Кривая Коха 5 порядка

2. L — системы и терл — графика

2.1 Снежинка Коха рекурсивный самоподобный фрактал программа Понятие «L-система» было введено А. Лидермайером в 1968 г. при изучении формальных языков. С их помощью оказывается возможным не только строить многие известные самоподобные фракталы, например, снежинку Коха, ковер Серпинского, кривые Пеано, Гильберта, Серпинского и др., но и создавать бесконечное разнообразие новых фракталов, укладывающихся в данную схему.

Терл-графика (от turtle — черепашка) является подсистемой вывода графического представления фрактального объекта. Основным исполнителем данной системы является «черепашка» (точка), которая перемещается по экрану дискретными шагами, прочерчивая или не прочерчивая свой след. «Мгновенное» положение «черепашки» задается тремя параметрами (x, y, б), где (x, y) — координаты «черепашки», б — направление следующего шага (угол, отсчитываемый от положительного направления оси x). Последовательность команд, определяющая направление перемещения и действия «черепашки», задается кодовым словом, буквы которого читаются слева направо. Кодовое слово, представляющее собой результат работы L-системы, может включать в себя следующие буквы:

· F — переместиться на один шаг вперед, прорисовывая след;

· b — переместиться на один шаг вперед, не прорисовывая след;

· [ — открыть ветвь;

· ] — закрыть ветвь;

· + — увеличить угол б на величину И;

·? — уменьшить угол б на величину И;

· X, Y — вспомогательные переменные.

Размер шага и величина приращения по углу И задаются изначально и остаются неизменными для всех перемещений «черепашки» .

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

Аксиома: F++F++F

Порождающее правило:

Newf =F — F + +F — F

И:= р/3

Очевидно, что графическим представлением аксиомы F + +F + +F является равносторонний треугольник. «Черепашка» делает один шаг вперед, затем угол увеличивается на 2р/3 и «черепашка» делает еще один шаг вперед, далее угол вновь увеличивается на 2р/3 и «черепашка» делает еще один шаг. На первом шаге каждая буква F в слове инициаторе заменяется на слово New:

F?F ++F ?F ++F? ++F?F ++F ?F ++F ?F

Повторяя этот процесс, на втором шаге получим:

F?F++F?F?F?F++F?F+F?F++F?F?F?F++F?F++F?F++F?F?F?F++F?F++F?F++F?F?F?F++F?F++F?F++F?F?F?F++F?F и т. д.

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

2.1.1 Код программы «RuleKoch.m» (возвращает функцию)

function z=RuleKoch (Lmax, Axiom, Newf, n, tmp);

% функция, возвращающая L-систему для снежинки Коха

% Входные параметры:

% Lmax — порядок снежинки

% Axiom — строка, содержащая аксиому

% Newf — строка, содержащая порождающее правило

% tmp — входная L-система

while n<=Lmax

if n==1

tmp=Axiom;

n=n+1;

else

tmp=strrep (tmp,'F', Newf); % замена всех букв F на

% порождающее правило

n=n+l;

tmp=RuleKoch (Lmax, Axiom, Newf, n, tmp); %

% рекурсия

end;

end;

z=tmp;

function [X, Y]=Koch (Lmax)

% функция, возвращающая изображение снежинки Коха

% Lmax — порядок снежинки

% порождающие правила

Axiom='F++F++F';

Newf='F-F++F-F';

teta=pi/3;

alpha=0;

p=[0;0]; % стартовая точка

p=Coord (p, Lmax, Axiom, Newf, alpha, teta, p); % обращение к функции,

% возвращающей координаты вершин

M=size (p, 2); % число вершин снежинки Коха

X=p (1:1,1:M); % создание вектора, содержащего Х-е координаты вершин снежинки

Y=p (2:2,1:M); % создание вектора, содержащего Y-e координаты вершин снежинки

figure (1); % инициализация графического окна

plot (X, Y,'Color','k'); % построение снежинки Коха

set (gca,'xtick',[],'ytick',[]);

set (gca,'XColor','w','YColor','w');

function z=Coord (p, Lmax, Axiom, Newf, alpha, teta)

% функция, возвращающая координаты вершин снежинки Коха

Rule=RuleKoch (Lmax, Axiom, Newf, 1,' '); % задание L-системы

M=length (Rule);

for i=l:M

Tmp=p (1:2, size (p, 2):size (p, 2));

if Rule (i)=='F' % шаг вперед

R=[cos (alpha);sin (alpha)];

R=R/(4^Lmax);

Tmp=Tmp+R;

p=cat (2,p, Tmp);

end;

if Rule (i)=='+' % увеличение угла, задающего направление движения

alpha=alpha+teta;

end;

if Rule (i)== '-' % уменьшение угла, задающего направление движения

alpha=alpha-teta; end; end;

z=p;

2.1.2 Код программы «Snowflake.m» (вывод изображения)

function [X, Y]=Snowflake (Lmax)

% функция, возвращающая изображение

% снежинки Коха

% Lmax — порядок снежинки

% порождающие правила

Axiom='F++F++F';

Newf='F-F++F-F';

teta=pi/3;

alpha=0;

p=[0;0]; % стартовая точка

p=Coord (p, Lmax, Axiom, Newf, alpha, teta);

% обращение к функции, возвращающей координаты вершин

M=size (p, 2); % число вершин снежинки Коха

X=p (1:1,1:M);

% создание вектора, содержащего Х-е координаты вершин снежинки

Y=p (2:2,1:M);

% создание вектора, содержащего Y-е координаты вершин снежинки

figure (1); % инициализация графического окна

plot (X, Y, 'Color', 'k'); % построение снежинки Коха

function z=Coord (p, Lmax, Axiom, Newf, alpha, teta)

% функциЯ, возвращающаЯ координаты вершин снежинки Коха

Rule=RuleKoch (Lmax, Axiom, Newf, 1, ' '); % задание L-системы

M=length (Rule);

for i=1:M

Tmp=p (1:2,size (p, 2):size (p, 2));

if Rule (i)== 'F ' % шаг вперед

R=[cos (alpha);sin (alpha)]; R=R/(4^Lmax);

Tmp=Tmp+R; p=cat (2,p, Tmp);

end

if Rule (i)== '+' % увеличение угла, задающего направление движения

alpha=alpha+teta;

end

if Rule (i)== '-' % уменьшение угла, задающего направление движения

alpha=alpha-teta;

end

end

z=p;

2.1.3 Изображение снежинки Коха Рисунок 19 — Снежинка Коха 4 порядка

2.2 Дракон Хартера-Хайтвея Однако при построении других фракталов, например, дракона Хартера-Хайтвея (рисунок 20), на некоторых шагах возникает необходимость в изменении направления чтения правила (не слева направо, а справа налево). Для решения данной проблемы вводят две дополнительные команды, обозначаемые X и Y, которые используются для создания соответствующей L-системы, но игнорируются «черепашкой» при перемещении. При использовании этих команд порождающее правило для дракона имеет вид:

Аксиома: FX

Порождающие правила:

Newf = FX

Newx = X +YF +

Newy = ?FX ?Y

В соответствии с данными правилами L-система имеет следующий вид:

1 шаг: FX+YF+

2 шаг: FX+YF++?FX?YF+

3 шаг: FX+YF ++?FX?YF++?FX+YF+??FX?YF+

4 шаг:

FX+YF++?FX?YF++?FX+YF+??FX?YF++?FX+YF++?FX?YF+??FX+YF+??FX?YF+

Ниже приводится листинг файла Dracon. m, содержащего описание функции, возвращающей изображение дракона, в соответствии с описанной выше L-системой.

2.2.1 Код программы «Dracon.m»

function [X, Y]=Dracon (Lmax)

% функция, возвращающая изображение дракона

% Lmax — порЯядок дракона

% порождающие правила

Axiom='FX';

Newf='F';

Newx='X+YF+';

Newy='-FX-Y';

teta=pi/2; alpha=0; p=[0;0];

p=Coord (Lmax, Axiom, Newf, Newx, Newy, alpha, teta);

% обращение к функции, возвращающей координаты угловых точек дракона

M=size (p, 2);

X=p (1:1,1:M);

% создание вектора, содержащего Х-е координаты угловых точек дракона

Y=p (2:2,1:M);

% создание вектора, содержащего Y-е координаты угловых точек дракона

plot (X, Y,'Color','k');

set (gca,'xtick', 'уtick',[]);

set (gca,'XColor','w','YColor','w');

function z=Coord (p, Lmax, Axiom, Newf, Newx, Newy, alpha, teta)

% функция, возвращающая координаты угловых точек дракона

Rule=DraconString (Lmax, Axiom, Newf, Newx, Newy, 1, ' '); % задание L-системы

M=length (Rule)

for i=l:M

Tmp=p (1:2, size (p, 2) :size (p, 2));

if Rule (i)=='F' % шаг вперед

R=[cos (alpha);sin (alpha)];

R=R/(2^Lmax);

p=cat (2,p, Tmp);

end;

if Rule (i) == ' + ' % увеличение угла, задающего направление движения

alpha=alpha+teta;

end;

if Rule (i) == '-' % уменьшение угла, задающего направление движения

alpha=alpha-teta;

end;

end;

z=p;

function z=DraconString (Lmax, Axiom, Newf, Newx, Newy, n, tmp);

% функция, возвращающая L-систему

a=1;

if n<=Lmax

if n==1

tmp=Axiom;

end;

M=length (tmp);

tmpl=' ';

for i=l:M

if tmp (i)=='X' tmpl=streat (tmpl, Newx); end;

if tmp (i)=='Y' tmpl=streat (tmpl, Newy); end;

if not (tmp (i)=='F') ¬(tmp (i)=='X') ¬(tmp (i)=='Y') tmpl=strcat (tmpl, tmp (i)); end;

tmp=tmpl; n=n+l; tmp=DraconString (Lmax, Axiom, Newf, Newx, Newy, n, tmp);

end;

end;

z=tmp;

2.2.2 Изображение дракона Хартера-Хайтвея Рисунок 20 — Дракон Хартера — Хайтвея 12 порядка

2.2.3 Изображение кривой Гильберта Заменяя в описанной выше программе порождающие правила, можно получить и другие фрактальные кривые, например, кривую Гильберта. Аксиома: X

Порождающие правила:

Newf = F

Newx = -YF + XFX + FY ;

Newy = +XFYFY — FX +

б = 0, и = р/2

Рисунок 21 — Кривая Гильберта 4 порядка

2.2.4 Изображение кривой Госпера Аксиома: XF

Порождающие правила:

Newx = X +YF + +YF? FX? ?FXFX ?YF +

Newy = ?FX +YFYF + +YF + FX? ?FX ?Y д б = 0, и = р/3

Рисунок 22 — Кривая Госпера 3 порядка

2.2.5 Изображение кривой Серпинского Аксиома: F + XF + F + XF

Порождающие правила:

Newf = F

Newx = XF? F + F? XF + F + XF? F + F? X

Newy = ' '

б = р 4, и = р/2.

Рисунок 23 — Кривая Серпинского 3 порядка

2.3 Ветвление Когда в L-системе встречается символ [ (открыть ветвь), необходимо запомнить координаты точки нахождения «черепашки» и направление ее движения, т. е. переменные (x, y, б). К сохраненным переменным следует вернуться после обнаружения символа ] (закрыть ветвь). Для хранения триплетов (x, y, б) в предлагается использовать стек:

в конец которого записываются новые данные. При закрытии ветви переменным (x, y, б) присваиваются значения, считанные из конца стека, затем эти значения из стека удаляются. В пакете MATLAB оказывается более удобным использовать матрицу с переменным числом столбцов:

причем координаты каждой новой точки ветвления добавляются в новый столбец матрицы M. После закрытия ветви переменным (x, y, б) присваиваются значения, считанные из последнего столбца матрицы M, затем этот столбец удаляется.

Таким образом, ветвление задается двумя символами:

[ — открыть ветвь: добавить вектор

] — закрыть ветвь: присвоить переменным

(x, y, б) значения координат вектора, являющегося последним столбцом матрицы M.

Пример фрактала, построенного с помощью операции ветвления, представлен на рисунке 24. Ниже приводится листинг файла Flower. m, содержащего описание функции, возвращающей изображение цветка, в соответствии с описанной выше L-системой.

2.3.1 Код программы «Flower.m»

% Листинг файла Flower. m

function [X, Y]=Flower (Lmax)

% функция, возвращающая изображение цветка

% Lmax — порядок цветка

% порождающие правила

Axiom='F[+F+F][-F-F][++F][—F]F';

Newf='FF[++F][+F][F][-F][—F] ';

teta=pi/16;

alpha=pi/2;

p=[0;0]; % начальная точка

Coord (p, Lmax, Axiom, Newf, alpha, teta);

% обращение к функции, возвращающей

% изображение цветка

function z=Coord (p, Lmax, Axiom, Newf, alpha, teta)

% функциЯ, возвращающаЯ изображение цветка

Rule=FlowerString (Lmax, Axiom, Newf, 1, ' ');

% задание L — системы

figure (1); hold on;

M=length (Rule);L=0; x0=p (1);y0=p (2);

for i=1:M

if Rule (i)== 'F' % шаг вперед

x1=x0+cos (alpha); y1=y0+sin (alpha);

X=[x0,x1]; Y=[y0,y1]; x0=x1; y0=y1;

plot (X, Y, 'Color', 'k'); end

if Rule (i)== '+' % увеличение угла, задающего направление движениЯ

alpha=alpha+teta;

end

if Rule (i)== '! ' % уменьшение угла, задающего направление движениЯ

alpha=alpha-teta;

end

if Rule (i)== ' ['; % открыть ветвь

if L==0 St=[x0;y0;alpha]; L=1; else St=cat (2,St,[x0;y0;alpha]);

end end

if Rule (i)== ']' % закрыть ветвь

M=size (St, 2); R=St (1:3,M:M);

x0=R (1); y0=R (2);

alpha=R (3); St=St (1:3,1:M-1);

end

end

hold off

function z=FlowerString (Lmax, Axiom, Newf, n, tmp)

% функция, возвращающая L — систему

while n<=Lmax

if n==1

tmp=Axiom; n=n+1; else

tmp=strrep (tmp, 'F', Newf);

n=n+1; tmp=FlowerString (Lmax, Axiom, Newf, n, tmp); % рекурсия

end

end

z=tmp;

2.3.2 Изображение цветка Рисунок 24 — Цветок 3 порядка

2.3.3 Изображение куста Заменяя в описанной программе порождающие правила, можно получить и другие ветвящиеся фрактальные объекты, например, куст.

Аксиома: F

Порождающие правила:

Newf = ?F + F +[+F? F?]?[?F + F + F]

б = р/2, и = р/8

Рисунок 25 — Куст 3 порядка

2.3.4 Изображение Снежинки Аксиома: [F]+[F]+[F]+[F]+[F]+[F]

Порождающие правила:

Newf = F[+ + F][?FF]FF[+F][?F]FF

б = 0, и = р/3

Рисунок 26 — Снежинка 3 порядка

3. Системы итерированных функций

3.1 Построение ковра Серпинского с помощью ДСИФ В общем случае для построения системы итерированных функций (СИФ) в рассмотрение вводится совокупность сжимающих отображений действующих в R. Эти m отображений используются для построения одного сжимающего отображения T в пространстве Щ всех непустых компактов из R Преобразование Хатчинсона T: Щ>Щ определяется следующим образом:

Данное преобразование ставит в соответствие «точкам» из Щ, под которыми здесь понимаются компактные множества, также «точки» из Щ. Системой итерированных функций называется совокупность приведенных выше отображений вместе с итерационной схемой Например, СИФ при построении ковра Серпинского задается тремя аффинными преобразованиями, которые в матричной форме имеют следующий вид:

Таким образом, для построения ковра Серпинского в пакете MATLAB с помощью ДСИФ можно использовать следующий алгоритм.

1. Задать порядок ковра Серпинского n.

2. Задать число точек начальной конфигурации m.

3. Задать координаты i точек (i = 1, 2,…, m), заполняющих начальное множество.

4. Перевести каждое из чисел 1, 2,…, 3n в троичную систему счисления.

5. Сформировать массив, состоящий из 3n строк длиной n символов.

6. Задать аффинные преобразования.

7. Для i-ой точки начальной конфигурации последовательно применить каждое из j = 1, 2,…, 3n итерационных правил и отобразить в графическом окне полученные образы каждой начальной точки.

Пример ковра Серпинского, построенного с помощью описанной выше модификации алгоритма ДСИФ, представлен на рисунке 27. Ниже приводится листинг файла SerpDSIF. m, содержащего описание соответствующей функции, возвращающей изображение ковра Серпинского.

3.1.1 Код программы «SerpDSIF.m»

% Листинг файла SerpDSIF. m

function z=SerpDSIF (Niter, NPoints)

% функциЯ, возвращающая изображение ковра Серпинского

% Niter — порядок ковра

% NPoints — число точек начальной конфигурации

x=zeros (NPoints, 1); y=zeros (NPoints, 1);

% задание координат точек начальной конфигурации

x1=0; y1=0; x2=1; y2=0; x3=½; y3=sin (pi/3);

j=1;

while j<=NPoints

tmpx=rand (1,1);

tmpy=sqrt (3)/2*rand (1,1);

if (-sqrt (3)*tmpx+tmpy<=0)&(sqrt (3)*tmpx+tmpy-sqrt (3)<=0)

x (j)=tmpx;

y (j)=tmpy;

j=j+1;

end;

end

% Формирование массива, содержащего правила итерации

for i=1:3^Niter

Tmp (i)=system3(i-1);

% перевод числа из десятичной в троичную систему счисления

end

n=1; s='0';

while n

s=strcat (s, '0'); n=n+1;

end

for i=1:3^Niter

tmp=num2str (Tmp (i)); tmp1=s;

for m=1:length (tmp)

tmp1(Niter-m+1)=tmp (length (tmp)-m+1);

end

Cod (i, 1: Niter)=tmp1;

end

a1=[0;0]; a2=[½;0]; a3=[¼;sqrt (3)/4]; A=[½, 0;0,½];

% задание аффинных преобразований

figure (1); hold on; set (gca, 'xtick',[], 'ytick',[]);

set (gca, 'XColor', 'w', 'YColor', 'w'); fill ([x1 x2 x3],[y1 y2 y3], 'w');

GosperDraw (Niter, NPoints, x, y, A, a1,a2,a3,Cod);

% построение ковра Серпинского

function z=GosperDraw (Niter, NPoints, x, y, A, a1,a2,a3,Cod)

% функциЯ, создающаЯ изображение ковра Серпинского

for m=1:3^Niter

X=x; Y=y; Rule=Cod (m:);

for i=1:Niter

tmp=Rule (Niter+1!i);

if tmp=='0 '

[X Y]=T (NPoints, X, Y, A, a1); % первое аффинное преобразование

end

if tmp=='1'

[X Y]=T (NPoints, X, Y, A, a2); % второе аффинное преобразование

end

if tmp=='2'

[X Y]=T (NPoints, X, Y, A, a3); % третье аффинное преобразование

end

end

plot (X, Y, '. ', 'MarkerSize', 1, 'MarkerEdgeColor', 'b');

% отображение результатов итерации

end

function [X, Y]=T (NPoints, x, y, A, a)

% ФункциЯ, возвращающаЯ результат аффинного преобразования

X=zeros (NPoints, 1); Y=zeros (NPoints, 1);

for i=1:NPoints

R=[x (i);y (i)]; R=A*R+a; X (i)=R (1); Y (i)=R (2);

end

function z=system3(D);

% функциЯ, возвращающаЯ значение целого числа

% в троичной системе счислениЯ

% D — число в десЯтичной системе счислениЯ

n=1;

while D>=3^n n=n+1; end

if n>1

a=floor (D/3^(n-1))*10^(n-1); b=mod (D, 3^(n-1));

if b>=3

b=system3(b); % рекурсия

end

z=a+b;

else

z=D;

end

3.1.2 Изображение ковра Серпинского построенного при помощи ДСИФ Рисунок 27 — Ковёр Серпинского 4 порядка (ДСИФ)

3.2 Кристалл построенный по алгоритму РСИФ В отличие от ДСИФ в рандомизированном алгоритме начальное множество S₀ состоит из одной точки (x₀, y₀), а правило, по которому точке ставится в соответствие точка (x_i, y_i), где i — номер правила, выбирается случайным образом из набора, содержащего все возможные правила аффинных преобразований. Например, применительно к ковру Серпинского это означает, что при построении ковра 2_го порядка преобразование должно случайным образом выбираться из следующего множества преобразований:

число элементов N которого равно

Таким образом, для построения ковра Серпинского в пакете MATLAB с помощью РСИФ можно использовать следующий алгоритм.

1. Задать порядок ковра Серпинского n.

2. Задать число испытаний NTrial.

3. Задать число аффинных преобразований m = 3.

4. Сформировать массив, содержащий набор правил для аффинных преобразований.

5. Задать координаты начальной точки (x₀, y₀).

6. Перевести каждое из чисел 1, 2,…, 3n в троичную систему счисления.

7. Задать аффинные преобразования.

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

9. Отобразить вычисленное множество точек в графическом окне.

Пример кристалла, построенного с помощью описанного выше алгоритма РСИФ, представлен на рисунке 28. Ниже приводится листинг файла Cristal. m, содержащего описание соответствующей функции.

3.2.1 Код программы «Cristal.m»

% Листинг файла Cristal. m

function z=Cristal (Niter, NTrial)

% функция, возвращающая изображение кристалла

% Niter — порядок кристалла

% NTrial — число испытаний

Na=4; % число аффинных преобразований

% Создание массива, содержащего набор

% правил для аффинных преобразований

k=1;

for m=1:Niter

for i=1:4^m Tmp (k)=system3(i-1,Na); k=k+1;

end

Q (1)=Na;

for m=2:Niter Q (m)=Q (m-1)+Na^m; end

n=1; s='0'; M=1;

while n<=length (Tmp)

m=1;

while n>Q (m) m=m+1; end

if m==1

S (n, 1:1)=s;

else

S (n, 1:1)=s;

for i=2:m S (n, 1: i)=strcat (S (n:), s); end

end

n=n+1;

end

for i=1:k-1 tmp=num2str (Tmp (i)); m=1;

while i>Q (m) m=m+1; end

tmp1(1:m)=S (i, 1: m);

for m=1:length (tmp)

tmp1(length (tmp1)-m+1:length (tmp1)!m+1)=tmp (length (tmp)!m+1:length (tmp)!m+1);

end

Cod (i, 1: length (tmp1))=tmp1;

end

x=0; y=0; % координаты начальной точки

% задание аффинных преобразований

A1=[0.2550,0.0000;0.0000,0.2550];

A2=[0.2550,0.0000;0.0000,0.2550];

A3=[0.2550,0.0000;0.0000,0.2550];

4=[0.3700,-0.6420;0.6420,0.3700];

a1=[0.3726;0.6714]; a2=[0.1146;0.2232];

a3=[0.6306;0.2232]; a4=[0.6356;!0.0061];

figure (1); hold on;

set (gca, 'xtick',[], 'ytick',[]);

set (gca, 'XColor ', 'w', 'YColor', 'w');

DrawFractal (Niter, NTrial, x, y, A1, A2,A3,A4,…

a1,a2,a3,a4,Cod); % визуализация фрактала

function z=DrawFractal (Niter, NTrial, x, y,…

A1,A2,A3,A4,a1,a2,a3,a4,Cod)

% функция, возвращающая изображение фрактала

X1=zeros (NTrial, 1); Y1=zeros (NTrial, 1);

X=x; Y=y;

for m=1:NTrial

Np=1+round ((size (Cod, 1)-1)*rand (1,1));

% выбор номера преобразования

Rule=Cod (Np:);

for i=1:length (Rule)

tmp=Rule (length (Rule)+1-i);

if tmp=='0' [X Y]=T (X, Y, A1,a1); end

if tmp=='1' [X Y]=T (X, Y, A2,a2); end

if tmp=='2 ' [X Y]=T (X, Y, A3,a3); end

if tmp=='3' [X Y]=T (X, Y, A4,a4); end

end

X1(m)=X; Y1(m)=Y;

end

plot (X1,Y1, '. ', 'MarkerSize', 1,'MarkerEdgeColor', 'b');

function [X, Y]=T (x, y, A, a)

% Функция, возвращающая результат

% аффинного преобразования

R=[x;y]; R=A*R+a;

X=R (1); Y=R (2);

function z=system3(D, m);

% функция, возвращающая значение

% числа в четверичной системе координат

n=1;

while D>=m^n n=n+1; end

if n>1

a=floor (D/m^(n-1))*10^(n-1);

b=mod (D, m^(n-1));

if b>=m b=system3(b, m); end

z=a+b;

else

z=D;

end

3.2.2 Изображение кристалла Рисунок 28 — Кристалл (РСИФ)

3.2.3 Код программы «Maple.m»

function z=Maple (Niter, NPoints)

x1=0; y1=0;

x2=1; y2=0;

x3=½; y3=sin (pi/3);

x=zeros (NPoints, 1); y=zeros (NPoints, 1);

j=1;

while j<=NPoints

tmpx=rand (1,1);

tmpy=sqrt (3)/2*rand (1,1);

if (-sqrt (3)*tmpx+tmpy<=0)&(sqrt (3)*…

tmpx+tmpy-sqrt (3)<=0)

x (j)=tmpx; y (j)=tmpy; j=j+1;

end

end

for i=1:2^Niter Tmp (i)=system2(i-1,2); end

n=1; s='0';

while n

for i=1:2^Niter

tmp=num2str (Tmp (i)); tmp1=s;

for m=1:length (tmp) tmp1(Niter-m+1)=…

tmp (length (tmp)-m+1); end

Cod (i, 1: Niter)=tmp1;

end

A1=[0.4,-0.3733;0.0600,0.6000];

A2=[-0.8000,-0.1867;0.1371,0.8000];

a1=[0.3533;0.000]; a2=[1.1000;0.1000];

figure (1); hold on;

set (gca,'xtick',[],'ytick',[]);

set (gca,'XColor','w','YColor','w');

FractalDraw (Niter, NPoints, x, y, A1, A2,…

a1,a2,Cod);

function z=FractalDraw (Niter, NPoints, x, y,…

A1,A2,a1,a2,Cod)

for m=1:2^Niter

X=x; Y=y;

Rule=Cod (m:);

for i=1:Niter

tmp=Rule (Niter+1-i);

if tmp=='0' [X Y]=T (NPoints, X, Y, A1, a1); end

if tmp=='1' [X Y]=T (NPoints, X, Y, A2, a2); end

end

plot (X, Y,'.','MarkerSize', 1,…

'MarkerEdgeColor','b'); end

function [X, Y]=T (NPoints, x, y, A, a)

X=zeros (NPoints, 1); Y=zeros (NPoints, 1);

for i=1:NPoints R=[x (i);y (i)]; R=A*R+a;

X (i)=R (1); Y (i)=R (2); end

function z=system2(D, m);

n=1;

while D>=m^n n=n+1; end

if n>1

a=floor (D/m^(n-1))*10^(n-1);

b=mod (D, m^(n-1));

if b>=m b=system2(b, m);

end

z=a+b;

else

z=D; end

3.2.4 Изображение Листа Рисунок 29 — Лист (ДСИФ)

3.2.5 Код программы «Paporotnic.m»

function z=Paporotnic (Niter, NPoints)

x1=0; y1=0; x2=1; y2=0; x3=½; y3=sin (pi/3);

x=zeros (NPoints, 1); y=zeros (NPoints, 1);

j=1;

while j<=NPoints

tmpx=rand (1,1); tmpy=sqrt (3)/2*rand (1,1);

if (-sqrt (3)*tmpx+tmpy<=0)&(sqrt (3)*tmpx+tmpy-sqrt (3)<=0)

x (j)=tmpx; y (j)=tmpy; j=j+1; end; end;

for i=1:4^Niter Tmp (i)=system2(i-1,4); end;

n=1; s='0';

while n

for i=1:4^Niter

tmp=num2str (Tmp (i)); tmp1=s;

for m=1:length (tmp)

tmp1(Niter-m+1)=tmp (length (tmp)-m+1);

end;

Cod (i, 1: Niter)=tmp1; end;

A1=[0.7000,0;0,0.7000]; A2=[0.1000,-0.4330;0.1732,0.2500];

A3=[0.1000,0.4330;-0.1732,0.2500]; A4=[0,0;0,0.3000];

a1=[0.1496;0.2962]; a2=[0.4478;0.0014];

a3=[0.4450;0.1559]; a4=[0.4987;0.0070];

figure (1); hold on; set (gca,'xtick',[],'ytick',[]);

set (gca,'XColor','w','YColor','w');

FractalDraw (Niter, NPoints, x, y, A1, A2,A3,A4,a1,a2,a3,a4,Cod);

function z=Simplex (Niter, NPoints, x, y, A1, A2,A3,A4,a1,a2,a3,a4,Cod)

for m=1:4^Niter

X=x; Y=y;

Rule=Cod (m:);

for i=1:Niter

tmp=Rule (Niter+1-i);

if tmp=='0' [X Y]=T (NPoints, X, Y, A1, a1); end;

if tmp=='1' [X Y]=T (NPoints, X, Y, A2, a2); end;

if tmp=='2' [X Y]=T (NPoints, X, Y, A3, a3); end;

if tmp=='3' [X Y]=T (NPoints, X, Y, A4, a4); end;

end;

plot (X, Y,'.','MarkerSize', 1,'MarkerEdgeColor','k') end;

function [X, Y]=T (NPoints, x, y, A, a)

X=zeros (NPoints, 1); Y=zeros (NPoints, 1);

for i=1:NPoints R=[x (i);y (i)]; R=A*R+a; X (i)=R (1); Y (i)=R (2); end;

function z=FractalDraw (Niter, NPoints, x, y, A1, A2,a1,a2,Cod)

for m=1:2^ Niter

X=x; Y=y; Rule=Cod (m:);

for i=1:Niter tmp=Rule (Niter+1-i);

if tmp=='0' [X Y]=T (NPoints, X, Y, A1, a1); end;

if tmp=='1' [X Y]=T (NPoints, X, Y, A2, a2); end;

end;

plot (X, Y,'.','MarkerSize', 1,'MarkerEdgeColor','k'); end;

function [X, Y]=T (NPoints, x, y, A, a)

X=zeros (NPoints, 1); Y=zeros (NPoints, 1);

for i=1:NPoints R=[x (i);y (i)]; R=A*R+a; X (i)=R (1); Y (i)=R (2); end;

function z=system2(D, m);

n=1;

while D>=m^n n=n+1; end;

if n>1

a=floor (D/m^(n-1))*10^(n-1); b=mod (D, m^(n-1));

if b>=m b=system2(b, m); end;

z=a+b;

else z=D; end;

3.2.6 Изображение папоротника Рисунок 30 — Папоротник (ДСИФ)

Заключение

В результате выполнения курсовой работы рассмотрены фракталы.

Фракталы — математические объекты дробной размерности, название которых было введено в математику Б. Мандельбротом, являются в настоящее время, как предметом самостоятельных математических исследований, так и инструментарием, используемым в целом ряде прикладных задач нелинейной динамики, теории хаоса, обработки сигналов. Однако только относительно недавно появилось первое полноценное учебное пособие по новой быстро развивающейся математической дисциплине, основой которого стал учебный курс, преподававшийся автором книги в течение ряда лет в университете Миссури Колумбия. Так как при изучении фракталов и хаоса большую роль играет компьютерное моделирование, в курсовой работе рассмотрено параллельное изучение теоретических вопросов и проведение компьютерных экспериментов.

В обобщенном виде подробно описаны известные алгоритмы построения фрактальных объектов (L-системы и терл-графика, аффинные преобразования, системы итерированных функций, случайные системы итерированных функций и др.).

Разработаны соответствующие программы, созданные в математическом пакете Matlab.

Список использованных источников

1. Mandelbrot B. B. Les object fractals: forme, hasard et dimantion.— Paris: Flamarion, 1975.

2. Каток А. Б., Хасселблат Б.

Введение

в современную теорию динамических систем.— М.: Факториал, 1999.

3. Кренкель Э. Т. Сжатие сигналов с применением теории фракталов.— М.: МТУСИ, 1996.

4. Кроновер Р. М. Фракталы и хаос в динамических системах. Основы теории.— М.: Постмаркер, 2000.

5. Лихтенберг А., Либерман М. Регулярная и стохастическая динамика.— М.: Меркурий_ПРЕСС, 2000.

6. Шустер Г. Детерминированный хаос.— М.: Мир, 1988.

7. Журнал Exponenta Pro. Математика в приложениях. № 3, 2003 г.