【第4章】分形图程序

（1）Koch曲线程序koch.mfunctionkoch(a1,b1,a2,b2,n)%koch(0,0,9,0,3)%a1,b1,a2,b2为初始线段两端点坐标,n为迭代次数a1=0;b1=0;a2=9;b2=0;n=3;%第i-1次迭代时由各条线段产生的新四条线段的五点横、纵坐标存储在数组A、B中[A,B]=sub_koch1(a1,b1,a2,b2);for

i=1:n

forj=1:length(A)/5;

w=sub_koch2(A(1+5*(j-1):5*j),B(1+5*(j-1):5*j));

fork=1:4

[AA(5*4*(j-1)+5*(k-1)+1:5*4*(j-1)+5*(k-1)+5),BB(5*4*(j-1)+5*(k-1)+1:5*4*(j-1)+5*(k-1)+5)]=sub_koch1(w(k,1),w(k,2),w(k,3),w(k,4));

end

end

A=AA;

B=BB;endplot(A,B)holdonaxisequal

%由以(ax,ay),(bx,by)为端点的线段生成新的中间三点坐标并把这五点横、纵坐标依次分别存%储在数组A,B中function[A,B]=sub_koch1(ax,ay,bx,by)cx=ax+(bx-ax)/3;cy=ay+(by-ay)/3;ex=bx-(bx-ax)/3;ey=by-(by-ay)/3;L=sqrt((ex-cx).^2+(ey-cy).^2);alpha=atan((ey-cy)./(ex-cx));if(ex-cx)<0

alpha=alpha+pi;enddx=cx+cos(alpha+pi/3)*L;dy=cy+sin(alpha+pi/3)*L;A=[ax,cx,dx,ex,bx];B=[ay,cy,dy,ey,by];

%把由函数sub_koch1生成的五点横、纵坐标A,B顺次划分为四组,分别对应四条折线段中%每条线段两端点的坐标,并依次分别存储在4*4阶矩阵k中,k中第i(i=1,2,3,4)行数字代表第%i条线段两端点的坐标functionw=sub_koch2(A,B)a11=A(1);b11=B(1);a12=A(2);b12=B(2);a21=A(2);b21=B(2);a22=A(3);b22=B(3);a31=A(3);b31=B(3);a32=A(4);b32=B(4);a41=A(4);b41=B(4);a42=A(5);b42=B(5);w=[a11,b11,a12,b12;a21,b21,a22,b22;a31,b31,a32,b32;a41,b41,a42,b42];

（2）Levy

曲线程序levy.mfunctionlevy(n)%levy(16),n为levy曲线迭代次数%x1,y1,x2,y2为初始线段两端点坐标,nn为迭代次数n=16;x1=0;y1=0;x2=1;y2=0;%第i-1次迭代时由各条线段产生的新两条线段的三端点横、纵坐标存储在数组X、Y中[X,Y]=levy1(x1,y1,x2,y2);fori=1:n

forj=1:length(X)/3

w=levy2(X(1+3*(j-1):3*j),Y(1+3*(j-1):3*j));

[XX(3*2*(j-1)+1:3*2*(j-1)+3),YY(3*2*(j-1)+1:3*2*(j-1)+3)]=levy1(w(1,1),w(1,2),w(1,3),w(1,4));

[XX(3*2*(j-1)+3+1:3*2*(j-1)+3+3),YY(3*2*(j-1)+3+1:3*2*(j-1)+3+3)]=levy1(w(2,1),w(2,2),w(2,3),w(2,4));

end

X=XX;

Y=YY;endplot(X,Y)holdonaxisequal

%由以(x1,y1),(x2,y2)为端点的线段生成新的中间点坐标并把(x1,y1),(x2,y2)连同新点横、纵坐%标依次分别存储在数组X,Y中function[X,Y]=levy1(x1,y1,x2,y2)x3=1/2*(x1+x2+y1-y2);y3=1/2*(-x1+x2+y1+y2);X=[x1,x3,x2];Y=[y1,y3,y2];

%把由函数levy1生成的三点横、纵坐标X,Y顺次划分为两组,分别对应两条折线段中每条线%段两端点的坐标,并依次分别存储在2*4阶矩阵w中,w中第i(i=1,2)行数字代表第i条线段%两端点的坐标functionw=levy2(X,Y)a11=X(1);b11=Y(1);a12=X(2);b12=Y(2);a21=X(2);b21=Y(2);a22=X(3);b22=Y(3);w=[a11,b11,a12,b12;a21,b21,a22,b22];

（3）

分形树程序tree.hfunctiontree(n,a,b)%tree(8,pi/8,pi/8),n为分形树迭代次数%a,b为分枝与竖直方向夹角%x1,y1,x2,y2为初始线段两端点坐标,nn为迭代次数n=8;a=pi/8;b=pi/8;x1=0;y1=0;x2=0;y2=1;plot([x1,x2],[y1,y2])holdon[X,Y]=tree1(x1,y1,x2,y2,a,b);holdonW=tree2(X,Y);w1=W(:,1:4);w2=W(:,5:8);%w为2^k*4维矩阵,存储第k次迭代产生的分枝两端点的坐标,%w的第i(i=1,2,…,2^k)行数字对应第i个分枝两端点的坐标w=[w1;w2];fork=1:n

fori=1:2^k

[X,Y]=tree1(w(i,1),w(i,2),w(i,3),w(i,4),a,b);

W(i,:)=tree2(X,Y);

end

w1=W(:,1:4);

w2=W(:,5:8);

w=[w1;w2];end

%由每个分枝两端点坐标(x1,y1),(x2,y2)产生两新点的坐标(x3,y3),(x4,y4),画两分枝图形,并把%(x2,y2)连同新点横、纵坐标分别存储在数组X,Y中function[X,Y]=tree1(x1,y1,x2,y2,a,b)L=sqrt((x2-x1)^2+(y2-y1)^2);if(x2-x1)==0

a=pi/2;

elseif(x2-x1)<0

a=pi+atan((y2-y1)/(x2-x1));

else

a=atan((y2-y1)/(x2-x1));

endendx3=x2+L*2/3*cos(a+b);y3=y2+L*2/3*sin(a+b);x4=x2+L*2/3*cos(a-b);y4=y2+L*2/3*sin(a-b);a=[x3,x2,x4];b=[y3,y2,y4];plot(a,b)axisequalholdonX=[x2,x3,x4];Y=[y2,y3,y4];

%把由函数tree1生成的X,Y顺次划分为两组,分别对应两分枝两个端点的坐标,并存储在一维%数组w中functionw=tree2(X,Y)a1=X(1);b1=Y(1);a2=X(2);b2=Y(2);a3=X(1);b3=Y(1);a4=X(3);b4=Y(3);w=[a1,b1,a2,b2,a3,b3,a4,b4];

（4）IFS算法画Sierpinski三角形程序sierpinski_ifs.hfunctionsierpinski_ifs(n,w1,w2,w3)%sierpinski_ifs(10000,1/3,1/3,1/3)%w1,w2,w3出现频率n=10000;w1=1/3;w2=1/3;w3=1/3;M1=[0.50

0

00.5

0];M2=[0.500.5

00.5

0];M3=[0.500.2500.5

0.5];x=0;y=0;%r为[0,1]区间内产生的n维随机数组r=rand(1,n);B=zeros(2,n);k=1;%当0<r(i)<1/3时,进行M1对应的压缩映射;%当1/3=<r(i)<2/3时,进行M2对应的压缩映射;%当2/3=<r(i)<1时,进行M3对应的压缩映射;fori=1:n

ifr(i)<w1

a=M1(1);b=M1(2);e=M1(3);c=M1(4);d=M1(5);f=M1(6);

elseif

r(i)<w1+w2

a=M2(1);b=M2(2);e=M2(3);c=M2(4);d=M2(5);f=M2(6);

elseif

r(i)<w1+w2+w3

a=M3(1);b=M3(2);e=M3(3);c=M3(4);d=M3(5);f=M3(6);

end

end

endx=a*x+b*y+e;y=c*x+d*y+f;B(1,k)=x;B(2,k)=y;k=k+1;end

plot(B(1,:),B(2,:),'.','markersize',0.1)

（5）IFS算法画Julia集程序julia_ifs.hfunctionjulia_ifs(n,cx,cy)%julia_ifs(100000,-0.77,0.08)%f(z)=z^2+c,cx=real(c);cy=image(c);n=10000;cx=-0.77;cy=0.08;%z^2+c=z0,x=real(z0);y=image(z0);x=1;y=1;

B=zeros(2,n);k=1;%A为产生的服从标准正态分布的n维随机数组A=randn(1,n);fori=1:n

wx=x-cx;

wy=y-cy;

ifwx>0

alpha=atan(wy/wx);

end

ifwx<0

alpha=pi+atan(wy/wx);

end

ifwx==0

alpha=pi/2;

end

alpha=alpha/2;

r=sqrt(wx^2+wy^2);

ifA(i)<0

r=-sqrt(r);

else

r=sqrt(r);

end

x=r*cos(alpha);

y=r*sin(alpha);

B(1,k)=x;

B(2,k)=y;

k=k+1;

endplot(B(1,:),B(2,:),'.','markersize',0.1)

(6)

逃逸时间算法画Sierpinski垫片程序sierpinski.hfunctionsierpinski(a,b,c,d,n,m,r)%sierpinski(0,0,1,1,12,200,200)%(a,b),(c,d)收敛区域左上角和右下角坐标,m为分辨率%n为逃逸时间,需要反复试探,r逃逸半径a=0;b=0;c=1;d=1;n=12;m=200;r=200;B=zeros(2,m*m);w=1;fori=1:m

x0=a+(c-a)*(i-1)/m;

forj=1:m

y0=b+(d-b)*(j-1)/m;

x=x0;

y=y0;

fork=1:n

ify>0.5

x=2*x;

y=2*y-1;

elseifx>=0.5

x=2*x-1;

y=2*y;

else

x=2*x;

y=2*y;

end

ifx^2+y^2>r

break;

end

end

ifk==n

B(1,w)=i;

B(2,w)=j;

w=w+1;

end

endendplot(B(1,:),B(2,:),'.','markersize',0.1)

(7)

元胞自动机算法画Sierpinski三角形程序

一维元胞自动机sierpinski_ca1.hfunctionsierpinski_ca1(m,n)%sierpinski_ca1(1000,3000)m=1000;n=3000;x=1;y=1;t=1;w=zeros(2,m*n);s=zeros(m,n);s(1,fix(n/3))=1;fori=1:m-1

forj=2:n-1

if(s(i,j-1)==1&s(i,j)==0&s(i,j+1)==0)|(s(i,j-1)==0&s(i,j)==0&s(i,j+1)==1)

s(i+1,j)=1;

w(1,t)=x+3+3*j;

w(2,t)=y+5*i;

t=t+1;

end

endendplot(w(1,:),w(2,:),'.','markersize',1)

²

二维元胞自动机sierpinski_ca2.hfunctionsierpinski_ca2(m,n)%sierpinski_ca2(400,400)m=400;n=400;t=1;w=zeros(2,m*n);s=zeros(m,n);s(m/2,n/2)=1;fori=[m/2:-1:2,m/2:m-1]

forj=[n/2:-1:2,n/2:n-1]

ifmod(s(i-1,j-1)+s(i,j-1)+s(i+1,j-1)+s(i-1,j)+s(i+1,j)+s(i-1,j+1)+s(i,j+1)+s(i+1,j+1),2)==1

s(i,j)=1;

w(1,t)=i;

w(2,t)=j;

t=t+1;

end

endendplot(w(1,:),w(2,:),'.','markersize',0.1)

(8)

IFS算法画Helix曲线程序helix_ifs.hfunctionhelix_ifs(n,w1,w2,w3)%helix_ifs(20000,0.9,0.05,0.05)%w1,w2,w3为出现频率n=20000;w1=0.9;w2=0.05;w3=0.05;M1=[0.787879

-0.424242

1.7586470.242424

0.859848

1.408065];M2=[-0.1212120.257576

-6.7216540.05303

0.05303

1.377236];M3=[0.181818

-0.136364

6.0861070.090909

0.181818

1.568035];x=0;y=0;%r为[0,1]区间内产生的n维随机数组r=rand(1,n);B=zeros(2,n);k=1;%当0<r(i)<1/3时,进行M1对应的压缩映射;%当1/3=<r(i)<2/3时,进行M2对应的压缩映射;%当2/3=<r(i)<1时,进行M3对应的压缩映射;fori=1:n

ifr(i)<w1

a=M1(1);b=M1(2);e=M1(3);c=M1(4);d=M1(5);f=M1(6);

elseif

r(i)<w1+w2

a=M2(1);b=M2(2);e=M2(3);c=M2(4);d=M2(5);f=M2(6);

elseif

r(i)<w1+w2+w3

a=M3(1);b=M3(2);e=M3(3);c=M3(4);d=M3(5);f=M3(6);

end

end

endx=a*x+b*y+e;y=c*x+d*y+f;B(1,k)=x;

B(2,k)=y;

k=k+1;

end

plot(B(1,:),B(2,:),'.','markersize',0.1)/link?url=CEfKyhT7lmb10gOVen5eSHpjztx2JuN1vJf02qoyK9oXrqfTYETbwxhPSE7QdgQrzqObbSGfbQixLyP3U09MokUvaST2CcsGP4uFqYi6YF_小宇宙

被誉为大自然的几何学的分形（Fractal）理论，是现代数学的一个新分支，但其本质却是一种新的世界观和方法论。它与动力系统的混沌理论交叉结合，相辅相成。它承认世界的局部可能在一定条件下，在某一方面（形态，结构，信息，功能，时间，能量等）表现出与整体的相似性，它承认空间维数的变化既可以是离散的也可以是连续的，因而拓展了视野。

分形几何的概念是美籍法国数学家曼德布罗（B.B.Mandelbrot）1975年首先提出的，但最早的工作可追朔到1875年，德国数学家维尔斯特拉斯（K.Weierestrass）构造了处处连续但处处不可微的函数，集合论创始人康托（G.Cantor，德国数学家）构造了有许多奇异性质的三分康托集。1890年，意大利数学家皮亚诺（G.Peano）构造了填充空间的曲线。1904年，瑞典数学家科赫（H.vonKoch）设计出类似雪花和岛屿边缘的一类曲线。1915年，波兰数学家谢尔宾斯基（W.Sierpinski）设计了象地毯和海绵一样的几何图形。这些都是为解决分析与拓朴学中的问题而提出的反例，但它们正是分形几何思想的源泉。1910年，德国数学家豪斯道夫（F.Hausdorff）开始了奇异集合性质与量的研究，提出分数维概念。1928年布利干（G.Bouligand）将闵可夫斯基容度应用于非整数维，由此能将螺线作很好的分类。1932年庞特里亚金（L.S.Pontryagin）等引入盒维数。1934年，贝塞考维奇（A.S.Besicovitch）更深刻地提示了豪斯道夫测度的性质和奇异集的分数维，他在豪斯道夫测度及其几何的研究领域中作出了主要贡献，从而产生了豪斯道夫-贝塞考维奇维数概念。以后，这一领域的研究工作没有引起更多人的注意，先驱们的工作只是作为分析与拓扑学教科书中的反例而流传开来。

真正令大众了解分形是从计算机的普及肇始，而一开始，分形图的计算机绘制也只是停留在二维平面，但这也足以使人们心驰神往。近来，一个分形体爱好者丹尼尔•怀特（英国一钢琴教师）提出一个大胆的方法，创造出令人称奇的3D分形影像，并将它们命名为芒德球（mandelbulb）。

在芒德球极其繁琐的外表下，这个集合实际上是由一种非常基础的算法得出的。那是一种利用复数的算法。就曼德布罗集而言，它是直接由最简单的乘方运算得出的——对复数进行乘方。但问题在于无法在三维空间恰当地扩展数的概念。与复数和平面点之间的关系不同，19世纪的数学家们曾证明，立体空间中的点是无法用适宜传统加法和乘法运算的代数工具来表示的。既然无法定义数字计算，自然也就无法勾画曼德布罗集的三维形象。解决方案之一是在四维空间中进行计算，然后将结果投射到三维空间中。四维空间中的每个点都可与“四元数”（quaternion）匹配，对它们可以进行传统算术操作。尽管四维空间无法用肉眼看到，但利用四元数便能轻而易举地列出与曼德布罗集相对应的算法，之后去掉一个分量，就能使结果显示成三维效果。但这个方案也令人失望，得到的画面比二维图像好不了多少。

为了避开这个难题，丹尼尔•怀特两年前冒出一个古怪的想法。彻底摆脱数学的羁绊，他在三维空间的点与点之间凭空构建出一种“伪分形”。尽管其处理手段算不上中规中矩的乘法，但至少将与曼德布罗集相对应的算法扩展到了三维空间中所有的点。丹尼尔•怀特对几百万个点进行了计算，之后又追加了光影和纹理以体现立体效果，终于，在他的屏幕上呈现出第一个芒德球，形状与严格的曼德布罗集十分近似。遗憾的是，这一结果没能满足他的期望：“图形令人惊叹，但我期望的是更精致的细节。”

尝试并未就此止步。丹尼尔•怀特在互联网上的一个分形体论坛上引起了美国一位年轻计算机专家保罗•尼兰德的注意。他接手怀特的研究，对算法进行稍事改动，把反复的平方操作换成更高次方（八次方），从而得到了一系列新的芒德球，指数越高，细节就越丰富。

这个芒德球引起了我的极大兴趣，下决心要学学分形体，于是乎决定从最简单的分形算法学起，希望与各位共勉。

以下开始介绍几例最简单的分形算法：

一、Cantor三分集的递归算法

选取一个欧氏长度的直线段，将该线段三等分，去掉中间一段，剩下两段。将剩下的两段分别再三等分，各去掉中间一段，剩下四段。将这样的操作继续下去，直到无穷，则可得到一个离散的点集。点数趋于无穷多，而欧氏长度趋于零。经无限操作，达到极限时所得到的离散点集称之为Cantor集。

1.给定初始直线两个端点的坐标（ax，ay）和（bx，by），按Cantor三分集的生成规则计算出个关键点的坐标如下：

cx=ax+(bx-ax)/3

cy=ay-d

dx=bx-(bx-ax)/3

dy=by-d

ay=ay-d

by=by-d

2.利用递归算法，将计算出来的新点分别对应于（ax，ay）和（bx，by），然后利用步骤1的计算关系计算出下一级新点（cx，cy）和（dx，dy），并压入堆栈。

3.给定一个小量c，当（bx，by）<c时，被压入堆栈中的值依次释放完毕，同时绘制直线段（ax，ay）-（bx，by），然后程序结束。

下面给出matlab程序：

functionf=cantor(ax,ay,bx,by)

c=0.005;d=0.005;

if(bx-ax)>c

x=[ax,bx];y=[ay,by];holdon;

plot(x,y,'LineWidth',2);holdoff;

cx=ax+(bx-ax)/3;

cy=ay-d;

dx=bx-(bx-ax)/3;

dy=by-d;

ay=ay-d;

by=by-d;

cantor(ax,ay,cx,cy);

cantor(dx,dy,bx,by);

end

运行cantor（0,5,5,5），出现图例如下：

二、Koch曲线的递归算法

在一单位长度的线段上对其三等分，将中间段直线换成一个去掉底边的等边三角形，再在每条直线上重复以上操作，如此进行下去直到无穷，就得到分形曲线Koch曲线。

1.给定初始直线（ax，ay）-（bx，by），按Koch曲线的构成原理计算出各关键点坐标如下：

cx=ax+(bx-ax)/3

cy=ay+(by-ay)/3

ex=bx-(bx-ax)/3

ey=by-(by-ay)/3

l=sqrt((ex-cx)^2+(ey-cy)^2)

alpha=atan((ey-cy)/(ex-cx))

dy=cy+sin(alpha+pi/3)*l

dx=cx+cos(alpha+pi/3)*l

2.利用递归算法，将计算出来的新点分别对应于（ax，ay）和（bx，by），然后利用步骤1中的计算公式计算出下一级新点（cx，cy），（dx，dy），（ex，ey），并压入堆栈。

3.给定一个小量c，当l<c时，被压入堆栈中的值依次释放完毕，同时绘制直线段（ax，ay）-（bx，by），然后结束程序。

下面给出matlab程序：

functionf=Koch(ax,ay,bx,by,c)

if(bx-ax)^2+(by-ay)^2<c

x=[ax,bx];y=[ay,by];

plot(x,y);holdon;

else

cx=ax+(bx-ax)/3;cy=ay+(by-ay)/3;

ex=bx-(bx-ax)/3;ey=by-(by-ay)/3;

l=sqrt((ex-cx)^2+(ey-cy)^2);

alpha=atan((ey-cy)/(ex-cx));

if(alpha>=0&(ex-cx)<0)|(alpha<=0&(ex-cx)<0)

alpha=alpha+pi;

end

dy=cy+sin(alpha+pi/3)*l;

dx=cx+cos(alpha+pi/3)*l;

Koch(ax,ay,cx,cy,c);

Koch(ex,ey,bx,by,c);

Koch(cx,cy,dx,dy,c);

Koch(dx,dy,ex,ey,c);

end

运行Koch(0,0,100,0,10)，出现图例如下：

三、生成填充Julia集

1.设定参数a，b以及一个最大的迭代步数N。

2.设定一个限界值R，即实数R≧max（2，sqrt（a^2+b^2）。

3.对于平面上以R为半径的圆盘内的每一点进行迭代，如果对于所有的n≦N，都有|x^2+y^2|≦R，那么，在屏幕上绘制出相应的起始点，否则不绘制。

下面给出matlab程序：

a=-0.11;b=0.65;r=2;

forx0=-1:0.01:1

fory0=-1:0.01:1

x=x0;y=y0;

ifx0^2+y0^2<1

forn=1:80

x1=x*x-y*y+a;

y1=2*x*y+b;

x=x1;

y=y1;

end

if(x*x+y*y)<r

plot(x0,y0);

end

holdon;

end

end

enda=-0.11,b=0.65a=-0.13,b=0.77a=-0.19,b=0.6557

四、牛顿迭代

牛顿迭代是在数值求解非线性方程（组）的时候经常使用的方法。有些牛顿迭代能够绘制出漂亮的图形来，所以现在也常用于设计图形。

Matlab程序如下：

首先编写newton函数：

functiony=newton(z)

if(z==0)

y=0;

return;

end

fori=1:1:2000

y=z-(z^3-1)/(3*z^2);

if(abs(y-z)<1.0e-7)

break;

end

z=y;

end

接着进入主程序：

clearall;clc;

A=1;B=0;C=1;

fora=-1:0.005:1

forb=-1:0.005:1

x0=a+b*i;

y=newton(x0);

ifabs(y-A)<1.0e-6

plot(a,b,'r');holdon;

elseifabs(y-B)<1.0e-6

plot(a,b,'g');holdon;

elseifabs(y-C)<1.0e-6

plot(a,b,'y');holdon;

end

end

end

五、迭代函数系IFS

IFS是分形的重要分支。它是分形图像处理中最富生命力而且最具有广阔应用前景的领域之一。这一工作最早可以追溯到Hutchinson于1981年对自相似集的研究。美国科学家M.F.Barnsley于1985年发展了这一分形构型系统，并命名为迭代函数系统（IteratedFunctionSystem，IFS），后来又由StephenDemko等人将其公式化，并引入到图像合成领域中。IFS将待生成的图像看做是由许多与整体相似的（自相似）或经过一定变换与整体相似的（自仿射）小块拼贴而成。

算法：

1.设定一个起始点（x0，y0）及总的迭代步数。

2.以概率P选取仿射变换W，形式为

X1=ax0+by0+e

Y1=cx0+dy0+f

3.以W作用点（x0，y0），得到新坐标（x1，y1）。

4.令x0=x1，y0=y1。

5.在屏幕上打出（x0，y0）。

6.重返第2步，进行下一次迭代，直到迭代次数大于总步数为止。

下面给出一些IFS植物形态的matlab程序：

a=[0.195-0.4880.3440.4330.44310.24520.25;

0.4620.414-0.2520.3610.25110.56920.25;

-0.058-0.070.453-0.1110.59760.09690.25;

-0.0350.07-0.469-0.0220.48840.50690.2;

-0.637000.5010.85620.25130.05];

x0=1;y0=1;

fori=1:10000

r=rand;

ifr<=0.25

x1=a(1,1)*x0+a(1,2)*y0+a(1,5);

y1=a(1,3)*x0+a(1,4)*y0+a(1,6);

end

ifr>0.25&r<=0.5

x1=a(2,1)*x0+a(2,2)*y0+a(2,5);

y1=a(2,3)*x0+a(2,4)*y0+a(2,6);

end

ifr>0.5&r<=0.75

x1=a(3,1)*x0+a(3,2)*y0+a(3,5);

y1=a(3,3)*x0+a(3,4)*y0+a(3,6);

end

ifr>0.75&r<=0.95

x1=a(4,1)*x0+a(4,2)*y0+a(4,5);

y1=a(4,3)*x0+a(4,4)*y0+a(4,6);

end

ifr>0.95&r<=1

x1=a(5,1)*x0+a(5,2)*y0+a(5,5);

y1=a(5,3)*x0+a(5,4)*y0+a(5,6);

end

x0=x1;y0=y1;

plot(x1,y1);holdon;

end

得到图例如下：

修改部分系数便可得到另一种形态：

六、三角形分形

functiontriangles(n);

clc;closeall;

ifnargin==0;

n=4;

end

rand('state',2);

C=rand(n+4,3);

figure;

axissquareequal;holdon;

a=-pi/6;

p=0;

r=1;

[p,r,n,a]=tritri(p,r,n,a,C);

function[p,r,n,a]=tritri(p,r,n,a,C);

%画一个三角形

%p是三角形中心

%r是三角形半径

%n是递归次数

%a是三角形角度

%C是颜色矩阵

z=p+r*exp(i*([0:3]*pi*2/3+a));

zr=p+r*exp(i*([0:3]*pi*2/3+a))/2;

pf=fill(real(z),imag(z),C(n+2,:));

set(pf,'EdgeColor',C(n+2,:));

ifn>0;

[p,r,n,a]=tritri(p,r/2,n-1,a+pi/3,C);

n=n+1;r=r*2;a=a-pi/3;

[zr(1),r,n,a]=tritri(zr(1),r/4,n-1,a,C);

n=n+1;r=r*4;

[zr(2),r,n,a]