bohua1208 发表于 2011-4-20 18:14

Lyapunov指数随参数变化的程序

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 = ;
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=;
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
    =ode45('DF1', 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 = 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=;
      yData=;            
   end
    % 重新定义起始时刻
   ic=X(rX,1:L);
      T1=T1+Iteration;
      T2=T2+Iteration;
      TSpan=;
      IC=;
   
   
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 = ;
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 发表于 2011-4-21 08:35

回复 1 # bohua1208 的帖子

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

lhy 发表于 2011-10-12 17:14

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

liliangbiao 发表于 2011-10-13 22:03

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

lhy 发表于 2011-10-18 09:10

回复 4 # liliangbiao 的帖子

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

chenshumeionly 发表于 2012-7-27 23:24

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

chenshumeionly 发表于 2012-8-3 20:06

回复 4 # liliangbiao 的帖子

你好，看来对Lyapunov指数的作图很有经验了，那我想问问duffing系统作Lyapunov指数图用哪种程序好呢？
页: [1]
查看完整版本: Lyapunov指数随参数变化的程序