请大家帮忙看看这个两个求最大Lyapunov指数的程序
有两种方法，都是根据指数工具箱来的，就是第二个降了一维，但是得到的结果就不一样！function OUT=duffing(t,X)
k=0.5;
f=2.804;
wd=2*pi;
%Rearrange input data in desired format
%Note: the input data is a column vector
x=X(1); y=X(2);z=X(3);
Q=[X(4),X(7),X(10);
X(5),X(8),X(11);
X(6),X(9),X(12)];
%Duffing's equation
dx=wd*y;
dy=wd*(-k*y+x-x^3+f*cos(wd*z));
dz=wd*1;%where z = t, this transformation is for changing
%the non-autonomous system to a autonomous one
DX1=; %Output data
%Linearized system
J=[ 0, wd*1, 0;
wd*(1-3*x^2), -k*wd, -f*wd*wd*sin(wd*z);
0, 0, 0];
%Variational equation
F=J*Q;
%Put output data in a column vector
OUT=;
function =lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
n1=n; n2=n1*(n1+1);
%Number of steps
nit = round((tend-tstart)/stept);
% Memory allocation
y=zeros(n2,1); cum=zeros(n1,1); y0=y;
gsc=cum; znorm=cum;
% Initial values
y(1:n)=ystart(:);
for i=1:n1 y((n1+1)*i)=1.0; end;
t=tstart;
% Main loop
for ITERLYAP=1:nit
% Solutuion of extended ODE system
= feval(fcn_integrator,rhs_ext_fcn,,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;
%
% construct new orthonormal basis by gram-schmidt
%
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;
% Output modification
if ITERLYAP==1
Lexp=lp;
Texp=t;
else
Lexp=;
Texp=;
end;
if (mod(ITERLYAP,ioutp)==0)
fprintf('t=%6.4f',t);
for k=1:n1 fprintf(' %10.6f',lp(k)); end;
fprintf('\n');
end;
for i=1:n1
for j=1:n1
y(n1*j+i)=y0(n1*i+j);
end;
end;
end;
=lyapunov(3,@duffing,@ode45,0,0.5,200,,10);
plot(T,Res);
title('Dynamics of Lyapunov exponents');
xlabel('Time'); ylabel('Lyapunov exponents');
得到的结果是：
方法二：
%计算Duffing振子的lyapunov指数
clear all;
clc;
global X V
yinit=;
orthmatrix=[1 0;
0 1];
y=zeros(6,1);
%初始化输入
IC(1:2)=yinit;
IC(3:6)=orthmatrix;
InitiaTime=0;%时间初始值 0
TimeStep=1e-2;%时间步长 0.02
wholetimes=10000;%总的循环次数 10000
UpdateStepNum=10;%每次演化的步数 10
iteratetimes =wholetimes/UpdateStepNum; %演化的次数 1000iterate 迭代的意思
DiscardItr=100; % 100
Iteration=UpdateStepNum*TimeStep; % 0.2
DiscardTime=DiscardItr*Iteration+InitiaTime; % 40
T1=InitiaTime; % T1=0
T2=T1+Iteration; % T2=0.2
Tspan=;% Tspan=
options=odeset('RelTol',1e-5,'AbsTol',1e-6);
mod=zeros(1,2);
d=2;
Sum=zeros(1,d);
n=0;
k=0;
xData=[];
yData=[];
for i=1:iteratetimes %1:1000
n=n+1
=ode45('twofun',Tspan,IC,options);
%取积分得到的最后一个时刻的值
=size(X);
L=cX-d*d;
for i=1:d
m1=L+1+(i-1)*d;
m2=m1+d-1;
A(:,i)=(X(rX,m1:m2))';
end
y0=TwoGS(A);
%取两个向量的模
mod(1)=sqrt(y0(:,1)'*y0(:,1));
mod(2)=sqrt(y0(:,2)'*y0(:,2));
y0(:,1)=y0(:,1)/mod(1);
y0(:,2)=y0(:,2)/mod(2);
if T2>DiscardTime
Q0=y0;
else
Q0=eye(d); %Y = eye(n) returns the n-by-n identity matrix.
end
if (T2>DiscardTime)
k=k+1;
T=k*Iteration;
TT=n*Iteration+InitiaTime;
Sum=Sum+log(abs(mod));
lambda=Sum/T;
Lambda=fliplr(sort(lambda));%B = sort(A) sorts the elements along different dimensions of an array, and
%arranges those elements in ascending order.
%B = fliplr(A) returns A with columns flipped in the left-right direction,
% that is, about a vertical axis.If A is a row vector, then fliplr(A) returns
%a vector of the same length with the order of its elements reversed. If A is
%a column vector, then fliplr(A) simply returns A.
=size(xData);
=size(yData);
xData=;
yData=;
end
ic=X(rX,1:L);
T1=T1+Iteration;
T2=T2+Iteration;
Tspan=;
IC=;
end
plot(xData,yData(:,1),xData,yData(:,2))
function dX=twofun(t,X)
k=0.5;
B=2.804;
wd=2*pi;
x=X(1);
y=X(2);
Y=[X(3) X(5);
X(4) X(6)];
dX=zeros(6,1);
dX(1)=wd*y;
dX(2)=wd*(-k*y+x-x^3+B*cos(wd*t));
J=[ 0,wd*1;
wd*(1-3*x^2), wd*(-k)];
dX(3:6)=J*Y;
function A=TwoGS(V)%V为3*3的向量
global a1 a2
v1=V(:,1);
v2=V(:,2);
a1=zeros(2,1);
a2=zeros(2,1);
a1=v1;
a2=v2-((a1'*v2)/(a1'*a1))*a1;
A=;
得到的结果是：
你犯了一个很典型的错误。
不知道你注意到没有，Lyapunov指数的定义，是基于自治系统，而不是非自治系统的
因此，我们在求LE的时候，都需要将原来的非自治系统转化为自治系统，在根据各种方法求解。 学习了，谢谢
页:
[1]