关于Lyapunov指数计算，欢迎大家到该帖提问

胡晓宇

是不是所有的计算Lyapunov指数的方法都需要用到系统的jacobian矩阵，如果系统的维数很多，方程也很复杂，不能直接写出Jacobian矩阵的，该采用哪一种方法计算Lyapunov指数？

houbeijunguan

octopussheng前辈：
    您好！我在求解超混沌系统的李雅普诺夫指数时遇到困难，恳请您帮忙，不胜感激！
    以下是超混沌系统的微分方程组：
dx/dt = - z - u
dy/dt = 2*y + z
dz/dt = 11*x - 12*y
du/dt = 92*x - 95*u -v + 101*(abs(u + 1) - abs(u - 1))
dv/dt = 15*z - 2*v
        这个状态方程组来自于《基于细胞神经网络的图像加密新算法》，文中指出该系统的李氏指数为：0．2 9 5 3 ， 0．5 2 8 5 ， 0．1 2 6 4 ，一 3．9 2 0 5 ，一 1 7．4 3 8。文中并未说明用何种方法求解的李氏指数。
        这个状态方程组中还存在绝对值（abs），我看了LET工具箱中的lorenzeq.m文件，要求用户然后按照readme.m文件中的方法仿照洛伦兹系统去构造一个ODE function.
        我的问题在于，洛伦兹系统中不存在绝对值，直接得到唯一的杰克比行列式。而我这个带有绝对值的系统方程组该如何构造杰克比行列式呢？而且是超混沌系统，需要得到两个或者两个以上的正的李雅普诺夫指数。我看了好多关于李氏指数的求法的文章和帖子，几乎都是关于混沌而不是超混沌的。看了您的很多帖子，我感觉看到了希望。恳请您能在百忙之中抽出时间指点一下，不胜感激！

houbeijunguan

您好，您看我上面提到的这个系统方程属于您所说的那种情况吗，我也是发愁如何构造杰克比行列式

liuhongjuan

版主，你好，我在用第一个程序（定义法）计算Rossler的Lyapunov指数时，得到的结果是一个Lyapunov指数谱，请问根据它如何计算出你给出的这个Lyapunov指数：0.063231  0.092635  -9.8924
Lyapunov指数谱的结果如下：

zxb1030

let工具箱 可以修改 计算lyapunov指数随某个参数变化的图像，请问怎么修改
我把程序贴上来
% Chapter 13 - Three-Dimensional Autonomous Systems and Chaos.
% Programs_13d - Lyapunov exponents of the Lorenz system.
% Copyright Birkhauser 2009. Stephen Lynch.

% Special thanks to Vasiliy Govorukhin for allowing me to use his M-files.
% For continuous and discrete systems see the Lyapunov Exponents Toolbox of
% Steve Siu at the mathworks/matlabcentral/fileexchange.

% Reference.
% A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, "Determining Lyapunov Exponents from a Time Series," Physica D,
% Vol. 16, pp. 285-317, 1985.
% You must read the above paper to understand how the program works.

% Lyapunov exponents for the Lorenz system below are:
% L_1 = 0.9022, L_2 = 0.0003, L_3 = -14.5691 when tend=10,000.

function [Texp,Lexp]=lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);

n=3;rhs_ext_fcn=@lorenz_ext;fcn_integrator=@ode45;
tstart=0;stept=0.5;tend=300;
ystart=[1 1 1];ioutp=10;
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.
  [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;

% 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;

%  Update running vector magnitudes.

  for k=1:n1 cum(k)=cum(k)+log(znorm(k)); end;

%  Normalize exponent.

  for k=1:n1
      lp(k)=cum(k)/(t-tstart);
  end;

% Output modification.

  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;

% Show the Lyapunov exponent values on the graph.
str1=num2str(Lexp(nit,1));str2=num2str(Lexp(nit,2));str3=num2str(Lexp(nit,3));
plot(Texp,Lexp);
title('Dynamics of Lyapunov Exponents');
text(235,1.5,'\lambda_1=','Fontsize',10);
text(250,1.5,str1);
text(235,-1,'\lambda_2=','Fontsize',10);
text(250,-1,str2);
text(235,-13.8,'\lambda_3=','Fontsize',10);
text(250,-13.8,str3);
xlabel('Time'); ylabel('Lyapunov Exponents');
% End of plot

function f=lorenz_ext(t,X);
%
% Values of parameters.
SIGMA = 10; R = 28; BETA = 8/3;

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)];

f=zeros(9,1);

%Lorenz equation.
f(1)=SIGMA*(y-x);
f(2)=-x*z+R*x-y;
f(3)=x*y-BETA*z;

%Linearized system.
Jac=[-SIGMA, SIGMA,     0;
         R-z,    -1,    -x;
           y,     x, -BETA];
  
%Variational equation.   
f(4:12)=Jac*Y;

%Output data must be a column vector.

% End of Programs_13d.

hillhaw

m和tau是用什么方法得到的？这组数据正确的结果应该是2.08，但这个程序是负的。

hillhaw

每条曲线的最后一个值便是Lyapunov指数了

hillhaw

每条曲线的最后一个值便是Lyapunov指数了。用Matlab，可以把三条曲线分开，看得更仔细

wh.gao

有没有jacobian法的详细介绍呀，matlab程序也行，急需，谢谢！

cqupenghao

回复 gurencai 的帖子


    照道理你这个应该是小数据量法求解最大Lyapunov指数的程序啊，我也用的这个求的，就是我有一个小问题。最后面用的是拟合求直线斜率得到最大Lyapunov指数，可是数据点差异那么大，我个人觉得有些点应该舍去，可是系统还是全部算进去进行拟合，我想问哈这样做得到的结果准确吗？？？

vicen

好的 一起来讨论！

hanweiw

回复 octopussheng 的帖子我得到一组数据，然后用C-C法计算其m和tau，可是为什么我这个程序跑到最后老是报错呢，这样就没办法用小数据量法进行下面的计算了。初学matlab，还请多指教。我的程序如下：
function [tau,m]=CC_Method(data)

x=data;
x=x';

X = [x-mean(x)]/[max(x)-min(x)];    % 归一化到均值为 0，振幅为 1

maxLags = 100;                  % 最大时延
m_vector = 2:5;                 % m 取值范围
sigma = std(X);
r_vector = sigma/2*[1:4];       % r 取值范围

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %初始化
S_mean = zeros(1,maxLags);
Sj = zeros(1,length(r_vector));
delta_S_mean = zeros(1,maxLags);
delta_S = zeros(length(m_vector),maxLags);

for t = 1:maxLags
    temp = 0;
    for i = 1:length(m_vector)
        for j = 1:length(r_vector)
            m = m_vector(i);
            r = r_vector(j);
            S = ccFunction0(m,X,r,t);       % 文献中公式（13）
            temp = temp + S;               % 文献中公式（17）
            Sj(j) = S;
        end
        delta_S(i,t) = max(Sj)-min(Sj);    % delta_S(m,t),文献中公式（15）
    end
    S_mean(t) = temp/(length(m_vector)*length(r_vector));   % 文献中公式（17）
    delta_S_mean = mean(delta_S);                           % 文献中公式（18）
end
S_cor = delta_S_mean + abs(S_mean);                         % 文献中公式（19）

for t=1:maxLags
    if S_mean(t)==0
        St=t;
        break
    end
end

for t=1:maxLags-1
    if delta_S_mean(t)<delta_S_mean(t+1)
        tau=t;
        break
    end
end
   
TscorV=min(S_cor);
for t=1:maxLags
    if TscorV==S_cor(t)
        Tscor=t;
        break
    end
end

m=round(Tscor/tau)+1;
============================
老是报错如下：
??? Undefined function or variable "tau".

Error in ==> CC_Method at 58
m=round(Tscor/tau)+1;
还有没有更好的计算方法呢？


   

mountain_77

在利用陆博士和skyshaw的程序时，总有很多参数无法确定，真是让人烦恼！

ranling

“LET工具箱的算法原理就是Jacobian方法。基本原理就是首先求解出连续系统微分方程的近似解，然后对系统的Jacobian矩阵进行QR分解，计算Jacobian矩阵特征值的乘积，最后计算出LE和分数维。”
请教一下上面那段话是什么意思。我看了很多文献，很多种说法，LE和Jacobian矩阵特征值的乘积有什么关系啊？谢谢！