Lyapunov指数随参数变化的程序

bohua1208

function dX = DF(t,X,flag,alafa)
global alafa;
a0=0.05;f1=5;f2=10;mu=0.1;w0=1;beta=0.1;ks=2;
b0=(f1-(f1^2-4*mu^2*w0^2*a0^2)^0.5)/(2*mu*w0);
omega=1;
x=X(1);y=X(2);z=X(3);
Y = [X(4), X(7), X(10);
    X(5), X(8), X(11);
    X(6), X(9), X(12)];
dX = zeros(12,1);
dX(1)=y;
dX(2)=-(w0^2*x-ks*beta*x^2+ks*alafa*x^3+ks*mu*y-ks*f1*cos(z)-ks*f2*cos(2*z));
dX(3)=omega;
J=[0                                  1        0 ;        
   -w0^2+2*ks*beta*x-3*ks*alafa*x^2   -mu*ks   ks*f1*sin(z)+2*ks*f2*sin(2*z);
   0                                  0        0];
dX(4:12) = J*Y;
% 计算duffing吸引子的Lyapunov指数谱，随参数alafa变化
clear all;
clc;
Z1=[ ];
Z2=[ ];
Z3=[ ];
global alafa;
% 初始化输入
a0=0.05;f1=5;f2=10;mu=0.1;w0=1;beta=0.1;ks=2;
b0=(f1-(f1^2-4*mu^2*w0^2*a0^2)^0.5)/(2*mu*w0);
yinit = [a0,b0,1];
orthmatrix = [1 0 0;
              0 1 0;
              0 0 1];
y = zeros(12,1);
IC(1:3) = yinit;
IC(4:12) = orthmatrix;
InitialTime = 0;                        % 时间初始值
TimeStep = 1e-2;                        % 时间步长
wholetimes = 1e2;                       % 总的循环时间
UpdateStepNum = 10;                     % 每次演化的步数
iteratetimes = wholetimes/UpdateStepNum; % 演化的次数
DiscardItr=200;                         %抛弃的不稳定迭代次数
Iteration=UpdateStepNum*TimeStep;
DiscardTime=DiscardItr*Iteration+InitialTime;
mod = zeros(1,3);
d=3;                          %维数
Sum=zeros(1,d);
n=0;   %总的迭代计数
k=0;   %对结果有作用的迭代计数
xData=[];
yData=[];
T1=InitialTime;
T2=T1+Iteration;
TSpan=[T1:TimeStep:T2];
options=odeset('RelTol',1e-5,'AbsTol',1e-6);
for alafa=linspace(0.8,2.4,200);
        lp=0;
for i=1:iteratetimes
     n=n+1
    [t,X] =ode45('DF1', TSpan, IC,options);
    [rX,cX]=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 = ThreeGS(A);                      %正交化
    % 取三个向量的模
    mod(1) = sqrt(y0(:,1)'*y0(:,1));
    mod(2) = sqrt(y0(:,2)'*y0(:,2));
    mod(3) = sqrt(y0(:,3)'*y0(:,3));
      
    y0(:,1) = y0(:,1)/mod(1);
    y0(:,2) = y0(:,2)/mod(2);
    y0(:,3) = y0(:,3)/mod(3);
     y(4:12) = y0';
   
   if T2>DiscardTime
         Q0=y0;
   else
         Q0=eye(d);
   end
      
  if (T2>DiscardTime )
      k=k+1;
      T=k*Iteration;
      TT=n*Iteration+InitialTime;
      %计算lyapunov指数
      Sum=Sum+log(abs(mod));
      lambda=Sum/T;
      %按降序排列指数
      Lambda=fliplr(sort(lambda));
        xData=[xData;TT];
        yData=[yData;lambda];            
     end
    % 重新定义起始时刻
     ic=X(rX,1:L);
      T1=T1+Iteration;
      T2=T2+Iteration;
      TSpan=[T1:TimeStep:T2];
      IC=[ic(:);y0(:)];  
   
   
end
end

alafa=linspace(0.8,2.4,200);
plot(alafa,yData(:,1),'-',alafa,yData(:,2),'-',alafa,yData(:,3),'-');
title('Lyapunov exponents spectrum of duffing')
xlabel('alafa'),ylabel('lyapunov exponents')
grid on
%G-S正交化
function A = ThreeGS(V)  % V 为3*3向量
v1 = V(:,1);
v2 = V(:,2);
v3 = V(:,3);
a1 = zeros(3,1);
a2 = zeros(3,1);
a3 = zeros(3,1);
a1 = v1;
a2 = v2-((a1'*v2)/(a1'*a1))*a1;
a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;
A = [a1,a2,a3];
matlab运行时出现??? Error using ==> plot
Vectors must be the same lengths.
Error in ==> canshu at 93
plot(alafa,yData(:,1),'-',alafa,yData(:,2),'-',alafa,yData(:,3),'-');
怎么回事，请高手指教

hsfy919

回复 1 # bohua1208 的帖子

绘图的矩阵必须是相同的长度

lhy

请问一下这个程序绘制出来的最大LE 指数谱对吗？

liliangbiao

凡是用这个程序计算指数图 我都没有得到正确的结果，我不知道错在哪里？

lhy

回复 4 # liliangbiao 的帖子

我更改了一下程序，做出来结果了，但是结果是这样的http://forum.vibunion.com/thread-106334-1-1.html 指数全部都是正数，怎么回事呢？

chenshumeionly

版主，你的问题怎么解决的？

chenshumeionly

回复 4 # liliangbiao 的帖子

你好，看来对Lyapunov指数的作图很有经验了，那我想问问duffing系统作Lyapunov指数图用哪种程序好呢？