求解微分方程时提示“index out of bounds because numel(u)=4”

chunshui2003

请大家看一下我建立的运动微分方程和求解过程。

运动微分方程如下：

1. 
   
2. 
   
3. function uu=diff_equ_elastic(t,u)
   
4. uu = zeros(8,1);
   
5. %%%%% 基本参数数值 %%%%%
   
6. mr = 600*10^3;            
   
7. mb = 1000;                 
   
8. De = 0.5*10^6;            
   
9. Ke = 0.5*10^10;            
   
10. delta_0 = 18*10^-3;        
    
11. e_0 = 2*10^-3;            
    
12. omega = 13.1;              
    
13. cb = 0.25*10^-3;           
    
14. 
    
15. %%%%% Fx,Fy 表达式 %%%%%%
    
16. Fx = mr*omega^2*e_0*cos(omega*t);
    
17. Fy = mr*omega^2*e_0*sin(omega*t);
    
18. global K11 K12 K21 K22 K_1 K_2
    
19. %%%% 关于轴颈偏心率和极角的参数设置 %%%
    
20. epsilon = sqrt(u(5)^2+u(7)^2);       %错误就出现在这一行，提示说数据溢出了！   
    
21. psi = atan(u(5)./(-u(7)));            
    
22. epsilon_de = (u(5).*u(6)+u(7).*u(8))./epsilon;        
    
23. psi_de = (u(5).*u(8)-u(7).*u(6))./epsilon.^2;         
    
24. G1_epsi=2.*epsilon.^2./(1-epsilon.^2).^2;
    
25. G2_epsi=pi*(1+2.*epsilon.^2)./(2*(1-epsilon.^2).^(2/5));
    
26. G3_epsi=pi*epsilon./(2*(1-epsilon).^(3/2));
    
27. G4_epsi=2*epsilon./(1-epsilon.^2).^2;
    
28. 
    
29. S_0 = 0.001;                        
    
30. L_vs_Rb = 1;                        
    
31. Coeff_oil = pi*S_0/(2*mb*omega*cb)*L_vs_Rb.^2;                %合成的一个系数，在方程中看着方便
    
32. 
    
33. P_epsi = (1-2.*psi_de).*G1_epsi + 2.*epsilon_de.*G2_epsi;
    
34. P_psi = (1-2.*psi_de).*G3_epsi + 2.*epsilon_de.*G4_epsi;
    
35. %%% 系统运动微分方程 %%%
    
36. 
    
37. uu(1) = u(2);
    
38. uu(2) = -De/(mr*omega).*u(2) -(K11+Ke)/(mr*omega^2).*u(1) -K12/(mr*omega^2).*u(3) +Ke*cb/(mr*omega^2*delta_0).*u(5) +e_0.*cos(omega*t)/delta_0 -K_1/(mr*omega^2*delta_0);
    
39. uu(3) = u(4);
    
40. uu(4) = -De/(mr*omega).*u(4) -K21/(mr*omega^2).*u(1) -(K22+Ke)/(mr*omega^2).*u(3) +Ke*cb/(mr*omega^2*delta_0).*u(7) +e_0.*sin(omega*t)/delta_0 -K_2/(mr*omega^2*delta_0);
    
41. uu(5) = u(6);
    
42. uu(6) = Ke*delta_0/(2*mb*omega^2*cb).*u(1) -Ke/(2*mb*omega^2).*u(5) -Coeff_oil.*P_epsi.*u(5)./epsilon - Coeff_oil.*P_psi.*u(7)./epsilon;
    
43. uu(7) = u(8);
    
44. uu(8) = Ke*delta_0/(2*mb*omega^2*cb).*u(3) -Ke/(2*mb*omega^2).*u(7) -Coeff_oil.*P_psi.*u(5)./epsilon + Coeff_oil.*P_epsi.*u(7)./epsilon;
    
45. 
    

复制代码

然后在主窗口中运行：

1. 
   
2. 
   
3. tic                                                                                                   
   
4. clear all                                                                                             
   
5. clc                                                                                                   
   
6. load K.mat K                    %从K.mat文件中将参数K提取出来                                         
   
7. load K.mat tt                                             
   
8. global K11 K12 K21 K22 K_1 K_2                                                                        
   
9. y0=[0.001 0.001 0.001 0.001];                                                                        
   
10. period=2*pi/13.1;               %%周期                                                               
    
11. step=period/100;                %%步长                                                               
    
12. YY=zeros(100,8);                %提取前100个点的结果            
    
13. for i=1:100                                                                                          
    
14.     i                                                                                                
    
15.     K11=K(1,i);                                                                                       
    
16.     K12=K(2,i);                                                                                       
    
17.     K21=K(3,i);                                                                                       
    
18.     K22=K(4,i);                                                                                       
    
19.     K_1=K(5,i);                                                                                       
    
20.     K_2=K(6,i);                                                                                       
    
21.     [t,Y]=ode45(@diff_equ_elastic,((i-1)*step:step:i*step),y0);           
    
22.     YY(i,:)=Y(end,:);                                   
    
23.     y0=Y(end,:);                                %将每一个步长计算后得到的Y值最终值作为下一次计算的初值
    
24. end                                                                                                   
    
25. toc     
    
26. 
    

复制代码

结果提示：

Attempted to access u(5); index out of bounds because numel(u)=4.

Error in ==> diff_equ_elastic at 26
epsilon = sqrt(u(5)^2+u(7)^2);      
   
Error in ==> odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> solu_equ_elas at 20
    [t,Y]=ode45(@diff_equ_elastic,((i-1)*step:step:i*step),y0);


其中，参数K是在不同时刻下的电磁刚度数值，每一个时间子步的K值都是不同的，具体数值我没有列出来。K的规模为K=zeros(6,101)


上面提示说访问u(5)的时候就溢出了，但我觉得不可能。之前我计算过其他的运动微分方程，也是将其中的某个参数用u(i)的形式表示，仍然可以计算，并且同样是u(1)到u(8)这8个变量。
请问大家我的function函数和计算程序是哪里出了问题才导致出现错误提示的。

messenger

应该是你的参数传递不对吧，diff_equ_elastic函数调用时，y0只有4个元素呀。

chunshui2003

非常感谢messenger，已经不是第一次帮助我了。这是我自己的大意造成的，因为之前曾经算过一个两自由度4个变量的方程，结果就把程序复制过来忘记改变初始值设置了。这样的错误以后不会再犯了！
真的非常感谢！