我的小数据量法计算Lyapunov指数的程序，大家看一下有没有问题

anqikeli

function lambda_1=Lyapunov(data,N,m,tau,P,delt_t)
clear
clc
%获得时间序列
[tt1,y1]=ode45(@Lorenz,[0,200],[0.1,0.1,0.1]);
xx1 = y1(:,1)';
delt_t = 0.01;
t1 = 0:delt_t:75;
x1 = spline(tt1, xx1, t1);
data = x1(1000:1:6000);
N=5001;
m=3;         %嵌入维数
tau=11;       %时延

P=ave_T(data,N);     %求平均周期
delt_t=0.01;
d_length=[];
d_content=[];
M=N-(m-1)*tau;
Y=reconstitution(data,M,m,tau);
idx_j=0;
for j=1:M
    d_min=10000;
    for jj=1:M                                             
        d_s=0;      %寻找相空间中每个点的最近距离点，并记下该点下标
        if abs(j-jj)>P                                      %限制短暂分离
            for i=1:1:m
                %d_s=d_s+abs(Y(j*M+i)-Y(jj*M+i));
                d_s=norm(abs(Y(j,:)-Y(jj,:)));
            end            
            if d_s<d_min
                d_min=d_s;
               idx_j=jj;
            end
        end
    end
max_i=min((M-j),(M-idx_j))
    d_length=[d_length,max_i];
    for k=1:max_i              %计算点j与其最近邻点在i个离散步后的距离
        d_j_i=0;        
        d_j_i=norm((Y(j+k,:)-Y(idx_j+k,:)),inf);
        d_content=[d_content,d_j_i];
    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=[y,S_j_i/(Count*delt_t)]; %对每个演化时间步长i，求所有的j平均
end
XX=1:length(y);
plot(XX,y(XX),'.','markersize',1);hold on;
YY=1:300;
plot(YY,y(YY),'.','markersize',1);%前300个点画图
linear=200；
XX1=1:linear;
pp=polyfit(XX1,y(XX1),1);  %前200个点拟合
lambda_1=pp(1)
yp=polyval(pp,XX1);
plot(XX1,yp,'r--');

function P=ave_T(x,N)
Y=fft(x);
N=length(Y);
L=round(N/2);%四舍五入取整
Y(1)=[];
power=abs(Y(1:L)).^2;%%幅度的平方为能量
nyquist=1/2;
freq=(1:L)/(L)*nyquist;
%subplot(2,1,1);plot(freq,power);
%title('freq-power gram');
period=1./freq;
%subplot(2,1,2);plot(period,power);
%title('period-power gram');
[mp,index]=max(power);
P=period(index);
return

function X=reconstitution(data,M,m,t)
%该函数用来重构相空间
% 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
xv=zeros(m,M);
for ii=1:m
    X1(ii,:)=data([((ii-1)*t+1):1:((ii-1)*t+M)]);
end
X=X1';
return
这是我根据论坛上的一些帖子整理出来的程序，运行的很慢，我估计是在哪里存在错误，但又实在是找不到，请各位高手帮忙看一下，不胜感激！

yufeiqun2008

语句d_content=[d_content,d_j_i];处于最内层循环,随着j的增大，d_content越来越大，当长度超过100000时(此时估计j还远小于M)，该语句执行耗费时间将非常长，建议修改该语句

无穷范数也最好换作max语句

[ 本帖最后由 yufeiqun2008 于 2008-1-8 13:03 编辑 ]

anqikeli

那我把d_count换成是一个M*（M-1）的矩阵，这样是不是会快呢，不过会浪费一些空间，你怎么做的啊？指点一下吧

anqikeli

还有个问题，限制短暂分离那里，会有jj取到M的情况，此时max_i=min(M-j,M-jj)就为0了，这个要怎么处理？？

octopussheng

能否考虑让jj只取到M-1呢？

yufeiqun2008

限制短暂分离,你的程序代码运行效率很低,我已经给你发了一份我做的程序,你可以看看,没有收到吗?

anqikeli

没有收到，再给我发一份好么？
153201909@qq.com
或anqikeli@126.com都可以，太感谢了，我是真的不知道怎么改了

gurencai

大哥，你能不能也发一份给我啊？万分感激你！
愿菩萨保佑好人一生平安！
:@) 我的邮箱是rencai2050@163.com

[ 本帖最后由 gurencai 于 2008-5-7 21:30 编辑 ]

swayfish

copuk321@163.com
能给我也发一份吗？谢谢！

雪山飞鸟

非常好，学习学习~~~~~~~~

165390292

yufeiqun2008 发表于 2008-1-8 23:35

                               
登录/注册后可看大图

限制短暂分离,你的程序代码运行效率很低,我已经给你发了一份我做的程序,你可以看看,没有收到吗?

能给我发一份吗？十分感谢！165390292@qq.com