对无阻尼两自由度自由振动的振动系统

对无阻尼两自由度自由振动的振动系统，质量块1和质量块2有初始位移x1=2，x2=2，初速度x3=0.8,x4=1.3。弹簧刚度k1=9,k2=12,k3=9。质量均为3kg。求位移与时间之间的关系。

syms k1k2k3m1m2abcdX1X2C1C2wehpsi1psi2r1r2tx1x2x3x4;

X1=C1*cos(w_1*t-psi1)+C2*cos(w_2*t-psi2);

X2=C1*r1*cos(w_1*t-psi1)+C2*r2*cos(w_2*t-psi2);

x1=2;

x2=2;

x3=0.8;

x4=1.3;

k=[9,12,9];

m=[3,3];

a=(k(1)+k(2))/m(1);

b=k(2)/m(1);

c=k(2)/m(2);

d=(k(2)+k(3))/m(2);

y1=w^2-(a+d)*w+(a*d-b*c);

y=solve('w^2 - 14*w + 33=0',w);

e=y(1);

h=y(2);

w=[e,h];

A=[(a-e^2)/b,(a-h^2)/b];

r1=simplify(A(1));

r2=simplify(A(2));

C1=(abs(r2-r1))^(-1)*sqrt((r2*x1-x2)^2+(r2*x3-x4)^2/e^2);

C2=(abs(r2-r1))^(-1)*sqrt((x2-r1*x1)+(x4-r1*x3)^2/h^2);

psi1=atan((r2*x3-x4)/(e*(r2*x1-x2)));

psi2=atan((r1*x3-x4)/(h*(r1*x1-x2)));

ts=0:0.01:10;

X1=C1*cos(e*ts-psi1)+C2*cos(h*ts-psi2);

X2=C1*r1*cos(e*ts-psi1)+C2*r2*cos(h*ts-psi2);

plot(ts,X1,'b',ts,X2,'r')