Numerics of PDE

Modeling with PDE produced a governing equation — the heat equation and its stationary form, the Laplace equation — but stopped short of solving it. For all but the simplest one-dimensional configurations there is no formula for the solution, so the temperature field has to be computed numerically. This page is about how that is done.

The first thing to settle is that the kind of PDE matters enormously. Whether an equation is elliptic, parabolic, or hyperbolic is not a cosmetic label: the three classes behave so differently that each demands its own numerical machinery, and a method built for one can fail outright on another. To keep the focus on the methods rather than on case analysis, we restrict throughout to the elliptic case, and more specifically to the simplest elliptic equation of all — the Laplace equation

Three Families of Discretization

Every numerical method for a PDE shares one opening move: replace the continuous problem, which asks for the solution at every one of infinitely many points, with a finite problem that asks for it at finitely many. What separates the major methods is how they perform that replacement, and three families dominate.

• The finite difference method (FD) approximates every derivative in the equation directly by a difference quotient. It is the most obvious approach, easy to implement, and needs little theoretical background — which is exactly why it is the one we build in full below.
• The finite volume method (FV) implements the continuum-mechanical conservation laws directly for a small control volume around each point, rather than approximating derivatives. It is the natural choice for flow problems, where conservation of mass and momentum is the physical heart of the model.
• The finite element method (FE) takes a variational route: instead of the PDE itself it enforces a slightly weakened version of it, expressed through integrals. It is more involved to implement, but it comes with a rich and powerful theory, and it is taken up in detail later on this page.

The Finite Difference Method

The grid is described by its mesh width (or grid width) $h = (h_x, h_y, h_z)$ — the spacing between neighboring points along each axis. A common simplification, used throughout the examples below, is an equidistant grid where the spacing is the same in every direction, $h = h_x = h_y = h_z$.

The derivatives are replaced one by one with the difference quotients already developed for ODEs.

A first derivative has three standard one-step approximations, differing only in which neighbors they reach for: the forward difference quotient reaches ahead of the point, its backward mirror image reaches behind it, and the symmetric central difference straddles the point. Any one of the three can stand in for the derivative — they are alternatives, not a sequence applied together:

A second derivative uses the second difference quotient, the standard three-point formula built from the point and one neighbor on each side:

The Laplace operator $\Delta u = u_{xx} + u_{yy} (+ u_{zz})$, being a sum of pure second derivatives, is approximated by combining several second difference quotients into a single formula spanning the point and its immediate neighbors.

What a stencil is

The combinations above all have the same shape: a value at a grid point is replaced by a fixed, weighted combination of the values at a small fixed pattern of nearby points. That pattern, weights and all, has a name.

The picture behind the word is a literal one: a stencil is a small template you lay over the grid, centered on the point you are currently working on, that tells you which surrounding values to take and in what proportion. Move to the next point, lay the same template down again — the pattern never changes, only where it sits.

The simplest example is the second difference quotient written above, whose stencil is the three-point stencil: the point itself and its two immediate neighbors, with weights $1, -2, 1$ (all divided by $h^2$). In two dimensions the Laplacian adds the same three-point pattern in the $y$-direction on top of the one in $x$. The two share the center point, so together they cover five points — the center and its four nearest neighbors:

This is the five-point stencil, the workhorse of the method in 2D: weight $-4$ on the center, $+1$ on each of the four neighbors. In three dimensions a third three-point pattern in $z$ joins in, giving the seven-point stencil — the center plus six neighbors, with weight $-6$ on the center. Richer stencils with more neighbors exist too, and they buy higher accuracy at a cost we return to at the end.

Setting up the equations

With the stencil chosen, the method writes one difference equation at each interior grid point: stamp the stencil down there and set the resulting weighted sum equal to the right-hand side of the PDE at that point. The right-hand side is whatever the differential operator is set equal to — the known data of the equation, with the unknown $u$ all on the left. For $\Delta u = 0$ it is simply zero; for an equation with a source term it is that prescribed function, evaluated at the grid point. The unknowns are the discrete values themselves.

Points on or near the boundary are handled by what the boundary condition prescribes there, and the two kinds split cleanly — mirroring how a boundary value problem for an ODE handles its endpoints:

• At a Dirichlet boundary the value is already given, so there is nothing to solve for: no difference equation is written and the point is not a degree of freedom.
• At a Neumann boundary the value is unknown but its normal derivative is prescribed, so the point keeps a difference equation — one of a special form that folds the boundary condition into the stencil.

Collecting one equation per interior degree of freedom turns the discretized PDE into a system of linear equations. Because each stencil reaches only a handful of neighbors, almost every entry of the system’s matrix is zero — the matrix is sparse (only a small, fixed number of nonzeros per row, regardless of how large the grid grows). Sparsity is what makes large grids tractable at all, and it is why fast iterative solution methods, rather than direct elimination, are essential — a theme that occupies much of the rest of this page.

Example: the Poisson equation on the unit square

The cleanest concrete case adds a source term to the Laplace equation.

Take it on the open unit square $\,]0, 1[^{\,2}$ (the outward-turned brackets mark an open interval — the same thing as $(0, 1)$, the endpoints $0$ and $1$ excluded; so $]0, 1[^{\,2}$ is the square’s interior without its edges) with an equidistant square grid of width $h = h_x = h_y = 1/N$. That leaves $N - 1$ interior points along each axis, so the number of degrees of freedom is

Applying the five-point stencil at every interior point produces a linear system $A x = b$. The construction is mechanical: each interior point contributes exactly one equation — its stencil, which ties the value there to its four neighbors — and that equation becomes one row of the system. Stacking all $M$ of them gives three objects:

• $x$ is the vector of unknowns: the $M$ discrete grid values themselves, the interior points collected into a single column in whatever order they are numbered.
• $A$ is the $M \times M$ coefficient matrix: row $i$ holds the stencil weights for point $i$ — the coefficient on its own value on the diagonal, and the coefficients on its neighbors in the columns those neighbors occupy in the numbering. Every other entry in the row is zero, so multiplying row $i$ by $x$ reproduces exactly the difference equation at point $i$.
• $b$ is the right-hand side: the known data, one entry per point. It carries the source term $f$ sampled at that point (so this part vanishes wherever $f$ does, and is entirely zero for pure Laplace $\Delta u = 0$), plus any known boundary values the stencil reaches, which are not unknowns and so move out of $A$ and over to the right.

The matrix $A$ is sparse — only a handful of nonzeros per row, no matter how large the grid grows — and, because the grid points are numbered row by row, has a regular band structure: its nonzeros sit on a few diagonals rather than scattered anywhere.

It is customary to multiply every equation through by $h^2$ to clear the denominators, leaving clean integer coefficients. The exact shape of the resulting rows depends on which boundary condition the problem carries:

• Pure Dirichlet. Each interior row has a $4$ on the diagonal and a $-1$ for each neighbor that is itself an interior unknown — between two and four such entries, fewer for rows whose stencil reaches a boundary point. The known boundary values that the stencil touches are moved to the right-hand side.
• Pure Neumann. Boundary points are now unknowns too, so they get their own rows; near the boundary the diagonal drops to $2$ or $3$ (a boundary point has fewer neighbors), while interior rows keep the $4$. The discretized flux condition contributes paired $(1, -1)$ entries along the boundary, with the prescribed fluxes carried to the right-hand side.

Discretization of a Dirichlet boundary

To see where the matrix rows come from, write the five-point stencil for $-\Delta u = f$ explicitly as the sum of the two three-point second differences, in $x$ and in $y$:

Multiplying by $h^2$ and reading off the coefficients gives a matrix row of the form $(\,\dots\ {-1}\ \dots\ {-1}\ \dots\ 4\ \dots\ {-1}\ \dots\ {-1}\ \dots\,)$: a $4$ on the diagonal and four $-1$‘s placed where the four neighbors sit in the global numbering. Where one of those neighbors lies on a Dirichlet boundary, its value $\psi$ is known; that term is not an unknown, so it is brought over to the right-hand side. Boundary-adjacent rows therefore carry fewer off-diagonal entries and a right-hand side enriched by the known boundary values.

That sorts the $M = (N-1)^2$ rows into three kinds. The great majority are deep-interior points, all four of whose neighbors are unknowns as well, and their rows keep the full complement of four $-1$‘s. Points sitting one step in from an edge have exactly one neighbor on the boundary, so their rows keep three; there are on the order of $4N$ of these, one band running along each side of the square. The four points tucked into the corners have two boundary neighbors apiece and so keep only two. The diagonal stays at $4$ in every row, because it comes from the stencil’s center weight, which no boundary condition touches.

Discretization of a Neumann boundary

A Neumann boundary is different in kind: the value there is not given, so the boundary point stays an unknown, but the normal derivative — the heat flow across the boundary — is prescribed. On the upper edge, for instance, that flux is a known one-sided difference,

and this relation is folded into the difference equation at the boundary point rather than replacing it. The consequence for the matrix is again a change of row shape near the edge: because a boundary point simply has fewer neighbors than an interior one, its rows carry a $3$ on the diagonal along an edge and a $2$ at a corner, instead of the interior $4$. The prescribed flux values are gathered, as always, on the right-hand side.

Accuracy and the curse of dimensionality

The standard five-point scheme is second-order accurate: the error between the computed and the exact solution shrinks quadratically with the mesh width,

Halving the grid spacing therefore cuts the error to a quarter — a good return on refinement. But refinement is not free, and in more than one dimension its cost grows alarmingly.

There are two standard ways to push back against this cost, each with its own catch:

• Higher-order stencils. Using more neighbors per stencil raises the order of accuracy — cubic, quartic, and beyond — so a coarser grid reaches the same accuracy. The price is that the matrix becomes denser: each row now couples to more points, eroding the sparsity that made the system cheap to solve.
• Locally refined (adaptive) grids. Rather than refining everywhere, place grid points densely only where the solution varies sharply and sparsely where it is smooth, economizing points overall. The difficulty is at the joints between coarse and fine regions: a fine-grid point may have no matching coarse-grid neighbor to reference, and deciding which value to use there takes extra care.

The Finite Element Method

The finite difference method attacks the derivatives head-on: every derivative in the equation is swapped for a difference quotient. The finite element method (FE) takes a route that looks, at first, nothing like that — it never approximates a derivative directly. Yet it arrives at the very same destination: a system of linear equations whose matrix is sparse, handed off to exactly the same solvers. The method is deep enough to fill courses of its own; what follows is only the rough idea, with attention to the mechanics that carry over to the finite difference systems we already built.

The whole method unfolds in five essential steps, each of which the rest of this section takes in turn:

• Substructuring and grid generation — decompose the domain into individual patches of a given pattern and of finite extent, the finite elements.
• Weak form — stop requiring the PDE to hold pointwise everywhere, and ask only that it hold in a weakened (inner-product) or averaged (integral) sense.
• Finite-dimensional ansatz space — replace the continuous solution in the weak form by a suitable finite-dimensional approximation.
• System of linear equations — use test functions to generate one equation per degree of freedom, and collect them into the associated linear system.
• Solution of the linear system — hand that system to a suitable iteration method, the subject of the remainder of this page.

Substructuring and grid generation

The decomposition step borrows its instinct straight from engineering statics: a complicated structure is understood by breaking it into standard components whose behavior is easy to describe, and then reassembling the behavior of the whole from the behavior of the parts. The finite element method does the same to a domain, cutting it into many small finite elements.

In three dimensions the result is a finite element net — a data structure with a tidy hierarchy of pieces:

• elements — the 3D atoms of the decomposition (cubes, tetrahedra, and the like),
• surfaces — their 2D faces (triangles, squares),
• edges — the 1D boundary structures of an element,
• grid points (or nodes) — the points where the unknowns live.

The choice of element shape is a trade-off: cubes are the simplest to work with, while triangles and tetrahedra are more accurate and bend far more readily around awkward geometry.

This is the sharpest contrast with finite differences. There the discretization was purely local: each grid point carried a difference formula relating it to its neighbors, and nothing tied those local relations into a single function over the whole domain. Here the trial functions assemble into one global function $u_n = \sum_k \alpha_k \varphi_k$ defined at every point of $\Omega$ — partitioned into local contributions, but a genuine function on the whole domain in the end.

The weak form

Let $L$ be the differential operator of the problem — for the Laplace equation, $L = \Delta$ (the Laplace operator). The finite element method does not ask for $Lu = f$ to hold at every point of $\Omega$. Instead it picks a finite set of test functions $\psi_l$, multiplies the equation through by each one, integrates over the whole domain, and asks only that

This is the weak form of the PDE, and the strategy is called the method of weighted residuals or the Galerkin approach.

The test functions span their own linear space, the test space $W_n$. Whether it coincides with the trial space matters enough to have a name:

• when test and trial space are the same, $W_n = V_n$, the method is the Ritz-Galerkin approach (the common case),
• when they differ, it is the Petrov-Galerkin approach.

Because $L$ is a linear operator, only the basis functions actually have to fulfill the weak form — agreement on a basis extends to every function built from it by linearity. It is customary to package the two sides of the weak form into a bilinear form $a(\cdot, \cdot)$ (the left-hand side, linear in each of its two arguments) and a linear form $b(\cdot)$ (the right-hand side), so the whole condition reads

These are $n$ linear equations, one per test function, where $n$ is the dimension of the test space.

Discrete approximation and the stiffness matrix

The weak form is $n$ equations, but $u$ is still a continuous function — infinitely many unknowns. The final move makes the unknowns finite too: replace the exact solution $u$ by a discrete approximation drawn from the trial space $V_n$, written as a linear combination of the trial functions with $n$ unknown coefficients $\alpha_k$:

Substituting $u_n$ for $u$ in the weak form and pulling the linear operator through the sum collapses everything onto the basis:

Read left to right: the bilinear form evaluated at the approximation $u_n$ equals the same form applied to its expansion $\sum_k \alpha_k \varphi_k$; linearity in the first argument lets each coefficient $\alpha_k$ come out front, leaving a weighted sum of the numbers $a(\varphi_k, \psi_l)$; and that sum must equal $b(\psi_l)$ for every test function $\psi_l$ in $W_n$. The doubly-indexed quantities $a(\varphi_k, \psi_l)$ assemble into a matrix, the singly-indexed $b(\psi_l)$ into a vector, and the unknowns are the coefficients $\alpha_k$.

The crucial observation is that every $a(\varphi_k, \psi_l)$ and every $b(\psi_l)$ depends only on the problem — the operator, the right-hand side, and the chosen basis functions — and not on the approximation of $u$ we are trying to find. They can therefore be computed once and for all, at the very beginning, producing a system of $n$ linear equations in the $n$ unknowns.

Solving the resulting system

The finite element route lands on a linear system just as the finite difference route did, and the two systems share the properties that govern how they are solved:

• With the Ritz-Galerkin approach ($V_n = W_n$), the stiffness matrix $A$ is often symmetric and positive definite (SPD) — the friendliest case for the iterative solvers below.

• The ideal constellation would be a diagonal $A$, which happens exactly when the trial and test functions are bi-orthogonal, so that

  $$a_{i,j} = a(\varphi_i, \varphi_j) = \int_\Omega L\varphi_i \cdot \varphi_j \, d\Omega = \delta_{i,j}.$$

  A diagonal system is trivial to solve, but such bases are rare and, where they exist at all, hardly accessible.

• The realistic guarantee is weaker but still decisive: because the trial and test functions have only local support, most pairs $(\varphi_k, \psi_l)$ never overlap, their integral vanishes, and $A$ comes out sparse — exactly as the five-point stencil produced a sparse matrix. With many unknowns and a sparse matrix, iterative solvers are again essential.

The strategy that follows from this is to choose trial and test spaces with good approximation properties, and then to construct bases for them that yield “nice” matrices — sparse, well-conditioned, fast to solve. The plain sample-point bases used here (one fixed basis function per grid point) are far from optimal; hierarchical bases, which resolve the solution at several scales at once, do considerably better, though their study belongs to a course of its own.

Solving the Linear System: Iterative Methods

Both routes — finite differences and finite elements — end at the same place: a large, sparse linear system $A x = b$. From here on it no longer matters where the matrix came from. The finite difference matrix is always sparse and the finite element matrix typically is too, and we simply have such a matrix in hand and want to solve it. This is one of the central problems of numerical simulation, and it shows up not only here but in the discretization of boundary value problems for ODEs as well.

The obvious idea — solve $A x = b$ exactly by a direct method like Gaussian elimination — turns out to be the wrong reflex for these systems, for three reasons:

• The number of unknowns is simply too large, especially for a PDE in 3D.
• Classical elimination fills in the matrix: the sparse structure that made the system cheap to store and apply is destroyed as elimination proceeds, and a dense system of this size is hopeless.
• Solving exactly is overkill when everything around it is already an approximation. The grid values only approximate the true solution to begin with — and in nonlinear problems a fresh linear system arises at every step of an outer iteration, so there is no point in pinning each one down to ten digits.

The alternative is an iterative method: start from a guess and improve it step by step, coming arbitrarily close to the solution without ever landing on it exactly. The aim is performance of the form “for 3 digits, one needs 10 steps” — and, crucially, a step count that stays fixed independently of the number of unknowns. The frustration that drives the rest of this page is that the classical iterative methods fall short of that aim: their speed of convergence degrades as the problem grows.

Basics of iterative methods

An iterative method begins at some starting vector $x^{(0)} \in \mathbb{R}^n$ and produces a sequence of ever-better approximations that, in the limit, reach the exact solution $x$ of $A x = b$:

The idea is that the starting point should not matter — from anywhere, the sequence homes in on $x$ — and of course we stop after finitely many steps rather than running to the limit. Throughout, a parenthesized superscript $x^{(i)}$ counts the iteration step, while a plain subscript $x_k$ picks out the $k$-th component of a vector.

How fast the sequence converges is captured by comparing successive errors:

for some $0 < \gamma < 1$ and a convergence order $s$. Larger $s$ is dramatically better: $s = 1$ is linear convergence, $s = 2$ quadratic, $s = 3$ cubic, and so on (the opposite sense to computational complexity, where a higher order is worse). The trouble is that a simple iteration for a linear system typically manages only

where $n$ is the number of grid points. With linear convergence and a contraction factor $\gamma$ this close to $1$, each step shaves off only a sliver of the error, and reaching a useful accuracy takes a great many steps.

This sets the strategy precisely. We want a method that costs only $O(n)$ arithmetic operations per step — obviously at least that much is needed, since each step must at minimum touch every unknown — while delivering a contraction factor bounded away from $1$, $\gamma < 1 - \text{const}$, with the constant independent of the grid. The two big families that chase this goal are the relaxation methods and the Krylov-subspace methods.

Relaxation methods

To design such a method, the natural thing to attack is the error of the current approximation — but the error is precisely what we do not know. The way out is to work with a closely related quantity that we can compute, the residual, and use it as a stand-in.

The chain above is worth reading slowly: the residual $b - A x^{(i)}$ measures how far the right-hand side is from being reproduced; substituting $b = A x$ (since $x$ is the exact solution) turns it into $A x - A x^{(i)} = -A\,e^{(i)}$, tying it directly to the unknown error through $A$. A small residual does not, in general, guarantee a small error — if $A$ has very large or very small entries the two can differ in size by a lot. But there is a reliable directional link between them: driving the residual down drives the error down too, roughly proportionally, so halving the residual roughly halves the error. That is enough to build a method on.

Each classical relaxation method is just a different answer to the question how do we exploit the residual to improve the approximation?

• Richardson uses the residual directly as the correction.
• Jacobi and Gauß-Seidel scale the correction so as to zero out one component of the residual at a time.
• SOR and damped methods take the same correction but deliberately overshoot or undershoot it.

The classical relaxation iterations

Richardson iteration

The simplest method takes the residual at face value and adds it straight onto the current approximation, component by component:

The two loops read just as their indices suggest: the outer loop counts the iteration step $i = 0, 1, \dots$, producing the sequence $x^{(0)}, x^{(1)}, \dots$, while the inner loop runs $k$ over the $n$ components $1, \dots, n$ of the vector, updating each one in turn. This same $i$-over-steps, $k$-over-components structure carries through all the methods below. This particular update is exactly the “residual as corrector” idea: since $r^{(i)} = -A\,e^{(i)}$ points opposite to the error, adding it nudges $x^{(i)}$ back toward $x$.

Jacobi iteration

The Jacobi iteration refines this with a normalization. Adding the raw residual is crude because the entries of $A$ set the scale: where the diagonal entry $a_{kk}$ is large the step is too big, where it is small the step is too timid. Dividing the $k$-th residual component by $a_{kk}$ corrects for this:

Using only the diagonal entry to rescale the $k$-th correction may look arbitrary, but it is well justified for the matrices we get here: they are heavily diagonally dominant (recall the five-point stencil — a $4$ on the diagonal against four $-1$‘s and a sea of zeros), so the diagonal already captures the bulk of each row, and $a_{kk}$ is the entry most responsible for component $k$.

Two details of the algorithm matter. First, the corrections $y_k$ are computed and stored in a first $k$-loop and only applied in a second $k$-loop, not immediately. The reason is subtle: if we applied $y_k$ the instant we computed it, the $k$-th component of the residual would momentarily drop to zero — equation $k$ would be solved exactly — but that gain would be undone in the very next substep $k+1$, because changing one component of $x$ changes the residual of the others (the residual was defined as $r^{(i)} = b - A x^{(i)}$, so it depends on the whole vector). Holding all corrections until the end of the step sidesteps the interference. Second, within each loop the components $k = 1, \dots, n$ do not depend on one another, so Jacobi is embarrassingly parallel — every component can be updated at once.

Gauß-Seidel iteration

The Gauß-Seidel iteration makes exactly the opposite choice: apply each correction immediately, so that when component $k$ is updated, the already-improved values of components $1, \dots, k-1$ are used in its residual. Writing the residual out with the updated entries spelled in gives

The two sums are the heart of it: the first runs over the components $j = 1, \dots, k-1$ already updated this step (so it uses $x_j^{(i+1)}$), the second over the components $j = k, \dots, n$ not yet reached (still $x_j^{(i)}$). The corrections are the same ones Jacobi computes; only the timing differs. This usually converges somewhat faster, but the immediate updates create a dependence between components and so cost the easy parallelism Jacobi enjoyed.

Damping and over-relaxation

Finally, the size of the correction can itself be tuned. Multiplying $y_k$ by a factor $\alpha$ before applying it,

gives damping when $0 < \alpha < 1$ (a deliberately cautious, too-small step) and over-relaxation when $1 < \alpha < 2$ (a deliberately bold, too-large step). Either can improve convergence in any of the three methods above. Over-relaxation applied to Gauß-Seidel is so common that it has its own name, the successive over-relaxation (SOR) method; for Jacobi, damping is the usual choice.

Convergence analysis: the additive decomposition

The algorithmic descriptions above are how the methods are run, but they are awkward to analyze. For a convergence analysis it pays to recast all four methods in a single algebraic form. They all rest on one simple idea: split the matrix as

choosing $M$ to meet two competing demands: a system $M x = b$ should be very easy to solve, and $M$ should be as close to $A$ as possible (so that the leftover $A - M$ stays small). The two pull against each other — $M = I$ makes $M x = b$ trivial but is a poor stand-in for $A$, whereas $M = A$ is a perfect stand-in but leaves $M x = b$ exactly as hard as the original problem — and a good method strikes a balance between them.

The algebraic trick is to take the trivial identity $M x + (A - M) x = b$ (which just restates $A x = b$) and split its two copies of the unknown across two successive iterates:

This is what turns an identity into an algorithm: the new approximation $x^{(i+1)}$ is computed from the old one $x^{(i)}$, and if the iteration ever settles down — $x^{(i+1)} = x^{(i)} = x$ — it reproduces the original identity exactly, so the resting point is the true solution. Solving for $x^{(i+1)}$ shows the residual reappearing:

Every relaxation method above is this one update with a particular choice of $M$. To name those choices, decompose $A$ additively into its diagonal part $D_A$, its strictly lower-triangular part $L_A$, and its strictly upper-triangular part $U_A$ (the “strictly” meaning the diagonal itself belongs to neither triangle):

The four methods then correspond to:

• Richardson: $M := I$
• Jacobi: $M := D_A$
• Gauß-Seidel: $M := D_A + L_A$
• SOR: $M := \tfrac{1}{\alpha} D_A + L_A$

The first two are immediate from the algorithms: Richardson used the residual directly as the correction, which is the update above with prefactor $M^{-1} = I$; Jacobi divided the residual by the diagonal, which is prefactor $M^{-1} = D_A^{-1}$. The Gauß-Seidel and SOR choices take a short derivation.

When the iteration converges: the spectral radius

Two consequences flow from the algebraic form $M x^{(i+1)} + (A - M) x^{(i)} = b$. The first we have already used: if the sequence $\left( x^{(i)} \right)$ converges at all, its limit is the exact solution $x$ of $A x = b$. The second is a sharp criterion for whether it converges, and it turns on a single matrix.

Whether repeated multiplication by the iteration matrix shrinks the error or lets it grow is governed by its eigenvalues — the eigenvalues and eigenvectors of $-M^{-1}(A - M)$.

The reasoning is direct from the error recurrence. Decomposed along the eigenvectors of the iteration matrix, each step scales every error component by the corresponding eigenvalue. If all eigenvalues are less than $1$ in absolute value — equivalently $\rho < 1$ — then every component shrinks each step and the error dies out. If $\rho > 1$, at least one component is amplified and the iteration diverges. The whole aim of constructing a good iterative method is therefore to make the spectral radius as small as possible, ideally near zero.

Convergence results, and why they are not enough

A handful of classical results pin down when the methods converge for the matrices that actually arise:

• A necessary condition for SOR to converge is $0 < \alpha < 2$ — over-relaxation past $2$ always diverges.
• If $A$ is positive definite, then both SOR (for $0 < \alpha < 2$) and Gauß-Seidel converge.
• If $A$ and $2 D_A - A$ are both positive definite, then Jacobi converges.
• If $A$ is strictly diagonally dominant — $a_{ii} > \sum_{j \neq i} |a_{ij}|$ for every row $i$ — then both Jacobi and Gauß-Seidel converge.
• In certain cases the optimal relaxation parameter $\alpha$ can be computed exactly, the one that minimizes $\rho$ and so maximizes the error reduction per step.

One natural guess turns out to be false: applying corrections immediately does not make Gauß-Seidel uniformly better than Jacobi. There are matrices on which Gauß-Seidel converges while Jacobi diverges, and others the other way around. In many practical cases, though, Gauß-Seidel reaches a given accuracy in about half as many steps.

The spectral radius decides not just whether the iteration converges but how fast: the smaller $\rho$, the more sharply every error component is cut each step. And here is where the classical methods disappoint. In practice $\rho$ often sits so close to $1$ that — convergence notwithstanding — the number of steps needed for acceptable accuracy is far too large. The discretized PDE is the cautionary example: its spectral radius depends on the problem size $n$ and hence on the grid resolution $h$, typically as

This is the worst possible coupling: the finer the grid — and so the more accurate the solution we are after — the closer $\rho$ creeps to $1$ and the more miserably the iteration crawls. Refining the grid to gain accuracy directly sabotages the solver meant to exploit it. Escaping this trap — building a solver whose convergence does not decay as the grid is refined — is the concern that drives the rest of this page.

Minimization Methods: Steepest Descent and Conjugate Gradients

The relaxation methods were the first of the two families named earlier. The second comes at the linear system from a completely different angle, and it begins with a change of question. Rather than solving $A x = b$ head-on, it steps up one level to an artificial minimization problem, solves that instead, and lands back on the solution of the linear system as a by-product.

Solving a linear system as minimizing a quadratic

The trick rests on a different way of looking at a positive definite matrix $A$ (one whose eigenvalues are all strictly positive). Attach to the system the scalar-valued quadratic function

where $A$ is the system matrix, $b$ the right-hand side, $c$ an arbitrary scalar constant, and the products are arranged so the whole expression returns a single number. This $f$ is artificial — nothing in the original problem asked for it; it is introduced purely as a solution device. Its point is the equivalence

and because $A$ is positive definite, that minimum is unique.

The method of steepest descent

Once solving is recast as finding the bottom of a bowl, the most natural descent strategy applies. From the gradient we know that $\nabla f$ points in the direction of steepest ascent and is orthogonal to the contour lines of $f$. To go downhill as fast as possible, then, step in the opposite direction — the negative gradient. And the negative gradient is an old friend: from $\nabla f(x) = A x - b$,

the residual. So at each step the search direction is simply the residual.

A direction alone is not enough; we also need a step size — how far to travel along it. The natural choice is the value that goes as low as possible along that one ray, a one-dimensional line search

and because $f$ is quadratic this minimization along a line has a closed-form answer. Together the direction and the optimal step size make up the method of steepest descent.

The three lines are the whole method. The first sets the step size to the line-search minimum. The second moves the approximation that far along the residual. The third updates the residual — and it is worth seeing that this is not a fresh computation but a cheap reuse: since $r^{(i+1)} = b - A x^{(i+1)} = b - A\bigl(x^{(i)} + \alpha_i r^{(i)}\bigr) = r^{(i)} - \alpha_i A r^{(i)}$, and the matrix-vector product $A r^{(i)}$ was already formed for the denominator of $\alpha_i$, the residual rides along for almost no extra work. The one strict requirement is that $A$ be positive definite, so that $f$ has a minimum to descend toward at all. It is a reasonable procedure, but a very heuristic one.

Why steepest descent is slow: the condition number

Before improving the method it pays to ask how it can possibly go wrong — and can it get any easier? One simplification is to search not along the gradient but along the coordinate axes in turn, one component at a time. Doing exactly that turns out to reproduce the Gauß-Seidel iteration from the relaxation family — a first sign that these methods, approached from opposite directions, are all related.

The trouble with steepest descent is its speed: it can be arbitrarily slow. The reason is that each step optimizes locally and in isolation, with no memory of the steps before it, so a new step routinely destroys part of what an earlier one had already achieved. On a long, narrow bowl the iterates zig-zag back and forth across the valley in tiny right-angle steps, creeping toward the minimum (this is the oscillation the line search cannot avoid: each direction is locally best but globally incoherent with the others).

How elongated that bowl is — and therefore how badly the method zig-zags — is measured by a single number built from the eigenvalues of $A$.

The condition number is not a quantity the algorithm computes or uses — it never appears in the update rules above. It is a diagnostic: it explains and predicts how fast the iteration will converge, and that is all. What it tells us is what to aim for. We want $\kappa(A)$ as small as possible, and — crucially — kept small even as the grid is refined, so the spectrum of $A$ does not spread out without bound as $n$ grows. Achieving that means designing algorithms that tame the eigenvalues: this time the eigenvalues of the system matrix $A$ itself, as opposed to the eigenvalues of the iteration matrix (the spectral radius) that governed the relaxation methods. The two cures developed below — conjugate directions and preconditioning — are two ways of doing exactly that.

Conjugate directions

The flaw to fix is the destruction: steepest descent keeps spoiling its own past progress. The cure is to choose search directions that never interfere — so that whatever a step accomplishes along its direction is never undone by a later step. Concretely, after $i$ steps the remaining error should be orthogonal to all the previous search directions, so that no earlier gain can leak back.

If the directions are chosen this way and the error component along each is eliminated exactly once, then in $\mathbb{R}^n$ the process must finish in at most $n$ steps — every direction handled, no error left. That makes it, in principle, a direct method: after $n$ steps it lands exactly on the optimum, no approximation about it.

The catch is that $n$ is enormous for a fine grid — millions or billions of unknowns in 3D — so running all $n$ steps is out of the question. In practice the method is stopped far earlier, the moment the approximation is accurate enough, and used exactly like an iterative method. A scheme that is a direct method in exact theory but is run as an iteration in practice is called semi-iterative: for a system with $n = 10^{10}$ unknowns it would reach the exact answer after $10^{10}$ steps, but no one ever takes $10^{10}$ steps — a handful suffices for the accuracy a simulation actually needs.

The update now steps along a dedicated search direction $d^{(i)}$ rather than the raw residual,

and the directions must be built so the no-interference property holds. The ideal condition would be that each direction is orthogonal to the new error, $0 = d^{(i)T} e^{(i+1)}$ — but the error $e^{(i+1)}$ is exactly the unknown we lack, so this cannot be tested directly. The way out is to replace plain orthogonality with orthogonality through the matrix $A$:

Conjugacy is the condition that can be enforced, precisely because $A\, e^{(i+1)} = -r^{(i+1)}$ is the (computable) residual — so the unknown error is once again accessed only through $A$, exactly as the residual gave us a handle on it before. The iteration starts from $d^{(0)} = r^{(0)}$, the steepest-descent direction. What remains is the recipe for producing the later conjugate directions $d^{(i)}$.

Conjugate gradients

The final method adds that recipe to the conjugate-direction iteration. The construction principle is Gram-Schmidt conjugation of the residuals — the standard orthogonalization procedure, but performed in the $A$-orthogonality sense to turn the successive residuals into conjugate directions. Carried out naively, Gram-Schmidt would be ruinously expensive: each new direction would have to be made conjugate against all its predecessors, and accumulating those corrections over $n$ steps costs on the order of $n$ work per step, hence $O(n^3)$ overall — no better than the direct elimination we were trying to avoid.

The remarkable discovery at the heart of the method is that, for this particular setup, almost all of that work is unnecessary: the conjugation against every earlier direction collapses to a correction against only the single previous one. Enforcing conjugacy against the last direction turns out to deliver conjugacy against all of them for free. That short recurrence is what makes the method practical, and it is the conjugate gradient method.

The iteration starts, as in the conjugate-direction scheme above, from $d^{(0)} = r^{(0)} = b - A x^{(0)}$ — the first search direction is just the residual at the starting guess $x^{(0)}$, the steepest-descent direction — after which the recurrence in the last line generates every subsequent $d^{(i+1)}$. The five lines of the loop body then read as follows:

• $\alpha_i := \dfrac{d^{(i)T} r^{(i)}}{d^{(i)T} A d^{(i)}}$ — the line-search step size, exactly as in steepest descent, but measured along the search direction $d^{(i)}$ instead of along the residual.
• $x^{(i+1)} := x^{(i)} + \alpha_i\, d^{(i)}$ — step the current approximation that far along the search direction.
• $r^{(i+1)} := r^{(i)} - \alpha_i\, A d^{(i)}$ — update the residual cheaply, reusing the matrix-vector product $A d^{(i)}$ already formed for $\alpha_i$ rather than recomputing $b - A x^{(i+1)}$ from scratch.
• $\beta_{i+1} := \dfrac{r^{(i+1)T} r^{(i+1)}}{r^{(i)T} r^{(i)}}$ — the conjugation coefficient: the ratio of the new residual’s squared length to the old one’s, a measure of how much the residual has shrunk this step.
• $d^{(i+1)} := r^{(i+1)} + \beta_{i+1}\, d^{(i)}$ — the new search direction: the fresh residual nudged by $\beta_{i+1}$ times the old direction.

The first three lines are steepest descent with $d^{(i)}$ in place of $r^{(i)}$; the last two are what is genuinely new. That single backward nudge by $\beta_{i+1} d^{(i)}$ is the whole of the cheap conjugation — it keeps the directions mutually conjugate while storing only the one previous direction, never the full history. CG is markedly faster than steepest descent — but its step count still depends on $n$.

There is a deeper structure to the directions it generates. The spaces spanned by the successive search directions form a so-called Krylov sequence:

Each iteration enlarges the search space by one more power of $A$ applied to the initial residual — and this is the structure that gives the second family the name promised earlier, the Krylov-subspace methods. CG is the most prominent member for symmetric positive definite systems; other widely used Krylov methods, built to handle matrices that are not symmetric positive definite, are GMRES and Bi-CGSTAB.

Preconditioning

Here the two families meet the same wall. The goal set out at the very start was a solver whose step count stays fixed no matter how many unknowns there are — “for 3 digits, ten steps,” independent of the grid resolution. Neither family reaches it. The relaxation methods stall because their spectral radius $\rho$ creeps toward $1$ as the grid is refined, so each step removes less and less of the error; CG is much faster, but its step count still grows with $n$. The decisive disadvantage they share is therefore exactly the one we were trying to avoid: the finer the grid, the more iterations the solver needs — convergence slows down precisely as the problem grows, which is precisely when we most want it to stay cheap.

This is not a small annoyance. Refining the grid is the whole point of the exercise — it is how the discrete solution is made accurate — so a solver that bogs down under refinement undercuts the very thing it is there to support. Escaping the wall is what the remainder of the page is about, and the cure splits along family lines:

• For the relaxation methods, the cure is an explicit application of the multigrid principle (the subject of the rest of this page).
• For the Krylov/CG methods, the cure is preconditioning — and the preconditioner is itself often multigrid-based.

Take preconditioning first. The quantity that governs CG’s convergence is the condition number $\kappa(A)$ of the matrix, and for a discretized PDE that condition number increases dramatically as $n$ grows — roughly like the inverse square of the mesh width, $\kappa(A) = O(h^{-2})$. The idea is blunt: if the matrix is the problem, change the matrix. Multiplying $A x = b$ from the left by some invertible matrix $M^{-1}$ leaves the solution $x$ untouched but replaces the system with an equivalent one,

and if $M$ is chosen so that $M^{-1} A$ has a far smaller condition number than $A$ — ideally close to $1$ — the minimization methods converge dramatically faster on the new system. The symbol $M^{-1}$ is just notation; the inverse is never formed explicitly, only applied.

There is a catch: $M^{-1} A$ need not be symmetric even when $A$ is, and CG relies on symmetric positive definiteness. The repair is to factor the (symmetric positive definite) preconditioner as $M = W W^T$ and transform symmetrically, which yields the fully expanded chain

where $M$ is symmetric positive definite, $W W^T = M$, and $y = W^T x$. The two preconditioned matrices $M^{-1} A$ and $W^{-1} A W^{-T}$ are similar (they share the same eigenvalues, so the conditioning improvement is identical), but the symmetric form $W^{-1} A W^{-T}$ has the advantage of preserving the symmetry of $A$ — keeping the system in the symmetric positive definite class that CG needs.

The two demands on $M$ pull in opposite directions in exactly the way the additive splitting $A = M + (A - M)$ did for the relaxation methods: the easier $M$ is to apply, the less it resembles $A$, and vice versa. The art is in the compromise.

Strategies for preconditioners

The choice of $M$ spans a spectrum between the two useless extremes:

• $M = I$ — the easiest possible choice and the cheapest to apply, but useless: it changes nothing, leaving the original unpreconditioned method.
• $M = A$ — the perfect choice for conditioning, since then $M^{-1} A = I$ exactly, but worthless in practice: applying $M^{-1} = A^{-1}$ means solving a system as hard as the original $A x = b$, so computing the inverse is even more expensive than the problem we started from.

Every practical preconditioner is a compromise between these poles, cheap to apply yet a decent stand-in for $A$:

• Diagonal (Jacobi) preconditioner — take $M = D_A$, the diagonal of $A$ (reusing the $A = L_A + D_A + U_A$ split from the relaxation analysis). It is the cheapest nontrivial choice and trivial to invert.

• No Gauß-Seidel or SOR — these are deliberately not used as preconditioners, because the matrix $M$ they would supply is not symmetric, which would break the symmetric positive definite structure CG depends on.

• SSOR preconditioner — the symmetric variant of SOR restores symmetry by combining a forward and a backward sweep. With the over-relaxation parameter $\alpha$ (the same parameter as in SOR, not the line-search step size $\alpha_i$ above), it is built from

  $$\begin{aligned} M^{(1/2)} &:= \alpha^{-1} D_A + L_A, \\ M^{(1)} &:= \alpha^{-1} D_A + U_A, \\ M &:= \frac{\alpha}{\alpha - 2}\,\bigl(M^{(1/2)}\bigr)^{-1} D_A^{-1} M^{(1)}. \end{aligned}$$

• Imperfect (incomplete) factorization — for example ILU (incomplete LU): compute only approximate triangular factors $L$ and $U$ instead of the exact ones a direct solver would use, cheap to form and sparse, and use their product $M = L U$ as the preconditioner. These are not the strictly-triangular parts $L_A$ and $U_A$ from the relaxation analysis, despite the similar names. That split was additive, $A = L_A + D_A + U_A$, with the diagonal held separately in $D_A$, so both triangles were strict. This one is multiplicative, $A \approx L U$, and each factor carries a diagonal of its own (by convention ones along $L$, the substantive diagonal in $U$).

• Sparse approximate inverse — construct a simple sparse matrix $B$ that approximates $A^{-1}$ directly, by minimizing the residual of the inverse relation

  $$\min_B\ \lVert I - A B \rVert^2, \qquad M^{-1} = B,$$

  so that applying the preconditioner is just a sparse multiplication by $B$.

• Multilevel (multigrid) preconditioner — build $M$ to follow the multigrid principle, the deepest of the cures and the subject of the rest of this page.