221: axis([uwmin,uwmax,ywmin,ywmax]);
222: axis off; hold on;
223: title(’Trace of Linearized Cable Motion’);
224:
225:
% Plot successive positions
226: for j=1:ntime
227: ut=u(j,:); plot(ut,y,’-’);
228: figure(gcf); pause(.5);
229:
230:
% Erase image before next one appears
231: if rubout & j < ntime, cla, end
232: end
7.2 Direct Integration Methods
Using stepwise integration methods to solve the structural dynamics equation pro-
vides an alternative to frequency analysis methods. If we invert the mass matrix and
save the result for later use, the n degree-of-freedom system can be expressed con-
cisely as a Þrst order system in 2n unknowns for a vector z =[x; v], where v is the
time derivative of x. The system can be solved by applying the variable step-size
differential equation integrator ode45 as indicated in the following function:
function [t,x]=strdynrk(t,x0,v0,m,c,k,functim)
% [t,x]=strdynrk(t,x0,v0,m,c,k,functim)
global Mi C K F n n1 n2
Mi=inv(m); C=c; K=k; F=functim;
n=size(m,1); n1=1:n; n2=n+1:2*n;
[t,z]=ode45(@sde,t,[x0(:);v0(:)]); x=z(:,n1);
%================================
function zp=sde(t,z)
global Mi C K F n n1 n2
zp=[z(n2); Mi*(feval(F,t)-C*z(n2)-K*z(n1))];

)
and

b
a
f(t)dt =
h
2
[f(a)+f(b)] +
h
2
12
[f

(a) − f

(b)] +
h
5
720
f
(4)
(
2
)
where a<
i
<band h = b −a. The Þrst formula, called the trapezoidal rule , gives
a zero truncation error term when applied to a linear function. Similarly, the second
formula, called the trapezoidal rule with end correction , has a zero Þnal term for a

t+h
t
Vdt.
For brevity we utilize a notation characterized by X(t)=X
0
, X(t + h)=X
1
,
˜
X = X
1
− X
0
. The trapezoidal rule immediately leads to

M +
h
2
C +
h
2
4
K

˜
V =

t+h
t
P (t)dt − h

which can be inverted once and used repeatedly if the step-size is not changed.
To integrate the forcing function we can use the midpoint rule [26] which states
that

b
a
P (t) dt = hP

a + b
2

+ O(h
3
).
© 2003 by CRC Press LLC
Solving for
˜
V yields
˜
V =

M +
h
2
C +
h
2
4
K


+
˜
V,X
1
= X
0
+
h
2
[V
0
+ V
1
]+O(h
3
).
A more accurate formula with truncation error of order h
5
can be developed from
the extended trapezoidal rule. This leads to
M
˜
V + C
˜
X + K

h
2
(
˜

˙
V
0
−
˙
V
1
]+O(h
5
).
Multiplying the last equation by M and employing the differential equation to reduce
the
˙
V
0
−
˙
V
1
terms gives
M
˜
X =
h
2
M[
˜
V +2V
0
]+

(M −
h
2
12
K)(C +
h
2
K)


˜
V
˜
X

=

hMV
0
+
h
2
12
(P
0
− P
1
)

Pdt−hKX

tion. Questions of computational efÞciency and numerical accuracy are examined for
two different step-sizes. Figures 7.7 and 7.8 present solution times as multiples of
the times needed for a modal response solution. The accuracy measures employed
© 2003 by CRC Press LLC
0 5 10 15 20 25 30 35 40 45 50
0
0.05
0.1
0.15
0.2
0.25
Solution Error For Implicit 2nd Order Integrator
time
solution error measure
h= 0.04, relative cputime= 34.6721
h= 0.08, relative cputime= 17.5615
Figure 7.7: Solution Error for Implicit 2nd Order Integrator
are described next. Note that the displacement response matrix has rows describ-
ing system positions at successive times. Consequently, a measure of the difference
between approximate and exact solutions is given by the vector
error_vector = \bsqrt(\bsum(((x_aprox-x_exact).ˆ2)’));
Typically this vector has small initial components (near t =0) and larger compo-
nents (near the Þnal time). The error measure is compared for different integrators
and time steps in the Þgures. Note that the fourth order integrator is more efÞcient
than the second order integrator because a larger integration step can be taken with-
out excessive loss in accuracy. Using h =0.4 for mckde4i achieved nearly the same
accuracy as that given by mckde2i with h =0.067. However, the computation time
for mckde2i was several times as large as that for mckde4i.
In the past it has been traditional to use only second order methods for solving
the structural dynamics equation. This may have been dictated by considerations on

8: %
9: % This program uses implicit second or fourth
10: % order integrators to compute the dynamical
11: % response of a cable which is suspended at
12: % one end and is free at the other end. The
13: % cable is given a uniform initial velocity.
14: % A plot of the solution error is given for
15: % two cases where approximate solutions are
16: % generated using numerical integration rather
17: % than modal response which is exact.
18: %
19: % User m functions required:
20: % mckde2i, mckde4i, cablemk, udfrevib,
21: % plterror
22:
23:
% Choose a model having twenty links of
24: % equal length
25:
26:
fprintf(
27: ’\nPlease wait: solution takes a while\n’)
28: clear all
29: n=20; gravty=1.; n2=1+fix(n/2);
30: masses=ones(n,1)/n; lengths=ones(n,1)/n;
31:
32:
% First generate the exact solution by
33: % modal superposition
34: [m,k]=cablemk(masses,lengths,gravty);

60: for j=1:j2
61: [t2,x2]=mckde2i(m,c,k,t0,dsp,vel,tfin,h2,i2);
62: end
63: tcp2=toc/j2; tr2=tcp2/tcpmr;
64:
65:
I2=5; H2=h/I2; tic;
66: for j=1:J2
67: [T2,X2]=mckde2i(m,c,k,t0,dsp,vel,tfin,H2,I2);
68: end
69: Tcp2=toc/J2; Tr2=Tcp2/tcpmr;
70:
71:
% Fourth order implicit results
72: i4=2; h4=h/i4; tic;
73: for j=1:j4
74: [t4,x4]=mckde4i(m,c,k,t0,dsp,vel,tfin,h4,i4);
75: end
76: tcp4=toc/j4; tr4=tcp4/tcpmr;
77:
78:
I4=1; H4=h/I4; tic;
79: for j=1:J4
80: [T4,X4]=mckde4i(m,c,k,t0,dsp,vel,tfin,H4,I4);
81: end
82: Tcp4=toc/J4; Tr4=Tcp4/tcpmr;
83:
84:
% Plot error measures for each solution
85: plterror(xmr,t2,h2,x2,T2,H2,X2,

113: % forcing function. This parameter
114: % should be omitted if no forcing
115: % function is used.
116: %
117: % Output:
118: %
119: % t time vector going from t0 to tmax
120: % in steps of
121: % x h*incout to yield a matrix of
122: % solution values such that row j
123: % is the solution vector at time t(j)
124: % tcp computer time for the computation
125: %
126: % User m functions called: none.
127: %
128:
129:
if (nargin > 9); force=1; else, force=0; end
130: if nargout ==3, tcp=clock; end
© 2003 by CRC Press LLC
131: hbig=h*incout;
132: t=(t0:hbig:tmax)’; n=length(t);
133: ns=(n-1)*incout; ts=t0+h*(0:ns)’;
134: xnow=x0(:); vnow=v0(:);
135: nvar=length(x0);
136: jrow=1; jstep=0; h2=h/2;
137:
138:
% Form the inverse of the effective
139: % stiffness matrix

168:
%=============================================
169:
170:
function [t,x,tcp] =
171: mckde4i(m,c,k,t0,x0,v0,tmax,h,incout,forc)
172: %
173: % [t,x,tcp]=
174: % mckde4i(m,c,k,t0,x0,v0,tmax,h,incout,forc)
175: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
© 2003 by CRC Press LLC
176: % This function uses a fourth order implicit
177: % integrator with fixed stepsize to solve the
178: % matrix differential equation
179: % m x’’ + c x’ + k x = forc(t)
180: % where m,c, and k are constant matrices and
181: % forc is an externally defined function.
182: %
183: % Input:
184: %
185: % m,c,k mass, damping and stiffness matrices
186: % t0 starting time
187: % x0,v0 initial displacement and velocity
188: % tmax maximum time for solution evaluation
189: % h integration stepsize
190: % incout number of integration steps between
191: % successive values of output
192: % forc externally defined time dependent
193: % forcing function. This parameter
194: % should be omitted if no forcing

221:
222:
% The forcing function is integrated using a
223: % 2 point Gauss rule
224: r3=sqrt(3); b1=h*(3-r3)/6; b2=h*(3+r3)/6;
225:
226:
% Initialize output matrix for x and other
227: % variables
228: xnow=x0(:); vnow=v0(:);
229: tnow=t0; zroforc=zeros(length(x0),1);
230:
231:
if force
232: fnow=feval(forc,tnow);
233: else
234: fnow=zroforc;
235: end
236: x=zeros(n,nvar); x(1,:)=xnow’; fnext=fnow;
237:
238:
% Main integration loop
239: for j=1:ns
240: tnow=t0+(j-1)*h; tnext=tnow+h;
241: if force
242: fnext=feval(forc,tnext);
243: di1=h12*(fnow-fnext);
244: di2=h2*(feval(forc,tnow+b1)+
245: feval(forc,tnow+b2));
246: z=mnv*[(di1+m*(h*vnow)); (di2-k*(h*xnow))];

272: %
273: % [m,k]=cablemk(masses,lngths,gravty)
274: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
275: % Form the mass and stiffness matrices for
276: % the cable.
277: %
278: % masses - vector of masses
279: % lngths - vector of link lengths
280: % gravty - gravity constant
281: % m,k - mass and stiffness matrices
282: %
283: % User m functions called: none.
284: %
285:
286:
m=diag(masses);
287: b=flipud(cumsum(flipud(masses(:))))*
288: gravty./lngths;
289: n=length(masses); k=zeros(n,n); k(n,n)=b(n);
290: for i=1:n-1
291: k(i,i)=b(i)+b(i+1); k(i,i+1)=-b(i+1);
292: k(i+1,i)=k(i,i+1);
293: end
294:
295:
%=============================================
296:
297:
function plterror(xmr,t2,h2,x2,T2,H2,X2,
298: t4,h4,x4,T4,H4,X4,tr2,Tr2,tr4,Tr4)

326: lg2=[’h= ’, num2str(H2),
327: ’, relative cputime= ’, num2str(Tr2)];
328: legend(lg1,lg2,2); figure(gcf);
329: disp(’Press [Enter] to continue’); pause
330: % print -deps deislne2
331:
332:
plot(t4,er4,’-’,T4,Er4,’ ’);
333: title([’Solution Error For Implicit ’,
334: ’4th Order Integrator’]);
335: xlabel(’time’);
336: ylabel(’solution error measure’);
337: lg1=[’h= ’, num2str(h4),
338: ’, relative cputime= ’, num2str(tr4)];
339: lg2=[’h= ’, num2str(H4),
340: ’, relative cputime= ’, num2str(Tr4)];
341: legend(lg1,lg2,2); figure(gcf);
342: % print -deps deislne4
343: disp(’ ’), disp(’All Done’)
344:
345:
%=============================================
346:
347:
% function [t,u,mdvc,natfrq]=
348: % udfrevib(m,k,u0,v0,tmin,tmax,nt)
349: % See Appendix B
350:
© 2003 by CRC Press LLC
Chapter 8

as infallible, and no numerical solution should be accepted as credible without ad-
equately investigating effects of parameter perturbation within uncertainty limits of
the parameters. To illustrate how sensitive a system can be to initial conditions, we

© 2003 by CRC Press LLC
might consider a very simple model concerning motion of a pendulum of length 
given an initial velocity v
0
starting from a vertically downward position. If v
0
ex-
ceeds 2
√
g, the pendulum will reach a vertically upward position and will go over
the top. If v
0
is less than 2
√
g, the vertically upward position is never reached.
Instead, the pendulum oscillates about the bottom position. Consequently, initial
velocities of 1.999
√
g and 2.001
√
g produce quite different system behavior with
only a tiny change in initial velocity. Other examples illustrating the difÞculties of
computing the response of nonlinear systems are cited below. These examples are
not chosen to discourage use of the powerful tools now available for numerical in-
tegration of differential equations. Instead, the intent is to encourage users of these
methods to exercise proper caution so that conÞdence in the reliability of results is

tion. Suppose a function y(x) satisÞes a differential equation of the form y

(x)=
f(x, y), subject to y(x
0
)=y
0
, where f is a known differentiable function. We
would like to compute an approximation of y(x
0
+ h) which agrees with a Taylor’s
series expansion up to a certain order of error. Hence,
y(x
0
+ h)=˜y(x
0
,h)+O(h
n+1
)
where O(h
n+1
) denotes a quantity which decreases at least as fast as h
n+1
for small
h. Taylor’s theorem allows us to write
y(x
0
+ h)=y(x
0
)+y

)+f
y
(x
0
,y
0
)f
0
]h
2
+ O(h
3
)
where f
0
= f(x
0
,y
0
). The last formula can be used to compute a second order
approximation ˆy(x
0
+h), provided the partial derivatives f
x
and f
y
can be evaluated.
However, this may be quite difÞcult since the function f(x, y) may not even be
known explicitly.
The idea leading to Runge-Kutta integration is to compute y(x

0
+ αh, y
0
+ βhf
0
)=f
0
+[f
x
(x
0
,y
0
)α + f
y
(x
0
,y
0
)f
0
β]h + O(h
2
),
we must have
˜y(x
0
+ h)=y
0
+ h[(k

)f
0
h +[f
x
(x
0
,y
0
)αk
1
+ f
y
(x
0
,y
0
)f
0
βk
1
]h
2
+ O(h
3
).
The last relation shows that
y(x
0
+ h)=˜y(x
0

0
+ h)=y(x
0
)+
1
2
[f
0
+ f(x
0
+ h, y
0
+ hf
0
)]h + O(h
3
).
Neglecting the truncation error O(h
3
) gives a difference approximation known as
Heun’s method [61], which is classiÞed as a second order Runge-Kutta method. Re-
ducing the step-size by h reduces the truncation error by about a factor of (
1
2
)
3
=
1
8
. Of course, the formula can be used recursively to compute approximations to

= f(x
0
,y
0
) ,k
2
= f(x
0
+
h
2
,y
0
+ k
1
h
2
),
k
3
= f(x
0
+
h
2
,y
0
+ k
2
h

time steps. This phenomenon, known as numerical instability, can be illustrated with
the simple differential equation
y

(t)=f(t, y)=λy
which has the solution y = ce
λt
. If the real part of λ is positive, the solution becomes
unbounded with increasing time. However, a pure imaginary λ produces a bounded
oscillatory solution, whereas the solution decays exponentially for real(λ) < 0.
Applying Heun’s method [43] gives
y(t + h)=y(t)

1+(λh)+
(λh)
2
2

.
This shows that at each integration step the next value of y is obtained by multiplying
the previous value by a factor
p =1+(λh)+
(λh)
2
2
,
which agrees with the Þrst three Taylor series terms of e
λh
. Clearly, the difference
relation leads to

= λy implies
y(t + h)=y(t)e
λh
and
e
λh
=
n

k=0
(λh)
k
k!
+ O(h
n+1
).
Consequently, points on the boundary of the stability region for a Runge-Kutta method
of order n are found by solving the polynomial
1 − e
ıθ
+
n

k=1
z
k
k!
=0
© 2003 by CRC Press LLC
for a dense set of θ-values ranging from 0 to 2π. Using MATLAB’s intrinsic function

3
4
real part of h*λ
imaginary part of h*λ
Stability Zone for Explicit Integrator of Order 6
Figure 8.2: Stability Zone for Explicit Integrator of Order 6
© 2003 by CRC Press LLC
MATLAB Example
Program rkdestab
1: % Example: rkdestab
2: % ~~~~~~~~~~~~~~~~~~
3: % This program plots the boundary of the region
4: % of the complex plane governing the maximum
5: % step size which may be used for stability of
6: % a Runge-Kutta integrator of arbitrary order.
7: %
8: % npts - a value determining the number of
9: % points computed on the stability
10: % boundary of an explicit Runge-Kutta
11: % integrator.
12: % xrang - controls the square window within
13: % which the diagram is drawn.
14: % [ -3, 3, -3, 3] is appropriate for
15: % the fourth order integrator.
16: %
17: % User m functions required: none
18:
19:
hold off; clf; close;
20: fprintf(’\nSTABILITY REGION FOR AN ’);

48: axis([-w,w,-w,w]); axis(’square’);
49: xlabel(’real part of h*\lambda’);
50: ylabel(’imaginary part of h*\lambda’);
51: ns=int2str(nordr);
52: st=[’Stability Zone for Explicit ’
53: ’Integrator of Order ’,ns];
54: title(st); grid on; figure(gcf);
55: % print -deps rkdestab
56: end
57:
58:
disp(’ ’); disp(’All Done’);
8.4 Discussion of Procedures to Maintain Accuracy by Varying
Integration Step-size
When we solve a differential equation numerically, our Þrst inclination is to seek
output at even increments of the independent variable. However, this is not the most
natural form of output appropriate to maintain integration accuracy. Whenever so-
lution components are changing rapidly, a small time step may be needed, whereas
using a small time step might be quite inefÞcient at times where the solution remains
smooth. Most modern ODE programs employ variable step-size algorithms which
decrease the integration step-size whenever some local error tolerance is violated and
conversely increase the step-size when the increase can be performed without loss of
accuracy. If results at even time increments are needed, these can be determined by
interpolation of the non-equidistant values. The differential equation integrators pro-
vide the capability to output results at an arbitrary vector of times over the integration
interval.
Although the derivation of algorithms to regulate step-size is an important topic,
development of these methods is not presented here. Several references [43, 46,
51, 61] discuss this topic with adequate detail. The primary objective in regulating
step-size is to gain computational efÞciency by taking as large a step-size as possible

nerable to accumulated errors from roundoff and arithmetic truncation. Such errors
usually render unreliable the results obtained sufÞciently far from the starting time.
This chapter concludes with the analysis of several realistic nonlinear problems
having certain properties of their exact solutions known. These known properties are
compared with numerical results to assess error growth. The Þrst problem involves
an inverted pendulum for which the loading function produces a simple exact dis-
placement function. Examples concerning top dynamics, a projectile trajectory, and
a falling chain are presented.
8.5 Example on Forced Oscillations of an Inverted Pendulum
The inverted pendulum in Figure 8.3 involves a weightless rigid rod of length l
which has a mass m attached to the end. Attached to the mass is a spring with
stiffness constant k and an unstretched length of γl. The spring has length l when
the pendulum is in the vertical position. Externally applied loads consist of a driv-
ing moment M (t), the particle weight, and a viscous damping moment cl
2
˙
θ. The
differential equation governing the motion of this system is
¨
θ = −(c/m)
˙
θ +(g/l)sin(θ)+M(t)/(ml
2
) − (2k/m)sin(θ)(1 − α/λ)
© 2003 by CRC Press LLC
k


mg
M(t)

solution for a nonlinear function. Let us assume that the driving moment M(τ)
produces a motion having the equation
θ
e
(τ)=θ
0
sin(ωτ)
for arbitrary θ
0
and ω. Then
˙
θ
e
(τ)=ωθ
0
cos(ωτ)
and
¨
θ
e
(τ)=−ω
2
θ
e
.
Consequently, the necessary driving moment is
P (τ)=−ω
2
θ
e

˙
θ). This is expressed as a Þrst
order matrix system by letting y
1
= θ, y
2
=
˙
θ, which gives
˙y
1
= y
2
, ˙y
2
= f(τ,y
1
,y
2
).
A function describing the system for solution by ode45 is provided at the end of this
section. Parameters θ
0
, ω
0
, α, ζ, and β are passed as global variables.
© 2003 by CRC Press LLC