clear
clc
t0=0;Y0=zeros(6,1);tN=20;h=0.01;%输入初始条件t0,y0,步长h，外加激励P0,周期w;
t=t0:h:tN;N=length(t);P0=1;w=3;i=1;
for i=1:N %主程序
t1=t0+h;
P=[P0*sin(w*t0);0.0;0.0];
K1=ZXC(t0,Y0,P);
P=[P0*sin(w*(t0+h/2));0.0;0.0];
K2=ZXC(t0+h/2,Y0+h*K1/2,P);
K3=ZXC(t0+h/2,Y0+h*K2/2,P);
P=[P0*sin(w*(t0+h));0.0;0.0];
K4=ZXC(t0+h,Y0+h*K3,P);
Y1=Y0+(h/6)*(K1+2*K2+2*K3+K4);
YY1(i)=Y1(1);YY2(i)=Y1(2);YY3(i)=Y1(3);
t0=t1;Y0=Y1;i=i+1;
end
%画出不同质量块的运动图像
subplot(131),plot(t,YY1);xlabel('{\itt}');ylabel('{\ity1}');title('R__K质量块1运动图');grid
subplot(132),plot(t,YY2);xlabel('{\itt}');ylabel('{\ity2}');title('R__K质量块2运动图');grid
subplot(133),plot(t,YY3);xlabel('{\itt}');ylabel('{\ity3}');title('R__K质量块3运动图');grid
function ydot=ZXC(t,Y,P) %子函数程序
M=[1,0,0;0,1,0;0,0,1];
K=[4,-2,0;-2,3,-1;0,-1,3];
C=zeros(size(M));
P1=[zeros(3,1);inv(M)*P];
A=[zeros(3,3),eye(3,3);-inv(M)*K,-inv(M)*C];
ydot=A*Y+P1;
end
