马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
本帖最后由 hit126 于 2016-4-22 11:39 编辑
近期为了把计算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),
- [indent] ;
- % 重新定义起始时刻
- 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 ; %若在min_eps～max_eps中找不到演化相点时，对max_eps的放宽幅度
- min_eps = min_d + dlt_eps / 2 ; %演化相点与当前相点距离的最小限
- max_eps = min_d + 2 * dlt_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
- [/indent]
- %*计算候选点与当前点的距离
- 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) %&&不是小于45度的角，跳过!
- continue;
- end
- 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,
- ,ln_C(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
- %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 以上的各种方法在实际应用的时候要根据具体情况来选择。 一般地，如果已知系统方程(当然系统不能太过复杂)时，则计算Lyapunov指数采用定义法、Jacobian方法要精确、简单些! 而如果系统方程比较复杂(如超维系统)、或者为一时间序列时，则推荐采样Wolf方法、小数据量方法。 Wolf方法的特点是时间序列无噪声，空间中小向量的演变高度非线性，而Jacobian方法则是噪声大，空间中小向量的演变接近线性。 小数据量方法的优点在于：(1)对小数据组的计算可靠;(2)计算量较小，比wolf方法快很多;(3)编程、操作较为容易。
而关于时间延迟、嵌入维数、平均周期的确定，还是推荐C-C方法和G-P算法，结果更为可靠一些!
|
发表于 2016-4-22 11:37
提升卡
变色卡
显身卡