|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 wyy20071127 于 2012-3-14 21:24 编辑
动力学方程为:
定义的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);
请问哪里不对啊...
|
|