|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
请大家看一下我建立的运动微分方程和求解过程。
运动微分方程如下:
-
- function uu=diff_equ_elastic(t,u)
- uu = zeros(8,1);
- %%%%% 基本参数数值 %%%%%
- mr = 600*10^3;
- mb = 1000;
- De = 0.5*10^6;
- Ke = 0.5*10^10;
- delta_0 = 18*10^-3;
- e_0 = 2*10^-3;
- omega = 13.1;
- cb = 0.25*10^-3;
- %%%%% Fx,Fy 表达式 %%%%%%
- Fx = mr*omega^2*e_0*cos(omega*t);
- Fy = mr*omega^2*e_0*sin(omega*t);
- global K11 K12 K21 K22 K_1 K_2
- %%%% 关于轴颈偏心率和极角的参数设置 %%%
- epsilon = sqrt(u(5)^2+u(7)^2); %错误就出现在这一行,提示说数据溢出了!
- psi = atan(u(5)./(-u(7)));
- epsilon_de = (u(5).*u(6)+u(7).*u(8))./epsilon;
- psi_de = (u(5).*u(8)-u(7).*u(6))./epsilon.^2;
- G1_epsi=2.*epsilon.^2./(1-epsilon.^2).^2;
- G2_epsi=pi*(1+2.*epsilon.^2)./(2*(1-epsilon.^2).^(2/5));
- G3_epsi=pi*epsilon./(2*(1-epsilon).^(3/2));
- G4_epsi=2*epsilon./(1-epsilon.^2).^2;
- S_0 = 0.001;
- L_vs_Rb = 1;
- Coeff_oil = pi*S_0/(2*mb*omega*cb)*L_vs_Rb.^2; %合成的一个系数,在方程中看着方便
- P_epsi = (1-2.*psi_de).*G1_epsi + 2.*epsilon_de.*G2_epsi;
- P_psi = (1-2.*psi_de).*G3_epsi + 2.*epsilon_de.*G4_epsi;
- %%% 系统运动微分方程 %%%
- uu(1) = u(2);
- 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);
- uu(3) = u(4);
- 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);
- uu(5) = u(6);
- 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;
- uu(7) = u(8);
- 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;
-
复制代码
然后在主窗口中运行:
-
- tic
- clear all
- clc
- load K.mat K %从K.mat文件中将参数K提取出来
- load K.mat tt
- global K11 K12 K21 K22 K_1 K_2
- y0=[0.001 0.001 0.001 0.001];
- period=2*pi/13.1; %%周期
- step=period/100; %%步长
- YY=zeros(100,8); %提取前100个点的结果
- for i=1:100
- i
- K11=K(1,i);
- K12=K(2,i);
- K21=K(3,i);
- K22=K(4,i);
- K_1=K(5,i);
- K_2=K(6,i);
- [t,Y]=ode45(@diff_equ_elastic,((i-1)*step:step:i*step),y0);
- YY(i,:)=Y(end,:);
- y0=Y(end,:); %将每一个步长计算后得到的Y值最终值作为下一次计算的初值
- end
- 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函数和计算程序是哪里出了问题才导致出现错误提示的。 |
-
|