Rossler 混沌系统Lyapunov指数Matlab实现

% Author: Thomas Lee振动资讯R!CZ{Q6rq

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

function dX = Rossler_ly(t,X)振动资讯5aE Vl~
%  Rossler吸引子，用来计算Lyapunov指数
7? |o H#`0M0%        a=0.15,b=0.20,c=10.0振动资讯C6}yS!la|(h2H!d+c"p
%        dx/dt = -y-z,振动资讯;\q~eB_e)W
%        dy/dt = x+ay,
Se#f'q-a,~b0%        dz/dt = b+z(x-c),振动资讯"W!qg| N4jV;u)@|
a = 0.15;振动资讯6KRC.e {)c`+g
b = 0.20;
Y2M,t{,f1G0c = 10.0;
Q|!I.tP"{0x=X(1); y=X(2); z=X(3);振动资讯0AcI!}3Ob
% Y的三个列向量为相互正交的单位向量
'bC4B7`b0Y = [X(4), X(7), X(10);
F0W/An/E/n;w!\0    X(5), X(8), X(11);
'?T {D(D*]MC0    X(6), X(9), X(12)];
h5J B:} {r0% 输出向量的初始化，必不可少
5L/PNtUA"[7e0dX = zeros(12,1);振动资讯&w2QLTA)s _ \$I
% Rossler吸引子
c6b#t%W5uE:Ih5Y0dX(1) = -y-z;
:]Qz!R,lg4ssE%y0dX(2) = x+a*y;振动资讯R-z*z0B_ Q,w
dX(3) = b+z*(x-c);振动资讯?0xRaw9\U6{&Q
% Rossler吸引子的Jacobi矩阵
ZX!XM8N5n K_9U0Jaco = [0 -1 -1;
4l*Ep$T6b0G5o0        1 a  0;
}6UJ _ j)u&K0        z 0  x-c];
h6J8Q4Y1PWi~0dX(4:12) = Jaco*Y;振动资讯&ju,]#R[hw

,y G/J;B [Q9p2i0求解LE代码：振动资讯b%l3[%R"f&T
% 计算Rossler吸引子的Lyapunov指数
:A/V5y oM~y0clear;振动资讯rVD*pn+t^,j/t-P n
yinit = [1,1,1];振动资讯 B`;[$Mv5`X
orthmatrix = [1 0 0;振动资讯[)At-bO}t:ey4E
              0 1 0;振动资讯#i-?X3p Ka;H
              0 0 1];振动资讯-d8cINo%Z/y3G
a = 0.15;
!d e3Ka6xkR*I!V`0b = 0.20;振动资讯opl.}8a} _
c = 10.0;
Z4z:Kuo)uF[FO'q0y = zeros(12,1);
Xu3M#H[#Hx:W"I0% 初始化输入振动资讯FVaO&}
y(1:3) = yinit;
u;c*{ q1c"u GnT u0y(4:12) = orthmatrix;
Pe#[4f;_k L L
FD0tstart = 0; % 时间初始值
he6d;tZT/@0tstep = 1e-3; % 时间步长
|(Q:e!Ol0wholetimes = 1e5; % 总的循环次数振动资讯
}$kc4w)qL6j
steps = 10; % 每次演化的步数 振动资讯%_4ms6UB0gu8{ T
iteratetimes = wholetimes/steps; % 演化的次数振动资讯]9o4g a U"i
mod = zeros(3,1);振动资讯8jhL @(kUwq$v
lp = zeros(3,1);
~0D8Ow5^6SmbDi7r0% 初始化三个Lyapunov指数振动资讯0`-Av:e#J/?KyR9L2kU#I
Lyapunov1 = zeros(iteratetimes,1);振动资讯nJ4`Tgx6kB-A
Lyapunov2 = zeros(iteratetimes,1);
D:CiKH3Pm0Lyapunov3 = zeros(iteratetimes,1);振动资讯W7~+~7U"w:i5r
for i=1:iteratetimes
SG];W;YX@ne0    tspan = tstart:tstep:(tstart + tstep*steps);   振动资讯-l1Br2E ew$N
    [T,Y] = ode45('Rossler_ly', tspan, y);
EaNL}qNj:|3s)Lw0    % 取积分得到的最后一个时刻的值
/Y4pA inz-jlIuT0    y = Y(size(Y,1),:);
z
y-y/|W#n3n
X{Lkv0    % 重新定义起始时刻振动资讯,O
h:R+l*olD5lM _5K
    tstart = tstart + tstep*steps;
$ok}/M|PWEW;B0    y0 = [y(4) y(7) y(10);
+aIb*nG'R%x:[m0          y(5) y(8) y(11);振动资讯1G;V.|MFT
          y(6) y(9) y(12)];
:^:I3{6d8`,onU0    %正交化
^4R `|:n$e#Y0    y0 = ThreeGS(y0);
"UM[.R$t+O^.V U0    % 取三个向量的模振动资讯d:@~5z'\d
    mod(1) = sqrt(y0(:,1)'*y0(:,1));
/}Tv!\K*B0    mod(2) = sqrt(y0(:,2)'*y0(:,2));
e4p!yGZct"|&\0    mod(3) = sqrt(y0(:,3)'*y0(:,3));
'N-R"N8n'fN:fe0    y0(:,1) = y0(:,1)/mod(1);振动资讯"vo2YIB
    y0(:,2) = y0(:,2)/mod(2);
[v f!\7xg1x#@j0    y0(:,3) = y0(:,3)/mod(3);振动资讯*u&n c8[$j-jZq@-Y
    lp = lp+log(abs(mod));
9r$J`
WX0    %三个Lyapunov指数
*r D$_E6_jj0    Lyapunov1(i) = lp(1)/(tstart);振动资讯Ng.~R5uWJ }
    Lyapunov2(i) = lp(2)/(tstart);振动资讯6v)qCLZ/j^+k
    Lyapunov3(i) = lp(3)/(tstart);
3Et!Kf+C#Ou?0o0        y(4:12) = y0';振动资讯6K'_2KK&Xl
end振动资讯6p+u;KDv)Y9x
% 作Lyapunov指数谱图振动资讯Lu g]4[h
i = 1:iteratetimes;振动资讯2X Ha!n&J:S;n)MO5}
plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)

R
}*{Zig({F0程序中用到的ThreeGS程序如下：振动资讯$m^0Z3n s3}bZA
%G-S正交化振动资讯C"A!i!Q!}
function A = ThreeGS(V)  % V 为3*3向量
|e9rtq9Y0v1 = V(:,1);振动资讯5P2?!x2s/X ` J?
v2 = V(:,2);
p8QI(Hs^
A0Vo0v3 = V(:,3);振动资讯q/u \Jw:e
a1 = zeros(3,1);
@yV B6@_6~0a2 = zeros(3,1);
3TG'zi6Y1Dk0a3 = zeros(3,1);
b:h-o o;CS0a1 = v1;振动资讯]/zv$w L
a2 = v2-((a1'*v2)/(a1'*a1))*a1;
(Gd/K*]OL(P0a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;振动资讯9qim {~c4m6j~
A = [a1,a2,a3];

计算得到的Rossler系统的LE为————  0.063231  0.092635  -9.8924振动资讯 bP rI:V