六个分形图形的matlab实现(6)

几个分形的matlab 实现
摘要：给出几个分形的实例，并用matlab 编程实现方便更好的理解分形，欣赏其带来的
数学美感
关键字：Koch 曲线 实验 图像
一、问题描述：
从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的两边代替，形成山丘形图形如下
图1
在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的两条边代替，再次形成新的图形如此迭代，形成Koch 分形曲线。

二、算法分析：
考虑由直线段（2个点）产生第一个图形（5个点）的过程。

图1中，设1P 和5P 分别为原始直线段的两个端点，现需要在直线段的中间依次插入三个点2P ，3P ，4P 。

显然2P 位于线段三分之一处，4P 位于线段三分之二处，3P 点的位置可看成是由4P 点以2P 点为轴心，逆时针旋转600
而得。

旋转由正交矩阵 ⎪⎪⎪⎪⎭⎫ ⎝
⎛-=)3cos()3sin()3sin()3cos(ππππA 实现。

算法根据初始数据（1P 和5P 点的坐标），产生图1中5个结点的坐标。

结点的坐标数组形成一个25⨯矩阵，矩阵的第一行为1P 的坐标，第二行为2P 的坐标……，第五行为5P 的坐标。

矩阵的第一列元素分别为5个结点的x 坐标，第二列元素分别为5个结点的y 坐标。

进一步考虑Koch 曲线形成过程中结点数目的变化规律。

设第k 次迭代产生的结点数为k n ，第1+k 次迭代产生的结点数为1+k n ，则k n 和1+k n 中间的递推关系为341-=+k k n n 。

三、实验程序及注释：
p=[0 0;10 0]; %P 为初始两个点的坐标,第一列为x 坐标,第二列为y 坐标
n=2; %n 为结点数
A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)]; %旋转矩阵
for k=1:4
d=diff(p)/3; %diff 计算相邻两个点的坐标之差,得到相邻两点确定的向量
%则d 就计算出每个向量长度的三分之一,与题中将线段三等分对应 m=4*n-3; %迭代公式
q=p(1:n-1,:); %以原点为起点,前n-1个点的坐标为终点形成向量
p(5:4:m,:)=p(2:n,:); %迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标 p(2:4:m,:)=q+d; %用向量方法计算迭代后处于4k+2位置上的点的坐标
p(3:4:m,:)=q+d+d*A'; %用向量方法计算迭代后处于4k+3位置上的点的坐标
p(4:4:m,:)=q+2*d; %用向量方法计算迭代后处于4k 位置上的点的坐标
n=m; %迭代后新的结点数目
end
plot(p(:,1),p(:,2)) %绘出每相邻两个点的连线
axis([0 10 0 10])
四、实验数据记录：
由第三部分的程序，可得到如下的Koch 分形曲线：
图2
五、注记：
１．参照实验方法，可绘制如下生成元的Koch 分形曲线：
图3
此时，旋转矩阵为：
⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎪⎭⎫ ⎝
⎛-=0110)2cos()2sin()2sin()2cos(ππππA 程序和曲线如下：
p=[0 0;10 0]; %P 为初始两个点的坐标,第一列为x 坐标,第二列为y 坐标
n=2; %n 为结点数
A=[0 -1;1 0]; %旋转矩阵
for k=1:4
d=diff(p)/3; %diff计算相邻两个点的坐标之差,得到相邻两点确定的向量
%则d就计算出每个向量长度的三分之一,与题中将线段三等分对应
m=5*n-4; %迭代公式
q=p(1:n-1,:); %以原点为起点,前n-1个点的坐标为终点形成向量
p(6:5:m,:)=p(2:n,:); %迭代后处于5k+1位置上的点的坐标为迭代前的相应坐标
p(2:5:m,:)=q+d; %用向量方法计算迭代后处于5k+2位置上的点的坐标
p(3:5:m,:)=q+d+d*A'; %用向量方法计算迭代后处于5k+3位置上的点的坐标
p(4:5:m,:)=q+2*d+d*A'; %用向量方法计算迭代后处于5k+4位置上的点的坐标
p(5:5:m,:)=q+2*d; %用向量方法计算迭代后处于5k位置上的点的坐标
n=m; %迭代后新的结点数目
end
plot(p(:,1),p(:,2)) %绘出每相邻两个点的连线
axis([0 10 0 10])
图4
由于中间三分之一部分是一个正方形时，有很多连接的部分。

所以我们将高度压缩到原来的0.7倍，即中间部分为一个长与宽之比为1:0.7的矩形时，得到程序和曲线如下：
p=[0 0;10 0]; %P为初始两个点的坐标,第一列为x坐标,第二列为y坐标
n=2; %n为结点数
A=[0 -1;1 0]; %旋转矩阵
for k=1:4
d=diff(p)/3; %diff计算相邻两个点的坐标之差,得到相邻两点确定的向量
%则d就计算出每个向量长度的三分之一,与题中将线段三等分对应
m=5*n-4; %迭代公式
q=p(1:n-1,:); %以原点为起点,前n-1个点的坐标为终点形成向量
p(6:5:m,:)=p(2:n,:); %迭代后处于5k+1位置上的点的坐标为迭代前的相应坐标
p(2:5:m,:)=q+d; %用向量方法计算迭代后处于5k+2位置上的点的坐标
p(3:5:m,:)=q+d+0.7*d*A'; %用向量方法计算迭代后处于5k+3位置上的点的坐标
p(4:5:m,:)=q+2*d+0.7*d*A'; %用向量方法计算迭代后处于5k+4位置上的点的坐标
p(5:5:m,:)=q+2*d; %用向量方法计算迭代后处于5k位置上的点的坐标
n=m; %迭代后新的结点数目
end
plot(p(:,1),p(:,2)) %绘出每相邻两个点的连线
axis([0 10 0 10])
图５
２．参照实验方法，我们由四边形的四个初始点出发，对于四边形的每条边，生成元如下：
图6
可得到火焰般的图形。

程序和曲线如下：
p=[0 10;10 0;0 -10;-10 0;0 10];
%P为四边形四个顶点的坐标,其中第五个点与第一个点重合,以便于绘图
%第一列为x坐标,第二列为y坐标
n=5; %n为结点数
A=[cos(-pi/3) -sin(-pi/3);sin(-pi/3) cos(-pi/3)]; %旋转矩阵,顺时针旋转60度
for k=1:5
d=diff(p)/3;m=4*n-3; %迭代公式
q=p(1:n-1,:);
p(5:4:m,:)=p(2:n,:);
p(2:4:m,:)=q+d;
p(3:4:m,:)=q+2*d+d*A';
p(4:4:m,:)=q+2*d;
n=m;
end
plot(p(:,1),p(:,2))
axis([-10 10 -10 10])
图7
３．参照实验方法，由下列的生成元，绘制Koch分形曲线：
图8
分析：为了绘图方便，我们将结点数处理一下，把第一次迭代产生的六个点看成十个点，即图中有五条线段（１－２，３－４，５－６，７－８，９－１０），我们将每条线段的每个端点看成新的两个结点，这样我们就可以很方便地用plot绘图了。

程序和曲线如下：
p=[0 0;10 10]; %P为初始两个点的坐标,第一列为x坐标,第二列为y坐标
n=2; %n为结点数
A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)];
B=[cos(-pi/3) -sin(-pi/3);sin(-pi/3) cos(-pi/3)];
%旋转矩阵A对应于第一次逆时针旋转60度,旋转矩阵B对应于第二次顺时针旋转60度
for k=1:4
d=diff(p)/3;
d1=d(1:2:n,:);%取每条线段对应的向量
m=5*n; %迭代公式
q1=p(1:2:n-1,:);
p(10:10:m,:)=p(2:2:n,:);
p(1:10:m,:)=p(1:2:n,:); %迭代后处于10k与10k+1位置上的点的坐标为迭代前的相应坐标p(2:10:m,:)=q1+d1;
%用向量方法计算迭代后处于10k+2,10k+3,10k+5位置上的点的坐标,都相同
p(3:10:m,:)=p(2:10:m,:);
p(4:10:m,:)=q1+d1+d1*A'; %用向量方法计算迭代后处于10k+4位置上的点的坐标
p(5:10:m,:)=p(2:10:m,:);
p(6:10:m,:)=q1+2*d1;
%用向量方法计算迭代后处于10k+6,10k+7,10k+9位置上的点的坐标,都相同
p(7:10:m,:)=p(6:10:m,:);
p(8:10:m,:)=q1+2*d1+d1*B';
p(9:10:m,:)=p(6:10:m,:);
n=m; %迭代后新的结点数目
end
plot(p(:,1),p(:,2)) %绘出每相邻两个点的连线
axis([0 10 0 10])
六，结束语
通过图形显示，更好的理解的分形同时也也加深对分形概念的进一步掌握参考文献：<<matlab从精通到入门>>。