Applications: Steady-State Heat

Learning Objectives

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

1. Derive Laplace's equation $\nabla^2 u = 0$ as the long-time limit of the heat equation.
2. Explain the relaxation rule $u_{i,j} \leftarrow \tfrac{1}{4}\bigl(u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1}\bigr)$ and why averaging neighbours encodes the Laplacian.
3. State and use the mean-value property and the maximum principle to predict, without computation, whether an interior point can ever be hotter than the boundary.
4. Solve a steady-state heat problem by hand on a small grid and on a large grid with Python / PyTorch.
5. Recognise where steady-state heat models appear in engineering — chip cooling, building insulation, planetary temperatures, image inpainting.

From Time Evolution to Equilibrium

"Wait long enough, and every temperature settles down. What remains is a frozen portrait of how heat balances itself."

Heat does not stand still. Pour hot coffee into a cold mug and the temperature changes everywhere in the system — the mug warms, the coffee cools, the air carries energy away. The transient stage of this story is told by the heat equation, $\partial_t u = \alpha\,\nabla^2 u$. But in many practical situations we do not care about the transient at all. We care about what the temperature looks like once nothing is changing any more: the engine block has been running for hours, the heat sink has reached its operating temperature, the planet has been receiving sunlight for billions of years.

When nothing changes any more, $\partial_t u = 0$. Plug that into the heat equation and the time derivative on the left collapses to zero, leaving

This is Laplace's equation. Every problem in this chapter — electrostatic potential, incompressible flow, gravitational potential in empty space — reduces to the same PDE because they all describe equilibrium states. Steady-state heat is simply the most physically intuitive entry point.

Warm-Up: The 1D Rod

Before going to two dimensions, settle the intuition in one. Take a thin uniform rod of length $L$, perfectly insulated along its sides, held at temperature $T_L$ on the left end and $T_R$ on the right. Heat flows from hot to cold along the rod until the temperature $u(x, t)$ stops changing.

In 1D the Laplacian $\nabla^2 u$ is just $u''(x)$, so the steady-state equation becomes

Integrate twice. The first integration gives $u'(x) = C_1$, a constant slope. The second gives $u(x) = C_1 x + C_2$, a straight line. The two boundary conditions pin $C_2 = T_L$ and $C_1 = (T_R - T_L)/L$, so

The 1D moral. "Steady state" means "the second derivative vanishes." In 1D that forces the profile to be a straight line — there is just one shape it can take. The heat flux $q = -k\,u'(x) = k\,(T_L - T_R)/L$ is also constant: just as much heat enters from the left as leaves from the right, every second.

Two dimensions are richer. The constraint $u_{xx} + u_{yy} = 0$ still forces smoothness, but the field has room to bend — it can be steep in one direction and shallow in the perpendicular one. The rest of this section explores how that extra freedom unfolds.

Why Laplace's Equation Governs Steady-State Heat

Let us be slightly more careful. The 2D heat equation is

Here $u(x, y, t)$ is the temperature at point $(x, y)$ at time $t$, and $\alpha > 0$ is the thermal diffusivity — large for copper, small for wood, ridiculous for ice cream.

Suppose the boundary temperatures $u\bigl|_{\partial \Omega}$ are held constant in time (a heater on one side, an ice bath on the other). Intuitively the interior temperature should approach a unique equilibrium — and it does. As $t \to \infty$, $u(x, y, t) \to u_\infty(x, y)$, and the limit obeys

with the same boundary data. This is a boundary value problem: no initial condition is needed because all memory of the past has decayed away.

There is also a direct, time-independent derivation that is worth carrying in your head. In equilibrium, the net heat flux into any small region must be zero (otherwise its temperature would change and we would not be in equilibrium). Fourier's law says the flux is $\mathbf{q} = -k\,\nabla u$. The divergence theorem turns "net flux through the boundary" into $\int \nabla \cdot \mathbf{q}\,dA$, which equals $-k \int \nabla^2 u\,dA$. Setting this to zero for every region forces $\nabla^2 u = 0$ pointwise. The Laplacian is just the "net flow out minus flow in" per unit area.

The Relaxation Rule and Why It Works

Discretise the plane on a square grid of spacing $h$. The standard second-order finite differences give

Add them and set the sum to zero. The $h^2$ cancels, and we are left with

which rearranges into the celebrated five-point stencil:

The intuition in one sentence: at equilibrium, every cell is the average of its four neighbours. The heat flowing in from above-plus-below must equal the heat flowing out to left-plus-right, otherwise the cell would be heating or cooling.

That algebraic identity is also an algorithm. Pick any initial guess. Sweep through every interior cell and overwrite it with the average of its four current neighbours. Repeat. The field relaxes to the unique steady state. This is called Jacobi iteration, and it is the cousin of gradient descent on the energy functional $E[u] = \tfrac{1}{2}\int |\nabla u|^2\,dA$.

The Mean-Value Property

The discrete rule above is a finite-difference shadow of an exact continuous statement called the mean-value property of harmonic functions:

for every disk of radius $r$ centred at $(x_0, y_0)$ that fits inside the domain. In words: the value at the centre equals the average of the values on the circle around it.

This is one of the most elegant facts in calculus. It says that harmonic functions are perfectly smooth and democratic: no single point can defy its neighbours, because it is forced to be their average. Slide the circle around the field below and watch the equality hold no matter where you put it.

Notice that the equality is exact — the slight numerical difference you see is purely discretisation error from sampling the continuous field on a finite grid.

The Maximum Principle: No Bumps Inside

From the mean-value property, an immediate corollary falls out that engineers rely on every day: the maximum principle.

Maximum principle (heat version). Suppose $u$ is harmonic on a domain $\Omega$. Then the maximum and minimum of $u$ are both attained on the boundary $\partial \Omega$, not in the interior.

Why? If $u$ had an interior local maximum at some point $p$, then $u(p)$ would be strictly greater than every nearby value — but the mean-value property forces $u(p)$ to equal the average of those nearby values. Contradiction. The same argument rules out interior local minima.

Physically: no interior point can be hotter than every single point on the boundary. If your chip is dissipating heat purely through its boundary, and the package surface is at 80°C, then nowhere in the silicon can exceed 80°C in steady state. The only way to get an interior hotspot is to add an internal heat source, which turns Laplace into Poisson's equation $\nabla^2 u = -f$ — the subject of Chapter 33.

Worked Example: 3×3 Grid by Hand

Before letting the computer iterate ten thousand times, do two passes by hand on a tiny grid. This builds enormous intuition for what the solver is actually doing.

Initial field:
 100  100  100  100  100
   0   25   25   25    0
   0   25   25   25    0
   0   25   25   25    0
   0    0    0    0    0

After iteration 1:
 100   100      100     100    100
   0    37.50   43.75   37.50    0
   0    18.75   25.00   18.75    0
   0    12.50   18.75   12.50    0
   0     0       0       0       0

After iteration 2:
 100    100      100      100    100
   0     40.625   50.000   40.625   0
   0     18.750   25.000   18.750   0
   0      9.375   12.500    9.375   0
   0      0        0         0      0

Interactive Solver

Here is the same algorithm scaled up to a 60×60 grid. Drag the four boundary sliders, press Play to watch the field relax, or step one iteration at a time. Click anywhere on the field to read off the temperature at that point.

A few experiments to do with this visualiser:

• Set the four sides to 100, 0, 100, 0. The diagonal symmetry forces the centre to be exactly 50 — confirm by clicking the middle.
• Set three sides to 0 and the fourth to 100. The centre converges to 25 (the rotation-symmetry argument from the worked example).
• Make all four sides equal. The solution is a constant equal to that value — the only harmonic function with constant Dirichlet data on a closed boundary.
• Watch the residual. It decays roughly geometrically, which reflects the slow eigenvalue of the Jacobi iteration matrix on a square grid (about $1 - \pi^2/(2N^2)$).

Python: Solving with Jacobi Iteration

Now let us write the same algorithm explicitly. Plain Python with NumPy makes every step transparent. Read every line carefully — the entire body of the function is just average your neighbours, repeated until nothing changes.

(100 + 0 + 0 + 0) / 4 = 25.0  →  every interior cell starts at 25.0

np.max(np.abs(u_new - u))  →  shrinks roughly like 0.97^k for this geometry

1import numpy as np
2
3def steady_state_heat(top, bottom, left, right, N=50, tol=1e-5, max_iter=20000):
4    """
5    Solve Laplace's equation  u_xx + u_yy = 0  on a unit square
6    with Dirichlet boundary conditions (constant temperature on each side)
7    via Jacobi iteration.
8
9    Returns the converged temperature field of shape (N, N).
10    """
11    u = np.full((N, N), (top + bottom + left + right) / 4, dtype=np.float64)
12
13    # Apply Dirichlet boundary conditions (these never change)
14    u[0, :]   = top      # row 0 is the top edge
15    u[-1, :]  = bottom
16    u[:, 0]   = left
17    u[:, -1]  = right
18
19    for k in range(max_iter):
20        # Average each interior cell against its four neighbours
21        u_new = u.copy()
22        u_new[1:-1, 1:-1] = 0.25 * (
23            u[:-2, 1:-1] + u[2:, 1:-1] + u[1:-1, :-2] + u[1:-1, 2:]
24        )
25        residual = np.max(np.abs(u_new - u))
26        u = u_new
27        if residual < tol:
28            print(f"Converged in {k + 1} iterations  (residual = {residual:.2e})")
29            return u
30
31    print(f"Stopped at max_iter ({max_iter}). Residual = {residual:.2e}")
32    return u
33
34
35# Try the classical "100° on top, 0° elsewhere" problem.
36u = steady_state_heat(top=100, bottom=0, left=0, right=0, N=51)
37
38# By a beautiful symmetry argument, the centre temperature is exactly 25°:
39# four superimposed problems (one hot side each) sum to a constant 100° field,
40# so by linearity each one contributes 25° at the centre.
41print(f"u(centre) = {u[25, 25]:.4f}      # should be very close to 25.0000")

Running this for $N = 51$ converges in about $2{,}500$ iterations to a residual of $10^{-5}$. The reported centre temperature is 25.0000 to four decimals — exactly the symmetry prediction.

PyTorch: The Same Idea, Tensorized

The five-point stencil is a tiny convolution. That observation makes the GPU port absurdly short — we describe the averaging rule as a 3×3 kernel and let torch.nn.functional.conv2d do the work. The boundary still needs to be re-imposed by hand because conv2d does not know it is special.

\frac{1}{4}\,\bigl(u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1}\bigr)

1import torch
2import torch.nn.functional as F
3
4def steady_state_heat_torch(top, bottom, left, right, N=128, tol=1e-5, max_iter=20000,
5                            device="cuda" if torch.cuda.is_available() else "cpu"):
6    """
7    Solve Laplace's equation on a unit square with constant Dirichlet BCs
8    using PyTorch tensors and a 3x3 averaging convolution.
9
10    Identical math to the NumPy version, but every iteration is one fused
11    GPU kernel — orders of magnitude faster for large N.
12    """
13    seed = (top + bottom + left + right) / 4.0
14    u = torch.full((1, 1, N, N), seed, dtype=torch.float32, device=device)
15
16    # Pin Dirichlet boundary conditions
17    u[..., 0,  :]  = top
18    u[..., -1, :]  = bottom
19    u[..., :,  0]  = left
20    u[..., :, -1]  = right
21
22    # Five-point stencil  =  3x3 conv kernel that picks N, S, E, W and averages
23    kernel = torch.tensor(
24        [[0.0, 0.25, 0.0],
25         [0.25, 0.0, 0.25],
26         [0.0, 0.25, 0.0]],
27        dtype=torch.float32, device=device,
28    ).view(1, 1, 3, 3)
29
30    for k in range(max_iter):
31        # Convolve with stride 1, padding 1 to keep the same shape
32        u_new = F.conv2d(u, kernel, padding=1)
33
34        # Re-apply boundary conditions (the conv would have rewritten them)
35        u_new[..., 0,  :]  = top
36        u_new[..., -1, :]  = bottom
37        u_new[..., :,  0]  = left
38        u_new[..., :, -1]  = right
39
40        residual = (u_new - u).abs().max().item()
41        u = u_new
42        if residual < tol:
43            print(f"Converged in {k + 1} iterations  (residual = {residual:.2e})")
44            return u.squeeze().cpu().numpy()
45
46    return u.squeeze().cpu().numpy()
47
48
49u = steady_state_heat_torch(top=100, bottom=0, left=0, right=0, N=129)
50print("u(centre) =", float(u[64, 64]))   # ≈ 25.0  (same 4-side superposition rule)

For $N = 1024$ the PyTorch version finishes in seconds on a consumer GPU while NumPy takes minutes. The algorithm is identical — the only thing that changed is who is doing the arithmetic.

Why this matters for ML practitioners. The same averaging kernel appears in image processing as a blur, in graph neural networks as a diffusion step, and in diffusion models as the score-matching denoiser. They are all discrete approximations of solutions to Laplace-like PDEs. Recognising the structure in one place teaches you the others.

Where This Shows Up in the Real World

Steady-state heat models are everywhere. A few representative examples:

The image-inpainting connection (a small detour)

Image inpainting — filling in a missing region of a photograph — is mathematically identical to a steady-state heat problem. Treat the known pixels as boundary conditions and let the missing pixels relax to the harmonic interpolant. The Jacobi iteration we just wrote, applied to a colour image with a hole cut out, is literally Poisson image editing without sources. Adding sources gives you Photoshop's seamless-cloning feature, which is Poisson's equation — the next chapter.

Summary

1. Steady-state heat distributions are exactly the harmonic functions: solutions to $\nabla^2 u = 0$ with prescribed boundary temperatures.
2. The discrete five-point rule $u_{i,j} = \tfrac{1}{4}\sum_{\text{neighbours}} u$ is both the discrete Laplacian-equals-zero condition and the iterative algorithm that finds it.
3. Harmonic functions obey the mean-value property — the centre of any disk equals the average around the rim — which forces all extrema to live on the boundary.
4. The same algorithm runs in NumPy (CPU) and as a 3×3 convolution in PyTorch (GPU). The math is identical; only the data plumbing changes.
5. Recognising harmonic interpolation in disguise — heat sinks, image inpainting, terrain shading, graph smoothing — turns a single PDE into a dozen practical tools.