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

动力系统中一些很好看的图形matlab实现

研究生课程考核试卷

（适用于课程论文、提交报告）

科 目： 动力系统 教 师： 舒永录 姓 名： 郑申海 学 号： 20110602024 专 业： 计算数学 类 别： 学术 上课时间： 2012 年 3 月至2012 年 6 月

考

生 成 绩：

阅卷评语：

阅卷教师 (签名

重庆大学研究生院制

动力系统中一些很好看的图形matlab实现

第一题 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 iterate=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分）

动力系统中一些很好看的图形matlab实现

Figure2.3(P51)

(a)

（a）、（b）对应的初始值分别是a=1.28、1.4 绘图程序：

(b)

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

[X Y ]=meshgrid(M); plot(X,Y,'.k') hold on

[Ii Jj]=size(X); R=zeros(Ii,Jj); for i=1:Ii for 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; end

if xm<Iter&&ym<Iter R(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);

动力系统中一些很好看的图形matlab实现

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); hold on;box on; x=0; y=0;

A=0:0.0001:1.25; for a=A

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

动力系统中一些很好看的图形matlab实现

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)

动力系统中一些很好看的图形matlab实现

(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、f for a=A

for 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

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分）

COMPUTER EXPERIMENT 2.2 P(76)

动力系统中一些很好看的图形matlab实现

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

x=0;y=0;%对比试验用x=0.5,y=0.5 for n=2:N; xm=x; ym=y;

x=ym+1-a*xm.*xm; y=b*xm; xn(n)=x;

动力系统中一些很好看的图形matlab实现

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) axis equal; figure(2)

plot(A,resulty,'.','markersize',1) axis equal;

第六题 Mandelbrot 集合(20分)

Figure4.10(P168)

绘图程序：

clc, clear

ITER=50; N=200; hold on for a=-1.5:0.005:0.6

for 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*x(n)*y(n)+b; end

if x(end)^2+y(end)^2<ITER plot(a,b,'.') end

动力系统中一些很好看的图形matlab实现

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); hold on tic

for i=1:C for 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; end

动力系统中一些很好看的图形matlab实现

if x(end)^2+y(end)^2<ITER plot(X(i,j),Y(i,j),'.') end end end

axis([-1.5 1.5 -1.5 1.5]) t=toc

第八题 计算机实验（20分）

COMPUTER EXPERIMENT 4.3 (P 170) （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 %画网格 hold on

plot(Y(:,1),X,'g'); plot(X,Y(:,1)','g'); z=X+Y*i;

[C,R]=size(z); ATTxx=[]; ATTyy=[]; tic flg=1;

动力系统中一些很好看的图形matlab实现

for i=2:C for j=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; 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)]; end

if x(end)^2+y(end)^2<ITER plot(x(1),y(1),'.b') if flg==1

ATTxx=[ATTxx;ATTX]; ATTyy=[ATTyy;ATTY]; flg=0; end else

plot(x(1),y(1),'.k') end end end

axis([-1.3 1.3 -1.3 1.3]) t=toc

（2）六周期点研究

通过查找相关资料，找到Mandelbrot 集合与Julia 集合对应的周期关系： (我通过循环查找六周期点，好几天都没找到结果，最后放弃了

)

然后结合上图与第六题结论，取c=0.4-0.215i, 应用上面的程序得到下图：

动力系统中一些很好看的图形matlab实现

第九题 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];

动力系统中一些很好看的图形matlab实现

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;

% 当0<r(i)<1/3时,进行M1对应的压缩映射; % 当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); 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);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)

本文来源：https://www.bwwdw.com/article/8bs1.html