分形幾何中一些經(jīng)典圖形的Matlab畫法(共14頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上分形幾何中一些經(jīng)典圖形的Matlab畫法（1）Koch曲線程序koch.mfunction koch(a1,b1,a2,b2,n)%koch(0,0,9,0,3)%a1,b1,a2,b2為初始線段兩端點坐標,n為迭代次數(shù)a1=0;b1=0;a2=9;b2=0;n=3;%第i-1次迭代時由各條線段產(chǎn)生的新四條線段的五點橫、縱坐標存儲在數(shù)組A、B中A,B=sub_koch1(a1,b1,a2,b2);for i=1:n for j=1:length(A)/5; w=sub_koch2(A(1+5*(j-1):5*j),B(1+5*(j-1):5*j); for k=1:4

2、 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)hold onaxis equal%由以(ax,ay),(bx,by)為端點的線段生成新的中間三點坐標并把這五點橫、縱坐標依次分別存%儲在數(shù)組A,B中function A,B=sub_koch1(ax,ay,bx,by)cx=ax+(bx-ax)/3;cy=ay+(by-ay)/3

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;%把由函數(shù)sub_koch1生成的五點橫、縱坐標A,B順次劃分為四組,分別對應(yīng)四條折線段中%每條線段兩端點的坐標,并依次分別存儲在4*4階矩陣k中,k中第i(i=1,2,3,4)行數(shù)字代表第%i

4、條線段兩端點的坐標function w=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;圖1 Von Koch曲線（2）Levy 曲線程序levy.mfunction levy(n)% levy(16

5、),n為levy曲線迭代次數(shù)%x1,y1,x2,y2為初始線段兩端點坐標,nn為迭代次數(shù)n=16;x1=0;y1=0;x2=1;y2=0;%第i-1次迭代時由各條線段產(chǎn)生的新兩條線段的三端點橫、縱坐標存儲在數(shù)組X、Y中X,Y=levy1(x1,y1,x2,y2);for i=1:n for j=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

6、*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)hold onaxis equal%由以(x1,y1),(x2,y2)為端點的線段生成新的中間點坐標并把(x1,y1),(x2,y2)連同新點橫、縱坐%標依次分別存儲在數(shù)組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=y

7、1,y3,y2;%把由函數(shù)levy1生成的三點橫、縱坐標X,Y順次劃分為兩組,分別對應(yīng)兩條折線段中每條線%段兩端點的坐標,并依次分別存儲在2*4階矩陣w中,w中第i(i=1,2)行數(shù)字代表第i條線段%兩端點的坐標function w=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;圖2 Levy 曲線（3） 分形樹程序tree.hfunction tree(n,a,b)% tree(8,pi/8,pi/8),n為

8、分形樹迭代次數(shù)%a,b為分枝與豎直方向夾角%x1,y1,x2,y2為初始線段兩端點坐標,nn為迭代次數(shù)n=8;a=pi/8;b=pi/8;x1=0;y1=0;x2=0;y2=1;plot(x1,x2,y1,y2)hold onX,Y=tree1(x1,y1,x2,y2,a,b);hold onW=tree2(X,Y);w1=W(:,1:4);w2=W(:,5:8);% w為2k*4維矩陣,存儲第k次迭代產(chǎn)生的分枝兩端點的坐標,% w的第i(i=1,2,2k)行數(shù)字對應(yīng)第i個分枝兩端點的坐標w=w1;w2;for k=1:n for i=1:2k X,Y=tree1(w(i,1),w(i,2),

9、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)產(chǎn)生兩新點的坐標(x3,y3),(x4,y4),畫兩分枝圖形,并把%(x2,y2)連同新點橫、縱坐標分別存儲在數(shù)組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; else if (x2-x1)<0 a=pi+atan(y2-y1)/(x2-x1); else

10、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)axis equalhold onX=x2,x3,x4;Y=y2,y3,y4;%把由函數(shù)tree1生成的X,Y順次劃分為兩組,分別對應(yīng)兩分枝兩個端點的坐標,并存儲在一維%數(shù)組w中function w=tree2(X,Y)a1=X(1);b1=Y(1);a2=X(2);b2=Y(2);a3=X(1);b3=Y(1

11、);a4=X(3);b4=Y(3);w=a1,b1,a2,b2,a3,b3,a4,b4;圖3 分形樹（4）IFS算法畫Sierpinski三角形程序sierpinski_ifs.hfunction sierpinski_ifs(n,w1,w2,w3)%sierpinski_ifs(10000,1/3,1/3,1/3)%w1,w2,w3出現(xiàn)頻率n=10000;w1=1/3;w2=1/3;w3=1/3;M1=0.5 0 0 0 0.5 0;M2=0.5 0 0.5 0 0.5 0;M3=0.5 0 0.25 0 0.5 0.5;x=0;y=0;% r為0,1區(qū)間內(nèi)產(chǎn)生的n維隨機數(shù)組r=rand(1

12、,n);B=zeros(2,n);k=1;% 當0<r(i)<1/3時,進行M1對應(yīng)的壓縮映射;% 當1/3=<r(i)<2/3時,進行M2對應(yīng)的壓縮映射;% 當2/3=<r(i)<1時,進行M3對應(yīng)的壓縮映射;for i=1:n if r(i)<w1 a=M1(1);b=M1(2);e=M1(3);c=M1(4);d=M1(5);f=M1(6); else if r(i)<w1+w2 a=M2(1);b=M2(2);e=M2(3);c=M2(4);d=M2(5);f=M2(6); else if r(i)<w1+w2+w3 a=M3(1)

13、;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)圖4 Sierpinski三角形（5）IFS算法畫Julia集程序julia_ifs.hfunction julia_ifs(n,cx,cy)% julia_ifs(,-0.77,0.08)% f(z)=z2+c,cx=real(c);cy=image(c);n=10000;

14、cx=-0.77;cy=0.08;% z2+c=z0,x=real(z0);y=image(z0);x=1;y=1; B=zeros(2,n);k=1;% A為產(chǎn)生的服從標準正態(tài)分布的n維隨機數(shù)組A=randn(1,n);for i=1:n wx=x-cx; wy=y-cy; if wx>0 alpha=atan(wy/wx); end if wx<0 alpha=pi+atan(wy/wx); end if wx=0 alpha=pi/2; end alpha=alpha/2; r=sqrt(wx2+wy2); if A(i)<0 r=-sqrt(r); else r=sq

15、rt(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)圖5 Julia 集(6) 逃逸時間算法畫Sierpinski墊片程序sierpinski.hfunction sierpinski(a,b,c,d,n,m,r)%sierpinski(0,0,1,1,12,200,200)%(a,b),(c,d)收斂區(qū)域左上角和右下角坐標,m為分辨率% n為逃逸時間,需要反復(fù)試探,r逃逸半徑a=0;b=0

16、;c=1;d=1;n=12;m=200;r=200;B=zeros(2,m*m);w=1;for i=1:m x0=a+(c-a)*(i-1)/m; for j=1:m y0=b+(d-b)*(j-1)/m; x=x0; y=y0; for k=1:n if y>0.5 x=2*x; y=2*y-1; else if x>=0.5 x=2*x-1; y=2*y; else x=2*x; y=2*y; end if x2+y2>r break; end end if k=n B(1,w)=i; B(2,w)=j; w=w+1; end endendplot(B(1,:),B(2

17、,:),'.','markersize',0.1)圖6 Sierpinski三角形墊片(7) 元胞自動機算法畫Sierpinski三角形程序² 一維元胞自動機sierpinski_ca1.hfunction sierpinski_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;for i=1:m-1 for j=2:n-1 if (s(i,j-1)=1&s(i,j)=0&s(i,j

18、+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) 圖7.1 一維元胞自動機畫Sierpinski三角形² 二維元胞自動機sierpinski_ca2.hfunction sierpinski_ca2(m,n)%sierpinski_ca2(400,400)m=400;n=400;t=1;w=zeros(2,m*n)

19、;s=zeros(m,n);s(m/2,n/2)=1;for i=m/2:-1:2,m/2:m-1 for j=n/2:-1:2,n/2:n-1 if mod(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)圖7.2 二維元胞自動機畫Sierpinski三角形(8) IFS算法畫Helix曲線程序helix_ifs.hfunction helix_ifs(n,w1,w2,w3)%helix_ifs(20000,0.9,0.05,0.