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

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

cablinea
```function cablinea
% Example: cablinea
% ~~~~~~~~~~~~~~~~~
% This program uses modal superposition to
% compute the dynamic response of a cable
% suspended at one end and free at the other.
% The cable is given a uniform initial
% velocity. Time history plots and animation
% of the motion are provided.
%
% User m functions required:
%    cablemk, udfrevib, canimate

% Initialize graphics
hold off; axis('normal'); close;

% Set physical parameters
n=30; gravty=1.; masses=ones(n,1)/n;
lengths=ones(n,1)/n;

% Obtain mass and stiffness matrices
[m,k]=cablemk(masses,lengths,gravty);

% for solution
dsp=zeros(n,1); vel=ones(n,1);
tmin=0; tmax=10; ntim=30;

% Compute the solution by modal superposition
[t,u,modvc,natfrq]=...
udfrevib(m,k,dsp,vel,tmin,tmax,ntim);

% Interpret results graphically
nt1=sum(t<=tmin); nt2=sum(t<=tmax);
u=[zeros(ntim,1),u];
y=cumsum(lengths); y=[0;y(:)];

% Plot deflection surface
disp(' '), disp('TRANSVERSE MOTION OF A CABLE')
surf(y,t,u); xlabel('y axis'); ylabel('time');
zlabel('transverse deflection');
title('Surface Showing Cable Deflection');
colormap(gray), brighten(.7)
view([30,30]); figure(gcf);
disp(['Press [Enter] to see the cable ',...
'position at two times'])
pause,  print -deps surface

% Show deflection configuration at two times
% Use closer time increment than was used
% for the surface plots.
mtim=4*ntim;
[tt,uu,modvc,natfrq]=...
udfrevib(m,k,dsp,vel,tmin,tmax,mtim);
uu=[zeros(mtim,1),uu];
tp1=.1*tmax; tp2=.2*tmax;
s1=num2str(tp1); s2=num2str(tp2);
np1=sum(tt<=tp1); np2=sum(tt<=tp2);
u1=uu(np1,:); u2=uu(np2,:);
yp=flipud(y(:)); ym=max(yp);
plot(u1,yp,'-k',u2,yp,'--k');
ylabel('distance from bottom');
xlabel('transverse displacement');
title(['Cable Transverse Deflection ' ...
'at t = ',s1,' and t = ',s2]);
legend('t = 1', 't = 2');

xm=.2*max([u1(:);u2(:)]);
ntxt=int2str(n); n2=1+fix(n/2);
str=strvcat(...
'The cable was initially vertical and was',...
'given a uniform transverse velocity.',...
['A ',ntxt,' link model was used.']);
text(xm,.9*ym,str), figure(gcf);
disp(['Press [Enter] to show the time ',...
'response at the middle and free end'])
pause, print -deps twoposn

% Plot time history for the middle and the end
clf; plot(tt,uu(:,n2),'--k',tt,uu(:,n+1),'-k');
xlabel('dimensionless time');
ylabel('transverse displacement');
title(['Position versus Time for the ' ...
'Cable Middle and End'])
legend('Midpoint','Lower end');
figure(gcf);
disp('Press [Enter] for a motion trace')
pause, print -deps 2timhist

% Plot animation of motion history
clf; canimate(y,u,t,0,.5*max(t),1);
print -deps motntrac
disp('Press [Enter] to finish'), pause, close;

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

function [m,k]=cablemk(masses,lngths,gravty)
%
% [m,k]=cablemk(masses,lngths,gravty)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Form the mass and stiffness matrices for
% the cable.
%
% masses     - vector of masses
% lngths     - vector of link lengths
% gravty     - gravity constant
% m,k        - mass and stiffness matrices
%
% User m functions called:  none.
%----------------------------------------------

m=diag(masses);
b=flipud(cumsum(flipud(masses(:))))* ...
gravty./lngths;
n=length(masses); k=zeros(n,n); k(n,n)=b(n);
for i=1:n-1
k(i,i)=b(i)+b(i+1); k(i,i+1)=-b(i+1);
k(i+1,i)=k(i,i+1);
end

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

function [t,u,mdvc,natfrq]=...
udfrevib(m,k,u0,v0,tmin,tmax,nt)
%
% [t,u,mdvc,natfrq]= ...
%              udfrevib(m,k,u0,v0,tmin,tmax,nt)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function computes undamped natural
% frequencies, modal vectors, and time response
% by modal superposition.  The matrix
% differential equation and initial conditions
% are
%
%    m u'' + k u = 0,  u(0) = u0, u'(0) = v0
%
% m,k       - mass and stiffness matrices
% u0,v0     - initial position and velocity
%             vectors
% tmin,tmax - time limits for solution
%             evaluation
% nt        - number of times for solution
% t         - vector of solution times
% u         - matrix with row j giving the
%             system response at time t(j)
% mdvc      - matrix with columns which are
%             modal vectors
% natfrq    - vector of natural frequencies
%
% User m functions called:  none.
%----------------------------------------------

% Call function eig to compute modal vectors
% and frequencies
[mdvc,w]=eig(m\k);
[w,id]=sort(diag(w)); w=sqrt(w);

% Arrange frequencies in ascending order
mdvc=mdvc(:,id); z=mdvc\[u0(:),v0(:)];

% Generate vector of equidistant times
t=linspace(tmin,tmax,nt);

% Evaluate the displacement as a
% function of time
u=(mdvc*diag(z(:,1)))*cos(w*t)+...
(mdvc*diag(z(:,2)./w))*sin(w*t);
t=t(:); u=u'; natfrq=w;

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

function canimate(y,u,t,tmin,tmax,norub)
%
% canimate(y,u,t,tmin,tmax,norub)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function draws an animated plot of
% data values stored in array u. The
% different columns of u correspond to position
% values in vector y. The successive rows of u
% correspond to different times. Parameter
% tpause controls the speed of the animation.
%
% u         - matrix of values for which
%             animated plots of u versus y
%             are required
% y         - spatial positions for different
%             columns of u
% t         - time vector at which positions
%             are known
% tmin,tmax - time limits for graphing of the
%             solution
% norub     - parameter which makes all
%             position images remain on the
%             screen. Only one image at a
%             time shows if norub is left out.
%             A new cable position appears each
%             time the user presses any key
%
% User m functions called:  none.
%----------------------------------------------

% If norub is input,
%   all images are left on the screen
if nargin < 6
rubout = 1;
else
rubout = 0;
end

% Determine window limits
umin=min(u(:)); umax=max(u(:)); udif=umax-umin;
uavg=.5*(umin+umax);
ymin=min(y); ymax=max(y); ydif=ymax-ymin;
yavg=.5*(ymin+ymax);
ywmin=yavg-.55*ydif; ywmax=yavg+.55*ydif;
uwmin=uavg-.55*udif; uwmax=uavg+.55*udif;
n1=sum(t<=tmin); n2=sum(t<=tmax);
t=t(n1:n2); u=u(n1:n2,:);
u=fliplr (u); [ntime,nxpts]=size(u);

hold off; cla; ey=0; eu=0; axis('square');
axis([uwmin,uwmax,ywmin,ywmax]);
axis off; hold on;
title('Trace of Linearized Cable Motion');

% Plot successive positions
for j=1:ntime
ut=u(j,:); plot(ut,y,'-k');
figure(gcf); pause(.5);

% Erase image before next one appears
if rubout & j < ntime, cla, end
end```