by the following property, for whose proof we refer to [GL89], pp. 475-477.

Property 5.12 Let A ∈ Rn—n be a symmetric matrix with eigenvalues

ordered as in (5.66) and let u1 , . . . , un be the corresponding orthonormal

eigenvectors. If ·1 , . . . , ·m denote the eigenvalues of Hm , with ·1 ≥ ·2 ≥

. . . ≥ ·m , then

(»1 ’ »n )(tan φ1 )2

»1 ≥ ·1 ≥ »1 ’ ,

(Tm’1 (1 + 2ρ1 ))2

234 5. Approximation of Eigenvalues and Eigenvectors

where cos φ1 = |(q(1) )T u1 |, ρ1 = (»1 ’ »2 )/(»2 ’ »n ) and Tm’1 (x) is the

Chebyshev polynomial of degree m ’ 1 (see Section 10.1.1).

A similar result holds of course for the convergence estimate of the eigen-

values ·m to »n

(»1 ’ »n )(tan φn )2

»n ¤ ·m ¤ »n + ,

(Tm’1 (1 + 2ρn ))2

where ρn = (»n’1 ’ »n )/(»1 ’ »n’1 ) and cos φn = |(q(n) )T un |.

A naive implementation of the Lanczos algorithm can be a¬ected by nu-

merical instability due to propagation of rounding errors. In particular, the

Lanczos vectors will not verify the mutual orthogonality relation, making

the extremal properties (5.67) false. This requires careful programming of

the Lanczos iteration by incorporating suitable reorthogonalization proce-

dures as described in [GL89], Sections 9.2.3-9.2.4.

Despite this limitation, the Lanczos method has two relevant features:

it preserves the sparsity pattern of the matrix (unlike Householder tridiag-

onalization), and such a property makes it quite attractive when dealing

with large size matrices; furthermore, it converges to the extremal eigen-

values of A much more rapidly than the power method does (see [Kan66],

[GL89], p. 477).

The Lanczos method can be generalized to compute the extremal eigen-

values of an unsymmetric matrix along the same lines as in Section 4.5

in the case of the solution of a linear system. Details on the practical im-

plementation of the algorithm and a theoretical convergence analysis can

be found in [LS96] and [Jia95], while some documentation of the latest

software can be found in NETLIB/scalapack/readme.arpack (see also the

MATLAB command eigs).

An implementation of the Lanczos algorithm is provided in Program 45.

The input parameter m is the size of the Krylov subspace in the tridiago-

nalization procedure, while toll is a tolerance monitoring the size of the

increment of the computed eigenvalues between two successive iterations.

The output vectors lmin, lmax and deltaeig contain the sequences of the

approximate extremal eigenvalues and of their increments between succes-

sive iterations. Program 42 is invoked for computing the eigenvalues of the

tridiagonal matrix Hm .

Program 45 - eiglancz : Extremal eigenvalues of a symmetric matrix

function [lmin,lmax,deltaeig,k]=eiglancz(A,m,toll)

n=size(A); V=[0*[1:n]™,[1,0*[1:n-1]]™];

beta(1)=0; normb=1; k=1; deltaeig(1)=1;

while k <= m & normb >= eps & deltaeig(k) < toll

5.12 Applications 235

vk = V(:,k+1); w = A*vk-beta(k)*V(:,k);

alpha(k)= w™*vk; w = w - alpha(k)*vk;

normb = norm(w,2); beta(k+1)=normb;

if normb ˜= 0

V=[V,w/normb];

if k==1

lmin(1)=alpha; lmax(1)=alpha;

k=k+1; deltaeig(k)=1;

else

d=alpha; b=beta(2:length(beta)-1);

[ak,bk,ck,nch,niter]=givsturm(d,b,1,toll);

lmax(k)=(ak(niter)+bk(niter))/2;

[ak,bk,ck,nch,niter]=givsturm(d,b,k,toll);

lmin(k)=(ak(niter)+bk(niter))/2;

deltaeig(k+1)=max(abs(lmin(k)-lmin(k-1)),abs(lmax(k)-lmax(k-1)));

k=k+1;

end

else

disp(™Breakdown™);

d=alpha; b=beta(2:length(beta)-1);

[ak,bk,ck,nch,niter]=givsturm(d,b,1,toll);

lmax(k)=(ak(niter)+bk(niter))/2;

[ak,bk,ck,nch,niter]=givsturm(d,b,k,toll);

lmin(k)=(ak(niter)+bk(niter))/2;

deltaeig(k+1)=max(abs(lmin(k)-lmin(k-1)),abs(lmax(k)-lmax(k-1)));

k=k+1;

end

end

k=k-1;

return

Example 5.18 Consider the eigenvalue problem for the matrix A∈ Rn—n with

n = 100, having diagonal entries equal to 2 and o¬-diagonal entries equal to -1

on the upper and lower tenth diagonal. Program 45, with m=100 and toll=eps,

takes 10 iterations to approximate the extremal eigenvalues of A with an absolute

•

error of the order of the machine precision.

5.12 Applications

A classical problem in engineering is to determine the proper or natural

frequencies of a system (mechanical, structural or electric). Typically, this

leads to solving a matrix eigenvalue problem. Two examples coming from

structural applications are presented in the forthcoming sections where the

buckling problem of a beam and the study of the free vibrations of a bridge

are considered.

236 5. Approximation of Eigenvalues and Eigenvectors

5.12.1 Analysis of the Buckling of a Beam

Consider the homogeneous and thin beam of length L shown in Figure

5.4. The beam is simply supported at the end and is subject to a normal

compression load P at x = L. Denote by y(x) the vertical displacement

of the beam; the structure constraints demand that y(0) = y(L) = 0. Let

L

P

x

y

FIGURE 5.4. A simply supported beam subject to a normal compression load

us consider the problem of the buckling of the beam. This amounts to

determining the critical load Pcr , i.e. the smallest value of P such that an

equilibrium con¬guration of the beam exists which is di¬erent from being

rectilinear. Reaching the condition of critical load is a warning of structure

instability, so that it is quite important to determine its value accurately.

The explicit computation of the critical load can be worked out under

the assumption of small displacements, writing the equilibrium equation

for the structure in its deformed con¬guration (drawn in dashed line in

Figure 5.4)

’E (J(x)y (x)) = Me (x), 0<x<L

(5.68)

y(0) = y(L) = 0,

where E is the constant Young™s modulus of the beam and Me (x) = P y(x)

is the momentum of the load P with respect to a generic point of the beam

of abscissa x. In (5.68) we are assuming that the momentum of inertia

J can be varying along the beam, which indeed happens if the beam has

nonuniform cross-section.

Equation (5.68) expresses the equilibrium between the external momen-

tum Me and the internal momentum Mi = ’E(Jy ) which tends to restore

the rectilinear equilibrium con¬guration of the beam. If the stabilizing re-

action Mi prevails on the unstabilizing action Me , the equilibrium of the

initial rectilinear con¬guration is stable. The critical situation (buckling of

the beam) clearly arises when Mi = Me .

Assume that J is constant and let ±2 = P/(EJ); solving the boundary value

problem (5.68), we get the equation C sin ±L = 0, which admits nontrivial

5.12 Applications 237

solutions ± = (kπ)/L, k = 1, 2, . . . . Taking k = 1 yields the value of the

2

critical load Pcr = π L2 .

EJ

To solve numerically the boundary value problem (5.68) it is conve-

nient to introduce for n ≥ 1, the discretization nodes xj = jh, with

h = L/(n + 1) and j = 1, . . . , n, thus de¬ning the vector of nodal ap-

proximate displacements uj at the internal nodes xj (where u0 = y(0) = 0,

un+1 = y(L) = 0). Then, using the ¬nite di¬erence method (see Section

12.2), the calculation of the critical load amounts to determining the small-

est eigenvalue of the tridiagonal symmetric and positive de¬nite matrix

A = tridiagn (’1, 2, ’1) ∈ Rn—n .

It can indeed be checked that the ¬nite di¬erence discretization of prob-

lem (5.68) by centered di¬erences leads to the following matrix eigenvalue

problem

Au = ±2 h2 u,

where u ∈ Rn is the vector of nodal displacements uj . The discrete coun-

terpart of condition C sin(±) = 0 requires that P h2 /(EJ) coincides with

the eigenvalues of A as P varies.

h

Denoting by »min and Pcr , the smallest eigenvalue of A and the (approx-

imate) value of the critical load, respectively, then Pcr = (»min EJ)/h2 .

h

Letting θ = π/(n + 1), it can be checked (see Exercise 3, Chapter 4) that

the eigenvalues of matrix A are

»j = 2(1 ’ cos(jθ)), j = 1, . . . , n. (5.69)

The numerical calculation of »min has been carried out using the Givens

algorithm described in Section 5.10.2 and assuming n = 10. Running the

Program 42 with an absolute tolerance equal to the roundo¬ unit, the

solution »min 0.081 has been obtained after 57 iterations.

It is also interesting to analyze the case where the beam has nonuniform

cross-section, since the value of the critical load, unlike the previous situa-

tion, is not exactly known a priori. We assume that, for each x ∈ [0, L], the

section of the beam is rectangular, with depth a ¬xed and height σ that

varies according to the rule

2

S x

’1 ’1 0 ¤ x ¤ L,

σ(x) = s 1 + ,

s L

where S and s are the values at the ends, with S ≥ s > 0. The momentum

of inertia, as a function of x, is given by J(x) = (1/12)aσ 3 (x); proceeding

similarly as before, we end up with a system of linear algebraic equations

of the form

˜

Au = (P/E)h2 u,

˜

where this time A = tridiagn (b, d, b) is a tridiagonal, symmetric and posi-

tive de¬nite matrix having diagonal entries di = J(xi’1/2 ) + J(xi+1/2 ), for

i = 1, . . . , n, and o¬-diagonal entries bi = ’J(xi+1/2 ), for i = 1, . . . , n ’ 1.

238 5. Approximation of Eigenvalues and Eigenvectors

Assume the following values of the parameters: a = 0.4 [m], s = a, S =

0.5 [m] and L = 10 [m]. To ensure a correct dimensional comparison, we

¯

have multiplied by J = a4 /12 the smallest eigenvalue of the matrix A in the

uniform case (corresponding to S = s = a), obtaining »min = 1.7283 · 10’4 .

Running Program 42, with n = 10, yields in the nonuniform case the value

»min = 2.243 · 10’4 . This result con¬rms that the critical load increases

for a beam having a wider section at x = 0, that is, the structure enters

the instability regime for higher values of the load than in the uniform

cross-section case.

5.12.2 Free Dynamic Vibration of a Bridge

We are concerned with the analysis of the free response of a bridge whose

schematic structure is shown in Figure 5.5. The number of the nodes of

the structure is equal to 2n while the number of the beams is 5n. Each

horizontal and vertical beam has a mass equal to m while the diagonal

√

beams have mass equal to m 2. The sti¬ness of each beam is represented

by the spring constant κ. The nodes labeled by “0” and “2n + 1” are

constrained to ground.

2n ’ 2 2n ’ 1 2n

n+2

n+1 n+3

0 2n + 1

n’2 n’1

3 n

1 2

FIGURE 5.5. Schematic structure of a bridge

Denoting by x and y the vectors of the 2n nodal horizontal and vertical

displacements the free response of the bridge can be studied by solving the

generalized eigenvalue problems

Mx = »Kx, My = »Ky, (5.70)

√

2, b = (β, . . . , β)T ∈

where M = mdiag2n (±, b, ±, γ, b, γ), where ± = 3 +

√ √

Rn’2 with β = 3/2 + 2 and γ = 1 + 2,

K11 K12

K=κ

K12 K11

for a positive constant κ and where K12 = tridiagn (’1, ’1, ’1), K11 =

tridiagn (’1, d, ’1) with d = (4, 5, . . . , 5, 4)T ∈ Rn . The diagonal matrix

M is the mass matrix while the symmetric and positive de¬nite matrix K

is the sti¬ness matrix.

5.12 Applications 239

100

90

80

70

60

50

40

30

20

10

0

0 10 20 30 40 50 60 70 80 90 100

FIGURE 5.6. Iterations number of the Lanczos method and of the inverse power

method versus the size 2n of matrix C. The solid and the dash-dotted curves

˜ ˜

refer to the inverse power method (for »2n and »2n’1 respectively), while the

˜

dashed and the dotted curves refer to the Lanczos method (still for »2n and

˜

»2n’1 , respectively)

For k = 1, . . . , 2n we denote by (»k , zk ) any eigenvalue/eigenvector pair

√

of (5.70) and call ωk = »k the natural frequencies and zk the modes of

vibration of the bridge. The study of the free vibrations is of primary im-

portance in the design of a structure like a bridge or a multi-story building.

Indeed, if the excitation frequency of an external force (vehicles, wind or,

even worse, an earthquake) coincides with one of the natural frequencies of

the structure then a condition of resonance occurs and, as a result, large

oscillations may dangerously arise.

Let us now deal with the numerical solution of the matrix eigenvalue

problem (5.70). For this purpose we introduce the change of variable z =

M1/2 x (or z = M1/2 y) so that each generalized eigenvalue problem in (5.70)

can be conveniently reformulated as

˜

Cz = »z

where » = 1/» and the matrix C = M’1/2 KM’1/2 is symmetric positive

˜

de¬nite. This property allows us to use the Lanczos method described in

Section 5.11 and also ensures quadratic convergence of the power iterations

(see Section 5.11).

˜ ˜

We approximate the ¬rst two subdominant eigenvalues »2n and »2n’1

of the matrix C (i.e., its smallest and second smallest eigenvalues) in the

case m = κ = 1 using the de¬‚ation procedure considered in Remark 5.3.

The inverse power iteration and the Lanczos method are compared in the

˜ ˜

computation of »2n and »2n’1 in Figure 5.6.

The results show the superiority of the Lanczos method over the inverse

iterations only when the matrix C is of small size. This is to be ascribed

to the fact that, as n grows, the progressive in¬‚uence of the rounding er-

240 5. Approximation of Eigenvalues and Eigenvectors

rors causes a loss of mutual orthogonality of the Lanczos vectors and, in

turn, an increase in the number of iterations for the method to converge.

Suitable reorthogonalization procedures are thus needed to improve the

performances of the Lanczos iteration as pointed out in Section 5.11.

We conclude the free response analysis of the bridge showing in Figure

5.7 (in the case n = 5, m = 10 and κ = 1) the modes of vibration z8

and z10 corresponding to the natural frequencies ω8 = 990.42 and ω10 =

2904.59. The MATLAB built-in function eig has been employed to solve

the generalized eigenvalue problems (5.70) as explained in Section 5.9.1.

1.5 1.5

1 1

0.5 0.5

0 0

’0.5 ’0.5

’1 0 1 2 3 4 5 6 7 ’1 0 1 2 3 4 5 6 7

FIGURE 5.7. Modes of vibration corresponding to the natural frequencies ω8

(left) and ω10 (right). The undeformed con¬guration of the bridge is drawn in

dotted line

5.13 Exercises

1. Using the Gershgorin theorems, localize the eigenvalues of the matrix A

which is obtained setting A = (P’1 DP)T and then a13 = 0, a23 = 0, where

D=diag3 (1, 50, 100) and

®

1 1 1

P = 10 20 30 .

° »

100 50 60

[Solution : σ(A) = {’151.84, 80.34, 222.5}.]

2. Localize the spectrum of the matrix

®

’1

12

A= 2 7 0 .

° »

’1 0 5

[Solution : σ(A) ‚ [’2, 9].]

5.13 Exercises 241

3. Draw the oriented graph of the matrix

®

13 0

A = 0 2 ’1 .

° »

’1 0 2

4. Check if the following matrices are reducible.

® ®