Rossler 混沌系统Lyapunov指数Matlab实现

% Author: Thomas Lee

H4j\/O2Ci0% E-mail: lixf1979@126.com % Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China

 振动资讯tjXAf&t,K

6B] vn I0function dX = Rossler_ly(t,X)振动资讯'njJ![iC
%  Rossler吸引子，用来计算Lyapunov指数
'|1`ue\ei1Z(H0%        a=0.15,b=0.20,c=10.0
1V7\lx:rZn5L'he,h0%        dx/dt = -y-z,
!Yn Y
{(f U(E$b(a0%        dy/dt = x+ay,振动资讯`%Le @'G7Fh
%        dz/dt = b+z(x-c),
A-m&\`F!j;@+S5v5k0a = 0.15;振动资讯,~Y!c0{-L.T|6L
b = 0.20;
1Sbr4G;e"h7D%b?0c = 10.0;振动资讯u T[.]4v6}
x=X(1); y=X(2); z=X(3);
H k0QJH
t0% Y的三个列向量为相互正交的单位向量振动资讯*C
~&W'z
OI-I
Y = [X(4), X(7), X(10);
,TS$rpq)yy;_}i0    X(5), X(8), X(11);振动资讯Prw9k9L;_;c
    X(6), X(9), X(12)];
Hiw U ~~?-Q+s0% 输出向量的初始化，必不可少
/G/vy `3l6a:L T0dX = zeros(12,1);
~ z-B H B0% Rossler吸引子
7\a!PNW0dX(1) = -y-z;振动资讯)h2x$ET'{s/Jg
dX(2) = x+a*y;振动资讯/U;O kG4Eb q/X2m%f
dX(3) = b+z*(x-c);
.G5RF X7Y/N0% Rossler吸引子的Jacobi矩阵
9iB;oo7J@%mh;y0Jaco = [0 -1 -1;振动资讯 L3_&~ T`'s0KxM
        1 a  0;振动资讯$X.N-z,U/KG
        z 0  x-c];
5S6O4Q;A7K7fDb6T0dX(4:12) = Jaco*Y;振动资讯+aV-z%mPT

求解LE代码：振动资讯da&U7XV;@
% 计算Rossler吸引子的Lyapunov指数振动资讯2COZj.|/z,`e$f7Gt
clear;
)@.fP p"m+u^0yinit = [1,1,1];
1xr0F:z D0orthmatrix = [1 0 0;振动资讯?7M
q3h S6z Su
              0 1 0;振动资讯I8z oF!d0P8Z0G iz
              0 0 1];
&ft}7@m5~"`3t0pi0a = 0.15;
{J+?Q^-Lj7M.G0b = 0.20;
I d/UJ T!a0c = 10.0;
j'jCl&d9z%E| }0y = zeros(12,1);
l FpY%R5Hpr9A0% 初始化输入
OAG,zy'iE{0y(1:3) = yinit;振动资讯M!a-oi&^G+b
y(4:12) = orthmatrix;振动资讯 uEf!~Z
tstart = 0; % 时间初始值
4u*q
Du;j0uaB'}0tstep = 1e-3; % 时间步长
aN_TxkUI0wholetimes = 1e5; % 总的循环次数振动资讯 FA^`$@H+C5z4Q
steps = 10; % 每次演化的步数 振动资讯iw&nE'R|_
iteratetimes = wholetimes/steps; % 演化的次数振动资讯0Pv|.bYQ
mod = zeros(3,1);振动资讯*a)h{p7s3K S4M\I,{%c
lp = zeros(3,1);
[Y1X#BO@'j(K@0% 初始化三个Lyapunov指数
{G sb7ji!n0Lyapunov1 = zeros(iteratetimes,1);
1U;GvD.~:y8cmD
eO0Lyapunov2 = zeros(iteratetimes,1);振动资讯YmL tb
Lyapunov3 = zeros(iteratetimes,1);振动资讯c@k/b*pF;gU1j
for i=1:iteratetimes
$wtB%J"g0    tspan = tstart:tstep:(tstart + tstep*steps);   振动资讯3z[4f8^ ug
    [T,Y] = ode45('Rossler_ly', tspan, y);
!a RlH3NUk(m"s0    % 取积分得到的最后一个时刻的值振动资讯/\,|.Sg'A0~3PeH
    y = Y(size(Y,1),:);振动资讯?Dok${2^D
    % 重新定义起始时刻
t[&Br4m\,P J,~ v0    tstart = tstart + tstep*steps;
7i.a-xzk&a8E0    y0 = [y(4) y(7) y(10);
%TF3CB"V0          y(5) y(8) y(11);振动资讯n;Zeg!p#VN
          y(6) y(9) y(12)];振动资讯4|;i&O@vFlFD1m8?Z
    %正交化
}Bx9`Qz:e[ v0    y0 = ThreeGS(y0);
'e#@\;W QNp$z0    % 取三个向量的模
lU'vl\3yEJ0    mod(1) = sqrt(y0(:,1)'*y0(:,1));振动资讯#@L/z AA2|
    mod(2) = sqrt(y0(:,2)'*y0(:,2));振动资讯l5^7DJ t(L0C|Q$Z
    mod(3) = sqrt(y0(:,3)'*y0(:,3));振动资讯"hwE1isAUux0?
^
    y0(:,1) = y0(:,1)/mod(1);振动资讯.Os'a"O*BCk0Z
    y0(:,2) = y0(:,2)/mod(2);振动资讯t n;y4W Q
E*vq0h8H
    y0(:,3) = y0(:,3)/mod(3);振动资讯}1Z;F9rgc#C%bN
    lp = lp+log(abs(mod));振动资讯`"B$X(zQ.o_2r T u
    %三个Lyapunov指数
jdP'D/r3n!~}U0    Lyapunov1(i) = lp(1)/(tstart);振动资讯gmxt T2d
    Lyapunov2(i) = lp(2)/(tstart);
m[M2_:G0    Lyapunov3(i) = lp(3)/(tstart);振动资讯7O`~&hPS"Q
        y(4:12) = y0';
WSS,Y0t0`IaS0end振动资讯O W4N6z V;Cp
% 作Lyapunov指数谱图
B#w?$q-i|hpcx#H"[0i = 1:iteratetimes;
e Ma:KaLqk0plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)

程序中用到的ThreeGS程序如下：振动资讯$v%v1in5~WC
%G-S正交化振动资讯!M0U(Ym/RZ9g3jI"b_}
function A = ThreeGS(V)  % V 为3*3向量振动资讯1^ m L`@ vR
v1 = V(:,1);振动资讯SE#C1R0n*l2x
v2 = V(:,2);
dXti"H7_0v3 = V(:,3);振动资讯"S0?E~w,MC
a1 = zeros(3,1);振动资讯$w#F0[ W?ow
a2 = zeros(3,1);振动资讯1m8d[0R5N)_
a3 = zeros(3,1);振动资讯!g S5xw/j
a1 = v1;
$uy*e5`leL eS0a2 = v2-((a1'*v2)/(a1'*a1))*a1;
bp,}e2X ^zxs0a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;振动资讯L0n.IYtH-f
A = [a1,a2,a3];

9JB-JCU;]8@'v0计算得到的Rossler系统的LE为————  0.063231  0.092635  -9.8924