贴出来，让大家看看 本帖最后由 tifanhuang 于 2012-8-3 20:08 编辑
求解的微分方程组（在附件中）
其中未知量为 , , .以下是我用的四阶龙格—库塔算法求解的程序：第一部分是建立函数：function dy = zunifun(t,y )
%此为两刚性支撑下悬臂转子模型
%以下为求解转轴刚度
E0=2.1e5;%转轴弹性模量
a=5.5e-1;
l=1.0e2;
d=1.0e1;%转轴直径
j=1/(64*pi*d^4);
x11=-(a^2)*(l+a)/(3*E0*j);
x12=a*(2*l+3*a)/(6*E0*j);
x21=x12;
x22=(l-3*a)/(3*E0*j);
A=;
K=inv(A);
k11=K(1,1);
k12=K(1,2);
k21=K(2,1);
k22=K(2,2);
Jp=2.7596e-11;%转子极转动惯量
Jd=0.5*Jp;%转子直径转动惯量
Cr=0.01;%转子阻尼系数
Ct=0.01;%转子阻尼系数
Cb=0.3;%外阻尼系数
m=8.64e-2;%转子质量
eu=5;%偏心距
f=20;%转子频率
Q=2*pi*f*t;
g=9.98;% 重力加速度
%转子运动微分方程组
dy=zeros(8,1)%行向量
dy(1)=y(2);
dy(2)=eu*(2*pi*f)^2*cos(Q)-((Ct+Cb)*y(2)+k11*y(1)+k12*y(6)+m*g)/m;
dy(3)=y(4);
dy(4)=eu*(2*pi*f)^2*sin(Q)-((Ct+Cb)*y(4)+k11*y(3)-k12*y(7))/m;
dy(5)=y(6);
dy(6)=(Jp*2*pi*f*y(8)-Cr*y(6)-k21*y(1)-k22*y(5))/Jd;
dy(7)=y(8);
dy(8)=-(Jp*2*pi*f*y(6)+Cr*y(8)-k21*y(3)+k22*y(7))/Jd;end第二部分是求解过程：clc;
t0=0;
tN=2;
y0=;
h=0.001;
t=t0:h:tN;
N=length(t);
j=1;
for i=1:N
t1=t0+h;
K1=zunifun(t0,y0);
K2=zunifun(t0+h/2,y0+h*K1/2);
K3=zunifun(t0+h/2,y0+h*K2/2);
K4=zunifun(t0+h,y0+h*K3);
y1=y0+(h/6)*(K1+2*K2+2*K3+K4);
yy1(j)=y1(1);
yy2(j)=y1(2);
yy3(j)=y1(3);
yy4(j)=y1(4);
yy5(j)=y1(5);
yy6(j)=y1(6);
yy7(j)=y1(7);
yy8(j)=y1(8);
t0=t1;
y0=y1;
j=j+1;
end出来的结果全是0，希望各高手能帮忙指出其中的错误，谢谢。 回复 17 # tifanhuang 的帖子
你为什么不用ode45来解呢、控制一下时间步长，也算是定步长啊。这样显得有点麻烦 我也用了，可一直出现电脑在跳动
如果ode45可以求解，那ode45跟四阶龙格—库塔算法也应该是差不多的，只是精度问题，而我的现在求的结果全为0
以上的程序有哪位大神帮帮忙呗，是在纠结呀 我看了两遍，怎么没有找到你的负载和驱动力矩啊
页:
1
[2]