Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

State Laplace's equation $\Delta u = 0$ and explain, in physical words, what it asserts about a function on a domain.
Identify whether a given function is harmonic by computing its Laplacian by hand and by code.
Use the mean value property as a self-contained definition of harmonicity, and explain why it forces the maximum principle.
Derive the discrete Laplace equation on a grid and recognise that Jacobi iteration is repeated averaging.
Solve a small boundary-value problem on a 4 × 4 grid by hand, then watch a 51 × 51 PyTorch implementation reproduce the same field at scale.

The Big Picture: Calculus of Equilibrium

Most equations in calculus describe change: velocity, growth, the propagation of heat over time. Laplace's equation is the opposite. It is the equation that survives once all the change has died down — the calculus of what remains after every wave has stilled, every diffusion has finished, every flow has reached its destination. It is the calculus of equilibrium.

The intuition. Imagine pressing a thin rubber sheet onto a wire frame bent into some shape. Once it stops vibrating, the sheet's height $u(x, y)$ at every interior point is determined entirely by where the frame is. The sheet has the smoothest possible shape consistent with those boundary heights. That shape is the harmonic function. Laplace's equation is the mathematical statement that no interior point is "pulling harder" than its neighbours.

The same equation governs the steady-state temperature in a metal plate, the electrostatic potential between charged plates, the velocity potential of a smooth fluid flow, the pressure in a slow liquid through porous rock, and the displacement of a stretched membrane. Five wildly different physical situations — same equation, same solutions. That kind of universality is rare and worth paying attention to.

The one-line takeaway

Laplace's equation says: at every interior point, the value of $u$ is the average of the values in a tiny ball around it. There are no peaks, no valleys, no local extrema in the interior — only the smoothest possible shape that connects the boundary to itself.

From the Heat Equation to Laplace

We met the heat equation in Chapter 26:

\dfrac{\partial u}{\partial t} \;=\; \alpha \, \Delta u, \qquad \Delta u \;=\; \dfrac{\partial^2 u}{\partial x^2} + \dfrac{\partial^2 u}{\partial y^2}.

It says that the rate of change of temperature at a point is proportional to how much the local temperature differs from the average of its immediate neighbours (the Laplacian measures exactly that difference). Hot spots cool, cold spots warm — and only if there is a difference.

Now ask: what happens after a very long time? If the boundary temperatures are held fixed, the interior will eventually settle into a state where nothing changes anymore. At that point:

\dfrac{\partial u}{\partial t} \;=\; 0 \quad \Longrightarrow \quad \alpha \, \Delta u \;=\; 0 \quad \Longrightarrow \quad \Delta u \;=\; 0.

That last equation is Laplace's equation. It is what you get when you ask "where does heat flow eventually stop?" It is the heat equation with the time derivative deleted. It is also what you would get by asking the same question of the wave equation, the diffusion equation, or any other equation whose right-hand side is the Laplacian.

Why this matters

Every parabolic equation (heat) and every hyperbolic equation (wave) has a Laplacian on the right-hand side. Their steady-state versions are all the same elliptic equation. That is why the small body of theory in this chapter covers so many physical phenomena: we are solving the universal steady-state problem, not any one of them in particular.

The Laplace Operator

Before stating the equation formally we need to know exactly what the symbol $\Delta$ means. In two dimensions:

\Delta u \;=\; \nabla \cdot \nabla u \;=\; \dfrac{\partial^2 u}{\partial x^2} + \dfrac{\partial^2 u}{\partial y^2}.

It is the sum of the unmixed second partial derivatives. In three dimensions you add a $\partial^2 u/\partial z^2$ term. In $n$ dimensions you sum $n$ terms. No cross derivatives like $u_{xy}$ ever appear — the Laplacian is invariant under rotation, and cross terms would break that.

Two readings of the Laplacian

Analytic: sum of the curvatures along each axis.

Geometric (more useful for intuition): $\Delta u(x_0) \;\approx\; \dfrac{2n}{r^2} \bigl( \overline{u}_{\text{ball}}(x_0, r) - u(x_0) \bigr)$ — proportional to the average value of u in a tiny ball around $x_0$ minus the value at the centre.

That second reading is the key to everything. The Laplacian measures the local discrepancy between a point's value and its neighbourhood average. Saying $\Delta u = 0$ is therefore the same as saying "every interior point equals the local average."

Laplace's equation: $\Delta u \;=\; 0 \quad \text{in } \Omega, \qquad u \;=\; g \quad \text{on } \partial \Omega.$ Find a function on the domain $\Omega$ whose Laplacian vanishes in the interior and that takes prescribed values $g$ on the boundary.

Harmonic Functions: The Solutions

A function $u$ on an open set is called harmonic if $\Delta u = 0$ at every point. Here are the most important harmonic functions:

Function	Laplacian	Why harmonic
u(x, y) = c (constant)	0 + 0 = 0	Trivially smooth
u(x, y) = ax + by + c	0 + 0 = 0	Planes have no curvature
u(x, y) = x² − y²	+2 + (−2) = 0	Curvatures cancel — the saddle
u(x, y) = xy	0 + 0 = 0	Rotated saddle
u(x, y) = e^x cos(y)	e^x cos(y) − e^x cos(y) = 0	Real part of e^z
u(x, y) = ln√(x² + y²)	0 (off origin)	The 2D fundamental solution

Notice the pattern. The polynomials are $1, x, y, x^2 - y^2, xy, x^3 - 3xy^2, 3x^2 y - y^3, \dots$ — these are exactly the real and imaginary parts of $(x + iy)^n$ . That is no coincidence: in 2D, every analytic complex function $f(z) = u + i v$ has $u$ and $v$ both harmonic, and conversely every harmonic $u$ arises as the real part of some complex analytic function. Complex analysis and 2D harmonic analysis are the same subject in disguise.

Why "harmonic"?

The name comes from acoustics. The normal modes of a vibrating drum on a domain $\Omega$ satisfy $-\Delta u = \lambda u$ . The simplest modes are the harmonics — pure tones. The static equilibrium case corresponds to $\lambda = 0$ : zero pitch, zero oscillation, just the silent steady shape. Those silent shapes are the "harmonics" we still call by that name.

Gallery: Harmonic, or Not?

The visual signature of a harmonic function is striking once you learn to see it. Below, six functions are plotted on the left and their Laplacians on the right. A harmonic function has a Laplacian that is uniformly the middle color of the colormap — flat zero everywhere. A non-harmonic one has hot or cold regions in its Laplacian, marking the points where the local average disagrees with the function's value.

Two takeaways. First: the harmonic surfaces all have a characteristic look — they are smooth, gently undulating, and almost "flat-looking" in the middle. They lack the dome-shape of the bowl $x^2 + y^2$ and the volcano-shape of the Gaussian. Second: the Laplacian of a non-harmonic function tells you, point by point, by how much the equation is violated. That residual is the central quantity in any numerical solver.

The Mean Value Property

We have asserted twice already that a harmonic function equals its local average. Now we state it precisely.

Mean Value Property

If $u$ is harmonic on an open set containing a closed disk $\overline{D}(x_0, r)$ , then:

u(x_0) \;=\; \dfrac{1}{2\pi r} \oint_{\partial D} u \, ds \;=\; \dfrac{1}{\pi r^2} \iint_{D} u \, dA.

The value at the centre equals both the average over the bounding circle and the average over the entire disk — for every centre and every radius for which the disk fits inside the harmonic domain.

Why is this true? Intuitively: harmonicity says the local average equals the value at the point in the limit of an infinitesimal ball. Integrating that infinitesimal statement out from radius $0$ to radius $r$ gives the finite-radius version. The clean proof uses the divergence theorem:

One-line proof sketch

Define $M(r) = \tfrac{1}{2\pi r}\oint_{\partial D(x_0, r)} u \, ds$ . Differentiating under the integral and applying the divergence theorem gives $M'(r) = \tfrac{1}{2\pi r} \iint_{D(x_0, r)} \Delta u \, dA$ . If $\Delta u = 0$ , then $M'(r) = 0$ , so $M(r)$ is constant in $r$ . Taking $r \to 0$ by continuity shows the constant value is $u(x_0)$ .

The converse is also true and almost as important: if $u$ is continuous and satisfies the mean value property on every small disk, then $u$ is harmonic. So you can use $\Delta u = 0$ and "value = local average" interchangeably as the definition of harmonic — they are the same condition wearing different clothes.

Interactive: The Mean Value Property

Drag the yellow centre, change the radius. The plotted field is $u(x, y) = x^2 - y^2$ — harmonic everywhere. The numerical average around the circle is computed live and compared with the value at the centre. They agree to round-off for every choice you make.

Things to try

Shrink the radius. The numerical match stays exact — radius is irrelevant for harmonic functions.
Move the centre near (1, 0). The point value goes up, the circle average rises with it.
Move along the line $y = x$ . Both values become zero — and stay zero — because $x^2 - y^2 = 0$ on that line.

The Maximum Principle

The mean value property has an immediate, deeply useful consequence:

Maximum (and Minimum) Principle

A non-constant harmonic function on a bounded domain attains its maximum and minimum only on the boundary. There are no interior peaks or valleys.

Why? Suppose a harmonic function had a maximum at some interior point $x_0$ . Then $u(x_0) \geq u(y)$ for every nearby point $y$ . But the mean value property says $u(x_0)$ equals the average of $u$ over a small circle around $x_0$ . The average can equal the maximum only if every value on the circle equals the maximum. By repeating the argument outward, $u$ would have to be constant on a connected open set — and by analytic continuation, everywhere. Contradiction.

Concrete consequence: if a metal plate's edges are all held between 0°C and 100°C, then no interior point is hotter than 100°C or colder than 0°C. Ever. The interior is squeezed by its boundary. This is the simplest non-trivial fact about Laplace's equation, and it underwrites almost everything else.

A direct corollary: uniqueness

If $u_1$ and $u_2$ are two harmonic functions satisfying the same boundary conditions, then $v = u_1 - u_2$ is harmonic and zero on the boundary. By the maximum/minimum principle, the maximum and minimum of $v$ on the closed domain are both 0. So $v \equiv 0$ , meaning $u_1 = u_2$ . The boundary uniquely determines the interior.

The Discrete Laplace Equation

To compute, we replace the continuum operator $\Delta$ by its finite-difference cousin. On a uniform grid with spacing $h$ :

\Delta u(x_i, y_j) \;\approx\; \dfrac{u_{i+1, j} + u_{i-1, j} + u_{i, j+1} + u_{i, j-1} - 4 \, u_{i, j}}{h^2}.

Setting this to zero and solving for $u_{i,j}$ :

u_{i, j} \;=\; \tfrac{1}{4} \bigl( u_{i+1, j} + u_{i-1, j} + u_{i, j+1} + u_{i, j-1} \bigr).

This is the discrete mean value property — and the basis of every elementary numerical method for Laplace's equation. It says: in the discrete picture, every interior cell's value must equal the average of its four neighbours. Solving Laplace's equation on a grid means making this true at every interior cell simultaneously.

The simplest way to make it true is to iterate. Pick any guess. Compute, at every interior cell, the average of its four current neighbours. Replace the guess with the averages. Do it again. And again. This is the Jacobi iteration:

u^{(k+1)}_{i, j} \;=\; \tfrac{1}{4} \bigl( u^{(k)}_{i+1, j} + u^{(k)}_{i-1, j} + u^{(k)}_{i, j+1} + u^{(k)}_{i, j-1} \bigr).

Each sweep replaces every cell by the average of its old neighbours. The boundary cells are held fixed. After enough sweeps, nothing changes anymore — and the discrete Laplace equation is satisfied everywhere.

Why does it converge?

The Jacobi update is a contractive averaging operator on the interior values. Each sweep reduces the maximum-norm error between the current iterate and the true solution by a factor like $\cos(\pi/N)$ on an $N \times N$ grid — so the error decays geometrically. The rate is slow (close to 1 for large $N$ ), which is why production solvers use multigrid or conjugate-gradient methods. But Jacobi is the one whose every step looks like the mean value property, so it is by far the most pedagogical.

Interactive: Solving Laplace on a Plate

Choose a boundary pattern. Hit play. Watch the interior fill in. Every frame is one or more Jacobi sweeps — each interior cell becomes the average of its four neighbours. The colormap tracks the temperature; the residual readout shows how close we are to the true steady state. Hit Jump to steady state to run thousands of sweeps in one click and see the final harmonic field.

What to look for

With Hot top edge, the field decays smoothly from top to bottom — no overshoots, no oscillations, no bumps. A harmonic field is the world's most boring interpolant.
With Two hot sides, the centre value approaches the average of the boundary: $u_{\text{centre}} \approx 0.5$ . The field bows outward toward the cold top and bottom.
With Sinusoidal top, the height of the wave drops off exponentially as you move down. This is the separation-of-variables eigenfunction we will derive in Section 28.3.
With Hot corners, the boundary has sharp jumps. The harmonic interior smooths every jump out — within a few cells, you cannot tell there was a discontinuity at all.

Worked Numerical Example

Let's solve the smallest non-trivial Laplace problem by hand, then watch Jacobi converge to the same answer. This is the kind of calculation you should be able to reproduce on paper.

The 4 × 4 plate problem — full hand calculation

Setup

A 4 × 4 grid with indices $(i, j) \in \{0, 1, 2, 3\}^2$ . Row $i = 0$ is the top of the plate, row $i = 3$ the bottom. Boundary values:

Top row ( $i = 0$ ): u = 100
Bottom row ( $i = 3$ ): u = 0
Left column ( $j = 0$ ): u = 0
Right column ( $j = 3$ ): u = 0

The four interior unknowns are $a = u_{1,1}, \; b = u_{1,2}, \; c = u_{2,1}, \; d = u_{2,2}$ .

Step 1 — write the discrete Laplace equation at each unknown

$4 a \;=\; u_{0,1} + u_{2,1} + u_{1,0} + u_{1,2} \;=\; 100 + c + 0 + b$

$4 b \;=\; u_{0,2} + u_{2,2} + u_{1,1} + u_{1,3} \;=\; 100 + d + a + 0$

$4 c \;=\; u_{1,1} + u_{3,1} + u_{2,0} + u_{2,2} \;=\; a + 0 + 0 + d$

$4 d \;=\; u_{1,2} + u_{3,2} + u_{2,1} + u_{2,3} \;=\; b + 0 + c + 0$

Step 2 — exploit symmetry

The boundary is symmetric under $j \mapsto 3 - j$ (left/right mirror), so the solution is too. That forces $a = b$ and $c = d$ . Two unknowns left.

Step 3 — substitute and solve

Using $a = b, c = d$ :

$4 a = 100 + c + a \;\;\Longrightarrow\;\; 3 a = 100 + c.$

$4 c = a + c \;\;\Longrightarrow\;\; 3 c = a \;\;\Longrightarrow\;\; c = a/3.$

Plug the second into the first:

$3 a = 100 + a/3 \;\;\Longrightarrow\;\; 9 a - a = 300 \;\;\Longrightarrow\;\; 8 a = 300 \;\;\Longrightarrow\;\; a = 37.5.$

$c = a/3 = 12.5.$

The exact answer

$u_{1,1} = u_{1,2} = 37.5, \quad u_{2,1} = u_{2,2} = 12.5.$ The top half of the interior is at 37.5°, the bottom half at 12.5°. The vertical decay is dramatic — 37.5 → 12.5 — because the cold left, right, and bottom boundaries pull the field down fast.

Sanity check — verify by plugging back

Equation at $a$ : $4 \cdot 37.5 = 150$ and $100 + 12.5 + 0 + 37.5 = 150$ . ✓

Equation at $c$ : $4 \cdot 12.5 = 50$ and $37.5 + 0 + 0 + 12.5 = 50$ . ✓

Step 4 — watch Jacobi converge to it

Starting from all zeros, the first four sweeps give:

Sweep	a = b	c = d	max change
0 (init)	0.000	0.000	—
1	25.000	0.000	25.000
2	31.250	6.250	6.250
3	34.375	9.375	3.125
4	35.938	10.938	1.562
10	37.451	12.451	0.049
20	37.500	12.500	5e−5
∞ (true)	37.5	12.5	0

Each sweep cuts the residual in half (roughly). After 20 sweeps the solution is correct to four decimal places. This is the entire algorithm — averaging, repeatedly — and it converges to the unique harmonic function with the prescribed boundary.

Python: Verifying a Function Is Harmonic

Two checks in one script. First we confirm that $u(x, y) = x^2 - y^2$ is harmonic by numerically computing $u_{xx} + u_{yy}$ . Then we verify the mean value property by averaging $u$ over a circle and comparing to the value at the centre. Read the explanations on the right — each annotation maps to exactly one line.

🐍python

Explanation(18)

Code(43)

1NumPy import

NumPy gives us vectorised arithmetic (np.cos, np.mean, np.linspace) so the circle-average integral becomes one short expression. No PyTorch yet — the goal here is to verify the definition of 'harmonic' by hand-derivable math, not to train anything.

5The candidate function u(x, y) = x² − y²

This is the prototypical 2D harmonic function — the 'hyperbolic saddle.' Its level sets are hyperbolas xy = const after a 45° rotation. We can compute its partials by hand: u_x = 2x, u_xx = 2, u_y = −2y, u_yy = −2. Their sum is 2 + (−2) = 0. We want NumPy to confirm this numerically before we trust anything more elaborate.

10Centred finite difference for u_xx

The exact second derivative is u_xx(x) = lim_{h→0} [u(x+h, y) − 2u(x, y) + u(x−h, y)] / h². We replace the limit by a small but nonzero h = 1e−3. Why centred? Because centred differences cancel the odd-order terms in the Taylor expansion, so the leading error is O(h²) rather than O(h). For this smooth polynomial that error is essentially floating-point noise (~1e−10), which is well below the values we care about (~2.0).

13Centred finite difference for u_yy

Same idea, but holding x fixed and varying y. Notice we are doing pure partials — no cross terms u_xy enter the Laplacian, only the two pure second derivatives along the axes. That is why Δ = ∂²/∂x² + ∂²/∂y² and not something more complicated.

16Sample point in the plane

We pick a generic interior point (0.5, 0.3) — not on any axis, not at the origin, nothing special. If a function is harmonic, the Laplacian must vanish at EVERY point in the open domain, so picking a 'random' point is a stronger test than picking a symmetric one.

18Compute u_xx at the sample point

Numerical value should be ≈ +2.0 (the exact second derivative of x² − y² in the x direction). Tiny deviations from 2.0 are the centred-difference truncation error and floating-point round-off, on the order of 1e−9.

19Compute u_yy at the sample point

Should be ≈ −2.0. The minus sign appears because y enters the function as −y². Geometrically, the surface curves upward as you walk in the x direction and downward as you walk in the y direction — those two curvatures are equal and opposite.

20Add them — the Laplacian

By definition Δu = u_xx + u_yy. Here 2 + (−2) = 0. The two curvatures cancel. The function passes through every point with zero net second-derivative — every point is the average of its neighborhood. That is precisely what 'harmonic' means.

22Print the three numerical values

We print all three so a reader can see the cancellation happen, not just the final 0. If the Laplacian is essentially zero (say, smaller than 1e−6) but the two pieces are not, that is the signature of a harmonic function — local curvatures exist, but they balance.

28circle_average — the integral as a Monte-Carlo-like sum

For the mean value property we need (1/2π) ∫₀^{2π} u(cx + r cos t, cy + r sin t) dt. A trapezoidal sum on n equispaced points is spectrally accurate for smooth periodic integrands like ours, so a thousand samples gives near-machine precision. We will compare this average to the value at the centre.

29n_samples = 1024 — why so many?

For a non-periodic integrand we would worry about discretisation error scaling like 1/n. But on a closed curve the periodic trapezoid rule converges exponentially fast in n (the error decays faster than any polynomial). 1024 is overkill — you can get the same accuracy with 64 — we just want zero ambiguity in the printed match.

30np.linspace(0, 2π, n, endpoint=False)

endpoint=False is important: a periodic integrand evaluated at both 0 and 2π would count the same point twice. Dropping the right endpoint gives a clean uniform partition of the circle.

31Parameterise the circle

A circle of radius r around (cx, cy) is (cx + r cos t, cy + r sin t) for t ∈ [0, 2π). Vectorising over t means we evaluate u at all 1024 boundary points in one NumPy call — much faster than a Python loop.

33Average via np.mean

np.mean computes (1/n) Σ u(sample_i). With equispaced samples on a closed curve this IS the trapezoidal estimate of (1/2π) ∮ u ds / r (after factoring out the constant arc-length element). It is the numerical face of the line-integral definition of 'mean value on a circle.'

36Pick a radius r = 0.7

Radius doesn't need to be tiny. The mean value property is exact for ANY radius, as long as the closed disk fits inside the domain where u is harmonic. We pick 0.7 to make it visible in the interactive viewer above (where the domain is the 2×2 square).

37centre_value = u(x0, y0)

By hand: u(0.5, 0.3) = 0.25 − 0.09 = 0.16. We will compare this to the circle average.

38Numerically compute the circle average

All 1024 evaluations of u on the boundary circle, then divide by 1024. Because the integrand x² − y² is a polynomial in cos t and sin t and the integral of cos t, sin t, cos t sin t over a full period is zero, the only surviving piece is (cx² − cy²) + r² · (½ − ½) = cx² − cy² = u(cx, cy). The match has to be exact analytically.

41Print and inspect

Expected output: u(0.5, 0.3) = +0.160000 and average ≈ +0.160000 with a difference around 1e−14 — pure round-off. The mean value property is now confirmed for this point and radius. It would be confirmed equally well for any other (cx, cy) and any r ≤ the distance to the boundary.

25 lines without explanation

1import numpy as np
2
3# Candidate function:  u(x, y) = x^2 - y^2
4# (a classic harmonic function — the "hyperbolic saddle".)
5def u(x, y):
6    return x ** 2 - y ** 2
7
8# Centred finite-difference second derivative.
9# u_xx(x, y) ~ [u(x + h, y) - 2 u(x, y) + u(x - h, y)] / h^2
10def u_xx(x, y, h=1e-3):
11    return (u(x + h, y) - 2 * u(x, y) + u(x - h, y)) / h ** 2
12
13def u_yy(x, y, h=1e-3):
14    return (u(x, y + h) - 2 * u(x, y) + u(x, y - h)) / h ** 2
15
16# Sample point in the interior of the plane.
17x0, y0 = 0.5, 0.3
18
19uxx = u_xx(x0, y0)
20uyy = u_yy(x0, y0)
21laplacian = uxx + uyy
22
23print(f"u_xx(x0, y0)        = {uxx:+.6f}")
24print(f"u_yy(x0, y0)        = {uyy:+.6f}")
25print(f"Laplacian (u_xx+u_yy) = {laplacian:+.6f}")
26
27# Mean value property — numerical check.
28# For ANY harmonic function, the value at the centre of a disk
29# equals the average of the function around the boundary circle.
30def circle_average(cx, cy, r, n_samples=1024):
31    t = np.linspace(0, 2 * np.pi, n_samples, endpoint=False)
32    xs = cx + r * np.cos(t)
33    ys = cy + r * np.sin(t)
34    return float(np.mean(u(xs, ys)))
35
36r = 0.7
37centre_value = u(x0, y0)
38average      = circle_average(x0, y0, r)
39
40print()
41print(f"u(x0, y0)              = {centre_value:+.6f}")
42print(f"average on circle r=0.7 = {average:+.6f}")
43print(f"difference              = {abs(centre_value - average):.2e}")

Expected output

u_xx(x0, y0)        = +2.000000
u_yy(x0, y0)        = -2.000000
Laplacian (u_xx+u_yy) = +0.000000

u(x0, y0)              = +0.160000
average on circle r=0.7 = +0.160000
difference              = 4.16e-15

The Laplacian is zero to round-off precision. The circle average equals the centre value to fourteen digits. Both fingerprints of a harmonic function are present.

PyTorch: Jacobi Relaxation on a Grid

Now we scale from a 4 × 4 hand calculation to a 51 × 51 grid. The whole Jacobi sweep becomes a single tensor expression — slice the array four ways, average. PyTorch handles 2401 simultaneous averages in one line and runs the same code on CPU or GPU.

🐍python

Explanation(20)

Code(40)

1torch import

PyTorch lets us write the entire Jacobi sweep as one tensor expression — no inner Python loop over (i, j). The same code will run on CPU or GPU just by changing dtype/device. We use float64 throughout because the convergence test compares values to a tolerance of 1e−5, and float32 round-off (~1e−7) is uncomfortably close to that for residuals that decay slowly.

6N = 51 — grid resolution

A 51 × 51 grid has 51 boundary points on each edge and 49 × 49 = 2401 interior unknowns. Why an odd N? Because then the geometric centre is exactly at (N // 2, N // 2) = (25, 25), and we can print 'u at the centre' without any half-cell interpolation.

7TOL = 1e−5 — convergence threshold

We stop when no grid point changed by more than 1e−5 in one sweep. The discrete Laplace residual is bounded by this number, so when TOL is hit the discrete equation Δ_h u = 0 holds at every interior point to five decimal places. Tightening to 1e−8 multiplies the iteration count by about 1000 — Jacobi is famously slow.

8MAX_ITER safety cap

A hard upper bound prevents infinite loops if convergence is slow. For a 51 × 51 grid Jacobi needs roughly O(N²) ≈ 2500 sweeps to reach 1e−5. Leaving headroom (20 000) means we will not abort mid-convergence by accident.

11Initial guess u = zeros

The Laplace equation has a UNIQUE solution given the boundary values, so the initial guess does not affect the final answer — only the speed of convergence. Zero is the simplest choice and matches our intuition that we know nothing about the interior yet.

14u[0, :] = 0.0 — bottom row

Convention: row index 0 is the bottom of the plate, row N − 1 is the top. The entire bottom edge is held at 0°. This is a Dirichlet boundary condition: prescribed value, not prescribed derivative. After this assignment, those values must NEVER be overwritten in the loop, which is why we update only the interior slice u[1:-1, 1:-1] below.

15u[-1, :] = 100.0 — hot top row

The top edge is held at 100°. This is the only boundary doing any 'pushing' — the other three edges are sinks at 0°. Intuitively, heat will flow inward from the top and trickle down toward the cold edges, with the steady state being a smooth interpolation. The exact shape is what the Laplace equation determines.

16u[:, 0] = 0.0 — left column

Cold edge. Together with the bottom and right edges, three of the four sides act as a heat sink at 0°. Symmetry note: the final solution will be symmetric under j → (N − 1 − j) (mirror across the vertical centre line) because the boundary conditions are.

17u[:, -1] = 0.0 — right column

Cold edge. With these four assignments, the boundary problem is fully specified. There is now exactly one harmonic function on the open interior that matches these values at every boundary point — that is the function we are about to compute.

23for step in range(MAX_ITER)

Outer loop: each iteration is one Jacobi sweep, which updates every interior point once. We will break out as soon as the residual falls below TOL. Note: this is NOT a time-stepping loop — there is no physical time in Laplace's equation. The iteration index is an algorithmic relaxation, not a clock.

24u_new = u.clone()

Why clone? Because Jacobi is meant to read ONLY the old values when computing the new ones. If we wrote into u in-place, half the neighbours would already be new and half would be old — that is the (different, also valid) Gauss-Seidel method. Jacobi is the pure 'average of old neighbours' update, so it requires two buffers.

25The big tensor expression — the heart of the solver

Read it as four shifted views of the array, added together: u[2:, 1:-1] is 'one cell up from each interior cell,' u[:-2, 1:-1] is 'one cell down,' u[1:-1, 2:] is 'one cell right,' u[1:-1, :-2] is 'one cell left.' Their sum, divided by 4, is the average of the 4 neighbours. We assign this average to the interior slice u_new[1:-1, 1:-1]. The boundaries remain untouched.

26Slicing: u[2:, 1:-1]

u[2:, ·] starts at row 2 and goes to the top — that is precisely the 'one row above' shift for the interior. The column slice 1:-1 keeps only the interior columns. So u[2:, 1:-1] has shape (N − 2, N − 2), the same shape as u_new[1:-1, 1:-1]. PyTorch broadcasts the addition cell-by-cell.

29diff = (u_new − u).abs().max() — the residual

For Laplace, the discrete equation is u_{i,j} = ¼(neighbours). If u solves it, then u_new = u exactly and diff = 0. The biggest cell-level change measures how badly the equation is currently violated. .item() pulls the scalar out of the 0-D tensor so we can compare it in pure Python.

30u = u_new — swap buffers

u is now the post-sweep state; u_new becomes a fresh blank in the next iteration via .clone(). For large grids you would re-use a pre-allocated buffer with .copy_() instead of allocating a new tensor every step — but for clarity we are not doing that here.

31if diff < TOL: break

Once no point has moved by more than 1e−5, the discrete Laplace equation is satisfied at every interior point to that precision. Continuing further would only sharpen agreement past the level we care about. For a 51 × 51 grid with this boundary, expect this to fire near step 2000.

33Convergence printout

The first two prints tell you whether the solver was healthy. If 'converged in' shows MAX_ITER, the residual never fell below TOL — usually means TOL was too tight, or there is a bug in the boundary handling that prevents the iteration from settling.

34u at the centre

For these four boundaries (top = 100, others = 0), symmetry alone tells us u(centre) = (100 + 0 + 0 + 0) / 4 = 25.0 exactly. This is a textbook consequence of the mean value property generalised to the centre of a square: any harmonic function on a domain symmetric under the dihedral group has its centre value equal to the AVERAGE of the boundary values (weighted by arc length). 25.0 is therefore the analytic 'correct answer' that the solver had to find on its own.

35u just below the top edge

Row N − 2 sits one cell below the hot edge. We expect this to be close to but less than 100 — the boundary is at 100, but the value of the harmonic interpolant decays as you move into the interior. Typical value on a 51 × 51 grid: about 75 at the column centre, dropping toward 50 as you approach the corners (where the cold side edges drag it down).

36u just above the bottom edge

Row 1 sits one cell above the cold bottom edge. Expected to be small but positive — the bottom edge is 0, but heat trickling down from the top is still felt faintly at the cell adjacent to the bottom. Typical value: about 4-5 at the column centre. The asymmetry between this number and 'u just below top edge' is the visual signature of a harmonic field — it decays exponentially away from its hot boundary.

20 lines without explanation

1import torch
2
3# Steady-state heat distribution on a square plate.
4# Top edge held at u = 100, all other edges at u = 0.
5# We want u(i, j) for every interior grid point such that
6# Laplacian(u) = 0 at every interior point.
7
8N = 51              # 51 x 51 grid → 49 x 49 interior points
9TOL = 1e-5          # stop when no point moves by more than this
10MAX_ITER = 20_000
11
12# Initial guess: zeros everywhere.
13u = torch.zeros(N, N, dtype=torch.float64)
14
15# Boundary conditions (held fixed every iteration).
16u[0, :]  = 0.0      # bottom row
17u[-1, :] = 100.0    # top row
18u[:, 0]  = 0.0      # left column
19u[:, -1] = 0.0      # right column
20
21# One Jacobi sweep is exactly the discrete Laplace equation
22# rearranged for u_{i,j}:
23#   4 u_{i,j} = u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1}
24# Solve for u_{i,j}, apply to every interior point at once.
25for step in range(MAX_ITER):
26    u_new = u.clone()
27    u_new[1:-1, 1:-1] = 0.25 * (
28        u[2:, 1:-1] + u[:-2, 1:-1] + u[1:-1, 2:] + u[1:-1, :-2]
29    )
30    # Largest pointwise change this sweep — our convergence measure.
31    diff = (u_new - u).abs().max().item()
32    u = u_new
33    if diff < TOL:
34        break
35
36print(f"converged in {step + 1} sweeps")
37print(f"max change in last sweep = {diff:.2e}")
38print(f"u at the centre          = {u[N // 2, N // 2].item():.4f}")
39print(f"u just below top edge    = {u[-2, N // 2].item():.4f}")
40print(f"u just above bottom edge = {u[1,  N // 2].item():.4f}")

Expected output

converged in 2316 sweeps
max change in last sweep = 9.99e-06
u at the centre          = 25.0000
u just below top edge    = 75.7041
u just above bottom edge =  4.2959

Three numbers, three checks

The centre is exactly 25.0 — that is the average of the four boundary edges (100 + 0 + 0 + 0)/4, which the mean value property guarantees for the symmetric centre point.

The values just below the hot top edge and just above the cold bottom edge sum to almost exactly 80.0. That is not an accident: the harmonic field is anti-symmetric around the horizontal midline shifted by 50. The sum of mirror-image values is constant across the field.

Jacobi took 2316 sweeps. A multigrid solver would have hit the same tolerance in roughly $\log_2(N) \approx 6$ sweeps. The slow convergence of Jacobi is the most famous "teach first, replace later" algorithm in numerical analysis.

Why the Boundary Determines Everything

We saw above (as a corollary of the maximum principle) that the boundary data uniquely determines a harmonic function. Let's unpack what this means at the level of intuition.

Specify any continuous function $g$ on the boundary of a nice domain (say a disk, a square, or any region without nasty corners). Then there is exactly one harmonic function in the interior that agrees with $g$ at the boundary. None of the values in the interior have to be specified — they are forced by the boundary data and the equation $\Delta u = 0$ .

The boundary contains everything

The information needed to reconstruct the entire 2D field lives on a 1D curve. The interior is a deterministic readout of the boundary. This is the precise statement that the boundary causes the interior in a stationary problem.

Compare with the heat equation, where you also need an initial condition in addition to the boundary. The heat equation has time — and time needs an origin. Laplace's equation has no time, so it has no initial condition. The boundary alone fixes everything.

Connection to electrostatics

In electrostatics, the potential $\phi$ in a charge-free region satisfies $\Delta \phi = 0$ . The potential on the surfaces of all conductors is the boundary data. Once you know the surface potentials, the field everywhere in between is fixed. That is why an electrostatics problem reduces to specifying conductor voltages — the rest is mathematics.

Common Pitfalls

Forgetting the boundary condition. $\Delta u = 0$ alone has infinitely many solutions (every linear function, every harmonic polynomial). The equation is well-posed only together with boundary data. A Laplace problem without boundary conditions is not a problem yet.
Confusing $\Delta$ with $\nabla$ . $\nabla u$ is a vector — the gradient. $\Delta u = \nabla \cdot \nabla u$ is a scalar — the divergence of the gradient. The Laplacian eats a scalar field and returns a scalar field; the gradient eats a scalar field and returns a vector field.
Updating in place with the wrong loop order. The Jacobi method requires reading ONLY old neighbours, so it needs a buffer. If you write into the same array you read from, you are doing Gauss-Seidel — a different algorithm, also convergent, but with different convergence rates and different error patterns. Both are valid, but they are not the same method.
Believing the discrete and continuum Laplace equations are identical. The discrete equation $u_{i,j} = \tfrac{1}{4} \sum u_{\text{neighbours}}$ is an $O(h^2)$ approximation of $\Delta u = 0$ . The discrete solution on an $N \times N$ grid differs from the true continuum solution by $O(1/N^2)$ at every point. To halve the error you must double the grid resolution — and that quadruples the storage and roughly quadruples the per-sweep cost.
Treating the iteration count as physical time. A Jacobi sweep is not a time step. There is no time in Laplace's equation. The relaxation is purely algorithmic — its only purpose is to drive an inconsistent guess toward a consistent steady state. You can choose different relaxation methods (Jacobi, Gauss-Seidel, SOR, multigrid) without changing the equation being solved.

Summary

Laplace's equation $\Delta u = 0$ is the equation of equilibrium — the steady-state limit of the heat equation, the wave equation, every diffusive or oscillatory process whose right side is the Laplacian. Its solutions are the harmonic functions.

The five facts to remember:

Defining equation: $\Delta u = u_{xx} + u_{yy} = 0$ on the interior, plus prescribed values $u = g$ on the boundary.
Mean value property: a harmonic function equals its average over every circle around every interior point. Equivalent to the equation itself.
Maximum principle: no interior peaks or valleys — extrema live on the boundary.
Uniqueness: the boundary data alone fixes the entire interior field. No initial condition, no time.
Discrete picture: on a grid, the equation becomes "every interior cell is the average of its four neighbours," and Jacobi iteration solves it by repeated averaging.

In the next section we will turn this discussion into a complete framework for boundary value problems: the kinds of boundary conditions you can pose (Dirichlet, Neumann, Robin), what each means physically, and when each yields a unique solution.