求助：小数据量方法求解le指数

竹园

小数据量方法计算le指数，数据采用lorenz序列10001-13000样本区间，cc法算出m=21，tau=10，le指数算出来是0.0122，但是在使用这个程序算其他的数据却得不到结果，望路过的大侠指教，我是刚刚动手学编程，程序参考了yuling的贴，自己作了修改，data0是我的数据，就是算不出le指数，在运行的时候总是出现没定义

主程序
clear all
t0=0;
tf=130;
[t,x]=ode45(@Lorenz,[t0:0.01:tf],[-1,0,1]);
Lorenz_data=x(10002:end,1);
m=21;
tau=10;
tic
lambda_1=largest_lyapunov_exponent(Lorenz_data,m,tau,t);
toc


function dy=Lorenz(t,y)
%方程的定义
a=16;
b=4.0;
c=45.92;
dy=zeros(3,1);
dy(1)=-a*(y(1)-y(2));
dy(2)=-y(1)*y(3)+c*y(1)-y(2);
dy(3)=y(1)*y(2)-b*y(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)            %由下标求出平均周期

function lambda_1=largest_lyapunov_exponent(data,m,tau,t)   
%本函数使用小数据量方法计算最大lyapunov指数
%data:时间序列
%m:嵌入维
%tau:时间延迟
%P:使用 FFT计算出的时间序列平均周期
%t:抽样时间间隔
%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(1)
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(2)
plot(I(1:(length(y)-1)),x);grid;
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);
plot(I1,YP,'r--');
data0.txt (40.32 KB, 下载次数: 9)

2012-12-20 11:51 上传

点击文件名下载附件
附件

Lorenz

lorenz