动力系统一些分形图像和matlab程序.docx
《动力系统一些分形图像和matlab程序.docx》由会员分享，可在线阅读，更多相关《动力系统一些分形图像和matlab程序.docx（18页珍藏版）》请在冰豆网上搜索。
动力系统一些分形图像和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;
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）
第二题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=[];
[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分〕
PUTEREXPERIMENT2.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;
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
for
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分〕
PUTEREXPERIMENT4.3（P170）
〔1〕
（a）（b）
i的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.4i,应用上面的程序得到如下图：
第九题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）
升级会员 