动力系统一些分形图像和matlab程序.docx

1、动力系统一些分形图像和matlab程序研究生课程考核试卷适用于课程论文、提交报告科 目： 动力系统 教 师： 舒永录 姓 名： X申海 学 号： 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;for a=A; x=0.5;for

2、 iterate=2:iterate_max x(iterate)=ga(a,x(iterate-1);endresult=result;x(iterate_max-plot_N+1):iterate_max);endplot(A,result,-)Figure1.7(P20)第二题 Henon映射初始条件10分Figure2.3(P51)(a)(b)绘图程序：a=1.4;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=;X Y =meshgrid(M);plot(X,Y,.k)hol

3、d onIi Jj=size(X);R=zeros(Ii,Jj);for i=1:Iifor j=1:Jj xm=X(i,j); ym=Y(i,j);for n=1:N x=a-xm.*xm+b*ym; y=xm; xm=x; ym=y;endif xmIter&ymIter R(i,j)=1; M_f=M_f,M(j); H_f=H_f,H(i);endendendm=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

4、(1,N);xn=zeros(1,N);yn=zeros(1,N);hold on;box on;x=0;y=0;A=0:0.0001:1.25;for a=Afor k=1:N; xm=x; ym=y; x=ym+1-a*xm.*xm; y=b*xm;end xn(1)=x;for n=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);endplot(A,result,.,markersize,1)xlim(0,a);第四题 Henon映射吸引子15分F

5、igure2.17P75(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);hold on;box on;x=0;y=0;A=0.9;%分别调节得到图a、b、c、d、e、ffor a=Afor k=1:N; xm=x; ym=y; x=ym+1-a*xm.*xm; y=b*xm;endxn(1)=x;for n=2:N; xm=x; ym=y; x=ym+1-a*xm.*xm; y=b*xm; xn(n)=x; yn(n)=y;

6、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.5 2 -0.8 0.8)第五题 计算机实验10分PUTER EXPERIMENT 2.2 P(76)(a) (b)(c)(d)如图是b=-0.3的henon映射分形图：(a)、(b)是初始值为(0,0)时对应x、y坐标的分析图，(c)、(d)是初始值为

7、(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);hold on;box on;A=0:0.0001:2.2;for a=A x=0;y=0;for n=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)

8、:N); resulty=resulty;yn(N-plot_N+1):N);endfigure(1)plot(A,resultx,.,markersize,1)axis equal;figure(2)plot(A,resulty,.,markersize,1)axis equal;第六题 Mandelbrot 集合(20分)Figure4.10(P168)绘图程序：clc, clear ITER=50; N=200; hold onforfor b=-1.0:0.005:1.0 x(1)=0; y(1)=0; for n=1:N x(n+1)=x(n)2-y(n)2+a; y(n+1)=2*

9、x(n)*y(n)+b; endif x(end)2+y(end)2ITER plot(a,b,.) endendend第七题 Julia 集合(20分)Figure4.11P169(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);%此三行用于画图bz=X+Y*i;C,R=size(

10、z);hold onticfor i=1:Cfor j=1:R x(1)=X(i,j); y(1)=Y(i,j); for n=1:N x(n+1)=x(n)2-y(n)2+a; y(n+1)=2*x(n)*y(n)+b; endif x(end)2+y(end)2ITER plot(X(i,j),Y(i,j),.) endendendaxis(-1.5 1.5 -1.5 1.5)t=toc第八题 计算机实验20分PUTER EXPERIMENT 4.3 (P 170)1(a)(b)i的Julia集合，其中蓝色的点是收敛的。图ab分别对应N=100、700的情况。蓝色和黑色的交界就是Julia

11、集合。绘图程序：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%画网格hold onplot(Y(:,1),X,g);plot(X,Y(:,1),g);z=X+Y*i;C,R=size(z);ATTxx=;ATTyy=;ticflg=1;for i=2:Cfor j=2:R ATTX=zeros(1,10); ATTY=zeros(1,10); x(1)=(X(i,j)+X(i,j-1)/2; y

12、(1)=(Y(i,j)+Y(i-1,j)/2; for n=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);endif x(end)2+y(end)2ITERplot(x(1),y(1),.b) if flg=1 ATTxx=ATTxx;ATTX; ATTyy=ATTyy;ATTY; flg=0;endelse plot(x(1),y(1),.k) endendendaxis(-1.3 1.3 -1.3 1.3)t=toc2六周期点研究通过查找相关资

13、料，找到Mandelbrot 集合与Julia 集合对应的周期关系：(我通过循环查找六周期点，好几天都没找到结果，最后放弃了)然后结合上图与第六题结论，取c=0.4i, 应用上面的程序得到如下图：第九题 Sierpinski三角形(20分)EXAMPLE 4.7(P159)绘图程序：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=rand(1,n);B=zeros(2,n);k=1;% 当0r(i)1/3时,进展M1对应的压缩映射;% 当

14、1/3=r(i)2/3时,进展M2对应的压缩映射;% 当2/3=r(i)1时,进展M3对应的压缩映射;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);elseif r(i)w1+w2 a=M2(1);b=M2(2);e=M2(3);c=M2(4);d=M2(5);f=M2(6); elseif r(i)w1+w2+w3 a=M3(1);b=M3(2);e=M3(3);c=M3(4);d=M3(5);f=M3(6);endendendx=a*x+b*y+e;y=c*x+d*y+f;B(1,k)=x;B(2,k)=y;k=k+1; endplot(B(1,:),B(2,:),.,markersize,0.1)