11.

DIRECT STEP-BY-STEP INTEGRATION METHODS

11.1. Introduction
In the previous three sections, the dynamic equilibrium equations for a multiple degree of freedom
system were solved by transforming them to a normal coordinate system based on the natural
vibrations modes of the structure. This transformation produced an uncoupled equilibrium
equation of motion for each vibration mode which could then be solved independently as a single
degree of freedom system.
However, the classical (i.e. non-complex) modal analysis method is based on a number of
fundamental assumptions. First, it is assumed that the damping satisfies modal orthogonality
requirements so that the modes can be uncoupled. Secondly, the superposition of modal
responses to obtain the total dynamic response is valid only for linear structural behavior. Finally,
the use of an eigenvalue analysis to determine the mode shapes and natural frequencies requires
that the mass and stiffness properties of the structure are frequency independent.
As an alternative, the coupled differential equations of motion expressed in the original (non-
modal) coordinate system can be integrated simultaneously in a step-by-step fashion. In this
approach, the equations of motion are satisfied at discrete time intervals rather than on a
continuous basis; the motion between time steps is assumed to follow some specific pattern.
Because of the discretization, the results from these methods are approximate, with the error
involved depending on the length of time step and the form of the motion assumed between
discrete intervals.
With the step-by-step solution method, there is no need to make any assumption regarding the
modal orthogonality of the damping matrix; however, the damping matrix must be formed
explicitly. Since modal response is not considered directly, it is also unnecessary to make any
assumption regarding the number of modes that make significant contributions. Possibly the
greatest advantage of the step-by-step time domain method, though, is that the structural property
matrices can be updated at each time step based on the current condition of the system; this
facilitates the consideration of nonlinear behavior. Since an eigenvalue solution is not required,
time step solution methods can also be advantageous in systems with a large number of degrees of
freedom which are subjected to a short duration transient load.
There are a large number of different direct integration schemes available, differing principally in
the manner in which the system is assumed to move between the discrete intervals. Only two of
the most common approaches will be discussed here.
11.2. Basic Finite Difference Schemes - Forward, Backward and Central Difference

11.2.1. Fundamental Concepts

Given a time-dependent function f (t ) , the first derivative of that function with respect to time is
defined as

df (t ) f ( t + ∆t ) − f ( t )
Eq. 11.1 = lim .
dt ∆t → 0 ∆t

11-1
 CE 804 - Structural Dynamics Page 11-2

For a finite time step of ∆t , Eq. 11.1 may be approximated by

df (t ) ∆f f ( t + ∆t ) − f ( t )
Eq. 11.2 ≈ = .
dt ∆t ∆t

Similarly, the second derivative of f ( t ) is approximately

d 2 f ( t ) 1  ∆f ( t + ∆t ) ∆f ( t ) 
≈ −
∆t  ∆t 
Eq. 11.3 .
dt 2 ∆t

In the finite difference approach, differential equations are solved by replacing the derivatives of a
function by expressions involving only the value of the function itself at discrete time steps and
subsequently solving for those values. For example, assume that the function f ( t ) is to be
determined at m discrete points that are equally spaced at intervals of ∆t , as shown in Fig. 11.1.
In the following discussion, expressions will be derived for the derivatives of the function at the ith
interval based on the value of the function at adjacent intervals.

f (t ) Point of Interest

f (t )

∆t

ti −1 ti ti +1 ti +2 t

Fig. 11.1 Function defined at discrete time intervals.

For all three schemes listed below, the error introduced into the solution by the discretization
process is in the order of ∆t 2 . Therefore, the accuracy of the solution can be increased by
reducing the length of the interval ∆t ; as ∆t becomes infinitely small, the approximate results
will converge on the exact solution.

11.2.2. Backward Difference Scheme

In the backward difference scheme, derivatives of the function at the ith interval are expressed in
terms of the value of the function at previous intervals. Adapting Eq. 11.2, the first derivative
becomes

∆f ( ti ) f (ti ) − f (ti −1 )
Eq. 11.4 = .
∆t ∆t
 CE 804 - Structural Dynamics Page 11-3

Likewise, the second derivative is obtained using Eq. 11.3 as

∆2 f ( ti ) 1  ∆f ( t i ) ∆f ( t i − 1 ) 
=  − 
∆t 2 ∆t  ∆t ∆t 
Eq. 11.5 .
[ ]
= 2 f (t i − 2 ) − 2 f ( t i − 1 ) + f (t i )
1
∆t

11.2.3. Forward Difference Scheme

In the forward difference scheme, derivatives of the function at the ith interval are expressed in
terms of the value of the function at the subsequent intervals. Using this approach, Eq. 11.2
becomes

∆f ( ti ) f ( ti +1 ) − f ( ti )
Eq. 11.6 =
∆t ∆t
while the second derivative is expressed as

∆2 f ( ti ) 1  ∆f (ti + 1 ) ∆f ( ti ) 
= −
∆t 2 ∆t  ∆t ∆t 
Eq. 11.7 .
[
= 2 f ( ti ) − 2 f ( ti + 1 ) + f ( ti + 2 )
1
∆t
]

11.2.4. Central Difference Scheme

In the central difference scheme, derivatives of the function at the ith interval are expressed in
terms of a time interval centered on the point in question. Consequently, Eq. 11.2 can be written

 ∆t   ∆t 
f  ti +  − f  ti − 
∆ f ( ti )  2  2
Eq. 11.8 = .
∆t ∆t
Since the values of the function are only determined at the intervals themselves and not mid way
in between, the values of the function at mid intervals are approximated by the expressions

 ∆t  1
Eq. 11.9a f  ti −  ≈
 2 2
[
f ( ti − 1 ) + f ( ti ) ]
and

 ∆t  1
Eq. 11.9b f  ti +  ≈
 2 2
[
f ( ti ) + f ( ti + 1 ) . ]
Substituting these expressions into Eq. 11.8 yields
 CE 804 - Structural Dynamics Page 11-4

∆ f ( ti )
Eq. 11.10
∆x
=
1
2 ∆t
[ ]
− f ( ti − 1 ) + f ( ti + 1 ) .

which, in turn, leads to

  ∆t   ∆t  
∆f  t +  ∆f  ti −  
∆2 f ( ti ) 1   i 2   2
=  − 
∆t 2
∆t  ∆t ∆t 
Eq. 11.11
 
1
[
= 2 f ( ti − 1 ) − 2 f ( ti ) + f ( ti + 1 )
∆t
]
Since the backward difference scheme provides the greatest accuracy for a given time step ∆t , it
is generally preferred as the basis for dynamic analysis. Consequently, subsequent discussions will
focus on the central difference scheme.

11.2.5. Dynamic Analysis Using the Central Difference Finite Difference Scheme*
Assume that the dynamic displacements of a structure at any time t are described by the vector
{ η( t )} where the jth element of this vector η j ( t ) corresponds to the displacement at the jth
degree of freedom for the structure. Using Eqs. 11.10 and 11.11, the velocity at the ith time
interval (i.e. t = ti ) can therefore be expressed as

Eq. 11.12 {η& (t )} = 2 1∆t [−{η(t )} + {η(t )}]

i i −1 i +1

while the acceleration would be

Eq. 11.13 {η&&(t )} = ∆1t [{η(t )} − 2 {η(t )} + {η(t )}] .

i 2 i −1 i i +1

The statement of dynamic equilibrium for the entire system at time ti can be expressed by the
matrix equation

Eq. 11.14 [m] {η&&(ti )} + [c] {η& (ti )} + [k ] {η(ti )} = {F (ti )}

which would include n coupled equations for a system with n degrees of freedom. Since this
statement of dynamic equilibrium is expressed in terms of conditions at the beginning of the
current time step (i.e. at t = ti ) when the displacements are known rather than at the end of the
time step (i.e. at t = ti + ∆t = ti +1 ) when the displacements are unknown, the central difference
scheme is referred to as an explicit finite difference method.

*
Reference: “Numerical Methods in Finite Element Analysis”, by Klause-Jürgen Bathe and Edward L. Wilson,
Prentice Hall, 1976
 CE 804 - Structural Dynamics Page 11-5

Substituting Eqs. 11.12 and 11.13 into Eq. 11.14, an expression can be obtained that has the
familiar form of

Eq. 11.15 [ K$ ] {η(t )} = {F$ }

i +1

[ ]
where K$ may be likened to a pseudo-stiffness matrix defined by

Eq. 11.16 [ K$ ] = ∆1t 2 [m] + 2 ∆t [c]

1

{}
and F$ may be thought of as a pseudo-force vector given by

Eq. 11.17 {F$ } = {F (t )} − [ k ] − ∆2t [m] {η(t )} −  ∆1t [m] − 2 1∆t [c] {η(t )} .
i 2 i 2 i −1

Assuming that the displacement vectors {η(t )} i and {η(t )}

i −1 are already known, the
displacement vector at the next time step (i.e. t = ti + ∆t = ti +1 ) can be solved for directly by

{η(t )} = [ K$ ] {F$ } .
−1
Eq. 11.18 i +1

It should be noted that the central difference method is only conditionally stable: that is, the
results will only approximate the true solution if the time step ∆t ≤ ∆tcrit where ∆tcrit is the
critical time step value given by

2
Eq. 11.19 ∆tcrit =
ωn

in which ω n is the highest natural frequency of a system with n degrees of freedom. If the time
step is larger than ∆tcrit , the error in the solution will grow progressively larger with each
successive time and will be unbounded.
• Analysis Procedure:
1. Start-Up Procedure: In practice, the initial conditions of the system (i.e. { η( 0)} , { η
& ( 0)} ,
and { η
&& ( 0)} ) are known or can be estimated. To calculate the displacements at the end of the
first time step (i.e. t = ∆t ), however, the displacements at two previous time steps are
required. To get around this problem, a set of fictitious displacements can be calculated for
t = − ∆t by solving for {η( ti −1 )} using Eqs. 11.12 and 13, giving

∆t 2 &&
Eq. 11.20 { η(− ∆t )} = { η( 0)} − ∆t { η& (0)} + 2
{ η( 0)} .
 CE 804 - Structural Dynamics Page 11-6

If only the initial displacements and velocities are known (i.e. { η( 0)} , { η& ( 0)} ), the initial
accelerations { η
&& ( 0)} can be determined from the dynamic equilibrium equation (Eq. 11.14).

2. General Procedure: To calculate the displacement vector at the end of the current time step,
(i.e. at t = ti +1 ) the following steps are performed:

• The external force vector at the start of the time step {F ( ti )} is defined.

• For a nonlinear system, the stiffness, damping, and possibly mass matrices can be updated
based on the current position of the system, {η(ti )} . For a linear system, the structural
property matrices need be defined only once at the start of the procedure.
• For the first time step, the a fictitious displacement vector {η( ti −1 )} is calculated as outlined
above. For subsequent time steps, the displacement vectors from the two previous time steps,
{η(ti )} and {η(ti −1 )} , will already be known at this stage.
• {}
The pseudo-force vector F$ is calculated using Eq. 11.17.

• [ ]
For a nonlinear system, the pseudo-stiffness matrix K$ is formed using Eq. 11.16. If the
mass and damping matrices do not change with time, [ K$ ] only has to formed once at the start
of the procedure.
• The displacements at the end of the time step, {η( ti +1 )} , are calculated using Eq. 11.18. Note
[ ]
that K$ only has to be inverted once if it remains constant, resulting in considerable
computational savings.
• If required, the velocity and acceleration vectors at the end of the time step can be calculated
using Eqs. 11.12 and 11.13. This step is not necessary unless these values are specifically
needed for some purpose.
11.3. Newmark’s Beta Method

11.3.1. Conventional Formulation*

Newmark’s method actually represents a whole series of solution methods that have a common
foundation. In this method, the acceleration is assumed to vary in a specified manner over the
time step. Dispalcement and velocity vectors at the beginning and end of each time step are then
related on the basis of the assumed acceleration pattern.
If the vectors describing the motion at the start of the time step (i.e at time t) are known and given
by {η} t , {η
& } , and {η
t
&&} , those at the end of the time step (i.e. at t = ti + ∆t ) may be expressed
t
as

*
Reference: “Numerical Methods in Finite Element Analysis”, by Klause-Jürgen Bathe and Edward L. Wilson,
Prentice Hall, 1976
 CE 804 - Structural Dynamics Page 11-7

 
Eq. 11.21a {η&&}t + ∆t =
1
β ( ∆t )
2 ({η} t + ∆t )
− {η}t −
1
β ( ∆t )
{η& } −  1 − 1 {η
t
 2β 
&&}
t

Eq. 11.21b {η& }t + ∆t = {η& }t + [(1 − γ ) {&&η}t + γ {&&η}t + ∆t ] (∆t )

 1 &&}  ( ∆t )2
Eq. 11.21c {η}t + ∆t = {η}t + {η& }t (∆t ) +  
− β {&&

η}t + β {η t + ∆t
 2 
The parameter β controls the manner in which the acceleration varies over the time step*. If
β = 0 , the acceleration at the beginning of the time step is assumed to remain constant over the
interval. For β = 1 4 , a constant acceleration is assumed that is equal to the average of the
acceleration values at the beginning and end of the time step. For β = 1 6 , a linearly varying
acceleration will be assumed.
While the linear acceleration method can produce superior accuracy for a given length of time
step, it is only conditionally stable; it will produce results that diverge from the true response if
∆t is too large. The constant-average acceleration method (β = 1 4 ), on the other hand, is
unconditionally stable, meaning that the errors involved will have an upper bound regardless of
the length of time step.
The fact that the constant-average acceleration is unconditionally stable does not imply that the
length of the time step is unimportant. If the time step gets too long, the periods of the resulting
motion will be elongated as compared to the true response. For reasonably accurate results,
therefore, the time step length should satisfy the criteria that ∆t ≤ Tmax 10 , where Tmax is the
period of the highest vibration mode which is likely to make a significant contribution to the
response.
The parameter γ controls the numerical (artificial) damping that is applied to the solution. The
numerical damping will be positive for γ > 1 2 , negative for γ < 1 2 , and absent for γ = 1 2 .
Numerical damping may be advisable in stiff, lightly damped structures to improve the
convergence rate of the solution. However, numerical damping will affect the amplitude of
calculated motion.
In Newmark’s method, the dynamic equilibrium equations are written in terms of the conditions at
the end of the current time step as follows:

Eq. 11.22 [m] {&&η}t + ∆t + [c] {η& }t + ∆t + [ k ] {η}t + ∆t = { F}t + ∆t .

Since the motion vectors at the end of the time step are initially unknown, this approach is
classified as an implicit direct integration method. Susbstituting the expressions in Eq. 11.21 into
Eq. 11.22 with considerable rearrangement yields a pseudo-static equation for the displacements
at the end of the time step given by

*
Reference: “Structural Dynamics by Finite Elements”, by William Weaver, Jr., and Paul R. Johnston, Prentice-
Hall, 1987, pg. 218
 CE 804 - Structural Dynamics Page 11-8

Eq. 11.23 [ K$ ] {η} t + ∆t {}

= F$

where the pseudo-stiffness matrix is defined as

γ
Eq. 11.24 [ K$ ] = [k ] + β ∆1t 2 [m] + β ∆t [c]

and the pseudo-load vector is given by

 1  1  
{F$ } = { F} t + ∆t
+ [ m] 
 β ∆t
2 { }t
η +
1
β ∆t
{η& }t +  − 1 {&&η}t  +
2β  
Eq. 11.25
 γ γ  ∆t  γ  
[ c]  {η}t +  − 1 {η& }t +  − 2 {&&η}t 
 β ∆t β  2 β  

[ ]
For a structure with constant structural properties, K$ only has to be formed and inverted once
at the start of the procedure. If the stiffness, mass, or damping matrix changes with position or
with time, it must be formed and inverted at every time step. The pseudo-load vector F$ must {}
be formed at every time step.
Solution Procedure:
• [ ]
If required for the current time step, form and invert the pseudo-stiffness matrix K$ based on
Eq. 11.24.

• {}
Form the pseudo-load vector F$ using Eq. 11.25.

• Calculate the displacements {η} t + ∆t at the end of the current time step using Eq. 11.23.

• Using these displacements, calculate the velocity and acceleration vectors at the end of the
time step using Eq. 11.21.

11.3.2. Incremental Formulation

For nonlinear systems it is often convenient to write the equations of motion in the incremental
form

Eq. 11.26 [m] {∆&&η} + [c] {∆η& } + [ k ] {∆η} = {∆F }

in which the incremental values over the period from time t to time t + ∆t are given by

Eq. 11.27a {∆η} = {η}t +∆t − {η}t

Eq. 11.27b {∆η& } = {η& }t +∆t − {η& }t

 CE 804 - Structural Dynamics Page 11-9

Eq. 11.27c {∆&&η} = {&&η}t +∆t − {&&η}t

and

Eq. 11.28 { ∆ F } = { F } t + ∆t − { F } t .
It is assumed that the displacement, velocity, and acceleration vectors are known at time t and that
the applied force vector is known at times t and t + ∆t . At this stage, the stiffness, mass and
damping matrices are based on the condition of the structure at time t; these properties are
initially assumed to remain constant over ∆t .
Substituting Eqs. 11.21 into Eq. 11.27 and rearranging, the incremental motion vectors
become
1 1 1
Eq. 11.29a {∆&&η} = {∆η} − {η& }t − {&&η}t
β ( ∆t )
2
β ( ∆t ) 2β

γ γ  γ 
Eq. 11.29b {∆η& } = {∆η} − {η}t −  − 1 (∆t ) {η&&}t .
β ( ∆t ) β  2β 

When these are put into the incremental equations of motion, Eq. 11.26 can be rewritten in the
following form:

Eq. 11.30 [ K$ ] {∆η} = {∆F$ }

where the pseudo-stiffness matirx is defined as before using Eq. 11.24 and the incremental
pseudo-force vector is given by

 
[∆F$ ] = [ ∆F ] + [m]  β (1∆t ) {η& } t
1
2β
{&&η} t  +

Eq. 11.31 .
γ  γ  
[c]t  β {η& }t +  2 β − 1 ( ∆t ) {&&η}t 
  N  

The procedure for obtaining an estimate for the motion vectors at time t + ∆t is

summarized below:

· the pseudo-static stiffness and incremental force vectors are determined using Eq. 11.31 and
5.9, respectively, based on the stiffness and damping matrices at the starting time t,
· the incremental displacement vector {∆η} is calculated from Eq. 11.30,
 CE 804 - Structural Dynamics Page 11-10

· the incremental acceleration and velocity vectors, {∆&&

η} and {∆η
& } , are determined from Eq.
11.29, and finally
· the motion vectors {η} t + ∆t , {η
& } t + ∆t , and {η
&&} t + ∆t are calculated from Eq. 11.27.

· Iterative Solution for Nonlinear Response: Because of the change in stiffness and damping
matrices during ∆t due to nonlinearities in the system, a state of dynamic equilibrium will
generally not be achieved at time t + ∆t using the motion vector estimates determined above.
The unbalanced forces existing within the system at the end of the ith iteration for time t + ∆t are
determined by the expression

Eq. 11.32 {δF }it + ∆t = {F }t + ∆t − [m] {η&&}it + ∆t − [c]it + ∆t {η& }it + ∆t − {FR }it + ∆t
where [ c]t + ∆t is the damping matrix based on the ith estimate of the displaced position of the mast
i

at time t + ∆t . The vector { FR } t + ∆t represents the elastic restoring forces corresponding to the
i

system degrees of freedom at the end of the ith iteration. The incremental motion vectors at time
t + ∆t required to achieve equilibrium may be written as

Eq. 11.33a {δη}i +1 = {η}it++1∆t − {η}it + ∆t

Eq. 11.33b {δη& } i +1 = {η& }it++1∆t − {η& }it + ∆t

Eq. 11.33c {δη&&} i +1 = {&&η}it++1∆t − {&&η}it + ∆t .

Substituting Eqs. 11.21a and 11.21b into Eqs. 11.33c and 11.33b, respectively, the incremental
acceleration and velocity vectors can be reduced to the following form:

1
Eq. 11.34a {δη&&} i +1 = {δη}i +1
β( ∆t )
2

γ
Eq. 11.34b {δη& } i +1 = {δη}i +1
β( ∆t )

The incremental equations of motion for the (i+1)th iteration are given by

Eq. 11.35a [m] {δη&&}i +1 + [c]it + ∆t {δη& }i +1 + [ k ]it + ∆t {δη}i +1 = {δF }i

or
 CE 804 - Structural Dynamics Page 11-11

[ K$ ]
j +1
Eq. 11.35b {δη} j +1 = {δF} j
where the updated pseudo- stiffness matrix is

[ K$ ] γ
j +1
= [ k ] t + ∆t + 2 [ ] [ c]t + ∆t .
1
m +
j j
Eq. 11.36
β ( ∆t ) β ( ∆t )

For each iteration after the initial estimate, therefore, the procedure is as follows:

· the updated pseudo- stiffness matrix is formed using Eq. 11.36,

the incremental displacement vector {δη}
j +1
· is calculated from Eq. 11.35,

· the incremental acceleration and velocity vectors, {δη

&&} i +1 and {δη
& } i +1 , are determined

from Eq. 11.34,

· the updated estimates for the motion vectors {η} it++1∆t , {η& } it++1∆t , and {&&η} it++1∆t are

calculated from Eq. 11.33, and

the new unbalanced force vector {δF }
i +1
· is determined from Eq. 11.32 based on the

current estimates for the motion vectors at time t + ∆t .

The iterative process is terminated when {δF } and {δη} fall within prescribed tolerances.
i +1 j +1