基于奇异值分解的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

切分段线性控制方法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;

切分段线性控制方法

% 连续时间系统混沌化切控制方法 % 《动力系统的混沌化》陈关荣 汪小帆 % % 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

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.m