请求无阻尼三自由度振动系统的的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=0

x100=0,x200=0;x300=0

M=[m1,0,0;0,m2,0;0,0,m3]

K=[k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3]  

A=inv(K)

[u,wn2]=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.5

wn1=data(3,3)

wn2=data(2,2)

wn3=data(1,1)

wn=[wn1 0 0;0 wn2 0;0 0 wn3 ]

U1=inv(U);

y0=U1*[x10,x20,x30]'

y00=U1*[x100,x200,x300]'

M1=U'*M*U

K1=U'*K*U  

j=1:3;

uz(:,j)=U(:,j)/(M1(j,j))^0.5

f1=0;f2=5;f3=0;

f=[f1;f2;f3]

Fz=uz'*f

w=0.2

b1=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)*b1

n2=Fz(2)/K1(2,2)*b1

n3=Fz(3)/K1(3,3)*b1

t=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 on

hold 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