分形几何中一些经典图形的Matlab画法
分形几何中一些经典图形的Matlab画法作者：日期：分形几何中一些经典图形的Matlab画法(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;%第i1次迭代时由各条线段产生的新四条线段的五点横、纵坐标存储在数组A、B中A,B=sub_koch1(al,b1,a2,b2);fori=1:nforj=1:length(A)/5;w=sub_koch2(A(1+5*(j-1)：5*j),B(1+5*(j1)：5*j);fork=1:4AA(5*4*(j-l)+5*(k-1)+1:5*4*(j-1)+5*(k-1)+5),BB(5*4*(j1)+5*(k1)+1:5*4*(j-1)+5*(k1)+5)=sub_koch1(w(k,1),w(k,2),w(k,3),w(k,4);endendA=AA;B=BB;endplot(A,B)ho1donaxisequal%由以(ax,ay),(bx,by)为端点的线段生成新的中间三点坐标并把这五点横、纵坐标依次分别存%储在数组A,B中functionA,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).A2);a1pha=atan(eycy)./(excx);if(ex-cx)<0a1pha=alpha+pi;enddx=cx+cos(alpha+pi/3)*L;dy=cy+sin(a1pha+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);bl2=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,b421；图1VonKoch曲线(2)Levy曲线程序levy.mfunctionlevy(n)%levy(16),n为lcvy曲线迭代次数%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:nforj=1:1ength(X)/3w=levy2(X(1+3*(j1):3*j),Y(1+3*(j-1):3*j);XX(3*2*(j-1)+1:3*2*(j1)+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);endX=XX;Y=YY;endplot(X,Y)holdonaxisequal%由以(x1,y1),(x2,y2)为端点的线段生成新的中间点坐标并把(x1,y1),(x2,y2)连同新点横、纵坐标依次分别存储在数组X,Y中functionX,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;%把由函数1evy1生成的三点横、纵坐标X,Y顺次划分为两组，分别对应两条折线段中每条线段两端点的坐标，并依次分别存储在2*4阶矩阵w中，w中第i(i=1,2)行数字代表第i条线段两端点的坐标functionw=1evy2(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;0.0060.4020-0.2-0,J.0 4-0200.204C,60 811214图2Levy曲线(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)holdonX,Y=tree1(x1,y1,x2,y2ab);holdonW=tree2(X,Y);w1=W(:,1:4);w2=W(:,5:8);%w为2人k*4维矩阵，存储第k次迭代产生的分枝两端点的坐标%w的第i(i=1,2,2")行数字对应第i个分枝两端点的坐标w=w1;w2;fork=1:nfori=1:2AkX,Y=tree1(w(i,1),w(i,2),w(i,3),w(i,4),a,b);W(i,:尸tree2(X,Y);endw1=W(:,1:4);w2=W(:,5:8);w=w1;w2;end%由每个分枝两端点坐标(x1,y1),(x2,y2)产生两新点的坐标(x3,y3),(x4,y4),画两分枝图形，并把(x2,y2)连同新点横、纵坐标分别存储在数组X,Y中functionX,Y=tree1(x1,y1,x2,y2,a,b)L=sqrt(x2-x1)A2+(y2-y1)人2);if(x2x1)=0a=pi/2;e1seif(x2-x1)v0a=pi+atan(y2y1)/(x2x1);elsea=atan(y2-y1)/(x2x1);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;%把由函数treel生成的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;252150.50-15-1-0500.511.6图3分形树(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.50000.50;M2=0.500.500.50;M3=0.500.2500.50.5;x=0;y=0;%r为0,1区间内产生的n维随机数组r=rand(1,n);B=zeros(2,n);k=1;%当0vr(i)v1/3时,进彳tM1对应的压缩映射；%当1/3=<r(i)<2/3时，进彳TM2对应的压缩映射；%当2/3=vr(i)<1时，进彳tM3对应的压缩映射；fori=1:nifr(i)<w1a=M1(1);b=M1(2);e=M1(3);c=M1(4);d=M1(5);f=M1(6);y()，B(2KM公工4s.hia0ln=10000;cx=-0.77;cy=0.08;y)3 8)k=k+1; end% juli % f(els eSia =M2on ju 1 ia _ i f _ifs(100000,-0. 7s(n,cx, cx=a*a=M3(1);b = M3 (2);e=M3( 3 ) ;c=M3(4);d=M3(5); f = M 3 (6) en dendendfun c t iel s e i f r(i)<w1+w2plot( B1(5 ) IFS算法画J uia集程序juli az )=z A 2 + c, c x=re ax +b* y + e ; c* x + d *y+f;0.)B(1,k)=x;B ( 2,k) = y! (1 );b=M 2 (2); e =M2( 3 ); c =M 2 (4); d = M2(5);f=M2(6 ); i f r(i) <w1+ w2+w3mark ersize' ,%zA2+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:nwx=x-cx;wy=y-cy;ifwx>0a1pha=atan(wy/wx);endifwx<0a1pha=pi+atan(wy/wx);endifwx=0aIpha=pi/2;enda1pha=aIpha/2;r=sqrt(wx人2+wyA2);ifA(i)<0r=sqrt(r);elser=sqrt(r);endx=r*cos(alpha);y=r*sin(alpha);B(1,k)=x;B(2,k)=y;k=k+1;endpIot(B(1,:),B(2,:),'.',/markersize',0.1)图5Julia集(6)逃逸时间算法画Sierpinski垫片程序sierpinski.hfunctionsierpinski(a,b,c,d,n,m,r)%sierpinski(0,0,1,1,12,200,200)%(a,b),(c,d)收敛区域左上角和右下角坐标
