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

octopussheng

如题，由于Lyapunov指数的那个帖子已经长达9页，不易于大家的讨论，因此新开一个帖子，请大家踊跃发言，讨论Lyapunov指数的计算问题，我会尽量回复。

wxl789

呵呵
我来抢个沙发。。。

http://forum.vibunion.com/forum/thread-73394-1-1.html

谢谢！

octopussheng

1. 数据量确实少了些，我的理解是可以取400个左右的数据点，这样应该可以比较好的反映系统状况了
2. 用小数据量法计算结果据陆振波博士的研究，有以下一些结论：
(1)小数据量法对嵌入维与时延的选取具有比较好的鲁棒性，可以用其它确定时延与嵌入维的方法代替自相关法和 G-P算法，求得的最大 Lyapunov 指数基本一致；
(2)在原混沌序列平均周期无法确定的情况下，不考虑限制短暂分离，也可以得到比较准确的计算结果；
(3)在曲线y(i)~i最终达到饱和之前， 曲线 y(i)-y(i-1)~i 随i的变化相对较小的区域即是理想
的线性区域，线性区域的下界应为(m-1)tau+1；
他的文章里面也谈到了线性区域的选择，你可以看一下。我这里截一部分出来供大家参考

wxl789

谢谢oct 的详细解答。
但有一个地方不明白：“线性区域的下界应为(m-1)tau+1”。这个是什么意思呢？

另外一个：
你截图里边的这些感觉都有一个很好的线性关系，如果做出来的图形没有这么好的线性关系，那么接下来该做如何处理呢？？

请指教。。

octopussheng

这个你最好还是看一下原文，呵呵，那样理解会更深刻。

另外，如果没有很好的线性关系的话，只有按照自己的判断了。

wxl789

谢谢  oct  
有问题继续请教哈。 。 。 。 :lol

tjyan07

最近在做非光滑系统的指数计算，很头痛……
怎样的系统称为耗散系统？耗散系统的指数和一定小于0？有没有指数是趋于负无穷的？另外三维以上的系统是不是至少有一个指数为0?
有一篇文章：At least on Lyapunov exponent vanishes if the trajectory of an attractor does not contain a fixed point(1983)，可能会提到指数为0的问题，但一直找不到……
问题可能有点肤浅，希望各位大侠指教啊……

gurencai

求教：小数据量法求解最大Lyapunov指数
最大Lyapunov指数的问题困惑我好久了，原本以为可用的程序今天验证了一下才发现得不到想要的结果。

我不知道是我没有正确使用还是程序本身有问题，希望各位高手能够帮帮我

我用的程序是：

function lambda_1=largest_lyapunov_exponent(data,N,m,tau,P)
%the function is used to calcultate largest lyapunov exponent with the
%mended algorithm,which put forward by lv jing hu.
%data:the time series
%N:the length of data
%m:enbedding dimention
%tau:time delay
%P:the mean period of the time series,calculated with FFT
%lambda_1:return the largest lyapunov exponent
%skyhawk

data=load('x.txt')
N=5000;
m=8;
tau=10;
P= 60;
delt_t=0.01;


Y=reconstitution(data,N,m,tau );%reconstitute state space
M=N-(m-1)*tau;%M is the number of embedded points in m-dimensional space
for j=1:M
    d_max=1e+100;
    for jj=1:M                                              %寻找相空间中每个点的最近距离点，并记下
        d_s=0;                                              %该点下标
        if abs(j-jj)>P                                      %限制短暂分离
            for i=1:m
                d_s=d_s+(Y(i,j)-Y(i,jj))*(Y(i,j)-Y(i,jj));
                d_min=d_max;
                if d_s<d_min
                   d_min=d_s;
                   idx_j=jj;
               end
            end
        end
    end
%     index(j)=idx_j;
    max_i=min((M-j),(M-idx_j));%计算点j的最大演化时间步长i
    for k=1:max_i              %计算点j与其最近邻点在i个离散步后的距离
        d_j_i=0;
        for kk=1:m
            d_j_i=d_j_i+(Y(kk,j+k)-Y(kk,idx_j+k))*(Y(kk,j+k)-Y(kk,idx_j+k));
            d(k,j)=d_j_i;
        end
    end
end

%对每个演化时间步长i，求所有的j的lnd(i,j)平均
[l_i,l_j]=size(d);
for i=1:l_i
    q=0;
    y_s=0;
    for j=1:l_j
        if d(i,j)~=0
            q=q+1;
            y_s=y_s+log(d(i,j));
        end
    end
    y(i)=y_s/(q*delt_t);
end
x=1:length(y);
pp=polyfit(x,y,1);
lambda_1=pp(1);
yp=polyval(pp,x);
plot(x,y,'-o',x,yp,'--')




function X=reconstitution(data,N,m,tau)
%该函数用来重构相空间
% m为嵌入空间维数
% tau为时间延迟
% data为输入时间序列
% N为时间序列长度
% X为输出,是m*n维矩阵

M=N-(m-1)*tau;%相空间中点的个数
for j=1:M           %相空间重构
    for i=1:m
        X(i,j)=data((i-1)*tau+j);
    end
end

zhuzhukaixin

请教oct：关于求解Rossler系统LE程序求解LE代码：
% 计算Rossler吸引子的Lyapunov指数
clear;
yinit = [1,1,1];
orthmatrix = [1 0 0;
              0 1 0;
              0 0 1];
a = 0.15;
b = 0.20;
c = 10.0;
y = zeros(12,1);
% 初始化输入
y(1:3) = yinit;
y(4:12) = orthmatrix;
tstart = 0; % 时间初始值
tstep = 1e-3; % 时间步长
wholetimes = 1e5; % 总的循环次数
steps = 10; % 每次演化的步数
iteratetimes = wholetimes/steps; % 演化的次数
mod = zeros(3,1);
lp = zeros(3,1);
% 初始化三个Lyapunov指数
Lyapunov1 = zeros(iteratetimes,1);
Lyapunov2 = zeros(iteratetimes,1);
Lyapunov3 = zeros(iteratetimes,1);
for i=1:iteratetimes
    tspan = tstart:tstep:(tstart + tstep*steps);   
    [T,Y] = ode45('Rossler_ly', tspan, y);
    % 取积分得到的最后一个时刻的值
    y = Y(size(Y,1),:);
    % 重新定义起始时刻
    tstart = tstart + tstep*steps;
    y0 = [y(4) y(7) y(10);
          y(5) y(8) y(11);
          y(6) y(9) y(12)];
    %正交化
    y0 = ThreeGS(y0);
    % 取三个向量的模
    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);
    lp = lp+log(abs(mod));
    %三个Lyapunov指数
    Lyapunov1(i) = lp(1)/(tstart);
    Lyapunov2(i) = lp(2)/(tstart);
    Lyapunov3(i) = lp(3)/(tstart);
        y(4:12) = y0';
end
% 作Lyapunov指数谱图
i = 1:iteratetimes;
plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)

程序中用到的ThreeGS程序如下：
%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];
这个程序同LE指数定义（oct总结过）对应吗？ 没看懂呢。。。而且计算结果同let工具箱（Jacobi方法）的计算结果相差甚大。能否就上述问题讲解下，先行谢过~~~

zhuzhukaixin

另求教oct
可否用求解Rossler方程LE的方法求解duffing和Lorenz方程的LE？致谢~

zhuzhukaixin

继续求教
基本原理就是首先求解出连续系统微分方程的近似解，然后对系统的Jacobian矩阵进行QR分解，计算Jacobian矩阵特征值的乘积，最后计算出LE和分数维
是对Jacobian方法求解LE的概括吗？毕竟数学描述过于抽象。计算出乘积后，后续计算LE工作又将如何进行呢？--求教

xiaocheng_2007

楼主，用wolf方法计算Lyapunov指数的方法我看不懂唉:@(
首先是重构相空间，然后取初始点，再找一个最近邻点，两者之间有一个距离。。。。
最近邻点是什么？从哪里得到的？请楼主帮忙详细解释一下:victory:

octopussheng

先从最后的回复吧
12楼——最近邻点我的理解就是相空间中和初始点间距离最短的点，在常见的wolf程序中，一般用这个代码实现：
      for k = 1 : m
                d = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
            end
保证了一次过程所有点和初始点之间的距离都可计算得到，从而开始下一步骤。

11楼——求解Rossler方程LE的方法求解duffing和Lorenz方程的LE，可以的
需要注意一个问题，Rossler是一个自治系统，目前如LET工具箱这些，用的算法都是雅克比算法，即要先对系统方程进行一些处理，如果方程是非自治的，需要增维，将原系统增加一个自由度变成自治系统，在套用Rossler方程LE的求解过程
另一个问题，我想你还是需要再好好读一下这方面的资料，呵呵，这里我有点讲不大清楚了！

9楼——这个程序的主要思想其实在wolf的那篇文章里面有，就是对方程的解进行施密特正交化，再进行求解。你可以看wolf文章里面相关的内容以及后面附录的fortran程序。
另外，不同方法之间肯定有差别，类似这个程序，可通过调节迭代步数等参数，应可以得到较好的结果

8楼——请将你遇到的问题贴出来吧，我不可能帮你算的，呵呵

7楼——不好意思，非光滑系统目前没有做过，不过你说的是有可能的。

xiaocheng_2007

楼主，请问求李雅谱诺夫指数时用的P，即平均周期是从哪里来的？

yuling

可以使用FFT计算。这是我编的的计算平均周期的程序，你可以参考一下。

function T_mean=period_mean_fft(data)
%该函数使用快速傅里叶变换FFT计算序列平均周期
%data：时间序列
%T_mean：返回快速傅里叶变换FFT计算出的序列平均周期

Y = fft(data);       %快速FFT变换
N = length(data);    %FFT变换后数据长度
Y(1) = [];           %去掉Y的第一个数据，它是data所有数据的和
power = abs(Y(1:N/2)).^2;  %求功率谱
nyquist = 1/2;

freq = (1:N/2)/(N/2)*nyquist; %求频率
subplot(121)
plot(freq,power); grid on     %绘制功率谱图
xlabel('频率')
ylabel('功率')
title('功率谱图')

period = 1./freq;                %计算周期
subplot(122)
plot(period,power); grid on  %绘制周期－功率谱曲线
ylabel('功率')
xlabel('周期')
title('周期—功率谱图')

[mp,index] = max(power);       %求最高谱线所对应的下标
T_mean=period(index);            %由下标求出平均周期