laoer180203 发表于 2010-3-31 22:35

李雅普诺夫指数计算的问题

选择了以下程序进行试计算,为了验证程序是否正确,我采用的测试数据是论坛上所提及的Lorenz方程算出的:

1.主程序
% 此程序用来测试CC_method



clear all
clear all

%利用方程获得
% 产生 Lorenz 时间序列
% dx/dt = sigma*(y-x)
% dy/dt = r*x - y - x*z
% dz/dt = -b*z + x*y
sigma=16;
% Lorenz 方程参数
b=4;

r=45.92;


y=[-1,0,1];
% 起始点 (1 x 3 的行向量)
h=0.01;
% 积分时间步长
k1=10000;
% 前面的迭代点数
k2=3000;
% 后面的迭代点数
Z=LorenzData(y,h,k1+k2,sigma,r,b);
X=Z(k1+1:end,1);


max_d=200;
% 最大延迟时间

% 调用C_CMethod_inf,求tau
tic
=C_CMethod_inf(X,max_d);
toc
tau_inf
tw_inf
% 相关作图
figure('name','CC法求时间延迟');
plot(1:max_d,Smean_inf,'-b');hold on;
plot(1:max_d,Sdeltmean_inf,'-*c');hold on;
plot(1:max_d,Scor_inf,'-m');hold on;
plot(1:max_d,zeros(1,max_d),'r');
title('C_CMethod_inf');xlabel('Lag');
legend('S(t)平均值','ΔS(t)平均值','Scor_inf');
% 将数据保持下来
fid=fopen('Smean_inf.txt','w');
fprintf(fid,'%f\n',Smean_inf);
fclose(fid);
fid=fopen('Sdeltmean_inf.txt','w');
fprintf(fid,'%f\n',Sdeltmean_inf);
fclose(fid);
fid=fopen('Scor_inf.txt','w');
fprintf(fid,'%f\n',Scor_inf);
fclose(fid);

2.子函数1
function =C_CMethod_inf(X,max_d)
% 用于求延迟时间tau
% X为输入时间序列
% max_d为最大时间延迟
% Smean,Sdeltmean,Scor为返回值
% tau为计算得到的延迟时间
% tw为时间窗口


N=length(X);
Smean=zeros(1,max_d);
Scmean=zeros(1,max_d);
Scor=zeros(1,max_d);
delt=std(X);

% 计算Smean,Sdeltmean,Scor
for t=1:max_d

S=zeros(4,4);

Sdelt=zeros(1,4);

for m=2:5

for j=1:4

r=delt*j/2;

Xdt=disjoint(X,N,t);

% 将时间序列X分解成t个不相交的时间序列

Xdt=Xdt';

s=0;

for tau=1:t

N_t=floor(N/t);
% 分成的子序列长度

Y=Xdt(:,tau);
% 每个子序列


Cs1(tau)=correlation_integral_inf(Y,N_t,r);% 计算C(1,N/t,r,t)


Z=reconstitution(Y,N_t,m,1);
% 相空间重构

Z=Z';

M=N_t-(m-1);

Cs(tau)=correlation_integral_inf(Z,M,r); % 计算C(m,N/t,r,t)

s=s+(Cs(tau)-Cs1(tau)^m);
% 对t个不相关的时间序列求和

end


S(m-1,j)=s/tau;


end

Sdelt(m-1)=max(S(m-1,:))-min(S(m-1,:));
% 差量计算

end

Smean(t)=mean(mean(S));
% 计算平均值

Sdeltmean(t)=mean(Sdelt);
% 计算平均值

Scor(t)=abs(Smean(t))+Sdeltmean(t);
end

% 寻找时间延迟tau:即Sdeltmean第一个极小值点对应的t
for i=2:length(Sdeltmean)-1

if Sdeltmean(i)<Sdeltmean(i-1)&Sdeltmean(i)<Sdeltmean(i+1)

tau=i;

break;

end
end
% 寻找时间窗口tw:即Scor最小值对应的t
for i=1:length(Scor)

if Scor(i)==min(Scor)

tw=i;

break;

end
end

3.子函数2
function data_d=disjoint(data,N,t)
% 此函数用于将时间序列分解成t个不相交的时间序列
% data:输入时间序列
% N:data的长度
% t:the index lag
% data_d:返回分解后的t个不相交的时间序列


for i=1:t

for j=1:(N/t)

data_d(i,j)=data(i+(j-1)*t);

end
end

4.子函数3
function Data=reconstitution(data,N,m,tau)
% 该函数用来重构相空间
% data为输入时间序列
% N为时间序列长度
% m为嵌入空间维数
% tau为时间延迟
% Y为输出,是M*m维矩阵


M=N-(m-1)*tau;
Data=zeros(m,M);
for i=1:m

Data(i,:)=data([((i-1)*tau+1):1:((i-1)*tau+M)]);
end

5.子函数4
function C=correlation_integral_inf(Y,M,r)
% 此函数用于计算关联积分,取无穷范数
% Y为重构的相空间
% M为相空间中点的个数
% r为搜索半径
% Y为输出,是M*m维矩阵


C=0;
for i=1:M-1

for j=i+1:M

d1=norm((Y(i,:)-Y(j,:)),inf);
% 计算状态空间中每两点之间的距离,取无穷范数

if r-d1>=0

C=C+1;


end

end
end
C=2*C/(M*(M-1));

6. 运行结果
tau=10,与书上的值吻合,但是tw算出来是148,而书上是100,通过对scor值的比较发现,确实是148处的值最小(0.0070341
),而在100时scor数值为0.016844,不知为何书上选择100处作为最小值?
另外85点处值为0.013411

7.对指数进行计算
由上cc法算出的嵌入维度是15,再进行lypunove指数计算,
lambda_1=lyapunov_wolf(X,3000,15,tau,20)


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; %&&要求最少搜索到的点数
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
               
                %*计算候选点与当前点的距离
                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);
算出的值是0.0084与A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, Determining Lyapunov exponents from a time series, Physica D 16 (1985) 285.上的1.5相差很远,不知何故,期望指点。
另外对于平均时间周期P的选值,究竟应该选择多少才更为合适,我试了一下,用P=20 和P=60的值基本一致。
再有就是gp法算嵌入维的话,在ln_r,ln_C曲线图上存在直线段,则该直线的斜率即为相关维度D。
反复以不同的嵌入维度m,计算各嵌入维度下之相关维度,即在不同嵌入维度下各有一条曲线与相关维度。若嵌入维度越大,且相关维度成收敛情形,即代表序列具有混沌现象,其收敛值即为此时间序列之相关维度。
究竟是怎么得出 嵌入维的?请大家指教

pestp 发表于 2010-4-2 09:26

好帖,详细,先顶了再说。我是初学的,要是我能用这个计算出来指数,再来谢你

pestp 发表于 2010-4-2 09:40

不会弄。。

laoer180203 发表于 2010-4-2 10:11

回复 沙发 pestp 的帖子

关键是我计算以后发现算出的数并不吻合啊

pestp 发表于 2010-4-4 21:20

同样。。计算出来的lorenz系统(a=16 b=4 c=45.92) 最大指数算出来太小了。。搞不懂

pestp 发表于 2010-4-9 10:48

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);
这个计算余弦值的 怎么跟余弦公式不一样。。是不是错了
cth=(b^2+c^2-a^2)/(2bc)
             CTH=( (Dnew))^2+(DK1)^2-(Y(k,j)-Y(k,old_Loc_Dk+1)^2 ) /(2Dnew*DK1)

wyk_1216 发表于 2010-4-15 15:19

1

Lorenzdata.m文件的内容能不能告诉我一下啊谢谢

yiqguan 发表于 2010-4-15 19:46

有点难:lol

guo2010 发表于 2010-4-22 11:10

很难,不明白% 相关作图
figure('name','CC法求时间延迟');
plot(1:max_d,Smean_inf,'-b');hold on;
plot(1:max_d,Sdeltmean_inf,'-*c');hold on;
plot(1:max_d,Scor_inf,'-m');hold on;
plot(1:max_d,zeros(1,max_d),'r');
title('C_CMethod_inf');xlabel('Lag');
legend('S(t)平均值','ΔS(t)平均值','Scor_inf');
% 将数据保持下来
fid=fopen('Smean_inf.txt','w');
fprintf(fid,'%f\n',Smean_inf);
fclose(fid);
fid=fopen('Sdeltmean_inf.txt','w');
fprintf(fid,'%f\n',Sdeltmean_inf);
fclose(fid);
fid=fopen('Scor_inf.txt','w');
fprintf(fid,'%f\n',Scor_inf);

mning 发表于 2010-5-7 19:57

回复 9楼 guo2010 的帖子

李雅普诺夫指数计算中P是什么意思?

cqupenghao 发表于 2010-8-21 16:06

回复 mning 的帖子


   P是 时间序列的平均轨道周期
页: [1]
查看完整版本: 李雅普诺夫指数计算的问题