动力系统程序及图像.rar_husbandn2k_动力系统_龙格——库塔法进行动力系统计算_龙格库塔_龙格库塔 步长

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