小数据量方法计算le指数

竹园

计算了实测数据的le指数，开始的时候采用了接近4000个数据运行时，老是提醒参数没定义，后来将数据缩短一半，结果出来了，平均轨道周期24，这和我采取的数据为小时数据基本吻合，但是数据大于4000的时候，算出来的平均轨道周期就是序列长度，不知道是什么原因，数据拟合时选取了前面200个点，最后le值为0.0018，有一些疑问，吕金虎的书中提到le>0就可以说明该系统处于混沌状态，但是我的结果波动比较大，不知道可不可以说明该系统处于混沌状态，请各位高手指点，下面附上我处理结果，谢谢了！

gghhjj

1. 你的结果好奇怪啊,最好附上代码和数据
2. 只要le>1就可以认为系统混沌

竹园

下面是我计算的结果和采用的数据，拟合选130个点，le指数=0.0069，请指教
function lambda_1=largest_lyapunov_exponent(data,m,tau)            
%本函数使用小数据量方法计算最大lyapunov指数
%data:时间序列
%m:嵌入维
%tau:时间延迟
%P:使用 FFT计算出的时间序列平均周期
%lambda_1:函数返回的最大lyapunov指数值
N=length(data);  %序列长度
P=period_mean_fft(data);
P=fix(P)%时间序列的平均周期
delt_t=1;
Y=reconstitution(data,m,tau );  %重构相空间
M=N-(m-1)*tau;  %M是m维重构相空间中总点数
d=zeros(M-1,M);
for j=1:M
    d_min=1e+100;
    for jj=1:M      %寻找相空间中每个点的最近距离点，并记下该点下标              
        if abs(j-jj)>P      %限制短暂分离
            d_s=norm(Y(:,j)-Y(:,jj));%计算分离后两点的距离
            if d_s<d_min
                d_min=d_s;
                idx_j=jj;
            end           
        end
    end
    max_i=min((M-j),(M-idx_j));    %计算点j的最大演化时间步长i
    for k=1:max_i                 %计算点j与其最近邻点在i个离散步后的距离
       d(k,j)=norm(Y(:,j+k)-Y(:,idx_j+k));
    end
end
%对每个演化时间步长i，求所有的j的lnd(i,j)平均
[l_i,l_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
    if q>0
      y(i)=y_s/(q*delt_t);
    end
end
I=1:length(y);
pp=polyfit(I,y,1);
lambda_1=pp(1);
yp=polyval(pp,I);
figure(2);
plot(I,y,'-');
xlabel('i');
ylabel('y(i)');hold on;
for i=1:length(y)-1
    x(i)=y(i+1)-y(i);
end
figure;
plot(I(1:(length(y)-1)),x);
xlabel('i');
ylabel('y(i)-y(i-1)');
hold on;
linear=input('请输入线型部分的长度')
I1=1:linear;
PP=polyfit(I1,y(I1),1);
lambda_1=PP(1)
YP=polyval(PP,I1);
figure(2);
plot(I1,YP,'r-','linewidth',3);


function T_mean=period_mean_fft(data)
%该函数使用快速傅里叶变换FFT计算序列平均周期
%data：时间序列
%T_mean：返回快速傅里叶变换FFT计算出的序列平均周期
Y = fft(data);       %快速FFT变换
N = length(data);    %FFT变换后数据长度
Y(1) = [];           %去掉Y的第一个数据，它是data所有数据的和
power = abs(Y(1:N/2)).^2;  %求功率谱
nyquist = 1/2;
freq = (1:N/2)/(N/2)*nyquist; %求频率
subplot(121)
plot(freq,power); grid on ;    %绘制功率谱图
xlabel('频率');
ylabel('功率');
title('功率谱图');
period = 1./freq;                %计算周期
subplot(122)
plot(period,power); grid on  %绘制周期－功率谱曲线
ylabel('功率');
xlabel('周期');
title('周期—功率谱图');
[mp,index] = max(power);       %求最高谱线所对应的下标
T_mean=period(index)            %由下标求出平均周期


data.txt (13.5 KB, 下载次数: 5)

2012-12-30 13:52 上传

点击文件名下载附件

竹园

gghhjj 发表于 2012-12-27 10:32

                               
登录/注册后可看大图

1. 你的结果好奇怪啊,最好附上代码和数据
2. 只要le>1就可以认为系统混沌

怎么个奇怪发呢，亲指教，下面是我的程序和结果

warden

gghhjj 发表于 2012-12-27 10:32

                               
登录/注册后可看大图

1. 你的结果好奇怪啊,最好附上代码和数据
2. 只要le>1就可以认为系统混沌

是大于1吗？不是说大于0就可以认为系统是混沌的么？