molshibi123 发表于 2009-10-20 09:01

用wolf法求解最大李氏指数请进来讨论

小弟找这个wolf最大李氏指数的程序已经好久了,一直没找到准确的,请高手指点下问题出在哪里:下面是我调用的程序,其中计算时间延迟和嵌入维数、周期都是准确的,我已经和吕金虎的书上对照过,但是调用wolf时怎么就不对呢请高手指点。
C-C法算时间延迟和嵌入维数:
clear all
clear all

=ode45(@LorenzDifEqn1,,);
xx1 = y1(:,1)';
delt_t1 = 0.01;
t1 = 0:delt_t1:75;
x = spline(tt1, xx1, t1);
X= x(1000:1:3999);

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);



function dx = LorenzDifEqn1(t,x);
% Lorenz系统微分方程
% 方程如下:
% dx = -σ*(x-y)
% dy = -x*z+r*x-y
% dz = x*y-b*z
dx_1 =35*(x(2)-x(1));
dx_2 = -7*x(1)-x(1)*x(3)+28*x(2);
dx_3 =x(1)*x(2)-3*x(3);
dx = ;


function =C_CMethod_inf(X,max_d)
% 用于求延迟时间tau
% X为输入时间序列
% max_d为最大时间延迟
% Smean,Sdeltmean,Scor为返回值
% tau为计算得到的延迟时间
% tw为时间窗口
% zhangli
% 2008-11-30

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




function data_d=disjoint(data,N,t)
% 此函数用于将时间序列分解成t个不相交的时间序列
% data:输入时间序列
% N:data的长度
% t:the index lag
% data_d:返回分解后的t个不相交的时间序列
% 2008-11-28
% zhangli
for i=1:t
    for j=1:(N/t)
      data_d(i,j)=data(i+(j-1)*t);
    end
end



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

% M=N-(m-1)*tau;                                  % 相空间中点的个数
% Data=zeros(M,m);
% for i=1:M                                    
%   for j=1:m
%         Data(i,j)=data(i+(j-1)*tau);             % 相空间重构
%   end
% end
% M=N-(m-1)*tau;                                  % 相空间中点的个数
% Data=zeros(M,m);
% for i=1:M                                    
%   for j=1:m
%         Data(i,j)=data(j+(i-1)*tau);             % 相空间重构
%   end
% end

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





function C=correlation_integral_inf(Y,M,r)
% 此函数用于计算关联积分,取无穷范数
% Y为重构的相空间
% M为相空间中点的个数
% r为搜索半径
% Y为输出,是M*m维矩阵
% 2008-11-30
% zhangli
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));


以上是Chen系统的时间延迟和嵌入维数的程序,运行得出的结果是时间延迟:9,时间窗口是68,计算得出嵌入维数8,和书上很接近。可以说是准确的。


以下是计算周期的




clear all
clear all

=ode45(@LorenzDifEqn1,,);
xx1 = y1(:,1)';
delt_t1 = 0.01;
t1 = 0:delt_t1:75;
x = spline(tt1, xx1, t1);
X= x(1000:1:3999);
%function TT=FFTcircle(X,N)%时间序列平均周期计算
Y=fft(X);
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/tongjishijianjiange')
title('Periodogram')
period = 1./freq;                     %年份(周期)
plot(period,power); % , axis(), grid on%绘制年份-功率谱曲线
ylabel('Power')
xlabel('Period(Years/Cycle)')
= max(power);       %求最高谱线所对应的年份下标
%求时间序列平均周期TT
TT=period(index);   


运行结果是60,书上是60.2













问题就是当我用以上的结果调用wolf时结果就不对了呢?请高手指教:


clc
clear all
close all
%--------------------------------------------------------------------------
% Henon映射
% x(n+1) = 1 - a * x(n)^2 + y(n);
% y(n+1) = b * x(n)
%a = 1.4;
%b = 0.3;
%x0 = 0.25;
%y0 = 0.25;
%k1 = 1000;                  % 前面的迭代点数
%k2 = 4000;                % 后面的迭代点数
%Z = zeros(k1+k2,2);
%for i = 1:k1+k2
%   x = 1 - a * x0^2 + y0 ;
%   y = b * x0;
%x0 = x;
%y0 = y;
   
   % Z(i,1) = x;
   % Z(i,2) = y;
%end
%X = Z(k1+1:end,1);
=ode45(@LorenzDifEqn1,,);
xx1 = y1(:,1)';
delt_t1 = 0.01;
t1 = 0:delt_t1:75;
x = spline(tt1, xx1, t1);
X= x(1000:1:3999);

N=length(X);
lambda_1=lyapunov_wolf(X,N,8,10,60)
%lambda_1=lyapunov_wolf(X,N,8,10,60)




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);






以上程序是别人所编,自己只是改了下系统,请高手指教。我这样调用chen系统的最大李氏指数是0..02了,结果明显不对,高手帮我看看,
附:我认为C-C法计算时间延迟和嵌入维数是对的,需要的可以下载用哦

molshibi123 发表于 2009-10-20 09:27

octopussheng 帮帮我吧,看该怎么改啊

octopussheng 发表于 2009-10-20 09:34

因为wolf法我用的不多,这里提一点小建议:wolf法中各个参数设置都要正确,不仅仅包括时间延迟和嵌入维数,平均周期也是关键的。此外,你的数据量同样重要,最好进行重新采样,再计算时间延迟、嵌入维数、平均周期等

molshibi123 发表于 2009-10-20 10:14

要是时间延迟和嵌入维数,平均周期都是准确的话,向这样调用wolf可以不?我就是怕这个wolf程序是错的,但我看到很多网站都是这个程序,

还有个问题,就是重新采样怎么理解?就是不能用ode45算出来的原始数据吗

molshibi123 发表于 2009-10-20 16:12

怎么没人理我啊,郁闷了几个月了

octopussheng 发表于 2009-10-20 16:54

为什么一定要wolf方法计算呢?
个人觉得小数据量方法要准确一些。

molshibi123 发表于 2009-10-21 08:08

回复 6楼 octopussheng 的帖子

你有小数据法算最大李氏指数的程序吗?我在论坛上下的小数据法求的也不对啊,尤其里面还有个算最大离散时间步进时间、线性长度我都不知该如何取值呢。 octopussheng你能发一份计算准确的小数据法程序吗?非常感谢,我的邮箱是:molshibi@163.com

octopussheng 发表于 2009-10-21 08:22

本版里面有不少讨论的,你可以搜索搜索看看。

molshibi123 发表于 2009-10-21 09:08

哦,论坛上我看多了,但没有一个算的准的,反正都是能运行,结果差了10万八千

molshibi123 发表于 2009-10-21 10:59

还有就是想问下octopussheng ,在用小数据量法时线性区域怎么选取呢,我看到讨论中有说到下界是(m-1)tau+1,那上界呢?看了小数据法上讲的是一个区域啊,还有就是最大离散步进时间怎么取的呢

lushuang8 发表于 2009-10-21 15:58

回复

版主给出的程序不能运行啊!!!!
>> chen
???In an assignment A(matrix,:) = B, the number of elements in the subscript of A and
the number of rows in B must be the same.

Error in ==> F:\wolf\C_CMethod_inf.m (reconstitution)
On line 106==>   Data(i,:)=data([((i-1)*tau+1):1:((i-1)*tau+M)]);

Error in ==> F:\wolf\C_CMethod_inf.m
On line 30==>               Z=reconstitution(Y,N_t,m,1);             % 相空间重构

Error in ==> F:\wolf\chen.m
On line 13==> =C_CMethod_inf(X,max_d);

>>

lushuang8 发表于 2009-10-21 16:07

回复

Henon映射:
lambda_1 =

    0.0075

>>

molshibi123 发表于 2009-10-21 17:06

能运行的,我都运行过,你用的matlab是哪个版本的啊

lushuang8 发表于 2009-10-21 19:20

回复

我的是6.5的

octopussheng 发表于 2009-10-21 21:44

小数据量法确实没怎么用,如果你不是很急的话,我近期可以帮你一起看看,如果很急的话,我也没办法了,呵呵,实在抱歉。
页: [1] 2 3
查看完整版本: 用wolf法求解最大李氏指数请进来讨论