求解微分方程时提示“index out of bounds because numel(u)=4”
请大家看一下我建立的运动微分方程和求解过程。运动微分方程如下：
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=;
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);
=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
=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函数和计算程序是哪里出了问题才导致出现错误提示的。 应该是你的参数传递不对吧，diff_equ_elastic函数调用时，y0只有4个元素呀。
回复 沙发 messenger 的帖子
非常感谢messenger，已经不是第一次帮助我了。这是我自己的大意造成的，因为之前曾经算过一个两自由度4个变量的方程，结果就把程序复制过来忘记改变初始值设置了。这样的错误以后不会再犯了！真的非常感谢！
页:
[1]