分形几何中一些经典图形的Matlab画法+[文档在线提供]

分形几何中一些经典图形的Matlab画法

（1）Koch曲线程序koch.m

function koch(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

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

[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;

end

plot(A,B)

hold on

axis equal

%由以(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;

end

dx=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条线段两端点的坐标

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 V on Koch曲线

（2）Levy 曲线程序levy.m

function levy(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);

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*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;

end

plot(X,Y)

hold on

axis equal

%由以(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条线段%两端点的坐标

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.h

function tree(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])

hold on

[X,Y]=tree1(x1,y1,x2,y2,a,b);

hold on

W=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];

for k=1:n

for i=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;

else if (x2-x1)<0

a=pi+atan((y2-y1)/(x2-x1));

else

a=atan((y2-y1)/(x2-x1));

end

end

x3=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 equal

hold on

X=[x2,x3,x4];

Y=[y2,y3,y4];

%把由函数tree1生成的X,Y顺次划分为两组,分别对应两分枝两个端点的坐标,并存储在一维%数组w中

function w=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];

图3 分形树

（4）IFS算法画Sierpinski三角形程序sierpinski_ifs.h

function sierpinski_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.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]区间内产生的n维随机数组

r=rand(1,n);

B=zeros(2,n);

k=1;

% 当0

% 当1/3=

% 当2/3=

for i=1:n

if r(i)

a=M1(1);b=M1(2);e=M1(3);c=M1(4);d=M1(5);f=M1(6);

else if r(i)

a=M2(1);b=M2(2);e=M2(3);c=M2(4);d=M2(5);f=M2(6);

else if r(i)

a=M3(1);b=M3(2);e=M3(3);c=M3(4);d=M3(5);f=M3(6);

end

end

end

x=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.h

function julia_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);

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(wx^2+wy^2);

if A(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;

end

plot(B(1,:),B(2,:),'.','markersize',0.1)

图5 Julia 集(6)逃逸时间算法画Sierpinski垫片程序sierpinski.h

function sierpinski(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;

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 x^2+y^2>r

break;

end

end

if k==n

B(1,w)=i;

B(2,w)=j;

w=w+1;

end

end

end

plot(B(1,:),B(2,:),'.','markersize',0.1)

图6 Sierpinski三角形垫片(7)元胞自动机算法画Sierpinski三角形程序

一维元胞自动机sierpinski_ca1.h

function 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+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

end

end

plot(w(1,:),w(2,:),'.','markersize',1)

图7.1 一维元胞自动机画Sierpinski三角形

二维元胞自动机sierpinski_ca2.h

function sierpinski_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;

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

end

end

plot(w(1,:),w(2,:),'.','markersize',0.1)

图7.2 二维元胞自动机画Sierpinski三角形

(8)IFS算法画Helix曲线程序helix_ifs.h

function helix_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.758647 0.242424 0.859848 1.408065];

M2=[-0.121212 0.257576 -6.721654 0.05303 0.05303 1.377236];

M3=[0.181818 -0.136364 6.086107 0.090909 0.181818 1.568035];

x=0;y=0;

% r为[0,1]区间内产生的n维随机数组

r=rand(1,n);

B=zeros(2,n);

k=1;

% 当0

% 当1/3=

% 当2/3=

for i=1:n

if r(i)

a=M1(1);b=M1(2);e=M1(3);c=M1(4);d=M1(5);f=M1(6);

else if r(i)

a=M2(1);b=M2(2);e=M2(3);c=M2(4);d=M2(5);f=M2(6);

else if r(i)

a=M3(1);b=M3(2);e=M3(3);c=M3(4);d=M3(5);f=M3(6);

end

end

end

x=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)

图8 Helix曲线