通过雅可比矩阵求解Rossler系统的李雅普诺夫指数指数

function main %求解LE代码： % 计算Rossler吸引子的Lyapunov指数 % Rossler吸引子，用来计算Lyapunov指数 % a=0.25,b=1,c=5.5 % dx/dt = -y-z, % dy/dt = x+ay, % dz/dt = b+z(x-c), clear; yinit = [1,1,1]; orthmatrix = [1 0 0; 0 1 0; 0 0 1]; a = 0.25; b = 1; c = 5.5; 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); % 取积分得到的最后一个时刻的值 function main %求解LE代码： % 计算Rossler吸引子的Lyapunov指数 % Rossler吸引子，用来计算Lyapunov指数 % a=0.25,b=1,c=5.5 % dx/dt = -y-z, % dy/dt = x+ay, % dz/dt = b+z(x-c), clear; yinit = [1,1,1]; orthmatrix = [1 0 0; 0 1 0; 0 0 1]; a = 0.25; b = 1; c = 5.5; 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) Le1=Lyapunov1(end),Le2=Lyapunov2(end),Le3=Lyapunov3(end) title('Lyapunov指数谱图') end %计算得到的Rossler系统的LE为———— 0.0992 0.0755 -5.2644 %计算得到的Rossler系统的LE为———— 0.0992 0.0755 -5.2644 function dX = Rossler_ly(t,X) a = 0.25; b = 1; c = 5.5; 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; %程序中用到的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];