杭州锐达数字技术有限公司
查看: 5953|回复: 26

[稳定性与分岔] 关于最大Lyapunov指数计算的几个疑问(非重复提问)

[复制链接]
发表于 2010-11-28 16:37 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
本帖最后由 chunshui2003 于 2010-11-28 16:45 编辑

首先要感谢无水、oct、咕噜噜、gghhjjliliangbiao等在求解非线性微分方程、制作poincare映射图和分岔图过程中给予的帮助,虽然有的问题没有在论坛中提问,但前辈的帖子和解释让我受益匪浅;其次,感谢 vickyvictoriahsfy919在进一步进行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的方法。看了gghhjjoctyufeiqun2008等前辈的帖子,现在有如下几个问题不是很清楚,希望得到解答:

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,其代码如下:



结合吕金虎老师的书,得知延迟时间taudS_t vs tmax图中第一个极小值点所对应的t值;而S_cor vs tmax 图中最小值点对应的t_min为时间窗口,嵌入维数按照网友的说法为m = t_min*delta_t / tau。(其中delta_t为采样时间,这个我也不是很清楚,在后面有说明)。不知道这样的计算taum方法是否正确。
b.
是否只有通过dS_t vs tmax S_cor vs tmax图才能确认taum
也看了其他的CC法帖子,但上面都没有给出直接的taum计算表达式,而是通过图形观察得到。那么,是不是意味着如果采用CC法,延迟时间和嵌入维数只能通过观察图形的方式获得?
c.
数据不同是否taum就不同
微分方程相同,如果更改了参数而得到不同的数据,那么是不是taum就不同了呢。如果相同,那么在采用小数据量法时,taum可以采用同一值;如果不同的话,是不是每次都要计算tau和m,或是有其他的方法。(在我的电脑上计算一次大概需要耗时30分钟)

data.txt

46.88 KB, 下载次数: 56

论坛中的计算数据

点评

赞成: 3.0
赞成: 3
yufeiqun2008的帖子打不开,而且“代码”也不见了……  发表于 2015-3-15 15:22

评分

1

查看全部评分

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

 楼主| 发表于 2010-11-28 16:39 | 显示全部楼层
本帖最后由 chunshui2003 于 2010-11-28 16:41 编辑

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


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


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

复制代码
 楼主| 发表于 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”。将二者贴出,希望指出是哪个环节导致如此:




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

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

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

  4. %data:the time series

  5. %N:the length of data

  6. %m:enbedding dimention

  7. %tau:time delay

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

  9. %lambda_1:return the largest lyapunov exponent

  10. %skyhawk

  11. load data.txt;



  12. N=length(data);

  13. P=38;

  14. m=5;

  15. tau=3;

  16. delt_t=1;

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

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

  19. for j=1:M

  20.     d_max=1e+100;

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

  22.         d_s=0;                                              %该点下标

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

  24.             for i=1:m

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

  26.                 d_min=d_max;

  27.                 if d_s<d_min

  28.                    d_min=d_s;

  29.                    idx_j=jj;

  30.                end

  31.             end

  32.         end

  33.     end

  34. %     index(j)=idx_j;

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

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

  37.         d_j_i=0;

  38.         for kk=1:m

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

  40.             d(k,j)=d_j_i;

  41.         end

  42.     end

  43. end



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

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

  46. for i=1:l_i

  47.     q=0;

  48.     y_s=0;

  49.     for j=1:l_j

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

  51.             q=q+1;

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

  53.         end

  54.     end

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

  56. end

  57. x=1:length(y);

  58. y

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

  60. lambda_1=pp(1);

  61. yp=polyval(pp,x);

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

  63. function X=reconstitution(data,N,m,tau)
  64. %该函数用来重构相空间
  65. % m为嵌入空间维数
  66. % tau为时间延迟
  67. % data为输入时间序列
  68. % N为时间序列长度
  69. % X为输出,是m*n维矩阵

  70. M=N-(m-1)*tau;%相空间中点的个数
  71. for j=1:M           %相空间重构
  72.     for i=1:m
  73.         X(i,j)=data((i-1)*tau+j);
  74.     end
  75. end

复制代码
 楼主| 发表于 2010-11-28 16:45 | 显示全部楼层


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

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

  44. end

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


  77. function X=reconstitution(data,N,m,tau)
  78. %该函数用来重构相空间
  79. % m为嵌入空间维数
  80. % tau为时间延迟
  81. % data为输入时间序列
  82. % N为时间序列长度
  83. % X为输出,是m*n维矩阵

  84. M=N-(m-1)*tau;%相空间中点的个数
  85. for j=1:M           %相空间重构
  86.     for i=1:m
  87.         X(i,j)=data((i-1)*tau+j);
  88.     end
  89. end

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

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

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

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

评分

1

查看全部评分

 楼主| 发表于 2010-11-29 08:44 | 显示全部楼层
回复 5 # hsfy919 的帖子

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

可能咱俩的研究有相似之处,我做的也含有非线性油膜力,不过我是12个自由度,但效果也挺好,呵呵,所以你先试试吧
 楼主| 发表于 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矩阵的时候有没有其他更好的方法。
发表于 2010-11-30 21:51 | 显示全部楼层
sign就是符号函数,可用if语句来实现。
发表于 2010-12-1 00:16 | 显示全部楼层
回复 8 # chunshui2003 的帖子

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

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

谢谢你的回答,我今天再试一下,尝试不同的方法。
发表于 2010-12-1 12:17 | 显示全部楼层
回复 11 # chunshui2003 的帖子

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

只要定义微分方程的程序能正确求解,问题就解决了。
 楼主| 发表于 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矩阵的。或者说我的理解有问题。
麻烦学长解释一下,谢谢。


 楼主| 发表于 2010-12-1 16:04 | 显示全部楼层
回复 13 # octopussheng 的帖子


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



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

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

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


  74. uu = zeros(90,1);

  75.   uu(1) = u(2);
  76.   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);
  77.   uu(3) = u(4);
  78.   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);
  79.   uu(5) = u(6);
  80.   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);
  81. uu(7) = u(8);
  82.   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);
  83.   uu(9) = 1;


  84. J = [];


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

复制代码



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

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表