关于最大Lyapunov指数计算的几个疑问(非重复提问)

chunshui2003

首先要感谢无水、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分钟）

chunshui2003

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

1. 
   
2. 
   
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;
    
41. 
    
42. 
    
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;
    
65. 
    

复制代码

chunshui2003

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. 
   
2. 
   
3. 
   
4. 
   
5. function lambda_1=largest_lyapunov_exponent(data,N,m,tau,P)
   
6. 
   
7. %the function is used to calcultate largest lyapunov exponent with the
   
8. 
   
9. %mended algorithm,which put forward by lv jing hu.
   
10. 
    
11. %data:the time series
    
12. 
    
13. %N:the length of data
    
14. 
    
15. %m:enbedding dimention
    
16. 
    
17. %tau:time delay
    
18. 
    
19. %P:the mean period of the time series,calculated with FFT
    
20. 
    
21. %lambda_1:return the largest lyapunov exponent
    
22. 
    
23. %skyhawk
    
24. 
    
25. load data.txt;
    
26. 
    
27. 
    
28. 
    
29. N=length(data);
    
30. 
    
31. P=38;
    
32. 
    
33. m=5;
    
34. 
    
35. tau=3;
    
36. 
    
37. delt_t=1;
    
38. 
    
39. Y=reconstitution(data,N,m,tau );%reconstitute state space
    
40. 
    
41. M=N-(m-1)*tau;%M is the number of embedded points in m-dimensional space
    
42. 
    
43. for j=1:M
    
44. 
    
45.     d_max=1e+100;
    
46. 
    
47.     for jj=1:M                                              %寻找相空间中每个点的最近距离点，并记下
    
48. 
    
49.         d_s=0;                                              %该点下标
    
50. 
    
51.         if abs(j-jj)>P                                      %限制短暂分离
    
52. 
    
53.             for i=1:m
    
54. 
    
55.                 d_s=d_s+(Y(i,j)-Y(i,jj))*(Y(i,j)-Y(i,jj));
    
56. 
    
57.                 d_min=d_max;
    
58. 
    
59.                 if d_s<d_min
    
60. 
    
61.                    d_min=d_s;
    
62. 
    
63.                    idx_j=jj;
    
64. 
    
65.                end
    
66. 
    
67.             end
    
68. 
    
69.         end
    
70. 
    
71.     end
    
72. 
    
73. %     index(j)=idx_j;
    
74. 
    
75.     max_i=min((M-j),(M-idx_j));%计算点j的最大演化时间步长i
    
76. 
    
77.     for k=1:max_i              %计算点j与其最近邻点在i个离散步后的距离
    
78. 
    
79.         d_j_i=0;
    
80. 
    
81.         for kk=1:m
    
82. 
    
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));
    
84. 
    
85.             d(k,j)=d_j_i;
    
86. 
    
87.         end
    
88. 
    
89.     end
    
90. 
    
91. end
    
92. 
    
93. 
    
94. 
    
95. %对每个演化时间步长i，求所有的j的lnd(i,j)平均
    
96. 
    
97. [l_i,l_j]=size(d);
    
98. 
    
99. for i=1:l_i
    
100. 
     
101.     q=0;
     
102. 
     
103.     y_s=0;
     
104. 
     
105.     for j=1:l_j
     
106. 
     
107.         if d(i,j)~=0
     
108. 
     
109.             q=q+1;
     
110. 
     
111.             y_s=y_s+log(d(i,j));
     
112. 
     
113.         end
     
114. 
     
115.     end
     
116. 
     
117.     y(i)=y_s/(q*delt_t);
     
118. 
     
119. end
     
120. 
     
121. x=1:length(y);
     
122. 
     
123. y
     
124. 
     
125. pp=polyfit(x,y,1);
     
126. 
     
127. lambda_1=pp(1);
     
128. 
     
129. yp=polyval(pp,x);
     
130. 
     
131. plot(x,y,'-o',x,yp,'--')
     
132. 
     
133. function X=reconstitution(data,N,m,tau)
     
134. %该函数用来重构相空间
     
135. % m为嵌入空间维数
     
136. % tau为时间延迟
     
137. % data为输入时间序列
     
138. % N为时间序列长度
     
139. % X为输出,是m*n维矩阵
     
140. 
     
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
     
147. 
     

复制代码

chunshui2003

1. 
   
2. 
   
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. %*************************************************************
    
18. 
    
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
    
47. 
    
48. end
    
49. 
    
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--');
    
82. 
    
83. 
    
84. function X=reconstitution(data,N,m,tau)
    
85. %该函数用来重构相空间
    
86. % m为嵌入空间维数
    
87. % tau为时间延迟
    
88. % data为输入时间序列
    
89. % N为时间序列长度
    
90. % X为输出,是m*n维矩阵
    
91. 
    
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
    
98. 
    

复制代码

c.        最小二乘法数据拟合

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

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

hsfy919

回复 4 # chunshui2003 的帖子

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

chunshui2003

回复 5 # hsfy919 的帖子

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

hsfy919

回复 6 # chunshui2003 的帖子

可能咱俩的研究有相似之处，我做的也含有非线性油膜力，不过我是12个自由度，但效果也挺好，呵呵，所以你先试试吧

chunshui2003

回复 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

oil force

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

因为你提到了也采用油膜力，而这个油膜力应用也较为广泛，所以想问一下你是如何处理的。或者在求解Jacobian矩阵的时候有没有其他更好的方法。

octopussheng

sign就是符号函数，可用if语句来实现。

hsfy919

回复 8 # chunshui2003 的帖子

oct说的对，你试试用if语句替换一下符号函数，改起来应该也不难。另外我不是用的这个油膜力表达式，是给予不同假设自己推导的表达式。

chunshui2003

回复 9 # octopussheng 的帖子

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

chunshui2003

回复 10 # hsfy919 的帖子

谢谢你的回答，我今天再试一下，尝试不同的方法。

octopussheng

回复 11 # chunshui2003 的帖子

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

只要定义微分方程的程序能正确求解，问题就解决了。

chunshui2003

回复 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

回复 13 # octopussheng 的帖子


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

1. 
   
2. 
   
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;            %大轴旋转角速度
    
30. 
    
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);
    
50. 
    
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.   ];
    
78. 
    
79. 
    
80. uu = zeros(90,1);
    
81. 
    
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;
    
91. 
    
92. 
    
93. J = [];
    
94. 
    
95. 
    
96. uu(10:90) = J*Q;
    
97. 
    

复制代码

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