数学与非线性科学 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

hohoo 发表于 2007-8-15 09:32

连续系统Lyapunov指数计算

根据LET工具箱的代码，降le指数正交化的方法进行了改进，以duffing系统为例，和LET的结果对比了一下，看看那个结果更精确
主程序如下：
% 计算duffing吸引子的Lyapunov指数
clear all;
clc;
% 初始化输入
yinit = ;
orthmatrix = [1 0 0;
         0 1 0;
         0 0 1];
y = zeros(12,1);
IC(1:3) = yinit;
IC(4:12) = orthmatrix;
InitialTime = 0;                      % 时间初始值
TimeStep = 1e-2;                      % 时间步长
wholetimes = 1e5;                   % 总的循环时间
UpdateStepNum = 10;                   % 每次演化的步数
iteratetimes = wholetimes/UpdateStepNum; % 演化的次数
DiscardItr=200;                      %抛弃的不稳定迭代次数
Iteration=UpdateStepNum*TimeStep;
DiscardTime=DiscardItr*Iteration+InitialTime;
mod = zeros(1,3);
d=3;                      %维数
Sum=zeros(1,d);
n=0; %总的迭代计数
k=0; %对结果有作用的迭代计数
xData=[];
yData=[];
T1=InitialTime;
T2=T1+Iteration;
TSpan=;
options=odeset('RelTol',1e-5,'AbsTol',1e-6);
for i=1:iteratetimes
n=n+1;
=ode45('duffingfun', TSpan, IC,options);
=size(X);
L=cX-d*d;

% 取积分得到的最后一个时刻的值
for i=1:d
   m1=L+1+(i-1)*d;
   m2=m1+d-1;
   A(:,i)=(X(rX,m1:m2))';
end
y0 = ThreeGS(A);                   %正交化
% 取三个向量的模
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);

if T2>DiscardTime
      Q0=y0;
else
      Q0=eye(d);
end

if (T2>DiscardTime )
   k=k+1;
   T=k*Iteration;
   TT=n*Iteration+InitialTime;
   %计算lyapunov指数
   Sum=Sum+log(abs(mod));
   lambda=Sum/T;
   %按降序排列指数
   Lambda=fliplr(sort(lambda));
   xData=;
   yData=;
end
% 重新定义起始时刻
ic=X(rX,1:L);
   T1=T1+Iteration;
   T2=T2+Iteration;
   TSpan=;
   IC=;
end
plot(xData,yData(:,1),xData,yData(:,2),xData,yData(:,3))

function dX = duffingfun(t,X)
k=0.1;
B=11;
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)];
dX = zeros(12,1);
dX(1)=y;
dX(2)=-k*y-x^3+B*cos(z);
dX(3)=1;
J=[ 0, 1,       0;
-3*x^2,-k,-B*sin(z);
   0, 0,       0];
dX(4:12) = J*Y;

function A = ThreeGS(V)
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 = ;

无水1324 发表于 2007-8-15 09:35

回复 #1 hohoo 的帖子

你的结论是什么样的？

hohoo 发表于 2007-8-15 10:07

计算结果是le1=0.11398975le2=0 le3=-0.21398975
LET的计算结果是0.10825 0 -0.20825

参考文献中的结果是0.1140 0 -0.2140

hohoo 发表于 2007-8-15 10:22

尝试了一下不改为自治系统的Duffing系统的le计算，结果和上面的有一点误差，应该是由options的设置有关，如果提高精度，应该能够得到和非自治系统一样的结论
程序如下：
% 计算duffing吸引子的Lyapunov指数
clear all;
clc;
yinit = ;
orthmatrix = [1 0 ;
         0 1 ];
y = zeros(6,1);
% 初始化输入
IC(1:2) = yinit;
IC(3:6) = orthmatrix;
InitialTime = 0; % 时间初始值
TimeStep = 1e-2; % 时间步长
wholetimes = 1e5; % 总的循环次数
UpdateStepNum = 10; % 每次演化的步数
iteratetimes = wholetimes/UpdateStepNum; % 演化的次数
DiscardItr=200;
Iteration=UpdateStepNum*TimeStep;
DiscardTime=DiscardItr*Iteration+InitialTime;

T1=InitialTime;
T2=T1+Iteration;
TSpan=;
options=odeset('RelTol',1e-5,'AbsTol',1e-6);
mod = zeros(1,2);
d=2;
Sum=zeros(1,d);

n=0;
k=0;
xData=[];
yData=[];

for i=1:iteratetimes
n=n+1;
= ode45('twofun', TSpan, IC,options);
% 取积分得到的最后一个时刻的值
=size(X);
L=cX-d*d;

for i=1:d
   m1=L+1+(i-1)*d;
   m2=m1+d-1;
   A(:,i)=(X(rX,m1:m2))';
end
y0 = TwoGS(A);

% 取三个向量的模
mod(1) = sqrt(y0(:,1)'*y0(:,1));
mod(2) = sqrt(y0(:,2)'*y0(:,2));
y0(:,1) = y0(:,1)/mod(1);
y0(:,2) = y0(:,2)/mod(2);
if T2>DiscardTime
      Q0=y0;
   else
      Q0=eye(d);
   end

if (T2>DiscardTime )
   k=k+1;
   T=k*Iteration;
   TT=n*Iteration+InitialTime;
   Sum=Sum+log(abs(mod));
   lambda=Sum/T;

   Lambda=fliplr(sort(lambda));
=size(xData);
   =size(yData);

      xData=;
      yData=;
end
ic=X(rX,1:L);
   T1=T1+Iteration;
   T2=T2+Iteration;
   TSpan=;

IC=;
end
plot(xData,yData(:,1),xData,yData(:,2))

function dX = twofun(t,X)
k=0.1;
B=11;
x=X(1); y=X(2);
Y=[X(3),X(5);
X(4),X(6)];
dX = zeros(6,1);
dX(1)=y;
dX(2)=-k*y-x^3+B*cos(t);
J=[ 0, 1;
-3*x^2,-k];
dX(3:6) = J*Y;

function A = TwoGS(V)% V 为3*3向量
v1 = V(:,1);
v2 = V(:,2);
a1 = zeros(2,1);
a2 = zeros(2,1);
a1 = v1;
a2 = v2-((a1'*v2)/(a1'*a1))*a1;
A = ;

结果为0.10953-0.20953

hohoo 发表于 2007-8-15 10:26

这里我想问一下大家计算le指数时，为什么通常要将非自治系统改为自治系统？有什么理论根据吗？

octopussheng 发表于 2007-8-15 10:39

回复 #5 hohoo 的帖子

目前的理论上推导的LE大多是基于自治系统的，但我们编程求解的时候，非自治系统还是作为非自治系统来求解的，如果直接改写肯定是不行的！
我不知道你怎么会有这个疑问的，毕竟相轨迹这些信息都是基于原系统求解得到的啊！

hohoo 发表于 2007-8-15 10:47

let工具箱是统统将非自治系统转换为自治系统求解的，就想问一下这样做的根据。

chuandong418 发表于 2007-8-15 10:54

想问下如果将系统的方程线性化后，如何调用let来运行程序！

hohoo 发表于 2007-8-15 11:02

将你的系统方程做成函数文件，在let的界面上输入你的初始参数即可

octopussheng 发表于 2007-8-15 11:05

回复 #7 hohoo 的帖子

不知道你发现没有，其实你函数的前三个方程仍然是非自治的呢？

hohoo 发表于 2007-8-15 11:12

我第一次发的程序是将duffing方程改为自治系统的计算程序，是三维的，有三个le指数第二次发的非自治的计算程序，是二维的，两个le指数

octopussheng 发表于 2007-8-15 11:14

其实程序里面将原方程进行线性化，其目的是为了得到Jacobian矩阵T，你看看下面的原理吧，应该没有问题了！

其实这些在我的那个总结帖子里面都有讲的，可能内容也是太多了点，呵呵！

chuandong418 发表于 2007-8-15 12:25

回复 #9 hohoo 的帖子

哦，谢谢

shenyongjun 发表于 2007-8-15 17:42

原帖由 hohoo 于 2007-8-15 09:32 发表 http://www.chinavib.com/forum/images/common/back.gif
根据LET工具箱的代码，降le指数正交化的方法进行了改进，以duffing系统为例，和LET的结果对比了一下，看看那个结果更精确
主程序如下：
% 计算duffing吸引子的Lyapunov指数
clear all;
clc;
% 初始化输入
...
你说的“改进”，或者说与LET不同之处在程序里哪一部分？

shenyongjun 发表于 2007-8-15 17:45

原帖由 hohoo 于 2007-8-15 10:07 发表 http://www.chinavib.com/forum/images/common/back.gif
计算结果是le1=0.11398975le2=0 le3=-0.21398975
LET的计算结果是0.10825 0 -0.20825

参考文献中的结果是0.1140 0 -0.2140
你的程序和LET中的UpdateStepNum不一样，改为9后，得到的结果为：
0.10974 -0.20974 0
不知道为什么LET中总要求UpdateStepNum为整数的平方？呵呵，也许是我看错了

页: [1] 2 3

声振论坛's Archiver

连续系统Lyapunov指数计算

回复 #1 hohoo 的帖子

回复 #5 hohoo 的帖子

回复 #7 hohoo 的帖子

回复 #9 hohoo 的帖子