14.

CALCULATION OF STIFFNESS AND

MASS ORTHOGONAL VECTORS
LDR Vectors are Always More Accurate than Using the
Exact Eigenvectors in a Mode Superposition Analysis

14.1 INTRODUCTION

The major reason to calculate mode shapes (or eigenvectors and eigenvalues) is
that they are used to uncouple the dynamic equilibrium equations for mode
superposition and/or response spectra analyses. The main purpose of a dynamic
response analysis of a structure is to accurately estimate displacements and
member forces in the real structure. In general, there is no direct relationship
between the accuracy of the eigenvalues and eigenvectors and the accuracy of
node point displacements and member forces.

In the early days of earthquake engineering, the Rayleigh-Ritz method of

dynamic analysis was used extensively to calculate approximate solutions. With
the development of high-speed computers, the use of exact eigenvectors replaced
the use of Ritz vectors as the basis for seismic analysis. It will be illustrated in
this book that Load-Dependent Ritz, LDR, vectors can be used for the dynamic
analysis of both linear and nonlinear structures. The new modified Ritz method
produces more accurate results, with less computational effort, than the use of
exact eigenvectors.
14-2 DYNAMIC ANALYSIS OF STRUCTURES

There are several different numerical methods available for the evaluation of the
eigenvalue problem. However, for large structural systems, only a few methods
have proven to be both accurate and robust.

14.2 DETERMINATE SEARCH METHOD

The equilibrium equation, which governs the undamped free vibration of a

typical mode, is given by:

[K - i2 M]v i = 0 or Ki v i = 0 (14.1)

Equation 14.1 can be solved directly for the natural frequencies of the structure
by assuming values for i and factoring the following equation:

T
Ki = Li Di Li (14.2)

From Appendix C the determinant of the factored matrix is defined by:

Det( i ) = D11 D22 - - - - DNN (14.3)

It is possible, by repeated factorization, to develop a plot of the determinant vs.

, as shown in Figure 14.1. This classical method for evaluating the natural
frequencies of a structure is called the determinant search method [1]. It should
be noted that for matrices with small bandwidths the numerical effort to factor
the matrices is very small. For this class of problem the determinant search
method, along with inverse iteration, is an effective method of evaluating the
undamped frequencies and mode shapes for small structural systems. However,
because of the increase in computer speeds, small problems can be solved by any
method in a few seconds. Therefore, the determinant search method is no longer
used in modern dynamic analysis programs.
EIGEN AND RITZ VECTOR EVALUATION 14-3

Det ( )
All Terms In D Positive

Two Neg. Dii

Six Neg. Dii

1 2 3 4 , 5 6

Five Neg. Dii

One Neg. Dii Three Neg. Dii

Figure 14.1 Determinant vs. Frequency for Typical System

14.3 STURM SEQUENCE CHECK

Figure 14.1 illustrates a very important property of the sequence of diagonal

terms of the factored matrix. One notes that for a specified value of i , one can
count the number of negative terms in the diagonal matrix and it is always equal
to the number of frequencies below that value. Therefore, it can be used to check
a method of solution that fails to calculate all frequencies below a specified
value. Also, another important application of the Sturm Sequence Technique is to
evaluate the number of frequencies within a frequency range. It is only necessary
to factor the matrix at both the maximum and minimum frequency points, and the
difference in the number of negative diagonal terms is equal to the number of
frequencies in the range. This numerical technique is useful in machine vibration
problems.

14.4 INVERSE ITERATION

Equation (14.1) can be written in an iterative solution form as:

K Vn(i) = n(i-1) M V(in -1) or LDLT V(i)

n =R
(i)
(14.4)
14-4 DYNAMIC ANALYSIS OF STRUCTURES

The computational steps required for the solution of one eigenvalue and eigenvector
can be summarized as follows:

1. Factor stiffness matrix into triangularized LD LT form during static load

solution phase.

2. For the first iteration, assume R (1) to be a vector of random numbers and
(1)
solve for initial vector V n .

3. Iterate with i = 1, 2 . . .

a. Normalize vector so that VT(i) (i)

n M Vn = I

T(i) (i)
b. Estimate eigenvalue (i)
n = Vn R

c. Check (i)
n for convergence - if converged, terminate

d. i = i + 1 and calculate R(i) = (in-1) MV (i-1)

e. Solve for new vector LD LT V(i)

n =R
(i)

f. Repeat Step 3

It can easily be shown that this method will converge to the smallest unique
eigenvalue.

14.5 GRAM-SCHMIDT ORTHOGONALIZATION

Additional eigenvectors can be calculated using the inverse iteration method if,
after each iteration cycle, the iteration vector is made orthogonal to all previously
calculated vectors. To illustrate the method, let us assume that we have an
approximate vector V that needs to be made orthogonal to the previously
calculated vector Vn . Or, the new vector can be calculated from:

V = V - Vn (14.5)

Multiplying Equation (14.3) by VTn M , we obtain:

EIGEN AND RITZ VECTOR EVALUATION 14-5

VTn MV = VTn MV - VTn M Vn = 0 (14.6)

Therefore, the orthogonality requirement is satisfied if:

VTn MV
= T
= VTn MV (14.7)
M
Vn Vn

If the orthogonalization step is inserted after Step 3.e in the inverse iteration
method, additional eigenvalues and vectors can be calculated.

14.6 BLOCK SUBSPACE ITERATION

Inverse iteration with one vector may not converge if eigenvalues are identical
and the eigenvectors are not unique. This case exists for many real three-
dimensional structures, such as buildings with equal stiffness and mass in the
principle directions. This problem can be avoided by iterating with a block of
orthogonal vectors [2]. The block subspace iteration algorithm is summarized in
Table 14.1 and is the method used in the modern versions of the SAP program.

Experience has indicated that the subspace block size b should be set equal to
the square root of the average bandwidth of the stiffness matrix, but, not less than
six. The block subspace iteration algorithm is relatively slow; however, it is very
accurate and robust. In general, after a vector is added to a block, it requires five
to ten forward reductions and back-substitutions before the iteration vector
converges to the exact eigenvector.
14-6 DYNAMIC ANALYSIS OF STRUCTURES

Table 14.1 Subspace Algorithm for the Generation of Eigenvectors

I. INITIAL CALCULATIONS

A. Triangularize Stiffness Matrix.

(0)
B. Use random numbers to form a block of b vectors V .

II. GENERATE L EIGENVECTORS BY ITERATION i = 1,2...

A. Solve for block of vectors, X(i) in, K X(i) = M V(i-1) .

(i)
B. Make block of vectors, X(i) , stiffness and mass orthogonal, V . Order
eigenvalues and corresponding vectors in ascending order.
(i)
C. Use Gram-Schmidt method to make V orthogonal to all previously
calculated vectors and normalized so that VT(i) M V(i) = I .

D. Perform the following checks and operations:

1. If first vector in block is not converged, go to Step A with i = i + 1 .

2. Save Vector n on Disk.

3. If n equals L , terminate iteration.

4. Compact block of vectors.

5. Add random number vector to last column of block.

Return to Step D.1 with n = n + 1

14.7 SOLUTION OF SINGULAR SYSTEMS

For a few types of structures, such as aerospace vehicles, it is not possible to use
inverse or subspace iteration directly to solve for mode shapes and frequencies.
This is because there is a minimum of six rigid-body modes with zero
frequencies and the stiffness matrix is singular and cannot be triangularized. To
EIGEN AND RITZ VECTOR EVALUATION 14-7

solve this problem, it is only necessary to introduce the following eigenvalue

shift, or change of variable:

n = n (14.8)

Hence, the iterative eigenvalue problem can be written as:

K Vn(i) = n(i-1) M V(in-1) or LDLT V(i)

n =R
(i)
(14.9)

The shifted stiffness matrix is now non-singular and is defined by:

K = K + M (14.10)

The eigenvectors are not modified by the arbitrary shift . The correct
eigenvalues are calculated from Equation (14.8).

14.8 GENERATION OF LOAD-DEPENDENT RITZ VECTORS

The numerical effort required to calculate the exact eigen solution can be
enormous for a structural system if a large number of modes are required.
However, many engineers believe that this computational effort is justifiable if
accurate results are to be obtained. One of the purposes of this section is to
clearly illustrate that this assumption is not true for the dynamic response
analyses of all structural systems.

It is possible to use the exact free-vibration mode shapes to reduce the size of
both linear and nonlinear problems. However, this is not the best approach for the
following reasons:

1. For large structural systems, the solution of the eigenvalue problem for the
free-vibration mode shapes and frequencies can require a significant amount
of computational effort.

2. In the calculation of the free-vibration mode shapes, the spatial distribution

of the loading is completely disregarded. Therefore, many of the mode
shapes that are calculated are orthogonal to the loading and do not participate
in the dynamic response.
14-8 DYNAMIC ANALYSIS OF STRUCTURES

3. If dynamic loads are applied at massless degrees-of-freedom, the use of all

the exact mode shapes in a mode superposition analysis will not converge to
the exact solution. In addition, displacements and stresses near the
application of the loads can be in significant error. Therefore, there is no
need to apply the static correction method as would be required if exact
eigenvectors are used for such problems.

4. It is possible to calculate a set of stiffness and mass orthogonal Ritz vectors,

with a minimum of computational effort, which will converge to the exact
solution for any spatial distribution of loading [2].

It can be demonstrated that a dynamic analysis based on a unique set of Load

Dependent Vectors yields a more accurate result than the use of the same number
of exact mode shapes. The efficiency of this technique has been illustrated by
solving many problems in structural response and in wave propagation types of
problems [4]. Several different algorithms for the generation of Load Dependent
Ritz Vectors have been published since the method was first introduced in 1982
[3]. Therefore, it is necessary to present in Table 14.2 the latest version of the
method for multiple load conditions.

Table 14.2 Algorithm for Generation of Load Dependent Ritz Vectors

I. INITIAL CALCULATIONS

A. Triangularize Stiffness Matrix K = LT DL .

B. Solve for block of b static displacement vectors us resulting from

spatial load patterns F ; or, K us = F .

C. Make block of vectors us , stiffness and mass orthogonal, V1 .

II. GENERATE BLOCKS OF RITZ VECTORS i = 2,....N

A. Solve for block of vectors, Xi , K Xi = M Vi-1 .

B. Make block of vectors, Xi ,stiffness and mass orthogonal, Vi .

EIGEN AND RITZ VECTOR EVALUATION 14-9

Table 14.2 Algorithm for Generation of Load Dependent Ritz Vectors

C. Use Modified Gram-Schmidt method (two times) to make Vi

orthogonal to all previously calculated vectors and normalized so that
VTi M Vi = I .

III. MAKE VECTORS STIFFNESS ORTHOGONAL

A. Solve Nb by Nb eigenvalue problem [K - 2 I]Z = 0 where

K = V KV .
T

B. Calculate stiffness orthogonal Ritz vectors, = VZ .

14.9 A PHYSICAL EXPLANATION OF THE LDR ALGORITHM

The physical foundation for the method is the recognition that the dynamic
response of a structure will be a function of the spatial load distribution. The
undamped, dynamic equilibrium equations of an elastic structure can be written
in the following form:

Mu
&&(t) + Ku(t) = R(t) (14.11)

In the case of earthquake or wind, the time-dependent loading acting on the

structure, R(t) , Equation (13.1), can be written as:

J
R (t ) = f j g(t ) j = F G (t ) (14.12)
j =1

Note that the independent load patterns F are not a function of time. For constant
earthquake ground motions at the base of the structure three independent load
patterns are possible. These load patterns are a function of the directional mass
distribution of the structure. In case of wind loading, the downwind mean wind
pressure is one of those vectors. The time functions G(t) can always be
expanded into a Fourier series of sine and cosine functions. Hence, neglecting
14-10 DYNAMIC ANALYSIS OF STRUCTURES

damping, a typical dynamic equilibrium equation to be solved is of the following

form:

&&(t) + Ku(t) = F sin t

Mu (14.13)

Therefore, the exact dynamic response for a typical loading frequency is of

the following form:

Ku = F + 2 Mu (14.14)

This equation cannot be solved directly because of the unknown frequency of the
loading. However, a series of stiffness and mass orthogonal vectors can be
calculated that will satisfy this equation using a perturbation algorithm. The first
block of vectors is calculated by neglecting the mass and solving for the static
response of the structure. Or:
Ku0 = F (14.15)

From Equation (14.14) it is apparent that the distribution of the error in the
solution, due to neglecting the inertia forces, can be approximated by:

F1 Mu0 (14.16)

Therefore, an additional block of displacement error, or correction, vectors can

be calculated from:
Ku1 = F1 (14.17)

In calculating u 1 the additional inertia forces are neglected. Hence, in continuing

this thought process, it is apparent the following recurrence equation exists:

Ku i = Mui 1 (14.18)

A large number of blocks of vectors can be generated by Equation (14.18).

However, to avoid numerical problems, the vectors must be stiffness and mass
orthogonal after each step. In addition, care should be taken to make sure that all
vectors are linearly independent. The complete numerical algorithm is
summarized in Table 14.2. After careful examination of the LDR vectors, one
can conclude that dynamic analysis is a simple extension of static analysis
EIGEN AND RITZ VECTOR EVALUATION 14-11

because the first block of vectors is the static response from all load patterns
acting on the structure. For the case where loads are applied at only the mass
degrees-of-freedom, the LDR vectors are always a linear combination of the
exact eigenvectors.

It is of interest to note that the recursive equation, used to generate the LDR
vectors, is similar to the Lanczos algorithm for calculating exact eigenvalues and
vectors, except that the starting vectors are the static displacements caused by the
spatial load distributions. Also, there is no iteration involved in the generation of
Load Dependent Ritz vectors.

14.10 COMPARISON OF SOLUTIONS USING EIGEN AND RITZ

VECTORS

The fixed-end beam shown in Figure 14.1 is subjected to a point load at the
center of the beam. The load varies in time as a constant unit step function.

100

10 @ 12 = 240

Modulus of Elasticity = 30,000,000

Moment of Inertia = 100
Mass per Unit Length = 0.1
Damping Ratio = 0.01
All units in Pounds and Inches

Figure 14.1 Dimensions, Stiffness and Mass for Beam Structure

The damping ratio for each mode was set at one percent and the maximum
displacement and moment occur at 0.046 second, as shown in Table 14.3.

The results clearly indicate the advantages of using load-dependent vectors. One
notes that the free-vibration modes 2, 4, 6 and 8 are not excited by the loading
because they are nonsymmetrical. However, the load dependent algorithm
14-12 DYNAMIC ANALYSIS OF STRUCTURES

generates only the symmetrical modes. In fact, the algorithm will fail for this
case, if more than five vectors are requested.

Table 14.3 Results from Dynamic Analyses of Beam Structure

Number of Free-Vibration Mode Shapes Load-Dependent Ritz Vectors

Vectors Displacement Moment Displacement Moment
0.004572 4178 0.004726 5907
1
(-2.41) (-22.8) (+0.88) (+9.2)
0.004572 4178 0.004591 5563
2
(-2.41) (-22.8) (-2.00) (+2.8)
0.004664 4946 0.004689 5603
3
(-0.46) (-8.5) (+0.08) (+3.5)
0.004664 4946 0.004688 5507
4
(-0.46) (-8.5) (+0.06) (+1.8)
0.004681 5188 0.004685 5411
5
(-0.08) (-4.1) (0.00) (0.0)
0.004683 5304
7
(-0.04) (-2.0)
0.004685 5411
9
(0.00) (0.0)

Note: Numbers is parentheses are percentage errors.

Both methods give good results for the maximum displacement. The results for
maximum moment, however, indicate that the load-dependent vectors give
significantly better results and converge from above the exact solution. It is clear
that free-vibration mode shapes are not necessarily the best vectors to be used in
mode-superposition dynamic response analysis. Not only is the calculation of the
exact free-vibration mode shapes computationally expensive, it requires more
vectors, which increases the number of modal equations to be integrated and
stored within the computer.

14.11 CORRECTION FOR HIGHER MODE TRUNCATION

In the analysis of many types of structures, the response of higher modes can be
significant. In the use of exact eigenvectors for mode superposition or response
EIGEN AND RITZ VECTOR EVALUATION 14-13

spectra analyses, approximate methods of analysis have been developed to

improve the results. The purpose of those approximate methods is to account for
missing mass or to add static response associated with higher mode
truncation. Those methods are used to reduce the number of exact eigenvectors
to be calculated, which reduces computation time and computer storage
requirements.

The use of Load Dependent Ritz, LDR, vectors, on the other hand, does not
require the use of those approximate methods because the static response is
included in the initial set of vectors. This is illustrated by the time history
analysis of a simple cantilever structure subjected to earthquake motions shown
in Figure 14.2. This is a model of a light-weight superstructure built on a massive
foundation supported on stiff piles that are modeled using a spring.

C o m p u te r M o d e l

Figure 14.2 Cantilever Structure on Massive Stiff Foundation

Only eight eigen or Ritz vectors can be used because the model has only eight
masses. The computed periods, using the exact eigen or Ritz method, are
summarized in Table 14.4. It is apparent that the eighth mode is associated with
the vibration of the foundation mass and the period is very short: 0.00517
seconds.
14-14 DYNAMIC ANALYSIS OF STRUCTURES

Table 14.4 Periods and Mass Participation Factors

PERIOD MASS PARTICIPATION

MODE NUMBER (Seconds) (Percentage)
1 1.27321 11.706
2 0.43128 01.660
3 0.24205 00.613
4 0.16018 00.310
5 0.11899 00.208
6 0.09506 00.100
7 0.07951 00.046
8 0.00517 85.375

The maximum foundation force using different numbers of eigen and LDR
vectors is summarized in Table 14.5. In addition, the total mass participation
associated with each analysis is shown. The integration time step is the same as
the earthquake motion input; therefore, no errors are introduced other than those
resulting from mode truncation. Five percent damping is used in all cases.

Table 14.5 Foundation Forces and Total Mass Participation

NUMBER FOUNDATION FORCE MASS PARTICIPATION

OF (Kips) (Total Percentage)
VECTORS EIGEN RITZ EIGEN RITZ
8 1,635 1,635 100.0 100.0
7 260 1,636 14.6 83.3
5 259 1,671 14.5 16.2
3 258 1,756 14.0 14.5
2 257 3,188 13.4 13.9

The solution for eight eigen or LDR vectors produces the exact solution for the
foundation force and 100 percent of the participating mass. For seven
eigenvectors, the solution for the foundation force is only 16 percent of the exact
EIGEN AND RITZ VECTOR EVALUATION 14-15

valuea significant error; whereas, the LDR solution is almost identical to the
exact foundation force. It is of interest to note that the LDR method
overestimates the force as the number of vectors is reduceda conservative
engineering result.

Also, it is apparent that the mass participation factors associated with the LDR
solutions are not an accurate estimate the error in the foundation force. In this
case, 90 percent mass participation is not a requirement if LDR vectors are used.
If only five LDR vectors are used, the total mass participation factor is only 16.2
percent; however, the foundation force is over-estimated by 2.2 percent.

14.12 VERTICAL DIRECTION SEISMIC RESPONSE

Structural engineers are required for certain types of structures, to calculate the
vertical dynamic response. During the past several years, many engineers have
told me that it was necessary to calculate several hundred mode shapes for a large
structure to obtain the 90 percent mass participation in the vertical direction. In
all cases, the "exact" free vibration frequencies and mode shapes were used in the
analysis.

To illustrate this problem and to propose a solution, a vertical dynamic analysis is

conducted of the two dimensional frame shown in Figure 14.3. The mass is
lumped at the 35 locations shown; therefore, the system has 70 possible mode
shapes.

Using the exact eigenvalue solution for frequencies and mode shapes, the mass
participation percentages are summarized in Table 14.6.

One notes that the lateral and vertical modes are uncoupled for this very simple
structure. Only two of the first ten modes are in the vertical direction. Hence, the
total vertical mass participation is only 63.3 percent.
14-16 DYNAMIC ANALYSIS OF STRUCTURES

Figure 14.3 Frame Structure Subjected to Vertical Earthquake Motions

Table 14.6 Mass Participation Percentage Factors for Exact Eigenvalues

LATERAL MASS VERTICAL MASS

PERIOD PARTICIPATION PARTICIPATION
MODE
(Seconds)
EACH MODE TOTAL EACH MODE TOTAL
1 1.273 79.957 79.957 0 0
2 0.421 11.336 91.295 0 0
3 0.242 4.172 95.467 0 0
4 0.162 1.436 96.903 0 0
5 0.158 0.650 97.554 0 0
6 0.148 0 97.554 60.551 60.551
7 0.141 0.031 97.584 0 60.551
8 0.137 0.015 97.584 0 60.551
9 0.129 0.037 97.639 0 60.551
10 0.127 0 97.639 2.775 63.326
EIGEN AND RITZ VECTOR EVALUATION 14-17

The first 10 Load Dependent Ritz vectors are calculated and the mass
participation percentages are summarized in Table 14.7. The two starting LDR
vectors were generated using static loading proportional to the lateral and vertical
mass distributions.

Table 14.7 Mass Participation Percentage Factors Using LDR Vectors

LATERAL MASS VERTICAL MASS

PERIOD PARTICIPATION PARTICIPATION
MODE
(Seconds)
EACH MODE TOTAL EACH MODE TOTAL
1 1.273 79.957 79.957 0 0
2 0.421 11.336 91.295 0 0
3 0.242 4.176 95.471 0 0
4 0.158 2.388 97.859 0 0
5 0.149 0 97.859 60.567 60.567
6 0.123 0 97.859 4.971 65.538
7 0.104 2.102 99.961 0 65.538
8 0.103 0 99.961 13.243 78.781
9 0.064 0 99.961 9.696 88.477
10 0.041 0 99.961 8.463 96.940

The ten vectors produced by the LDR method more than satisfy the 90 percent
code requirement. It would require the calculation of 34 eigenvectors for the
exact eigenvalue approach to obtain the same mass participation percentage. This
is just one additional example of why use of the LDR method is superior to the
use of the exact eigenvectors for seismic loading.

The reason for the impressive accuracy of the LDR method compared to the
exact eigenvector method is that only the mode shapes that are excited by the
seismic loading are calculated.
14-18 DYNAMIC ANALYSIS OF STRUCTURES

14.13 SUMMARY

There are three different mathematical methods for the numerical solution of the
eigenvalue problem. They all have advantages for certain types of problems.

First, the determinant search method, which is related to finding the roots of a
polynomial, is a fundamental traditional method. It is not efficient for large
structural problems. The Sturm sequence property of the diagonal elements of the
factored matrix can be used to determine the number of frequencies of vibration
within a specified range.

Second, the inverse and subspace iteration methods are subsets of a large number of
power methods. The Stodola method is a power method. However, the use of a
sweeping matrix to obtain higher modes is not practical because it eliminates the
sparseness of the matrices. Gram-Schmidt orthogonalization is the most effective
method to force iteration vectors to converge to higher modes.

Third, transformation methods are very effective for the calculation of all eigenvalues
and eigenvectors of small dense matrices. Jacobi, Givens, Householder, Wilkinson
and Rutishauser are all well-known transformation methods. The author prefers to use
a modern version of the Jacobi method in the ETABS and SAP programs. It is not the
fastest; however, we have found it to be accurate and robust. Because it is only used
for problems equal to the size of the subspace, the computational time for this phase
of the solution is very small compared to the time required to form the subspace
eigenvalue problem. The derivation of the Jacobi method is given in Appendix D.

The use of Load Dependent Ritz vectors is the most efficient approach to solve for
accurate node displacements and member forces within structures subjected to
dynamic loads. The lower frequencies obtained from a Ritz vector analysis are always
very close to the exact free vibration frequencies. If frequencies and mode shapes are
missed, it is because the dynamic loading does not excite them; therefore, they are of
no practical value. Another major advantage of using LDR vectors is that it is not
necessary to be concerned about errors introduced by higher mode truncation of a set
of exact eigenvectors.

All LDR mode shapes are linear combinations of the exact eigenvectors;
therefore, the method always converges to the exact solution. Also, the
EIGEN AND RITZ VECTOR EVALUATION 14-19

computational time required to calculate the LDR vectors is significantly less

than the time required to solve for eigenvectors.

14.14 REFERENCES

1. Bathe, K. J., and E. L. Wilson. 1972. "Large Eigenvalue Problems in

Dynamic Analysis," Proceedings, American Society of Civil Engineers,
Journal of the Engineering Mechanics Division, EM6. December. pp.
1471-1485.

2. Wilson, E. L., and T. Itoh. 1983. "An Eigensolution Strategy for Large
Systems," in J. Computers and Structures. Vol. 16, No. 1-4. pp. 259-265.

3. Wilson, E. L., M. Yuan and J. Dickens. 1982. Dynamic Analysis by Direct

Superposition of Ritz Vectors, Earthquake Engineering and Structural
Dynamics. Vol. 10. pp. 813-823.

4. Bayo, E. and E. L. Wilson. 1984. "Use of Ritz Vectors in Wave Propagation

and Foundation Response," Earthquake Engineering and Structural
Dynamics. Vol. 12. pp. 499-505.
14-20 DYNAMIC ANALYSIS OF STRUCTURES