动力系统一些分形图像和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;

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

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

[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

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

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;

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

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)

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

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

plot(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);

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

if x(end)^2+y(end)^2

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;

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

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, 应用上面的程序得到下图：

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

% 当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)