毕业了,本博客今后不再进行维护!勿发邮件,请谅解。

 

Friday, December 29, 2006

基于奇异值分解的Lyapunov指数计算


前几星期一个学友和我探讨了《一类基于奇异值分解的Lyapunov指数计算方法》一文
的方法,他把Henon映射的程序做出了,但Lorenz的程序却得不到正确的值。我试了
并得到了比较好的结果。现在贴出和大家分享。

henon_le.m




% 奇异值分解求Lyapunov法
% 差分henon系统

x=0.6;y=0.4;
V=eye(2);
S=V;b=0;
k=200;
for i=1:k
J=[-2.8*x 0.3;1 0];
B=J*V*S;
[V,S,U]=svd(B);
a_max=max(diag(S));
S=(1/a_max)*S;
b=b+log(a_max);
x_next=1+0.3*y-1.4*x*x;
y_next=x;
x=x_next;
y=y_next;
end
Lyapunov=(log(diag(S))+b)/k



Lyapunov =

0.4120
-1.6160


lorenz_euler.m



% Lorenz图形(欧拉方法)
%
% Author:yujunjie
% Author's email: ustb03-07@yahoo.com.cn
%

clear
h=0.006;a=16;b=4;c=49.52;
x=20;y=20;z=50;
Y=[];
for i=1:8000
x1=x+h*a*(y-x);
y1=y+h*(c*x-x*z-y);
z1=z+h*(x*y-b*z);
x=x1;y=y1;z=z1;
Y(i,:)=[x y z];
end
plot3(Y(:,1),Y(:,2),Y(:,3));





lorenz_le_eu.m



% 奇异值分解求Lyapunov法
% 微分lorenz系统
%
% Author:yujunjie
% Author's email: ustb03-07@yahoo.com.cn
%

x=20;y=20;z=50;h=0.002;
a=16;b=4;c=49.52;
V=eye(3);
S=V;b1=0;
k=4000;
for i=1:k
x1=x+h*a*(y-x);
y1=y+h*(c*x-x*z-y);
z1=z+h*(x*y-b*z);
x=x1;y=y1;z=z1;
J=[-a a 0
c-z -1 -x
y x -b];
J=eye(3)+h*J;
B=J*V*S;
[V,S,U]=svd(B);
a_max=max(diag(S));
S=(1/a_max)*S;
b1=b1+log(a_max);
end
Lyapunov=(log(diag(S))+b1)/(k*h)



Lyapunov =

1.5172
0.0016
-22.5679


lorenz_le_rk.m



% 奇异值分解求Lyapunov法
% 微分lorenz系统
%
% dx/dt = SIGMA*(y - x)
% dy/dt = BETA*x - y -x*z
% dz/dt= x*y - GAMA*z
% In this demo, SIGMA = 16,BETA = 45.92, GAMA = 4
% Initial conditions: x(0) = 1, y(0) = 1, z(0) = 1;
% Reference values:
% LE1 = 1.497, LE2 = 0.00, LE3 = -22.46, LD = 2.07

sigma=16;beta=45.92;gama=4;
x(1)=1;y(1)=1;z(1)=1;V=diag(ones(1,3));S=V;I=V;b=0;
h=0.005;k=5000;
for i=1:k
k11=sigma*(-x(i)+y(i));
k21=beta*x(i)-x(i)*z(i)-y(i);
k31=x(i)*y(i)-gama*z(i);

k12=sigma*(-(x(i)+0.5*h*k11)+(y(i)+0.5*h*k21));
k22=beta*(x(i)+0.5*h*k11)-(x(i)+0.5*h*k11)...
*(z(i)+0.5*h*k31)-(y(i)+0.5*h*k21);
k32=(x(i)+0.5*h*k11)*(y(i)+0.5*h*k21)-gama*(z(i)+0.5*h*k31);

k13=sigma*(-(x(i)+0.5*h*k12)+(y(i)+0.5*h*k22));
k23=beta*(x(i)+0.5*h*k12)-(x(i)+0.5*h*k12)*...
(z(i)+0.5*h*k32)-(y(i)+0.5*h*k22);
k33=(x(i)+0.5*h*k12)*(y(i)+0.5*h*k22)-gama*(z(i)+0.5*h*k32);

k14=sigma*(-(x(i)+h*k13)+(y(i)+h*k23));
k24=beta*(x(i)+h*k13)-(x(i)+h*k13)*(z(i)+h*k33)-(y(i)+h*k23);
k34=(x(i)+h*k13)*(y(i)+h*k23)-gama*(z(i)+h*k33);

J1=[-sigma sigma 0;beta-z(i) -1 -x(i);y(i) x(i) -gama];
J2=[-sigma sigma 0;
beta-(z(i)+0.5*h*k31) -1 -(x(i)+0.5*h*k11);
y(i)+0.5*h*k21 x(i)+0.5*h*k11 -gama];
J3=[-sigma sigma 0;
beta-(z(i)+0.5*h*k32) -1 -(x(i)+0.5*h*k12);
y(i)+0.5*h*k22 x(i)+0.5*h*k12 -gama];
J4=[-sigma sigma 0;
beta-(z(i)+h*k33) -1 -(x(i)+h*k13);
y(i)+h*k23 x(i)+h*k13 -gama];

J=I+h*(J1+2*J2*(I+0.5*h*J1)+2*J3*(I+0.5*h*J2*(I+0.5*h*J1))...
+J4*(I+0.5*h*J3*(I+0.5*h*J2*(I+0.5*h*J1))))/6;

B=J*V*S;
[V,S,U]=svd(B);
am=max(diag(S));
S=S/am;
b=b+log(am);

x(i+1)=x(i)+h*(k11+2*k12+2*k13+k14)/6;
y(i+1)=y(i)+h*(k21+2*k22+2*k23+k24)/6;
z(i+1)=z(i)+h*(k31+2*k32+2*k33+k34)/6;

end
for i=1:3
le(i)=(log(diag(S(i,i)))+b)/(k*h);
end
le'



1.3452
0.1502
-22.5502


Saturday, November 25, 2006

切分段线性控制方法2




% 连续时间系统混沌化切控制方法
% 《动力系统的混沌化》陈关荣 汪小帆
%
% Example(函数图象):
% [T,Y]=ode45('chaos2',40,[2;2;1]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
% hold on
% [T,Y]=ode45('chaos2',40,[2;2;-1]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
%
function dx=chaos2(t,x);
a=3;b=20;c=-20;
k=4;d=10;
m=4;e=-10;
A= [ a b 0
-b a 0
0 0 c];
if x(3)+norm(x(1:2))>k & x(3)>0
u=k*[-x(1) -x(2) d]';
elseif x(3)-norm(x(1:2))<-m & x(3)<0
u=m*[-x(1) -x(2) e]';
else
u=0;
end
dx=A*x+u;




连续时间系统混沌化切控制方法41

连续时间系统混沌化切控制方法42


Thursday, November 23, 2006

切分段线性控制方法





% 连续时间系统混沌化切控制方法
% 《动力系统的混沌化》陈关荣 汪小帆
%
% Example(函数图象):
% [T,Y]=ode45('chaos1',40,[2;2;1]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
%
function dx=chaos1(t,x);
a=3;
%global a;
b=20;c=-20;
k=4;d=10;
A= [ a b 0
-b a 0
0 0 c];
if x(3)+norm(x(1:2))>k
u=k*[-x(1) -x(2) d]';
else
u=0;
end
dx=A*x+u;









% 求最大Lyapunov指数谱
% 连续时间系统混沌化切控制方法
% 《动力系统的混沌化》陈关荣 汪小帆
%
Z=[];
global a;
d0=1e-8;
for a=linspace(0.1,6,150);
y0=[2;2;1];
y=[2+d0;2;1];
lsum=0;
p=0;
for k=1:60
[T,Y0] = ode45('chaos1',.4,y0);
[T,Y] = ode45('chaos1',.4,y);
y0 = Y0(size(Y0,1),:);
y = Y(size(Y,1),:);
d1=norm(y-y0);
if d1~=0
y=y0+(d0/d1)*(y-y0);
if k>10
lsum=lsum+log(d1/d0);
end
else
p=p+1;
y(1)=y0(1)+d0;
end
end
le=lsum/(k-10-p);
Z=[Z,a+le*i];
end
plot(Z,'-')
title('Lyapunov exponents')
xlabel('parameter a'),ylabel('Maximum lyapunov exponents')
grid on




切分段线性控制方法
切分段线性控制方法lp

Sunday, October 22, 2006

Cat映射的周期性





function cattest2
% Cat 加密图象
clear all
f=imread('cat.jpg'); % 读曲图象
A=double(f)/255;
[m,n,l]=size(A);

Cat=[321 40;8 1];N=124;

m=m-rem(m,N);
n=n-rem(n,N);
A=A(1:m,1:n,:);
subplot(2,3,1)
imagesc(A,[0 1]);
mk=m/N;nk=n/N;

for k=1:5
% 将图象分割成124*124的数据块
for mki=0:mk-1
for nki=0:nk-1
% 处理N*N的数据块
for i=N*mki:N*mki+123
for j=N*nki:N*nki+123
c=[i j]';
c=Catmap(c,Cat,N);
B(mki*N+c(1)+1,nki*N+c(2)+1,:)=A(i+1,j+1,:);
end
end
end
end
subplot(2,3,k+1)
imagesc(B,[0 1]);
A=B;
end




Cat映射的周期性



相关程序 Cat.m

Cat 映射


Cat.m



function newx=Catmap(x,A,N)
% Cat映射
% x(n+1)=A*x(n) (mod N)
% x是m*1向量,A是m*m矩阵,N是整数
%
% Example:
% global Cat;
% global N;
% A=[1 1;1 2];N=1;
% newx=Catmap([.4 .5]')
%
% Author:LDYU
% Author's email: ustb03-07@yahoo.com.cn
%

mewx=A*x;
newx=mod(mewx,N);



混沌 Cat 映射




function newx=Catmap(x)
% Cat映射
% x(n+1)=x(n)+y(n) (mod 1)
% y(n+1)=x(n)+2*y(n) (mod 1)
%
%
newx(1)=mod(x(1)+x(2),1);
newx(2)=mod(x(1)+2*x(2),1);







function test
clear all
% 读曲图象(路径按自己情况改)
f=imread('cat.jpg');
A=double(f)/255;
[m,n,l]=size(A);
subplot(3,2,1)
imagesc(A,[0 1]);

for k=1:5
S=zeros(m,n);
for i=1:m
for j=1:n
c=[(i-1)/m (j-1)/n];
c=Catmap(c);
ib=fix(m*c(1))+1;
jb=fix(n*c(2))+1;
if S(ib,jb)==0
B(ib,jb,:)=A(i,j,:);
S(ib,jb)=1;
% 如有变换后的点重叠,移到邻近单元,减少信息损失
elseif jb<n & S(ib,jb+1)==0
B(ib,jb+1,:)=A(i,j,:);
S(ib,jb+1)=1;
elseif ib<m & S(ib+1,jb)==0
B(ib+1,jb,:)=A(i,j,:);
S(ib+1,jb)=1;
end
end
end
subplot(3,2,k+1)
imagesc(B,[0 1]);
A=B;
end



变换的原图


变换的原图




变换图


Saturday, October 21, 2006

庞加莱截面和倍周期分岔


以前的关于微分方程组的庞加莱截面和倍周期分岔程序有误,现在更正如下:



Z=[];
[T,Y]=ode45('Lorenz',[0,1000],[1;1;1;16;4;49.52]);
Y(:,3)=Y(:,3)-30;
for k=1:(length(Y)-1)
if Y(k,3)<0
if Y(k+1,3)>0
% 用直线连接位于截面两侧的相邻离散点,近似得到穿过截面的点
x=Y(k,1)-Y(k,3)*(Y(k+1,1)-Y(k,1))/(Y(k+1,3)-Y(k,3));
y=Y(k,2)-Y(k,3)*(Y(k+1,2)-Y(k,2))/(Y(k+1,3)-Y(k,3));
Z=[Z x+i*y];
end

% 如果只画从一侧穿过截面的图象,以下的if-end语句可删去
else
if Y(k+1,3)<0
x=Y(k,1)-Y(k,3)*(Y(k+1,1)-Y(k,1))/(Y(k+1,3)-Y(k,3));
y=Y(k,2)-Y(k,3)*(Y(k+1,2)-Y(k,2))/(Y(k+1,3)-Y(k,3));
Z=[Z x+i*y];
end

end
end
plot(Z,'.','markersize',2)
title('Lorenz 系统的 Poincare 映像 z=30')
xlabel('x'),ylabel('y')










Z=[];
[T,Y]=ode45('Chen',[0,1000],[1;1;1;35;3;28]);
for k=1:(length(Y)-1)
if Y(k,1)<0
if Y(k+1,1)>0
y=Y(k,2)-Y(k,1)*(Y(k+1,2)-Y(k,2))/(Y(k+1,1)-Y(k,1));
z=Y(k,3)-Y(k,1)*(Y(k+1,3)-Y(k,3))/(Y(k+1,1)-Y(k,1));
Z=[Z y+i*z];
end
else
if Y(k+1,1)<0
y=Y(k,2)-Y(k,1)*(Y(k+1,2)-Y(k,2))/(Y(k+1,1)-Y(k,1));
z=Y(k,3)-Y(k,1)*(Y(k+1,3)-Y(k,3))/(Y(k+1,1)-Y(k,1));
Z=[Z y+i*z];
end
end
end
plot(Z,'.','markersize',2)
title('Chen 系统的 Poincare 映像 x=0')
xlabel('y'),ylabel('z')











Z=[];
[T,Y]=ode45('Lorenz3',[0,5000],[0.5;0.5;2;.9;.1;2;0]);
Y(:,2)=Y(:,2)-.2;
for k=1:(length(Y)-1)
if Y(k,2)<0
if Y(k+1,2)>0
x=Y(k,1)-Y(k,2)*(Y(k+1,1)-Y(k,1))/(Y(k+1,2)-Y(k,2));
z=Y(k,3)-Y(k,2)*(Y(k+1,3)-Y(k,3))/(Y(k+1,2)-Y(k,2));
Z=[Z x+i*z];
end
else
if Y(k+1,2)<0
x=Y(k,1)-Y(k,2)*(Y(k+1,1)-Y(k,1))/(Y(k+1,2)-Y(k,2));
z=Y(k,3)-Y(k,2)*(Y(k+1,3)-Y(k,3))/(Y(k+1,2)-Y(k,2));
Z=[Z x+i*z];
end
end
end
plot(Z,'.','markersize',2)
title('模拟 Lorenz 系统的 Poincare 映像 y=0.2')
xlabel('x'),ylabel('z')











Z=[];
for r=linspace(1,500,250);
% 舍弃前面迭带的结果,用后面的结果画图
[T,Y]=ode45('Lorenz',1,[1;1;1;10;8/3;r]);
[T,Y]=ode45('Lorenz',50,Y(length(Y),:));
Y(:,1)=Y(:,2)-Y(:,1);
for k=2:length(Y)
f=k-1;
if Y(k,1)<0
if Y(f,1)>0
y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
Z=[Z r+abs(y)*i];
end
else
if Y(f,1)<0
y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
Z=[Z r+abs(y)*i];
end
end
end
end
plot(Z,'.','markersize',1)
title('Lorenz映射分岔图')
xlabel('r'),ylabel('|y| where x=y')








Saturday, October 14, 2006

Chen 系统的 Poincare 映像



Z=[];
[T,Y]=ode45('Chen',[0,5000],[1;1;1;35;3;28]);
for k=1:length(Y)
if abs(Y(k,1))<1e-2
Z=[Z Y(k,2)+i*Y(k,3)];
end
end
plot(Z,'.','markersize',2)
title('Chen 系统的 Poincare 映像 x=0')
xlabel('y'),ylabel('z')




Chen 系统的 Poincare 映像 x=0





Z=[];
[T,Y]=ode45('Chen',[0,5000],[1;1;1;35;3;28]);
for k=1:length(Y)
if abs(Y(k,3)-28)<1e-2
Z=[Z Y(k,1)+i*Y(k,2)];
end
end
plot(Z,'.','markersize',2)
title('Chen 系统的 Poincare 映像 z=c')
xlabel('x'),ylabel('y')




Chen 系统的 Poincare 映像 z=c


部分放大
Chen 系统的 Poincare 映像 z=c 部分放大


相关程序 Chen.m

模拟 Lorenz 系统的 Poincare 映像



Z=[];
[T,Y]=ode45('Lorenz3',[0,6000],[0.5;0.5;2;.9;.1;2;0]);
for k=1:length(Y)
if (abs(Y(k,2)-0.2)<1e-2)
Z=[Z Y(k,1)+i*Y(k,3)];
end
end
plot(Z,'.','markersize',2)
title('模拟 Lorenz 系统的 Poincare 映像 y=0.2')
xlabel('x'),ylabel('z')




模拟 Lorenz 系统的 Poincare 映像


相关程序 Lorenz3.m

Friday, October 06, 2006

第一阶段个人工作小结


 

第一阶段个人工作小结



第一阶段的任务是查资料及上机实现Lyapunov指数的算法,在这些日子里我查阅
了大量的资料并且上机实现了一些求Lyapunov指数的算法。在查资料的过程中虽然这
方面的东西有很多,但却没有什么现成的程序,网上好不容易找到一些却不是我需要的,
有时找到可以用的程序,运行结果却和理论值有较大的偏差。书上和论文上几乎没给什
么程序,有伪代码的也是十分简略,增加了实现的难度。不过还是发现一些有价值的程
序,这些不全是用Matlab写的,有些还是用Basic写的,这些程序给了我很大的启发。
以前按算法写的程序总得不到正确的结果,现在回头看来很大程度是我理解的意思和书
上要表达的意思不一致,如果书上能提供这些算法的程序,无论是用什么语言实现,也
可以帮助读者理解,而不至于让人走弯路。
虽然遇到了些困难,但实现的程序结果和书上,论文上的结果比较吻合。总之是完
成了第一阶段的任务了,另外今天是中秋加十一放假期间,不熄灯,庆祝一下:)

 

Lyapunov exponents of Chen


Chen2.m




function dX = Chen2(t,X)
%
% Chen吸引子,用来计算Lyapunov指数
% dx/dt=a*(y-x)
% dy/dt=(c-a)*x+c*y-x*z
% dz/dt=x*y-b*z

global a; % 变量不放入参数表中
global b;
global c;

x=X(1); y=X(2); z=X(3);
% Y的三个列向量为相互正交的单位向量
Y = [X(4), X(7), X(10);
X(5), X(8), X(11);
X(6), X(9), X(12)];
% 输出向量的初始化
dX = zeros(12,1);

% Lorenz吸引子
dX(1) = a*(y-x);
dX(2) = (c-a)*x+c*y-x*z;
dX(3) = x*y-b*z;

% Lorenz吸引子的Jacobi矩阵
Jaco = [-a a 0;
c-a-z c -x;
y x -b];

dX(4:12) = Jaco*Y;







Z1=[];
Z2=[];
Z3=[];
global a;
global b;
global c;
b=3;c=28;
for a=linspace(32,40,100);
y=[1;1;1;1;0;0;0;1;0;0;0;1];
lp=0;
for k=1:200
[T,Y] = ode45('Chen2', 1, y);
y = Y(size(Y,1),:);
y0 = [y(4) y(7) y(10);
y(5) y(8) y(11);
y(6) y(9) y(12)];
y0=GS(y0);
mod(1)=norm(y0(:,1));
mod(2)=norm(y0(:,2));
mod(3)=norm(y0(:,3));
lp = lp+log(abs(mod));
y0(:,1)=y0(:,1)/mod(1);
y0(:,2)=y0(:,2)/mod(2);
y0(:,3)=y0(:,3)/mod(3);
y(4:12) = y0';
end
lp=lp/200;
Z1=[Z1 lp(1)];
Z2=[Z2 lp(2)];
Z3=[Z3 lp(3)];
end
a=linspace(32,40,100);
plot(a,Z1,'-',a,Z2,'-',a,Z3,'-');
title('Lyapunov exponents of Chen')
xlabel('b=3,c=28,parameter a'),ylabel('lyapunov exponents')
grid on











Z=[];
d0=1e-8;
for a=linspace(32,40,80)
lsum=0;
x=1;y=1;z=1;
x1=1;y1=1;z1=1+d0;
for i=1:100
[T1,Y1]=ode45('Chen',1,[x;y;z;a;3;28]);
[T2,Y2]=ode45('Chen',1,[x1;y1;z1;a;3;28]);
n1=length(Y1);n2=length(Y2);
x=Y1(n1,1);y=Y1(n1,2);z=Y1(n1,3);
x1=Y2(n2,1);y1=Y2(n2,2);z1=Y2(n2,3);
d1=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);
x1=x+(d0/d1)*(x1-x);
y1=y+(d0/d1)*(y1-y);
z1=z+(d0/d1)*(z1-z);
if i>50
lsum=lsum+log(d1/d0);
end
end
Z=[Z lsum/(i-50)];
end
a=linspace(32,40,80);
plot(a,Z,'-');
title('Chen 系统最大lyapunov指数')
xlabel('parameter a'),ylabel('lyapunov exponents')







相连的程序有Chen.m

Thursday, October 05, 2006

模拟 Lorenz 系统最大lyapunov指数谱


Z=[];
d0=1e-8;
for a=linspace(0.3,2,70)
lsum=0;
x=1;y=1;z=1;
x1=1;y1=1;z1=1+d0;
for i=1:100
[T1,Y1]=ode45('Lorenz3',1,[x;y;z;a;.1;2;0]);
[T2,Y2]=ode45('Lorenz3',1,[x1;y1;z1;a;.1;2;0]);
n1=length(Y1);n2=length(Y2);
x=Y1(n1,1);y=Y1(n1,2);z=Y1(n1,3);
x1=Y2(n2,1);y1=Y2(n2,2);z1=Y2(n2,3);
d1=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);
x1=x+(d0/d1)*(x1-x);
y1=y+(d0/d1)*(y1-y);
z1=z+(d0/d1)*(z1-z);
if i>50
lsum=lsum+log(d1/d0);
end
end
Z=[Z lsum/(i-50)];
end
a=linspace(0.3,2,70);
plot(a,Z,'-');
title('模拟 Lorenz 系统最大lyapunov指数谱')
xlabel('parameter a')
ylabel('lyapunov exponents')




相连的程序有Lorenz3.m

模拟洛仑滋系统

Lorenz3.m




function dx=Lorenz3(t,x);
% 模拟洛仑滋系统[不显含时间t的自治系统]
% dx=Lorenz(t,[x;y;z;a;b;c;m])
% t-可以取任何数字,x,y,z-为自变量,a,b,r-为如下方程所示的参数
% eg: dx=Lorenz(0,[1;1;1;.9;.1;2;0])
%
% 方程如下:
% dx/dt=a*(y-x)
% dy/dt=sign(x)*(c-z)+m
% dz/dt=|x|-b*z
%
%
% Author:yujunjie
% Author's email: ustb03-07@yahoo.com.cn
%

dx(1,1)=x(4)*(x(2)-x(1));
dx(2,1)=sign(x(1))*(x(6)-x(3))+x(7);
dx(3,1)=abs(x(1))-x(5)*x(3);
dx(4,1)=0;
dx(5,1)=0;
dx(6,1)=0;
dx(7,1)=0;



>> [T,Y]=ode45('Lorenz3',[0,10],[1;1;1;.9;.1;2;0]);
x=Y(length(Y),:);
[T,Y]=ode45('Lorenz3',[0,500],x);
plot3(Y(:,1),Y(:,2),Y(:,3));




模拟洛仑滋系统 m=0



>> [T,Y]=ode45('Lorenz3',[0,10],[1;1;1;.9;.1;2;1]);
x=Y(length(Y),:);
[T,Y]=ode45('Lorenz3',[0,500],x);
plot3(Y(:,1),Y(:,2),Y(:,3));




模拟洛仑滋系统 m=1



>> [T,Y]=ode45('Lorenz3',[0,10],[1;1;1;.9;.1;2;-1]);
x=Y(length(Y),:);
[T,Y]=ode45('Lorenz3',[0,500],x);
plot3(Y(:,1),Y(:,2),Y(:,3));




模拟洛仑滋系统 m=-1

系统对初始条件敏感依赖性





[T1,Y1]=ode45('Lorenz',20,[1;1;1;16;4;45.92]);
plot(T1,Y1(:,1),'-');
hold on
[T2,Y2]=ode45('Lorenz',20,[1+.001;1;1;16;4;45.92]);
plot(T2,Y2(:,1),'r--');


相连的程序有Lorenz.m


Lorenz系统对初始条件敏感依赖性

Wednesday, October 04, 2006

Lorenz指数谱


Lorenlp.m




% Lyapunov exponents of Lorenz
%
% Author:yujunjie
% Author's email: ustb03-07@yahoo.com.cn
%
function lorenzlp
Z1=[];
Z2=[];
Z3=[];
global a;
global b;
global r;
a=10;b=8/3;
for r=linspace(0,500,80);
y=[1;1;1;1;0;0;0;1;0;0;0;1];
lp=0;
for k=1:100
[T,Y] = ode45('Lorenz2', 1, y);
y = Y(size(Y,1),:);
y0 = [y(4) y(7) y(10);
y(5) y(8) y(11);
y(6) y(9) y(12)];
y0=GS(y0);
mod(1)=norm(y0(:,1));
mod(2)=norm(y0(:,2));
mod(3)=norm(y0(:,3));
lp = lp+log(abs(mod));
y0(:,1)=y0(:,1)/mod(1);
y0(:,2)=y0(:,2)/mod(2);
y0(:,3)=y0(:,3)/mod(3);
y(4:12) = y0';
end
lp=lp/100;
Z1=[Z1 lp(1)];
Z2=[Z2 lp(2)];
Z3=[Z3 lp(3)];
end
r=linspace(0,500,80);
plot(r,Z1,'-',r,Z2,'-',r,Z3,'-');
title('Lyapunov exponents of Lorenz')
xlabel('parameter r'),ylabel('lyapunov exponents')
grid on




>> lorenzlp


相连的程序有GS.m and Lorenz2.m



Lyapunov exponents of Lorenz

计算Lorenz系统的Lyapunov指数


Lorenz2.m




function dX = Lorenz2(t,X)
%
% Lorenz吸引子,用来计算Lyapunov指数
% dx/dt = a(y-x),
% dy/dt = x(r-z)-y,
% dz/dt = xy-bz,

global a; % 变量不放入参数表中
global b;
global r;

x=X(1); y=X(2); z=X(3);
% Y的三个列向量为相互正交的单位向量
Y = [X(4), X(7), X(10);
X(5), X(8), X(11);
X(6), X(9), X(12)];
% 输出向量的初始化
dX = zeros(12,1);

% Lorenz吸引子
dX(1) = a*(y-x);
dX(2) = x*(r-z)-y;
dX(3) = x*y-b*z;

% Lorenz吸引子的Jacobi矩阵
Jaco = [-a a 0;
r-z -1 -x;
y x -b];

dX(4:12) = Jaco*Y;



llp.m



function llp
global a;
global b;
global r;
a=16;b=4;r=45.92;
y=[1;1;1;1;0;0;0;1;0;0;0;1];
lp=0;
for k=1:300
[T,Y] = ode45('Lorenz2', 1, y);
y = Y(size(Y,1),:);
y0 = [y(4) y(7) y(10);
y(5) y(8) y(11);
y(6) y(9) y(12)];
y0=GS(y0);
mod(1)=norm(y0(:,1));
mod(2)=norm(y0(:,2));
mod(3)=norm(y0(:,3));
lp = lp+log(abs(mod));
y0(:,1)=y0(:,1)/mod(1);
y0(:,2)=y0(:,2)/mod(2);
y0(:,3)=y0(:,3)/mod(3);
y(4:12) = y0';
end
lp=lp/300



>> llp

lp =

1.5096 -0.5995 -21.8105


相连的程序有GS.m

Tuesday, October 03, 2006

Lyapunov exponents of Henon


要用的相关程序 Henon.mGS.m




% Lyapunov exponents of Henon
%
% Author:yujunjie
% Author's email: ustb03-07@yahoo.com.cn
%
function henonlp
Zmax=[];
Zmin=[];
for p=linspace(0,1.4,300);
x=[1;0;p;.3];
e=[1 0
0 1];
lp=0;
for k=1:500
x=Henon(x);
J=[-2*p*x(1) 1
.3 0];
e=GS(J*e);
mod(1)=norm(e(:,1));
mod(2)=norm(e(:,2));
lp = lp+log(abs(mod));
e(:,1)=e(:,1)/mod(1);
e(:,2)=e(:,2)/mod(2);
end
lp=lp/500;
Zmax=[Zmax lp(1)];
Zmin=[Zmin lp(2)];
end
p=linspace(0,1.4,300);
plot(p,Zmax,'-',p,Zmin,'-');
title('Lyapunov exponents of Henon')
xlabel('parameter p'),ylabel('lyapunov exponents')
grid on

>> henonlp




Lyapunov exponents of Henon

GS 正交化算法

GS.m




% G-S正交化
% A = GS(V)
% V 为n*n向量
% A 为n*n正交向量
%
% Example:
% A=GS([1 0 1;2 2 0;3 1 1])
%
% Author:yujunjie
% Author's email: ustb03-07@yahoo.com.cn
%
function A = GS(V)
n=length(V);
for i=1:n
k=0;
for j=1:i-1
k=k+(A(:,j)'*V(:,i)/(A(:,j)'*A(:,j)))*A(:,j);
end
A(:,i)=V(:,i)-k;
end


>> A=GS([1 0 1;2 2 0;3 1 1])

A =

1.0000 -0.5000 0.3810
2.0000 1.0000 0.0952
3.0000 -0.5000 -0.1905


Logistic系统


一个网上的程序




% This file draws the bifurcation diagram and Lyapunov exponent
% for the Logistic Iterated Map x_n+1 = r*x_n.*(1-x_n); 0 < x < 1.
% The parameter r is in interval [0, 4]
% adapted from a program by G Lindblad gli@theophys.kth.se
% http://www.theophys.kth.se/5a1352/mat.html
rect = [200 80 700 650];
set(0, 'defaultfigureposition',rect);

n1=200; %% no of lattice points in coordinate and parameter r
n2=500; %% no of iterations to reach attractor
n3=250; %% no of iterations for bifurcation diagram
n4=4000; %% no of iterations for Lyapunov exponent
k1=[]; kk=[]; q1=[];
home,
if isempty(kk) %> setting default
kk=[3 4];
end

q3=input('> Choose a r-interval [a b] or >>> ');
q=isempty(q3);
if q==0
disp('> The r-interval is set = ');
disp(q3);
kk = q3;
end
disp('> The calculations could take a few seconds .... ');
k0=linspace(kk(1), kk(2) ,n1)';
seed=rand(1);
x1=seed;
ww=[];
for iter1=1:n2 %% this is just to reach the attractor
x1=k0.*x1.*(1-x1);
end
for it3=1:n3 %% points on the attractor, we hope
x1=k0.*x1.*(1-x1);
ww=[ww, x1];
end

% We pick a lattice of n1 values of the parameter k1
k0=linspace(kk(1), kk(2) ,n1)';
% We start the iteration at a random point
seed=rand(1);x1=seed;
% The first part of the iteration is discarded -
% this is just to reach the attractor
for it2=1:n1 %
x1=k0.*x1.*(1-x1);
end
% Now we generate a sequence of approximations to the Lyapunov exponent
% as a function of the vector k0!
aa=log(abs(1-2*x1));
for it3=1:n4 %%
x1=k0.*x1.*(1-x1);
y1=log(abs(1-2*x1));
aa=(y1+it3*aa)/(1+it3);
end
aa = aa + log(k0+eps); % adding a term coming from the k0 factor

subplot(211)
plot(k0,ww, '. k','MarkerSize',4);
axis([kk(1) kk(2) 0 1]);
%axis tight;
title('A bifurcation diagram of the logistic equation');

subplot(212)
plot(k0,[aa, zeros(size(k0))]);
axis([kk(1) kk(2) -1 1]),
title('Lyapunov exponent of the logistic equation');







logistic



Monday, October 02, 2006

微分混沌方程组最大Lyapunov指数算法

Lorenz系统
在参数a=16,b=4,r=46.92下的最大Lyapunov指数




d0=1e-8;
le=0;
lsum=0;
x=1;y=1;z=1;
x1=1;y1=1;z1=1+d0;
for i=1:500
[T1,Y1]=ode45('Lorenz',[0,1],[x;y;z;16;4;46.92]);
[T2,Y2]=ode45('Lorenz',[0,1],[x1;y1;z1;16;4;46.92]);
n1=length(Y1);n2=length(Y2);
x=Y1(n1,1);y=Y1(n1,2);z=Y1(n1,3);
x1=Y2(n2,1);y1=Y2(n2,2);z1=Y2(n2,3);
d1=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);
x1=x+(d0/d1)*(x1-x);
y1=y+(d0/d1)*(y1-y);
z1=z+(d0/d1)*(z1-z);
if i>100
lsum=lsum+log(d1/d0);
end
end
le=lsum/(i-100)




le =

1.4954

Chen系统和Rossler系统也可以用以上的算法只不过将Lorenz换成Chen或Rossler即可。

Lorenz系统

Lorenz.m




function dx=Lorenz(t,x);
% 洛仑滋方程[不显含时间t的自治系统]
% dx=Lorenz(t,[x;y;z;a;b;r])
% t-可以取任何数字,x,y,z-为自变量,a,b,r-为如下方程所示的参数
% eg: dx=Lorenz(0,[1;1;1;16;4;49])
%
% 方程如下:
% dx/dt=a*(y-x)
% dy/dt=x*(r-z)-y
% dz/dt=x*y-b*z
%
% Example(函数图象):
% [T,Y]=ode45('Lorenz',[0,50],[1;1;1;16;4;49]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
%
% Example(分岔图):
%
% Z=[];
% for r=linspace(1,500,350);
% [T,Y]=ode45('Lorenz',[0,50],[1;1;1;10;8/3;r]);
% n=length(Y);
% for k=round(n/2):n
% if abs(Y(k,1))<1
% Z=[Z,r+abs(Y(k,2))*i];
% end
% end
% end
% plot(Z,'.','markersize',1)
% title('Lorenz映射分岔图')
% xlabel('r'),ylabel('|y| where x=0')
%
% Author:yujunjie
% Author's email: ustb03-07@yahoo.com.cn
%

dx(1,1)=x(4)*(x(2)-x(1));
dx(2,1)=x(1)*(x(6)-x(3))-x(2);
dx(3,1)=x(1)*x(2)-x(5)*x(3);
dx(4,1)=0;
dx(5,1)=0;
dx(6,1)=0;






Lorenz函数图象


Lorenz映射分岔图


Chen and Rossler

Chen.m




function dx=Chen(t,x);
% Chen方程[不显含时间t的自治系统]
% dx=Chen(t,[x;y;z;a;b;c])
% t-可以取任何数字,x,y,z-为自变量,a,b,c-为如下方程所示的参数
% eg: dx=Chen(0,[1;1;1;35;3;28])
%
% 方程如下:
% dx/dt=a*(y-x)
% dy/dt=(c-a)*x+c*y-x*z
% dz/dt=x*y-b*z
%
% Example(函数图象):
% [T,Y]=ode45('Chen',[0,50],[1;1;1;35;3;28]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
%
% Author's email: ustb03-07@yahoo.com.cn
%
dx(1,1)=x(4)*(x(2)-x(1));
dx(2,1)=(x(6)-x(4))*x(1)+x(6)*x(2)-x(1)*x(3);
dx(3,1)=x(1)*x(2)-x(5)*x(3);
dx(4,1)=0;
dx(5,1)=0;
dx(6,1)=0;






Chen




Rossler.m




function dx=Rossler(t,x);
% Rossler方程[不显含时间t的自治系统]
% dx=Rossler(t,[x;y;z;a;b;c])
% t-可以取任何数字,x,y,z-为自变量,a,b,c-为如下方程所示的参数
% eg: dx=Rossler(0,[1;1;1;.15;.2;10])
%
% 方程如下:
% dx/dt=-y-z
% dy/dt=x+a*y
% dz/dt=b+z*(x-c)
%
% Example(函数图象):
% [T,Y]=ode45('Rossler',[0,400],[1;1;1;.15;.2;10]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
%
% Author's email: ustb03-07@yahoo.com.cn
%
dx(1,1)=-x(2)-x(3);
dx(2,1)=x(1)+x(4)*x(2);
dx(3,1)=x(5)+x(3)*(x(1)-x(6));
dx(4,1)=0;
dx(5,1)=0;
dx(6,1)=0;







Rossler

Henon系统 (2)

Henon.m






function newx=Henon(x);
% Henon方程[差分方程]
% newx=Henon([x;y;p;q])
%
% 方程如下:
% x(k+1)=-p*x(k)*x(k)+y(k)+1
% y(k+1)=q*x(k)
%
% Example(函数图象):
%
% X=[1;0;1.4;.3];Y=[];
% for i=1:1000
% X=feval(@Henon,X);
% Y(i,:)=X(1:2,1);
% end
% plot(Y(:,1),Y(:,2),'.');
% title('Henon映射图'),xlabel('x'),ylabel('y')
%
% Example(分岔图):
%
% Z=[];
% for p=linspace(0,1.4,300);
% x=[1;0;p;.3];
% for k=1:300;
% x=Henon(x);
% if k>60
% Z=[Z,p+x(1)*i];
% end
% end
% end
% plot(Z,'.','markersize',1)
% title('Henon映射分岔图'),xlabel('p'),ylabel('x')
%
% Example(最大Lyapunov指数谱图):
%
% d0=1e-8;
% Z=[];
% for p=linspace(0,1.4,300)
% le=0;
% lsum=0;
% x=[1;0;p;.3];
% x1=[1;d0;p;.3];
% for k=1:800
% x=Henon(x);
% x1=Henon(x1);
% d1=sqrt((x(1)-x1(1))^2+(x(2)-x1(2))^2);
% x1=x+(d0/d1)*(x1-x);
% if k>100
% lsum=lsum+log(d1/d0);
% end
% end
% le=lsum/(k-100);
% Z=[Z,p+le*i];
% end
% plot(Z,'-')
% title('Henon最大Lyapunov指数谱图'),xlabel('p'),ylabel('lyapunov')
% grid on
%
% Author:yujunjie
% Author's email: ustb03-07@yahoo.com.cn


newx(1,1)=-x(3)*x(1)*x(1)+x(2)+1;
newx(2,1)=x(4)*x(1);
newx(3,1)=x(3);
newx(4,1)=x(4);







Henon映射图


Henon映射分岔图


henon最大Lyapunov指数谱图

Lozi系统

Lozi.m




function newx=Lozi(x);
% Lozi方程[差分方程]
% newx=Lozi([x;y;p;q]);
% 方程如下:
% x(k+1)=-p*x(k)+y(k)+1
% y(k+1)=q*x(k)
%
% Example(函数图象):
%
% X=[1;0;1.7;.5];Y=[];
% for i=1:1000
% X=feval(@Lozi,X);
% Y(i,:)=X(1:2,1);
% end
% plot(Y(:,1),Y(:,2),'.');
% title('Lozi映射图'),xlabel('x'),ylabel('y')
%
% Example(分岔图):
%
% Z=[];
% for p=linspace(0,1.7,300);
% x=[1;0;p;.5];
% for k=1:300;
% x=Lozi(x);
% if k>60
% Z=[Z,p+x(1)*i];
% end
% end
% end
% plot(Z,'.','markersize',1)
% title('Lozi映射分岔图'),xlabel('p'),ylabel('x')
%
% Example(最大Lyapunov指数谱图):
%
% d0=1e-8;
% Z=[];
% for p=linspace(0,1.7,300)
% le=0;
% lsum=0;
% x=[1;0;p;.5];
% x1=[1;d0;p;.5];
% for k=1:800
% x=Lozi(x);
% x1=Lozi(x1);
% d1=sqrt((x(1)-x1(1))^2+(x(2)-x1(2))^2);
% x1=x+(d0/d1)*(x1-x);
% if k>100
% lsum=lsum+log(d1/d0);
% end
% end
% le=lsum/(k-100);
% Z=[Z,p+le*i];
% end
% plot(Z,'-')
% title('Lozi最大Lyapunov指数图'),xlabel('p'),ylabel('lyapunov')
% grid on
%
% Author:yujunjie
% Author's email: ustb03-07@yahoo.com.cn

newx(1,1)=-x(3)*abs(x(1))+x(2)+1;
newx(2,1)=x(4)*x(1);
newx(3,1)=x(3);
newx(4,1)=x(4);



Lozi映射图


Lozi映射分岔图


Lozi最大Lyapunov指数图

Friday, September 22, 2006

推荐一个学习混沌的网站

这上面有大量的关于混沌和分形matlab程序,大家可以上去看看。

Thursday, September 21, 2006

Henon系统

function newx=Henon(x);
% Henon方程[差分方程]
% newx=Henon(x)
%
% 方程如下:
%   x(k+1)=-p*x(k)*x(k)+y(k)+1
%   y(k+1)=q*x(k)
%
% Example:
%   Y=Iterative(1000,[1;0]);
%   plot(Y(:,1),Y(:,2),'.');
%
% Author's email: ustb03-07@yahoo.com.cn
%
p=1.4;
q=.3;
newx=zeros(2,1);
 
newx(1)=-p*x(1)*x(1)+x(2)+1;
newx(2)=q*x(1);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
function Y=Iterative(n,X);
% 方程迭代
% 参数如下:
%   Y=Iterative(n,X);
%   n-迭代次数
%   X-函数的初值
%
% Example:
%   Y=Iterative(1000,[1;0]);
%
% Author's email: ustb03-07@yahoo.com.cn
%
for i=1:n
   
% 对Henon方程迭代
  X=feval(@Henon,X);
  
  Y(i,:)=X;
end

Rossler系统

function dx=Rossler(t,x);
% Rossler方程[不显含时间t的自治系统]
% 方程如下:
% dx/dt=-y-z
% dy/dt=x+a*y
% dz/dt=b+z*(x-c)
%
% Example:
% [T,Y]=ode45('Rossler',[0,400],[1;1;1]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
%
% Author's email: ustb03-07@yahoo.com.cn
%
a=.15;
b=.2;
c=10;
dx=zeros(3,1);

dx(1)=-x(2)-x(3);
dx(2)=x(1)+a*x(2);
dx(3)=b+x(3)*(x(1)-c);

Chen系统

function dx=Chen(t,x);
% Chen方程[不显含时间t的自治系统]
% 方程如下:
% dx/dt=a*(y-x)
% dy/dt=(c-a)*x+c*y-x*z
% dz/dt=x*y-b*z
%
% Example:
% [T,Y]=ode45('Chen',[0,50],[1;1;1]);
% plot3(Y(:,1),Y(:,2),Y(:,3));
%
% Author's email: ustb03-07@yahoo.com.cn
%
a=35;
b=3;
c=28;
dx=zeros(3,1);

dx(1)=a*(x(2)-x(1));
dx(2)=(c-a)*x(1)+c*x(2)-x(1)*x(3);
dx(3)=x(1)*x(2)-b*x(3);

Copyright © 2006 LDYU (USTB OF CHINA)