Linear systems of equations appear throughout mathematics, science, engineering, economics, and data analysis. In calculus and linear algebra courses, students learn to solve these systems using Gaussian elimination, matrix inverses, Cramer’s rule, and other symbolic techniques. For small systems with tidy coefficients, these methods work well. But in real computational practice, linear systems are often large — hundreds or thousands of unknowns — and their coefficients carry rounding errors from measurement, modeling, or prior computation. In these settings, the question is not merely whether an exact symbolic answer exists, but whether a reliable numerical answer can be found efficiently.

Chapter 9 introduces numerical linear algebra: the study of how to solve linear systems in ways that are computationally efficient, numerically stable, and sensitive to conditioning. Students will revisit Gaussian elimination from a numerical perspective, study pivoting strategies that reduce error, develop LU decomposition as a structured factorization tool, and explore iterative methods that approximate solutions through repeated refinement. Throughout the chapter, the emphasis is on understanding not just how to perform computations, but when to trust them.

Why This Chapter Matters. Large linear systems are among the most common computational problems in applied mathematics. They arise when discretizing differential equations, fitting data models, analyzing networks, computing stresses in structural engineering, modeling economic systems, training machine learning models, and simulating physical processes. No matter how carefully a system is set up, numerical errors enter the picture. Gaussian elimination, which students may regard as a clean symbolic procedure, must be performed in floating-point arithmetic — and floating-point arithmetic introduces rounding at every step. Without careful strategy, those rounding errors can grow dramatically. Conditioning tells us whether a system is inherently sensitive to input errors. Stability tells us whether an algorithm manages those errors well. This chapter teaches students how to solve linear systems numerically with awareness of both.

Linear algebra provides a powerful symbolic language for solving systems of equations. Given a system

where \(\mathbf{A}\) is an \(n \times n\) matrix and \(\mathbf{b}\) is a vector of \(n\) known values, the symbolic solution is \(\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}\), provided \(\mathbf{A}\) is invertible. For small systems, students can compute this directly or use row reduction to find exact rational solutions.

But real computational problems rarely stay small. In modern applications, \(n\) might be thousands, millions, or more. The matrix \(\mathbf{A}\) might come from a physical model, a finite difference approximation, or a data fitting procedure — and its entries are often decimal approximations that already carry rounding error. In this context, several new concerns arise.

Numerical linear algebra addresses all of these concerns. It does not abandon the symbolic framework of linear algebra — it works within that framework, but with careful attention to how computations are actually performed in finite-precision arithmetic.

To appreciate why numerical methods are needed, consider how often linear systems arise in practice.

In all of these settings, the challenge is not formulating the system — it is solving it accurately, efficiently, and with honest awareness of potential error.

Gaussian elimination is the foundational direct method for solving \(\mathbf{A}\mathbf{x} = \mathbf{b}\). Students have likely encountered it in algebra or linear algebra courses. This section reviews the procedure and sets up the numerical analysis that follows.

The goal of Gaussian elimination is to transform the augmented matrix \([\mathbf{A} \mid \mathbf{b}]\) into upper triangular form using elementary row operations, then recover the solution by back substitution.

The three elementary row operations are:

• Swap two rows.
• Multiply a row by a nonzero scalar.
• Add a multiple of one row to another.

None of these operations change the solution set of the system.

The elimination process. In the first stage, the first column below the diagonal is eliminated by subtracting multiples of the first row from each subsequent row. In the second stage, the second column below the diagonal is eliminated using the second row as a pivot. This continues until the matrix is upper triangular.

In the basic formulation of Gaussian elimination, the pivot element at each stage is the diagonal entry of the current row. This works well when all diagonal entries are nonzero and have reasonably large magnitudes. But in practice, two problems can occur.

Zero pivots. If the pivot element is zero, the algorithm breaks down because we would be dividing by zero. A row swap can resolve this — we swap the current row with a lower row that has a nonzero entry in the pivot column.

Small pivots. Even more insidious is the case of small but nonzero pivots. When the multiplier \(m_{ij} = a_{ij} / a_{jj}\) is computed with a very small \(a_{jj}\), the multiplier becomes very large. Large multipliers amplify rounding errors that were introduced in earlier steps, potentially causing the computed solution to be wildly inaccurate.

The standard remedy is partial pivoting: at each stage of elimination, before using the current diagonal entry as a pivot, scan the entries in the pivot column at or below the diagonal and swap the current row with the row containing the largest absolute value in that column.

This strategy ensures that the pivot element is as large as possible in magnitude, which keeps the multipliers bounded in absolute value by 1. Bounded multipliers mean that rounding errors do not grow explosively during elimination.

A refinement of partial pivoting is scaled partial pivoting: before comparing entries in the pivot column, scale each row by dividing by the largest absolute value in that row. This accounts for the possibility that different equations have very different scales, which can cause partial pivoting alone to give suboptimal results. Most professional numerical software libraries implement scaled partial pivoting or an equivalent strategy automatically.

Even with partial pivoting, Gaussian elimination is not perfectly immune to roundoff error. Understanding how errors enter helps students interpret the reliability of computed solutions.

Each floating-point multiplication and addition introduces a relative error on the order of machine precision \(\epsilon_{\text{mach}}\). Over the course of \(n\) elimination steps, these errors accumulate. In the worst case, the accumulated error in the computed solution \(\hat{\mathbf{x}}\) satisfies a bound of the form:

where \(\kappa(\mathbf{A})\) is the condition number of the matrix \(\mathbf{A}\) (developed in Section 9.8) and \(\epsilon_{\text{mach}}\) is the smallest floating-point number that, when added to 1, produces a result distinguishable from 1 on the computing system.

This bound has a critical implication: even if the algorithm is perfectly implemented and the arithmetic is as accurate as the hardware allows, the accuracy of the solution is fundamentally limited by how well-conditioned or ill-conditioned the system is.

A key practical tool is the residual: after computing a solution \(\hat{\mathbf{x}}\), compute

If \(\mathbf{r}\) is close to zero, the computed solution satisfies the equations closely. But a small residual does not guarantee that \(\hat{\mathbf{x}}\) is close to the true solution \(\mathbf{x}\) — for ill-conditioned systems, a solution can have a small residual and a large error. This subtlety is exactly what the condition number captures, as Section 9.8 will explain.

Gaussian elimination transforms \(\mathbf{A}\) into an upper triangular matrix \(\mathbf{U}\). The multipliers used in this process can be organized into a lower triangular matrix \(\mathbf{L}\) with ones on the diagonal. The result is:

This is called the LU decomposition of \(\mathbf{A}\).

If we need to solve \(\mathbf{A}\mathbf{x} = \mathbf{b}\) for many different right-hand sides \(\mathbf{b}\) while \(\mathbf{A}\) stays the same, LU decomposition allows us to factor \(\mathbf{A}\) once and then solve two triangular systems for each new \(\mathbf{b}\):

Step 1. Forward substitution. Solve \(\mathbf{L}\mathbf{y} = \mathbf{b}\) for \(\mathbf{y}\). Since \(\mathbf{L}\) is lower triangular, this proceeds from the top row downward.

Step 2. Back substitution. Solve \(\mathbf{U}\mathbf{x} = \mathbf{y}\) for \(\mathbf{x}\). Since \(\mathbf{U}\) is upper triangular, this proceeds from the bottom row upward.

The total cost of solving for \(k\) right-hand sides after factoring is \(O(n^3) + k \cdot O(n^2)\). Without LU decomposition, each new right-hand side would require another full \(O(n^3)\) elimination. The savings for large \(k\) are substantial.

When Gaussian elimination is performed without row swaps, the multipliers \(m_{ij}\) used to eliminate entry \((i,j)\) become the entries of \(\mathbf{L}\) below the diagonal:

where \(a_{ij}^{(k)}\) denotes the entry in position \((i,j)\) after \(k\) stages of elimination. The diagonal entries of \(\mathbf{L}\) are all 1.

To measure the size of errors in vectors and matrices, numerical linear algebra uses norms. Students are likely familiar with the Euclidean norm for vectors from calculus and linear algebra. Here we introduce the concepts needed to discuss condition number and error bounds.

A vector norm \(\|\mathbf{v}\|\) assigns a nonnegative real number to each vector \(\mathbf{v}\) such that:

• \(\|\mathbf{v}\| \geq 0\), with equality only when \(\mathbf{v} = \mathbf{0}\),
• \(\|c\mathbf{v}\| = |c| \cdot \|\mathbf{v}\|\) for any scalar \(c\),
• \(\|\mathbf{u} + \mathbf{v}\| \leq \|\mathbf{u}\| + \|\mathbf{v}\|\) (triangle inequality).

The three most common vector norms in numerical work are:

In order: the 1-norm (sum of absolute values), the 2-norm (Euclidean length), and the infinity-norm (maximum absolute value).

A matrix norm \(\|\mathbf{A}\|\) measures the “size” of a matrix. The most useful matrix norms in error analysis are the induced norms, defined by:

This says: what is the most a matrix can stretch a vector? The induced infinity-norm, the one most often used in hand calculations, is the maximum absolute row sum:

For our purposes, students do not need to compute sophisticated matrix norms for every problem. The key conceptual point is that norms allow us to state rigorous bounds on errors in linear system solutions.

The condition number of a matrix is the single most important indicator of how sensitive a linear system is to perturbations in its data.

For an invertible matrix \(\mathbf{A}\), the condition number with respect to a given matrix norm is:

The condition number is always \(\geq 1\).

Suppose the right-hand side \(\mathbf{b}\) has a small perturbation \(\delta \mathbf{b}\), producing a perturbed solution \(\mathbf{x} + \delta \mathbf{x}\) where \(\mathbf{A}(\mathbf{x} + \delta \mathbf{x}) = \mathbf{b} + \delta \mathbf{b}\). Then it can be shown that:

In words: the relative change in the solution is bounded by the condition number times the relative change in the right-hand side. A condition number of \(\kappa(\mathbf{A}) = 1000\) means a relative perturbation of \(10^{-6}\) in \(\mathbf{b}\) could cause a relative change of up to \(10^{-3}\) in the solution.

More importantly for floating-point computation: since \(\mathbf{b}\) itself may be rounded to machine precision \(\epsilon_{\text{mach}} \approx 10^{-16}\) (for double precision), the solution can be expected to lose approximately \(\log_{10}(\kappa(\mathbf{A}))\) digits of accuracy. If \(\kappa(\mathbf{A}) \approx 10^{12}\), then from a 16-digit number, only about 4 reliable digits remain in the solution.

Direct methods like Gaussian elimination and LU decomposition solve \(\mathbf{A}\mathbf{x} = \mathbf{b}\) in a finite, predictable number of steps. They work well for dense systems of modest size. But for large sparse systems — where \(\mathbf{A}\) has many zero entries — direct methods are often impractical because:

• The fill-in problem: elimination introduces nonzero entries into positions that were originally zero, requiring far more memory than the sparse matrix.
• Computational cost: \(O(n^3)\) operations becomes prohibitive for \(n\) in the millions.

Iterative methods offer an alternative. Starting from an initial guess \(\mathbf{x}^{(0)}\), an iterative method produces a sequence of approximations \(\mathbf{x}^{(1)}, \mathbf{x}^{(2)}, \mathbf{x}^{(3)}, \ldots\) that converge to the true solution \(\mathbf{x}\) — if the method converges.

The key ideas are:

• Each iteration is inexpensive — often \(O(n)\) or \(O(n \cdot \text{nnz})\), where nnz is the number of nonzero entries.
• Convergence may require many iterations, but for practical tolerances the total cost can be far less than \(O(n^3)\).
• Iterative methods are well-suited for large sparse systems arising from finite difference and finite element discretizations.

The two classical iterative methods introduced here are the Jacobi method and the Gauss-Seidel method.

The Jacobi method is based on a simple idea: if we know all components of \(\mathbf{x}\) except \(x_i\), we can solve the \(i\)-th equation for \(x_i\) exactly. Since we do not know all components exactly, we substitute the current best approximations for all other components.

Write the \(i\)-th equation of \(\mathbf{A}\mathbf{x} = \mathbf{b}\) as:

Solving for \(x_i\):

The Jacobi iteration replaces \(x_j\) on the right-hand side by the current approximation \(x_j^{(k)}\):

Crucially, all components of the new approximation \(\mathbf{x}^{(k+1)}\) are computed using only values from the previous iteration \(\mathbf{x}^{(k)}\). This makes Jacobi naturally parallelizable — all components can be updated simultaneously.

The Jacobi method is guaranteed to converge if \(\mathbf{A}\) is strictly diagonally dominant: for every row \(i\),

That is, the diagonal entry in each row is larger in absolute value than the sum of all other entries in that row. Many linear systems arising from physical models satisfy this condition.

If \(\mathbf{A}\) is not diagonally dominant, Jacobi may still converge, but it may also diverge. There is no guarantee.

The Gauss-Seidel method improves on Jacobi by using the most recently computed component values immediately, rather than waiting until the end of an iteration.

The update formula for Gauss-Seidel is:

When computing \(x_i^{(k+1)}\), the method uses the already-updated values \(x_1^{(k+1)}, \ldots, x_{i-1}^{(k+1)}\) for components computed earlier in this iteration, and the previous iteration values \(x_{i+1}^{(k)}, \ldots, x_n^{(k)}\) for components not yet updated.

This uses fresh information as soon as it is available, typically improving convergence speed compared to Jacobi.

Like Jacobi, Gauss-Seidel is guaranteed to converge for strictly diagonally dominant matrices. In practice, it tends to converge roughly twice as fast as Jacobi in terms of iteration count, though it is inherently sequential (components must be updated in order) and therefore less parallelizable.

Large linear systems arise throughout applied mathematics and engineering. This section illustrates three important settings.

When engineers analyze a bridge, building, or aircraft, they use the finite element method (FEM) to divide the structure into small elements and write equilibrium equations at each node. The result is a system \(\mathbf{K}\mathbf{u} = \mathbf{f}\), where \(\mathbf{K}\) is the stiffness matrix, \(\mathbf{u}\) is the displacement vector, and \(\mathbf{f}\) is the force vector. For a realistic structure, this system might have millions of unknowns. The stiffness matrix is sparse, symmetric, and (for a stable structure) positive definite — properties that enable efficient iterative solvers.

In an electrical circuit, resistor network, or traffic flow model, each node corresponds to an unknown voltage or flow rate, and each edge corresponds to a conservation equation. Applying Kirchhoff’s current law at each node produces a system of linear equations. The coefficient matrix is a version of the Laplacian matrix of the network graph. For large networks — power grids, internet routers, supply chains — these systems are large and sparse, but well-structured.

Fitting a linear model \(y = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p\) to \(n\) observations leads to the least squares problem of minimizing \(\|\mathbf{A}\boldsymbol{\beta} - \mathbf{y}\|^2\), where \(\mathbf{A}\) is the \(n \times (p+1)\) matrix of predictor values. The normal equations \(\mathbf{A}^T \mathbf{A} \boldsymbol{\beta} = \mathbf{A}^T \mathbf{y}\) form a \((p+1) \times (p+1)\) linear system. When \(p\) is large, the condition number of \(\mathbf{A}^T \mathbf{A}\) can be very large — sometimes catastrophically so — requiring regularization or more stable solution techniques.

Numerical linear algebra is a subject where subtle errors can lead to dramatically wrong answers. The following mistakes are among the most common and consequential.

Chapter 9 has developed the numerical tools for solving linear systems: Gaussian elimination with pivoting, LU decomposition, matrix norms, condition numbers, and iterative methods. These foundations appear again in Chapter 10, which addresses a related but distinct problem: computing eigenvalues and eigenvectors numerically.

Eigenvalue problems are central in structural vibrations, stability analysis, dimensionality reduction, network analysis, and machine learning. The power method introduced in Chapter 10 is itself an iterative process, and understanding why it works requires the same thinking about convergence, error, and matrix properties that Chapter 9 has developed.

Students should carry forward from Chapter 9 several key ideas:

• Numerical linear algebra is not just symbolic row reduction done slowly on a computer. It is a discipline that must account for rounding, stability, conditioning, and efficiency.
• The condition number of a matrix determines how sensitive a problem is to perturbations — and no algorithm can overcome severe ill-conditioning.
• Direct methods (LU decomposition) are reliable and efficient for moderate-sized dense systems. Iterative methods are preferred for large sparse systems.
• Always verify solutions by checking the residual and, for sensitive problems, estimating the condition number.

Transform \([\mathbf{A} \mid \mathbf{b}]\) to upper triangular form. For column \(j\):

Then back substitute.

Solve \(\mathbf{L}\mathbf{y} = \mathbf{b}\) (forward), then \(\mathbf{U}\mathbf{x} = \mathbf{y}\) (back).