设为首页收藏本站联系我们

声振论坛

 找回密码
 我要加入

QQ登录

只需一步，快速开始

声振论坛 › 展示 ›振动理论› 非线性振动 › 查看内容

关于最大Lyapunov指数计算的几个疑问

2015-10-23 10:51| 发布者: aspen| 查看: 1616| 评论: 25|原作者: chunshui2003|来自: 声振论坛

摘要: 首先要感谢无水、oct、咕噜噜、gghhjj和liliangbiao等在求解非线性微分方程、制作poincare映射图和分岔图过程中给予的帮助，虽然有的问题没有在论坛中提问，但前辈的帖子和解释让我受益匪浅；其次，感谢 vickyvictor ...
首先要感谢无水、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求解方程，如[t,u] = 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分钟）

data.txt

【本文来自声振论坛热议帖子，详细内容见后续评论】
收藏 分享 邀请
上一篇：非自治系统分岔图绘制实例下一篇：庞加莱截面图绘制matlab代码
发表评论

最新评论

引用 chunshui2003 2010-11-28 16:39
本帖最后由 chunshui2003 于 2010-11-28 16:41 编辑

不知道为什么程序会变成那个样子，重新贴一下,按照yufeiqun2008的帖子http://forum.vibunion.com/forum-viewthread-tid-57387-highlight-Lyapunov.html：


  3. function [St,S_t,dS_m_t,dS_t,S_cor]=CC(x)
  4. %% C-C方法提取延迟时间和延迟时间窗口
  5. % by Yu Feiqun
  6. % 2007-12-09
  7. N=length(x);
  8. % 计算给定序列的标准差 deta
  9. deta=std(x);
  10. tmax=200;
  11. t=1:tmax;
  12. m=2:5;
  13. r=1:4;
  14. St=zeros(length(t),length(m),length(r));
  15. %dSt=zeros(length(t),length(m),length(r));
  16. for tt=1:length(t)
  17.     % 计算时间序列分为t个不相交的时间序列时，每一时间序列的长度
  18.     L=floor(N/t(tt));
  19.     for mm=1:length(m)
  20.         for rr=1:length(r)
  21.             rjj=deta*r(rr)/2;
  22.             smrt=0;
  23.             for ss=1:t(tt)
  24.                 % 提取第ss个子序列
  25.                 xsub=x([ss:t(tt):N]);
  26.                 xsub=xsub([1:L]);
  27.                 c1=getCValue(xsub,m(mm),rjj,t(tt));
  28.                 c2=getCValue(xsub,1,rjj,t(tt));
  29.                 smrt=smrt+c1-c2.^m(mm);
  30.                
  31.             end
  32.             [tt,mm,rr]
  33.             St(tt,mm,rr)=smrt/t(tt);
  34.         end
  35.     end
  36. end
  37. S_t=sum(sum(St,3),2)/16;
  38. dS_m_t=max(St,[],3)-min(St,[],3);
  39. dS_t=sum(dS_m_t,2)/4;
  40. S_cor=abs(S_t)+dS_t;


  43. function y=getCValue(x,m,r,t)
  44. % 计算时间序列x的关联积分
  45. % m为嵌入维数 t为时间延迟 N为序列长度 r为半径
  46. % 由于序列x已经经过划分,其相邻两个数据已具有时间延迟t
  47. N=length(x);
  48. M=N-(m-1);
  49. % 重构相空间 嵌入维数为m 时间延迟为t
  50. xv=zeros(m,M);
  51. for ii=1:m
  52.     xv(ii,:)=x([((ii-1)+1):1:((ii-1)+M)]);
  53. end
  54. y=0;
  55. xt=xv;
  56. for ii=2:M
  57.     xt=circshift(xt,[0,-1]);
  58.     y=y+sum(Heaviside(r-max(abs(xv-xt),[],1)));
  59. end
  60. y=y/M/(M-1);
  61.         
  62. function y=Heaviside(x)
  63. y=zeros(size(x));   
  64. y(find(x>=0))=1;

复制代码
引用 chunshui2003 2010-11-28 16:44
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”。将二者贴出，希望指出是哪个环节导致如此：




  5. function lambda_1=largest_lyapunov_exponent(data,N,m,tau,P)

  7. %the function is used to calcultate largest lyapunov exponent with the

  9. %mended algorithm,which put forward by lv jing hu.

  11. %data:the time series

  13. %N:the length of data

  15. %m:enbedding dimention

  17. %tau:time delay

  19. %P:the mean period of the time series,calculated with FFT

  21. %lambda_1:return the largest lyapunov exponent

  23. %skyhawk

  25. load data.txt;



  29. N=length(data);

  31. P=38;

  33. m=5;

  35. tau=3;

  37. delt_t=1;

  39. Y=reconstitution(data,N,m,tau );%reconstitute state space

  41. M=N-(m-1)*tau;%M is the number of embedded points in m-dimensional space

  43. for j=1:M

  45.     d_max=1e+100;

  47.     for jj=1:M                                              %寻找相空间中每个点的最近距离点，并记下

  49.         d_s=0;                                              %该点下标

  51.         if abs(j-jj)>P                                      %限制短暂分离

  53.             for i=1:m

  55.                 d_s=d_s+(Y(i,j)-Y(i,jj))*(Y(i,j)-Y(i,jj));

  57.                 d_min=d_max;

  59.                 if d_s<d_min

  61.                    d_min=d_s;

  63.                    idx_j=jj;

  65.                end

  67.             end

  69.         end

  71.     end

  73. %     index(j)=idx_j;

  75.     max_i=min((M-j),(M-idx_j));%计算点j的最大演化时间步长i

  77.     for k=1:max_i              %计算点j与其最近邻点在i个离散步后的距离

  79.         d_j_i=0;

  81.         for kk=1:m

  83.             d_j_i=d_j_i+(Y(kk,j+k)-Y(kk,idx_j+k))*(Y(kk,j+k)-Y(kk,idx_j+k));

  85.             d(k,j)=d_j_i;

  87.         end

  89.     end

  91. end



  95. %对每个演化时间步长i，求所有的j的lnd(i,j)平均

  97. [l_i,l_j]=size(d);

  99. for i=1:l_i

  101.     q=0;

  103.     y_s=0;

  105.     for j=1:l_j

  107.         if d(i,j)~=0

  109.             q=q+1;

  111.             y_s=y_s+log(d(i,j));

  113.         end

  115.     end

  117.     y(i)=y_s/(q*delt_t);

  119. end

  121. x=1:length(y);

  123. y

  125. pp=polyfit(x,y,1);

  127. lambda_1=pp(1);

  129. yp=polyval(pp,x);

  131. plot(x,y,'-o',x,yp,'--')

  133. function X=reconstitution(data,N,m,tau)
  134. %该函数用来重构相空间
  135. % m为嵌入空间维数
  136. % tau为时间延迟
  137. % data为输入时间序列
  138. % N为时间序列长度
  139. % X为输出,是m*n维矩阵

  141. M=N-(m-1)*tau;%相空间中点的个数
  142. for j=1:M           %相空间重构
  143.     for i=1:m
  144.         X(i,j)=data((i-1)*tau+j);
  145.     end
  146. end

复制代码
引用 chunshui2003 2010-11-28 16:45


  3. function lambda_1=largest_lyapunov_p3(data,N,m,tau,P,delt_t)
  4. m=8;tau=9;
  5. delt_t=1;
  6. load data.txt
  7. %**********************************************************
  8. x=data;
  9. xPower=abs(fft(x)).^2;
  10. NN=length(xPower);
  11. xPower(1)=[];%去除直流分量
  12. NN=floor(NN/2);
  13. xPower=xPower(1:NN);
  14. freq=(1:NN)/NN*0.5;
  15. [mP,index]=max(xPower);
  16. P=index;  
  17. %*************************************************************

  19. d_length=[];
  20. d_content=[];
  21. Y=reconstitution(data,N,m,tau);
  22. M=N-(m-1)*tau;
  23. idx_j=0;
  24. for j=1:M
  25.     d_min=10000;
  26.     for jj=1:M                                             
  27.         d_s=0;      %寻找相空间中每个点的最近距离点，并记下该点下标
  28.         if abs(j-jj)>P                                      %限制短暂分离
  29.             d_s=sum(abs(Y(j)-Y(jj)));
  30.             if d_s<d_min
  31.                 d_min=d_s;
  32.                idx_j=jj;
  33.             end
  34.         end
  35.     end
  36.     if ((M-j)>(M-idx_j));%计算点j的最大演化时间步长i
  37.         max_i=M-idx_j;
  38.     else
  39.         max_i=M-j;
  40.     end
  41.     d_length=[d_length,max_i];
  42.     for k=1:max_i              %计算点j与其最近邻点在i个离散步后的距离
  43.         d_j_i=0;        
  44.         d_j_i=norm(Y(j+k)-Y(idx_j+k));
  45.         d_content=[d_content,d_j_i];
  46.     end

  48. end

  50. %对每个演化时间步长i，求所有的j的lnd(i,j)平均
  51. y=[];
  52. for i=1:max(d_length)
  53.    S_j_i=0;
  54.    sum_former=0;
  55.    Count=0;
  56.    for j=1:M
  57.        if j==1
  58.           former=0;
  59.        else
  60.           former=d_length(j-1);
  61.        end   
  62.        sum_former=sum_former+former;
  63.        if i<=(d_length(j))
  64.            if d_content(sum_former+i)~=0
  65.               S_j_i=S_j_i+log(d_content(sum_former+i));
  66.               Count=Count+1;
  67.            end
  68.        end
  69.    end
  70.    y=[y,S_j_i/(Count*delt_t)]; %对每个演化时间步长i，求所有的j平均
  71. end
  72. XX=1:length(y);
  73. figure;
  74. plot(XX,y(XX),'-','markersize',1);hold on;%——————————【上面的图是从这儿画出来的】
  75. linear=input('请输入线形部分的长度');
  76. XX1=1:linear;
  77. pp=polyfit(XX1,y(XX1),1);
  78. lambda_1=pp(1)
  79. yp=polyval(pp,XX1);
  80. figure;
  81. plot(XX1,yp,'r--');


  84. function X=reconstitution(data,N,m,tau)
  85. %该函数用来重构相空间
  86. % m为嵌入空间维数
  87. % tau为时间延迟
  88. % data为输入时间序列
  89. % N为时间序列长度
  90. % X为输出,是m*n维矩阵

  92. M=N-(m-1)*tau;%相空间中点的个数
  93. for j=1:M           %相空间重构
  94.     for i=1:m
  95.         X(i,j)=data((i-1)*tau+j);
  96.     end
  97. end

复制代码
c.        最小二乘法数据拟合

一些帖子中提到在进行数据拟合时，要选用线性部分较好的一段计算。那么从我的观点理解，判断好坏最直接的办法就是画图观察。而这样就出现了一个问题：每次都观察的话，如果绘制参数变化的最大LE曲线，岂不是循环计算无法实施，只能看一次重新拟合计算一次。也许我的理解有误，或者大家有更好的方法，希望指点。

帖子有些长，但这是我在看了许多参考资料后总结出来的问题，也许有一些很简单，但提出来至少也说明我用心了，相信论坛中也有很多像我一样对这些问题不甚了解的同学。希望各位前辈和同学给予指点，谢谢。
引用 hsfy919 2010-11-28 22:25
回复 4 # chunshui2003 的帖子

首先，我很佩服你的这种学习精神，另外，也非常意外的看到你在前面能提到我，呵呵。
至于你说的LE的计算，我大概看了一下你的帖子，你好像是在做连续系统的理论分析时需要计算Lyapunov指数吧，如果是这样为什么不用Jacbain方法呢。你所提到的方法应该是需要相空间重构，而你已经有了控制系统的微分方程，这样的话我认为就没必要解出时间序列再重构厢空间了，我建议你试试Jacbain的方法，根据LE的定义，线性化后可以得到在某个参数下的所有LE，进而也方便编程，实现参数变化时的最大LE曲线。
引用 chunshui2003 2010-11-29 08:44
回复 5 # hsfy919 的帖子

再次感谢hsfy919的回答。之前参看oct的帖子后，下载了一个let的工具箱，里面用的就是Jacobian法（应该是这样吧）。之后oct说了小数据量法的优点，并且我也看了吕金虎老师的书，都说这种方法不错，精度也可以。
我的系统运动方程包括x1,y1,x2,y2四个自由度，换算成系统运动微分方程是8个。我不知道你的研究中是否有油膜力或者不平衡磁拉力等因素，因为这些表达式非常复杂，在求解Jacobian矩阵的时候很复杂，所以一开始就没有考虑这种方法。
不过经过你的提醒后我决定再把这个方法好好看一下，然后编程试试。
非常感谢你的帮助。
引用 hsfy919 2010-11-29 09:51
回复 6 # chunshui2003 的帖子

可能咱俩的研究有相似之处，我做的也含有非线性油膜力，不过我是12个自由度，但效果也挺好，呵呵，所以你先试试吧
引用 chunshui2003 2010-11-29 16:07
回复 7 # hsfy919 的帖子

hsfy919，你好。
按照你说的方法，我将非自治系统改为自治系统然后参照let工具箱中的思路确定了Q，但是在求解Jacobian矩阵的时候遇到了一个问题：
我的打算是 Jac =jacobian (expression    [u1 u2 u3 u4 u5 u6 u7 u8 u9])

系统运动微分方程中所采用的油膜力来自于   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

alpha

oil force

oil force


其中的alpha中含有函数“sign”，如果是数值计算则没有问题，但是因为要确定Jacobian矩阵，所以将系统中的变量用符号表达式进行了处理，然而在应用时，却出现了提示错误“Undefined function or method 'sign' for input arguments of type 'sym'.”，导致后续工作无法进行。

因为你提到了也采用油膜力，而这个油膜力应用也较为广泛，所以想问一下你是如何处理的。或者在求解Jacobian矩阵的时候有没有其他更好的方法。
引用 octopussheng 2010-11-30 21:51
sign就是符号函数，可用if语句来实现。
引用 hsfy919 2010-12-1 00:16
回复 8 # chunshui2003 的帖子

oct说的对，你试试用if语句替换一下符号函数，改起来应该也不难。另外我不是用的这个油膜力表达式，是给予不同假设自己推导的表达式。
引用 chunshui2003 2010-12-1 08:03
回复 9 # octopussheng 的帖子

谢谢oct的回答，昨天rocwoods（吴鹏老师）也给予了帮助。吴老师的意见是：
sign(x)   等价于 x./(abs(x) + (x==0))  ，我在matlab中试验了一下，Jacobian矩阵可以求出，但是表达式过于复杂，最终提示超出屏幕显示范围（25000字节）。看来这种方法还是不行，我再试一下。
引用 chunshui2003 2010-12-1 08:06
回复 10 # hsfy919 的帖子

谢谢你的回答，我今天再试一下，尝试不同的方法。
引用 octopussheng 2010-12-1 12:17
回复 11 # chunshui2003 的帖子

其实Jacobian矩阵也是数值计算，没有必要一定要得到它的符号表达。

只要定义微分方程的程序能正确求解，问题就解决了。
引用 chunshui2003 2010-12-1 15:38
回复 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 a  0;
        z 0  x-c];
dX(4:12) = Jaco*Y;

这部分内容是不是就是“微分方程的定义程序”呢。但里面依然有Jaco的求解，如果按照学长的说法，不进行符号预算，好像是无法求出Jacobian矩阵的。或者说我的理解有问题。
麻烦学长解释一下，谢谢。


引用 chunshui2003 2010-12-1 16:04
回复 13 # octopussheng 的帖子


另外，学长，将我的定义微分方程列出来，你看一下过程是否有误：



  3. function  uu = equ_elastic(t,u)
  4. m1 = 60;               %转子质量 kg
  5.   m2 = 25;               %轴颈质量 kg
  6.   c1 = 4000;             %转子阻尼 N.s/m
  7.   c2 = 1200;             %轴承阻尼 N.s/m
  8.   Ke = 6.2*10^6;          %转轴刚度 N/m
  9.   
  10.   e0 = 0.6*10^-3;        %转子质量偏心 m
  11.   
  12.   Rr = 60*10^-3;         %转子半径 m
  13.   Lr = 150*10^-3;        %转子长 m
  14.   delta_0 = 4.5*10^-3;   %均匀气隙大小 m
  15.   
  16.   cz = 0.2*10^-3;        %轴颈间隙 m
  17.   
  18.   Rb = 50*10^-3;         %轴承半径 m
  19.   
  20.       
  21.   Lb = 20*10^-3;          %轴承长 m
  22.   
  23.   mu = 18*10^-3;         %绝对润滑油粘度 无单位
  24.   mu_0 = 4*pi*10^-7;     %空气磁导系数 H/m  
  25.   kj = 5.2;              %气隙基波磁动势系数 无单位
  26.   Ij = 4;            %励磁电流
  27.   Fj = kj^2 * Ij^2;              %励磁电流的基波磁动势幅值
  28.   
  29.   omega = 13;            %大轴旋转角速度

  31.   
  32.   
  33.   %%% 不平衡磁拉力表达式 Fx  Fy   %%%%
  34.   e = sqrt(u(1).^2 + u(3).^2);                           %转子轴心轨迹
  35.   F_coeffi = Rr * Lr *pi *mu_0 *Fj^2 /delta_0.^2;
  36.   Fx = F_coeffi * (1/2 * e + 5/8 * e.^3) * u(1)./e;
  37.   Fy = F_coeffi * (1/2 * e + 5/8 * e.^3) * u(3)./e;       %不平衡磁拉力
  38.   
  39.   %%% 非线性油膜力表达式 fx  fy  %%%%
  40.   
  41.   sigma = mu * omega * Rb * Lb * (Rb/cz)^2 * (Lb/2/Rb)^2;        %Sommerfeld修正系数
  42.   
  43.   A1 = u(7) + 2*u(6);
  44.   A2 = u(5) - 2*u(8);
  45.   
  46.   alpha = atan(A1./A2) - pi/2* sign(A1./A2) - pi/2* sign(A1);
  47.   
  48. E = sqrt(u(5).^2 +u(7).^2);                     %轴颈轴心轨迹
  49. E_minus = sqrt(1 -E.^2);

  51. B1 = u(7)*cos(alpha) - u(5)*sin(alpha);
  52. B2 = u(5)*cos(alpha) + u(7)*sin(alpha);
  53.   
  54.   G = 2./E_minus * ( pi/2 + atan(B1./ E_minus) );
  55.   V = (2 + B1 * G) ./  E_minus.^2;
  56.   S = B2 ./ ( 1 - B2.^2 );
  57.   
  58.   F_oil_same = -sqrt( A2.^2 + A1.^2 ) ./ E_minus.^2;
  59.   
  60.   f_x = F_oil_same * ( 3* u(5) *V - sin(alpha) * G - 2 *cos(alpha)*S );
  61.   f_y = F_oil_same * ( 3* u(7) *V + cos(alpha) * G - 2 *sin(alpha)*S );     %油膜力的无量纲表达式
  62.   
  63.   fx = sigma * f_x;
  64.   fy = sigma * f_y;              %非线性油膜力
  65.   
  66.   
  67.   
  68.   Q = [u(10) u(19) u(28) u(37) u(46) u(55) u(64) u(73) u(82);
  69.       u(11) u(20) u(29) u(38) u(47) u(56) u(65) u(74) u(83);
  70.       u(12) u(21) u(30) u(39) u(48) u(57) u(66) u(75) u(84);
  71.       u(13) u(22) u(31) u(40) u(49) u(58) u(67) u(76) u(85);
  72.       u(14) u(23) u(32) u(41) u(50) u(59) u(68) u(77) u(86);
  73.       u(15) u(24) u(33) u(42) u(51) u(60) u(69) u(78) u(87);
  74.       u(16) u(25) u(34) u(43) u(52) u(61) u(70) u(79) u(88);
  75.       u(17) u(26) u(35) u(44) u(53) u(62) u(71) u(80) u(89);
  76.       u(18) u(27) u(36) u(45) u(54) u(63) u(72) u(81) u(90);
  77.   ];


  80. uu = zeros(90,1);

  82.   uu(1) = u(2);
  83.   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);
  84.   uu(3) = u(4);
  85.   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);
  86.   uu(5) = u(6);
  87.   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);
  88. uu(7) = u(8);
  89.   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);
  90.   uu(9) = 1;


  93. J = [];


  96. uu(10:90) = J*Q;

复制代码



这里Jacobian矩阵没有列出来，因为实在是太长了，尤其涉及到油膜力求导的部分已经超过了matlab屏幕显示范围。

引用 octopussheng 2010-12-1 20:54
确实，因为有sign函数的存在，Jacobian矩阵的推导很麻烦。

不过没办法，要用这个方法的话，就必须得推导Jacobian矩阵了。
耐心点吧，呵呵！
引用 hsfy919 2010-12-1 23:39
回复 15 # chunshui2003 的帖子

你试试不要让matlab显示Jacbain矩阵，直接传递变量，调用Jacbain程序。
另外你的轴承参数是实际参数吗？单位是多少啊？想借鉴一下
引用 chunshui2003 2010-12-2 09:18
回复 17 # hsfy919 的帖子

恩，你说的这个方法我倒是可以试试。另外，由于sign函数最后确定的是数值而不是变量，所以我做了一个假设，将

  2. alpha = atan(A1./A2) - pi/2* sign(A1./A2) - pi/2* sign(A1)
复制代码

直接简化为

  2. alpha = atan(A1./A2)
复制代码


其余参数保持不变，这样通过let工具箱得到的一个LE曲线如图所示：

最大LE

最大LE



我觉得这应该是不对的，因为在这个参数条件下，得到的转子轴心轨迹和poincare映射图显示为拟周期运动，最大LE应该小于0我想alpha改变的影响应该不至于这么大。

轨迹图

轨迹图

另外，参数的选取是根据上海和西安交通大学的老师和前辈文章中的计算和实验参数分析比较后得到，应该很可信。单位在后面都有标明的。具体可以参照
Cheng, M., G. Meng, and J. Jing, Numerical and experimental study of a rotor-bearing-seal system. Mechanism and Machine Theory, 2007. 42(8): p. 1043-1057.
以及
Shen, X. and M. Zhao, Effect of the Seal Force on Nonlinear Dynamics and Stability of the Rotor-bearing-seal System. Journal of Vibration and Control, 2009. 15(2): p. 197-217.
and
李松涛, et al. , 迷宫密封不平衡转子动力系统的稳定性与分岔. 应用数学和力学, 2003(11).

poincare图

poincare图














引用 chunshui2003 2010-12-2 09:24
回复 16 # octopussheng 的帖子

otc前辈，如果这个方法效果不是很理想的话，是不是就要采用其他诸如小数据量法了呢。那么又回到之前的问题：
1数据是计算后得到的连续点还是周期点
2延迟时间tau和嵌入维数m，对于同一个微分方程不同参数得到的不同数据，是不是每次通过CC法计算均会得到不同的数值
最后，请问学长是采用什么样的方法进行LE的计算的。
引用 octopussheng 2010-12-3 08:24
小数据量法可以参考陆振波博士的一些工作，他做的不错。

1）数据要进行采样；
2）时间延迟和嵌入维数，是需要根据相应程序计算的。

查看全部评论(25)

相关分类

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-26 12:37 , Processed in 0.073790 second(s), 26 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部