DOC《1. 关于连续系统Lyapunov指数的计算方法》

1.

关于连续系统Lyapunov指数的计算方法 连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法,或者通过求解系统的微分方程,得到微分方程解的时间序列,然后利用时间序列(即离散系统)的LE求解方法来计算得到. 关于连续系统LE的计算,主要以定义方法、Jacobian方法做主要介绍内容. (1)定义法 关于定义法求解的程序,和matlab板块的"连续系统LE求解程序"差不多.以Rossler系统为例 Rossler系统微分方程定义程序 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:tsteptstart + 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为DDDD 0.063231 0.092635 -9.8924 Wolf文章中计算得到的Rossler系统的LE为DDDD0.09

0 -9.77 需要注意的是DD定义法求解的精度有限,对有些系统的计算往往出现计果和理论值有偏差的现象. 正交化程序可以根据上面的扩展到N*N向量,这里就不加以说明了,对matlab用户来说应该还是比较简单的! (2)Jacobian方法 通过资料检索,发现论坛中用的较多的LET工具箱的算法原理就是Jacobian方法. 基本原理就是首先求解出连续系统微分方程的近似解,然后对系统的Jacobian矩阵进行QR分解,计算Jacobian矩阵特征值的乘积,最后计算出LE和分数维. 经过计算也证明了这种方法精度较高,对目前常见的混沌系统,如Lorenz、Henon、Duffing等的Lyapunov指数的计算精度都很高,而且程序编写有一定的规范,个人很推荐使用.(虽然我自己要做的系统并不适用) LET工具箱可以在网络上找到,这里就不列出了!关于LET工具箱如果有问题,欢迎加入本帖讨论! 对离散动力系统,或者说是非线性时间序列,往往不需要计算出所有的Lyapunov指数,通常只需计算出其最大的Lyapunov指数即可."1983年,格里波基证明了只要最大Lyapunov指数大于零,就可以肯定混沌的存在". 目前常用的计算混沌序列最大Lyapunov指数的方法主要有一下几种: (1)由定义法延伸的Nicolis方法 (2)Jacobian方法 (3)Wolf方法 (4)P-范数方法 (5)小数据量方法 其中以Wolf方法和小数据量方法应用最为广泛,也最为普遍. 下面对Nicolis方法、Wolf方法以及小数据量方法作一一介绍. (1)Nicolis方法 这种方法和连续系统的定义方法类似,而且目前应用很有限制,因此只对其理论进行介绍,编程应用方面就省略了 (2)Wolf方法 Wolf方法的Matlab程序如下: function lambda_1=lyapunov_wolf(data,N,m,tau,P) % 该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法 % m: 嵌入维数 % tau:时间延迟 % data:时间序列 % N:时间序列长度 %P:时间序列的平均周期,选择演化相点距当前点的位置差,即若当前相点为I,则演化相点只能在|I-J|>P的相点中搜寻 % lambda_1:返回最大lyapunov指数值 min_point=1 ;