夷红志 发表于 2014-12-1 15:27

请求无阻尼三自由度振动系统的的matlab实现程序对不对?


         无阻尼三自由度振动系统。已知系统质量与刚度参数m1、m2、m3、k1、k2、k3如表1所示,初始条件取为x10=0.5m,x20=x30=0m,0m/s,0m/s。假定质量m2作用简谐激振力f=5sin0.2t,求振动总响应。function shiyan02(m1,m2,m3,k1,k2,k3x10,x20,x30,f1,f2,f3,w,x100,x200,x300)m1=1;m2=2;m3=3;k1=1;k2=2;k3=3;x10=0.5;x20=0;x30=0x100=0,x200=0;x300=0M=K=A=inv(K) =eig(M^(-1)*K) u11(:,1)=u(:,3)u11(:,2)=u(:,2)u11(:,3)=u(:,1)U(:,1)=u11(:,1)/u11(1,1)U(:,2)=u11(:,2)/u11(1,2)U(:,3)=u11(:,3)/u11(1,3)data=wn2^0.5wn1=data(3,3)wn2=data(2,2)wn3=data(1,1)wn=U1=inv(U);y0=U1*'y00=U1*'M1=U'*M*UK1=U'*K*Uj=1:3;uz(:,j)=U(:,j)/(M1(j,j))^0.5f1=0;f2=5;f3=0;f=Fz=uz'*fw=0.2b1=1/abs(1-w^2/wn1^2)b2=1/abs(1-w^2/wn2^2)b3=1/abs(1-w^2/wn3^2)n1=Fz(1)/K1(1,1)*b1n2=Fz(2)/K1(2,2)*b1n3=Fz(3)/K1(3,3)*b1t=0:0.001:200;Y1=y0(1)*cos(wn1*t)+((y00(1)/wn1)-n1*(w/wn1))*sin(wn1*t)+n1*sin(w*t);Y2=y0(2)*cos(wn2*t)+((y00(2)/wn2)-n2*(w/wn2))*sin(wn2*t)+n2*sin(w*t);Y3=y0(3)*cos(wn3*t)+((y00(3)/wn3)-n3*(w/wn3))*sin(wn3*t)+n3*sin(w*t);figure(1)plot(t,Y1,'k:',t,Y2,'b',t,Y3,'r--');xlabel('t(s)')set(get(gca,'xlabel'),'fontsize',20)ylabel('Y(m)')set(get(gca,'ylabel'),'fontsize',20)set(gca,'fontsize',20)title('金宇','fontsize',35)legend('Y1虚线','Y2实线','Y3点划线')grid onhold on y1=n1*sin(w*t);y2=n2*sin(w*t);y3=n3*sin(w*t);X1=Y1*uz(1,1)+Y2*uz(1,2)+Y3*uz(1,3);X2=Y1*uz(2,1)+Y2*uz(2,2)+Y3*uz(2,3);X3=Y1*uz(3,1)+Y2*uz(3,2)+Y3*uz(3,3);figure(2)plot(t,X1,'k:',t,X2,'b',t,X3,'r--');xlabel('t(s)')set(get(gca,'xlabel'),'fontsize',20)ylabel('X(m)')set(get(gca,'ylabel'),'fontsize',20)set(gca,'fontsize',20)title('jinyu','fontsize',35)legend('X1虚线','X2实线','X3点划线')file:///C:/Users/ADMINI~1/AppData/Local/Temp/msohtml1/01/clip_image002.gifgrid on
hold on

页: [1]
查看完整版本: 请求无阻尼三自由度振动系统的的matlab实现程序对不对?