分形图程序、分形图基本图形以及源程序、分形图形学实验指导

（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]=tritri(zr(2),r/4,n-1,a,C);

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

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

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

end

七、曼德布罗集合

Mandelbrotset是在复平面上组成分形的点的集合。Mandelbrot集合可以用复二次多项式f(z)=z^2+c来定义。其中c是一个复参数。对于每一个c，从z=0开始对f(z)进行迭代序列(0,f(0),f(f(0)),f(f(f(0))),.......)的值或者延伸到无限大，或者只停留在有限半径的圆盘内。曼德布罗集合就是使以上序列不延伸至无限大的所有c点的集合。从数学上来讲，曼德布罗集合是一个复数的集合。一个给定的复数c或者属于曼德布罗集合M，或者不是。

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

2.设定一个限界值R，不妨设实数R=2。

3.对于参数平面上的每一点c（a，b），使用以R为半径的圆盘内的每一点进行迭代，如果对于所有的n≤N，都有|x*x+y*y|≤R*R，那么，在屏幕上绘制出相应的起始点c（a，b），否则不绘制。

下面给出matlab程序：

r=4;%限界值

fora=-2:0.002:1

forb=-2:0.002:1%参数a，b取到一个范围

x=a;y=b;%初始的复数c

forn=1:20

x1=x*x-y*y+a;%复数平方加一个c的运算

y1=2*x*y+b;

x=x1;%迭代

y=y1;

end

if(x*x+y*y)<r%限界

plot(a,b);

end

holdon;

end

end

八、脑分形

作为IFS的一种应用

a=[0.03000.45000.05;

-0.0300-0.4500.40.15;

0.56-0.560.560.5600.40.4;

0.560.56-0.560.5600.40.4];

x0=1;y0=1;

fori=1:100000

r=rand;

ifr<=0.05

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.05&r<=0.2

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.2&r<=0.6

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.6&r<=1

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

x0=x1;y0=y1;

plot(x1,y1);

holdon;

end最简单的brainfractal分形图基本图形以及源程序第一部分本人新手，如有错误请指正。程序完成于2011/6/17晚间到2011/6/18。很多变量名称采用的是同学的姓名拼音，为的是告诉大家这些都是可以随意命名的变量或函数名，一般大写字母开头的是系统定义的变量不可以随意更改。(*雪花*)源程序lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=60Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,e,T={{ca,-sa},{sa,ca}}},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*2/3+zhengguojie[[i+1]]/3;e=zhengguojie[[i]]/3+zhengguojie[[i+1]]*2/3;d=c+T.(e-c);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,e,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1/2,Sqrt[3]/2},{1,0},{0,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,0]],AspectRatioSqrt[3]/2]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,5]],AspectRatioSqrt[3]/2]]基本生成元分形图(*上三角下三角1*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=60Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,g,f,e,T={{ca,-sa},{sa,ca}},S={{ca,sa},{-sa,ca}}},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*3/4+zhengguojie[[i+1]]/4;e=zhengguojie[[i]]/2+zhengguojie[[i+1]]/2;f=zhengguojie[[i]]/4+zhengguojie[[i+1]]*3/4;d=c+T.(e-c);g=e+S.(f-e);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,e,g,f,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatioSqrt[3]/2]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,4]],AspectRatioSqrt[3]/2]]基本生成元分形图(*上三角下三角2*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=60Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,g,f,e,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*3/4+zhengguojie[[i+1]]/4;e=zhengguojie[[i]]/2+zhengguojie[[i+1]]/2;f=zhengguojie[[i]]/4+zhengguojie[[i+1]]*3/4;d=c+T.(e-c);S=Transpose[T];g=e+S.(f-e);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,e,g,f,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatioSqrt[3]/2]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,4]],AspectRatioSqrt[3]/2]]图形同上(*上三角下三角3*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=90Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,d1,g1,g,f,e,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*3/4+zhengguojie[[i+1]]/4;e=zhengguojie[[i]]/2+zhengguojie[[i+1]]/2;f=zhengguojie[[i]]/4+zhengguojie[[i+1]]*3/4;d=c+T.(e-c);d1=d+e-c;S=Transpose[T];g=e+c-d;g1=g+f-e;weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,d1,e,g,g1,f,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatio1/2]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,4]],AspectRatio1/2]]图形同上(*上正方下正方形*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=90Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,d1,g1,g,f,e,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*3/4+zhengguojie[[i+1]]/4;e=zhengguojie[[i]]/2+zhengguojie[[i+1]]/2;f=zhengguojie[[i]]/4+zhengguojie[[i+1]]*3/4;d=c+T.(e-c);d1=e+T.(f-e);S=Transpose[T];g=e+S.(f-e);g1=f+S.(f-e);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,d1,e,g,g1,f,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatio1/2]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,4]],AspectRatio1/2]]基本生成元分形图(*下正方形上正方形*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=90Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,d1,g1,g,f,e,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*3/4+zhengguojie[[i+1]]/4;e=zhengguojie[[i]]/2+zhengguojie[[i+1]]/2;f=zhengguojie[[i]]/4+zhengguojie[[i+1]]*3/4;S=Transpose[T];d=c+S.(e-c);d1=e+S.(f-e);g=e+T.(f-e);g1=f+T.(f-e);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,d1,e,g,g1,f,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatio1/2]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,4]],AspectRatio1/2]]基本生成元分形图(*下三角上三角*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=60Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,g,f,e,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*3/4+zhengguojie[[i+1]]/4;e=zhengguojie[[i]]/2+zhengguojie[[i+1]]/2;f=zhengguojie[[i]]/4+zhengguojie[[i+1]]*3/4;S=Transpose[T];d=c+S.(e-c);g=e+T.(f-e);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,e,g,f,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatioSqrt[3]/2]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,5]],AspectRatioSqrt[3]/2]]基本生成元分形图(*单个上正方*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=90Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,d1,g1,g,f,e,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*2/3+zhengguojie[[i+1]]/3;e=zhengguojie[[i]]/3+zhengguojie[[i+1]]*2/3;S=Transpose[T];d=c+T.(e-c);d1=e+T.(e-c);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,d1,e,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,7]],AspectRatio1/3]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatio1/3]]基本生成元分形图(*一个正方形向外长大*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=90Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,d1,g1,g,f,e,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*2/3+zhengguojie[[i+1]]/3;e=zhengguojie[[i]]/3+zhengguojie[[i+1]]*2/3;S=Transpose[T];d=c+T.(e-c);d1=e+T.(e-c);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,d1,e,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0},{1,-1},{0,-1},{0,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,4]],AspectRatio1/1]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatio1/1]]基本生成元分形图(*一个正方形向内长大*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=90Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,d1,g1,g,f,e,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*2/3+zhengguojie[[i+1]]/3;e=zhengguojie[[i]]/3+zhengguojie[[i+1]]*2/3;S=Transpose[T];d=c+T.(e-c);d1=e+T.(e-c);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,d1,e,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0},{1,1},{0,1},{0,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,5]],AspectRatio1/1]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatio1/1]]基本生成元分形图(*一个M形状图形*)lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=60Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,g,f,e,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*3/4+zhengguojie[[i+1]]/4;e=zhengguojie[[i]]/2+zhengguojie[[i+1]]/2;f=zhengguojie[[i]]/4+zhengguojie[[i+1]]*3/4;S=Transpose[T];d=c+T.(e-c);g=e+T.(f-e);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,e,g,f,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatioSqrt[3]/8]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,5]],AspectRatioSqrt[3]/8]]基本生成元分形图两个上三角形横线lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=60Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,g,f,e,e1,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*4/5+zhengguojie[[i+1]]/5;e=zhengguojie[[i]]*3/5+zhengguojie[[i+1]]*2/5;e1=zhengguojie[[i]]*2/5+zhengguojie[[i+1]]*3/5;f=zhengguojie[[i]]/5+zhengguojie[[i+1]]*4/5;S=Transpose[T];d=c+T.(e-c);g=e1+T.(f-e1);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,e,e1,g,f,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatioSqrt[3]/10]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,5]],AspectRatioSqrt[3]/10]]生成元分形图上三角形横线下三角形lovelyduwangen[zhengguojie_List]:=Block[{weihuayan={},i,wuxiaonan=Length[zhengguojie],gengping=60Degree,sa=Sin[gengping],ca=Cos[gengping],c,d,g,f,e,e1,T={{ca,-sa},{sa,ca}},S},For[i=1,i<wuxiaonan,i++,c=zhengguojie[[i]]*4/5+zhengguojie[[i+1]]/5;e=zhengguojie[[i]]*3/5+zhengguojie[[i+1]]*2/5;e1=zhengguojie[[i]]*2/5+zhengguojie[[i+1]]*3/5;f=zhengguojie[[i]]/5+zhengguojie[[i+1]]*4/5;S=Transpose[T];d=c+T.(e-c);g=e1+S.(f-e1);weihuayan=Join[weihuayan,{zhengguojie[[i]],c,d,e,e1,g,f,zhengguojie[[i+1]]}]];weihuayan]dongquanfa={{0,0},{1,0}};Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,1]],AspectRatioSqrt[3]/5]]Show[Graphics[Line[Nest[lovelyduwangen,dongquanfa,5]],AspectRatioSqrt[3]/5]]生成元分形图分形图基本图形以及源程序第二部分该程序写的太复杂的,至少看着程序段太多,在本文档系列第三部分我会给一个简单的程序实现它.注意,从word直接拷贝到mathematica,箭头会乱码，请自己改一改如这里的箭头AspectRatio1/GoldenRatio挖空一个黑色三角形sierpinski[tris_List]:=Block[{tmp={},i,p=Length[tris]/3,a,b,c,d,e,f},For[i=0,i<p,i=i+1,a=tris[[3i+1]];b=tris[[3i+2]];c=tris[[3i+3]];d=(a+b)/2;e=(a+c)/2;f=(b+c)/2;tmp=Join[tmp,{a,d,e,d,b,f,e,f,c}]];tmp]Showsierpinski[pts_List]:=Block[{tmp={},i,p=Length[pts]/3},For[i=0,i<p,i=i+1,AppendTo[tmp,Polygon[{pts[[3i+1]],pts[[3i+2]],pts[[3i+3]]}]]];Show[Graphics[tmp],AspectRatio1/GoldenRatio]]triangle={{-1,0},{1,0},{0,Sqrt[3]}}p0=Showsierpinski[Nest[sierpinski,triangle,0]]p1=Showsierpinski[Nest[sierpinski,triangle,1]]p2=Showsierpinski[Nest[sierpinski,triangle,2]]p3=Showsierpinski[Nest[sierpinski,triangle,3]]p4=Showsierpinski[Nest[sierpinski,triangle,4]]Show[GraphicsArray[{{p1,p2},{p3,p4}}]];生成元分形图挖空一个彩色的三角形sierpinski[tris_List]:=Block[{tmp={},i,p=Length[tris]/4,a,b,c,d,e,f},For[i=0,i<p,i=i+1,a=tris[[4i+1]];b=tris[[4i+2]];c=tris[[4i+3]];d=(a+b)/2;e=(a+c)/2;f=(b+c)/2;tmp=Join[tmp,{a,d,e,a,d,b,f,d,e,f,c,e}]];tmp]Showsierpinski[pts_List]:=Block[{tmp={},i,p=Length[pts]/4},For[i=0,i<p,i=i+1,(*AppendTo[tmp,{RGBColor[1,0,0],Polygon[{pts[[4i+1]],pts[[4i+2]],pts[[4i+3]]}]}],*)AppendTo[tmp,{Thickness[.02],RGBColor[0,0,1],Line[{pts[[4i+1]],pts[[4i+2]],pts[[4i+3]],pts[[4i+1]]}]}]];Show[Graphics[tmp]]]Showsierp[pts_List]:=Block[{tmp={},i,p=Length[pts]/4},For[i=0,i<p,i=i+1,AppendTo[tmp,{RGBColor[1,0,0],Polygon[{pts[[4i+1]],pts[[4i+2]],pts[[4i+3]],pts[[4i+1]]}]}]];Show[Graphics[tmp]]]triangle={{-1,0},{1,0},{0,Sqrt[3]},{-1,0}}p0=Showsierpinski[Nest[sierpinski,triangle,2]]p1=Showsierp[Nest[sierpinski,triangle,2]]Show[{p0,p1}]生成元分形图第一题的逆,填充挖去的部分彩色swwwierpinski[wtris_List]:=Block[{tmp={},duu={},mm={{},{}},i,j,p=Length[wtris[[1]]]/4,dd=Length[wtris[[2]]],a,b,c,d,e,f},For[j=1,j<dd+1,j=j+1,mm[[2]]=Join[mm[[2]],{wtris[[2]][[j]]}]];For[i=0,i<p,i=i+1,a=wtris[[1]][[4i+1]];b=wtris[[1]][[4i+2]];c=wtris[[1]][[4i+3]];d=(a+b)/2;e=(a+c)/2;f=(b+c)/2;mm[[1]]=Join[mm[[1]],{a,d,e,a,d,b,f,d,e,f,c,e}];mm[[2]]=Join[mm[[2]],{d,e,f,d}]];mm]Shwwowsierp[pts_List]:=Block[{tmp={},i,p=Length[pts[[2]]]/4},For[i=0,i<p,i=i+1,AppendTo[tmp,{RGBColor[1,0,0],Polygon[{pts[[2]][[4i+1]],pts[[2]][[4i+2]],pts[[2]][[4i+3]],pts[[2]][[4i+4]]}]}]];Show[Graphics[tmp]]]triangle={{{-1,0},{1,0},{0,Sqrt[3]},{-1,0}},{}};p9=Shwwowsierp[Nest[s