OCT前辈提供的关于LE指数程序的疑问
本帖最后由 王瑛 于 2013-7-3 22:51 编辑1. 关于连续系统Lyapunov指数的计算方法关于定义法求解的程序，和matlab板块的“连续系统LE求解程序”差不多。以Rossler系统为例
Rossler系统微分方程定义程序
function dX = Rossler_ly(t,X)
%Rossler吸引子，用来计算Lyapunov指数
% a=0.15,b=0.20,c=10.0
% dx/dt = -y-z,
% dy/dt = x+ay,
% dz/dt = b+z(x-c),
a = 0.15;
b = 0.20;
c = 10.0;
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) = -y-z;
dX(2) = x+a*y;
dX(3) = b+z*(x-c);
% Rossler吸引子的Jacobi矩阵
Jaco = [0 -1 -1;
1 a0;
z 0x-c];
dX(4:12) = Jaco*Y;
求解LE代码：
% 计算Rossler吸引子的Lyapunov指数
clear;
yinit = ;
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);
= 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 = ;
计算得到的Rossler系统的LE为————0.0632310.092635-9.8924
Wolf文章中计算得到的Rossler系统的LE为————0.09 0 -9.77
需要注意的是——定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。
正交化程序可以根据上面的扩展到N*N向量，这里就不加以说明了，对matlab用户来说应该还是比较简单的！
本帖最后由 王瑛 于 2013-7-3 22:51 编辑
这个程序计算得到的Rossler系统的LE为————0.0632310.092635-9.8924，第一LE指数应该是最大LE指数对吗？但是程序算出来的大过了第二LE指数。并且混沌系统中，应该要有一个指数为0，也就是第二LE指数为0.我想这个计算结果应该是有问题的。我自己看了一下程序，实在找不出问题在哪。并且这个程序的帖子在论坛这么多年了，许多前辈也都讨论过。我想如果真有问题，应该会早就提出来了。但我对这个结果还是有些质疑。想请明白的前辈指点一下。谢谢 每次遇到问题就跑来问，最后发现都是比较基础的问题。可能是我刚接触这方面，很多知识都不扎实。调一下步长和计算次数，可以得到比较好的结果。
页:
[1]