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)

cantbfrq
function cantbfrq      
% Example:  cantbfrq
% ~~~~~~~~~~~~~~~~
% This program computes approximate natural 
% frequencies of a uniform depth cantilever 
% beam using finite difference and finite 
% element methods. Error results are presented 
% which demonstrate that the finite element 
% method is much more accurate than the finite 
% difference method when the same matrix orders
% are used in computation of the eigenvalues.
%
% User m functions required: 
%   cbfrqnwm, cbfrqfdm, cbfrqfem, frud, 
%   examplmo, animate, plotsave, inputv

clear;
fprintf('\n\n');
fprintf('CANTILEVER BEAM FREQUENCIES BY ');
fprintf('FINITE DIFFERENCE AND');
fprintf(...
'\n           FINITE ELEMENT APPROXIMATION\n');

fprintf('\nGive the number of frequencies ');
fprintf('to be computed');
fprintf('\n(use an even number greater ');
fprintf('than 2)\n'), n=input('? > ');
if rem(n,2) ~= 0, n=n+1; end

% Exact frequencies from solution of 
% the frequency equation
wex = cbfrqnwm(n,1e-12);

% Frequencies for the finite 
% difference solution
wfd = cbfrqfdm(n);

% Frequencies, modal vectors, mass matrix, 
% and stiffness matrix from the finite 
% element solution.
nelts=n/2; [wfe,mv,mm,kk] = cbfrqfem(nelts);
pefdm=(wfd-wex)./(.01*wex); 
pefem=(wfe-wex)./(.01*wex);

nlines=17; nloop=round(n/nlines);
v=[(1:n)',wex,wfd,pefdm,wfe,pefem];
disp(' '); lo=1;
t1=[' freq.     exact.        fdif.' ...
    '       fd. pct.'];
t1=[t1,'    felt.       fe. pct.'];
t2=['number     freq.         freq.' ...
    '        error  '];
t2=[t2,'    freq.        error  '];
while lo < n
  disp(t1),disp(t2); hi=min(lo+nlines-1,n);
  for j=lo:hi
    s1=sprintf('\n %4.0f %13.5e %13.5e', ...
               v(j,1),v(j,2),v(j,3));
    s2=sprintf(' %9.3f %13.5e %9.3f', ...
               v(j,4),v(j,5),v(j,6));
    fprintf([s1,s2]);
  end
  fprintf('\n\nPress [Enter] to continue\n\n');
  pause; lo=lo+nlines;
end
plotsave(wex,wfd,pefdm,wfe,pefem);
nfe=length(wfe); nmidl=nfe/2;
if rem(nmidl,2)==0, nmidl=nmidl+1; end
x0=zeros(nfe,1); v0=x0; w=0;
f1=zeros(nfe,1); f2=f1; f1(nfe-1)=1; 
f1(nmidl)=-5;
xsav=examplmo(mm,kk,f1,f2,x0,v0,wfe,mv);
close; fprintf('All Done\n');

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

function z=cbfrqnwm(n,tol)
%
% z=cbfrqnwm(n,tol)
% ~~~~~~~~~~~~~~~~~
% Cantilever beam frequencies by Newton's 
% method.  Zeros of 
%        f(z) = cos(z) + 1/cosh(z)  
% are computed.
%
% n   - Number of frequencies required
% tol - Error tolerance for terminating 
%       the iteration
% z   - Dimensionless frequencies are the 
%       squares of the roots of f(z)=0
%
% User m functions called:  none
%----------------------------------------------

if nargin ==1, tol=1.e-5; end

% Base initial estimates on the asymptotic 
% form of the frequency equation
zbegin=((1:n)-.5)'*pi; zbegin(1)=1.875; big=10;

% Start Newton iteration
while big > tol
  t=exp(-zbegin); tt=t.*t; 
  f=cos(zbegin)+2*t./(1+tt);
  fp=-sin(zbegin)-2*t.*(1-tt)./(1+tt).^2; 
  delz=-f./fp;
  z=zbegin+delz; big=max(abs(delz)); zbegin=z;
end
z=z.*z;

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

function [wfindif,mat]=cbfrqfdm(n)
%
% [wfindif,mat]=cbfrqfdm(n)
% ~~~~~~~~~~~~~~~~~~~~~~~~~
% This function computes approximate cantilever
% beam frequencies by the finite difference 
% method. The truncation error for the 
% differential equation and boundary 
% conditions are of order h^2.
%
% n       - Number of frequencies to be 
%           computed
% wfindif - Approximate frequencies in 
%           dimensionless form
% mat     - Matrix having eigenvalues which 
%           are the square roots of the 
%           frequencies
%
% User m functions called:  none
%----------------------------------------------

% Form the primary part of the frequency matrix
mat=3*diag(ones(n,1))-4*diag(ones(n-1,1),1)+...
    diag(ones(n-2,1),2); mat=(mat+mat');

% Impose left end boundary conditions 
% y(0)=0 and y'(0)=0
mat(1,[1:3])=[7,-4,1]; mat(2,[1:4])=[-4,6,-4,1];

% Impose right end boundary conditions
% y''(1)=0 and y'''(1)=0
mat(n-1,[n-3:n])=[1,-4,5,-2]; 
mat(n,[n-2:n])=[2,-4,2];

% Compute approximate frequencies and 
% sort these values
w=eig(mat); w=sort(w); h=1/n; 
wfindif=sqrt(w)/(h*h);

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

function [wfem,modvecs,mm,kk]= ...
                     cbfrqfem(nelts,mas,len,ei)
%
% [wfem,modvecs,mm,kk]=
%                    cbfrqfem(nelts,mas,len,ei)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Determination of natural frequencies of a 
% uniform depth cantilever beam by the Finite 
% Element Method.
%
%  nelts   - number of elements in the beam
%  mas     - total beam mass
%  len     - total beam length
%  ei      - elastic modulus times moment 
%            of inertia
%  wfem    - dimensionless circular frequencies
%  modvecs - modal vector matrix
%  mm,kk   - reduced mass and stiffness 
%            matrices
%
% User m functions called:  none
%----------------------------------------------

if nargin==1, mas=1; len=1; ei=1; end
n=nelts; le=len/n; me=mas/n; 
c1=6/le^2; c2=3/le; c3=2*ei/le;

% element mass matrix
masselt=me/420* ...
        [   156,   22*le,     54,  -13*le
          22*le,  4*le^2,  13*le, -3*le^2
             54,   13*le,    156,  -22*le
         -13*le, -3*le^2, -22*le,  4*le^2];

% element stiffness matrix
stifelt=c3*[ c1,  c2,  -c1,  c2
             c2,   2,  -c2,   1
            -c1, -c2,   c1, -c2
             c2,   1,  -c2,   2];

ndof=2*(n+1); jj=0:3; 
mm=zeros(ndof);  kk=zeros(ndof);

% Assemble equations
for i=1:n
  j=2*i-1+jj; mm(j,j)=mm(j,j)+masselt;
  kk(j,j)=kk(j,j)+stifelt;
end

% Remove degrees of freedom for zero 
% deflection and zero slope at the left end.
mm=mm(3:ndof,3:ndof); kk=kk(3:ndof,3:ndof);

% Compute frequencies
if nargout ==1
  wfem=sqrt(sort(real(eig(mm\kk))));
else
  [modvecs,wfem]=eig(mm\kk); 
  [wfem,id]=sort(diag(wfem));
  wfem=sqrt(wfem); modvecs=modvecs(:,id);
end

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

function [t,x]= ...
        frud(m,k,f1,f2,w,x0,v0,wn,modvc,h,tmax)
%
% [t,x]=frud(m,k,f1,f2,w,x0,v0,wn,modvc,h,tmax)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function employs modal superposition 
% to solve
%
%    m*x'' + k*x = f1*cos(w*t) + f2*sin(w*t)
%
% m,k    - mass and stiffness matrices
% f1,f2  - amplitude vectors for the forcing 
%          function
% w      - forcing frequency not matching any
%          natural frequency component in wn
% wn     - vector of natural frequency values
% x0,v0  - initial displacement and velocity 
%          vectors
% modvc  - matrix with modal vectors as its 
%          columns
% h,tmax - time step and maximum time for 
%          evaluation of the solution
% t      - column of times at which the 
%          solution is computed
% x      - solution matrix in which row j 
%          is the solution vector at 
%          time t(j)
%
% User m functions called:  none
%----------------------------------------------

t=0:h:tmax; nt=length(t); nx=length(x0); 
wn=wn(:); wnt=wn*t;

% Evaluate the particular solution.
x12=(k-(w*w)*m)\[f1,f2]; 
x1=x12(:,1); x2=x12(:,2);
xp=x1*cos(w*t)+x2*sin(w*t);

% Evaluate the homogeneous solution.
cof=modvc\[x0-x1,v0-w*x2]; 
c1=cof(:,1)'; c2=(cof(:,2)./wn)';
xh=(modvc.*c1(ones(1,nx),:))*cos(wnt)+...
   (modvc.*c2(ones(1,nx),:))*sin(wnt);

% Combine the particular and 
% homogeneous solutions.
t=t(:); x=(xp+xh)';

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

function x=examplmo(mm,kk,f1,f2,x0,v0,wfe,mv)
%
% x=examplmo(mm,kk,f1,f2,x0,v0,wfe,mv)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Evaluate the response caused when a downward
% load at the middle and an upward load at the 
% free end is applied.
%
% mm, kk - mass and stiffness matrices
% f1, f2 - forcing function magnitudes
% x0, v0 - initial position and velocity
% wfe    - forcing function frequency
% mv     - matrix of modal vectors
%
% User m functions called:  frud, animate, inputv
%----------------------------------------------

w=0; n=length(x0); t0=0; x=[];
s1=['\nEvaluate the time response from two',...
  '\nconcentrated loads. One downward at the',...
  '\nmiddle and one upward at the free end.'];
while 1
  fprintf(s1); fprintf('\n\n'); 
  fprintf('Input the time step and ');
  fprintf('the maximum time ');
  fprintf('\n(0.04 and 5.0) are typical.');
  fprintf(' Use 0,0 to stop\n');
  [h,tmax]=inputv; disp(' ');
  if norm([h,tmax])==0 | isnan(h), return; end

  [t,x]= ...
     frud(mm,kk,f1,f2,w,x0,v0,wfe,mv,h,tmax);
  x=x(:,1:2:n-1); x=[zeros(length(t),1),x];
  [nt,nc]=size(x); hdist=linspace(0,1,nc);

  clf; plot(t,x(:,nc),'k-');
  title('Position of the Free End of the Beam');
  xlabel('dimensionless time');
  ylabel('end deflection'); shg;
  disp('Press [Enter] for a surface plot of');
  disp(' transverse deflection versus x and t');
  pause
  print -deps endpos1 
  xc=linspace(0,1,nc); zmax=1.2*max(abs(x(:)));

  clf; surf(xc,t,x); view(30,35); 
  axis([0,1,0,tmax,-zmax,zmax]);
  xlabel('x axis'); ylabel('time'); 
  zlabel('deflection');
  title(['Cantilever Beam Deflection ' ...
       'for Varying Position and Time']); shg
  print -deps endpos2
  disp(['Press [Enter] to animate',...
        ' the beam motion']); pause
    
  titl='Cantilever Beam Animation'; 
  xlab='x axis'; ylab='displacement';
  animate(hdist,x,0.1,titl,xlab,ylab); close;
end

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

function animate(x,u,tpause,titl,xlabl,ylabl)
%
% animate(x,u,tpause,titl,xlabl,ylabl)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function draws an animated plot of data 
% values stored in array u.  The different 
% columns of u correspond to position values 
% in vector x.  The successive rows of u 
% correspond to different times. Parameter 
% tpause controls the speed of animation.
%
%  u      - matrix of values to animate plots 
%           of u versus x
%  x      - spatial positions for different 
%           columns of u
%  tpause - clock seconds between output of 
%           frames. The default is .1 secs 
%           when tpause is left out. When 
%           tpause=0, a new frame appears 
%           when the user presses any key.
%  titl   - graph title
%  xlabl  - label for horizontal axis
%  ylabl  - label for vertical axis
%
% User m functions called:  none
%----------------------------------------------
 
if nargin<6, ylabl=''; end; 
if nargin<5, xlabl=''; end
if nargin<4, titl=''; end; 
if nargin<3, tpause=.1; end;

[ntime,nxpts]=size(u); 
umin=min(u(:)); umax=max(u(:));
udif=umax-umin; uavg=.5*(umin+umax); 
xmin=min(x); xmax=max(x); 
xdif=xmax-xmin; xavg=.5*(xmin+xmax);
xwmin=xavg-.55*xdif; xwmax=xavg+.55*xdif;
uwmin=uavg-.55*udif; uwmax=uavg+.55*udif; clf;
axis([xwmin,xwmax,uwmin,uwmax]); title(titl);
xlabel(xlabl); ylabel(ylabl); hold on;

for j=1:ntime
  ut=u(j,:); plot(x,ut,'k-'); axis('off'); shg
  if tpause==0 
    pause; 
  else 
    pause(tpause); 
  end
  if j==ntime, break, else, cla; end
end
print -deps cntltrac
hold off; clf;

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

function plotsave(wex,wfd,pefd,wfe,pefem)
%
% function plotsave(wex,wfd,pefd,wfe,pefem)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function plots errors in frequencies 
% computed by two approximate methods.
%
% wex        - exact frequencies
% wfd        - finite difference frequencies
% wfe        - finite element frequencies
% pefd,pefem - percent errors by both methods
%
% User m functions called:  none
%----------------------------------------------

% plot results comparing accuracy 
% of both frequency methods
w=[wex(:);wfd(:);wfd]; 
wmin=min(w); wmax=max(w); n=length(wex);
wht=wmin+.001*(wmax-wmin); j=1:n;
 
semilogy(j,wex,'k-',j,wfe,'k--',j,wfd,'k:')
title('Cantilever Beam Frequencies');
xlabel('frequency number');
ylabel('frequency values');
legend('Exact freq.','Felt. freq.', ...
       'Fdif. freq.',2); figure(gcf); 
disp(['Press [Enter] for a frequency ',...
      'error plot']); pause 
print -deps beamfrq1 

plot(j,abs(pefd),'k:',j,abs(pefem),'k--');
title(['Cantilever Beam Frequency ' ...
       'Error Percentages']);
xlabel('frequency number'); 
ylabel('percent frequency error'); 
legend('Fdif. pct. error','Felt. pct. error',4);
figure(gcf);
disp(['Press [Enter] for a transient ',...
'response calculation']); pause
print -deps beamfrq2

Contact us