lyapunov指数总结


http://blog.sina.com.cn/s/blog_61c41eb80100gfkh.html 【总结】Lyapunov 指数(LE)的计算方法 转自“百思论坛”帖子: http://www.baisi.net/viewthread.php?tid=797646 近期为了把计算 LE 的一些问题弄清楚,看了有 7~9 本书!下面以吕金虎《混沌时间序 列分析及其应用》 、马军海《复杂非线性系统的重构技术》为主线,把目前已有的 LE 计算 方法做一个汇总!

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: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 Wolf 文章中计算得到的 Rossler 系统的 LE 为――――0.09 0 -9.77

需要注意的是――定义法求解的精度有限, 对有些系统的计算往往出现计果和理论值有偏差 的现象。 正交化程序可以根据上面的扩展到 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:时间序列长度 % :时间序列的平均周期,选择演化相点距当前点的位置差,即若当前相点为 I,则演化相 的相点中搜寻

点只能在|I-J|>

% lambda_1:返回最大 lyapunov 指数值 min_point=1 ; %&&要求最少搜索到的点数 MAX_CISHU=5 ; %&&最大增加搜索范围次数 %FLYINGHAWK

% 求最大、最小和平均相点距离 max_d = 0; %最大相点距离 min_d = 1.0e+100; %最小相点距离 avg_dd = 0; Y=reconstitution(data,N,m,tau); %相空间重构 M=N-(m-1)*tau; %重构相空间中相点的个数 for i = 1 : (M-1) for j = i+1 : M d = 0; for k = 1 : m d = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j)); end d = sqrt(d); if max_d < d max_d = d; end if min_d > d min_d = d; end avg_dd = avg_dd + d; end end avg_d = 2*avg_dd/(M*(M-1)); %平均相点距离 dlt_eps = (avg_d - min_d) * 0.02 ; 对 max_eps 的放宽幅度 min_eps = min_d + dlt_eps / 2 ; max_eps = min_d + 2 * dlt_eps ; % %若在 min_eps~max_eps 中找不到演化相点时, %演化相点与当前相点距离的最小限 %&&演化相点与当前相点距离的最大限

从 P+1~M-1 个相点中找与第一个相点最近的相点位置(Loc_DK)及其最短距离 DK DK = 1.0e+100; %第 i 个相点到其最近距离点的距离 Loc_DK = 2; %第 i 个相点对应的最近距离点的下标 for i = (P+1) M-1) %限制短暂分离,从点 P+1 开始搜索

d = 0; for k = 1 : m d = d + (Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1)); end d = sqrt(d); if (d < DK) & (d > min_eps) DK = d; Loc_DK = i; end end

% 以下计算各相点对应的李氏数保存到 lmd()数组中 % i 为相点序号,从 1 到(M-1),也是 i-1 点的演化点;Loc_DK 为相点 i-1 对应最短距离 的相点位置,DK 为其对应的最短距离 % Loc_DK+1 为 Loc_DK 的演化点,DK1 为 i 点到 Loc_DK+1 点的距离,称为演化距离 % 前 i 个 log2(DK1/DK)的累计和用于求 i 点的 lambda 值 sum_lmd = 0 ; % 存放前 i 个 log2(DK1/DK)的累计和 for i = 2 : (M-1) % 计算演化距离 DK1 = 0; for k = 1 : m DK1 = DK1 + (Y(k,i)-Y(k,Loc_DK+1))*(Y(k,i)-Y(k,Loc_DK+1)); end DK1 = sqrt(DK1); old_Loc_DK = Loc_DK ; % 保存原最近位置相点 old_DK=DK; % 计算前 i 个 log2(DK1/DK)的累计和以及保存 i 点的李氏指数 if (DK1 ~= 0)&( DK ~= 0) sum_lmd = sum_lmd + log(DK1/DK) /log(2); end lmd(i-1) = sum_lmd/(i-1); % 以下寻找 i 点的最短距离:要求距离在指定距离范围内尽量短,与 DK1 的角度最小 point_num = 0 ; % &&在指定距离范围内找到的候选相点的个数 cos_sita = 0 ; %&&夹角余弦的比较初值 ――要求一定是锐角 zjfwcs=0 ;%&&增加范围次数 while (point_num == 0) % * 搜索相点 for j = 1 : (M-1) if abs(j-i) <=(P-1) %&&候选点距当前点太近,跳过! continue; end %*计算候选点与当前点的距离 dnew = 0; for k = 1 : m dnew = dnew + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j)); end dnew = sqrt(dnew); if (dnew < min_eps)|( dnew > max_eps ) %&&不在距离范围,跳过! continue; end %*计算夹角余弦及比较 DOT = 0; for k = 1 : m

DOT = DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_Loc_DK+1)); end CTH = DOT/(dnew*DK1); if acos(CTH) > (3.14151926/4) continue; end %&&不是小于 45 度的角,跳过!

if CTH > cos_sita %&&新夹角小于过去已找到的相点的夹角,保留 cos_sita = CTH; Loc_DK = j; DK = dnew; end point_num = point_num +1; end if point_num <= min_point max_eps = max_eps + dlt_eps; zjfwcs =zjfwcs +1; if zjfwcs > MAX_CISHU %&&超过最大放宽次数,改找最近的点 DK = 1.0e+100; for ii = 1 : (M-1) if abs(i-ii) <= (P-1) %&&候选点距当前点太近,跳过! continue; end d = 0; for k = 1 : m d = d + (Y(k,i)-Y(k,ii))*(Y(k,i)-Y(k,ii)); end d = sqrt(d); if (d < DK) & (d > min_eps) DK = d; Loc_DK = ii; end end break; end point_num = 0 cos_sita = 0; end end end ; %&&扩大距离范围后重新搜索

%取平均得到最大李雅普诺夫指数 lambda_1=sum(lmd)/length(lmd);

程序中用到的 reconstitution 函数如下: function X=reconstitution(data,N,m,tau) %该函数用来重构相空间 % m 为嵌入空间维数 % tau 为时间延迟 % data 为输入时间序列 % N 为时间序列长度 % X 为输出,是 m*n 维矩阵 M=N-(m-1)*tau;%相空间中点的个数 for j=1:M %相空间重构 for i=1:m X(i,j)=data((i-1)*tau+j); end end

这里声明一下,这些程序并非我自己编写的,均是转载,其使用我已经验证过,绝对可以运 行! (3)小数据量方法 说小数据量方法是目前最实用、应用最广泛的方法应该不为过吧,呵呵!

*** 上面两种方法,即 Wolf 方法和小数据量方法,在计算 LE 之前,都要求对时间序列进 行重构相空间,重构相空间的优良对于最大 LE 的计算精度影响非常大!因此重构相空间的 几个参数的确定就非常重要。 (1)时间延迟 主要推荐两种方法――自相关函数法、C-C 方法 自相关函数法――对一个混沌时间序列, 可以先写出其自相关函数, 然后作出自相关函数关 于时间 t 的函数图像。根据数值试验结果,当自相关函数下降到初始值的 1-1/e 时,所得的 时间 t 即为重构相空间的时间延迟。 C-C 方法――可以同时计算出时间延迟和时间窗口,个人推荐使用这种方法! (2)平均周期 平均周期的计算可以采用 FFT 方法。 matlab 帮助中有一个太阳黑子的例子,现摘录如下: 在 load sunspot.dat %装载数据文件 year = sunspot(:,1); %读取年份信息 wolfer = sunspot(:,2); %读取黑子活动数据 plot(year,wolfer) %绘制原始数据图

title('Sunspot Data') Y = fft(wolfer); %快速 FFT 变换

N = length(Y); %FFT 变换后数据长度 Y(1) = []; %去掉 Y 的第一个数据,它是所有数据的和 power = abs(Y(1:N/2)).^2; %求功率谱 nyquist = 1/2; freq = (1:N/2)/(N/2)*nyquist;%求频率 plot(freq,power), grid on %绘制功率谱图 xlabel('cycles/year') title('Periodogram') period = 1./freq; %年份(周期) plot(period,power), axis([0 40 0 2e7]), grid on %绘制年份-功率谱曲线 ylabel('Power') xlabel('Period(Years/Cycle)') [mp,index] = max(power); period(index) %求最高谱线所对应的年份下标 %由下标求出平均周期

(3)嵌入维数 目前嵌入维数的主要计算方法是采用 Grassberger 和 Procaccia 提出的 G-P 算法计算出序列 的关联维数 d,然后利用嵌入维数 m>=2d+1,选取合适的嵌入维数。 G―P 算法程序如下: function [ln_r,ln_C]=G_P(data,N,tau,min_m,max_m,ss) % the function is used to calculate correlation dimention with G-P algorithm % 计算关联维数的 G-P 算法 % data:the time series 时间序列 % N: the length of the time series 时间序列长度 % tau: the time delay 时间延迟 % min_m:the least embedded dimention m 最小的嵌入维数 % max_m:the largest embedded dimention m 最大的嵌入维数 % ss:the stepsize of r r 的步长 %skyhawk for m=min_m:max_m Y=reconstitution(data,N,m,tau);%reconstitute state space M=N-(m-1)*tau;%the number of points in state space for i=1:M-1 for j=i+1:M d(i,j)=max(abs(Y(:,i)-Y(:,j)));%calculate the distance of each two end %points in state space 计算状态空间中每两点之间的距离 end max_d=max(max(d));%the max distance of all points 得到所有点之间的最大距离

d(1,1)=max_d; min_d=min(min(d));%the min distance of all points 得到所有点间的最短距离 delt=(max_d-min_d)/ss;%the stepsize of r 得到 r 的步长 for k=1:ss r=min_d+k*delt; C(k)=correlation_integral(Y,M,r);%calculate the correlation integral ln_C(m,k)=log(C(k));%lnC(r) ln_r(m,k)=log(r);%lnr fprintf('%d/%d/%d/%d\n',k,ss,m,max_m); end plot(ln_r(m, hold on; end fid=fopen('lnr.txt','w'); fprintf(fid,'%6.2f %6.2f\n',ln_r); fclose(fid); fid = fopen('lnC.txt','w'); fprintf(fid,'%6.2f %6.2f\n',ln_C); fclose(fid); 程序中的 correlation_integral 函数如下: function C_I=correlation_integral(X,M,r) %the function is used to calculate correlation integral %C_I:the value of the correlation integral %X:the reconstituted state space,M is a m*M matrix %m:the embedding demention %M:M is the number of embedded points in m-dimensional sapce %r:the radius of the Heaviside function,sigma/2<r<2sigma %calculate the sum of all the values of Heaviside %skyhawk sum_H=0; for i=1:M % fprintf('%d/%d\n',i,M); for j=i+1:M d=norm((X(:,i)-X(:,j)),inf);%calculat the distances of each two points in matris M with sup-norm sita=heaviside(r,d);%calculate the value of the heaviside function sum_H=sum_H+sita; end end C_I=2*sum_H/(M*(M-1));%the value of correlation integral ,ln_C(m, );

以上的各种方法在实际应用的时候要根据具体情况来选择。 一般地,如果已知系统方程(当然系统不能太过复杂)时,则计算 Lyapunov 指数采用定义 法、Jacobian 方法要精确、简单些! 而如果系统方程比较复杂(如超维系统) 、或者为一时间序列时,则推荐采样 Wolf 方法、小 数据量方法。 Wolf 方法的特点是时间序列无噪声,空间中小向量的演变高度非线性,而 Jacobian 方法则 是噪声大,空间中小向量的演变接近线性。 小数据量方法的优点在于: (1)对小数据组的计算可靠; (2)计算量较小,比 wolf 方法快 很多; (3)编程、操作较为容易。

而关于时间延迟、嵌入维数、平均周期的确定,还是推荐 C-C 方法和 G-P 算法,结果更 为可靠一些!


相关文档

更多相关文档

Lyapunov 指数
-Lyapunov指数的计算方法
matlab编写的Lyapunov指数计算程序
Rossler 混沌系统Lyapunov指数Matlab实现代码
不光滑混沌系统的Lyapunov指数计算
电脑版