连续系统Lyapunov指数计算

hohoo

根据LET工具箱的代码，降le指数正交化的方法进行了改进，以duffing系统为例，和LET的结果对比了一下，看看那个结果更精确
主程序如下：
% 计算duffing吸引子的Lyapunov指数
clear all;
clc;
% 初始化输入
yinit = [1,1,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 = 1e5;                       % 总的循环时间
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 i=1:iteratetimes
     n=n+1;  
    [t,X] =ode45('duffingfun', 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);
   
   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
plot(xData,yData(:,1),xData,yData(:,2),xData,yData(:,3))

function dX = duffingfun(t,X)
k=0.1;
B=11;
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)=-k*y-x^3+B*cos(z);
dX(3)=1;  
J=[    0,    1,          0;
   -3*x^2,  -k,  -B*sin(z);
        0,   0,          0];
dX(4:12) = J*Y;

function A = ThreeGS(V)
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];

无水1324

你的结论是什么样的？

hohoo

计算结果是le1=0.11398975  le2=0   le3=-0.21398975
LET的计算结果是0.10825     0      -0.20825

参考文献中的结果是0.1140   0   -0.2140

hohoo

尝试了一下不改为自治系统的Duffing系统的le计算，结果和上面的有一点误差，应该是由options的设置有关，如果提高精度，应该能够得到和非自治系统一样的结论
程序如下：
% 计算duffing吸引子的Lyapunov指数
clear all;
clc;
yinit = [1,1];
orthmatrix = [1 0 ;
              0 1 ];            
y = zeros(6,1);
% 初始化输入
IC(1:2) = yinit;
IC(3:6) = orthmatrix;
InitialTime = 0; % 时间初始值
TimeStep = 1e-2; % 时间步长
wholetimes = 1e5; % 总的循环次数
UpdateStepNum = 10; % 每次演化的步数
iteratetimes = wholetimes/UpdateStepNum; % 演化的次数
DiscardItr=200;
Iteration=UpdateStepNum*TimeStep;       
DiscardTime=DiscardItr*Iteration+InitialTime;

T1=InitialTime;
T2=T1+Iteration;
TSpan=[T1:TimeStep:T2];
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
     n=n+1;  
    [t,X] = ode45('twofun', 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 = 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);
      end
   
  if (T2>DiscardTime )
      k=k+1;
      T=k*Iteration;
      TT=n*Iteration+InitialTime;
      Sum=Sum+log(abs(mod));
      lambda=Sum/T;

      Lambda=fliplr(sort(lambda));      
    [rxD,cxD]=size(xData);
      [ryD,cyD]=size(yData);
     
         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
plot(xData,yData(:,1),xData,yData(:,2))


function dX = twofun(t,X)
k=0.1;
B=11;
x=X(1); y=X(2);
Y=[X(3),X(5);
   X(4),X(6)];   
dX = zeros(6,1);
dX(1)=y;
dX(2)=-k*y-x^3+B*cos(t);       
J=[    0,    1;         
   -3*x^2,  -k];
dX(3:6) = J*Y;

function A = TwoGS(V)  % V 为3*3向量
v1 = V(:,1);
v2 = V(:,2);
a1 = zeros(2,1);
a2 = zeros(2,1);
a1 = v1;
a2 = v2-((a1'*v2)/(a1'*a1))*a1;
A = [a1,a2];

结果为0.10953  -0.20953

hohoo

这里我想问一下大家 计算le指数时，为什么通常要将非自治系统改为自治系统？ 有什么理论根据吗？

octopussheng

目前的理论上推导的LE大多是基于自治系统的，但我们编程求解的时候，非自治系统还是作为非自治系统来求解的，如果直接改写肯定是不行的！
我不知道你怎么会有这个疑问的，毕竟相轨迹这些信息都是基于原系统求解得到的啊！

hohoo

let工具箱是统统将非自治系统转换为自治系统求解的，就想问一下这样做的根据。

chuandong418

想问下如果将系统的方程线性化后，如何调用let来运行程序！

hohoo

将你的系统方程做成函数文件，在let的界面上输入你的初始参数即可

octopussheng

不知道你发现没有，其实你函数的前三个方程仍然是非自治的呢？

hohoo

我第一次发的程序是将duffing方程改为自治系统的计算程序，是三维的，有三个le指数  第二次发的非自治的计算程序，是二维的，两个le指数

octopussheng

其实程序里面将原方程进行线性化，其目的是为了得到Jacobian矩阵T，你看看下面的原理吧，应该没有问题了！

其实这些在我的那个总结帖子里面都有讲的，可能内容也是太多了点，呵呵！

chuandong418

哦，谢谢

shenyongjun

原帖由 hohoo 于 2007-8-15 09:32 发表

                               
登录/注册后可看大图

根据LET工具箱的代码，降le指数正交化的方法进行了改进，以duffing系统为例，和LET的结果对比了一下，看看那个结果更精确
主程序如下：
% 计算duffing吸引子的Lyapunov指数
clear all;
clc;
% 初始化输入
...

你说的“改进”，或者说与LET不同之处在程序里哪一部分？

shenyongjun

原帖由 hohoo 于 2007-8-15 10:07 发表

                               
登录/注册后可看大图

计算结果是le1=0.11398975  le2=0   le3=-0.21398975
LET的计算结果是0.10825     0      -0.20825

参考文献中的结果是0.1140   0   -0.2140

你的程序和LET中的UpdateStepNum不一样，改为9后，得到的结果为：
0.10974   -0.20974   0
不知道为什么LET中总要求UpdateStepNum为整数的平方？呵呵，也许是我看错了