经典图形的Matlab.docx
《经典图形的Matlab.docx》由会员分享，可在线阅读，更多相关《经典图形的Matlab.docx（17页珍藏版）》请在冰豆网上搜索。
经典图形的Matlab
分形几何中一些经典图形的Matlab画法
来源：
杨昊的日志
（1）Koch曲线程序koch.m
functionkoch（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
forj=1:
length（A）/5;
w=sub_koch2（A（1+5*（j-1）:
5*j）,B（1+5*（j-1）:
5*j））;
fork=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）
holdon
axisequal
%由以（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条线段两端点的坐标
functionw=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];
（2）Levy曲线程序levy.m
functionlevy（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）;
fori=1:
n
forj=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）
holdon
axisequal
%由以（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条线段%两端点的坐标
functionw=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];
（3）分形树程序tree.h
functiontree（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]）
holdon
[X,Y]=tree1（x1,y1,x2,y2,a,b）;
holdon
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];
fork=1:
n
fori=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;
elseif（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）
axisequal
holdon
X=[x2,x3,x4];
Y=[y2,y3,y4];
%把由函数tree1生成的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];
（4）IFS算法画Sierpinski三角形程序sierpinski_ifs.h
functionsierpinski_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.50 0 00.5 0];
M2=[0.500.5 00.5 0];
M3=[0.500.2500.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=fori=1:
n
ifr（i） a=M1
（1）;b=M1
（2）;e=M1（3）;c=M1（4）;d=M1（5）;f=M1（6）;
elseif r（i） a=M2
（1）;b=M2
（2）;e=M2（3）;c=M2（4）;d=M2（5）;f=M2（6）;
elseif 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）
（5）IFS算法画Julia集程序julia_ifs.h
functionjulia_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）;
fori=1:
n
wx=x-cx;
wy=y-cy;
ifwx>0
alpha=atan（wy/wx）;
end
ifwx<0
alpha=pi+atan（wy/wx）;
end
ifwx==0
alpha=pi/2;
end
alpha=alpha/2;
r=sqrt（wx^2+wy^2）;
ifA（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）
（6） 逃逸时间算法画Sierpinski垫片程序sierpinski.h
functionsierpinski（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;
fori=1:
m
x0=a+（c-a）*（i-1）/m;
forj=1:
m
y0=b+（d-b）*（j-1）/m;
x=x0;
y=y0;
fork=1:
n
ify>0.5
x=2*x;
y=2*y-1;
elseifx>=0.5
x=2*x-1;
y=2*y;
else
x=2*x;
y=2*y;
end
ifx^2+y^2>r
break;
end
end
ifk==n
B（1,w）=i;
B（2,w）=j;
w=w+1;
end
end
end
plot（B（1,:
）,B（2,:
）,'.','markersize',0.1）
（7） 元胞自动机算法画Sierpinski三角形程序
一维元胞自动机sierpinski_ca1.h
functionsierpinski_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;
fori=1:
m-1
forj=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）
² 二维元胞自动机sierpinski_ca2.h
functionsierpinski_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;
fori=[m/2:
-1:
2,m/2:
m-1]
forj=[n/2:
-1:
2,n/2:
n-1]
ifmod（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）
（8） IFS算法画Helix曲线程序helix_ifs.h
functionhelix_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.7586470.242424 0.859848 1.408065];
M2=[-0.1212120.257576 -6.7216540.05303 0.05303 1.377236];
M3=[0.181818 -0.136364 6.0861070.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=fori=1:
n
ifr（i） a=M1
（1）;b=M1
（2）;e=M1（3）;c=M1（4）;d=M1（5）;f=M1（6）;
elseif r（i） a=M2
（1）;b=M2
（2）;e=M2（3）;c=M2（4）;d=M2（5）;f=M2（6）;
elseif 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）
升级会员 