关于最大Lyapunov指数计算的几个疑问(非重复提问)
本帖最后由 chunshui2003 于 2010-11-28 16:45 编辑首先要感谢无水、oct、咕噜噜、gghhjj和liliangbiao等在求解非线性微分方程、制作poincare映射图和分岔图过程中给予的帮助，虽然有的问题没有在论坛中提问，但前辈的帖子和解释让我受益匪浅；其次，感谢 vickyvictoria和hsfy919在进一步进行poincare分析和分岔图制作中提供的帮助和解惑（http://forum.vibunion.com/thread-96554-1-1.html）。当我完成了分岔图的制作后，发现在某些区间内结合poincare图，系统出现了混沌现象，按照文献Chang-Jian, C., Non-linear dynamic analysis of dual flexible rotors supported by long journal bearings. 2010. 45(6): p. 844-866上的描述，如果加入最大Lyapunov指数图，并配合poincare映射图，可以更好地说明系统的运动现象。所以，最近几天一直在论坛中搜索相关内容，试图找到一个适合计算最大LE的方法。看了gghhjj、oct和yufeiqun2008等前辈的帖子，现在有如下几个问题不是很清楚，希望得到解答：
1
计算的原始数据
这个问题我也看到有人问过，但没有回答，可能很基础，但我确实不清楚。我的问题是，在采用CC法和小数据量法中都要涉及到数据，根据吕老师的书上写长度为3000。如果通过ode求解方程，如 = ode45(@fangcheng,(0:period/100:10000*period),y0)，那么这些数据在过滤掉瞬态后，是像poincare那样每一个周期选取一个点，达到3000个周期，还是连续选取3000个积分数值点。2 关于CC法a.
按照yufeiqun2008的帖子http://forum.vibunion.com/forum-viewthread-tid-57387-highlight-Lyapunov.html，其代码如下：
结合吕金虎老师的书，得知延迟时间tau为dS_t vs tmax图中第一个极小值点所对应的t值；而S_cor vs tmax 图中最小值点对应的t_min为时间窗口，嵌入维数按照网友的说法为m = t_min*delta_t / tau。（其中delta_t为采样时间，这个我也不是很清楚，在后面有说明）。不知道这样的计算tau和m方法是否正确。
b.
是否只有通过dS_t vs tmax图 和 S_cor vs tmax图才能确认tau和m也看了其他的CC法帖子，但上面都没有给出直接的tau和m计算表达式，而是通过图形观察得到。那么，是不是意味着如果采用CC法，延迟时间和嵌入维数只能通过观察图形的方式获得？c.
数据不同是否tau和m就不同微分方程相同，如果更改了参数而得到不同的数据，那么是不是tau和m就不同了呢。如果相同，那么在采用小数据量法时，tau、m可以采用同一值；如果不同的话，是不是每次都要计算tau和m，或是有其他的方法。（在我的电脑上计算一次大概需要耗时30分钟） 本帖最后由 chunshui2003 于 2010-11-28 16:41 编辑
不知道为什么程序会变成那个样子，重新贴一下,按照yufeiqun2008的帖子http://forum.vibunion.com/forum-viewthread-tid-57387-highlight-Lyapunov.html：
function =CC(x)
%% C-C方法提取延迟时间和延迟时间窗口
% by Yu Feiqun
% 2007-12-09
N=length(x);
% 计算给定序列的标准差 deta
deta=std(x);
tmax=200;
t=1:tmax;
m=2:5;
r=1:4;
St=zeros(length(t),length(m),length(r));
%dSt=zeros(length(t),length(m),length(r));
for tt=1:length(t)
% 计算时间序列分为t个不相交的时间序列时，每一时间序列的长度
L=floor(N/t(tt));
for mm=1:length(m)
for rr=1:length(r)
rjj=deta*r(rr)/2;
smrt=0;
for ss=1:t(tt)
% 提取第ss个子序列
xsub=x();
xsub=xsub();
c1=getCValue(xsub,m(mm),rjj,t(tt));
c2=getCValue(xsub,1,rjj,t(tt));
smrt=smrt+c1-c2.^m(mm);
end
St(tt,mm,rr)=smrt/t(tt);
end
end
end
S_t=sum(sum(St,3),2)/16;
dS_m_t=max(St,[],3)-min(St,[],3);
dS_t=sum(dS_m_t,2)/4;
S_cor=abs(S_t)+dS_t;
function y=getCValue(x,m,r,t)
% 计算时间序列x的关联积分
% m为嵌入维数 t为时间延迟 N为序列长度 r为半径
% 由于序列x已经经过划分,其相邻两个数据已具有时间延迟t
N=length(x);
M=N-(m-1);
% 重构相空间 嵌入维数为m 时间延迟为t
xv=zeros(m,M);
for ii=1:m
xv(ii,:)=x([((ii-1)+1):1:((ii-1)+M)]);
end
y=0;
xt=xv;
for ii=2:M
xt=circshift(xt,);
y=y+sum(Heaviside(r-max(abs(xv-xt),[],1)));
end
y=y/M/(M-1);
function y=Heaviside(x)
y=zeros(size(x));
y(find(x>=0))=1;
3． 关于小数据量法
a. delta_t的意义
delta_t是采样时间，在帖子中看到有取值为1也有其他值，书上也没有太明确的说明。可否将其理解为是计算微分方程时的时间步长step呢，另外书中提到“只要delta_t小到能确定一个吸引子的轨道上的最小点的数目，大约在n到10n之间”就能得到满意结果，是否可以认为delta_t越小越好。
b. 小数据量法程序
在帖子中主要找到两个程序，分别是lambda_1=largest_lyapunov_exponent(data,N,m,tau,P)以及lambda_1=largest_lyapunov_p3(data,N,m,tau,P,delt_t)，其中应用前者计算一组从帖子中得到的数据耗时大约为2秒；而后者在10分钟后依然提示“busy”。将二者贴出，希望指出是哪个环节导致如此：
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
load data.txt;
N=length(data);
P=38;
m=5;
tau=3;
delt_t=1;
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)平均
=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);
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
function lambda_1=largest_lyapunov_p3(data,N,m,tau,P,delt_t)
m=8;tau=9;
delt_t=1;
load data.txt
%**********************************************************
x=data;
xPower=abs(fft(x)).^2;
NN=length(xPower);
xPower(1)=[];%去除直流分量
NN=floor(NN/2);
xPower=xPower(1:NN);
freq=(1:NN)/NN*0.5;
=max(xPower);
P=index;
%*************************************************************
d_length=[];
d_content=[];
Y=reconstitution(data,N,m,tau);
M=N-(m-1)*tau;
idx_j=0;
for j=1:M
d_min=10000;
for jj=1:M
d_s=0; %寻找相空间中每个点的最近距离点，并记下该点下标
if abs(j-jj)>P %限制短暂分离
d_s=sum(abs(Y(j)-Y(jj)));
if d_s<d_min
d_min=d_s;
idx_j=jj;
end
end
end
if ((M-j)>(M-idx_j));%计算点j的最大演化时间步长i
max_i=M-idx_j;
else
max_i=M-j;
end
d_length=;
for k=1:max_i %计算点j与其最近邻点在i个离散步后的距离
d_j_i=0;
d_j_i=norm(Y(j+k)-Y(idx_j+k));
d_content=;
end
end
%对每个演化时间步长i，求所有的j的lnd(i,j)平均
y=[];
for i=1:max(d_length)
S_j_i=0;
sum_former=0;
Count=0;
for j=1:M
if j==1
former=0;
else
former=d_length(j-1);
end
sum_former=sum_former+former;
if i<=(d_length(j))
if d_content(sum_former+i)~=0
S_j_i=S_j_i+log(d_content(sum_former+i));
Count=Count+1;
end
end
end
y=; %对每个演化时间步长i，求所有的j平均
end
XX=1:length(y);
figure;
plot(XX,y(XX),'-','markersize',1);hold on;%——————————【上面的图是从这儿画出来的】
linear=input('请输入线形部分的长度');
XX1=1:linear;
pp=polyfit(XX1,y(XX1),1);
lambda_1=pp(1)
yp=polyval(pp,XX1);
figure;
plot(XX1,yp,'r--');
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
c. 最小二乘法数据拟合
一些帖子中提到在进行数据拟合时，要选用线性部分较好的一段计算。那么从我的观点理解，判断好坏最直接的办法就是画图观察。而这样就出现了一个问题：每次都观察的话，如果绘制参数变化的最大LE曲线，岂不是循环计算无法实施，只能看一次重新拟合计算一次。也许我的理解有误，或者大家有更好的方法，希望指点。
帖子有些长，但这是我在看了许多参考资料后总结出来的问题，也许有一些很简单，但提出来至少也说明我用心了，相信论坛中也有很多像我一样对这些问题不甚了解的同学。希望各位前辈和同学给予指点，谢谢。
回复 4 # chunshui2003 的帖子
首先，我很佩服你的这种学习精神，另外，也非常意外的看到你在前面能提到我，呵呵。
至于你说的LE的计算，我大概看了一下你的帖子，你好像是在做连续系统的理论分析时需要计算Lyapunov指数吧，如果是这样为什么不用Jacbain方法呢。你所提到的方法应该是需要相空间重构，而你已经有了控制系统的微分方程，这样的话我认为就没必要解出时间序列再重构厢空间了，我建议你试试Jacbain的方法，根据LE的定义，线性化后可以得到在某个参数下的所有LE，进而也方便编程，实现参数变化时的最大LE曲线。 回复 5 # hsfy919 的帖子
再次感谢hsfy919的回答。之前参看oct的帖子后，下载了一个let的工具箱，里面用的就是Jacobian法（应该是这样吧）。之后oct说了小数据量法的优点，并且我也看了吕金虎老师的书，都说这种方法不错，精度也可以。
我的系统运动方程包括x1,y1,x2,y2四个自由度，换算成系统运动微分方程是8个。我不知道你的研究中是否有油膜力或者不平衡磁拉力等因素，因为这些表达式非常复杂，在求解Jacobian矩阵的时候很复杂，所以一开始就没有考虑这种方法。
不过经过你的提醒后我决定再把这个方法好好看一下，然后编程试试。
非常感谢你的帮助。 回复 6 # chunshui2003 的帖子
可能咱俩的研究有相似之处，我做的也含有非线性油膜力，不过我是12个自由度，但效果也挺好，呵呵，所以你先试试吧 回复 7 # hsfy919 的帖子
hsfy919，你好。
按照你说的方法，我将非自治系统改为自治系统然后参照let工具箱中的思路确定了Q，但是在求解Jacobian矩阵的时候遇到了一个问题：
我的打算是 Jac =jacobian (expression )
系统运动微分方程中所采用的油膜力来自于 Adiletta, G., A.R. Guido, and C. Rossi, Chaotic motions of a rigid rotor in short journal bearings. Nonlinear Dynamics, 1996. 10(3): p. 251-269
其中的alpha中含有函数“sign”，如果是数值计算则没有问题，但是因为要确定Jacobian矩阵，所以将系统中的变量用符号表达式进行了处理，然而在应用时，却出现了提示错误“Undefined function or method 'sign' for input arguments of type 'sym'.”，导致后续工作无法进行。
因为你提到了也采用油膜力，而这个油膜力应用也较为广泛，所以想问一下你是如何处理的。或者在求解Jacobian矩阵的时候有没有其他更好的方法。
sign就是符号函数，可用if语句来实现。 回复 8 # chunshui2003 的帖子
oct说的对，你试试用if语句替换一下符号函数，改起来应该也不难。另外我不是用的这个油膜力表达式，是给予不同假设自己推导的表达式。 回复 9 # octopussheng 的帖子
谢谢oct的回答，昨天rocwoods（吴鹏老师）也给予了帮助。吴老师的意见是：
sign(x) 等价于 x./(abs(x) + (x==0))，我在matlab中试验了一下，Jacobian矩阵可以求出，但是表达式过于复杂，最终提示超出屏幕显示范围（25000字节）。看来这种方法还是不行，我再试一下。 回复 10 # hsfy919 的帖子
谢谢你的回答，我今天再试一下，尝试不同的方法。 回复 11 # chunshui2003 的帖子
其实Jacobian矩阵也是数值计算，没有必要一定要得到它的符号表达。
只要定义微分方程的程序能正确求解，问题就解决了。 回复 13 # octopussheng 的帖子
oct学长：
你好。
我是在看了你的总结帖并且参照LET工具箱之后才这样求解的。因为里面是求出了运动微分方程的Jacobian表达式之后才进行LE的计算的。
可能我对这个求解方法的理解还不够，对于学长所提到的“只要定义微分方程的程序能正确求解，问题就解决了。”这句话的意思不是很理解。
在帖子http://forum.vibunion.com/forum-viewthread-tid-49671-highlight-Lyapunov.html中
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;
这部分内容是不是就是“微分方程的定义程序”呢。但里面依然有Jaco的求解，如果按照学长的说法，不进行符号预算，好像是无法求出Jacobian矩阵的。或者说我的理解有问题。
麻烦学长解释一下，谢谢。
回复 13 # octopussheng 的帖子
另外，学长，将我的定义微分方程列出来，你看一下过程是否有误：
functionuu = equ_elastic(t,u)
m1 = 60; %转子质量 kg
m2 = 25; %轴颈质量 kg
c1 = 4000; %转子阻尼 N.s/m
c2 = 1200; %轴承阻尼 N.s/m
Ke = 6.2*10^6; %转轴刚度 N/m
e0 = 0.6*10^-3; %转子质量偏心 m
Rr = 60*10^-3; %转子半径 m
Lr = 150*10^-3; %转子长 m
delta_0 = 4.5*10^-3; %均匀气隙大小 m
cz = 0.2*10^-3; %轴颈间隙 m
Rb = 50*10^-3; %轴承半径 m
Lb = 20*10^-3; %轴承长 m
mu = 18*10^-3; %绝对润滑油粘度 无单位
mu_0 = 4*pi*10^-7; %空气磁导系数 H/m
kj = 5.2; %气隙基波磁动势系数 无单位
Ij = 4; %励磁电流
Fj = kj^2 * Ij^2; %励磁电流的基波磁动势幅值
omega = 13; %大轴旋转角速度
%%% 不平衡磁拉力表达式 FxFy %%%%
e = sqrt(u(1).^2 + u(3).^2); %转子轴心轨迹
F_coeffi = Rr * Lr *pi *mu_0 *Fj^2 /delta_0.^2;
Fx = F_coeffi * (1/2 * e + 5/8 * e.^3) * u(1)./e;
Fy = F_coeffi * (1/2 * e + 5/8 * e.^3) * u(3)./e; %不平衡磁拉力
%%% 非线性油膜力表达式 fxfy%%%%
sigma = mu * omega * Rb * Lb * (Rb/cz)^2 * (Lb/2/Rb)^2; %Sommerfeld修正系数
A1 = u(7) + 2*u(6);
A2 = u(5) - 2*u(8);
alpha = atan(A1./A2) - pi/2* sign(A1./A2) - pi/2* sign(A1);
E = sqrt(u(5).^2 +u(7).^2); %轴颈轴心轨迹
E_minus = sqrt(1 -E.^2);
B1 = u(7)*cos(alpha) - u(5)*sin(alpha);
B2 = u(5)*cos(alpha) + u(7)*sin(alpha);
G = 2./E_minus * ( pi/2 + atan(B1./ E_minus) );
V = (2 + B1 * G) ./E_minus.^2;
S = B2 ./ ( 1 - B2.^2 );
F_oil_same = -sqrt( A2.^2 + A1.^2 ) ./ E_minus.^2;
f_x = F_oil_same * ( 3* u(5) *V - sin(alpha) * G - 2 *cos(alpha)*S );
f_y = F_oil_same * ( 3* u(7) *V + cos(alpha) * G - 2 *sin(alpha)*S ); %油膜力的无量纲表达式
fx = sigma * f_x;
fy = sigma * f_y; %非线性油膜力
Q = [u(10) u(19) u(28) u(37) u(46) u(55) u(64) u(73) u(82);
u(11) u(20) u(29) u(38) u(47) u(56) u(65) u(74) u(83);
u(12) u(21) u(30) u(39) u(48) u(57) u(66) u(75) u(84);
u(13) u(22) u(31) u(40) u(49) u(58) u(67) u(76) u(85);
u(14) u(23) u(32) u(41) u(50) u(59) u(68) u(77) u(86);
u(15) u(24) u(33) u(42) u(51) u(60) u(69) u(78) u(87);
u(16) u(25) u(34) u(43) u(52) u(61) u(70) u(79) u(88);
u(17) u(26) u(35) u(44) u(53) u(62) u(71) u(80) u(89);
u(18) u(27) u(36) u(45) u(54) u(63) u(72) u(81) u(90);
];
uu = zeros(90,1);
uu(1) = u(2);
uu(2) = -c1/(m1*omega) * u(2) - Ke/(m1*omega^2) * u(1) + Ke * cz/(m1*omega^2*delta_0)*u(5) + e0 * cos(u(9))/delta_0 + Fx/(m1*omega^2*delta_0);
uu(3) = u(4);
uu(4) = -c1/(m1*omega) * u(4) - Ke/(m1*omega^2) * u(3) + Ke * cz/(m1*omega^2*delta_0)*u(7) + e0 * sin(u(9))/delta_0 + Fy/(m1*omega^2*delta_0);
uu(5) = u(6);
uu(6) = -c2/(m2*omega) * u(6) + Ke * delta_0/(2*m2*omega^2*cz) * u(1) - Ke/(2*m2*omega^2) * u(5) + fx/(m2*omega^2 *cz);
uu(7) = u(8);
uu(8) = -c2/(m2*omega) * u(8) + Ke * delta_0/(2*m2*omega^2*cz) * u(3) - Ke/(2*m2*omega^2) * u(7) + fy/(m2*omega^2 *cz);
uu(9) = 1;
J = [];
uu(10:90) = J*Q;
这里Jacobian矩阵没有列出来，因为实在是太长了，尤其涉及到油膜力求导的部分已经超过了matlab屏幕显示范围。
页:
[1]
2