请教一个动力学 系统方程解法的问题

wyy20071127

动力学方程为:

定义的function程序为
function ydot=ZCX_sub(t,y)
n=300;w=2*pi*n/60;x=w*t;
m1=13300;m2=4300;m3=11583;
c1=7.87E+5;c2=0;c3=1.117E+5;
k1=2.91E+10;k2=2.37E+10;k3=6.737E+6;
M=[m1,0,0;0,m2,0;0,0,m3];
C=[c1+c2,-c2,-c1;-c2,c2,0;-c1,0,c1+c3];
K=[k1+k2,-k2,-k1;-k2,k2,0;-k1,0,k1+k3];
A=[zeros(3,3),eye(3,3);-inv(M)*K,-inv(M)*C];
F=[1;0;0]*500*cos(x)+[0;1;0]*300*cos(x);
P=[zeros(3,1);inv(M)*F];
ydot=A*y+P;
分别使用ode45和RK方法编程得到的结果都不对
ode45程序
clc
clear
t0=0;
y0=zeros(6,1);
h=0.02;
[t,y]=ode45('ZCX_sub',0,0.8,y0,1e-10);
plot(t,y(:,1));
RK方法程序
clear
clc
t0=0;h=0.02;
j=1;
t=0:h:0.8;n=length(t);
for i=1:n;
   y0=zeros(6,1);
   l1=ZCX_sub(t0,y0);
   l2=ZCX_sub(t0+h/2,y0+(h/2)*l1);
   l3=ZCX_sub(t0+h/2,y0+(h/2)*l2);
   l4=ZCX_sub(t0+h,y0+h*l3);     
   t1=t0+h;
   y1=y0+(h/6)*(l1+2*l2+2*l3+l4);
   u1(j)=y1(1);u2(j)=y1(2);u3(j)=y1(3);
   t0=t1;y0=y1;
   j=j+1;
end
plot(t,u1);
请问哪里不对啊...

dw04116

请教下，怎么往上贴图片啊。。你的东西太多了，实在看不下去，貌似是振动的数值算法，我前面也用过纽马克直接积分法，，也有很多问题没解决。。

wyy20071127

贴图片要用高级模式

笑石头

拉普拉斯变换后使用SIMULINK建模求解还方便一些

wyy20071127

使用simulink仿真修改参数不方便啊