动力系统一些分形图像和matlab程序.docx
《动力系统一些分形图像和matlab程序.docx》由会员分享,可在线阅读,更多相关《动力系统一些分形图像和matlab程序.docx(20页珍藏版)》请在冰豆网上搜索。
动力系统一些分形图像和matlab程序
研究生课程考核试卷
(适用于课程论文、提交报告)
科目:
动力系统教师:
舒永录
姓名:
郑申海学号:
20110602024
专业:
计算数学类别:
学术
上课时间:
2012年3月至2012年6月
考生成绩:
卷面成绩
平时成绩
课程综合成绩
阅卷评语:
阅卷教师(签名)
重庆大学研究生院制
第一题Logistic映射(15分)
Figure1.6(P19)
绘图程序:
ga=inline('a*x*(1-x)');
plot_N=100;
iterate_max=200;
result=[];
A=1:
0.001:
4;
fora=A;
x=0.5;
foriterate=2:
iterate_max
x(iterate)=ga(a,x(iterate-1));
end
result=[result;x((iterate_max-plot_N+1):
iterate_max)];
end
plot(A,result,'-')
Figure1.7(P20)
注:
这两个图是Figure1.6的局部放大图
第二题Henon映射初始条件(10分)
Figure2.3(P51)
(a)(b)
(a)、(b)对应的初始值分别是a=1.28、1.4
绘图程序:
a=1.4;%a=1.28
b=-0.3;
N=200;
Iter=3;
M=linspace(-2.5,2.5,500);
M_f=[];
H=linspace(-2.5,2.5,500);
H_f=[];
[XY]=meshgrid(M);
plot(X,Y,'.k')
holdon
[IiJj]=size(X);
R=zeros(Ii,Jj);
fori=1:
Ii
forj=1:
Jj
xm=X(i,j);
ym=Y(i,j);
forn=1:
N
x=a-xm.*xm+b*ym;
y=xm;
xm=x;
ym=y;
end
ifxmR(i,j)=1;
M_f=[M_f,M(j)];
H_f=[H_f,H(i)];
end
end
end
m=size(M_f);
h=size(H_f);
plot(M_f,H_f,'.w')
第三题Henon映射分叉图(15分)
Figure2.16(P74)
绘图程序:
b=0.4;
N=200;
plot_N=150;
result=[];
an=ones(1,N);
xn=zeros(1,N);
yn=zeros(1,N);
holdon;boxon;
x=0;
y=0;
A=0:
0.0001:
1.25;
fora=A
fork=1:
N;
xm=x;
ym=y;
x=ym+1-a*xm.*xm;
y=b*xm;
end
xn
(1)=x;
forn=2:
N;
xm=x;
ym=y;
x=ym+1-a*xm.*xm;
y=b*xm;
xn(n)=x;
yn(n)=y;
end
result=[result;xn((N-plot_N+1):
N)];
end
plot(A,result,'.','markersize',1)
xlim([0,a]);
第四题Henon映射吸引子(15分)
Figure2.17(P75)
(a)(b)
(c)(d)
(e)(f)
绘图程序:
b=0.4;
N=2000;
plot_N=1500;
resultx=[];
resulty=[];
an=ones(1,N);
xn=zeros(1,N);
yn=zeros(1,N);
holdon;boxon;
x=0;
y=0;
A=0.9;%分别调节得到图a、b、c、d、e、f
fora=A
fork=1:
N;
xm=x;
ym=y;
x=ym+1-a*xm.*xm;
y=b*xm;
end
xn
(1)=x;
forn=2:
N;
xm=x;
ym=y;
x=ym+1-a*xm.*xm;
y=b*xm;
xn(n)=x;
yn(n)=y;
end
resultx=[resultx;xn((N-plot_N+1):
N)]
resulty=[resulty;yn((N-plot_N+1):
N)]
end
plot(resultx,resulty,'+','markersize',8)%a、b、d用
%plot(resultx,resulty,'.','markersize',3)%c、e、f用此
axis([-1.52-0.80.8])
第五题计算机实验(10分)
COMPUTEREXPERIMENT2.2P(76)
(a)(b)
(c)(d)
如图是b=-0.3的henon映射分形图:
(a)、(b)是初始值为(0,0)时对应x、y坐标的分析图,(c)、(d)是初始值为(0.5,0.5)时对应x、y坐标的分析图。
从图中可以看出,henon分形图与初始值有关,x、y坐标对应的分形图周期相同,但是轨迹不同。
绘图程序:
b=-0.3;
N=200;
plot_N=100;
resultx=[];
resulty=[];
an=ones(1,N);
xn=zeros(1,N);
yn=zeros(1,N);
holdon;boxon;
A=0:
0.0001:
2.2;
fora=A
x=0;y=0;%对比试验用x=0.5,y=0.5
forn=2:
N;
xm=x;
ym=y;
x=ym+1-a*xm.*xm;
y=b*xm;
xn(n)=x;
yn(n)=y;
end
resultx=[resultx;xn((N-plot_N+1):
N)];
resulty=[resulty;yn((N-plot_N+1):
N)];
end
figure
(1)
plot(A,resultx,'.','markersize',1)
axisequal;
figure
(2)
plot(A,resulty,'.','markersize',1)
axisequal;
第六题Mandelbrot集合(20分)
Figure4.10(P168)
绘图程序:
clc,clear
ITER=50;N=200;holdon
fora=-1.5:
0.005:
0.6
forb=-1.0:
0.005:
1.0
x
(1)=0;y
(1)=0;
forn=1:
N
x(n+1)=x(n)^2-y(n)^2+a;y(n+1)=2*x(n)*y(n)+b;
end
ifx(end)^2+y(end)^2plot(a,b,'.')
end
end
end
第七题Julia集合(20分)
Figure4.11(P169)
(a)(b)
(b)(d)
绘图程序:
a=0.29;
b=0.54;%调节参数a、b得到相应的图
c=a+b*i;
N=100;
ITER=100;
[X,Y]=meshgrid(linspace(-1.5,1.5,700));
%xx=linspace(-0.19,0.01,500);
%yy=linspace(0.89,1.09,500);
%[X,Y]=meshgrid(xx,yy);%此三行用于画图b
z=X+Y*i;
[C,R]=size(z);
holdon
tic
fori=1:
C
forj=1:
R
x
(1)=X(i,j);y
(1)=Y(i,j);
forn=1:
N
x(n+1)=x(n)^2-y(n)^2+a;y(n+1)=2*x(n)*y(n)+b;
end
ifx(end)^2+y(end)^2plot(X(i,j),Y(i,j),'.')
end
end
end
axis([-1.51.5-1.51.5])
t=toc
第八题计算机实验(20分)
COMPUTEREXPERIMENT4.3(P170)
(1)
(a)(b)
图为c=0.29+0.54i的Julia集合,其中蓝色的点是收敛的。
图(a)(b)分别对应N=100、700的情况。
蓝色和黑色的交界就是Julia集合。
绘图程序:
a=0.29;
b=0.54;
c=a+b*i;
N=100;
ITER=50;
[X,Y]=meshgrid(linspace(-1.3,1.3,100));
%[X,Y]=meshgrid(linspace(-1.3,1.3,700));%700的时候要运行490s
%画网格
holdon
plot(Y(:
1),X,'g');
plot(X,Y(:
1)','g');
z=X+Y*i;
[C,R]=size(z);
ATTxx=[];
ATTyy=[];
tic
flg=1;
fori=2:
C
forj=2:
R
ATTX=zeros(1,10);
ATTY=zeros(1,10);
x
(1)=(X(i,j)+X(i,j-1))/2;y
(1)=(Y(i,j)+Y(i-1,j))/2;
forn=1:
N
x(n+1)=x(n)^2-y(n)^2+a;y(n+1)=2*x(n)*y(n)+b;
ATTX=[ATTX(2:
end),x(n+1)];
ATTY=[ATTY(2:
end),y(n+1)];
end
ifx(end)^2+y(end)^2plot(x
(1),y
(1),'.b')
ifflg==1
ATTxx=[ATTxx;ATTX];
ATTyy=[ATTyy;ATTY];
flg=0;
end
else
plot(x
(1),y
(1),'.k')
end
end
end
axis([-1.31.3-1.31.3])
t=toc
(2)六周期点研究
通过查找相关资料,找到Mandelbrot集合与Julia集合对应的周期关系:
(我通过循环查找六周期点,好几天都没找到结果,最后放弃了)
然后结合上图与第六题结论,取c=0.4-0.215i,应用上面的程序得到下图:
第九题Sierpinski三角形(20分)
EXAMPLE4.7(P159)
绘图程序:
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=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);
elseifr(i)a=M2
(1);b=M2
(2);e=M2(3);c=M2(4);d=M2(5);f=M2(6);
elseifr(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)