# Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

strngrun(rundemo)
```function strngrun(rundemo)
%
% strngrun(rundemo)
% ~~~~~~~~~~~~~~~~
% This function illustrates propagation of
% waves in a tightly stretched string having
% given initial deflection. Calling strngrun
% with no input argument causes data to be
% executes a sample data case.
%
% User m functions called: strngwav animate

pltsav=0; % flag to save or not save graphs

close, disp(' ')
disp('WAVE PROPAGATION IN A STRING'), disp(' ')
if nargin==0 % Input data interactively
[a,len]=inputv(['Input wave speed (a) and ',...
'string length (len) > ? ']);
disp(' ')
disp(['Enter the number of interior ',...
'data points (the fixed'])
disp(['end point coordinates are ',...
n=input('? '); if isempty(n), return, end
xd=zeros(n+2,1); xd(n+2)=len;
yd=zeros(n+2,1); disp(' ')
disp(['The string stretches between ',...
'fixed endpoints at'])
disp(['x=0 and x=',num2str(len),'.']),disp(' ')
disp(['Enter ',num2str(n),...
' sets of x,y to specify interior'])
disp(['initial deflections ',...
'(one pair per line)'])
for j=2:n+1,[xd(j),yd(j)]=inputv; end;
disp(' ')
disp('Input tmax and the number of time steps')
[tmax,nt]=inputv('(Try len/a and 40) > ? ');
disp(' ')
disp('Specify position x=x0 where the time')
x0=input(...
'history is to be evaluated (try len/4) > ? ');
disp(' ')
disp('Specify time t=t0 when the deflection')
t0=input('curve is to be plotted > ? ');
disp(' ')
titl=input('Input a graph title > ? ','s');

else % Example for triangular initial deflection
a=1; len=1; tmax=len/a; nt=40;
xd=[0,.33,.5,.67,1]*len; yd=[0,0,-1,0,0];

% Different example for a truncated sine curve
% xd=linspace(0,len,351); yd=sin(3*pi/len*xd);
% k=find(yd<=0); xd=xd(k); yd=yd(k);

x0=0.25*len; t0=0.33*len/a;
titl='TRANSLATING WAVE OVER HALF A PERIOD';
end

nx=80; x=0:len/nx:len; t=0:tmax/nt:tmax;

h=max(abs(yd)); xplot=linspace(0,len,201);
tplot=linspace(0,max(t),251)';

[Y,X,T]=strngwav(xd,yd,x,t,len,a);
plot3(X',T',Y','k'); xlabel('x axis')
ylabel('time'), zlabel('y(x,t)'), title(titl)
if pltsav, print(gcf,'-deps','strngplot3'); end
drawnow, shg, disp(' ')

disp(['when t = ',num2str(t0)]), pause

[yt0,xx,tt]=strngwav(xd,yd,xplot,t0,len,a);
close; plot(xx(:),yt0(:),'k')
xlabel('x axis'), ylabel('y(x,t0)')
title(['DEFLECTION WHEN T = ',num2str(t0)])
axis([min(xx),max(xx),-h,h])
if pltsav, print(gcf,'-deps','strngyxt0'); end
drawnow, shg

disp(' ')
disp(['at x = ',num2str(x0)]), pause

yx0=strngwav(xd,yd,x0,tplot,len,a);
plot(tplot,yx0,'k')
xlabel('time'), ylabel('y(x0,t)')
title(...
['DEFLECTION HISTORY AT X = ',num2str(x0)])
axis([0,max(t),-h,h])
if pltsav, print(gcf,'-deps','strngyx0t'); end
drawnow, shg

disp(' ')
disp('over two periods of motion'), pause
x=linspace(0,len,101); t=linspace(0,4*len/a,121);
[Y,X,T]=strngwav(xd,yd,x,t,len,a);
titl='MOTION OVER TWO PERIODS';
animate(X(1,:),Y',titl,.1), pause(2)

if pltsav, print(gcf,'-deps','strnganim'); end

disp(' '), disp('All Done')

%===============================================

function [Y,X,T]=strngwav(xd,yd,x,t,len,a)
%
% [Y,X,T]=strngwav(xd,yd,x,t,len,a)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function computes the dynamic response of
% a tightly stretched string released from rest
% with a piecewise linear initial deflection. The
% string ends are fixed.
%
% xd,yd - data vectors defining the initial
%         deflection as a piecewise linear
%         function. xd values should be increasing
%         and lie between 0 and len
% x,t   - position and time vectors for which the
%         solution is evaluated
% len,a - string length and wave speed

if nargin<6, a=1; end; if nargin <5, len=1; end
xd=xd(:); yd=yd(:);  p=2*len;

% If end values are not zero, add these points
if xd(end)~=len, xd=[xd;len]; yd=[yd;0]; end
if xd(1)~=0, xd=[0;xd]; yd=[0;yd]; end
nd=length(xd);

% Eliminate any repeated abscissa values
k=find(diff(xd)==0); tiny=len/1e6;
if length(k)>0, xd(k)=xd(k)+tiny; end

% Extend the data definition for len < x < 2*len
xd=[xd;p-xd(nd-1:-1:1)]; yd=[yd;-yd(nd-1:-1:1)];
[X,T]=meshgrid(x,t); xp=X+a*T; xm=X-a*T;
shape=size(xp); xp=xp(:); xm=xm(:);

% Compute the general solution for a piecewise
% linear initial deflection
Y=(sign(xp).*interp1(xd,yd,rem(abs(xp),p),...
'linear','extrap')+sign(xm).*interp1(xd,yd,...
rem(abs(xm),p),'linear','extrap'))/2;
Y=reshape(Y,shape);

%===============================================

function animate(x,y,titl,tim,trace)
%
% animate(x,y,titl,tim,trace)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function performs animation of a 2D curve
% x,y - arrays with columns containing curve positions
%       for successive times. x can also be a single
%       vector if x values do not change. The animation
%       is done by plotting (x(:,j),y(:,j)) for
%       j=1:size(y,2).
% titl- title for the graph
% tim - the time in seconds between successive plots

if nargin<5, trace=0; else, trace=1; end;
if nargin<4, tim=.05; end
if nargin<3, trac=''; end; [np,nt]=size(y);
if min(size(x))==1, j=ones(1,nt); x=x(:);
else, j=1:nt; end; ax=newplot;
if trace, XOR='none'; else, XOR='xor'; end
r=[min(x(:)),max(x(:)),min(y(:)),max(y(:))];
%axis('equal') % Needed for an undistorted plot
axis(r), % axis('off')
curve = line('color','k','linestyle','-',...
'erase',XOR, 'xdata',[],'ydata',[]);
xlabel('x axis'), ylabel('y axis'), title(titl)
for k = 1:nt
set(curve,'xdata',x(:,j(k)),'ydata',y(:,k))
if tim>0, pause(tim), end, drawnow, shg
end

%===============================================

% function varargout=inputv(prompt)
% See Appendix B```