3. Numerical Linear Algebra and Eigenvalue Problems

Introduction and Context

Many engineering models become, after discretization, systems of algebraic equations. Finite-element analysis, finite-difference schemes, computational fluid dynamics, structural optimization, heat transfer, and control-system simulation commonly lead to a linear system

\[

A\mathbf{x}=\mathbf{b},

\]

where \(A\) is an \(n\times n\) matrix, \(\mathbf{x}\) is the vector of unknowns, and \(\mathbf{b}\) is the load, source, or excitation vector. In dynamic and stability analyses, the central problem often changes form to an eigenvalue equation such as

\[

A\mathbf{v}=\lambda\mathbf{v}

\]

or, for mechanical vibration,

\[

K\boldsymbol{\phi}=\omega^2M\boldsymbol{\phi}.

\]

Here \(K\) and \(M\) are stiffness and mass matrices, \(\omega\) is a natural frequency, and \(\boldsymbol{\phi}\) is a mode shape.

Numerical linear algebra provides the algorithms, error analysis, and implementation principles required to solve these problems reliably. At postgraduate level, it is not sufficient to know that a software package returns a solution. The engineer must understand conditioning, sparsity, stability, convergence, preconditioning, and the relationship between mathematical results and physical interpretation.

Learning Objectives

After studying this chapter, the reader should be able to:

1. Formulate engineering models as dense or sparse linear systems.

2. Distinguish direct, iterative, stationary, and Krylov-subspace methods.

3. Explain conditioning, stability, residuals, and backward error.

4. Select and apply LU, Cholesky, QR, Jacobi, Gauss–Seidel, CG, and GMRES methods appropriately.

5. Explain the purpose and operation of preconditioners.

6. Formulate standard, generalized, symmetric, and nonsymmetric eigenvalue problems.

7. Compute and interpret eigenvalues, eigenvectors, natural frequencies, and vibration modes.

8. Diagnose numerical difficulties and make informed solver choices for engineering applications.

Algebraic Structure of Engineering Systems

A matrix is not merely a rectangular collection of numbers. Its structure often reflects the physics and discretization of the problem. A matrix may be dense, sparse, symmetric, positive definite, banded, diagonally dominant, or indefinite.

A matrix \(A\) is symmetric if \(A=A^T\). It is positive definite if \(\mathbf{z}^TA\mathbf{z}>0\) for every nonzero vector \(\mathbf{z}\). Stiffness matrices for stable, properly constrained elastic systems are frequently symmetric positive definite. Such structure permits efficient and robust algorithms. By contrast, convection-dominated flow problems may produce nonsymmetric matrices, while constrained mechanics and saddle-point formulations may produce indefinite matrices.

The condition number in a chosen norm is

\[

\kappa(A)=\|A\|\|A^{-1}\|.

\]

A large condition number means that small perturbations in data or round-off errors may cause relatively large changes in the solution. This is a property of the problem, not necessarily a defect in the algorithm. An algorithm can be numerically stable while solving an ill-conditioned problem whose answer is intrinsically sensitive.

The residual of an approximate solution \(\hat{\mathbf{x}}\) is

\[

\mathbf{r}=\mathbf{b}-A\hat{\mathbf{x}}.

\]

A small residual indicates that \(\hat{\mathbf{x}}\) nearly satisfies the equations. However, a small residual does not always imply a small solution error when \(A\) is ill-conditioned. For this reason, residual monitoring should be combined with scaling, condition estimates, and, where possible, backward-error analysis.

Direct Solution of Linear Systems

Direct methods produce a solution through a finite sequence of algebraic operations, although floating-point arithmetic introduces approximation. They are often preferred for small and medium systems, multiple right-hand sides, and problems where predictable behavior is important.

Gaussian Elimination and LU Factorization

Gaussian elimination transforms \(A\mathbf{x}=\mathbf{b}\) into an upper triangular system \(U\mathbf{x}=\mathbf{y}\). The solution is then obtained by back substitution. In matrix form,

\[

PA=LU,

\]

where \(P\) is a permutation matrix, \(L\) is lower triangular, and \(U\) is upper triangular. Partial pivoting selects a suitably large pivot in each column and reduces the risk of division by a small number. In practical computation, pivoting is normally essential unless special matrix properties guarantee stability.

For one matrix with several load vectors, factorization is performed once and triangular substitutions are repeated. This is particularly valuable in structural analysis with multiple load cases. The approximate arithmetic cost for a dense factorization is proportional to \(n^3\), while storage is proportional to \(n^2\). These costs become prohibitive for very large systems.

Cholesky Factorization

If \(A\) is symmetric positive definite, it can be written as

\[

A=LL^T,

\]

where \(L\) is lower triangular. Cholesky factorization requires roughly half the arithmetic and storage of general LU factorization and avoids general pivoting. It is a natural method for many finite-element stiffness systems after appropriate treatment of boundary conditions.

The method must not be used merely because a matrix appears symmetric. Numerical symmetry and positive definiteness should be verified from the model and, when appropriate, checked computationally. An unconstrained structural model may possess rigid-body modes, making its stiffness matrix singular rather than positive definite.

QR Factorization

QR factorization expresses a matrix as

\[

A=QR,

\]

where \(Q\) is orthogonal and \(R\) is upper triangular. Orthogonal transformations preserve the Euclidean norm and are generally more numerically stable than normal-equation approaches for least-squares problems. QR is therefore preferred when data fitting or rank deficiency makes stability important, although it is usually more expensive than normal equations or LU.

Sparse Direct Methods

Discretized partial differential equations and finite-element models often generate sparse matrices: most entries are zero. Sparse direct solvers exploit this structure, but elimination can create new nonzero entries, called fill-in. Ordering strategies such as approximate minimum degree and nested dissection reduce fill-in and improve computational efficiency. The choice of ordering may matter as much as the factorization itself.

Iterative Methods for Large Systems

Iterative methods begin with an initial estimate \(\mathbf{x}_0\) and generate a sequence of approximations. They are attractive when matrices are very large and sparse, because the dominant operation can be a matrix–vector product rather than explicit factorization. They also permit early termination when the engineering accuracy requirement has been reached.

Stationary Iterations

For Jacobi iteration, write \(A=D+L+U\), where \(D\) is diagonal. The update is

\[

\mathbf{x}^{(k+1)}=D^{-1}[\mathbf{b}-(L+U)\mathbf{x}^{(k)}].

\]

Jacobi is simple and highly parallelizable, but often converges slowly.

Gauss–Seidel uses newly computed components immediately:

\[

(D+L)\mathbf{x}^{(k+1)}=\mathbf{b}-U\mathbf{x}^{(k)}.

\]

It may converge faster than Jacobi, but its sequential dependence limits straightforward parallel implementation. Both methods require suitable spectral conditions; diagonal dominance or positive-definite structure often supports convergence, but neither method is universally reliable.

Successive over-relaxation introduces a parameter \(\omega_r\). Values between zero and two may accelerate convergence for suitable problems, but a poor parameter can cause stagnation or divergence.

Conjugate Gradient Method

The conjugate gradient (CG) method is designed for symmetric positive-definite matrices. Rather than making residuals orthogonal in the ordinary sense alone, it constructs search directions that are conjugate with respect to \(A\):

\[

\mathbf{p}_i^TA\mathbf{p}_j=0\quad (i\ne j).

\]

In exact arithmetic, CG reaches the solution in at most \(n\) steps. In floating-point arithmetic, practical convergence depends strongly on the eigenvalue distribution and condition number. A standard bound is

\[

\|\mathbf{x}_k-\mathbf{x}_\|_A\leq2\left(\frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1}\right)^k\|\mathbf{x}_0-\mathbf{x}_\|_A.

\]

This bound is conservative, but it explains why reducing the effective condition number improves convergence.

GMRES and Other Krylov Methods

For nonsymmetric systems, GMRES constructs an approximation in the Krylov subspace

\[

\mathcal{K}_m(A,\mathbf{r}_0)=\operatorname{span}\{\mathbf{r}_0,A\mathbf{r}_0,\ldots,A^{m-1}\mathbf{r}_0\}

\]

and minimizes the residual norm over that subspace. Full GMRES normally requires increasing storage and orthogonalization, so restarted GMRES(m) is commonly used. Restarting reduces memory but may slow convergence or cause stagnation.

BiCGSTAB and related methods use less memory and can be effective for nonsymmetric problems, but their convergence may be less smooth and can be sensitive to breakdown or loss of accuracy. The method should therefore be selected with reference to matrix structure and observed convergence, not by popularity alone.

Preconditioning: Making the Problem Easier

A preconditioner is a matrix or operator \(M\) that approximates \(A\) but is cheaper to apply or factorize. Instead of solving \(A\mathbf{x}=\mathbf{b}\) directly, an iterative method solves an equivalent preconditioned system. Left preconditioning gives

\[

M^{-1}A\mathbf{x}=M^{-1}\mathbf{b}.

\]

Right preconditioning gives

\[

AM^{-1}\mathbf{y}=\mathbf{b},\qquad \mathbf{x}=M^{-1}\mathbf{y}.

\]

An effective preconditioner reduces the number of iterations while keeping setup and application costs reasonable. It need not approximate \(A\) accurately in every entry; its purpose is to improve the spectral or algebraic properties relevant to convergence.

Common choices include diagonal scaling, incomplete LU (ILU), incomplete Cholesky for symmetric positive-definite systems, block preconditioners, and multigrid methods. Incomplete factorizations limit fill-in through a prescribed sparsity pattern or dropping rule. Algebraic multigrid is particularly powerful for many elliptic partial differential equations because it reduces both high- and low-frequency error components at different levels.

Preconditioning changes the residual that an algorithm minimizes, so stopping criteria must be interpreted carefully. A robust implementation monitors the original residual as well as the preconditioned residual when practical.

Eigenvalue Problems and Modal Analysis

An eigenvalue \(\lambda\) and nonzero eigenvector \(\mathbf{v}\) satisfy

\[

A\mathbf{v}=\lambda\mathbf{v}.

\]

Eigenvalues characterize amplification, decay, resonance, stability, and dominant modes. For a general matrix, eigenvalues may be real or complex, and eigenvectors may be nonorthogonal. For a symmetric matrix, eigenvalues are real and eigenvectors can be chosen orthogonal. If \(A\) is symmetric positive definite, all eigenvalues are positive.

The generalized eigenproblem is

\[

A\mathbf{v}=\lambda B\mathbf{v}.

\]

When \(B\) is symmetric positive definite and \(A\) is symmetric, eigenvectors can be chosen to satisfy \(B\)-orthogonality:

\[

\mathbf{v}_i^TB\mathbf{v}_j=0,\quad i\ne j.

\]

Mechanical Vibration

For an undamped linear mechanical system,

\[

M\ddot{\mathbf{u}}+K\mathbf{u}=\mathbf{0}.

\]

Assuming \(\mathbf{u}(t)=\boldsymbol{\phi}e^{i\omega t}\) produces

\[

(K-\omega^2M)\boldsymbol{\phi}=\mathbf{0}.

\]

Thus \(\lambda=\omega^2\). The natural frequency is \(f=\omega/(2\pi)\). Mode shapes are defined only up to a scalar factor, so normalization may be selected for convenience. Mass normalization satisfies

\[

\boldsymbol{\phi}_i^TM\boldsymbol{\phi}_j=\delta_{ij}.

\]

With this normalization, modal coordinates decouple the undamped equations. This is central to modal superposition, model reduction, response prediction, and experimental correlation.

Computing Eigenvalues

For dense matrices, the QR algorithm is a standard general-purpose method. For symmetric matrices, tridiagonal reduction followed by a symmetric QR or divide-and-conquer method is efficient and reliable. For very large sparse systems, it is usually unnecessary to compute the complete spectrum. Krylov methods such as Lanczos for symmetric problems and Arnoldi for nonsymmetric problems target a selected portion of the spectrum.

The power method is simple and useful for estimating the dominant eigenvalue in magnitude. Its convergence is slow when the largest and second-largest eigenvalues are close. Shift-and-invert methods transform the problem so that eigenvalues near a chosen shift become dominant. For a standard problem,

\[

(A-\sigma I)^{-1}\mathbf{v}=\mu\mathbf{v},

\qquad \mu=\frac{1}{\lambda-\sigma}.

\]

This approach is powerful for extracting low-frequency vibration modes, but each iteration requires solving a shifted linear system. Consequently, eigenvalue computation and linear-system solution are closely connected.

Practical Solver Selection

A useful decision process begins with matrix properties rather than problem size alone.

• A dense, moderate-size general system commonly uses pivoted LU.
• A symmetric positive-definite system commonly uses Cholesky or preconditioned CG.
• A very large sparse nonsymmetric system commonly uses GMRES or BiCGSTAB with ILU or a block preconditioner.
• A symmetric generalized vibration problem commonly uses a Lanczos-type method with mass-aware orthogonalization.
• A least-squares problem with possible rank deficiency is usually better served by QR or singular-value methods than by normal equations.

Scaling is frequently overlooked. If rows or columns differ by many orders of magnitude, diagonal equilibration can improve numerical behavior and preconditioner performance. Boundary conditions must also be imposed consistently. Incorrect elimination of constrained degrees of freedom can produce singular matrices, spurious modes, or physically meaningless reactions.

Stopping tolerances should reflect engineering accuracy. A tolerance of \(10^{-12}\) is not automatically superior to \(10^{-6}\); excessive iteration may waste resources and expose the computation to accumulated round-off error. Conversely, a small residual may still be inadequate for an ill-conditioned quantity of interest. Verification should include mesh refinement, balance checks, symmetry checks, and comparison with analytical or experimental benchmarks where available.

Illustrative Case Study: Modal Analysis of a Flexible Support

Background

Consider an illustrative finite-element model of a flexible machine support. The model has several thousand displacement degrees of freedom. The objective is to determine the first six natural frequencies so that operating excitation frequencies can be kept away from structural resonances.

Problem and Initial Diagnosis

The assembled matrices satisfy

\[

K\boldsymbol{\phi}=\omega^2M\boldsymbol{\phi}.

\]

The stiffness matrix is sparse and symmetric. The mass matrix is symmetric positive definite after removal of constrained degrees of freedom. A dense eigensolver would require excessive memory, while computing the full spectrum would provide information that is not needed.

Computational Approach

A Lanczos-based sparse eigensolver is selected to target the lowest eigenvalues. Because low frequencies are required, a shift-and-invert transformation is used with a shift near zero. Each iteration requires solving a system involving \(K-\sigma M\). A sparse Cholesky factorization is used because the shifted matrix retains a suitable symmetric positive-definite structure for the selected shift. Nested-dissection ordering limits fill-in.

The eigenvectors are orthogonalized in the mass inner product. After convergence, each computed pair is checked using the generalized residual

\[

\mathbf{r}_i=K\boldsymbol{\phi}_i-\omega_i^2M\boldsymbol{\phi}_i.

\]

The frequencies are normalized and compared with a reduced hand calculation and a refined mesh.

Outcome and Interpretation

The computed modes are accepted only after residuals, mass orthogonality, and mesh sensitivity meet predetermined engineering criteria. The first mode is identified as a global bending mode, while higher modes involve torsional and local deformation. An operating speed that lies close to the second natural frequency is rejected, and a design modification is evaluated to increase stiffness in the region associated with that mode.

Learning Points

The case illustrates that eigenvalue analysis is not simply an invocation of an eigensolver. Matrix structure determines the algorithm; the eigenvalue region determines the transformation; factorization cost determines feasibility; and physical interpretation requires examination of mode shapes, not frequencies alone. It also demonstrates the dual role of linear algebra: the eigenvalue method depends on repeated solution of linear systems.

Challenges, Errors, and Professional Responsibilities

Common mistakes include applying CG to a nonsymmetric matrix, using Cholesky without confirming positive definiteness, ignoring pivoting in LU, forming an explicit inverse, treating a small residual as proof of accuracy, and computing many eigenvalues when only a few are required. Another serious error is failing to distinguish rigid-body modes from numerical singularities in structural models.

Floating-point cancellation, loss of orthogonality, poor scaling, and excessive fill-in can all degrade performance. Iterative methods may stagnate because the preconditioner is ineffective, the tolerance is inappropriate, or the matrix has unfavorable spectral properties. Eigenvectors associated with clustered or repeated eigenvalues can rotate substantially under small perturbations; the invariant subspace may be more meaningful than an individual vector.

In professional practice, numerical results must be traceable. Solver versions, tolerances, scaling, ordering, convergence histories, matrix checks, and boundary-condition treatments should be recorded. Safety-critical decisions should not rely on an unexplained software default. Verification, validation, independent calculations, and sensitivity studies are part of engineering governance.

Digital Implementation and Emerging Directions

Modern scientific-computing libraries provide sparse storage formats, threaded factorization, GPU-enabled matrix–vector products, automatic differentiation, and distributed Krylov solvers. The main computational kernel in many iterative methods is sparse matrix–vector multiplication, whose performance is often limited by memory bandwidth rather than arithmetic throughput. Matrix-free methods avoid assembling \(A\) explicitly by evaluating its action on a vector; they are valuable in very large discretizations and high-order methods.

Adaptive preconditioning, hierarchical matrices, domain decomposition, and multilevel solvers continue to expand the feasible problem size. Machine-learning techniques are being explored for predicting solver parameters and preconditioner choices, but such approaches require careful validation. A learned decision rule cannot replace residual checks, physical verification, or mathematical guarantees.

Best Practices

1. Inspect symmetry, definiteness, sparsity, scaling, and conditioning before choosing a solver.

2. Exploit physical and mathematical structure rather than treating every matrix as general and dense.

3. Prefer factorization reuse for multiple right-hand sides.

4. Use preconditioning deliberately and monitor the original residual.

5. Compute only the eigenvalues and modes required by the engineering question.

6. Check residuals, orthogonality, sensitivity, and physical plausibility.

7. Treat boundary conditions and units as part of numerical correctness.

8. Document tolerances, ordering, stopping criteria, and convergence evidence.

Chapter Summary

Engineering discretizations produce linear systems and eigenvalue problems whose reliable solution depends on matrix structure, conditioning, and numerical stability. Direct methods such as pivoted LU, Cholesky, and QR are robust choices for appropriate dense or moderate-size problems, while iterative Krylov methods are essential for very large sparse systems. Preconditioning accelerates convergence by transforming the algebraic problem without changing its exact solution. Eigenvalue methods provide stability information, natural frequencies, and mode shapes; generalized symmetric problems are especially important in structural dynamics. Effective practice combines algorithm selection with scaling, residual analysis, convergence monitoring, verification, and physical interpretation. Modern sparse, parallel, matrix-free, and multilevel technologies extend these ideas to large-scale computational engineering.

Key Terms

Condition number: A measure of the sensitivity of a linear-system solution to perturbations.

Residual: The vector \(\mathbf{b}-A\hat{\mathbf{x}}\) measuring equation imbalance.

Backward error: The size of a perturbation to the data for which a computed result is an exact solution.

Pivoting: Reordering rows or columns to obtain numerically safer pivots.

LU factorization: Decomposition of a matrix into lower and upper triangular factors, usually with permutation.

Cholesky factorization: Factorization \(A=LL^T\) for symmetric positive-definite matrices.

QR factorization: Decomposition into an orthogonal factor and an upper triangular factor.

Sparsity: A matrix property in which most entries are zero.

Fill-in: New nonzero entries created during sparse factorization.

Krylov subspace: The span generated by successive matrix actions on a starting vector.

Preconditioner: An operator that improves the convergence properties of an iterative method.

Conjugate gradient: A Krylov method for symmetric positive-definite systems.

GMRES: A Krylov method that minimizes the residual for general nonsymmetric systems.

Eigenvalue: A scalar satisfying \(A\mathbf{v}=\lambda\mathbf{v}\) for a nonzero vector.

Eigenvector: A nonzero vector associated with an eigenvalue.

Generalized eigenproblem: A problem of the form \(A\mathbf{v}=\lambda B\mathbf{v}\).

Mode shape: A spatial deformation pattern associated with a natural frequency.

Shift-and-invert: An eigenvalue transformation used to target eigenvalues near a selected shift.

Mass orthogonality: Orthogonality defined by the inner product \(\mathbf{u}^TM\mathbf{v}\).

Review Questions

Short Answer Questions

1. What information does the residual of an approximate linear-system solution provide?

2. Why is pivoting normally used in Gaussian elimination?

3. State the matrix properties required by the conjugate gradient method.

4. What is the purpose of a preconditioner?

5. How are natural frequencies related to generalized eigenvalues in an undamped mechanical system?

Long Answer Questions

1. Explain LU, Cholesky, and QR factorization, including the matrix structures and engineering problems for which each is appropriate.

2. Compare Jacobi, Gauss–Seidel, conjugate gradient, and GMRES with respect to convergence, storage, and matrix requirements.

3. Discuss conditioning, residuals, backward error, and the distinction between problem sensitivity and algorithmic stability.

4. Develop the generalized eigenvalue formulation for finite-element vibration analysis and explain mass normalization.

5. Explain sparse direct methods, fill-in, ordering, and preconditioning in the solution of large engineering systems.

Analytical / Application Questions

1. A symmetric stiffness matrix has several zero eigenvalues before boundary conditions are applied. Explain the likely physical cause, the numerical consequence, and the appropriate corrective action.

2. A GMRES computation reaches a small preconditioned residual but the original residual remains large. Diagnose possible causes and propose a verification procedure.

3. An engineer must compute the first four natural frequencies of a very large sparse model. Select an eigensolver, transformation, linear solver, and validation checks, and justify each choice.

22. Multiple Choice Questions with Answers

Question 1

Which method is specifically designed for a symmetric positive-definite linear system?

A. Conjugate gradient

B. GMRES

C. Power iteration

D. Unsymmetric QR only

Correct Answer: A. Conjugate gradient

Question 2

What is the principal purpose of pivoting in LU factorization?

A. To make the matrix sparse

B. To reduce the risk associated with small pivots

C. To compute eigenvectors directly

D. To eliminate all round-off error

Correct Answer: B. To reduce the risk associated with small pivots

Question 3

A small residual may fail to indicate a small solution error primarily when the matrix is:

A. Diagonal

B. Well scaled

C. Ill-conditioned

D. Positive definite

Correct Answer: C. Ill-conditioned

Question 4

Which factorization is most directly suited to a symmetric positive-definite matrix?

A. Cholesky

B. Incomplete QR only

C. Schur without reduction

D. Fourier

Correct Answer: A. Cholesky

Question 5

In a generalized vibration problem Kφ = λMφ, what does λ represent for an undamped system?

A. The damping ratio

B. The square of the natural circular frequency

C. The modal mass only

D. The applied force

Correct Answer: B. The square of the natural circular frequency

Question 6

What is the main reason for using a preconditioner?

A. To change the physical model

B. To reduce iterative convergence difficulty

C. To guarantee zero round-off error

D. To convert every matrix into a diagonal matrix

Correct Answer: B. To reduce iterative convergence difficulty

Question 7

Which method minimizes a residual over a Krylov subspace for a general nonsymmetric system?

A. GMRES

B. Cholesky

C. Jacobi factorization

D. Classical modal normalization

Correct Answer: A. GMRES

Question 8

In sparse direct factorization, fill-in refers to:

A. Rounding entries to zero

B. New nonzero entries created during elimination

C. Additional physical loads

D. Repeated eigenvalues

Correct Answer: B. New nonzero entries created during elimination

Question 9

Why is shift-and-invert useful in eigenvalue computation?

A. It targets eigenvalues near a selected shift

B. It removes the need for linear solves

C. It guarantees all eigenvalues are positive

D. It applies only to diagonal matrices

Correct Answer: A. It targets eigenvalues near a selected shift

Question 10

For a mechanical mode shape, mass normalization commonly imposes:

A. φᵀφ = 0 for every mode

B. φᵀMφ = 1 for each normalized mode

C. Kφ = 0 for every mode

D. Mφ = λφ regardless of K

Correct Answer: B. φᵀMφ = 1 for each normalized mode