Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 8

Linear Systems and Eigenvalue Problems

The companion programs show how elimination, relaxation, tridiagonal solves, and eigenvalue iterations work on small matrices. For a larger physics calculation, first identify whether the matrix is dense or sparse, symmetric or general, and whether one solution or many are needed. The resources below connect those choices to tested routines, provide independently checked results for the companion examples, and explain what the numerical checks can and cannot establish.

Chapter 8 codes Chapter 8 exercises & projects Chapter 8 updates & errata Chapter 8 overview

Software and library resources

The companion implementations expose the algorithms. Library routines offer pivoting, structure-aware factorizations, sparse storage, and convergence reporting for practical calculations. The routines in a row solve the same kind of problem but need not use the identical algorithm. For example, a sparse direct solver does not necessarily perform the Thomas algorithm. Check the routine’s status and its assumptions before interpreting the result. Appendix B of the book discusses additional routines.

On a narrow screen, scroll horizontally to view all language columns.

Numerical task MATLAB Python C++ Fortran
One dense linear system Solve A*x = b directly without forming the inverse. A\b numpy.linalg.solve Eigen PartialPivLU for an invertible matrix LAPACK dgesv
Tridiagonal system Retain the three-diagonal structure instead of storing a full matrix. spdiags to build a sparse matrix, then A\b scipy.linalg.solve_banded((1,1),ab,b) Eigen SparseLU for a general sparse solve LAPACK dgtsv
Several right-hand sides with the same matrix Factor once, then reuse the factors for each solve. dA = decomposition(A), then dA\b lu_factor and lu_solve Reuse an Eigen PartialPivLU object LAPACK dgetrf and dgetrs
Large sparse, symmetric positive-definite system Use a preconditioned iterative method when direct storage is costly. pcg scipy.sparse.linalg.cg Eigen ConjugateGradient PETSc KSPCG
All eigenvalues and eigenvectors of a dense symmetric matrix A tested full-spectrum solver for normal modes or energy levels. eig numpy.linalg.eigh Eigen SelfAdjointEigenSolver LAPACK dsyevd
A few eigenvalues and eigenvectors of a large sparse symmetric matrix Request selected modes rather than the full spectrum. eigs scipy.sparse.linalg.eigsh Spectra SymEigsSolver SLEPc EPS

Conjugate gradient requires a symmetric positive-definite matrix; the companion SOR example also requires strict diagonal dominance under its input check. For a general sparse system, choose a suitable sparse direct or Krylov solver, such as PETSc KSP. PETSc is a larger external package with Fortran bindings; it is listed for practical large-scale work rather than as a small built-in Fortran routine. A solver for a selected eigenvalue may need a shift to reach an interior eigenvalue, and a shift close to an eigenvalue makes the associated linear solve ill-conditioned. Numerical-library tolerances and returned status codes should be checked separately from the tolerances in the companion programs.

Benchmark matrices and reference results

The linear-system values below were checked with exact rational elimination. The eigenvalue formula follows from the characteristic polynomial of the displayed matrix and was checked numerically with an independent symmetric eigensolver. These are the companion programs’ default inputs; several differ from the book’s printed examples. Verify computed results with the original equations, before rounding the displayed values.

Matrix problem Reference result Why it is useful
Gaussian elimination: A = [0 1 2; 2 1 4; 2 4 6], b = [4,3,7] Exact solution x = [−5/2, 0, 2]. The first pivot is zero, so row pivoting is necessary. Check A*x − b.
SOR: A = [8 2 4; 2 6 1; 1 1 8], b = [2,6,4] Exact solution x = [−1/5, 1, 2/5]. Compare the reported residual and iteration count after varying the relaxation parameter.
Thomas algorithm: eight rows, diagonal 4, adjacent diagonals −1, and b = [0,0,0,0,0,0,0,16] x1 = 16/40545 ≈ 0.000394623; x8 = 173824/40545 ≈ 4.287187076. Compare the first and last components, then test all eight original equations.
Matrix inverse: A = [3 2 4; 2 −3 1; 1 1 2] det(A) = −7; A−1 = [1 0 −2; 3/7 −2/7 −5/7; −5/7 1/7 13/7]. Check A*A−1 against the identity, using the original matrix.
Eigenvalue programs: A = [1 2 3; 2 2 −2; 3 −2 4] Eigenvalues (1 − √37)/2 ≈ −2.541381265, (1 + √37)/2 ≈ 3.541381265, and 6. The power method targets 6; shift 3.4 targets 3.541381265; Jacobi returns all three.

For the eigenvalue matrix, det(λI − A) = (λ − 6)(λ2 − λ − 9). Eigenvector signs are arbitrary. Compare A*v − λv, and for a full symmetric eigensystem also compare VT*V with the identity.

Choosing a linear or eigenvalue method

Problem Good method to try Notes
Moderate dense system with one right-hand side Pivoted Gaussian elimination or a library dense solver Inspect scaling and conditioning. A small residual need not imply a small solution error.
Same dense matrix with many right-hand sides Factor once and reuse LU or Cholesky factors Use Cholesky when the matrix is symmetric positive definite. Do not compute an explicit inverse merely to obtain solutions.
Tridiagonal system from nearest-neighbor coupling Thomas algorithm or a tested tridiagonal solver The companion Thomas code does not pivot; its example is strictly diagonally dominant. A library routine can handle a wider range of matrices.
Large sparse symmetric positive-definite system Preconditioned conjugate gradient Use a tested routine, preserve sparse storage, monitor convergence and the residual. The chapter discusses CG; there is no separate CG companion code.
Dominant eigenvalue in magnitude Power method for a small example It finds the dominant eigenvalue and a corresponding eigenvector when the dominant magnitude is isolated and the initial vector has a component in that direction.
Eigenvalue near a physically motivated estimate Shifted inverse iteration Factor A − μI once. An exact or nearly exact eigenvalue as shift makes this solve singular or ill-conditioned.
All eigenvalues and eigenvectors of a small real symmetric matrix Jacobi rotations Check ||A*v − λv|| for each eigenvalue and eigenvector, as well as eigenvector orthogonality. For larger dense matrices, use a symmetric library eigensolver.
Only a few modes of a large sparse symmetric matrix Lanczos-based library solver Avoid computing the full spectrum; check the requested portion of the spectrum and each mode’s residual.

The meaning of a residual depends on the scales in the physical equations. For A*x = b, compare ||A*x − b|| with ||A|| ||x|| + ||b||, and vary the input data if the matrix may be ill-conditioned. For each eigenvalue and eigenvector, check ||A*v − λv|| and whether the eigenvalue and mode fall in the physically allowed range. Refining a tolerance without examining the model and matrix conditioning cannot establish physical accuracy.

Selected references and documentation

← Back to Chapter 8