数学实验课件：迭代与分形

迭代与分形

在我们所处的世界里，存在着许多极不规则的复杂现象，比如：弯弯曲曲的海岸线、变化的云朵、宇宙中星系的分布、金融市场上价格的起伏图等，为了获得解释这些极端复杂现象的数学模型，我们需要认识其中蕴涵的特性，构造出相应的数学规则.曼德尔布罗特（Mandelbrot）在研究英国的海岸线形状等问题时，总结出自然界中很多现象从标度变换角度表现出对称性.

Mandelbrot将这类几何形体称为分形(fractal)，意思就是不规则的、分数的、支离破碎的，并对它们进行了系统的研究，创立了分形几何这一新的数学分支.Mandelbrot认为海岸、山峦、云彩和其他很多自然现象都具有分形的特性，因此可以说：分形是大自然的几何学.

分形几何体一般来说都具有无限精细的自相似的层次结构，即局部与整体的相似性，图形的每一个局部都可以被看作是整体图形的一个缩小的复本.

早在19世纪就已经出现了一些具有自相似特性的分形图形，比如：瑞典数学家科赫（vonKoch）设计的类似雪花和岛屿边缘的一类曲线，即Koch曲线（见图14-1）；图14-1Koch曲线

英国植物学家布朗通过观察悬浮在水中的花粉的运动轨迹，提出来的布朗运动轨迹（见图14-2）.

分形几何把自然形态看作是具有无限嵌套的层次结构，并且在不同尺度下保持某种相似的属性，于是，简单的迭代过程，就是描述复杂的自然形态的有效方法.图14-2布朗运动轨迹

本章简单介绍迭代和分形的基础知识，利用MATLAB绘制Koch曲线、Sierpinski地毯、Mandelbrot集和分形树枝，通过观察这些图形来了解数学之美.CONTENTS14.1迭代14.2.分形10.3MATLAB求解14.1迭代

迭代法是常用的一种数学方法，就是将一种规则反复作用在某个对象上，它可以产生非常复杂的行为.我们这里介绍图形迭代和函数迭代两种方式.

(1)图形迭代.给定初始图形F0，以及一个替换规则R,

将R反复作用在初始图形F0上,产生一个图形序列：

R(F0)＝F1，R(F1)＝F2，R(F2)＝F3，…

科赫（Koch）曲线是通过图形迭代的方式产生的，其迭代规则是：对一条线段，首先将它分成三等份，然后将中间的一份替换成以此为底边的等边三角形的另外两条边.

无限次迭代下去，最终形成的曲线就是Koch曲线（见图14-3）.a）迭代1次c）迭代5次b）迭代2次图14-3Koch曲线

Sierpinski地毯也是通过图形迭代的方式产生的，其迭代规则是：对一个正方形，首先将它分成九个小正方形，然后挖掉中间的一个.

无限次迭代下去，最终形成的图形就是Sierpinski地毯（见图14-4）.a）迭代1次c）迭代4次b）迭代2次图14-4Sierpinski地毯

计算Sierpinski地毯的面积：

初始面积为1；迭代一次后减少的面积为

，剩余

；以后每次迭代，面积减少

，剩余面积是原来的

倍；故经过

n

次迭代后，剩余面积为

.

故Sierpinski地毯面积趋于零..

(2)函数迭代

给定初始值x0，以及一个函数f(x),将f(x)反复作用在初始值x0上，产生一个数列：

f(x0)＝x1，f(x1)＝x2，f(x2)＝x3，…称为迭代序列.

设迭代关系为zk+1=zk2+c，其中zk(k=0,1,2,...)、c为复数.给定初值z0和c可得复数迭代列{zk}，(k=0,1,2,...).若固定，可得Julia集；固定z0，可得Mandelbrot集.14.2分形

分形几何的概念是美籍法国数学家曼德尔布罗特（Mandelbrot）于1975年首先提出的，但最早的工作可追朔到1875年，德国数学家维尔斯特拉斯（Weierestrass）构造了处处连续但处处不可微的函数，集合论创始人康托尔（Cantor，德国数学家）构造了有许多奇异性质的康托尔三分集.

1890年，意大利数学家皮亚诺（Peano）构造了填充空间的曲线.1915年，波兰数学家谢尔宾斯基（Sierpinski）设计了像地毯和海绵一样的几何图形.这些都是为解决分析与拓扑学中的问题而提出的反例，但它们正是分形几何思想的源泉.

分形树枝是为了模拟自然界中树木花草的形状，其迭代方式与前面两个问题不同，主要原理是设定基本的绘图规则，然后让计算机根据这些规则进行反复迭代，最终生成分形图.

下面介绍一种迭代规则：对一条线段进行三等分，得到内部两个分点，在这两个分点处（以分点为起点）分别以旋转角

左右生长出新线段（新线段长度可以为等分后线段的长度）.

无限次迭代下去，最终形成的图形就像一棵树.14.3MATLAB求解例14.1

利用MATLAB绘制科赫（Koch）曲线.解

以下的代码是绘制迭代k次的Koch曲线图形的函数：functionplotkoch(k)%显示迭代k次后的Koch曲线图p=[00;100];%存放结点坐标，每行一个点，初始值为两结点的坐标n=1;%存放线段的数量，初始值为1A=[cos(pi/3)-sin(pi/3);sin(pi/3)cos(pi/3)];%旋转矩阵，用于计算新的结点fors=1:k%实现迭代过程，计算所有的结点的坐标j=0;%以下根据线段两个结点的坐标，计算迭代后它们之间增加的三个%结点的坐标，并且将这些点的坐标按次序存暂时放到r中fori=1:n%每条边计算一次

q1=p(i,:);%目前线段的起点坐标q2=p(i+1,:);%目前线段的终点坐标d=(q2-q1)/3;j=j+1;r(j,:)=q1;%原起点存入rj=j+1;r(j,:)=q1+d;%新1点存入rj=j+1;r(j,:)=q1+d+d*A';%新2点存入rj=j+1;r(j,:)=q1+2*d;%新3点存入rend%原终点作为下条线段的起点，在迭代下条线段时存入rn=4*n;%全部线段迭代一次后，线段数量乘4

clearp%清空p，注意：最后一个终点q2不在r中p=[r;q2];%重新装载本次迭代后的全部结点endplot(p(:,1);p(:,2))%显示各结点的连线图axisequal%各坐标轴同比例

在命令行窗口调用plotkoch(1)，图像为图14-3a；调用plotkoch(2)，图像为图14-3b；调用plotkoch(5)，图像为图14-3c.例14.2利用MATLAB绘制谢尔宾斯基（Sierpinski）地毯.解以下的代码是绘制迭代n次的Sierpinski地毯图形的函数：functionplotSierpinski(x,y,d,n)%x为正方形的顶点的横坐标，可取0（一个顶点代表一个小正方形）%y为正方形的顶点的纵坐标，可取0%d为初始正方形边长，可取1%n为迭代次数，可取4forp=1:n;%实现迭代过程,计算所有的顶点坐标a1=[];%保存迭代后所有顶点的x坐标b1=[];%保存迭代后所有顶点的y坐标%根据小正方形的顶点坐标，%计算迭代后形成的8个新的小正方形的顶点坐标forq=1:length(x);%每个小正方形计算一次x1=x(q)+[0,d/3,2*d/3,0,2*d/3,0,d/3,2*d/3];%新的x坐标y1=y(q)+[0,0,0,d/3,d/3,2*d/3,2*d/3,2*d/3];%新的y坐标

a1=[a1,x1];%所有顶点x坐标存入a1b1=[b1,y1];%所有顶点y坐标存入b1endd=d/3;%迭代一次，边长缩小x=a1;%全部的x坐标重新放入xy=b1;%全部的y坐标重新放入yendholdon%在同一个图形窗口显示forq=1:length(x);%用蓝色注满多边形区域fill(x(q)+[0,d,d,0,0],y(q)+[0,0,d,d,0],'b')endholdoffaxisoff%不要坐标轴axisequal%各坐标轴同比例%不显示这些正方形的边界set(findobj(gcf,'type','patch'),'edgecolor','none')

在命令行窗口调用plotSierpinski(0,0,1,1)，图像为图14-4a；调用plotSierpinski(0,0,1,2)，图像为图14-4b；调用plotSierpinski(0,0,1,4)，图像为图14-4c.例14.3绘制分形——Mandelbrot集.解首先编写函数文件Mandelbrot.m.MATLAB代码如下：functionMandelbrot(a,b,M,n,cx,cy,zm)%b*a为图像网格，M为阈值，n是迭代次数，（cx,cy）是图像中心，zm是放缩倍数delta=2/zm;xl=cx-delta;xr=cx+delta;yd=cy-delta;yu=cy+delta;x=linspace(xl,xr,a);y=linspace(yd,yu,b);[X,Y]=meshgrid(x,y);Z=X+Y*i;P=zeros(b,a);C=Z;fork=1:nZ=Z.^2+C;P(abs(Z)>M)=k;Z(abs(Z)>M)=0;C(abs(Z)>M)=0;endimshow(P,[])在命令行窗口调用Mandelbrot(500,500,4,200,0,0,1)，可得图形如图14-5所示.图14-5Mandelbrot集例14.4利用MATLAB绘制分形树枝.解首先编写函数文件plottree.m.MATLAB代码如下：functionplottree(p,theta,n)m=2;plot(p(:,1),p(:,2),'k')holdonA=[cos(theta)sin(theta);-sin(theta)cos(theta)];%变换矩阵fork=1:ni=1;forj=1:2:mp1=p(j,:);p2=p(j+1,:);d=(p2-p1)/3;w(i,:)=p1;i=i+1;q1=p1+d;w(i,:)=q1;i=i+1;w(i,:)=q1;i=i+1;q2=q1+d*A;w(i,:)=q2;i=i+1;w(i,:)=q1;i=i+1;q3=p1+2*d;w(i,:)=q3;i=i+1;w(i,:)=q3;i=i+1;q4=q3+d*A';w(i,:)=q4;i=i+1;w(i,:)=q3;i=i+1;w(i,:)=p2;i=i+1;point=[q1;q2];plot(point(:,1),point(:,2),'k');point=[q3;q4];plot(point(:,1),