Lyapunov指数的计算

chenwenjuan86

T=[1.49 2 1;-2 1.7 0;4 -4 2];
x=[x1;x2; x3];
f=-x+T*tanh(x);

这个系统的lyapunov指数怎么计算呢？
高手请指教一下，谢谢！

chenwenjuan86

T=[1.49 2 1;-2 1.7 0;4 -4 2];
x=[x1;x2; x3];
x'=-x+T*tanh(x);

tt1341

用程序计算的

[ 本帖最后由 xinyuxf 于 2006-12-23 09:50 编辑 ]

octopussheng

楼上的大哥，有源码吗？？？

chenwenjuan86

能把源代码贴上来吗？
谢谢！

chenwenjuan86

哪位高手把算法贴上来吧，急用！！！
谢谢！！！

chinacer

最讨厌那些张口就要源码的人了。
这些东西只有自己去亲自编写程序了，才能有更深的了解。
想当初我为了写一个程序，花费了两个月时间

chenwenjuan86

论坛是交流和学习的地方，不是水吧。
我是新手，看了论坛里的一些程序，不是很明白。
我只是希望高手能够帮忙解答一下。

chenwenjuan86

function dX = chaos(t,X)
x=X(1); y=X(2); z=X(3);

% Y的三个列向量为相互正交的单位向量
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) = -X(1)+149/100*tanh(X(1))+2*tanh(X(2))+tanh(X(3))
dX(2) = -X(2)-2*tanh(X(1))+17/10*tanh(X(2))
dX(3) = -X(3)+4*tanh(X(1))-4*tanh(X(2))+2*tanh(X(3))

% 吸引子的Jacobi矩阵
Jaco = [ 49/100-149/100*tanh(X(1))^2  2-2*tanh(X(2))^2               1-tanh(X(3))^2;
        -2+2*tanh(X(1))^2             7/10-17/10*tanh(X(2))^2        0;
        4-4*tanh(X(1))^2              -4+4*tanh(X(2))^2              1-2*tanh(X(3))^2];
        
dX(4:12) = Jaco*Y;




yinit = [0,0,0.1];
orthmatrix = [1 0 0;
              0 1 0;
              0 0 1];

y = zeros(12,1);
% 初始化输入
y(1:3) = yinit;
y(4:12) = orthmatrix;

tstart = 0; % 时间初始值
tstep = 1e-3; % 时间步长
wholetimes = 1e4; % 总的循环次数
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('chaos', 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 = orth(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)

chenwenjuan86

算出来的结果是
lp =1.0e-013 *
   -0.3331
   -0.5040
   -0.4030
但正确结果应该是三楼给出的，哪里出错了呢？

dc1hawk

楼上的，你的问题解决了没？我最近也在做这个东东，你好像%正交化    y0 = orth(y0); 这个有点早了

sizhiyuan2006

谢谢楼主分享