以下是我自己编的程序,可能运行速度不是很理想,但经过验证,应该不存在问题,仅供参考!
function [I,x]=biggest_lyapunov_P2(data,m,tau,t)
%data代表时间序列
%m代表嵌入维数
%tau代表时间延迟
%t代表抽样时间间隔
N=length(data);
M=N-(m-1)*tau;
P=cycle(data);
P=fix(P);%时间序列的平均周期
X=reconstitution(data,M,m,tau); %相空间重构
di=zeros(M-1,M-1);
for i=1:M-1
d_min=1e+100;
for j=1:M
d=0;
if abs(i-j)>P
for k=1:m
d=d+(X(i,k)-X(j,k))^2;
end
d=sqrt(d);
if d<d_min
d_min=d;
DK=j;
end
end
end %找出相空间中每个点X(t)的最近邻点
dt=zeros(min(M-i,M-DK));
for ii=1:min(M-i,M-DK)
for kk=1:m
dt(ii)=dt(ii)+(X(i+ii,kk)-X(DK+ii,kk))^2;
end
dt(ii)=sqrt(dt(ii));
di(i,ii)=dt(ii);
end
end %对相空间中每个点X(t),计算出该邻域点对的i个离散时间步后的距离
jj=0;
for l=1:M-1
q=0;
ln_d=0;
for n=1:M-1
if di(n,l)~=0
ln_d=ln_d+log(di(n,l));
q=q+1;
end
end
if ln_d~=0
jj=jj+1;
x(jj)=ln_d/(q*t);
end
end %对每个i求出所有t的lnd(i)平均x(i)
I=1:length(x);
subplot(211)
plot(I,x);grid;
xlabel('i');
ylabel('y(i)');
hold on;
for i=1:length(x)-1
y(i)=x(i+1)-x(i);
end
subplot(212)
plot(I(1:(length(x)-1)),y);grid;
xlabel('i');
ylabel('y(i)-y(i-1)');
hold on
你可以用它试试!
下面是我用它做Lorenz系统得到的图形! |