论坛 › 数学与非线性科学 › 关于小数据量法求混沌lyapunov指数的问题

darcykuzy 发表于 2012-1-1 16:56

关于小数据量法求混沌lyapunov指数的问题

本帖最后由 darcykuzy 于 2012-1-1 17:11 编辑

我从网上找了2个小数据量法求最大lyapunov指数的程序，但运行后得到两个完全不同的结果，一个为正，一个为负，不确定哪个程序是正确的，恳请各位帮忙看看，相应程序及结果如下：
（1）function lambda_1=largest_lyapunov_exponent(data,N,m,tau,P)
clc
clear all
%the function is used to calcultate largest lyapunov exponent with the
%mended algorithm,which put forward by lv jing hu.
%data:the time series
%N:the length of data
%m:enbedding dimention
%tau:time delay
%P:the mean period of the time series,calculated with FFT
%lambda_1:return the largest lyapunov exponent
%skyhawk

data=load('E:\ma\an.txt');
N=length(data);
m=9;
tau=4;
P= 52; %%% or P=1;
delt_t=1;


Y=reconstitution(data,N,m,tau );%reconstitute state space
M=N-(m-1)*tau;%M is the number of embedded points in m-dimensional space
for j=1:M
    d_max=1e+100;
    for jj=1:M                                              %寻找相空间中每个点的最近距离点，并记下
      d_s=0;                                              %该点下标
      if abs(j-jj)>P                                    %限制短暂分离
            for i=1:m
                d_s=d_s+(Y(i,j)-Y(i,jj))*(Y(i,j)-Y(i,jj));
                d_min=d_max;
                if d_s<d_min
                   d_min=d_s;
                   idx_j=jj;
               end
            end
      end
    end
%   index(j)=idx_j;
    max_i=min((M-j),(M-idx_j));%计算点j的最大演化时间步长i
    for k=1:max_i            %计算点j与其最近邻点在i个离散步后的距离
      d_j_i=0;
      for kk=1:m
            d_j_i=d_j_i+(Y(kk,j+k)-Y(kk,idx_j+k))*(Y(kk,j+k)-Y(kk,idx_j+k));
            d(k,j)=d_j_i;
      end
    end
end

%对每个演化时间步长i，求所有的j的lnd(i,j)平均
=size(d);
for i=1:l_i
    q=0;
    y_s=0;
    for j=1:l_j
      if d(i,j)~=0
            q=q+1;
            y_s=y_s+log(d(i,j));
      end
    end
    y(i)=y_s/(q*delt_t)
end
x=1:length(y);
pp=polyfit(x,y,1);
lambda_1=pp(1);
yp=polyval(pp,x);
plot(x,y,'-o',x,yp,'--')

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

运行后的结果为：0.0058

（2）function lambda_1=largest_lyapunov_p3(data,N,m,tau,P,delt_t)
clc
clear all
m=9;tau=4;
delt_t=1;
data=load('E:/ma/an.txt');
N=length(data);
P=1;%%%平均周期
d_length=[];
d_content=[];
Y=reconstitution(data,N,m,tau);
M=N-(m-1)*tau;
idx_j=0;
for j=1:M
    d_min=1000;
    for jj=1:M                                             
      d_s=0;      %寻找相空间中每个点的最近距离点，并记下该点下标
      if abs(j-jj)>P                                    %限制短暂分离
            d_s=sum(abs(Y(j)-Y(jj)));
            if d_s<d_min
                d_min=d_s;
               idx_j=jj;
            end
      end
    end
    if ((M-j)>(M-idx_j));%计算点j的最大演化时间步长i
      max_i=M-idx_j;
    else
      max_i=M-j;
    end
    d_length=;
    for k=1:max_i            %计算点j与其最近邻点在i个离散步后的距离
      d_j_i=0;      
      d_j_i=norm(Y(j+k)-Y(idx_j+k));
      d_content=;
    end
end
%对每个演化时间步长i，求所有的j的lnd(i,j)平均
y=[];
for i=1:max(d_length)
   S_j_i=0;
   sum_former=0;
   Count=0;
   for j=1:M
       if j==1
          former=0;
       else
          former=d_length(j-1);
       end   
       sum_former=sum_former+former;
       if i<=(d_length(j))
         if d_content(sum_former+i)~=0
            S_j_i=S_j_i+log(d_content(sum_former+i));
            Count=Count+1;
         end
       end
   end
   y=; %对每个演化时间步长i，求所有的j平均
end
XX=1:length(y);
figure;
plot(XX,y,'-','markersize',1);
hold on;
linear=input('请输入线形部分的长度');
XX1=1:linear;
pp=polyfit(XX1,y(XX1),1);
lambda_1=pp(1)
yp=polyval(pp,XX1);
figure;
plot(XX1,yp,'r--')

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
运行后得到的图为：

输入线性部分长度后（100:150，自己随意取的）得到的结果及图形为：-0.0040

darcykuzy 发表于 2012-1-1 17:26

问题补充：
（1）所得的图是一个周期的，且直线的斜率是对全部数据进行最小二乘拟合，我看一些资料里说的是只求直线部分的斜率，是不是得修改相应代码？
（2）中所得图形包含所有数据，但与(1)所得图形相比，应如何确定线性部分?且（1）（2）中y轴值正、负不同的原因是什么？
问题比较弱，希望各位前辈不要见笑，谢谢！

darcykuzy 发表于 2012-1-4 11:06

等的好辛苦啊！！各位大侠赶紧现身吧！{:3_56:}
{:{13}:}

live101 发表于 2012-1-6 10:53

我也在学习这一部分，有很多问题也不懂。斜率是取前面线性部分的

后来1206 发表于 2012-5-13 19:32

不知道怎么联系楼主。我的邮箱czf1206@163.com希望向楼主学习

165390292 发表于 2012-9-26 19:42

这个问题我也想知道！同求。

liguangzhigong 发表于 2014-5-16 23:43

自然是第二种对了。但是你的线性区选择的不太准确。但看起来叶差不多。第一种方法没有结束。只要你根据想吐找出线性区再继续算就可以了。呵呵。

原来是小飞 发表于 2014-5-19 22:04

liguangzhigong 发表于 2014-5-16 23:43
自然是第二种对了。但是你的线性区选择的不太准确。但看起来叶差不多。第一种方法没有结束。只要你根据想吐 ...

用第二种方法，当数据量很大的时候（20000多），计算时间是不是很长啊？

liguangzhigong 发表于 2014-5-20 06:28

现在有计算机的高运算速度，大家更加关注的是计算的精度问题，计算时间还可以吧。我用的是C_C方法求tau和关联维，数据选取5000点左右。计算时间25分钟左右吧。毕竟是小数据的方法，从20000点中挑出5000点计算就可以了。我是从120000的数据中挑的5000点。呵呵

原来是小飞 发表于 2014-6-3 16:17

liguangzhigong 发表于 2014-5-20 06:28
现在有计算机的高运算速度，大家更加关注的是计算的精度问题，计算时间还可以吧。我用的是C_C方法求tau和关 ...

为什么我用第二种方法计算4253个数据，使用7.1和7.6版本，计算时间会很长，出不来结果呢？

liguangzhigong 发表于 2014-6-13 22:03

你说的“使用7.1和7.6版本”是matlab的版本吗？你所说的很长时间是多长？超过30分钟吗？我认为我的方法已经很耗时间了，也就是大概25分钟左右。出不来结果是什么意思?有错误提示吗？

原来是小飞 发表于 2014-6-14 10:47

liguangzhigong 发表于 2014-6-13 22:03
你说的“使用7.1和7.6版本”是matlab的版本吗？你所说的很长时间是多长？超过30分钟吗？我认为我的方法已经 ...

“使用7.1和7.6版本是指的matlab版本。出不来结果指的是一直运行一个小时左右都没有结果，一直运行中，不知道咋回事。

原来是小飞 发表于 2014-6-14 10:48

liguangzhigong 发表于 2014-6-13 22:03
你说的“使用7.1和7.6版本”是matlab的版本吗？你所说的很长时间是多长？超过30分钟吗？我认为我的方法已经 ...

“使用7.1和7.6版本是指的matlab版本。出不来结果指的是一直运行一个小时左右都没有结果，一直运行中，不知道咋回事。

liguangzhigong 发表于 2014-7-17 18:38

对不起，长时间没来了，你说的问题。你说取4253个点，为什么非要选取这个数，一般取3000或者5000最好。因为在相空间重构中，你取4253这个数会出现好多个多余的点，就会出现错误，即便不会提示，在运行过程中也会出现问题，毕竟我们编的程序都是短平快的方式编的，不像官方程序那样考虑得比较全面。

iwanna2008 发表于 2014-7-22 13:14

能不能把数据上传上来呢 分享下

查看完整版本: 关于小数据量法求混沌lyapunov指数的问题