<< . .

. 35
( : 95)



. . >>

sequences {»1 (Hm )} and {»m (Hm )} to »1 and »n , respectively, is governed
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.
®  ® 

<< . .

. 35
( : 95)



. . >>