结构动力学及其工程应用

Structural Dynamics and Its Applications

Solving Eigenvalue Problem

阳佳桦，Ph.D.，副教授
广西大学土木建筑工程学院，工程力学研究中心
Email: javayang@gxu.edu.cn
Dr. Jia-Hua Yang, Ph.D., Associate Professor
College of Civil Engineering and Architecture
Scientific Research Center of Engineering Mechanics
Guangxi University
Contents
 Fundamental Mode Analysis

 Power Spectral Density (PSD)

 Fourier Series, Fourier Transform and PSD

2
Fundamental Mode Analysis

Inertial force

displacements resulting from these forces may be calculated by solving

the static deflection problem

Define dynamic matrix

It summarizes the dynamic
properties of the structure
So,
3
To start the iteration, assume a trial displacement vector
Inertial force by the Reasonable guess of
assumed shape the 1-st mode shape

A better approximation of the first mode shape than was the initial
vector using the previous equation

We only care about the “shape”, so the frequency scaling factor can be
dropped from the above equation:

Iteration equation

Then the improved iteration vector is obtained finally by normalizing

this shape, dividing it by an arbitrary reference element of the vector

4
Using the iteration equation repeatedly to get better approximation of
the mode shape
(𝑠) (𝑠−1)
𝐯ത1 = 𝐃𝐯1
Suppose that after s iterations the vector is accurate enough to approximate
the mode shape, i.e., it has converged. The s-th vector is obtained from the (s-
1)-th vector by
The frequency scaling factor is
added back to get the actual
vector for approximation

(𝑠−1)
We think that 𝐯1 is accurate
enough to be mode shape 𝛟1

The natural frequency can be approximated using any component of the

two successive vectors
(𝑠−1)
𝑣𝑘1
𝜔12 ≈ (𝑠)
5
𝑣𝑘1
ҧ
The range of the natural frequency
𝑠−1 𝑠−1
𝑣𝑘1 𝑣𝑘1
𝑠
< 𝜔12 < 𝑠
𝑣𝑘1
ҧ 𝑚𝑖𝑛
𝑣𝑘1
ҧ 𝑚𝑎𝑥

Two frequency approximation methods. It is normalized such

𝑠−1 that the maximum
max 𝐯1 1 component is 1.
𝜔12 ≈ 𝑠
= 𝑠
max 𝐯ത1 max 𝐯ത1
Averaged way using all
components: Pre- 𝑠 𝑇 𝑠−1
𝐯ത1 𝐦𝐯1
multiplying numerator and 𝜔12 ≈
𝑠 𝑇 𝑠
𝑠 𝑇 𝐯ത1 𝐦ത
𝐯1
ത1 𝐦
denominator by 𝐯

6
Clough and Penzien, 2000. Dynamics of Structures. CSI.
 Example

Clough and Penzien, 2000. Dynamics of Structures. CSI.

  Iteration

Note that the factor 1/3,600 has

not been considered in this phase
of the analysis because only the
relative shape is important. The
shapes have been normalized by
dividing by the largest
displacement component.

8
Clough and Penzien, 2000. Dynamics of Structures. CSI.
Using the largest displacement component, the first mode frequency is
found to be

in which it will be noted that the factor 1/3,600 has now been included.
It also is of interest to determine the range of frequencies obtained after
one cycle Range is large, and the frequency approximation is not good.

The averaged formula gives an acceptable approximation

9
Clough and Penzien, 2000. Dynamics of Structures. CSI.
Assumed shape in modal coordinates:

Inertial force corresponding to the assumed shape:

Construct

The deflections derived from these inertial forces are

10
Clough and Penzien, 2000. Dynamics of Structures. CSI.
Re-write

𝐤𝛟𝑛 = 𝜔𝑛2 𝐦𝛟𝑛

𝛟𝑛 = 𝜔𝑛2 𝐤 −1 𝐦𝛟𝑛
𝛟𝑛 = 𝜔𝑛2 𝐃𝛟𝑛

Final step of the 1-st iteration

11
Clough and Penzien, 2000. Dynamics of Structures. CSI.
The second iteration

Prove it in detail.
The s-th iteration:

Will be zero after enough iterations.

Therefore, this iteration converges to the first mode

Clough and Penzien, 2000. Dynamics of Structures. CSI.

13
Analysis of The Second Mode

Purify the first mode Substitute

component from the
assumed shape
Sweeping matrix

14
Clough and Penzien, 2000. Dynamics of Structures. CSI.
 Similar iteration method with the one for the fundamental mode.

Iterate after taking out

the first mode component

Consider purifying and iterating together:

Iteration will continue using the above equation. By using 𝐃2 , it is

equivalent to the first-mode iteration method. Each iteration will take
out the first-mode component by the sweeping matrix 𝐒1 .
Frequency approximation after one iteration:
1
𝐯ത2 denotes the vector without
the frequency scaling factor.

15
Clough and Penzien, 2000. Dynamics of Structures. CSI.
 Continuous iteration
1 2 1 2 1
𝐯 = 𝐃 𝐯
2 2 , 𝑜𝑟 ത
𝐯 = 𝐃 𝐯
2 2
𝜔22 2 2

3 2 4 3
𝐯ത2 = 𝐃2 𝐯2 , 𝐯ത2 = 𝐃2 𝐯2 , ⋯
Because each iteration takes out the first-mode component, this
iteration will automatically converge to the second mode.
 Example. Use the vector iteration method to analyze the second
mode.

Clough and Penzien, 2000. Dynamics of Structures. CSI.

17
Clough and Penzien, 2000. Dynamics of Structures. CSI.
Clough and Penzien, 2000. Dynamics of Structures. CSI.
After 12 iterations, the approximate mode shape is accurate enough

The frequency of the second mode vibration derived from the top story
displacement after the first cycle of iteration is given by

1 2 1
𝜔2 = 1
max 𝐯ത2
Many iterations are needed to get the
accurate frequency and mode shape for
 Comparison the 2nd mode. Even more iterations are
needed for higher modes.
After 4 iterations, the approximation using the top displacement is
𝜔22 = 30.90 𝑟𝑎𝑑/𝑠
The average formula (Rayleigh quotient): 𝜔22 = 31.06 𝑟𝑎𝑑/𝑠
2
Correct:
19
𝜔2 = 31.048 𝑟𝑎𝑑/𝑠
Clough and Penzien, 2000. Dynamics of Structures. CSI.
Analysis of Third and Higher Modes
For the 3-rd mode, take out the first- and second- mode components
(0) (0) (0) (0)
𝐯3 = 𝛟1 𝑌1 + 𝛟2 𝑌2 + 𝛟3 𝑌3 + ⋯
(0) (0) (0)
𝛟1𝑇 𝐦𝐯3 = 𝛟1𝑇 𝐦𝛟1 𝑌1 + 𝛟1𝑇 𝐦𝛟2 𝑌2 + ⋯
𝑇 (0) 𝑇 (0) 𝑇 (0)
𝛟2 𝐦𝐯3 = 𝛟2 𝐦𝛟1 𝑌1 + 𝛟2 𝐦𝛟2 𝑌2 + ⋯

20
Clough and Penzien, 2000. Dynamics of Structures. CSI.
The sweeping matrix that takes out the first- and second- mode
components:

Clough and Penzien, 2000. Dynamics of Structures. CSI.

 Iteration that converges to the 3-rd mode can be established.
Purify the 1-st and 2-nd modes
Iterate equation

 We have the similar procedure for the 4-th mode.

Sweeping matrix for the 4-th mode

Purify modes 1 to 3
Iteration matrix
The matrices suitable for calculating any mode can be obtained easily by
analogy from these

Clough and Penzien, 2000. Dynamics of Structures. CSI.

Clearly the most important limitation of this procedure is that all the
lower mode shapes must be calculated before any given higher mode can
be evaluated. Also, it is essential to evaluate these lower modes with great
precision if the sweeping matrix for the higher modes is to perform
effectively. Generally this process is used directly for the calculation of
no more than four or five modes.

Clough and Penzien, 2000. Dynamics of Structures. CSI.

Analysis of Highest Mode
Based on
Pre-multiply to make use of the stiffness matrix. Compared with
the fundamental mode analysis, the inverse of the stiffness matrix
(flexibility matrix) is used.
𝐦−1 𝐤𝐯ො𝑛 = 𝜔𝑛2 𝐦−1 𝐤 𝐤 −1 𝐦ො 𝐯𝑛

If a trial shape for the highest (N-th) mode of vibration is introduced, the
iteration is

1 2
Notice that the frequency scaling has changed from 2 to 𝜔𝑁
𝜔𝑁
24
Clough and Penzien, 2000. Dynamics of Structures. CSI.
 Frequency approximation using one component of the vector (usually
the maximal component)

 Frequency approximation using the average formula

Compared with the fundamental
mode case, the numerator and
denominator have changed
locations.

Clough and Penzien, 2000. Dynamics of Structures. CSI.

Proof of Convergence to The Highest Mode
Expanding the assumed initial vector and using the iteration equation
gives the vector after the 1-st iteratoin 0
𝐯𝑛
1 1 0 0 0
𝐯ത𝑁 = 2 𝐄 𝛟1 𝑌1 + 𝛟2 𝑌2 + 𝛟3 𝑌3 + ⋯
𝜔𝑁
2
1 𝜔12 0 1 𝜔 2
1 0 1 𝜔 𝑁 0
= 2 2 𝐄𝛟1 𝑌1 + 2 2 𝐄𝛟2 𝑌2 + ⋯ + 2 2 𝐄𝛟𝑁 𝑌𝑁
𝜔1 𝜔𝑁 𝜔2 𝜔𝑁 𝜔𝑁 𝜔𝑁
For each term on the right-hand side, the eigenvalue equation gives
1 𝐄 1
2 −1
𝐤𝛟𝑛 = 𝜔𝑛 𝐦𝛟𝑛 ⟹ 2 𝐦 𝐤𝛟𝑛 = 𝛟𝑛 ⟹ 𝛟𝑛 = 2 𝐄𝛟𝑛
𝜔𝑛 𝜔𝑛
So the first iteration becomes
𝑁
1 𝜔12 0 𝜔12 0 0 𝜔𝑛2 0
𝐯ത𝑁 = 2 𝛟1 𝑌1 + 2 𝛟2 𝑌2 + ⋯ + 𝛟𝑁 𝑌𝑁 = ෍ 2 𝛟2 𝑌2
26 𝜔𝑁 𝜔𝑁 𝜔𝑁
Clough and Penzien, 2000. Dynamics of Structures. CSI. 𝑛=1
The iteration continues
𝑁
1 1 𝜔 2
2 1 𝑛 0
𝐯ത𝑁 = 2 𝐄ത
𝐯𝑁 = 2 𝐄 ෍ 2 𝛟𝑛 𝑌𝑛
𝜔𝑁 𝜔𝑁 𝜔𝑁
𝑛=1

The previous derivation is repeated, and finally we can construct

Only the term of the highest mode in the above series remains, so this
iteration will converge to the highest mode. Compare this proof with the
one for fundamental mode analysis.

Clough and Penzien, 2000. Dynamics of Structures. CSI.

 Example

The stiffness form of the dynamic

matrix is

It is evident that this iteration

process converges toward the highest
mode shape much more slowly than
the convergence toward the lowest
mode.

Clough and Penzien, 2000. Dynamics of Structures. CSI.

Inverse Iteration - Recommended
Noting that the effect of the frequency will be removed subsequently by
the normalization step, in this formulation the frequency is assumed to be
unity and the resulting inertial forces are denoted by
No frequency scaling factor

The improved displacement vector

(1) Define (2) Solve

Solve
(3) Normalize
29
Clough and Penzien, 2000. Dynamics of Structures. CSI.
It is important to note that the narrow-banded character of the stiffness
matrix k is retained in the triangular matrices L and U, consequently the
efficiency of this inverse displacement analysis is greatly enhanced
relative to the flexibility matrix formulation used with direct iteration.

Clough and Penzien, 2000. Dynamics of Structures. CSI.

Inverse Iteration With Shifts

The essential concept of shifting is the representation of each

eigenvalue as the sum of a shift plus a residual

Consider the entire diagonal matrix of eigenvalues

The shift can be visualized as a displacement of the origin in a plot of

the eigenvalues. Its effect is to transform the eigenvalue problem to the
analysis of the residuals rather than the actual eigenvalues

31
Clough and Penzien, 2000. Dynamics of Structures. CSI.
It is apparent that the above eigenvalue problem is entirely equivalent to
the original one and that the shifted matrix has the same eigenvectors as
E.

Clough and Penzien, 2000. Dynamics of Structures. CSI.

The solution of this new eigenproblem may be carried out by inverse
iteration, because it is in the stiffness form. The result of the first cycle
of iteration can be expressed as follows

After s iteration cycles, This can be proved using

the similar procedures
shown before.

where 𝛿𝑘 represents the smallest residual eigenvalue, that is,

Clough and Penzien, 2000. Dynamics of Structures. CSI.

This iteration formula will converge to the smallest residual
eigenvalue, working in a similar mechanism that the iteration
converges to the first mode.

𝛿𝑘
Will be zero if s is large enough, because 𝛿𝑘 is the smallest and <1
𝛿𝑛

Thus it is evident that the two summations will become negligibly

small after a sufficient number of iteration cycles, so the computed
mode shape converges to
This analysis therefore shows that the process of inverse iteration with
eigenvalue shift converges to the mode shape for which the eigenvalue is
closest to the shift position

𝛿2 is the smallest.

Will converge to this eigenvalue that is closest to the shifted origin.

Clough and Penzien, 2000. Dynamics of Structures. CSI.

Similarly, the residual eigenvalue can be estimated as

The actual eigenvalue is thus

Because the speed of convergence can be accelerated by shifting to a

value very close to the eigenvalue that is sought, it is good practice to
shift at intervals during iteration, i.e., change the shifted point for every
𝑁𝑠 iteration cycles using the Rayleigh quotient:

36
Clough and Penzien, 2000. Dynamics of Structures. CSI.
 The efficient inverse iteration procedures are also applicable for
inverse iteration with shifts.

Introducing shifts

Assume an initial trial vector for mode 𝑘, to obtain a trial inertial

load vector . The first iteration becomes

This gives the improved displacements . In order to take advantage of

the narrow banding of , it is reduced by Gauss elimination to upper and
lower triangular form as before

Efficient solution as follows

37
Clough and Penzien, 2000. Dynamics of Structures. CSI.
 (1) Define (2) Solve

Solve
(3) Normalize

Very often, the requirement is to determine the lowest few eigenpairs.

In such cases, inverse iteration with shift must be combined with
iteration without shift.
The first few, say five, frequencies and mode shapes are calculated by
the standard inverse iteration technique along with Gram–Schmidt
orthogonalization. The same process is then employed for the sixth
eigenvalue, but instead of completing the iteration, it is discontinued
after a few passes, at which time an estimate of the sixth frequency is
obtained by use of the Rayleigh quotient This estimated frequency is
now used as the new shift point. The next five frequencies and mode
shapes are then calculated using inverse iteration with the shifted
origin.
38
Clough and Penzien, 2000. Dynamics of Structures. CSI.
Subspace Method
The subspace method developed by Bathe KJ is an efficient method for
the eigensolution of large systems when only the lower modes are of
interest. The number of trial vectors m is usually much less than N, the
order of the matrices involved, but is larger than p, the number of
modes to be determined. In practice, the method is found to be most
efficient when m is chosen to be the smaller of 2p and p + 8.

Solving the first updated shape vector with inertial load constructed
using the initial shape matrix 𝐕0 :
𝐊𝐕෩1 = 𝐌𝐕0 ⟹ 𝐕 ෩1 = 𝐊 −1 𝐌𝐕0

39
Two requirements: 1. 𝐕 ෩1 must be orthogonalized so that they will
converge to different eigenvectors or mode shapes and not all to the
lowest one; 2. The eigenvectors should be normalized in some way so
that the numbers remain within reasonable bounds.
These can be done by solving a new eigenvalue problem

is mass (not the original mass matrix) normalized so that

The improved trial vectors for the next iteration

40
 The improved trial vectors are orthogonal with respect to the
stiffness and mass matrices of the original system and will therefore
converge to different mode shapes. Proof:

41
 Summary of Subspace Method.
 Given , find the next iteration vector by solving the
following equation (note the inertial load with the shape):

 Find the projection of the stiffness and mass matrices in the subspace
of . The size of this subspace is 𝑚 × 𝑚, which is much smaller
than the original eigenvector space.

 Solve the following eigenvalue problem in the reduced space

 Normalize the eigenvectors so that

42
 Find an improved estimate of the eigenvectors

When k becomes large the vectors in 𝐕𝑘+1 will converge to the first m
eigenvectors and the diagonal elements of 𝚪𝑘+1 will converge to the
first m eigenvalues.

43
Lanczos Method For Standard Eigenvalue Problem
The Lanczos Method transforms a symmetric matrix 𝐀 to a tridiagonal
matrix 𝐓 using a series of mutually orthogonal vectors contained in the
columns of 𝐗, i.e.,
or
The Lanczos Method is considered to be as effective a tool as the
Subspace Method for the solution of large vibration problems (ANSYS
uses Lanczos Method).
When written in its expanded form, the above transformation becomes
Our purpose is to determine
these coefficients.

44
Or,

The process of finding the Lanczos vectors begins with an arbitrary

selection for the initial vector , and then update it to that is
normalized with unit length (for our generalized eigenvalue problem, it
will be mass normalized), i.e., orthonormal.
Unit length because
of orthonormality

Scaling factor

45
Next, we determine 𝐱 2 and select 𝛼1 . Bear in mind that it should be
orthonormal to 𝐱1 . Multiplication of the first equation by 𝐱1𝑇

If we select

and notice , the first equation becomes

𝛼1 = 𝛼1 + 𝛽2 𝐱1𝑇 𝐱 2 ⟹ 𝐱1𝑇 𝐱 2 = 0
So this choice of 𝛼1 ensures that 𝐱1 and 𝐱 2 are mutually orthogonal.
Using this 𝛼1 , we can get the next shape

Using the normalization by vector length gives

46
Next, we determine 𝐱 3 and select 𝛼2 using the similar procedures.
Multiplication of the second equation by 𝐱 2𝑇

If we select

and notice and , the second equation becomes

𝛼2 = 𝛼2 + 𝛽3 𝐱 2𝑇 𝐱 3 ⟹ 𝐱 2𝑇 𝐱 3 = 0
So this choice of 𝛼2 ensures that 𝐱 2 and 𝐱 3 are mutually orthogonal.
Using this 𝛼2 , we can get the next shape 𝐱෤ 3
𝐱෤ 3 = 𝛽3 𝐱 3 = 𝐀𝐱 2 − 𝛽2 𝐱1 − 𝛼2 𝐱 2
1/2
𝛽3 = 𝐱෤ 3𝑇 𝐱෤ 3
𝐱෤ 3
𝐱3 =
𝛽3
47
 The complete iteration process can be summarized as

1. The two eigenvalue problems

imply that the eigenvalues of
matrix T are the same as those of
matrix A, while the eigenvectors of
the two matrices are related.
2. The smaller size eigenvalue
problem of matrix T can be solved
by QR or LR method.
The iteration begins with 𝑗 = 1 and 𝐱 0 = 𝟎 and it ends when 𝐱 𝑁 has
been obtained.
Pre-multiply 𝐗 𝑇 : 𝐗 𝑇 𝐀𝐪 = 𝜆𝐗 𝑇 𝐪

where 𝐓 = 𝐗 𝑇 𝐀𝐗, 𝐪
෥ = 𝐗𝑇 𝐪
Use

48 Construct the target

tridiagonal matrix T
Theoretically, the Lanczos method vectors are mutually orthogonal to
each other. However, because of round-off errors during computations,
orthogonality relationship breaks down when vectors are sufficiently
separated from each other. For large systems, this source of error
makes Lanczos method unstable. It is therefore necessary to re-
orthogonalize a newly determined Lanczos vector by sweeping off any
contribution from vectors determined previously. Purification is carried
out by using the Gram–Schmidt process. The purified vector:

To verify whether convergence has been achieved for the eigenvalues,

it may be necessary to solve the eigenvalues for 𝐓𝑚−1 as well, and to
compare the magnitude of the desired eigenvalues obtained from 𝐓𝑚−1
and 𝐓𝑚 .

49
Lanczos Method For The Generalized
Eigenvalue Problem
In the case of structural dynamics, the normalization of the vectors 𝐗 in
a lower subspace is chosen to be mass normalization

The generalized eigenvalue problem needs to be transformed to the

standard one using another set of eigenvectors in a lower subspace.

Using another set of

orthonormal vectors to
represent the mode shape.

1 Pre-multiply
𝐊 −1 𝐌𝐗𝜙෨ = 𝐗𝜙෨ (Note that before we multiply 𝐗 𝑇
𝜆 for standard eigenvalue problem)
50
 𝐈
Standard eigenvalue problem

Because ,

𝐓 is a tridiagonal matrix as the standard eigenvalue problem before, so

we have the following

51
To begin the iteration, an arbitrary vector 𝐱෤ 1 is selected and is
normalized with respect to the mass matrix so that .

52
Next, update the shape and choose 𝛼1 . To do this, first update the
shape 𝐱ത 2 using the inertial load based on the normalized 𝐱1 , i.e., solve
the following equation
𝐊ത𝐱 2 = 𝐌𝐱1 → 𝐱ത 2 = 𝐊 −1 𝐌𝐱1
We then choose 𝛼1 = 𝐱1𝑇 𝐌ത𝐱 2 = 𝐱1𝑇 𝐌𝐊 −1 𝐌𝐱1 , this is because we find
that if we pre-multiply 𝐱1𝑇 𝐌 at both sides of the first equation

𝐱1𝑇 𝐌𝐊 −1 𝐌𝐱1 = 𝛼1 𝐱1𝑇 𝐌𝐱1 + 𝛽2 𝐱1𝑇 𝐌𝐱 2

With the above choice =1
we know that this is 𝛼1

Therefore, and this choice of 𝛼1 automatically gives the

orthonormality for 𝐱1 and 𝐱 2 .
After the updating to get 𝐱ത 2 and 𝛼1 , we can update 𝐱ത 2 to 𝐱෤ 2 using the
first equation
53
Then, normalize the vector to finish one cycle of iteration

The complete iteration process is summarized as

54
 Example of Lanczos Method

55
56
57