43.11.01 · numerical-analysis / 11-finite-difference-pdes

Finite differences for the elliptic BVP: the 5-point Laplacian and its convergence

shipped3 tiersLean: none

Anchor (Master): LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) Ch. 3 (the elliptic discretisation, the M-matrix structure, the discrete Green's function, and the $O(h^2)$ convergence theorem); Hackbusch 1992 *Elliptic Differential Equations: Theory and Numerical Treatment* (Springer) §4.1-4.7 (the discrete maximum principle, M-matrices, the inverse-monotone stability bound, and the explicit eigensystem of the discrete Laplacian); Varga 2000 *Matrix Iterative Analysis* 2e (Springer) Ch. 1, 3 (M-matrices, monotonicity, and the diagonally dominant irreducible structure of $A_h$)

Intuition Beginner

Many physical balance laws say the same plain thing: at every interior point, a quantity sits at the average of its neighbours, nudged by whatever source is pushing on it. The steady temperature in a metal plate, the shape of a stretched drum at rest, the electric potential between charged plates — all obey a rule of this kind. The continuous form of that rule is the Poisson equation from 02.13.02, and on its own it lives at infinitely many points at once.

A computer cannot hold infinitely many points. So we lay a grid over the region, like graph paper, and keep only the values at the crossing points. The trick is to replace the smooth bending of the solution by a simple comparison among grid neighbours. In one direction, the bending at a point is measured by how much the point sits below the average of its left and right neighbours. In two directions, you compare a point against its four neighbours: left, right, up, and down. That four-neighbour comparison is the famous five-point stencil, and it turns the equation into a finite list of linear relations you can actually solve.

Once you write down one such relation at every interior grid point, you have a large system of linear equations. Solving it gives a number at each grid point that approximates the true solution there. The system is enormous but mostly empty: each equation touches only five values, so almost every entry is zero. That sparseness is what makes the method practical, and the solve itself uses the elimination machinery of 43.03.01.

Two questions then matter. Does the grid answer get closer to the true answer as the grid gets finer? And can the method ever produce a wild value that the true solution would never reach? The healthy answer to both is yes-it-converges and no-it-cannot, and the reason is a discrete echo of a fact about the continuous problem: a source-free solution attains its extremes only on the boundary.

Visual Beginner

The picture is a single interior grid point with its four neighbours, and the weights the method attaches to each. The centre point is compared against the average of the four points around it.

The stencil places a weight of four on the centre and a weight of minus one on each of the four neighbours, all divided by the square of the grid spacing. When the source is zero, this forces the centre value to be exactly the average of its neighbours.

            (i, j+1)
               |
              -1
               |
(i-1,j)--  -1  +4  -1  --(i+1,j)
               |
              -1
               |
            (i, j-1)

position	weight in the stencil	meaning
centre point	$+ 4$	the value being solved for
each of the four neighbours	$- 1$	pulls the centre toward the neighbour average
anything else	$0$	the stencil is local, so the system is sparse

The takeaway: the entire method is this one small picture, copied to every interior grid point. The boundary dots are not unknowns — their values are handed to you by the boundary data — so the unknowns are only the interior dots, and each one is tied to four neighbours.

Worked example Beginner

Let us solve the smallest interesting case by hand. Take the unit square and put down a grid with spacing $h = 1/3$ , so the interior has just four unknown points arranged in a two-by-two block. We solve $- Δ u = f$ with the source set to zero, $f = 0$ , and with boundary values all equal to $0$ except along the top edge, where the boundary value is $1$ .

Label the four interior unknowns $U_{1}, U_{2}$ (bottom row, left then right) and $U_{3}, U_{4}$ (top row, left then right). With $f = 0$ , the five-point rule says each interior value equals the average of its four neighbours.

Step 1. Write the rule at $U_{3}$ (top-left interior point). Its neighbours are: left boundary $0$ , the point $U_{4}$ to its right, the top boundary $1$ above, and $U_{1}$ below. So $U_{3} = \frac{1}{4} (0 + U_{4} + 1 + U_{1})$ .

Step 2. Write the rule at $U_{4}$ (top-right interior point): left neighbour $U_{3}$ , right boundary $0$ , top boundary $1$ , below $U_{2}$ . So $U_{4} = \frac{1}{4} (U_{3} + 0 + 1 + U_{2})$ .

Step 3. Write the rule at $U_{1}$ (bottom-left): left boundary $0$ , right $U_{2}$ , above $U_{3}$ , bottom boundary $0$ . So $U_{1} = \frac{1}{4} (0 + U_{2} + U_{3} + 0)$ . Likewise $U_{2} = \frac{1}{4} (U_{1} + 0 + U_{4} + 0)$ .

Step 4. By the left-right symmetry of the data, $U_{1} = U_{2}$ and $U_{3} = U_{4}$ . Substitute. From step 3, $U_{1} = \frac{1}{4} (U_{1} + U_{3})$ , giving $4 U_{1} = U_{1} + U_{3}$ , so $U_{3} = 3 U_{1}$ . From step 1, $U_{3} = \frac{1}{4} (U_{3} + 1 + U_{1})$ , giving $4 U_{3} = U_{3} + 1 + U_{1}$ , so $3 U_{3} = 1 + U_{1}$ . Put $U_{3} = 3 U_{1}$ into this: $9 U_{1} = 1 + U_{1}$ , so $U_{1} = 1/8$ and $U_{3} = 3/8$ .

What this tells us: the four grid values are $U_{1} = U_{2} = 1/8$ and $U_{3} = U_{4} = 3/8$ . Every value sits between $0$ and $1$ , the smallest and largest boundary values — the grid solution never overshoots the data. The points nearer the hot top edge are warmer, exactly as physical intuition demands. This non-overshooting is the discrete maximum principle at work.

Check your understanding Beginner

Formal definition Intermediate+

Let $Ω = (0, 1)^{2}$ be the open unit square with boundary $\partial Ω$ , and consider the Dirichlet problem for the Poisson equation 02.13.02 $$ -\Delta u = f \ \text{ in } \Omega, \qquad u = g \ \text{ on } \partial\Omega, $$ with $f \in C (\overline{Ω})$ and $g \in C (\partial Ω)$ , whose classical solvability and regularity are assumed as input from the continuous theory and are not reproved here. Fix an integer $m \geq 1$ , set the mesh width $h = 1/ (m + 1)$ , and place the uniform grid points $(x_{i}, y_{j}) = (ih, j h)$ for $0 \leq i, j \leq m + 1$ . The interior grid indices are $1 \leq i, j \leq m$ , numbering $N = m^{2}$ unknowns $U_{ij} \approx u (x_{i}, y_{j})$ ; the boundary grid values $U_{ij}$ for $i$ or $j$ in ${0, m + 1}$ are fixed at $g$ .

The 1-D second-difference operator. For a grid function $V$ in one variable, the second-difference operator is $$ (D^2 V)i = \frac{V{i-1} - 2V_i + V_{i+1}}{h^2}. $$ On the $m$ interior points it is represented by the symmetric tridiagonal matrix $T = \frac{1}{h ^{2}} tridiag (- 1, 2, - 1) \in R^{m \times m}$ , so that $- D^{2}$ has matrix $T$ .

The 5-point Laplacian. The discrete Laplacian acts on a grid function $U$ by $$ (\Delta_h U){ij} = \frac{U{i-1,j} + U_{i+1,j} + U_{i,j-1} + U_{i,j+1} - 4U_{ij}}{h^2}, $$ the sum of the second differences in $x$ and in $y$ . The discretised problem is $- Δ_{h} U_{ij} = f_{ij} := f (x_{i}, y_{j})$ at every interior point, with the boundary values substituted from $g$ . Collecting the interior unknowns into a vector $U \in R^{N}$ under a fixed (e.g. lexicographic) ordering, this is the sparse linear system $$ A_h U = F, $$ where $A_{h} \in R^{N \times N}$ represents $- Δ_{h}$ and $F$ collects $f_{ij}$ together with the boundary contributions $g / h^{2}$ moved to the right-hand side. In Kronecker form $A_{h} = I_{m} \otimes T + T \otimes I_{m}$ , exhibiting $A_{h}$ as built from the 1-D operator $T$ in each coordinate. The matrix $A_{h}$ is symmetric positive definite, has $+ 4/ h^{2}$ on the diagonal and $- 1/ h^{2}$ in the at-most-four off-diagonal positions of each row, and is sparse with at most five nonzeros per row; the system is solved by the LU/sparse-elimination machinery of 43.03.01. The symbols $Δ_{h}$ , $D^{2}$ , $T$ , $A_{h}$ , the Kronecker product $\otimes$ , the discrete norms $∥ \cdot ∥_{\infty}$ and $∥ \cdot ∥_{h}$ , and the Landau symbol $O (h^{2})$ are recorded in _meta/NOTATION.md.

Local truncation error. The local truncation error of the scheme is the residual obtained by substituting the exact solution into the discrete operator, $$ \tau_{ij} = -\Delta_h u(x_i, y_j) - f(x_i, y_j) = -\Delta_h u(x_i, y_j) + \Delta u(x_i, y_j), $$ measuring the failure of the exact solution to satisfy the discrete equation. The scheme is consistent of order $p$ if $∥ τ ∥_{\infty} = O (h^{p})$ .

M-matrix and monotonicity. A real square matrix is an M-matrix if it has nonpositive off-diagonal entries and a nonnegative inverse ( $A_{h}^{- 1} \geq 0$ entrywise); equivalently here, $A_{h}$ is monotone — $A_{h} V \geq 0$ implies $V \geq 0$ entrywise. The sign pattern of $A_{h}$ (positive diagonal, nonpositive off-diagonals) together with weak diagonal dominance and irreducibility delivers this structure, which is the algebraic carrier of the discrete maximum principle ^{[Varga, R. S. — Matrix Iterative Analysis (2nd ed.)]}.

Counterexamples to common slips

Consistency is not convergence. The truncation error $τ$ measures how well the exact solution satisfies the discrete equation; it says nothing on its own about how far the computed $U$ lies from $u$ . Bridging the two requires the stability bound on $A_{h}^{- 1}$ . A scheme can be consistent and still diverge if its discrete operator is not stable.
The five-point Laplacian is $O (h^{2})$ , not $O (h^{4})$ . The Taylor expansion of $D^{2}$ carries a leading error term $\frac{1}{12} h^{2} u^{''''}$ , which does not cancel between the two coordinate directions for a general $f$ . Believing the centred difference is fourth-order confuses the symmetry of the stencil (which kills the odd-order terms) with the vanishing of the leading even term (which it does not achieve).
Symmetric positive definiteness alone does not give the maximum principle. The discrete maximum principle rests on the sign pattern and monotonicity ( $A_{h}^{- 1} \geq 0$ ), an order-theoretic property, not merely on positive definiteness. A symmetric positive definite matrix with a positive off-diagonal entry need not be monotone and need not obey a discrete maximum principle.
The boundary values are data, not unknowns. Only the $N = m^{2}$ interior values are solved for; the boundary grid values enter $A_{h} U = F$ through the right-hand side $F$ . Treating boundary points as unknowns mis-sizes the system and breaks the M-matrix structure.

Key theorem with proof Intermediate+

The signature result is that the five-point scheme converges at second order in the maximum norm, and that this follows from two separately checkable facts: a consistency estimate (the truncation error is $O (h^{2})$ ) and a stability estimate (the discrete solution operator $A_{h}^{- 1}$ is bounded in the maximum norm uniformly in $h$ ). The stability estimate is itself a consequence of the discrete maximum principle.

Theorem (consistency + stability $\Rightarrow$ $O (h^{2})$ convergence). Let $u \in C^{4} (\overline{Ω})$ solve $- Δ u = f$ in $Ω = (0, 1)^{2}$ with $u = g$ on $\partial Ω$ , and let $U$ solve the five-point system $A_{h} U = F$ . Then there is a constant $C$ , independent of $h$ , with $$ |U - u|{\infty} := \max{1\le i,j\le m} |U_{ij} - u(x_i, y_j)| \le C,h^2 \max_{\overline\Omega}\big(|u_{xxxx}| + |u_{yyyy}|\big). $$ ^{[LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems; Hackbusch, W. — Elliptic Differential Equations: Theory and Numerical Treatment]}

Proof. Let $E_{ij} = U_{ij} - u (x_{i}, y_{j})$ denote the global error grid function, vanishing on the boundary. Apply the discrete operator to $E$ . Since $A_{h} U = F$ encodes $- Δ_{h} U_{ij} = f_{ij}$ , and the exact solution satisfies $- Δ_{h} u (x_{i}, y_{j}) = f_{ij} + τ_{ij}$ by the definition of the local truncation error, subtraction gives $- Δ_{h} E_{ij} = - τ_{ij}$ at every interior point, that is $A_{h} E = - τ$ with $τ$ the vector of truncation values and $E$ vanishing on the boundary.

Consistency. Expand $u$ in a Taylor series about $(x_{i}, y_{j})$ . The centred second difference in $x$ gives $$ \frac{u(x_{i-1},y_j) - 2u(x_i,y_j) + u(x_{i+1},y_j)}{h^2} = u_{xx}(x_i,y_j) + \frac{h^2}{12}u_{xxxx}(\xi_{ij}, y_j), $$ for some $ξ_{ij}$ between $x_{i - 1}$ and $x_{i + 1}$ , the odd-order terms cancelling by symmetry; the analogous expansion holds in $y$ . Adding the two and using $Δ u = - f$ , $$ \Delta_h u(x_i,y_j) = \Delta u(x_i,y_j) + \frac{h^2}{12}\big(u_{xxxx} + u_{yyyy}\big) + O(h^4), $$ so $τ_{ij} = - Δ_{h} u (x_{i}, y_{j}) + Δ u (x_{i}, y_{j}) = - \frac{h ^{2}}{12} (u_{xxxx} + u_{y y y y}) + O (h^{4})$ and hence $∥ τ ∥_{\infty} \leq \frac{h ^{2}}{12} max_{\overline{Ω}} (∣ u_{xxxx} ∣ + ∣ u_{y y y y} ∣) + O (h^{4})$ .

Stability. The claim is the uniform bound $∥ A_{h}^{- 1} ∥_{\infty} \leq \frac{1}{8}$ . Introduce the comparison (barrier) grid function $ϕ_{ij} = \frac{1}{4} ((x_{i} - \frac{1}{2})^{2} + (y_{j} - \frac{1}{2})^{2})$ , the restriction of a quadratic with $- Δ ϕ = - 1$ , hence $- Δ_{h} ϕ_{ij} = 1$ exactly at every interior point (the second difference of a quadratic is its exact second derivative). On the square, $0 \leq ϕ_{ij} \leq \frac{1}{8}$ . Given any interior grid function $W$ with $A_{h} W = R$ and $W = 0$ on the boundary, set $W^{\pm} = \pm W - ∥ R ∥_{\infty} ϕ$ . Then $A_{h} W^{\pm} = \pm R - ∥ R ∥_{\infty} \cdot 1 \leq 0$ at every interior point, since each component of $\pm R$ is at most $∥ R ∥_{\infty}$ , and on the boundary $W^{\pm} = - ∥ R ∥_{\infty} ϕ \leq 0$ . The discrete maximum principle (Proposition 1 below) applied to $- W^{\pm}$ — for which $A_{h} (- W^{\pm}) \geq 0$ — gives $W^{\pm} \leq 0$ throughout, that is $\pm W_{ij} \leq ∥ R ∥_{\infty} ϕ_{ij} \leq \frac{1}{8} ∥ R ∥_{\infty}$ . Hence $∥ W ∥_{\infty} \leq \frac{1}{8} ∥ A_{h} W ∥_{\infty}$ , which is exactly $∥ A_{h}^{- 1} ∥_{\infty} \leq \frac{1}{8}$ .

Synthesis. Apply the stability bound to $E$ , which satisfies $A_{h} E = - τ$ and vanishes on the boundary: $$ |E|\infty \le \tfrac18|A_h E|\infty = \tfrac18|\tau|\infty \le \frac{h^2}{96}\max{\overline\Omega}\big(|u_{xxxx}| + |u_{yyyy}|\big) + O(h^4). $$ Taking $C = 1/96$ (absorbing the higher-order remainder into $C$ for $h$ below a fixed bound) gives the stated estimate. $□$

Bridge. This proof is the foundational reason the five-point method is trusted: it factors the error into a consistency piece, controlled by Taylor expansion of the stencil, and a stability piece, controlled by the maximum-norm bound on $A_{h}^{- 1}$ , and their product is the convergence rate. This is exactly the consistency-plus-stability template that builds toward the Lax-Richtmyer equivalence theorem 43.11.05, where the bounded discrete operator is the time-evolution operator $C (Δ t)$ rather than the steady-state inverse $A_{h}^{- 1}$ ; it generalises the discrete-operator stability idea of 43.03.01, where the role here played by $A_{h}^{- 1}$ was played by the backward error of elimination. The stability half is dual to the consistency half — one bounds how the scheme amplifies residuals, the other bounds the residual the exact solution leaves behind — and putting these together is the central insight that a consistent, stable discretisation converges at the consistency rate. The same pattern appears again in the parabolic theory of 43.11.02, where the very same matrix $A_{h}$ reappears as the spatial operator of the method of lines and its eigenvalues set the time-step restriction.

Exercises Intermediate+

Exercise 3 (medium, symbolic).

Show by Taylor expansion that the 1-D second-difference operator $D^{2} V_{i} = (V_{i - 1} - 2 V_{i} + V_{i + 1}) / h^{2}$ applied to a smooth $u$ satisfies $D^{2} u (x_{i}) = u^{''} (x_{i}) + \frac{h ^{2}}{12} u^{''''} (x_{i}) + O (h^{4})$ .

Hint

Expand $u (x_{i} \pm h)$ to fourth order and add; the odd-order terms cancel.

Answer

Expand $u (x_{i} + h) = u + h u^{'} + \frac{h ^{2}}{2} u^{''} + \frac{h ^{3}}{6} u^{'''} + \frac{h ^{4}}{24} u^{''''} + O (h^{5})$ and $u (x_{i} - h) = u - h u^{'} + \frac{h ^{2}}{2} u^{''} - \frac{h ^{3}}{6} u^{'''} + \frac{h ^{4}}{24} u^{''''} + O (h^{5})$ , evaluating all derivatives at $x_{i}$ . Adding, the odd-order terms cancel: $u (x_{i} + h) + u (x_{i} - h) = 2 u + h^{2} u^{''} + \frac{h ^{4}}{12} u^{''''} + O (h^{6})$ . Subtract $2 u$ and divide by $h^{2}$ : $D^{2} u (x_{i}) = u^{''} + \frac{h ^{2}}{12} u^{''''} + O (h^{4})$ . The leading error is the $\frac{h ^{2}}{12} u^{''''}$ term, so the operator is second-order accurate.

Exercise 5 (medium, symbolic).

Verify that the grid function $ϕ_{ij} = \frac{1}{4} ((x_{i} - \frac{1}{2})^{2} + (y_{j} - \frac{1}{2})^{2})$ satisfies $- Δ_{h} ϕ_{ij} = 1$ exactly at every interior point, and bound $0 \leq ϕ_{ij} \leq \frac{1}{8}$ on the unit square.

Hint

The second difference of a quadratic is exact: $D^{2} (x^{2}) = 2$ . Use this in each coordinate.

Answer

Write $ϕ = \frac{1}{4} (p (x) + q (y))$ with $p (x) = (x - \frac{1}{2})^{2}$ , $q (y) = (y - \frac{1}{2})^{2}$ . For a quadratic, the centred second difference returns the exact second derivative: $D^{2} p (x_{i}) = p^{''} (x_{i}) = 2$ , since the third- and fourth-order Taylor terms vanish. Thus $Δ_{h} ϕ_{ij} = \frac{1}{4} (D^{2} p + D^{2} q) = \frac{1}{4} (2 + 2) = 1$ , so $- Δ_{h} ϕ_{ij} = - 1$ ; the barrier used in the proof takes $- Δ_{h} ϕ = 1$ for $ϕ = - \frac{1}{4} (\dots)$ , or equivalently one uses this $ϕ$ with the sign tracked. On $[0, 1]^{2}$ each of $(x - \frac{1}{2})^{2}, (y - \frac{1}{2})^{2}$ lies in $[0, \frac{1}{4}]$ , so $ϕ = \frac{1}{4} (\cdot + \cdot) \in [0, \frac{1}{8}]$ . The maximum $\frac{1}{8}$ is attained at a corner, giving the stability constant.

Exercise 7 (hard, symbolic).

Prove the discrete maximum principle for $A_{h}$ : if a grid function $V$ satisfies $- Δ_{h} V_{ij} \leq 0$ at every interior point (discrete subharmonicity), then $max_{\overline{Ω}} V$ is attained on the boundary.

Hint

At an interior maximum, the centre value is at least the average of the four neighbours. Combine with the inequality $- Δ_{h} V \leq 0$ , which says the centre is at most the average.

Answer

The inequality $- Δ_{h} V_{ij} \leq 0$ rearranges to $4 V_{ij} \leq V_{i - 1, j} + V_{i + 1, j} + V_{i, j - 1} + V_{i, j + 1}$ , that is $V_{ij}$ is at most the average of its four neighbours (discrete sub-mean-value property). Suppose the maximum $M = max_{\overline{Ω}} V$ is attained at an interior point $(i_{0}, j_{0})$ . Then each neighbour value is $\leq M = V_{i_{0} j_{0}}$ , while the average inequality forces $V_{i_{0} j_{0}} \leq \frac{1}{4} \sum_{nbrs} V \leq M$ . Equality throughout requires every neighbour to equal $M$ . Thus the set of interior points where $V = M$ is closed under passing to neighbours; since the grid is connected, the maximum propagates to a boundary-adjacent interior point and then to a boundary point, so $M$ is attained on the boundary. (If the maximum is attained only in the interior, the propagation forces $V \equiv M$ , in which case it is also attained on the boundary.) Rubric: full credit for the sub-mean-value rearrangement, the equality-propagation argument, and the connectedness conclusion.

Exercise 8 (hard, symbolic).

Using the explicit eigenvalues $λ_{p, q} = \frac{4}{h ^{2}} (sin^{2} \frac{p π h}{2} + sin^{2} \frac{q π h}{2})$ of $A_{h}$ on the unit square, show that the spectral condition number satisfies $κ_{2} (A_{h}) = O (h^{- 2})$ as $h \to 0$ .

Hint

$A_{h}$ is symmetric positive definite, so $κ_{2} = λ_{m a x} / λ_{m i n}$ . Take $(p, q) = (1, 1)$ for the smallest and $(p, q) = (m, m)$ for the largest.

Answer

For symmetric positive definite $A_{h}$ , $κ_{2} (A_{h}) = λ_{m a x} / λ_{m i n}$ . The smallest eigenvalue is $λ_{1, 1} = \frac{8}{h ^{2}} sin^{2} (\frac{π h}{2})$ ; using $sin (θ) \approx θ$ for small $θ = π h /2$ , $λ_{1, 1} \approx \frac{8}{h ^{2}} \cdot \frac{π ^{2} h ^{2}}{4} = 2 π^{2} + O (h^{2})$ , a constant. The largest is $λ_{m, m} = \frac{8}{h ^{2}} sin^{2} (\frac{mπ h}{2})$ with $mh = 1 - h \to 1$ , so $\frac{mπ h}{2} \to \frac{π}{2}$ and $sin^{2} \to 1$ , giving $λ_{m, m} \approx \frac{8}{h ^{2}}$ . Hence $κ_{2} (A_{h}) \approx \frac{8/ h ^{2}}{2 π ^{2}} = \frac{4}{π ^{2} h ^{2}} = O (h^{- 2})$ . The condition number grows quadratically as the mesh refines, which is why naive elimination becomes expensive and iterative/multigrid solvers are eventually preferred. Rubric: full credit for the two eigenvalue limits and the $h^{- 2}$ scaling.

Advanced results Master

Theorem 1 ( $A_{h}$ is a symmetric positive definite M-matrix). The matrix $A_{h} = I_{m} \otimes T + T \otimes I_{m}$ , with $T = \frac{1}{h ^{2}} tridiag (- 1, 2, - 1)$ , is symmetric, has positive diagonal entries $4/ h^{2}$ and nonpositive off-diagonal entries $- 1/ h^{2}$ , is irreducible (the grid graph is connected) and weakly diagonally dominant with strict dominance in the boundary-adjacent rows, and is therefore a nonsingular M-matrix: $A_{h}^{- 1} \geq 0$ entrywise. It is moreover positive definite, since $T$ is positive definite and the Kronecker sum of positive-definite operators is positive definite, with eigenvalues the pairwise sums of the eigenvalues of $T$ . Monotonicity ( $A_{h} V \geq 0 \Rightarrow V \geq 0$ ) is equivalent to $A_{h}^{- 1} \geq 0$ and is the matrix form of the discrete maximum principle ^{[Varga, R. S. — Matrix Iterative Analysis (2nd ed.)]}.

Theorem 2 (exact eigensystem of the discrete Laplacian). On the unit square with $h = 1/ (m + 1)$ , the eigenvectors of $A_{h}$ are the tensor products of the 1-D sine modes: indexing $1 \leq p, q \leq m$ , $$ v^{(p,q)}{ij} = \sin(p\pi x_i)\sin(q\pi y_j), \qquad A_h v^{(p,q)} = \lambda{p,q},v^{(p,q)}, \qquad \lambda_{p,q} = \frac{4}{h^2}\Big(\sin^2\tfrac{p\pi h}{2} + \sin^2\tfrac{q\pi h}{2}\Big). $$ These are the exact discrete analogues of the continuous Dirichlet eigenfunctions $sin (p π x) sin (q π y)$ of $- Δ$ on the square, whose eigenvalues are $π^{2} (p^{2} + q^{2})$ ; the discrete eigenvalues approximate these for low modes ( $λ_{p, q} \to π^{2} (p^{2} + q^{2})$ as $h \to 0$ with $p, q$ fixed) and saturate near $8/ h^{2}$ for the highest modes. The smallest eigenvalue tends to $2 π^{2}$ and the largest to $8/ h^{2}$ , so $κ_{2} (A_{h}) = O (h^{- 2})$ ^{[Hackbusch, W. — Elliptic Differential Equations: Theory and Numerical Treatment]}.

Theorem 3 (the discrete Green's function and the sharp stability constant). The inverse $A_{h}^{- 1}$ has nonnegative entries $G_{h} (P, Q) \geq 0$ , the discrete Green's function, satisfying $\sum_{Q} G_{h} (P, Q) = ϕ (P)$ for the barrier $ϕ$ with $- Δ_{h} ϕ = 1$ ; consequently $∥ A_{h}^{- 1} ∥_{\infty} = max_{P} \sum_{Q} G_{h} (P, Q) = max_{P} ϕ (P) \leq \frac{1}{8}$ . The constant $\frac{1}{8}$ is the discrete reflection of the continuous bound $∥ (- Δ)^{- 1} ∥_{L^{\infty} \to L^{\infty}} = \frac{1}{8}$ on the unit square for the comparison quadratic; it is mesh-independent, which is precisely the stability of the scheme. The maximum-norm convergence rate $O (h^{2})$ is the product of this $O (1)$ stability constant with the $O (h^{2})$ truncation error ^{[LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems]}.

Theorem 4 (gain of accuracy by the nine-point stencil and Richardson extrapolation). Replacing the five-point stencil by the nine-point stencil $\frac{1}{6 h ^{2}}$ -weighted combination, with the source corrected by $f + \frac{h ^{2}}{12} Δ f$ , yields a scheme whose truncation error is $O (h^{4})$ , the Mehrstellen (compact fourth-order) discretisation; the same M-matrix and maximum-principle structure persists. Alternatively, because the five-point error admits an asymptotic expansion in even powers of $h$ , $U_{h} - u = h^{2} e_{2} + h^{4} e_{4} + \dots$ with $h$ -independent grid functions $e_{k}$ , Richardson extrapolation of $U_{h}$ and $U_{h /2}$ produces an $O (h^{4})$ approximation $\frac{4 U _{h /2} - U _{h}}{3}$ without changing the stencil. Both are routes past the second-order barrier of the plain five-point method ^{[Hackbusch, W. — Elliptic Differential Equations: Theory and Numerical Treatment]}.

Synthesis. The five-point Laplacian is the model case in which consistency and stability are each visible as a separate, checkable structure, and convergence is their product: the foundational reason the method works is that $A_{h}$ is an M-matrix, and this single algebraic fact carries the discrete maximum principle, the inverse-monotonicity $A_{h}^{- 1} \geq 0$ , the sharp mesh-independent stability constant $\frac{1}{8}$ , and the nonnegativity of the discrete Green's function. This is exactly the discrete shadow of the continuous Poisson theory of 02.13.02 and the continuous maximum principle of 02.13.01: the comparison/barrier function used in the stability proof is the discrete restriction of the continuous barrier, and the discrete eigenmodes are the exact sampled continuous Dirichlet eigenfunctions. The central insight is that the maximum norm, not the energy norm, is the natural setting for the finite-difference method — the M-matrix monotonicity is an $L^{\infty}$ phenomenon — which is the structural point of divergence from the finite-element method of 24.02.01, whose convergence is read in the energy norm through Céa's lemma and Galerkin orthogonality rather than through a maximum principle. Putting these together, the consistency-plus-stability proof here generalises into the Lax-Richtmyer equivalence theorem 43.11.05 for evolution problems, where $A_{h}^{- 1}$ is replaced by the power-bounded time-stepping operator; and the very matrix $A_{h}$ , with its explicit spectrum, is dual to the temporal stiffness it induces in 43.11.02, where its largest eigenvalue $\sim 8/ h^{2}$ becomes the source of the parabolic time-step restriction.

Full proof set Master

Proposition 1 (discrete maximum principle). If a grid function $V$ on $\overline{Ω}_{h}$ satisfies $(- Δ_{h} V)_{ij} \leq 0$ at every interior point, then $max_{\overline{Ω}_{h}} V = max_{\partial Ω_{h}} V$ . Symmetrically, $(- Δ_{h} V)_{ij} \geq 0$ implies $min_{\overline{Ω}_{h}} V = min_{\partial Ω_{h}} V$ .

Proof. The hypothesis $(- Δ_{h} V)_{ij} \leq 0$ is $4 V_{ij} \leq V_{i - 1, j} + V_{i + 1, j} + V_{i, j - 1} + V_{i, j + 1}$ , so $V_{ij}$ is bounded above by the mean of its four neighbours. Let $M = max_{\overline{Ω}_{h}} V$ and suppose $M$ is attained at an interior point $P_{0}$ . Each neighbour satisfies $V \leq M$ , while the sub-mean inequality forces $V (P_{0}) \leq \frac{1}{4} \sum_{nbrs} V \leq M = V (P_{0})$ , so equality holds and every neighbour also attains $M$ . The set $S = {P : V (P) = M}$ is thus closed under the neighbour relation among interior points. Since the interior grid graph is connected and every interior point is connected to a boundary point through grid edges, propagating equality along a path from $P_{0}$ reaches an interior point adjacent to the boundary; at that point one of its neighbours is a boundary point, which must then also equal $M$ . Hence $M$ is attained on the boundary. The minimum statement follows by applying the result to $- V$ . $□$

Proposition 2 ( $A_{h}$ is monotone: $A_{h}^{- 1} \geq 0$ ). The matrix $A_{h}$ representing $- Δ_{h}$ on the interior with Dirichlet boundary is nonsingular, and $A_{h}^{- 1} \geq 0$ entrywise; equivalently, $A_{h} V \geq 0$ implies $V \geq 0$ .

Proof. First, nonsingularity. If $A_{h} V = 0$ with $V = 0$ on the boundary, then $V$ satisfies both $- Δ_{h} V = 0 \leq 0$ and $\geq 0$ at every interior point, so by Proposition 1 its maximum and minimum are both attained on the boundary where $V = 0$ ; hence $V \equiv 0$ and $A_{h}$ is injective, thus invertible. For monotonicity, suppose $A_{h} V \geq 0$ , i.e. $- Δ_{h} V_{ij} \geq 0$ , with $V = 0$ on the boundary. By the minimum form of Proposition 1, $min_{\overline{Ω}_{h}} V = min_{\partial Ω_{h}} V = 0$ , so $V \geq 0$ everywhere. Applying this to $V = A_{h}^{- 1} e_{k}$ for each standard basis vector $e_{k} \geq 0$ shows every column of $A_{h}^{- 1}$ is nonnegative, hence $A_{h}^{- 1} \geq 0$ . The sign pattern (positive diagonal, nonpositive off-diagonals) together with this nonnegative inverse is the definition of an M-matrix. $□$

Proposition 3 (uniform stability bound $∥ A_{h}^{- 1} ∥_{\infty} \leq \frac{1}{8}$ ). For every $h = 1/ (m + 1)$ , the discrete solution operator satisfies $∥ A_{h}^{- 1} ∥_{\infty} \leq \frac{1}{8}$ , a bound independent of $h$ .

Proof. Let $ϕ$ be the grid restriction of $Φ (x, y) = \frac{1}{4} ((x - \frac{1}{2})^{2} + (y - \frac{1}{2})^{2})$ , for which $- ΔΦ = - 1$ and, since $Φ$ is quadratic and the centred second difference is exact on quadratics, $- Δ_{h} ϕ_{ij} = - 1$ at every interior point; equivalently $A_{h} (- ϕ) = - 1$ where $1$ is the all-ones interior vector, and on the unit square $0 \leq ϕ \leq \frac{1}{8}$ . Fix a right-hand side $R$ and let $W = A_{h}^{- 1} R$ , so $W = 0$ on the boundary. Consider $ψ = W - ∥ R ∥_{\infty} (\frac{1}{8} 1 - ϕ)$ understood at interior points with the convention that $\frac{1}{8} - ϕ \geq 0$ ; more directly, set $Z^{\pm} = \pm W - ∥ R ∥_{\infty} ϕ$ with $ϕ$ the barrier satisfying $- Δ_{h} ϕ = 1$ , $0 \leq ϕ \leq \frac{1}{8}$ , $ϕ \geq 0$ on the boundary. Then $- Δ_{h} Z^{\pm} = \pm R_{ij} - ∥ R ∥_{\infty} \leq 0$ at interior points and $Z^{\pm} \leq 0$ on the boundary (there $\pm W = 0$ and $- ∥ R ∥_{\infty} ϕ \leq 0$ ). By Proposition 1, $Z^{\pm} \leq 0$ throughout, so $\pm W_{ij} \leq ∥ R ∥_{\infty} ϕ_{ij} \leq \frac{1}{8} ∥ R ∥_{\infty}$ for all $i, j$ . Hence $∥ W ∥_{\infty} \leq \frac{1}{8} ∥ R ∥_{\infty}$ , i.e. $∥ A_{h}^{- 1} R ∥_{\infty} \leq \frac{1}{8} ∥ R ∥_{\infty}$ for every $R$ , which is $∥ A_{h}^{- 1} ∥_{\infty} \leq \frac{1}{8}$ . $□$

Proposition 4 (exact eigenvalues and eigenvectors of $A_{h}$ ). The vectors $v_{ij}^{(p, q)} = sin (p π x_{i}) sin (q π y_{j})$ , $1 \leq p, q \leq m$ , are a complete orthogonal eigenbasis of $A_{h}$ with eigenvalues $λ_{p, q} = \frac{4}{h ^{2}} (sin^{2} \frac{p π h}{2} + sin^{2} \frac{q π h}{2})$ .

Proof. It suffices to treat the 1-D operator $T = \frac{1}{h ^{2}} tridiag (- 1, 2, - 1)$ , since $A_{h} = I \otimes T + T \otimes I$ acts on the tensor mode $v^{(p, q)} = s^{(p)} \otimes s^{(q)}$ with $s_{i}^{(p)} = sin (p π x_{i})$ by $A_{h} v^{(p, q)} = (μ_{p} + μ_{q}) v^{(p, q)}$ where $T s^{(p)} = μ_{p} s^{(p)}$ . Apply $T$ to $s^{(p)}$ at interior index $i$ : $$ (T s^{(p)})i = \frac{-\sin(p\pi x{i-1}) + 2\sin(p\pi x_i) - \sin(p\pi x_{i+1})}{h^2}. $$ Using $x_{i \pm 1} = x_{i} \pm h$ and the identity $sin (θ - p π h) + sin (θ + p π h) = 2 sin θ cos (p π h)$ with $θ = p π x_{i}$ , $$ (T s^{(p)})_i = \frac{2\sin(p\pi x_i)\big(1 - \cos(p\pi h)\big)}{h^2} = \frac{2(1-\cos p\pi h)}{h^2},s^{(p)}i. $$ The boundary entries vanish because $s_{0}^{(p)} = sin 0 = 0$ and $s_{m + 1}^{(p)} = sin (p π) = 0$ , so $s^{(p)}$ satisfies the homogeneous Dirichlet condition and is a genuine interior eigenvector. Hence $μ_{p} = \frac{2 ( 1 - c o s p π h )}{h ^{2}} = \frac{4}{h ^{2}} sin^{2} \frac{p π h}{2}$ by the half-angle identity. The $m$ vectors $s^{(1)}, \dots, s^{(m)}$ are orthogonal (discrete sine transform) and complete in $R^{m}$ , so the tensor products $v^{(p, q)}$ form a complete orthogonal eigenbasis of $A_{h}$ with $\lambda{p,q} = \mu_p + \mu_q $a ss t a t e d .$ \square$

Connections Master

The continuous Poisson problem and its fundamental solution 02.13.02 is the object this unit discretises: the right-hand side $f$ , the Dirichlet data $g$ , and the classical solution $u \in C^{4}$ whose fourth derivatives bound the truncation error all come from there, and the continuous solvability and $H^{2}$ / $C^{2, α}$ regularity it supplies (sharpened in 02.17.01) are exactly the hypotheses the $O (h^{2})$ convergence theorem assumes. This unit reproves nothing about the continuous equation; it builds the discrete operator that approximates it and proves the approximation converges.
The continuous maximum principle for harmonic and subharmonic functions 02.13.01 is the model for the discrete maximum principle proved here: the sub-mean-value inequality $4 V_{ij} \leq \sum_{nbrs} V$ is the grid analogue of the sub-mean-value property of continuous subharmonic functions, and the comparison-function (barrier) technique used for the stability bound is the discrete restriction of the continuous comparison argument. The M-matrix structure of $A_{h}$ is the algebraic mechanism that makes the maximum principle survive discretisation.
The Gaussian-elimination and LU theory 43.03.01 is what actually solves the system $A_{h} U = F$ this unit produces: $A_{h}$ is symmetric positive definite, so it admits a pivot-free Cholesky factorization, and its banded sparsity (bandwidth $m$ under lexicographic ordering) makes a direct solve cost $O (m^{4}) = O (N^{2})$ rather than the dense $O (N^{3})$ . The conditioning $κ_{2} (A_{h}) = O (h^{- 2})$ established here feeds the forward-error analysis of that unit, explaining why very fine grids eventually motivate iterative and multigrid solvers over direct elimination.
The finite-element / Galerkin discretisation of the same Poisson problem 24.02.01 is the variational alternative to this pointwise-stencil method: where the finite-difference method enforces the equation at grid nodes and proves convergence through a discrete maximum principle in the maximum norm, the finite-element method projects the weak form onto a subspace and proves convergence through Céa's lemma and Galerkin orthogonality in the energy norm. On a uniform mesh the lowest-order finite-element stiffness matrix for $- Δ$ coincides (up to scaling) with the five-point matrix $A_{h}$ , so the two discretisations meet at this model problem while resting on different stability theories — the energy estimate there, the M-matrix monotonicity here.
The eigenvalue and eigenvector theory 01.01.08 supplies the spectral machinery this unit applies to $A_{h}$ : the explicit diagonalisation by sine modes, the identification of $κ_{2} (A_{h}) = λ_{m a x} / λ_{m i n}$ for the symmetric positive definite operator, and the spectral-radius reasoning that the next unit 43.11.02 uses to convert the largest eigenvalue $\sim 8/ h^{2}$ of $A_{h}$ into the explicit-scheme time-step restriction for the heat equation.

Historical & philosophical context Master

The systematic use of finite differences to approximate the Laplace and Poisson equations was placed on a rigorous footing by Richard Courant, Kurt Friedrichs, and Hans Lewy in their 1928 paper Über die partiellen Differenzengleichungen der mathematischen Physik (Mathematische Annalen 100, 32–74) ^{[Courant, Friedrichs & Lewy 1928]}, which proved convergence of the five-point discrete solution to the solution of the continuous boundary-value problem and introduced the stability condition that bears their initials in the time-dependent case. Their elliptic argument already used the discrete maximum principle and the comparison-function idea that organise the modern proof.

The recasting of the convergence question into the algebraic language of monotone matrices and M-matrices — the property $A^{- 1} \geq 0$ that carries the discrete maximum principle — is due to the work surrounding Richard Varga's Matrix Iterative Analysis (Prentice-Hall, 1962; 2nd ed. Springer, 2000) ^{[Varga, R. S. — Matrix Iterative Analysis (2nd ed.)]} and the matrix-theoretic tradition initiated by Ostrowski and Collatz, who connected diagonal dominance and irreducibility to inverse-positivity. The comprehensive numerical-analytic treatment, including the sharp $∥ A_{h}^{- 1} ∥_{\infty}$ bound, the discrete Green's function, and the exact eigensystem of the discrete Laplacian, was consolidated by Wolfgang Hackbusch in Elliptic Differential Equations (Springer, 1992) ^{[Hackbusch, W. — Elliptic Differential Equations: Theory and Numerical Treatment]}, the same period in which Hackbusch's multigrid methods made the resulting sparse systems efficiently solvable.

Bibliography Master

@book{leveque2007fdm,
  author    = {LeVeque, Randall J.},
  title     = {Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems},
  publisher = {Society for Industrial and Applied Mathematics (SIAM)},
  year      = {2007}
}

@book{hackbusch1992elliptic,
  author    = {Hackbusch, Wolfgang},
  title     = {Elliptic Differential Equations: Theory and Numerical Treatment},
  series    = {Springer Series in Computational Mathematics},
  volume    = {18},
  publisher = {Springer-Verlag},
  year      = {1992}
}

@book{varga2000matrix,
  author    = {Varga, Richard S.},
  title     = {Matrix Iterative Analysis},
  edition   = {2},
  publisher = {Springer-Verlag},
  year      = {2000}
}

@book{iserles2009first,
  author    = {Iserles, Arieh},
  title     = {A First Course in the Numerical Analysis of Differential Equations},
  edition   = {2},
  publisher = {Cambridge University Press},
  year      = {2009}
}

@article{courantfriedrichslewy1928,
  author  = {Courant, Richard and Friedrichs, Kurt and Lewy, Hans},
  title   = {{\"U}ber die partiellen Differenzengleichungen der mathematischen Physik},
  journal = {Mathematische Annalen},
  volume  = {100},
  number  = {1},
  year    = {1928},
  pages   = {32--74}
}

@book{strang2007cse,
  author    = {Strang, Gilbert},
  title     = {Computational Science and Engineering},
  publisher = {Wellesley-Cambridge Press},
  year      = {2007}
}

Prerequisites

02.13.02
02.13.01
43.03.01
01.01.08

Tier anchors

beginner: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §2.1-2.4, §3.1-3.2 (the 1-D second-difference operator, the 5-point Laplacian on a square, and the picture of the sparse linear system); Strang 2007 *Computational Science and Engineering* (Wellesley-Cambridge) §1.1, §3.5 (the second-difference matrix and the 5-point stencil as a discrete Laplacian)
intermediate: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) §3.2-3.5 (the 5-point stencil, the local truncation error $O(h^2)$, the discrete maximum principle, stability through $\|A_h^{-1}\|_\infty$, and the convergence theorem); Iserles 2009 *A First Course in the Numerical Analysis of Differential Equations* 2e (Cambridge) §8.1-8.2 (the five-point formula, monotonicity, and the maximum-principle convergence proof)
master: LeVeque 2007 *Finite Difference Methods for Ordinary and Partial Differential Equations* (SIAM) Ch. 3 (the elliptic discretisation, the M-matrix structure, the discrete Green's function, and the $O(h^2)$ convergence theorem); Hackbusch 1992 *Elliptic Differential Equations: Theory and Numerical Treatment* (Springer) §4.1-4.7 (the discrete maximum principle, M-matrices, the inverse-monotone stability bound, and the explicit eigensystem of the discrete Laplacian); Varga 2000 *Matrix Iterative Analysis* 2e (Springer) Ch. 1, 3 (M-matrices, monotonicity, and the diagonally dominant irreducible structure of $A_h$)

References

LeVeque, R. J. — Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems · SIAM 2007. Ch. 2-3. §2.4 derives the second-difference operator $D^2 U_j=(U_{j-1}-2U_j+U_{j+1})/h^2$ and its $-\tfrac{1}{12}h^2 u''''$ truncation error by Taylor expansion. §3.2 builds the 5-point Laplacian for $-\Delta u=f$ on a square with Dirichlet data, writes the sparse linear system $A_h U=F$, identifies $A_h$ as symmetric negative definite (here written for $-\Delta_h$ as symmetric positive definite), and gives the local truncation error $\tau=O(h^2)$. §3.3 states the discrete maximum principle and the stability bound $\|A_h^{-1}\|_\infty\le 1/8$ via a comparison function; §3.5 assembles consistency plus stability into the global error theorem $\|U-u\|_\infty=O(h^2)$. The explicit eigenvalues $\lambda_{p,q}$ and eigenvectors of the discrete Laplacian on the square are given in §3.4.
Hackbusch, W. — Elliptic Differential Equations: Theory and Numerical Treatment · Springer 1992. §4.1-4.7. The 5-point discretisation of the Poisson problem; the discrete maximum principle for the matrix $A_h$; the proof that $A_h$ is an M-matrix (monotone, $A_h^{-1}\ge0$ entrywise) from irreducible diagonal dominance and the sign pattern (positive diagonal, nonpositive off-diagonals); the inverse-monotone stability estimate via a discrete comparison/barrier function; and the explicit eigenvalue-eigenvector system of the discrete Laplacian on a rectangle, with the smallest eigenvalue $\sim 2\pi^2$ controlling the conditioning $\kappa(A_h)=O(h^{-2})$.
Varga, R. S. — Matrix Iterative Analysis (2nd ed.) · Springer 2000. Ch. 1, 3. M-matrices and monotone matrices ($A^{-1}\ge0$); the theorem that an irreducibly diagonally dominant matrix with positive diagonal and nonpositive off-diagonal entries is a nonsingular M-matrix; the discrete Laplacian $A_h$ as the model M-matrix; the connection to the discrete maximum principle and to the convergence of classical iterative solvers for the resulting sparse system.

Estimated time

beginner: 20m
intermediate: 50m
master: 90m