Rossler_混沌系统Lyapunov指数Matlab实现代码

水面清圆

Rossler 混沌系统Lyapunov指数Matlab实现代码  分类： 转载 2010-04-11 10:49    混沌系统的基本特点就是系统对初始值的极端敏感性，两个相差无几的初值所产生的轨迹，随着时间的推移按指数方式分离，lyapunov指数就是定量的描述这一现象的量。

　　Lyapunov指数是衡量系统动力学特性的一个重要定量指标,它表征了系统在相空间中相邻轨道间收敛或发散的平均指数率。对于系统是否存在动力学混沌, 可以从最大Lyapunov指数是否大于零非常直观的判断出来: 一个正的Lyapunov指数,意味着在系统相空间中,无论初始两条轨线的间距多么小,其差别都会随着时间的演化而成指数率的增加以致达到无法预测,这就是混沌现象。

　　Lyapunov指数的和表征了椭球体积的增长率或减小率,对Hamilton 系统,Lyapunov指数的和为零; 对耗散系统,Lyapunov指数的和为负。如果耗散系统的吸引子是一个不动点,那么所有的Lyapunov指数通常是负的。如果是一个简单的m维流形(m = 1或m = 2分别为一个曲线或一个面) ,那么,前m 个Lyapunov指数是零,其余的Lyapunov指数为负。不管系统是不是耗散的,只要λ1 > 0就会出现混沌。

　　微分动力系统L yapunov指数的性质

　　对于一维(单变量) 情形,吸引子只可能是不动点(稳定定态) 。此时λ是负的。对于二维情形, 吸引子或者是不动点或者是极限环。对于不动点,任意方向的δxi , 都要收缩, 故这时两个Lyapunov指数都应该是负的, 即对于不动点, (λ1 ,λ2 ) = ( - , - ) 。至于极限环,如果取δxi 始终是垂直于环线的方向,它一定要收缩,此时λ < 0;当取δxi沿轨道切线方向,它既不增大也不缩小,可以想像,这时λ = 0。事实上,所有不终止于定点而又有界的轨道(或吸引子) 都至少有一个Lyapunov指数等于零,它表示沿轨线的切线方向既无扩展又无收缩的趋势。所以极限环的Lyapunov指数是(λ1 ,λ2 ) = (0, - ) 。

　　在三维情形下有

　　(λ1 ,λ2 ,λ3 ) = ( - , - , - ) :稳定不动点;

　　(λ1 ,λ2 ,λ3 ) = (0, - , - ) :极限环;

　　(λ1 ,λ2 ,λ3 ) = (0, 0, - ) :二维环面;

　　(λ1 ,λ2 ,λ3 ) = ( +, +, 0) :不稳极限环;

　　(λ1 ,λ2 ,λ3 ) = ( +, 0, 0) :不稳二维环面;

　　(λ1 ,λ2 ,λ3 ) = ( +, 0, - ) :奇怪吸引子。

　　李雅谱诺夫指数小于零，则意味着相邻点最终要靠拢合并成一点，这对应于稳定的不动点和周期运动;若指数大于零，则意味着相邻点最终要分离，这对应于轨道的局部不稳定，如果轨道还有整体的稳定因素（如整体有界、耗散、存在捕捉区域等），则在此作用下反复折叠并形成混沌吸引子。指数越大，说明混沌特性越明显，混沌程度越高.   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Thomas Lee
function dX = Rossler_ly(t,X)
% Rossler吸引子，用来计算Lyapunov指数
%        a=0.15,b=0.20,c=10.0
%        dx/dt = -y-z,
%        dy/dt = x+ay,
%        dz/dt = b+z(x-c),
a = 0.15;
b = 0.20;
c = 10.0;
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);
% Rossler吸引子
dX(1) = -y-z;
dX(2) = x+a*y;
dX(3) = b+z*(x-c);
% Rossler吸引子的Jacobi矩阵
Jaco = [0 -1 -1;
        1 a 0;
        z 0 x-c];
dX(4:12) = Jaco*Y;

求解LE代码：
% 计算Rossler吸引子的Lyapunov指数
clear;
yinit = [1,1,1];
orthmatrix = [1 0 0;
              0 1 0;
              0 0 1];
a = 0.15;
b = 0.20;
c = 10.0;
y = zeros(12,1);
% 初始化输入
y(1:3) = yinit;
y(4:12) = orthmatrix;
tstart = 0; % 时间初始值
tstep = 1e-3; % 时间步长
wholetimes = 1e5; % 总的循环次数
steps = 10; % 每次演化的步数
iteratetimes = wholetimes/steps; % 演化的次数
mod = zeros(3,1);
lp = zeros(3,1);
% 初始化三个Lyapunov指数
Lyapunov1 = zeros(iteratetimes,1);
Lyapunov2 = zeros(iteratetimes,1);
Lyapunov3 = zeros(iteratetimes,1);
for i=1:iteratetimes
    tspan = tstart:tstep:(tstart + tstep*steps);   
    [T,Y] = ode45('Rossler_ly', tspan, y);
    % 取积分得到的最后一个时刻的值
    y = Y(size(Y,1),:);
    % 重新定义起始时刻
    tstart = tstart + tstep*steps;
    y0 = [y(4) y(7) y(10);
          y(5) y(8) y(11);
          y(6) y(9) y(12)];
    %正交化
    y0 = ThreeGS(y0);
    % 取三个向量的模
    mod(1) = sqrt(y0(:,1)'*y0(:,1));
    mod(2) = sqrt(y0(:,2)'*y0(:,2));
    mod(3) = sqrt(y0(:,3)'*y0(:,3));
    y0(:,1) = y0(:,1)/mod(1);
    y0(:,2) = y0(:,2)/mod(2);
    y0(:,3) = y0(:,3)/mod(3);
    lp = lp+log(abs(mod));
    %三个Lyapunov指数
    Lyapunov1(i) = lp(1)/(tstart);
    Lyapunov2(i) = lp(2)/(tstart);
    Lyapunov3(i) = lp(3)/(tstart);
        y(4:12) = y0';
end
% 作Lyapunov指数谱图
i = 1:iteratetimes;
plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)
程序中用到的ThreeGS程序如下：
%G-S正交化
function A = ThreeGS(V) % V 为3*3向量
v1 = V(:,1);
v2 = V(:,2);
v3 = V(:,3);
a1 = zeros(3,1);
a2 = zeros(3,1);
a3 = zeros(3,1);
a1 = v1;
a2 = v2-((a1'*v2)/(a1'*a1))*a1;
a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;
A = [a1,a2,a3];
计算得到的Rossler系统的LE为———— 0.063231 0.092635 -9.8924

水面清圆

牛贴 不错。

wxfx1988

好贴，先收藏了

猫头鹰先生

谢谢分享啊