基于MATLAB实现分形图形的绘制论文.doc

目 录前 言1第一章 MATLAB介绍21.1 MATLAB简介21.2 MATLAB语言21.2.1 创建向量、向量元素的访问：21.2.2 创建矩阵、矩阵元素的访问31.2.3 流程控制41.3 MATLAB语言的传统优点5第二章 分形入门知识62.1 分形理论62.2 分形几何观及其应用7第三章 Koch雪花的绘制83.1 von Koch曲线简介83.2 Koch雪花算法设计9第四章 Frac_tree绘制11第五章 Mandelbort集的绘制135.1 Mandelbort集简介135.2 Mandelbort集算法设计13第六章 Julia集的绘制186.1 Julia集简介186.2 Julia集的算法设计186.3 Julia集与Mandelbort集20第七章 花篮簇的绘制22总 结23主要参考文献：23 前 言分形是描述不规则几何形态的有力工具。不言而喻，不规则的几何形态在我们的周围处处可见，诸如花草、山脉、烟云、火焰等举目皆是。至于微观世界的复杂物质结构，宏观世界浩瀚天体的演变，更展现出了层出不穷的不规则几何形态，它们往往都是分形几何的研究对象。大自然向人类展示其美丽多变形态的同时，也提出了难以回答的询问：怎样描述复杂的自然表象？恰恰是分形几何学，它把自然形态看作是具有无限嵌套层次的精细结构，并且在不同的尺度下保持某种相似的属性，于是在变换与迭代中得到描述自然形态的有效方法。分形的研究离不开计算机。如果不是计算机图形图像处理功能的增强，不能想象怎样才能直观地看到Julia集和Mandelbort集的精细结构，更不能想象可以产生具有无限细结的自然景物和高度真实感的三维动画。反过来，分形理论与方法又极大地丰富了计算机图形学内容，甚至分形的思想会在计算机科学的发展上产生一定的影响。下面我们将以MATLAB为平台，讲述用M语言实现分形图形的绘制。在毕业设计中我主要负责绘制图形的程序部分，由负责其界面制作、程序与界面接口调试工作。在设计过程中，遇到了许多问题、困难，但是，在我的导师蒲俊教授和我的搭档黄娟同学一如既往的帮助和鼓励下，我们终于完成了毕业设计。在此，特别向我的导师及搭档致以诚挚的谢意。第一章 MATLAB介绍1.1 MATLAB简介MATLAB是由美国的Math Work公司推出的一个为科学和工程计算而专门设计的高级交互式软件包。它是一种高性能的用于工程计算的编程软件，它把科学计算、编程和结果的可视化都集中在一个使用非常方便的环境中。MATLAB最早是做为矩阵实验室（Matrix Laboratory），20世纪70年代后期，时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库程序的矩阵软件接口，此即用FORTRAN编写的萌芽状态的MATLAB。当其以商品形式出现后，短短几年，就以其良好的开放性和运行的可靠性，使原先控制领域里的封闭式软件包纷纷淘汰，而改以MATLAB为平台加以重建。进入20世纪90年代， MATLAB已经成为国际控制界公认的标准计算软件。Math Works公司于1993年推出MATLAB4.0版本，从此告别DOS版。从1997年，MATLAB5.0版问世到现在最新版本MATLAB6.5.1，MATLAB的功能不断得到加强，使其在科学计算、控制系统设计与分析、数学信号处理、数字图像处理、通讯系统仿真与设计等领域的应用独领风骚。其典型的应用范围包括以下方面：(1) 数学计算(2) 算法推导(3) 建模和仿真模拟(4) 工程图形的绘制(5) 应用程序的开发，包括用户图形界面的建立1.2 MATLAB语言1.2.1 创建向量、向量元素的访问：向量的创建方法主要有下面两种：第一种:逐个输入元素，彼此元素之间可以使用空格或者逗号作为间隔符。如x=1 3 pi 3+5i或者x=1,3,pi,3+5i形式；第二种:利用冒号运算符创建向量,基本语法为X=J:INC:K，其中J为向量的第一个元素，而K为向量的最后一个元素，INC为向量元素递增的步长。若在表达式中忽略INC，则默认的递增步长为1。INC可以为正数也可以为负数，如果INC为正数，则必须JK，否则创建的为空向量。访问向量的元素只要使用相应元素的索引即可。1. A(n) 访问向量的第n个元素；2. A(i j k)或者 A(i,j,k) 访问向量的第i、j、k个元素；3. A(j:k,k:-1:j) 重复访问向量中的元素。1.2.2 创建矩阵、矩阵元素的访问这里，由于输入矩阵的元素、行、列较少，我们采用直接输入矩阵元素法来创建矩阵。如输入A=1 2 3;4 5 6;7 8 9则在创建过程中，矩阵的元素行与行之间需要用分号“;”间隔，也可以在需要分行的地方用回车键间隔，矩阵的元素之间可以使用逗号“,”或者空格间隔。其实创建上面的矩阵还可以这么做A=1:3;4:6;7:9。我们可以将矩阵的每一行或者每一列看作一个向量，矩阵就是由若干向量组合而成的。那么矩阵元素的访问则类似于向量元素的访问，使用索引即可。表 1.1 使用索引访问矩阵元素的方法矩阵元素的访问说 明A(i,j)访问矩阵A的第i行第j列上的元素，其中i和j为标量A(I,J)访问有向量I和J指定的矩阵A中的元素A(i,:)访问矩阵A中第i行的所有元素A(:,j)访问矩阵A中第j列的所有元素A(:)访问矩阵A中的所有元素，将矩阵看作一个向量A(l)使用单下标的方式访问矩阵元素，其中l为标量A(L)访问有向量L指定的矩阵A的元素，向量L中的元素为矩阵元素的单下标数值矩阵元素的单下标与全坐标之间的转换关系如下，以m*n的矩阵为例，该矩阵的第i行第j列的元素全下标表示为单下标l=(j-1)*m+i。为了方便全下标和单下标之间的转换，MATLAB提供了两个函数分别完成两者之间的相互转换：sub2ind：根据全下标计算单下标。ind2sub：根据单下标计算全下标。1.2.3 流程控制程序流程控制包含控制程序流程的基本结构和语法，例如应用程序的选择和循环结构，这也是结构化编程的基本基本结构，使用结构化的应用程序设计方法可以使设计的程序结构清晰、可读性强，能够提高应用程序的设计效率，增强程序的可维护性。结构化的程序主要有三种基本的程序结构：顺序结构、选择结构、循环结构。本节主要介绍选择结构中的if语句，循环结构中的for循环。1. if语句的基本语法结构有三种，分别如下： (1) if(关系运算表达式) MATLAB语句 end这种形式的选择结构表示，当关系运算表达式计算的结果为逻辑真的时候，执行MATLAB语句，这里的MATLAB语句也可以是多个MATLAB表达式。(2) if(关系运算表达式) MATLAB语句AelseMATLAB语句Bend这种形式的选择结构表示，当关系运算表达式的计算结果为逻辑真的时候，执行MATLAB语句A，否则执行MATLAB语句B。(3) if(关系表达式a) MATLAB语句Aelseif(关系运算表达式b) MATLAB语句Belse(关系表达式c)end这种选择结构可以判断多条关系运算表达式的计算结果，然后按照执行的逻辑关系执行相应的语句。2. for语句构成循环的循环是最灵活、简洁的方法，不过使用for语句循环需要预先知道循环体执行的次数。在MATLAB中for循环的基本结构如下：for index=start:increment:end MATLAB语句end其中，index的取值取决于start和end的值，increment为步长，默认的数值为1。和其他高级语言类似，MATLAB的循环结构也可以进行嵌套使用，在使用时需要注意for关键字和end关键字之间的配对使用。1.3 MATLAB语言的传统优点MATLAB自问世起，就以数值计算称雄。MATLAB进行数值计算的基本处理单位是复数数组，并且数组维数是自动按照规则确定的。这一方面使MATLAB程序可以被高度“向量化”，另方面使用户易读易写。比如已知t的采样数据是(n*m)维数组，要计算。对一般的计算语言来说，必须采用两层循环才能得到结果，不但程序复杂，而且那些讨厌的循环十分费时。MATLAB处理这类问题则简洁快捷的多，只需直接了当的一条指令，就可获得同样结果。这就是所谓的“数组运算”。这种运算在下面的分形图绘制过程中，将被频繁使用。第二章 分形入门知识2.1 分形理论谁创立了分形几何学1973年，美籍法国数学家B.Mandelbrot在法兰西学院讲课时，首次提出了分维和分形几何的设想。分形（Fractal）一词，是Mandelbort教授创造出来的，其原意具有不规则、支离破碎等意义，分形几何学是一门以非规则几何形态为研究对象的几何学。由于不规则现象在自然界是普遍存在的，因此分形几何又称为描述大自然的几何学。分形几何建立以后，很快就引起了许多学科的关注，这是由于它不仅在理论上，而且在实用上都具有重要价值。分形几何与传统几何相比有何特点(1) 从整体上看，分形几何图形是处处不规则的。例如，海岸线和山川形状，从远距离观察，其形状是极不规则的。(2) 在不同尺度上，图形的规则性又是相同的。上述的海岸线和山川形状，从近距离观察，其局部形状又和整体形态相似，它们从整体到局部，都是自相似的。当然，也有一些分形几何图形，它们并不完全是自相似的。其中一些是用来描述一般随即现象的，还有一些是用来描述混沌和非线性系统的。Fractal(分形)一词的由来据Mandelbort教授自己说，fractal一词是1975年夏天的一个寂静夜晚，他在冥思苦想之余偶翻他儿子的拉丁文字典时突然想到的。此词源于拉丁文形容词fractus，对应的拉丁文动词是frangere（破碎、产生无规碎片）。此外与英文的fraction（碎片、分数）及fragment（碎片）具有相同的词根。在70年代中期以前，Mandelbort一直使用英文fractional一词来表示他的分形思想。因此，取拉丁词之头，撷英文之尾的fractal，本意是不规则的、破碎的、分数的。Mandelbort是想用此词来描述自然界中传统欧几里德几何学所不能描述的一大类复杂无规的几何对象。例如，弯弯曲曲的海岸线、起伏不平的山脉，粗糙不堪的断面，变幻无常的浮云，九曲回肠的河流，纵横交错的血管，令人眼花僚乱的满天繁星等。它们的特点是，极不规则或极不光滑。直观而粗略地说，这些对象都是分形。分形的定义Mandelbort曾经为分形下过两个定义： (1) 满足下式条件的集合A，称为分形集。其中，Dim(A)为集合A的Hausdoff维数（或分维数），dim(A)为其拓扑维数。一般说来，Dim(A)不是整数，而是分数。(2) 部分与整体以某种形式相似的形，称为分形。然而，经过理论和应用的检验，人们发现这两个定义很难包括分形如此丰富的内容。实际上，对于什么是分形，到目前为止还不能给出一个确切的科学定义，正如生物学中对“生命”也没有严格明确的定义一样，人们通常是列出生命体的一系列特性来加以说明。对分形的定义也可同样的处理。a. 分形集都具有任意小尺度下的比例细节，或者说它具有精细的结构。b. 分形集不能用传统的几何语言来描述，它既不是满足某些条件的点的轨迹，也不是某些简单方程的解集。c. 分形集具有某种自相似形式，可能是近似的自相似或者统计的自相似。d. 一般，分形集的“分形维数”，严格大于它相应的拓扑维数。e. 在大多数令人感兴趣的情形下，分形集由非常简单的方法定义，可能以变换的迭代产生。2.2 分形几何观及其应用平面上决定一条直线或圆锥曲线只需数个条件。那么决定一片蕨叶需要多少条件？如果把蕨叶看成是由线段拼合而咸，那么确定这片蕨叶的条件数相当可现，然而当人们以分形的眼光来看这片蕨叶时，可以把它认为是一个简单的迭代函数系统的结果，而确定该系统所需的条件数相比之下要少得多这说明用特定的分形拟合蕨叶比用折线拟合蕨叶更为有效。分形观念的引入并非仅是一个描述手法上的改变，从根本上讲分形反映了自然界中某些规律性的东西，以植物为例，植物的生长是植物细胞按一定的遗传规律不断发育、分裂的过程，这种按规律分裂的过程可以近似地看做是递归、迭代过程，这与分形的产生极为相似。在此意义上，人们可以认为一种植物对应一个迭代函数系统，人们甚至可以通过改变该系统中的某些参数来模拟植物的变异过程。分形几何还被用于海岸线的描绘及海图制作、地震预报、图象编码理论、信号处理等领域，并在这些领域内取得了个人注目的成绩。作为多个学科的交叉，分形几何对以往欧氏几何不屑一顾（或说是无能为力）的“病态”曲线的全新解释是人类认识客体不断开拓的必然结果。当前，人们迫切需要一种能够更好地研究、描述各种复杂自然曲线的几何学：而分形几何恰好可以堪当此用。所以说，分形几何也就是自然几何，以分形或分形的组合的眼光来看待周围的物质世界就是自然几何观。 第三章 Koch雪花的绘制3.1 von Koch曲线简介自从有了函数曲线的连续与可微性质及其关系以后，是否存在一个处处连续而点点不可微的函数曲线成了研究的热门。首先解决这个问题的是大数学家Weierstrass，他于1872年设计了如下的一个函数：图 3.1 W(x) (a=0.5,b=3)Weierstrass证明了对某些a和b的值，该函数无处可微。1916年，Hardy证明了对满足上列条件的所有a和b的值，W(x)都是无处可微的。W(x)的缺点是极难绘图，顾不够直观。到1904，瑞典数学家von Koch设计了一条被称之为Koch曲线的图形，起设计步骤如下：设E为单位区间0，1，第一步，即n=1，以E的中间1/3线段为底，向上作一个等边三角形，然后去掉区间(1/3,2/3)，得一条4折线段的多边形E。第二步，既n=2，对E的四条折线段重复上述过程，得到一条16折线段多边形E。如图3.2a b 图 3.2 von Koch曲线 3.2 Koch雪花算法设计不难想象，如果改变生成元(如图 3.3 n=0)，便可导至另外的曲线，图3.3便是一例。利用这一生成元迭代的作下去，迭代步数记为n，可以得到图 3.3所示的结果。当n时，生成的则是著名的Koch雪花。假设生成元的顶点为p0、p1、p2、p3，其中p0=p3，这样构成的等边三角形则是闭合的。我们现在用P=p0 p1 p2 p3来记录迭代过程中产生的新节点。首先由p0、p1计算出p2、p3，p4=ptemp*。这里我们将ptemp分解为ptemp(1)和ptemp(2)。其中主要算法如下：取点p0=(0,0),点p1=(1,0)如果令P=p0,p1，则有p0=P(:,i);p1=P(:,i+1);下面由p0、p1计算其它点 p2=2*p0/3+p1/3;p3=p0/3+2*p1/3;ptemp=p3-p2;p4=ptemp(1)*cos(pi/3)-ptemp(2)*sin(pi/3);ptemp(1)*sin(pi/3)+ptemp(2)*cos(pi/3)+p2;由下面的M语句完成对过程中产生的新节点的保护工作。Ptemp(:,flag)=p0;Ptemp(:,flag+1)=p2;Ptemp(:,flag+2)=p4;Ptemp(:,flag+3)=p3;n=0 n=1 n=2 n=10图 3.3 Koch曲线如果改变变换矩阵中的旋转角度，将会产生另一张图形，见图3.4。这里我们将变换矩阵变为：图 3.4 改变变换矩阵后的Koch雪花图第四章 Frac_tree绘制任意选定一个二维平面上的初始点坐标(x,y)。假设我们可以生成一个在0，1区间上均匀分布的随机数，那么根据其取值的大小，可以按下面的公式生成一个新的坐标点(x,y)。令新生成的点(x,y)为初始点(x,y)，可以再生成一个新的点，我们可以多次重复这样的过程，这样就可以生成一族坐标点。假设我们想根据这样的方式产生N个点，则可以由下面的MATLAB语句编写一个函数。function x,y=frac_tree(x0,y0,v,N)x=x0;zeros(N-1,1);y=y0;zeros(N-1,1);for i=2:Nvv=v(i);if vv0.05, y(i)=0.5*y(i-1);elseif vv0.45,x(i)=0.42*(x(i-1)-y(i-1);y(i)=0.2+0.42*(x(i-1)+y(i-1);elseif vv N=10000；v=rand(N,1);x,y=frac_tree(0,0,v,N);h=plot(x(1:10000),y(1:10000),.) set(h,markersize,4)将markersize设置为5； set(h,markersize,5)markersize=4 markersize=5图4.2 分形树第五章 Mandelbort集的绘制5.1 Mandelbort集简介1980年，Mandelbort给世人提供了一幅无与伦比的杰作：Mandelbort集。其创作过程如下：令,z是复变量，是复常数。从上述的复反馈过程得到如下迭代公式 (5.1.1)式中=(0,0)，为计算机荧屏位于(p,q)位置的像素。于是上式成了给定N为一个正整数。当像素位于(p,q)且n=N时,仍然小于给定的一个值K,则在(p,q)位置着色为1,否则当nM,则选择颜色k，转至步骤4;如果k=K，则选择颜色0，转至步骤4;如果rM，且kK，转至步骤2;步骤4：对点()着颜色k，转至下一个点(步骤1)。利用上面的过程，灵活地选取画面窗口，就可以观察到Mandelbort集的精细结构。在下面的程序中，用v1，v2，来决定一个平面，用offset来移动这个平面。图如5.1(a)% 由v1、v2和offset来决定一个平面v1=1 0 0 0; v2=0 1 0 0; offset=0 0 0 0;% p、q的起点、终点，决定画图的范围pmin=-2.25; pmax=0.75;qmin=-1.50; qmax=1.50;% 设收敛上界为100M=100;% 假设显示器的分辨率是a*b点，则a、b便决定了图的精细度a=800; b=600; delta_p=(pmax-pmin)/a; delta_q=(qmax-qmin)/b;% 用(a+1)*(b+1)零矩阵来记录每个发散的点out=zeros(a+1,b+1);for i=1:a+1;for j=1:b+1;u=(pmin+(i-1)*delta_p).*v1+(qmin+(j-1)*delta_q).*v2+offset;n=0; k=u; c=u;% 矩阵的模M且由c定义的数列的前64项都M，则认为c(由i、j决定)收敛，n=64，n越小收敛越快while norm(k)M & n x=linspace(-1.5,1.5,300);y=x;X,Y=meshgrid(x,y); n_iter=30;c=-0.7; W=Julia(X,Y,c,n_iter); 这样得出的z矩阵中有很多元素为NaN。MATLAB在解决这样问题中的优势就是它能有效地处理计算过程中的NaN，所以用MATLAB语言来计算这样的问题是合适的。我们可以调用pcolor()函数用为彩色描述得出的z矩阵。如果我们将色调设置成prism(),则最终可以由下面的MATLAB语句绘制出Julia分形图，如图6.1(a)所示。 pcolor(X,Y,W),shading flat; axis(square);colormap prism(256);(a) n_iter=30 (b) n_iter=300图 6.1 不同映射次数下的Julia图选择300次映射，则可以绘制出如图6.1(b)所示的结果。 n_iter=300; W=Julia(X,Y,c,n_iter);pcolor(X,Y,W),shading flat; axis(square);colormap prism(256); max(max(abs(W-W1)ans = 0.2305仔细比较这两个图，则可以看出它们之间是有很大区别的。在两种影射次数下的结果中的矩阵之间也有较大的偏差30次映射是没有完全收敛的。选择不同的参数和着色方案，我们可以看出分别由下面的MATLAB命令绘制出Julia集，如图6.2(a)所示。 n_iter=300; c=0.481762-0.321657j; W=Julia(X,Y,c,n_iter); pcolor(X,Y,W1),shading flat; axis(square); colormap prism(256);给定c=-0.766227+0.096990j，选用summer色调。如图6.2 (b) colormap summer(256); (a) c=0.481762-0.321657j (b) c=-0.766227+0.096990j图6.2 参数变化时的Julia6.3 Julia集与Mandelbort集 一幅Julia集的图完全依赖于常数c，所以