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

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

cbfreq
```function cbfreq
% Example:  cbfreq
% ~~~~~~~~~~~~~~~~
% 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, beamanim, plotsave, inputv

clear, close; 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
% 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, beamanim, 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;
if norm([h,tmax])==0 | isnan(h), return, end
disp(' ')

[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'), figure(gcf)
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)
colormap([1 1 1])
axis([0,1,0,tmax,-zmax,zmax])
xlabel('x axis'); ylabel('time')
zlabel('deflection')
title(['Cantilever Beam Deflection ' ...
'for Varying Position and Time'])
figure(gcf);
print -deps endpos2
disp(' '), disp(['Press [Enter] to animate',...
' the beam motion'])
pause

titl='Cantilever Beam Animation';
xlab='x axis'; ylab='displacement';
beamanim(hdist,x,0.1,titl,xlab,ylab), close
end

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

% function beamanim(x,u,tpause,titl,xlabl,ylabl)
% See Appendix B

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

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

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

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