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

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

vdb
```function vdb
% Example: vdb
% ~~~~~~~~~~~~
%
% This program calculates the shear, moment,
% slope, and deflection of a variable depth
% indeterminate beam subjected to complex
% input data are defined in the program
% statements below.
%
% User m functions required:
%   sngf, trapsum

clear all; Problem=1;
if Problem == 1
Title=['Problem from Arbabi and Li'];
Printout=10;  % Output frequency
BeamLength=3; % Beam length
NoSegs=301;   % # of beam divisions for
% integration
% External concentrated loads and location
ExtForce= [-1]; ExtForceX=[.5];
% External ramp loads and range
%        q1  q2  x1  x2
ExtRamp=[-1  -1   1   2];
% Interior supports: initial displacement
% and location
IntSupX=    [1; 2]; IntSupDelta=[0; 0];
% End (left and right) conditions
EndCondVal= [0; 0; 0; 0];   % magnitude
% 1=shear,2=moment,3=slope,4=delta
EndCondFunc=[3; 4; 3; 4];
% 1=left end,2=right end
EndCondEnd= [1; 1; 2; 2];
% EI or beam depth specification
EIorDepth=1;  % 1=EI values specified
% 2=depth values specified
if EIorDepth == 1
% Discretize the parabolic haunch for the
% three spans
Width=1; E=1; a=0.5^2; Npts=100;
h1=0.5; k1=1; x1=linspace(0,1,Npts);
h2=1.5; k2=1; x2=linspace(1,2,Npts);
h3=2.5; k3=1; x3=linspace(2,3,Npts);
y1=(x1-h1).^2/a+k1; y2=(x2-h2).^2/a+k2;
y3=(x3-h3).^2/a+k3;
EIx=[x1 x2 x3]'; h=[y1 y2 y3]';
EIvalue=E*Width/12*h.^3;
mn=min(EIvalue); EIvalue=EIvalue./mn;
else
% Beam width and Young's modulus
BeamWidth=[]; BeamE=[]; Depth=[]; DepthX=[];
end
elseif Problem == 2
Title=['From Timoshenko and Young,', ...
' p 434, haunch beam'];
Printout=12; NoSegs=144*4+1; BeamLength=144;
ExtForce=[]; ExtForceX=[];
ExtRamp=[-1 -1 0 108];
IntSupX=[36; 108]; IntSupDelta=[0; 0];
EndCondVal=[0; 0; 0; 0];
EndCondFunc=[2; 4; 2; 4];
EndCondEnd= [1; 1; 2; 2]; EIorDepth=2;
if EIorDepth == 1
EIvalue=[]; EIx=[];
else
BeamWidth=[1]; BeamE=[1];
% Discretize the parabolic sections
a=36^2/5; k=2.5; h1=0; h2=72; h3=144;
N1=36; N2=72; N3=36;
x1=linspace(  0, 36,N1); y1=(x1-h1).^2/a+k;
x2=linspace( 36,108,N2); y2=(x2-h2).^2/a+k;
x3=linspace(108,144,N3); y3=(x3-h3).^2/a+k;
Depth=[y1 y2 y3]'; DepthX=[x1 x2 x3]';
% Comparison values
I=BeamWidth*Depth.^3/12; Imin=min(I); L1=36;
k1=BeamE*Imin/L1; k2=k1/2; k3=k1;
t0=10.46/k1; t1=15.33/k1; t2=22.24/k1;
t3=27.95/k1;
fprintf('\n\nValues from reference');
fprintf('\n  Theta (x=  0): %12.4e',t0);
fprintf('\n  Theta (x= 36): %12.4e',t1);
fprintf('\n  Theta (x=108): %12.4e',t2);
fprintf('\n  Theta (x=144): %12.4e\n',t3);
end
end

% Load input parameters into matrices
Force=[ExtForce,ExtForceX];
NoExtForce=length(ExtForce);
[NoExtRamp,ncol]=size(ExtRamp);
IntSup=[IntSupDelta,IntSupX];
NoIntSup=length(IntSupX);
EndCond=[EndCondVal,EndCondFunc,EndCondEnd];
if EIorDepth == 1
BeamProp=[]; NoEIorDepths=length(EIx);
EIdata=[EIvalue EIx];
else
BeamProp=[BeamWidth BeamE];
NoEIorDepths=length(DepthX);
EIdata=[Depth DepthX];
end

% more on

% Output input data
label1=['shear     ';'moment    '; ...
'slope     ';'deflection'];
label2=['left   ';'right  '];
fprintf('\n\nAnalysis of a Variable Depth ');
fprintf('Elastic Beam');
fprintf('\n--------------------------------');
fprintf('---------');
fprintf('\n\n');
disp(['Title: ' Title]);
fprintf...
('\nBeam Length:                    %g', ...
BeamLength);
fprintf...
('\nNumber of integration segments: %g', ...
NoSegs);
fprintf...
('\nPrint frequency for results:    %g', ...
Printout);
fprintf('\n\nInterior Supports: (%g)', ...
NoIntSup);
if NoIntSup > 0
fprintf('\n  |   #   X-location   Deflection');
fprintf('\n  | --- ------------ ------------');
for i=1:NoIntSup
fprintf('\n  |%4.0f %12.4e %12.4e', ...
i,IntSup(i,2),IntSup(i,1));
end
end
fprintf('\n\nConcentrated Forces: (%g)', ...
NoExtForce);
if NoExtForce > 0
fprintf('\n  |   #   X-location        Force');
fprintf('\n  | --- ------------ ------------');
for i=1:NoExtForce
fprintf('\n  |%4.0f %12.4e %12.4e', ...
i,Force(i,2),Force(i,1));
end
end
if NoExtRamp > 0
fprintf('\n  | --- ------------ ------------');
fprintf(' ------------ ------------');
for i=1:NoExtRamp
fprintf('\n  |%4.0f %12.4e %12.4e ', ...
i,ExtRamp(i,3),ExtRamp(i,1));
fprintf('%12.4e %12.4e', ...
ExtRamp(i,4),ExtRamp(i,2));
end
end
fprintf('\n\nEnd conditions:');
fprintf('\n  | End    Function           Value');
fprintf('\n  ');
fprintf('| ------ ----------  ------------\n');
for i=1:4
j=EndCond(i,3); k=EndCond(i,2);
strg=sprintf('  %12.4e',EndCond(i,1));
disp(['  | ' label2(j,:) label1(k,:) strg]);
end
if EIorDepth == 1
fprintf('\nEI values are specified');
fprintf('\n  |   #      X-start     EI-value')
fprintf('\n  | --- ------------ ------------');
for i=1:NoEIorDepths
fprintf('\n  |%4.0f %12.4e %12.4e', ...
i,EIdata(i,2),EIdata(i,1));
end
else
fprintf('\nDepth values are specified for ');
fprintf('rectangular cross section');
fprintf('\n  | Beam width:      %12.4e', ...
BeamProp(1));
fprintf('\n  | Young''s modulus: %12.4e', ...
BeamProp(2));
fprintf('\n  |');
fprintf('\n  |   #      X-start        Depth')
fprintf('\n  | --- ------------ ------------');
for i=1:NoEIorDepths
fprintf('\n  |%4.0f %12.4e %12.4e', ...
i,EIdata(i,2),EIdata(i,1));
end
end
disp(' ');

% Begin analysis
x=linspace(0,BeamLength,NoSegs)'; t=clock;
[V,M,Theta,Delta,Reactions]= ...
bmvardep(NoSegs,BeamLength,Force,ExtRamp, ...
EndCond,IntSup,EIdata,BeamProp);
t=etime(clock,t);

% Output results
disp(' ');
disp(['Solution time was ',num2str(t),' secs.']);
if NoIntSup > 0
fprintf('\nReactions at Internal Supports:');
fprintf('\n  |   X-location     Reaction');
fprintf('\n  | ------------ ------------');
for i=1:NoIntSup
fprintf('\n  | %12.8g %12.4e', ...
IntSup(i,2),Reactions(i));
end
end
fprintf('\n\nTable of Results:');
fprintf('\n  |  X-location        Shear');
fprintf('       Moment');
fprintf('        Theta        Delta');
fprintf('\n  | ----------- ------------ ');
fprintf('------------');
fprintf(' ------------ ------------');
if Printout > 0
for i=1:Printout:NoSegs
fprintf('\n  |%12.4g %12.4e %12.4e', ...
x(i),V(i),M(i));
fprintf(' %12.4e %12.4e',Theta(i),Delta(i));
end
disp(' ');
else
i=1; j=NoSegs;
fprintf('\n  |%12.4g %12.4e %12.4e', ...
x(i),V(i),M(i));
fprintf(' %12.4e %12.4e',Theta(i),Delta(i));
fprintf('\n  |%12.8g %12.4e %12.4e', ...
x(j),V(j),M(j));
fprintf(' %12.4e %12.4e',Theta(j),Delta(j));
end
fprintf('\n\n');
subplot(2,2,1);
plot(x,V,'k-'); grid; xlabel('x axis');
ylabel('Shear'); title('Shear Diagram');
subplot(2,2,2);
plot(x,M,'k-'); grid; xlabel('x axis');
ylabel('Moment'); title('Moment Diagram')
subplot(2,2,3);
plot(x,Theta,'k-'); grid; xlabel('x axis');
ylabel('Slope'); title('Slope Curve');
subplot(2,2,4);
plot(x,Delta,'k-'); grid; xlabel('y axis');
ylabel('Deflection');
title('Deflection Curve'); subplot
drawnow; subplot ; figure(gcf)
%print -deps vdb

% more off

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

function [V,M,Theta,Delta,Reactions]= ...
bmvardep(NoSegs,BeamLength,Force,ExtRamp, ...
EndCond,IntSup,EIdata,BeamProp)
% [V,M,Theta,Delta,Reactions]=bmvardep ...
% (NoSegs,BeamLength,Force,ExtRamp,EndCond, ...
% IntSup,EIdata,BeamProp)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function computes the shear, moment,
% slope, and deflection in a variable depth
% elastic beam having specified end conditions,
% intermediate supports with given
% allowing concentrated loads and linearly
%
% NoSegs     - number of beam divisions for
%              integration
% BeamLength - beam length
% Force      - matrix containing the magnitudes
%              and locations for concentrated
% ExtRamp    - matrix containing the end
%              magnitudes and end locations
% EndCond    - matrix containing the type of
%              end conditions, the magnitudes,
%              and whether values are for the
%              left or right ends
% IntSup     - matrix containing the location
%              and delta for interior supports
% EIdata     - either EI or depth values
% BeamProp   - either null or beam widths
%
% V          - vector of shear values
% M          - vector of moment values
% Theta      - vector of slope values
% Delta      - vector of deflection values
% Reactions  - reactions at interior supports
%
% User m functions required:
%----------------------------------------------

if nargin < 8, BeamProp=[]; end
% Evaluate function value coordinates and 1/EI
x=linspace(0,BeamLength,NoSegs)';
kk=oneovrei(x,EIdata,BeamProp);

% External load contributions to shear and
% moment interior to span and at right end

% Deflections and position of interior supports
ns=size(IntSup,1);
if ns > 0
ysprt=IntSup(:,1); r=IntSup(:,2);
snf=sngf(x,r,1);
else
ysprt=[]; r=[]; snf=zeros(NoSegs,0);
end

% Form matrix governing y''(x)
smat=kk(:,ones(1,ns+3)).* ...
[x,ones(NoSegs,1),snf,me];

% Integrate twice to get slope and deflection
% matrices
smat=trapsum(0,BeamLength,smat);
ymat=trapsum(0,BeamLength,smat);

% slope/deflection at the right end
ss=smat(NoSegs,ns+3); yy=ymat(NoSegs,ns+3);

% Equations to solve for left end conditions
% and internal reactions
ns4=ns+4; j=1:4; a=zeros(ns4,ns4);
b=zeros(ns4,1); js=1:ns; js4=js+4;

% Account for four independent boundary
% conditions.  Usually two conditions will be
% imposed at each end.
for k=1:4
val=EndCond(k,1); typ=EndCond(k,2);
wchend=EndCond(k,3);
if wchend==1
b(k)=val; row=zeros(1,4); row(typ)=1;
a(k,j)=row;
else
if typ==1     % Shear
a(k,j)=[1,0,0,0]; b(k)=val-vv;
if ns>0
a(k,js4)=sngf(BeamLength,r,0);
end
elseif typ==2 % Moment
a(k,j)=[BeamLength,1,0,0]; b(k)=val-mm;
if ns>0
a(k,js4)=sngf(BeamLength,r,1);
end
elseif typ==3 % Slope
a(k,j)=[smat(NoSegs,1:2),1,0];
b(k)=val-ss;
if ns>0
a(k,js4)=smat(NoSegs,3:ns+2);
end
else          % Deflection
a(k,j)=[ymat(NoSegs,1:2),BeamLength,1];
b(k)=val-yy;
if ns>0
a(k,js4)=ymat(NoSegs,3:ns+2);
end
end
end
end

% Interpolate to assess how support deflections
% are affected by end conditions, external
if ns>0
a(js4,1)=interp1(x,ymat(:,1),r);
a(js4,2)=interp1(x,ymat(:,2),r);
a(js4,3)=r; a(js4,4)=ones(ns,1);
for j=1:ns-1
a(j+5:ns+4,j+4)= ...
interp1(x,ymat(:,j+2),r(j+1:ns));
end
end
b(js4)=ysprt-interp1(x,ymat(:,ns+3),r);

% Solve for unknown reactions and end conditions
c=a\b; v0=c(1); m0=c(2); s0=c(3); y0=c(4);
Reactions=c(5:ns+4);

% Compute the shear, moment, slope, deflection
% for all x
if ns > 0
V=v0+ve+sngf(x,r,0)*Reactions;
M=m0+v0*x+me+sngf(x,r,1)*Reactions;
Theta=s0+smat(:,ns+3)+smat(:,1:ns+2)* ...
[v0;m0;Reactions];
Delta=y0+s0*x+ymat(:,ns+3)+ ...
ymat(:,1:ns+2)*[v0;m0;Reactions];
else
Reactions=[]; V=v0+ve; M=m0+v0*x+me;
Theta=s0+smat(:,ns+3)+smat(:,1:2)*[v0;m0];
Delta=y0+s0*x+ymat(:,ns+3)+ ...
ymat(:,1:2)*[v0;m0];
end

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

(x,Force,ExtRamp)
%                       (x,Force,ExtRamp)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function computes the shear, moment,
% slope, and deflection in a uniform depth
% Euler beam which is loaded by a series of
% of shear, moment, slope and deflection all
% equal zero when x=0.
%
% x       - location along beam
% Force   - concentrated force matrix
% ExtRamp - distributed load matrix
%
% V       - shear
% M       - moment
% EITheta - slope
% EIDelta - deflection
%
% User m functions required: sngf
%----------------------------------------------

nf=size(Force,1); nr=size(ExtRamp,1);
nx=length(x); V=zeros(nx,1); M=V;
EITheta=V; EIDelta=V;
if nf > 0
F=Force(:,1); f=Force(:,2); V=V+sngf(x,f,0)*F;
M=M+sngf(x,f,1)*F;
if nargout > 2
EITheta=EITheta+sngf(x,f,2)*(F/2);
EIDelta=EIDelta+sngf(x,f,3)*(F/6);
end
end
if nr > 0
P=ExtRamp(:,1); Q=ExtRamp(:,2);
p=ExtRamp(:,3); q=ExtRamp(:,4);
S=(Q-P)./(q-p); sp2=sngf(x,p,2);
sq2=sngf(x,q,2); sp3=sngf(x,p,3);
sq3=sngf(x,q,3); sp4=sngf(x,p,4);
sq4=sngf(x,q,4);
V=V+sngf(x,p,1)*P-sngf(x,q,1)* ... % Shear
Q+(sp2-sq2)*(S/2);
M=M+sp2*(P/2)-sq2*(Q/2)+ ...       % Moment
(sp3-sq3)*(S/6);
if nargout > 2
EITheta=EITheta+sp3*(P/6)- ...  % EI*Theta
sq3*(Q/6)+(sp4-sq4)*(S/24);
EIDelta=EIDelta+sp4*(P/24)- ... % EI*Delta
sq4*(Q/24)+(sngf(x,p,5)- ...
sngf(x,q,5))*(S/120);
end
end

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

function val=oneovrei(x,EIdata,BeamProp)
% [val]=oneovrei(x,EIdata,BeamProp)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function computes 1/EI by piecewise
% linear interpolation through a set of data
% values.
%
% x        - location along beam
% EIdata   - EI or depth values
% BeamProp - null or width values
%
% val      - computed value for 1/EI
%
% User m functions required: none
%----------------------------------------------

if size(EIdata,1) < 2  % uniform depth case
v=EIdata(1,1);
EIdata=[v,min(x);v,max(x)];
end
if ( nargin > 2 ) & ( sum(size(BeamProp)) > 0)
% Compute properties assuming the cross
% section is rectangular and EIdata(:,1)
% contains depth values
width=BeamProp(1); E=BeamProp(2);
EIdata(:,1)=E*width/12*EIdata(:,1).^3;
end
val=1./lintrp(EIdata(:,2),EIdata(:,1),x);

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

function y=sngf(x,x0,n)
% y=sngf(x,x0,n)
% ~~~~~~~~~~~~~~
%
% This function computes the singularity
% function defined by
%    y=<x-x0>^n for n=0,1,2,...
%
% User m functions required: none
%----------------------------------------------

if nargin < 3, n=0; end
x=x(:); nx=length(x); x0=x0(:)'; n0=length(x0);
x=x(:,ones(1,n0)); x0=x0(ones(nx,1),:); d=x-x0;
s=(d>=zeros(size(d))); v=d.*s;
if n==0
y=s;
else
y=v;
for j=1:n-1; y=y.*v; end
end

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

function v=trapsum(a,b,y,n)
%
% v=trapsum(a,b,y,n)
% ~~~~~~~~~~~~~~~~~~
%
% This function evaluates:
%
%   integral(a=>x, y(x)*dx) for a<=x<=b
%
% by the trapezoidal rule (which assumes linear
% function variation between succesive function
% values).
%
% a,b - limits of integration
% y   - integrand that can be a vector-valued
%       function returning a matrix such that
%       function values vary from row to row.
%       It can also be input as a matrix with
%       the row size being the number of
%       function values and the column size
%       being the number of components in the
%       vector function.
% n   - the number of function values used to
%       perform the integration.  When y is a
%       matrix then n is computed as the number
%       of rows in matrix y.
%
% v   - integral value
%
% User m functions called:  none
%----------------------------------------------

if isstr(y)
% y is an externally defined function
x=linspace(a,b,n)'; h=x(2)-x(1);
Y=feval(y,x); % Function values must vary in
% row order rather than column
% order or computed results
% will be wrong.
m=size(Y,2);
else
% y is column vector or a matrix
Y=y; [n,m]=size(Y); h=(b-a)/(n-1);
end
v=[zeros(1,m); ...
h/2*cumsum(Y(1:n-1,:)+Y(2:n,:))];

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

% function y=lintrp(xd,yd,x)
% See Appendix B
```