回复 2 # hsfy919 的帖子
该程序如下
function [Texp,Lexp]=lyapunov2(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
n1=n; n2=n1*(n1+1);
nit = round((tend-tstart)/stept);
y=zeros(n2,1); cum=zeros(n1,1); y0=y;
gsc=cum; znorm=cum;
y(1:n)=ystart(:);
for i=1:n1 y((n1+1)*i)=1.0;
end;
t=tstart;
for ITERLYAP=1:nit
[T,Y] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y);
t=t+stept;
y=Y(size(Y,1),:);
for i=1:n1
for j=1:n1 y0(n1*i+j)=y(n1*j+i); end;
end;
znorm(1)=0.0;
for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)^2; end;
znorm(1)=sqrt(znorm(1));
for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end;
for j=2:n1
for k=1:(j-1)
gsc(k)=0.0;
for l=1:n1
gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k);
end;
end;
for k=1:n1
for l=1:(j-1)
y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);
end;
end;
znorm(j)=0.0;
for k=1:n1
znorm(j)=znorm(j)+y0(n1*k+j)^2;
end;
znorm(j)=sqrt(znorm(j));
for k=1:n1
y0(n1*k+j)=y0(n1*k+j)/znorm(j);
end;
end;
for k=1:n1
cum(k)=cum(k)+log(znorm(k));
end;
for k=1:n1
lp(k)=cum(k)/(t-tstart);
end;
if ITERLYAP==1
Lexp=lp;
Texp=t;
else
Lexp=[Lexp; lp];
Texp=[Texp; t];
end;
for i=1:n1
for j=1:n1
y(n1*j+i)=y0(n1*i+j);
end;
end;
end;
使用的时候,只要定义好函数,再调用以上文件就行了,不知道能不能改成离散系统使用的。
先谢谢大家啦^_^ |