Boundary Value Problems

Learning Objectives

By the end of this section you will be able to:

1. State the three classical types of boundary conditions for Laplace's equation (Dirichlet, Neumann, Robin) and recognise which one a physical setup demands.
2. Explain why a PDE alone — even one as rigid as $\nabla^2 u = 0$ — has infinitely many solutions, and how boundary data picks out exactly one.
3. Use the maximum principle to argue uniqueness and to bound the interior of a harmonic function in terms of its boundary.
4. Check the compatibility condition $\oint_{\partial\Omega} g\,ds = 0$ that a pure Neumann problem must satisfy.
5. Implement a Jacobi relaxation solver from scratch in Python, and a differentiable version in PyTorch that you can use inside a learning loop.

The Big Picture

In section 28.1 we introduced Laplace's equation $\nabla^2 u = u_{xx} + u_{yy} = 0$ and saw that its solutions — harmonic functions — model equilibrium: steady-state heat, electrostatic potential, ideal-fluid streamfunctions, gravitational potential outside masses. Equilibrium is what is left after every transient has died out.

A natural question follows immediately: given that the interior is in equilibrium, what is it in equilibrium with? The answer is the boundary. A drumhead at rest is in equilibrium with the rim it is glued to. A metal plate's steady temperature is in equilibrium with whatever the edges are doing. A capacitor's electric potential is in equilibrium with the voltages on its plates.

Slogan. Laplace's equation is the law of the interior. Boundary conditions are the law of the edge. A Boundary Value Problem (BVP) is the two of them welded together — and only then does the problem have a unique answer.

Why a PDE Needs Boundary Conditions

For an ordinary differential equation, you fix the solution with initial conditions: tell me where the particle is and how fast it is moving at $t = 0$, and I will tell you $x(t)$ for all later time.

For Laplace's equation there is no time. Time has already finished — we are sitting in the equilibrium that remains. So we cannot pin the answer down with initial conditions. Instead we pin it down with spatial conditions on the edge of the region.

Here is the most striking demonstration. The function $u(x, y) = 0$ is harmonic. So is $u(x, y) = x$. So is $u(x, y) = x^2 - y^2$. So is $u(x, y) = e^{x} \cos(y)$. All four solve the same PDE $u_{xx} + u_{yy} = 0$ identically. The PDE alone cannot distinguish them — they look different only at the boundary.

The Three Classical Boundary Conditions

On the boundary $\partial\Omega$ of our region $\Omega \subset \mathbb{R}^2$ we can prescribe the value of $u$, the outward normal derivative $\partial u / \partial n$, or a linear combination of the two. These give the three classical boundary conditions.

Dirichlet: Prescribed Value

A Dirichlet boundary condition prescribes the value of the solution on the boundary:

Physical picture. A metal plate whose edges are welded to ice baths and steam pipes: the temperature on the edge is whatever the bath dictates, no matter what is happening inside. Electrostatically, a Dirichlet boundary is a conductor held at a fixed voltage by a battery.

Neumann: Prescribed Flux

A Neumann boundary condition prescribes the outward normal derivative:

where $\mathbf{n}$ is the outward unit normal to $\partial\Omega$.

Physical picture. Now we are not specifying the temperature on the edge — we are specifying the rate at which heat flows through it. If $g = 0$, the boundary is perfectly insulated: no heat enters or leaves. If $g > 0$ at some point, heat is being pumped out there. Electrostatically, prescribing $\partial u / \partial n$ on a surface is the same as prescribing the surface charge density (up to a constant).

Robin: A Mix of Both

A Robin (or "third-kind") boundary condition is a linear combination:

Physical picture. Newton's law of cooling. Suppose the boundary is a thin wall exchanging heat with surroundings at temperature $T_\infty$. Heat flux through the wall is proportional to the temperature difference:

Rearranged, this is exactly the Robin form. As $h_{\text{conv}} \to \infty$ Robin degenerates to Dirichlet (the wall instantly matches the outside temperature). As $h_{\text{conv}} \to 0$ Robin degenerates to Neumann with $g = 0$ (the wall is a perfect insulator).

Mixed vs. Robin — Don't Confuse Them

Two phrases sound similar but mean different things:

What "Well-Posed" Really Means

Jacques Hadamard introduced three criteria a sensible BVP should meet. We call the problem well-posed if all three hold:

1. Existence. A solution exists for any reasonable boundary data.
2. Uniqueness. No two distinct solutions share the same boundary data.
3. Continuous dependence. Small changes in the boundary data produce small changes in the solution. Tiny measurement errors cannot blow up into wildly different temperatures.

The Dirichlet problem for Laplace's equation on a bounded well-behaved domain satisfies all three. The Neumann problem satisfies all three once the compatibility condition holds and we agree to consider solutions modulo an additive constant. Robin is well-posed whenever $\alpha, \beta > 0$.

The Maximum Principle (the Secret Weapon)

Harmonic functions have a remarkable shape constraint. The strong maximum principle says:

Theorem. Let $\Omega$ be a bounded connected open set in $\mathbb{R}^n$ and suppose $u$ is harmonic in $\Omega$ and continuous on $\overline{\Omega}$. Then $u$ attains its maximum and minimum on the boundary $\partial\Omega$. Moreover, if $u$ attains its max or min at any interior point, then $u$ is constant.

Why does this happen? Because of the mean-value property: $u(\mathbf{x}_0)$ is the average of $u$ on any small disk around $\mathbf{x}_0$. If $\mathbf{x}_0$ were strictly larger than every value in its neighbourhood, it could not equal their average — a contradiction. So a strict interior maximum is impossible.

Uniqueness from the Maximum Principle

The maximum principle gives a one-line proof of uniqueness for the Dirichlet BVP. Suppose $u_1$ and $u_2$ both solve $\nabla^2 u = 0$ in $\Omega$ with the same boundary data $u = f$ on $\partial\Omega$. Let $w = u_1 - u_2$. Then $w$ is harmonic in $\Omega$ with $w = 0$ on $\partial\Omega$. By the maximum principle, $\max w = \min w = 0$, so $w \equiv 0$ in $\Omega$ — meaning $u_1 \equiv u_2$. The Dirichlet problem has at most one solution.

Interactive: Maximum Principle Demo

Below is a numerical solver for Laplace's equation on the unit square with four piecewise-constant Dirichlet boundaries. Drag the probe around the interior and watch the inequality $\min_{\partial\Omega} u \leq u(\mathbf{x}) \leq \max_{\partial\Omega} u$ hold at every point. You can never find an interior point hotter than the hottest edge or colder than the coldest edge.

The Neumann Compatibility Condition

For the pure Neumann problem

the data $g$ cannot be arbitrary. Take the PDE and integrate over $\Omega$:

The leftmost integral is zero because $\nabla^2 u = 0$ in $\Omega$. So:

Physical interpretation. For a steady temperature field, the net heat flux through the entire boundary must vanish — otherwise heat would be accumulating or being lost, and the interior could not be steady. If you try to pump heat into a fully sealed room without letting any out, there is no steady-state temperature; the room just keeps heating up.

Interactive: Compatibility Check

Slide the four normal-derivative values. The widget shows the net flux $\oint g\,ds$ and turns green only when the data is compatible.

Interactive: Solve a BVP Yourself

The widget below is a full Jacobi relaxation solver on the unit square. You can switch each of the four edges between Dirichlet and Neumann, adjust the boundary values, and press Run to watch the field relax to equilibrium. The red border marks Dirichlet edges; the blue border marks Neumann edges. The iteration counter tells you how many Jacobi sweeps have been applied so far.

Things to try. Pick the first preset (hot top, cold bottom, cold sides) and watch a thermal "dome" develop near the top edge. Switch to the third preset (insulated sides) and notice that the solution becomes purely vertical: $u = y$. Try the all-Neumann case (set every side to Neumann) — if the four fluxes do not sum to zero, the solver still iterates, but the compatibility warning lights up and the field drifts without settling.

Worked Example (Step by Step)

Let's trace the Jacobi solver by hand on a tiny grid you can fit on a napkin. Take the unit square with $N = 5$ grid points per axis (so $3 \times 3 = 9$ interior unknowns). Use Dirichlet boundaries: top edge $u = 1$, the other three edges $u = 0$. We start with $u = 0$ everywhere in the interior.

row i=4   [1   1   1   1   1]   <- top boundary (hot)
row i=3   [0   0   0   0   0]
row i=2   [0   0   0   0   0]
row i=1   [0   0   0   0   0]
row i=0   [0   0   0   0   0]   <- bottom (cold)
            j=0 j=1 j=2 j=3 j=4

row i=3   [0  0.25 0.25 0.25  0]
row i=2   [0  0    0    0     0]
row i=1   [0  0    0    0     0]

row i=3   [0  0.3125 0.375  0.3125 0]
row i=2   [0  0.0625 0.0625 0.0625 0]
row i=1   [0  0      0      0      0]

row i=3   [0  0.359375 0.421875 0.359375 0]
row i=2   [0  0.09375  0.125    0.09375  0]
row i=1   [0  0.015625 0.015625 0.015625 0]

u(centre) = 0.25   (exact, by symmetry + linearity)
Jacobi after 2048 iterations on a 33x33 grid: 0.24998

Python Implementation (Jacobi Relaxation)

Below is the same algorithm we just traced by hand, written in plain Python plus NumPy. Hover the cards on the left to highlight the corresponding code line.

south neighbours

north neighbours

west  neighbours

east  neighbours

1import numpy as np
2
3def laplace_bvp_jacobi(bc, N=33, tol=1e-6, max_iter=20000):
4    """
5    Solve Laplace's equation  div(grad u) = 0  on the unit square (0,1)x(0,1)
6    by Jacobi relaxation on an N x N grid.
7
8    Parameters
9    ----------
10    bc : dict with keys "top", "bottom", "left", "right".
11         Each value is a tuple ("dirichlet", value)  or  ("neumann", flux).
12    N        : number of grid points along each axis (including boundary).
13    tol      : convergence tolerance on max |u_new - u_old|.
14    max_iter : safety cap on iterations.
15
16    Returns
17    -------
18    u    : the converged solution as an (N, N) array, u[i, j] = u(x_j, y_i).
19    iters: number of iterations actually used.
20    """
21    h = 1.0 / (N - 1)
22    u = np.zeros((N, N))
23
24    def apply_bc(u):
25        # top row  (y = 1)
26        kind, v = bc["top"]
27        if kind == "dirichlet":
28            u[-1, :] = v
29        else:                                  # outward normal = +y
30            u[-1, :] = u[-2, :] + h * v        # ghost-point form
31
32        # bottom row (y = 0)
33        kind, v = bc["bottom"]
34        if kind == "dirichlet":
35            u[0, :] = v
36        else:                                  # outward normal = -y
37            u[0, :] = u[1, :] + h * v
38
39        # left column (x = 0)
40        kind, v = bc["left"]
41        if kind == "dirichlet":
42            u[:, 0] = v
43        else:                                  # outward normal = -x
44            u[:, 0] = u[:, 1] + h * v
45
46        # right column (x = 1)
47        kind, v = bc["right"]
48        if kind == "dirichlet":
49            u[:, -1] = v
50        else:                                  # outward normal = +x
51            u[:, -1] = u[:, -2] + h * v
52
53    apply_bc(u)
54
55    for it in range(1, max_iter + 1):
56        u_old = u.copy()
57        # Jacobi sweep:  u_{i,j} <- 1/4 (north + south + east + west)
58        u[1:-1, 1:-1] = 0.25 * (
59            u_old[:-2, 1:-1] +     # south neighbor
60            u_old[2:,  1:-1] +     # north neighbor
61            u_old[1:-1, :-2] +     # west  neighbor
62            u_old[1:-1, 2:]        # east  neighbor
63        )
64        apply_bc(u)
65        err = np.max(np.abs(u - u_old))
66        if err < tol:
67            return u, it
68    return u, max_iter
69
70
71# --- Demo: hot top, cold bottom, insulated sides ---
72bc = {
73    "top":    ("dirichlet", 1.0),
74    "bottom": ("dirichlet", 0.0),
75    "left":   ("neumann",   0.0),
76    "right":  ("neumann",   0.0),
77}
78u, n_iter = laplace_bvp_jacobi(bc, N=33, tol=1e-7)
79print(f"Converged in {n_iter} iterations.")
80print(f"u(center)       = {u[16, 16]:.6f}")  # expected ~ 0.5 by symmetry
81print(f"max |u| on boundary = {max(u[0].max(), u[-1].max()):.3f}")

PyTorch Implementation (Differentiable)

Rewriting the solver in PyTorch unlocks something the NumPy version cannot do: backprop through the BVP. If we treat the boundary values as learnable parameters, the gradient of a downstream loss (e.g. desired interior temperature) flows all the way back to the boundary. This is the foundation of physics-informed learning and shape/topology optimisation.

1import torch
2
3def laplace_bvp_torch(top, bottom, left, right, N=33, n_iter=2000):
4    """
5    Differentiable Jacobi solver for Laplace's equation on the unit square.
6    All four boundary values are torch tensors that may carry gradients —
7    so we can backprop a downstream loss (e.g. desired interior pattern)
8    back through the solver to the boundary data.
9
10    All four sides are Dirichlet here for simplicity.
11    """
12    device = top.device
13    dtype  = top.dtype
14    u = torch.zeros(N, N, dtype=dtype, device=device)
15
16    # Paint the boundary
17    u[-1, :] = top
18    u[ 0, :] = bottom
19    u[:,  0] = left
20    u[:, -1] = right
21
22    for _ in range(n_iter):
23        # discrete Laplacian average — runs on GPU if u is on GPU
24        u_new = u.clone()
25        u_new[1:-1, 1:-1] = 0.25 * (
26            u[:-2, 1:-1] + u[2:, 1:-1] +
27            u[1:-1, :-2] + u[1:-1, 2:]
28        )
29        # re-paint boundary so it stays glued to the input tensors
30        u_new[-1, :] = top
31        u_new[ 0, :] = bottom
32        u_new[:,  0] = left
33        u_new[:, -1] = right
34        u = u_new
35
36    return u
37
38
39# --- We want u(centre) = target.  Optimise the top edge to achieve it. ---
40target_centre = 0.7
41
42top    = torch.full((33,), 1.0, requires_grad=True)
43bottom = torch.zeros(33)
44left   = torch.zeros(33)
45right  = torch.zeros(33)
46
47opt = torch.optim.Adam([top], lr=0.05)
48for step in range(200):
49    u = laplace_bvp_torch(top, bottom, left, right, N=33, n_iter=400)
50    loss = (u[16, 16] - target_centre) ** 2
51    opt.zero_grad()
52    loss.backward()
53    opt.step()
54    if step % 40 == 0:
55        print(f"step {step:3d}  loss={loss.item():.5f}  u_centre={u[16,16].item():.4f}")
56print(f"final mean(top) = {top.detach().mean().item():.4f}")

Real-World Applications

Image inpainting is a fun one: replace a missing region with the solution of Laplace's equation whose Dirichlet data is the surrounding known pixels. The mean-value property smoothly interpolates the hole — colours blend, gradients match, no jarring edges. Modern deep-learning inpainters are more powerful (they can hallucinate structure), but the Laplace-equation inpainter is still the right thing to use for smooth backgrounds.

Common Pitfalls

• Forgetting the compatibility condition. If you set up a pure Neumann problem with $\oint g\,ds \neq 0$, no solution exists. The solver will diverge or drift; the discretised linear system is singular.
• Forgetting the constant freedom in pure Neumann. Even compatible Neumann data leaves $u$ defined only up to $+C$. Add one additional pin like "average of u over $\Omega$ is zero" to make it unique.
• Mixing up ∂u/∂n and ∂u/∂x. The Neumann condition uses the outward normal derivative. On the right edge of a rectangle, this is $+\partial u/\partial x$; on the left edge it is $-\partial u/\partial x$.
• First-order Neumann in code. The ghost-point trick $u_{\text{bd}} = u_{\text{nbr}} + h g$ is only first-order accurate. For high precision, use second-order one-sided differences $u_{\text{bd}} = \tfrac{4}{3} u_{\text{nbr,1}} - \tfrac{1}{3} u_{\text{nbr,2}} + \tfrac{2}{3} h g$ (Taylor expand a quadratic to confirm).
• Jacobi is slow. Convergence is $O(N^2)$ iterations. For large grids use multigrid, conjugate gradient, or a direct sparse solver.
• Re-entrant corners. On non-convex domains (think of an L-shape) the solution can have a singularity at the inward corner: $u \sim r^{2/3}$. The PDE still well-posed in the variational sense, but pointwise convergence near the corner is slow and grid-refinement is essential.

Test Your Understanding

Summary

Key Takeaways

1. A PDE on its own is under-determined; boundary conditions turn it into a well-posed problem.
2. The three classical types — Dirichlet, Neumann, Robin — cover essentially every physical scenario you will meet in heat conduction, electrostatics, fluid flow, and elastic membranes.
3. The maximum principle is the secret weapon: it gives uniqueness, monotonicity, and pointwise bounds, all in one stroke.
4. Pure Neumann problems need compatibility: $\oint_{\partial\Omega} g\,ds = 0$. Net flux must vanish, otherwise no steady state exists.
5. Jacobi relaxation is the discrete embodiment of the mean-value property — and writing it in PyTorch makes the entire BVP differentiable, which is the foundation of modern physics-informed deep learning.

Coming next. In section 28.3 we leave abstract boundary types behind and solve a concrete problem: Laplace's equation on a rectangle, by separation of variables and Fourier series. The maximum principle and compatibility conditions you have built here will keep guiding us — but the algebra will get specific.