请大家帮忙分析定义法计算chen’s系统的Lyapunov指数有什么错误？

Jessie5663

我在用定义法计算chen’s系统的Lyapunov指数过程中发现一个非常奇怪的现象。
我参考论坛中的《【总结】Lyapunov指数的计算方法》一篇，把Rossler系统的方程改为chen’s系统的方程。

由语句
[T,Y] = ode45('chen_ly', tspan, y);
可知，程序中Y存储了每一次演化的结果。则Y(:,1:3)为：
5.5302    7.7102    21.929
5.5378    7.7158    21.927
5.5454    7.7214    21.925
5.553      7.727      21.922
5.5606    7.7326    21.92
……
6.2377    8.2289    21.75
6.2447    8.234      21.748
6.2517    8.2391    21.747
6.2586    8.2442    21.746
但是参考吕金虎《混沌时间序列分析及其应用》一书
-11.6126    -12.5833    22.9856
-11.6160    -12.5837    22.9933
-11.6194    -12.5841    23.0010
-11.6228    -12.5845    23.0088
-11.6261    -12.5848    23.0165
……
-11.8868    -12.5364    23.7326
-11.8891    -12.5349    23.7403
-11.8913    -12.5335    23.7481
-11.8936    -12.5320    23.7559
-11.8958    -12.5305    23.7637

都是
yinit = [0,1,0];
tstart = 0;                      % 时间初始值
tstep = 0.0001;             % 时间步长
wholetimes =500000;   % 总的循环次数
为什么会有差别呢？
谢谢！

[ 本帖最后由 Jessie5663 于 2008-3-12 10:20 编辑 ]

octopussheng

你这个结果确实有问题，用的是哪些代码？贴一下看看！

Jessie5663

修改的：

function dX = chen_ly(t,X)
a = 35;
b = 3;
c = 28;
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);
% Rossler吸引子
dX(1) = a*(y-x);
dX(2) = (c-a)*x-x*z+c*y;
dX(3) = x*y-b*z;
% Rossler吸引子的Jacobi矩阵
Jaco = [-a a 0;
        c-a-z c  -x;
        y x -b];
dX(4:12) = Jaco*Y;



% 计算Rossler吸引子的Lyapunov指数
clear;
yinit = [0,1,0];
orthmatrix = [1 0 0;
              0 1 0;
              0 0 1];
a = 35;
b = 3;
c = 28;
y = zeros(12,1);
% 初始化输入
y(1:3) = yinit;
y(4:12) = orthmatrix;
tstart = 0; % 时间初始值
tstep = 0.0001; % 时间步长
wholetimes =500000; % 总的循环次数
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('chen_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程序未做改动。

Jessie5663

% chen's attractor
% 模型
% chen方程            dx=a*(y-x);
%                     dy=(c-a)*x-x*z+c*y;                                       
%                     dz=x*y-b*z;
clear all;
clc;
global a;
global b;
global c;
% 参数设定
a=35;
b=3;
c=28;
% 数据初始化
x0=[0;1;0];                     % 初始值
td=0.0001;                      % 迭代步长
t_time_abandon=0;           % 放弃的迭代次数
t_time=500000;                  % 选择的迭代次数
tspan=[0:td:td*(t_time+t_time_abandon-1)];
[t,y]=ode45('chen_de',tspan,x0);    % 求解方程
% 放弃前t_time_abandon个不稳定的迭代结果
y=y';
y(1:1:t_time_abandon*3)=[];
y=reshape(y,3,t_time);
t=t';
t(1:1:t_time_abandon)=[];



微分方程
function dy=chen_de(t,y)
global a;
global b;
global c;
dy=[a*(y(2)-y(1));(c-a)*y(1)-y(1)*y(3)+c*y(2);y(1)*y(2)-b*y(3)];

[ 本帖最后由 Jessie5663 于 2008-3-13 09:07 编辑 ]

octopussheng

问题1：你在楼主贴中问的问题应该是，Chen's系统随时间演化时计算得到的数据为什么和吕金虎一书中的数据差别那么大？对不？

问题2：初始条件设置是否和书中一致？

问题3：t_time_abandon=0;           % 放弃的迭代次数，你设置为0，就是不舍弃求解的结果，这个有些问题!

Jessie5663

感谢octopussheng！

我在修改的程序中设定：
yinit = [0,1,0];
tstart = 0; % 时间初始值
tstep = 0.0001; % 时间步长
wholetimes =500000; % 总的循环次数
steps = 10; % 每次演化的步数
iteratetimes = wholetimes/steps; % 演化的次数


在参考吕金虎编程序中设置：
x0=[0;1;0];                     % 初始值
td=0.0001;                      % 迭代步长
t_time_abandon=0;           % 放弃的迭代次数
t_time=500000;                  % 选择的迭代次数
tspan=[0:td:td*(t_time+t_time_abandon-1)];

因为吕金虎除去了前10000个点，所以我在编程时，变量定义t_time_abandon来表示出去的点数。但是在修改的程序中，没有考虑去除点的问题，故设t_time_abandon=0。

我始终想不明白初始条件一样、步长一样、总步数一样，为什么每次演化10步，演化50000次和一口气演化500000步不一样呢？

octopussheng

我看是否因为在循环计算过程中，演化次数的不同造成的呀？

xw263

我按照 octopussheng 提供的程序，也做了一下Chen系统，但是无论如何也得不到吕金虎的结果。
真是郁闷呀。既使去掉了前面10000个点，也不行。在吕金虎的论文中，明显有一个LE指数是为0的。倒底是怎么回事呢？？？？


还有：请教  octopussheng  在程序中

Jaco = [-a a 0;
        c-a-z c  -x;
        y x -b];
dX(4:12) = Jaco*Y;


是什么意思？线性化方程组吗？？？？？

octopussheng

吕金虎论文中有些结果我觉得不大可信，做个参考就可以了。这些经典系统最好还是参考其他的一些经典文献。

这两句话的目的是得到原系统的线性化系统。

xw263

y(4:12) = y0';

本人觉得这条语句有错误,应该改为:  y(4:12) = y0;

octopussheng

这句话应是没有问题的。