|
ly1=[];ly2=[];C=[];x=0;y=0;w=eye(2,2);
for a=0:0.001:1.4
b=0.3;
N=1000; % NUMBER OF ITERATIONS
sl1=0; sl2=0;
for i=1:N
xprev=x;
yprev=y;
x=a-xprev.*xprev+b*yprev;
y=xprev;
jac=[-2*x b; 1 0]; %%% JACOBIAN OF THE HENON MAP
F=jac*w;
[w,r]=qr(F);
sl1 = sl1 + log(abs(diag(r)));
l1=sl1/N;
end
if N==1000
ly1=[ly1;l1(1)];ly2=[ly2;l1(2)];
C=[C;a];
%fprintf(1,'l1=%f\n',ly1);
end
end
subplot(2,1,1)
hold on
a=0:0.001:1.4;
line(a,0,'Color','k','LineWidth',4)
plot(C,ly1,'k',C,ly2,'k')
grid
%画Henon系统的分岔图程序!
b=0.3;
n=2000;
x=0;
y=0;
a=0:0.001:1.4;
for i=1:n/100
xprev=x;
yprev=y;
x=a-xprev.*xprev+b*yprev;
y=xprev;
end
for i=1:n
xprev=x;
yprev=y;
x=a-xprev.*xprev+b*yprev;
y=xprev;
if i>1990
hold on
subplot(2,1,2)
plot(a,x,'k.','markersize',1);
title('Henon Bifurcation');
end
end
grid
%%%%%%%求取最大LE的程序
clear all
d0=1e-12;C=[];Le=[];
for i=1:639
c=i/320;
x1=0;y1=0;
x2=0;y2=d0;
lsum=0;
for j=1:500
x3=1-c*x1*x1+0.2*y1;
y1=x1;
x1=x3;
x4=1-c*x2*x2+0.2*y2;
y2=x2;
x2=x4;
d1=sqrt((x2-x1)^2+(y2-y1)^2);
x2=x1+(d0/d1)*(x2-x1);
y2=y1+(d0/d1)*(y2-y1);
if j>100
lsum=lsum+log(d1/d0);
end
end
le=lsum/(j-100);
C=[C;c];Le=[Le;le];
end
plot(C,Le,'k')
|
|