No BSD License  

Highlights from
Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

image thumbnail

Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

by

Howard Wilson

 

14 Oct 2002 (Updated )

Companion Software (amamhlib)

[t,x,y,theta,cptim]=cablenl
function [t,x,y,theta,cptim]=cablenl   
% [t,x,y,theta,cptim]=cablenl
% Example: cablenl
% ~~~~~~~~~~~~~~~~
% Numerical integration of the matrix 
% differential equations for the nonlinear 
% dynamics of a cable of rigid links with
% the outer ends of the cable fixed. 
%
% t     - time vector for the solution
% x,y   - matrices with nodal coordinates
%         stored in the columns. The time
%         history of point j is in x(:,j)
%         and y(:,j) 
% theta - matrix with inclination angles
%         stored in the columns
% cptim - number of seconds to integrate
%         the equations of motion
%
% User m functions required: 
%         plotmotn, equamo

clear all; close;

% Make variables global for use by 
% function equamo
global first_ n_ m_ len_ grav_ b_ mas_ py_

fprintf('\nNONLINEAR DYNAMICS OF A ')
fprintf('FALLING CABLE\n')
fprintf(...
'\nNote: The calculations take awhile\n')

% Set up data for a cable of n_ links, 
% initially arranged in a triangular 
% deflection configuration.

% parameter controlling initialization of 
% auxiliary parameters used in function 
% equamo
first_=1; 
% number of links in the cable
n_=8; n=n_; nh=n_/2;           
% vector of lengths and gravity constant
len_=1/n*ones(n,1); grav_=1;   
% vector of mass constants
m_=ones(1,n_)/n_;              

% initial position angles
th0=pi/4*[ones(nh,1);-ones(nh,1)]; 
td0=zeros(size(th0)); z0=[th0;td0];

% time limits, integration tolerances, 
% and the number of solution points
tmin=0; tmax=25; nt=201;
t=linspace(0,tmax,nt)';
tolrel=1e-6; tolabs=1e-8; len=len_;

% Perform the numerical integration using a 
% variable stepsize Runge-Kutta integrator
tic; 
odetol=odeset('RelTol',tolrel,'AbsTol',tolabs);
[t,w]=ode45(@equamo,t,z0,odetol);
theta=w(:,1:n); cptim=toc; 

% Compute node point coordinates
Z=[zeros(nt,1),repmat(len',nt,1).*exp(i*theta)];
Z=cumsum(Z.').'; x=real(Z); y=imag(Z); 

% Plot the horizontal position of the midpoint
clf; plot(t,x(:,1+n_ /2));  
ylabel('x coordinate'); xlabel('time')
title(['Horizontal Position of the ' ...
       'Cable Midpoint'])
grid on; figure(gcf); 
% print -deps xmidl

disp(' '), disp(...
'Press [Enter] to see the error growth curve');
pause, close

% Show error growth indicated by symmetry 
% loss of the vertical deflection symmetry.
% An approximately linear trend on the semilog
% plot indicates exponential growth of the error.
unsymer=sqrt(sum((y-y(:,end:-1:1)).^2,2));
hold off; axis('normal'); clf;
semilogy(t,unsymer); 
xlabel('time'); ylabel('asymmetry error');
title(['Growing Loss of Symmetry in ' ...
       'Vertical Deflection']);
grid on; figure(gcf);
% print -deps unsymerr

disp(' '), disp(...
'Press [Enter] to see the response animation');

% Show animation of the cable response 
disp(' ')
disp('The  motion can be animated or a trace')
disp('can be shown for successive positions')
disp(['between t = ',num2str(tmin),...
      ' and t = ',num2str(tmax)])

% Plot the position for different times limits
titl='CABLE MOTION FOR T = ';
while 1
  disp(' '), disp(...
  ['Choose a plot option (1 <=> animate, ',...
   ' 2 <=> trace,'])
  opt=input('3 <=> stop)  > ? ');
  if opt==3, break, end
  disp(' '), disp(...
  'Give a time vector such as 0:.1:15')
  Tp=input('Time vector > ? ','s'); 
  if isempty(Tp), break, end
  tp=eval(Tp); tp=tp(:); T=[titl,Tp];
  xp=interp1q(t,x,tp); yp=interp1q(t,y,tp);
  if opt ==1, plotmotn(xp,yp,T)
  else, plotmotn(xp,yp,T,1), end
end
fprintf('\nAll Done\n')

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

function plotmotn(x,y,titl,isave)
%
% plotmotn(x,y,titl,isave)
% ~~~~~~~~~~~~~~~~~~~~
% This function plots the cable time 
% history described by coordinate values 
% stored in the rows of matrices x and y.
%
% x,y   - matrices having successive rows 
%         which describe position 
%         configurations for the cable
% titl  - a title shown on the plots
% isave - parameter controlling the form 
%         of output. When isave is not input, 
%         only one position at a time is shown
%         in rapid succession to animate the
%         motion. If isave is given a value,
%         then successive are all shown at
%         once to illustrate a kinematic 
%         trace of the motion history.
%
% User m functions called:  none
%----------------------------------------------

% Set a square window to contain all 
% possible positions
[nt,n]=size(x); 
if nargin==4, save =1; else, save=0; end
xmin=min(x(:)); xmax=max(x(:));
ymin=min(y(:)); ymax=max(y(:)); 
w=max(xmax-xmin,ymax-ymin)/2;
xmd=(xmin+xmax)/2; ymd=(ymin+ymax)/2;  
hold off; clf; axis('normal'); axis('equal'); 
range=[xmd-w,xmd+w,ymd-w,ymd+w];
title(titl)
xlabel('x axis'); ylabel('y axis')
if save==0
  for j=1:nt
    xj=x(j,:); yj=y(j,:);
    plot(xj,yj,'-k',xj,yj,'ok');
    axis(range), axis off
    title(titl)
    figure(gcf), drawnow, pause(.1)
  end
  pause(2)
else
  hold off; close
  for j=1:nt
    xj=x(j,:); yj=y(j,:);
    plot(xj,yj,'-k',xj,yj,'ok');
    axis(range), axis off, hold on
  end
  title(titl)
  figure(gcf), drawnow, hold off, pause(2)
end

% Save plot history for subsequent printing 
% print -deps plotmotn

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

function zdot=equamo(t,z)
%
% zdot=equamo(t,z)
% ~~~~~~~~~~~~~~~~
% Equation of motion for a cable fixed at 
% both ends and loaded by gravity forces only
%
% t        current time value
% z        column vector defined by
%          [thet(t);theta'(t)]
% zdot     column vector defined by 
%          the concatenation
%          z'(t) = [theta'(t);theta''(t)]
%
% User m functions called:  none.
%----------------------------------------------

% Values accessed as global variables
global first_ n_ m_ len_ grav_ b_ mas_ py_

% Initialize parameters first time 
% function is called
if first_==1, first_=0;
% mass parameters
  m_=m_(:); b_=flipud(cumsum(flipud(m_))); 
  mas_=b_(:,ones(n_,1)); 
  mas_=tril(mas_)+tril(mas_,-1)';
% load effects from gravity forces
  py_=-grav_*(b_-b_(n_)); 
end

% Solve for zdot = [theta'(t); theta''(t)];
n=n_; len=len_; 
th=z(1:n); td=z(n+1:2*n); td2=td.*td;
x=len.*cos(th); y=len.*sin(th);

% Matrix of mass coefficients and 
% constraint conditions
amat=[[mas_.*(x*x'+y*y'),x,y];
      [x,y;zeros(2,2)]'];

% Right side vector involves applied forces 
% and inertial terms
bmat=x*y'; bmat=mas_.*(bmat-bmat');

% Solve for angular acceleration. 
% Most computation occurs here.
soln=amat\[x.*py_+bmat*td2; y'*td2; -x'*td2];

% Final result needed for use by the  
% numerical integrator
zdot=[td; soln(1:n)];

Contact us