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

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

runimpv
```function runimpv
% Example:  runimpv
% ~~~~~~~~~~~~~~~~~
% This is a driver program for the
% earthquake example.
%
% User m functions required:
%    fouaprox, imptp, hsmck,
%    shkbftss, lintrp

% Make the undamped period about one
% second long
m=1; k=36;

% Use damping equal to 5 percent of critical
c=.05*(2*sqrt(m*k));

% Choose a period equal to length of
% Imperial Valley earthquake data
prd=53.8;

nft=1024; tmin=0; tmax=prd;
ntimes=200; nsum=80; % ntimes=501; nsum=200;
tplt=linspace(0,prd,ntimes);
y20trm=fouaprox('imptp',prd,tplt,20); close
plot(tplt,y20trm,'-',tplt,imptp(tplt),'--');
xlabel('time, seconds');
ylabel('unitized displacement');
title('Result from a 20-Term Fourier Series')
figure(gcf);
disp('Press [Enter] to continue');
dumy=input('','s');
% print -deps 20trmplt

% Show how magnitudes of Fourier coefficients
% decrease with increasing harmonic order

fcof=fft(imptp((0:1023)/1024,1))/1024;
clf; plot(abs(fcof(1:100)));
xlabel('harmonic order');
ylabel('coefficient magnitude');
title(['Coefficient Magnitude in Base ' ...
'Motion Expansion']); figure(gcf);
disp('Press [Enter] to continue');
dumy=input('','s');
% print -deps coefsize

% Compute forced response
[t,ys,ys0,vs0,as]= ...
shkbftss(m,c,k,'imptp',prd,nft,nsum, ...
tmin,tmax,ntimes);

% Compute homogeneous solution
[t,yh,ah]= ...
hsmck(m,c,k,-ys0,-vs0,tmin,tmax,ntimes);

% Obtain the combined solution
y=ys(:)+yh(:); a=as(:)+ah(:);
clf; plot(t,y,'-',t,yh,'--');
xlabel('time'); ylabel('displacement');
title('Total and Homogeneous Response');
legend('Total response','Homogeneous response');
figure(gcf);
disp('Press [Enter] to continue');
dumy=input('','s');
print -deps displac;

clf; plot(t,a,'-');
xlabel('time'); ylabel('acceleration')
title('Acceleration Due to Base Oscillation')
figure(gcf); print -deps accel

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

function y=fouaprox(func,per,t,nsum,nft)
%
% y=fouaprox(func,per,t,nsum,nft)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Approximation of a function by a Fourier
% series.
%
% func   - function being expanded
% per    - period of the function
% t      - vector of times at which the series
%          is to be evaluated
% nsum   - number of terms summed in the series
% nft    - number of function values used to
%          compute Fourier coefficients. This
%          should be an integer power of 2.
%          The default is 1024
%
% User m functions called:  none.
%----------------------------------------------

if nargin<5, nft=1024; end;
nsum=min(nsum,fix(nft/2));
c=fft(feval(func,per/nft*(0:nft-1)))/nft;
c(1)=c(1)/2; c=c(:); c=c(1:nsum);
w=2*pi/per*(0:nsum-1);
y=2*real(exp(i*t(:)*w)*c);

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

function ybase=imptp(t,period)
%
% ybase=imptp(t,period)
% ~~~~~~~~~~~~~~~~~~~~~
% This function defines a piecewise linear
% function resembling the ground motion of
% the earthquake which occurred in 1940 in
% the Imperial Valley of California. The
% maximum amplitude of base motion is
% normalized to equal unity.
%
% period - period of the motion
%          (optional argument)
% t      - vector of times between
%          tmin and tmax
% ybase  - piecewise linearly interpolated
%          base motion
%
% User m functions called:  lintrp
%----------------------------------------------

tft=[ ...
0.00    1.26    2.64    4.01    5.10 ...
5.79    7.74;   8.65    9.74   10.77 ...
13.06   15.07   21.60   25.49;  27.38 ...
31.56   34.94   36.66   38.03   40.67 ...
41.87;  48.40   51.04   53.80    0    ...
0       0       0 ]';
yft=[ ...
0       0.92   -0.25    1.00   -0.29 ...
0.46   -0.16;  -0.97   -0.49   -0.83 ...
0.95    0.86   -0.76    0.85;  -0.55 ...
0.36   -0.52   -0.38    0.02   -0.19 ...
0.08;  -0.26    0.24    0.00    0    ...
0       0       0 ]';
tft=tft(:); yft=yft(:);
tft=tft(1:24); yft=yft(1:24);
if nargin == 2
tft=tft*period/max(tft);
end
ybase=lintrp(tft,yft,t);

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

function [t,ys,ys0,vs0,as]=...
shkbftss(m,c,k,ybase,prd,nft,nsum, ...
tmin,tmax,ntimes)
%
% [t,ys,ys0,vs0,as]=...
%   shkbftss(m,c,k,ybase,prd,nft,nsum, ...
%            tmin,tmax,ntimes)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function determines the steady-state
% solution of the scalar differential equation
%
%    m*y''(t) + c*y'(t) + k*y(t) =
%                  k*ybase(t) + c*ybase'(t)
%
% where ybase is a function of period prd
% which is expandable in a Fourier series
%
% m,c,k     - Mass, damping coefficient, and
%             spring stiffness
% ybase     - Function or vector of
%             displacements equally spaced in
%             time which describes the base
%             motion over a period
% prd       - Period used to expand xbase in a
%             Fourier series
% nft       - The number of components used
%             in the FFT (should be a power
%             of two). If nft is input as
%             zero, then ybase must be a
%             vector and nft is set to
%             length(ybase)
% nsum      - The number of terms to be used
%             to sum the Fourier series
%             expansion of ybase. This should
%             not exceed nft/2.
% tmin,tmax - The minimum and maximum times
%             for which the solution is to
%             be computed
% t         - A vector of times at which
%             the solution is computed
% ys        - Vector of steady-state solution
%             values
% ys0,vs0   - Position and velocity at t=0
% as        - Acceleration ys''(t), if this
%             quantity is required
%
% User m functions called:  none.
%----------------------------------------------

if nft==0
nft=length(ybase); ybft=ybase(:)
else
tbft=prd/nft*(0:nft-1);
ybft=fft(feval(ybase,tbft))/nft;
ybft=ybft(:);
end
nsum=min(nsum,fix(nft/2)); ybft=ybft(1:nsum);
w=2*pi/prd*(0:nsum-1);
t=tmin+(tmax-tmin)/(ntimes-1)*(0:ntimes-1)';
etw=exp(i*t*w); w=w(:);
ysft=ybft.*(k+i*c*w)./(k+w.*(i*c-m*w));
ysft(1)=ysft(1)/2;
ys=2*real(etw*ysft); ys0=2*real(sum(ysft));
vs0=2*real(sum(i*w.*ysft));
if nargout > 4
ysft=-ysft.*w.^2; as=2*real(etw*ysft);
end

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

function [t,yh,ah]= ...
hsmck(m,c,k,y0,v0,tmin,tmax,ntimes)
%
% [t,yh,ah]=hsmck(m,c,k,y0,v0,tmin,tmax,ntimes)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Solution of
%     m*yh''(t) + c*yh'(t) + k*yh(t) = 0
% subject to initial conditions of
%     yh(0) = y0 and yh'(0) = v0
%
% m,c,k      -  mass, damping and spring
%               constants
% y0,v0      -  initial position and velocity
% tmin,tmax  -  minimum and maximum times
% ntimes     -  number of times to evaluate
%               solution
% t          -  vector of times
% yh         -  displacements for the
%               homogeneous solution
% ah         -  accelerations for the
%               homogeneous solution
%
% User m functions called:  none.
%----------------------------------------------

t=tmin+(tmax-tmin)/(ntimes-1)*(0:ntimes-1);
r=sqrt(c*c-4*m*k);
if r~=0
s1=(-c+r)/(2*m); s2=(-c-r)/(2*m);
g=[1,1;s1,s2]\[y0;v0];
yh=real(g(1)*exp(s1*t)+g(2)*exp(s2*t));
if nargout > 2
ah=real(s1*s1*g(1)*exp(s1*t)+ ...
s2*s2*g(2)*exp(s2*t));
end
else
s=-c/(2*m);
g1=y0; g2=v0-s*g1; yh=(g1+g2*t).*exp(s*t);
if nargout > 2
ah=real(s*(2*g2+s*g1+s*g2*t).*exp(s*t));
end
end

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

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