Systems of Linear Equations

Systems of Linear Equations

Linear systems Ax = b arise in nearly every area of numerical computing: curve fitting, network flow, numerical PDEs, economics, and optimization. Solving them efficiently and accurately is a core topic of numerical analysis.

Direct Methods

Direct methods produce exact answers in finite steps (ignoring round-off). The two principal direct methods are Gaussian elimination and LU decomposition. They are used when the matrix is moderately sized and dense.

Gaussian Elimination

Gaussian elimination reduces A to upper triangular form by row operations, then solves the triangular system by back substitution. Complexity is O(n³) for an n×n matrix. The method is simple but may fail if a pivot element becomes zero.

Pivoting Strategies

To avoid zero or small pivots, partial pivoting swaps rows so the largest magnitude element in the current column becomes the pivot. Full pivoting swaps both rows and columns. Pivoting dramatically improves numerical stability; without it, small round-off errors can catastrophically amplify.

Gauss–Jordan Elimination

Gauss–Jordan continues elimination until the coefficient matrix becomes identity; the right-hand side then contains the solution directly. Computing matrix inverses by Gauss–Jordan on [A | I] yields [I | A^-1]. It uses roughly 50% more operations than Gaussian elimination and is mainly used when the inverse is explicitly needed.

LU Decomposition

LU decomposition factors A as the product of a lower triangular L and an upper triangular U. Once computed, any right-hand side b is solved by forward substitution (Ly = b) and back substitution (Ux = y). This is efficient when solving multiple systems with the same matrix, such as computing A^-1 column by column or handling time-stepping problems. Variants include Doolittle, Crout, and Cholesky (for symmetric positive-definite matrices).

Iterative Methods

For large sparse systems, iterative methods produce successive approximations x^(k) that converge to the solution. They require much less memory than direct methods and exploit sparsity. Common methods include Jacobi, Gauss–Seidel, and Successive Over-Relaxation (SOR).

Jacobi Method

Jacobi iteration updates each x_i using the previous iteration's values: x_i^(k+1) = (b_i − Σ_j≠i a_ij x_j^(k)) / a_ii. Convergence is guaranteed for diagonally dominant matrices. Jacobi is easy to parallelize because every component is computed independently.

Gauss–Seidel Method

Gauss–Seidel updates components in place, using the most recent values. It typically converges in half as many iterations as Jacobi but is harder to parallelize. Successive Over-Relaxation (SOR) introduces a relaxation parameter ω to further accelerate convergence; optimal ω lies in (1, 2) for most problems.

Condition Number and Accuracy

The condition number κ(A) = ‖A‖ · ‖A^-1‖ quantifies sensitivity of x to perturbations in b. A large condition number signals an ill-conditioned system where small round-off can produce large solution errors. Scaling, equilibration, and iterative refinement are techniques to improve computed accuracy.

Summary

Linear system solvers are workhorses of scientific computing. Direct methods are standard for moderate dense systems; iterative methods dominate for large sparse ones. The choice depends on matrix size, structure, sparsity, and required accuracy.