声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3155|回复: 2

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

[复制链接]
发表于 2010-7-25 08:47 | 显示全部楼层 |阅读模式

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

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

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

运动微分方程如下:



  1. function uu=diff_equ_elastic(t,u)
  2. uu = zeros(8,1);
  3. %%%%% 基本参数数值 %%%%%
  4. mr = 600*10^3;            
  5. mb = 1000;                 
  6. De = 0.5*10^6;            
  7. Ke = 0.5*10^10;            
  8. delta_0 = 18*10^-3;        
  9. e_0 = 2*10^-3;            
  10. omega = 13.1;              
  11. cb = 0.25*10^-3;           

  12. %%%%% Fx,Fy 表达式 %%%%%%
  13. Fx = mr*omega^2*e_0*cos(omega*t);
  14. Fy = mr*omega^2*e_0*sin(omega*t);
  15. global K11 K12 K21 K22 K_1 K_2
  16. %%%% 关于轴颈偏心率和极角的参数设置 %%%
  17. epsilon = sqrt(u(5)^2+u(7)^2);       %错误就出现在这一行,提示说数据溢出了!   
  18. psi = atan(u(5)./(-u(7)));            
  19. epsilon_de = (u(5).*u(6)+u(7).*u(8))./epsilon;        
  20. psi_de = (u(5).*u(8)-u(7).*u(6))./epsilon.^2;         
  21. G1_epsi=2.*epsilon.^2./(1-epsilon.^2).^2;
  22. G2_epsi=pi*(1+2.*epsilon.^2)./(2*(1-epsilon.^2).^(2/5));
  23. G3_epsi=pi*epsilon./(2*(1-epsilon).^(3/2));
  24. G4_epsi=2*epsilon./(1-epsilon.^2).^2;

  25. S_0 = 0.001;                        
  26. L_vs_Rb = 1;                        
  27. Coeff_oil = pi*S_0/(2*mb*omega*cb)*L_vs_Rb.^2;                %合成的一个系数,在方程中看着方便

  28. P_epsi = (1-2.*psi_de).*G1_epsi + 2.*epsilon_de.*G2_epsi;
  29. P_psi = (1-2.*psi_de).*G3_epsi + 2.*epsilon_de.*G4_epsi;
  30. %%% 系统运动微分方程 %%%

  31. uu(1) = u(2);
  32. 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);
  33. uu(3) = u(4);
  34. 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);
  35. uu(5) = u(6);
  36. 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;
  37. uu(7) = u(8);
  38. 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;

复制代码


然后在主窗口中运行:



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

复制代码


结果提示:

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函数和计算程序是哪里出了问题才导致出现错误提示的。
弹性支撑运动微分方程.jpg
回复
分享到:

使用道具 举报

发表于 2010-7-25 18:32 | 显示全部楼层
应该是你的参数传递不对吧,diff_equ_elastic函数调用时,y0只有4个元素呀。

评分

1

查看全部评分

 楼主| 发表于 2010-7-26 08:52 | 显示全部楼层

回复 沙发 messenger 的帖子

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

本版积分规则

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

GMT+8, 2024-11-14 15:56 , Processed in 0.077560 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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