clc % Graphics parameters xbox1=[20 30 30 20 20];ybox1=[0 0 10 10 0 ]; xbox2=[50 60 60 50 50];ybox2=ybox1; xspring1=[0 2 4 8 12 16 18 20];yspring=[5 5 7 3 7 3 5 5]; xspring2=xspring1+30*ones(1,8); xspring3=xspring1+60*ones(1,8); basex=[0 -5 1 18 28 52 75 84 80 80 0 0 ]; basey=[15 3 -8 -4 -10 -3 -7 -2 15 -.3 -.3 15]; arrowx=[30 44 41 41 44];arrowy=[15 15 16 14 15]; % Draw system figure(1);clf axis off axis([-10 90 -20 30]);hold on text(15,25,'Two Degree of Freedom System'); patch(basex,basey,'r'); patch(xbox1,ybox1,'r'); patch(xbox2,ybox2,'r'); plot(xspring1,yspring,'k','LineWidth',[2]); plot(xspring2,yspring,'k','LineWidth',[2]); plot(xspring3,yspring,'k','LineWidth',[2]); plot(arrowx,arrowy,'k');plot(arrowx+30*ones(1,5),arrowy,'k'); text(8,8,'k_1');text(38,8,'k_2');text(68,8,'k_3'); text(23,8,'m_1');text(53,8,'m_2'); text(32,13,'u_1');text(62,13,'u_2') text(15,-12,'Press Enter to Continue') pause ndof=2; k1=2;k2=2;k3=2; m1=1;m2=1; a=1; b=-(((k1+k2)/m1)+((k2+k3)/m2)); c=((k1*k2)+(k2*k3)+(k3*k1))/(m1*m2); omega1sq=(-b-sqrt((b^2)-4*a*c))/(2*a); omega2sq=(-b+sqrt((b^2)-4*a*c))/(2*a); disp('The radian natural frequencies are') omega1=sqrt(omega1sq) omega2=sqrt(omega2sq) mu1=(k1+k2-m1*omega1sq)/k2; mu2=(k1+k2-m1*omega2sq)/k2; mode1=[1;mu1] mode2=[1;mu2] omegan=[omega1 omega2]; disp(' The Hertzian undamped natural frequencies are') fn=omegan/(2*pi) disp(' The undamped natural periods are') Tn=ones(1,ndof)./fn% natural periods in sec. disp(' The modal matrix is') PHI=[mode1 mode2] u1zero=[];u2zero=[]; disp(' Input initial displacements') while isempty(u1zero) u1zero=input(' u1(0) = '); end while isempty(u2zero) u2zero=input(' u2(0) = '); end uzero=[u1zero;u2zero]; disp(' Press Enter to Continue') pause A=inv(PHI)*uzero; Ts=2*pi/(40*omegan(1,2)); t=0:Ts:400*Ts; c=[A(1,1)*cos(omegan(1,1)*t);A(2,1)*cos(omegan(1,2)*t)]; u=PHI*c; %Plot displacements figure(2);clf Umax=max(max(abs(u))); axis([0 400*Ts -1.2*Umax 1.2*Umax]); plot(t,u) legend('u_1(t)','u_2(t)'); xlabel('Time t, sec.') ylabel('Displacements, u_1(t) and u_2(t)'); text(t(1,201),Umax/2,'Press Enter to Continue'); pause % Set up animation figure(1);clf axis off axis([-10 90 -20 30]);hold on; text(20,25,'Two Degree of Freedom System'); patch(basex,basey,'r');hold on Mass1=patch(xbox1,ybox1,'r','EraseMode','xor','LineWidth',[2]);hold on Mass2=patch(xbox2,ybox2,'r','EraseMode','xor','LineWidth',[2]);hold on Spring1=plot(xspring1,yspring,'k','EraseMode','xor','LineWidth',[2]);hold on Spring2=plot(xspring2,yspring,'k','EraseMode','xor','LineWidth',[2]);hold on Spring3=plot(xspring3,yspring,'k','EraseMode','xor','LineWidth',[2]);hold on txt1='Press Enter to Set Initial Displacement'; Txt1=text(15,-12,txt1); pause % Set initial displacement configuration U=u*7/Umax; figure(1);clf axis off;hold on axis([-10 90 -20 30]);hold on; text(15,25,'Two Degree of Freedom System'); patch(basex,basey,'r');hold on Mass1=patch(xbox1+U(1,1)*ones(1,5),ybox1,'r','EraseMode','xor','LineWidth',[2]);hold on Mass2=patch(xbox2+U(2,1)*ones(1,5),ybox2,'r','EraseMode','xor','LineWidth',[2]);hold on xspring11=xspring1+[0 0 U(1,1)*.125 U(1,1)*.375 U(1,1)*.625 U(1,1)*.875 U(1,1) U(1,1)]; Spring1=plot(xspring11,yspring,'k','EraseMode','xor','LineWidth',[2]);hold on xrel=U(2,1)-U(1,1); xspring22=xspring2+[U(1,1) U(1,1) U(1,1)+.125*xrel U(1,1)+.375*xrel U(1,1)+.625*xrel U(1,1)+.875*xrel U(2,1) U(2,1)]; Spring2=plot(xspring22,yspring,'k','EraseMode','xor','LineWidth',[2]);hold on xspring33=xspring3+[U(2,1) U(2,1) .875*U(2,1) .625*U(2,1) .375*U(2,1) .125*U(2,1) 0 0]; Spring3=plot(xspring33,yspring,'k','EraseMode','xor','LineWidth',[2]);hold on txt1=['Press Enter to Animate umax =' num2str(Umax)]; Txt1=text(15,-12,txt1); Txt2=text(15,-16,['t = ' num2str(t(1,1)) ' sec.']); pause % Animate system for k=1:400 set(Mass1,'Xdata',xbox1+U(1,k)*ones(1,5)); set(Mass2,'Xdata',xbox2+U(2,k)*ones(1,5)); xspring11=xspring1+[0 0 .125*U(1,k) .375*U(1,k) .625*U(1,k) .875*U(1,k) U(1,k) U(1,k)]; set(Spring1,'Xdata',xspring11) xrel=U(2,k)-U(1,k); xspring22=xspring2+[U(1,k) U(1,k) U(1,k)+.125*xrel U(1,k)+.375*xrel U(1,k)+.625*xrel U(1,k)+.875*xrel U(2,k) U(2,k)]; set(Spring2,'Xdata',xspring22); xspring33=xspring3+[U(2,k) U(2,k) .875*U(2,k) .625*U(2,k) .375*U(2,k) .125*U(2,k) 0 0]; set(Spring3,'Xdata',xspring33) set(Txt2,'String',['t = ' num2str(t(1,k)) ' sec.']); pause(.04) end