% forcedstring.m Last modified 8/1/2008 % String driven by a sinusoidal point force of amplitude P0 at x/L=a/L % with radian frequency omega. Inputs are the ratio of the % forcing frequency to the first natural frequency and the nondimensional % forcing location, a/L. % This m-file was written at the University of Wyoming in the Electrical % and Computer Engineering Department and is to be distributed without % cost. clear all set(0,'DefaultAxesFontSize',12); set(0,'DefaultTextFontSize',12); fratio=-1; while fratio<=0|fratio==1|fratio==2|fratio==3 string1=' The ratio of excitation freq. to first natural freq. = '; fratio=input(string1); if fratio==1|fratio==2|fratio==3 disp(' ') disp(' This is the resonant case choose another freq. ratio'); disp(' ') end end disp(' ') aoverL=-1; while aoverL<=0|aoverL>=1 aoverL=input('The location of the force is a/L = '); end d=2*fix((aoverL*100)+.001)+1; xoverL=linspace(0,1,201); x1=xoverL(1,1:d);x2=xoverL(1,d+1:201); betaL=pi*fratio; delta=sin(betaL); y1=(1/betaL)*(sin(betaL*(1-aoverL))*sin(betaL*x1))/delta; y2=(1/betaL)*(sin(betaL*(ones(1,201-(d))-x2))*sin(betaL*aoverL))/delta; y=[y1 y2]; lim=max(abs(y)); figure(1);clf; axis([0 1 -.2*lim 1.2*lim]); hold on plot(xoverL,abs(y)) hold on plot([0 1],[0 0]); text(.8,1.1*lim,['f/f_1 = ',num2str(fratio)]) xlabel('Dimensionless distance from the fixed end, x/L') ylabel('Dimensionless Amplitude, y(x)T/P_0L') hold off text(.1,1.1*lim,'Press Enter to Continue') xp1=.5+.5*[1 1 1.1 1.3 1.2 1.32 1.2 1 1]; yp1=1.5*lim*[0 .4 .3 .2 .1 0 -.2 -.4 0]; pause figure(2);clf; axis([-.3 1.3 -2*lim 2*lim]); hold on set(gca,'Box','on') patch(xp1,yp1,'r') xp2=[0 0 -.1 -.3 -.2 -.32 -.2 0 0]*.5; patch(xp2,yp1,'r'); xlabel('Distance, x/L') ylabel('Displacement, y(x,t)T/P_0L') hold on x1=[0 1];y1=[0 0]; plot(x1,y1) hold on texthandle2=text(-.1,-1.6*lim,['f/f_1 = ',num2str(fratio)]); ys=y*0; plot(xoverL,ys,'k','LineWidth',[1.5]); hold on ydata=[ys]; texthandle2=text(0,1.5*lim,'Press Enter to Animate one Frame at a Time'); xlabel('Distance , x/L'); ylabel('Displacement, y(x,t)/Y_0'); pause t=linspace(0,8*pi,257); sine=sin(t); for k=1:2:65 ys=y*sine(1,k); plot(xoverL, ys,'k','LineWidth',[1.5]); hold on set(texthandle2,'String', 'Press Enter to Animate One Frame at a Time'); pause end hold off set(texthandle2,'String','Press Enter to Continue'); pause figure(3);clf; axis([-.3 1.3 -2*lim 2*lim]); hold on set(gca,'Box','on') patch(xp1,yp1,'r'); hold on patch(xp2,yp1,'r'); hold on xlabel('Distance, x/L') ylabel('Displacement, y(x,t)T/P_0L') hold on x1=[0 1];y1=[0 0]; plot(x1,y1) hold on texthandle2=text(-.1,-1.6*lim,['f/f_1 = ',num2str(fratio)]); ys=y*0; L=plot(xoverL,ys,'k','EraseMode','xor','LineWidth',[2.5]); ydata=[ys]; texthandl=text(0,1.5*lim,'Press Enter to Animate'); xlabel('Distance , x/L') ylabel('Displacement, y(x,t)T/P_0L') pause set(texthandl,'String',' '); t=linspace(0,8*pi,257); sine=sin(t); for k=2:257 ys=y*sine(1,k);; set(L,'Ydata',ys) pause(.05) ydata=[ydata ;ys]; end hold off set(texthandl,'String','Press Enter to Continue'); pause figure(4);clf;%3-d plot [X,T]=meshgrid(xoverL,t); mesh(X,T,ydata) axis([0 1 0 30 -1.2*lim 1.2*lim]); text(1,5,-lim,'Dimensionless Time, \omegat','Rotation',-10) text(1,30,-lim,'Dimensionless Distance, x/L','Rotation',32) zlabel('Displacement, y(x,t)T/P_0L') text(0,4,1.1*lim,'Press Enter to Continue','Rotation',-10) view(120,30) pause