Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Recognize Laplace's equation $\nabla^2 u = 0$ as the law of perfect balance — every interior value is the local average of its neighbours.
Rewrite the equation in polar coordinates and explain why the geometry of a disk demands it.
Derive the separated solutions $r^n\cos n\theta$ and $r^n\sin n\theta$ and explain why the companion solutions $r^{-n}$ must be discarded inside the disk.
Match any boundary function $f(\theta)$ to a Fourier series and assemble the unique harmonic extension.
Apply the Poisson integral formula and the mean-value theorem to read off interior values without solving any equation.
Verify the theory numerically — first with plain Python, then with a vectorized PyTorch solver that handles the whole disk in one matrix multiply.

The Big Picture

"Laplace's equation is the mathematics of a system that has run out of reasons to change."

Imagine a thin circular metal plate. You clamp a soldering iron and a block of ice around its rim in an intricate pattern of hot and cold spots — the boundary temperature $f(\theta)$ . You wait. Heat flows around the disk, settling. After enough time, the interior temperature stops changing entirely. What pattern does it settle into?

The answer obeys a single, supremely elegant equation:

\displaystyle \nabla^2 u \;=\; u_{xx} + u_{yy} \;=\; 0.

This is Laplace's equation. Its solutions are called harmonic functions. In any region of the plane where it holds, the function has a stunning property:

The mean-value property — the soul of harmonic functions

The value of $u$ at any point equals the average of $u$ on any circle drawn around that point (so long as the disk inside is also inside the domain). No peaks, no valleys, no surprises — every interior value is dictated by the boundary alone.

For a circular domain, this property is not just an abstract theorem — it is the recipe. In this section we will turn it into an explicit formula:

\displaystyle u(r,\theta) \;=\; \frac{1-r^2}{2\pi}\int_{0}^{2\pi} \frac{f(\varphi)\,d\varphi}{1 - 2r\cos(\theta - \varphi) + r^2}.

Every symbol in this Poisson integral formula will be derived geometrically, interpreted intuitively, and verified by hand. By the end of the section it should feel inevitable.

Why Polar Coordinates Are the Right Language

Cartesian coordinates are flat and unbiased — they treat every direction equally. That is exactly the wrong feature here. The domain is a disk, and a disk has a centre and a boundary; it cares deeply about radius and angle. So we rewrite $u(x,y)$ as $u(r,\theta)$ where

x = r\cos\theta, \qquad y = r\sin\theta.

The boundary becomes the simple set $r = 1$ (taking the disk radius to be 1 without loss of generality). The boundary data becomes a function of one variable, $f(\theta)$ . And — crucially — the differential operator $\nabla^2$ takes a form whose pieces separate cleanly into radial and angular halves.

A small but powerful change of language

The same equation that looks like a coupled second-order PDE in $x, y$ almost falls apart into two independent ODEs in $r, \theta$ . The geometry of the disk is doing the work.

Laplace's Equation in Polar Form

A short calculation with the chain rule (you can do it in two pages, or look it up) converts $u_{xx} + u_{yy} = 0$ into:

\displaystyle u_{rr} \;+\; \frac{1}{r}\, u_{r} \;+\; \frac{1}{r^2}\, u_{\theta\theta} \;=\; 0.

Three terms, each with a clean interpretation:

$u_{rr}$ — pure curvature in the radial direction. Tells you how the temperature bows as you walk straight out from the centre.
$\tfrac{1}{r} u_{r}$ — a geometric correction: as you move out, you sweep through more circumference, so a uniform radial flux must be balanced by a slope. The $1/r$ factor encodes the cylindrical area element.
$\tfrac{1}{r^2} u_{\theta\theta}$ — angular curvature. The factor $1/r^2$ is there because arc length along a circle of radius $r$ is $r\,d\theta$ , so a second derivative with respect to $\theta$ must be rescaled by $r^{-2}$ to be a true second derivative along the surface.

Multiply through by $r^2$ to clear the fractions:

\displaystyle r^2\, u_{rr} \;+\; r\, u_{r} \;+\; u_{\theta\theta} \;=\; 0.

This form is the launching pad for separation of variables — the two halves now have matching degrees in $r$ .

Separation of Variables on the Disk

We try a solution that factors into a radial piece and an angular piece:

u(r,\theta) \;=\; R(r) \cdot \Theta(\theta).

Substitute into the polar Laplace equation and divide by $R\Theta$ :

\displaystyle \frac{r^2 R'' + r R'}{R} \;=\; -\,\frac{\Theta''}{\Theta}.

The left side depends only on $r$ ; the right side depends only on $\theta$ . The only way two such functions can be equal is if both are constant. Call the constant $n^2$ . Then we get two ordinary differential equations:

\displaystyle \Theta'' + n^2 \Theta = 0, \qquad r^2 R'' + r R' - n^2 R = 0.

Angular equation

$\Theta'' + n^2 \Theta = 0$ has the familiar oscillatory solutions $\Theta(\theta) = \cos n\theta$ and $\Theta(\theta) = \sin n\theta$ . The domain wraps around — $\theta$ and $\theta + 2\pi$ are the same physical point — so periodicity forces $n$ to be a non-negative integer.

Why integer n?

If $n$ were anything else (say $n = 1.7$ ), $\cos(1.7\,\theta)$ would not return to its original value as $\theta$ went around the circle. The function would be multi-valued — physically nonsense for a temperature on a disk.

Radial equation — Cauchy–Euler form

$r^2 R'' + r R' - n^2 R = 0$ is a Cauchy–Euler equation. The natural ansatz is $R(r) = r^\alpha$ . Substituting:

\alpha(\alpha-1) r^\alpha + \alpha r^\alpha - n^2 r^\alpha = 0 \;\Longrightarrow\; \alpha^2 = n^2 \;\Longrightarrow\; \alpha = \pm n.

So the general radial solution for $n \ge 1$ is

R(r) = A_n r^n + B_n r^{-n}.

The special case $n = 0$ gives a double root; its solutions are

R_0(r) = A_0 + B_0 \ln r.

Boundedness at the Origin Removes Half the Solutions

The disk contains the point $r = 0$ . The functions $r^{-n}$ for $n \ge 1$ and $\ln r$ for $n = 0$ blow up there. Real physical temperature does not blow up — there is no point source at the centre.

So we set every $B_n = 0$ . The surviving radial functions are $R_n(r) = r^n$ for $n = 0, 1, 2, \dots$ .

In an annulus, the story changes

If the domain were a ring $1 \le r \le 2$ instead of a disk, the origin would not be inside, and the $r^{-n}$ and $\ln r$ terms would survive. The general solution would then be a sum of both sets of modes — and you would need two boundary conditions to fix the coefficients.

Combining the radial and angular pieces, the bounded separated solutions on the disk are:

\displaystyle u_n(r,\theta) \;=\; r^n\left(a_n \cos n\theta + b_n \sin n\theta\right), \quad n = 0, 1, 2, \dots

Each of these is a perfectly good harmonic function on the disk. The general solution is their superposition:

\displaystyle u(r,\theta) \;=\; a_0 \;+\; \sum_{n=1}^{\infty} r^n\bigl(a_n \cos n\theta + b_n \sin n\theta\bigr).

Two ideas in one formula

The angular part is a Fourier series — the same tool you used for vibrating strings and heat flow on a rod. The radial part attaches a damping factor $r^n$ in front of every mode. Higher modes decay faster as you move inward, so the centre of the disk sees only the gentlest features of the boundary.

Coupling to the Boundary via Fourier Series

At $r = 1$ the formula must match the prescribed boundary data $f(\theta)$ :

\displaystyle f(\theta) \;=\; a_0 + \sum_{n=1}^{\infty}\bigl(a_n \cos n\theta + b_n \sin n\theta\bigr).

This is just the Fourier series of $f$ . The coefficients are computed by the usual projection integrals:

\displaystyle a_0 = \frac{1}{2\pi}\int_0^{2\pi} f(\theta)\,d\theta, \quad a_n = \frac{1}{\pi}\int_0^{2\pi} f(\theta)\cos n\theta\,d\theta, \quad b_n = \frac{1}{\pi}\int_0^{2\pi} f(\theta)\sin n\theta\,d\theta.

Notice that $a_0$ is the average of $f$ on the boundary — exactly the boundary mean.

Once those numbers are computed, the interior solution is fixed:

\displaystyle u(r,\theta) \;=\; a_0 + \sum_{n=1}^{\infty} r^n\bigl(a_n \cos n\theta + b_n \sin n\theta\bigr).

Mode n	Angular shape	Damping factor	Reaches the centre?
0	constant	r⁰ = 1	Yes — survives at r = 0
1	cos θ, sin θ	r	Falls off linearly
2	cos 2θ, sin 2θ	r²	Falls off quadratically
n	cos nθ, sin nθ	rⁿ	Decays like rⁿ

The intuition baked into the damping

A sharp feature on the boundary (a hot spike) has lots of high-frequency content — large $n$ components. Those modes are crushed by $r^n$ as you move inward. By the time you reach the centre only the boundary average survives. Diffusion in space, baked into a single algebraic factor.

Interactive: Basis Modes of the Disk

Here is the basis you just derived. Pick a mode $n$ and the cosine / sine mix; watch the radial damping $r^n$ in action.

Two things to notice:

The pattern always vanishes at the centre except for $n = 0$ .
There are exactly $2n$ nodal radii — rays where the mode is zero. Between any two consecutive rays the sign flips. This is the same combinatorics as the zeros of $\cos n\theta$ on the unit circle.

The Poisson Integral Formula

The Fourier-series solution is exact, but it asks us to compute infinitely many coefficients. Poisson's trick is to sum the series in closed form, swapping the infinite sum for a single integral against a magical kernel.

Substitute the projection formulas for $a_n, b_n$ into the series and swap the sum with the integral:

\displaystyle u(r,\theta) = \frac{1}{2\pi}\int_0^{2\pi} f(\varphi)\left[1 + 2\sum_{n=1}^{\infty} r^n \cos n(\theta - \varphi)\right] d\varphi.

The bracket is a geometric series in disguise. Using $\cos n\alpha = \tfrac{1}{2}(e^{in\alpha} + e^{-in\alpha})$ and summing the resulting geometric series gives

\displaystyle 1 + 2\sum_{n=1}^{\infty} r^n \cos n\alpha \;=\; \frac{1-r^2}{1 - 2r\cos\alpha + r^2}.

Putting it back into the integral with $\alpha = \theta - \varphi$ produces the Poisson integral formula:

\displaystyle u(r,\theta) \;=\; \frac{1-r^2}{2\pi}\int_0^{2\pi} \frac{f(\varphi)\,d\varphi}{1 - 2r\cos(\theta - \varphi) + r^2}.

The quantity in the integrand is the Poisson kernel:

\displaystyle P_r(\alpha) \;=\; \frac{1-r^2}{1 - 2r\cos\alpha + r^2}.

Geometric reading: the denominator $1 - 2r\cos\alpha + r^2$ is, by the law of cosines, the squared distance from the interior point at radius $r$ to the boundary point at angular offset $\alpha$ . So

\displaystyle P_r(\alpha) \;=\; \frac{1 - r^2}{|\text{interior point} - \text{boundary point}|^2}.

Boundary points near the interior probe (small distance) get a large weight; boundary points on the opposite side (large distance) get a small weight. The Poisson kernel is an inverse-square influence function — exactly what physical intuition (think electrostatics, think gravity) would suggest.

A different way to read the formula

$u(r,\theta)$ is a weighted average of the boundary data. The weights are positive, sum to 1 over the circle, and concentrate near the boundary point closest to the probe. The interior temperature is just "what you see when you smear the boundary by these weights."

Interactive: The Poisson Kernel

Slide $r$ from 0 toward 1 and watch the kernel transform.

Three limits worth noting:

At $r = 0$ , the kernel is the constant $1$ . The integral becomes a plain average of $f$ over the circle — this is the mean-value theorem.
At $r \to 1$ , the kernel becomes a spike at $\alpha = 0$ . The integral becomes $f(\theta)$ itself — the interior solution recovers the boundary as you approach it.
The peak height is $P_r(0) = (1+r)/(1-r)$ . The kernel remains positive for every $r < 1$ and $\alpha \in [-\pi, \pi]$ — this positivity is the deep reason harmonic functions obey a maximum principle (no interior peaks).

The Mean-Value Theorem

Set $r = 0$ in either the Fourier sum or the Poisson integral. Every $r^n$ with $n \ge 1$ evaluates to zero, and the Poisson kernel collapses to $P_0(\alpha) = 1$ . Both formulas say the same thing:

\displaystyle u(0,\theta) \;=\; \frac{1}{2\pi}\int_0^{2\pi} f(\varphi)\,d\varphi \;=\; a_0.

The value of a harmonic function at the centre of a disk equals its average on the boundary. And by translation invariance: the value at any point equals its average on any circle around that point, as long as the surrounding disk lies inside the domain.

A startling consequence: no interior extrema

If $u$ had a local maximum at some interior point, the value at that point would be strictly greater than the average on a small surrounding circle — contradicting the mean-value property. The same argument rules out interior minima. So every local extremum of a harmonic function lies on the boundary. This is the maximum principle, and it is one of the cornerstones of elliptic PDE theory.

Interactive: The Full Disk Explorer

Time to play. Pick a boundary preset; the entire interior heatmap is the unique harmonic extension. Drag the probe to read off $u(r,\theta)$ ; turn on the Poisson kernel halo to see how a single boundary point influences each interior pixel. Confirm the mean-value theorem by checking that the centre value always equals the boundary average.

Two experiments worth trying:

Pick cos(φ) as the boundary. The interior solution is exactly $u = r\cos\theta = x$ — a perfectly linear function in $x$ . The mean is zero, so the centre is zero. Verify both.
Pick cos(5φ) and watch the colour fade to a near-uniform grey by $r = 0.5$ . The damping factor $r^5$ at $r = 0.5$ is just $0.03125$ — high-frequency boundary features are essentially invisible in the middle.

Worked Example (Collapsible, By Hand)

Let us solve one problem completely by hand. Boundary: $f(\theta) = 50 + 50\cos\theta$ . Find $u(r,\theta)$ on the unit disk and evaluate at three interior points.

Show full hand-worked solution

Step 1. Read off the Fourier coefficients

The boundary $f(\theta) = 50 + 50\cos\theta$ already is its own Fourier series — no integrals needed. By inspection,

$a_0 = 50$
$a_1 = 50$
All other $a_n, b_n = 0$ .

Step 2. Attach the radial damping

The mode-1 term gets a factor of $r^1 = r$ . So

u(r,\theta) \;=\; 50 \;+\; 50\,r\,\cos\theta.

In Cartesian coordinates this is just $u = 50 + 50 x$ — a linear function, which is automatically harmonic. Notice how the polar machinery reproduces the obvious Cartesian answer.

Step 3. Evaluate at three interior points

Point	(r, θ)	x = r cos θ	u = 50 + 50 r cos θ
Centre	(0, 0)	0	50.00
Right midpoint	(0.5, 0)	+0.50	75.00
Left midpoint	(0.5, π)	−0.50	25.00
Top midpoint	(0.5, π/2)	0	50.00
NE diagonal	(0.7, π/4)	≈ +0.495	≈ 74.75

Step 4. Verify the mean-value theorem

The boundary average is

\displaystyle \frac{1}{2\pi}\int_0^{2\pi}(50 + 50\cos\theta)\,d\theta \;=\; 50 \;+\; 50 \cdot 0 \;=\; 50.

And indeed $u(0,0) = 50$ . Done.

Step 5. Sanity-check against the Poisson integral

Take any interior point, say $(r,\theta) = (0.5, 0)$ . The Poisson integral becomes

\displaystyle u(0.5, 0) = \frac{1 - 0.25}{2\pi}\int_0^{2\pi} \frac{50 + 50\cos\varphi}{1 - \cos\varphi + 0.25}\,d\varphi.

Direct numerical evaluation (or a contour integral if you are feeling ambitious) gives exactly $75$ — matching the much simpler Fourier-series answer. The two formulas must agree; they are two different windows onto the same harmonic extension.

Step 6. Try a second case to keep your hand in

Repeat with $f(\theta) = \cos(3\theta)$ . The answer is $u(r,\theta) = r^3 \cos(3\theta)$ . Three sets of two nodal radii each, vanishing at the centre, inheriting a single Fourier mode from the boundary.

Python: Fourier Coefficients and the Truncated Series

Plain Python first — no fancy libraries, no shortcuts. We will compute the Fourier coefficients of a square-wave boundary, attach the radial damping $r^n$ to each mode, and read off the interior value at a probe point. Every line is annotated.

Truncated harmonic extension of a square-wave boundary

🐍python

Explanation(18)

Code(35)

1NumPy import

We need np.cos, np.sin, np.where, and array arithmetic to evaluate the boundary samples and assemble the Fourier sum. PyTorch comes later; we keep this first pass in plain NumPy so the algorithm is bare and easy to hand-check.

4Boundary function f(φ)

The boundary data is a square wave: 100 on the top half of the circle (where sin φ > 0) and 0 on the bottom half. This is a deliberately rough boundary — it has a jump discontinuity at φ = 0 and φ = π. Laplace's equation will smooth this jump out as we move inward; that smoothing is the entire point of harmonic functions.

10M = 720 boundary samples

We sample the boundary at 720 equally-spaced angles. 720 is overkill for a smooth f but it gives an essentially exact estimate of the Fourier coefficients even for our discontinuous square wave. The trapezoidal rule on a 2π-periodic function is spectrally accurate, so even very large M is cheap.

11Sample angles φ in [0, 2π)

endpoint=False makes the last sample land at 2π·(M−1)/M, NOT at 2π. We do this because f is 2π-periodic — including 2π would double-count the φ = 0 sample and bias the integrals.

12fb = f(φ)

The boundary values themselves, as a length-M array. This is the raw boundary data we will project onto sines and cosines.

15Fourier truncation order N

We will keep modes n = 0, 1, ..., 40. The square wave has spectrum that falls off only like 1/n — so 40 modes is plenty for the centre but the boundary itself will still show Gibbs overshoot. Try lowering N to 5 and see the boundary become a smooth blob.

18a₀ = boundary mean

The zeroth Fourier coefficient is just the average of f over the circle. For our square wave it equals (100 + 0) / 2 = 50. By the mean-value theorem, this is also u at the centre — a fact we will verify two lines below.

19Loop over the modes

For every n ≥ 1 we project the boundary onto cos(nφ) and sin(nφ). aₙ measures how much of f looks like cos(nφ) on the circle; bₙ measures how much looks like sin(nφ). Together (aₙ, bₙ) IS the discrete Fourier transform of fb.

20aₙ projection

The continuous formula is aₙ = (1/π) ∫₀^{2π} f(φ) cos(nφ) dφ. Replacing the integral with a Riemann sum on M equally-spaced points and using dφ = 2π/M gives aₙ ≈ (2/M) · Σ fb·cos(nφ). The constant 2/M absorbs both the 1/π and the dφ.

21bₙ projection

Same formula with sin(nφ). For our square wave the symmetry f(π − φ) = f(φ) does NOT hold, but f(−φ) = −f(φ) + 100 holds (an odd-around-zero offset), so the cosine coefficients are mostly zero except a₀ and the sine coefficients carry most of the spectrum.

24Pick an interior probe point

r = 0.5 is halfway out to the boundary; θ = π/2 points to the top of the disk. Intuitively the top is in the hot region, so we expect u > 50 there. The mean-value theorem alone tells us u(0) = 50, but moving off-centre toward the hot half MUST push the value up.

25Start with the constant term

Every harmonic function inside the disk inherits the boundary average as its DC offset. Without this term the entire solution would be wrong by a constant.

26Mode-by-mode summation

The harmonic extension is u(r, θ) = a₀ + Σₙ rⁿ · (aₙ cos(nθ) + bₙ sin(nθ)). The factor rⁿ is the key: it DAMPS each mode the further inside the disk we go. High-frequency boundary features (large n) decay quickly toward the centre — they cannot reach the middle of the disk.

27Single line that assembles every term

The same separation-of-variables solution we derived by hand, written out in numpy. (r**n) is the radial factor; a[n]·cos(n·θ) + b[n]·sin(n·θ) is the angular factor. Together they form the n-th separated solution, weighted by the projection of the boundary onto that mode.

30u(0) = a₀ (the mean-value theorem)

At r = 0 every rⁿ term with n ≥ 1 evaluates to 0. Only a₀ survives — and a₀ is the boundary average. This is the mean-value theorem for harmonic functions, popping straight out of the Fourier sum. No integration required.

32Print boundary mean

Expected output: 50.00000 — exactly the average of 100 (top half) and 0 (bottom half).

33Print u at (0.5, π/2)

Expected output: about 79.516 — well above the mean. The truncated sum at N = 40 is already within 10⁻⁴ of the limit value 79.5167; truncation at lower N would still hit the same neighbourhood because (0.5)ⁿ kills the tail very quickly.

34Print u at the centre

Expected output: 50.00000 — identical to a₀. This is the mean-value theorem expressed as one print statement.

17 lines without explanation

1import numpy as np
2
3# Boundary data on the unit circle: a square wave.
4# f(phi) = 100  on the top half  (sin(phi) > 0)
5# f(phi) =   0  on the bottom half
6def f(phi):
7    return np.where(np.sin(phi) > 0, 100.0, 0.0)
8
9# Sample the boundary uniformly in phi.
10M  = 720
11phi = np.linspace(0.0, 2*np.pi, M, endpoint=False)
12fb  = f(phi)
13
14# Compute the Fourier coefficients a_n, b_n  (n = 0, 1, ..., N).
15N  = 40
16a  = np.zeros(N + 1)
17b  = np.zeros(N + 1)
18a[0] = fb.mean()                                # the boundary average
19for n in range(1, N + 1):
20    a[n] = (2.0 / M) * np.sum(fb * np.cos(n * phi))
21    b[n] = (2.0 / M) * np.sum(fb * np.sin(n * phi))
22
23# Evaluate the harmonic extension u(r, theta) at one interior point.
24r, theta = 0.5, np.pi / 2
25u = a[0]
26for n in range(1, N + 1):
27    u += (r ** n) * (a[n] * np.cos(n * theta) + b[n] * np.sin(n * theta))
28
29# Also evaluate at the centre — by the mean-value theorem this must
30# equal the boundary average a[0].
31u_centre = a[0]
32
33print(f"a_0 (boundary mean)       = {a[0]:.5f}")
34print(f"u(r=0.5, theta=pi/2)      = {u:.5f}    (N = {N})")
35print(f"u(r=0, theta=*) = a_0     = {u_centre:.5f}")

Run it. The boundary mean is $50$ , the centre value is $50$ , and the probe value at $(0.5, \pi/2)$ is about $79.5167$ . Two numbers worth burning into memory — they will come back from a completely different formula in the next subsection.

Python: The Poisson Integral Directly

Now the rival method: skip the Fourier series and evaluate the Poisson integral by Riemann sums. One function, four lines of active computation.

Poisson integral as a single quadrature

🐍python

Explanation(13)

Code(25)

1Import

Same dependencies as before. Plain NumPy is enough; the Poisson formula is just a single one-dimensional integral.

4Boundary function

We keep the same square-wave boundary so the two approaches — Fourier-series sum versus Poisson integral — can be compared on the same numbers.

8M = 4000 quadrature points

The Poisson kernel has a sharp peak near α = 0 when r is close to 1, so we need a finer grid than the Fourier method needed. M = 4000 is enough that even at r = 0.99 the quadrature error is well below 10⁻⁵.

9Boundary angles

Same uniform grid on [0, 2π). The integration measure dφ = 2π / M is constant, which lets us use plain numpy sums instead of a trapezoidal-rule helper.

10fb = boundary samples

Vectorized — fb[k] is the boundary temperature at angle phi[k]. We do not need to recompute it inside the integral loop.

11dphi

The infinitesimal in the Riemann sum approximation of the integral. Pulling it out as a constant means the integrand inside the sum is just K · fb, multiplied by dphi at the end.

14Poisson kernel P_r(α)

The single most important function in this section. P_r(α) = (1 − r²)/(1 − 2r cos α + r²). Geometric interpretation: 1 − 2r cos α + r² is the squared distance from the interior point at radius r to the boundary point at angle α (relative to θ), by the law of cosines. So the kernel is large when α is small (close boundary point, high influence) and small when α is large (far boundary point, low influence).

18Vectorized kernel evaluation

θ − phi is a length-M array of angular offsets. poisson_kernel returns a length-M array of kernel values, one per boundary sample. The sum K · fb is a weighted average of boundary values, with the Poisson kernel as the weight function.

19Riemann-sum approximation of the integral

The continuous formula is u(r, θ) = (1/2π) ∫₀^{2π} P_r(θ − φ) · f(φ) dφ. Discretizing with M equally-spaced points and dφ = 2π/M turns it into (1/2π) · Σ K · fb · dφ. This is the workhorse equation of the entire section.

21Sample at the top of the disk

Expected 79.516713 — within 4 ppm of the Fourier-series answer. Two completely different formulas, same number.

22Sample at the bottom of the disk

Expected 20.483287. Notice this is exactly 100 − 79.516713 — a manifestation of the symmetry f(φ + π) = 100 − f(φ) for our boundary, which makes u(r, θ + π) = 100 − u(r, θ).

23Sample near the hot boundary

Expected 96.652458 — almost 100. Moving from r = 0.5 to r = 0.9 toward the top boundary, the Poisson kernel becomes a sharp spike at the source angle θ = π/2 directly under the probe, so u closely tracks the local boundary value.

24Centre value (mean-value theorem)

At r = 0 the Poisson kernel reduces to P_0(α) = (1 − 0)/(1 − 0 + 0) = 1. The integral is then simply the average of f. Expected output: 50.000000 — the mean of the boundary, identical to a₀.

12 lines without explanation

1import numpy as np
2
3# Same boundary data as before.
4def f(phi):
5    return np.where(np.sin(phi) > 0, 100.0, 0.0)
6
7# A high-resolution boundary grid for the integral.
8M   = 4000
9phi = np.linspace(0.0, 2*np.pi, M, endpoint=False)
10fb  = f(phi)
11dphi = 2*np.pi / M
12
13# Poisson kernel evaluated at a single interior point (r, theta).
14def poisson_kernel(r, alpha):
15    return (1 - r**2) / (1 - 2*r*np.cos(alpha) + r**2)
16
17# Compute u(r, theta) by direct quadrature of the Poisson integral.
18def u(r, theta):
19    K = poisson_kernel(r, theta - phi)
20    return (1.0 / (2*np.pi)) * np.sum(K * fb) * dphi
21
22print(f"u(0.5, pi/2)   = {u(0.5,  np.pi/2):.6f}")
23print(f"u(0.5, -pi/2)  = {u(0.5, -np.pi/2):.6f}")
24print(f"u(0.9, pi/2)   = {u(0.9,  np.pi/2):.6f}")
25print(f"u(0, anything) = {u(0.0, 1.234):.6f}    # mean value theorem")

The output matches the Fourier solver to four decimal places at every probe point — and uses exactly zero Fourier coefficients to do it. Two windows, one harmonic function.

PyTorch: A Vectorized Disk Solver

The Poisson integral evaluated at a single point is a one-dimensional integral. Evaluated on a $G \times G$ grid of interior points, it becomes a single matrix–vector product:

\displaystyle u_i \;\approx\; \frac{1}{2\pi}\sum_{k=0}^{M-1} P_{r_i}(\theta_i - \varphi_k)\, f(\varphi_k)\, \Delta\varphi.

With $P$ a $(G^2) \times M$ matrix and $f$ a length- $M$ vector, this is one PyTorch operation. The same code runs on a GPU without a single tweak.

Whole-disk Poisson solver: one matmul

🐍python

Explanation(22)

Code(38)

1PyTorch + math

We switch to PyTorch tensors. The same Poisson integral now becomes a single matrix–vector product, which means we get the disk solver for free on a GPU. The math module provides math.pi as a Python scalar (PyTorch does not export a tensor-friendly pi by default).

5Boundary samples on the circle

M = 720 angle samples. We slice off the last element of linspace(0, 2π, M+1) because the function is periodic — we do not want phi to include both 0 and 2π, which would double-count the same physical boundary point.

6fb — the boundary temperature tensor

torch.where is the tensor analogue of np.where. For each angle phi[k] we pick 100.0 if sin(phi[k]) > 0 (top half) and 0.0 otherwise. The wrapping of the scalars in torch.tensor is needed because torch.where requires its 'true' and 'false' branches to be tensors of compatible type.

11Interior Cartesian grid

We use a G×G Cartesian grid covering the bounding square [−1, 1]². A polar grid would also work but Cartesian makes plotting trivial (one imshow), and PyTorch's atan2 and sqrt are vectorized so the conversion is cheap.

12Stack ys

Same construction in y. xs and ys are 1-D; we will broadcast them into a 2-D grid with meshgrid in the next line.

13Build the 2D grid

indexing="xy" makes X[i, j] = xs[j] and Y[i, j] = ys[i] — the convention every plotting library expects. X and Y are both (G, G).

14Radial coordinate R

R[i, j] = √(X² + Y²). Vectorized over the whole grid — no Python loop. We will use R both as the harmonic-extension radius and as a mask predicate (R < 1 = inside the disk).

15Angular coordinate Θ

atan2(Y, X) returns the polar angle in [−π, π]. It handles the four quadrants correctly — atan(Y/X) alone cannot do that. Together (R, Θ) is the polar coordinate of every Cartesian grid point.

16Disk mask

Boolean tensor of shape (G, G), True for grid cells inside the unit disk. We will use it later to zero out the kernel for points outside the disk so the solution does not bleed into the bounding square.

20Flatten Θ for broadcasting

We want to broadcast the (G*G,) angle vector against the (M,) boundary-phi vector to form an (G*G, M) matrix of differences. reshape(-1, 1) gives Θ a trailing length-1 axis so broadcasting works.

21Flatten R for broadcasting

Same trick for the radius. Now both R_flat and Theta_flat have shape (G*G, 1), ready to broadcast against the (M,) boundary angles.

22Pairwise angle differences alpha

phi.unsqueeze(0) reshapes the boundary angles to (1, M). Subtracting from Theta_flat broadcasts to (G*G, M). alpha[i, k] = θᵢ − φₖ, the angular offset between every interior point and every boundary point — exactly what the Poisson kernel eats.

23Build the kernel matrix P

Same formula as in the Python version, but now applied to the entire (G*G, M) tensor at once. Every row of P is the kernel weight from one interior point to all M boundary samples. Conceptually this matrix IS the Poisson integral operator.

24Numerical clamp

Near r = 1 and α near 0 the kernel can become a huge spike; tiny floating-point noise can occasionally make the denominator negative by a few ULPs. clamp(min=0) suppresses those nonphysical specks without affecting correct values.

25Mask outside the disk

For grid cells with R ≥ 1 we set the kernel weights to zero. The Poisson formula is only defined inside the disk; we want the solution outside to be visibly blank, not a polynomial extrapolation of the kernel.

28Discrete integration step dphi

Same constant as before — 2π / M. We multiply once at the end rather than inside the kernel to keep the matmul as clean as possible.

29One matrix multiply solves the entire disk

P @ fb has shape (G*G,). Element i of the result is Σₖ P[i, k] · fb[k] — exactly the Riemann sum of the Poisson integral at interior point i. Multiplying by dphi / (2π) finishes the normalization. On a GPU this single line replaces a doubly-nested Python loop and runs in microseconds.

30Reshape back to the 2D grid

We flatten before the matmul to keep the algebra clean, then put the result back into (G, G) so it can be plotted with imshow. u[i, j] is now the solution at Cartesian grid cell (xs[j], ys[i]).

33Centre index

G is even, so the 'centre' is between two grid cells. G // 2 is the nearest cell to (0, 0); for G = 256 the cell is at (xs[128], ys[128]) ≈ (0.0039, 0.0039) — essentially the origin.

34Centre value (mean-value theorem)

Expected ≈ 50.0 — the boundary average. Identical (to plotting precision) to the a₀ value from the Fourier solver. The same theorem keeps reappearing because it is a corollary of P_0(α) = 1.

35Top probe value

Row G/4 in our 'xy' convention corresponds to y ≈ +0.5 — the top half. Expected ≈ 79.51 — matches the NumPy answer to four decimals.

36Bottom probe value

Row 3G/4 corresponds to y ≈ −0.5 — the bottom half. Expected ≈ 20.48, the antipode of 79.51 under the symmetry f(φ + π) = 100 − f(φ).

16 lines without explanation

1import torch
2import math
3
4# 1.  Boundary data on the unit circle.
5M   = 720
6phi = torch.linspace(0.0, 2*math.pi, M + 1)[:-1]   # length M, no double 2pi
7fb  = torch.where(torch.sin(phi) > 0,
8                  torch.tensor(100.0),
9                  torch.tensor(0.0))               # square wave
10
11# 2.  Interior grid (Cartesian) — every point (x, y) inside the unit disk.
12G  = 256
13xs = torch.linspace(-1.0, 1.0, G)
14ys = torch.linspace(-1.0, 1.0, G)
15X, Y = torch.meshgrid(xs, ys, indexing="xy")
16R     = torch.sqrt(X**2 + Y**2)
17Theta = torch.atan2(Y, X)
18mask  = R < 1.0                                     # inside the disk
19
20# 3.  Poisson kernel as a (G*G, M) matrix.
21#     P[i, k]  =  P_{r_i}( theta_i  -  phi_k )
22Theta_flat = Theta.reshape(-1, 1)                   # (G*G, 1)
23R_flat     = R.reshape(-1, 1)
24alpha      = Theta_flat - phi.unsqueeze(0)          # (G*G, M)
25P = (1 - R_flat**2) / (1 - 2*R_flat*torch.cos(alpha) + R_flat**2)
26P = P.clamp(min=0.0)                                # ignore numerical specks
27P[~mask.reshape(-1), :] = 0.0                       # zero outside disk
28
29# 4.  ONE matrix–vector product evaluates u on the entire grid.
30dphi = 2*math.pi / M
31u    = (P @ fb) * dphi / (2*math.pi)                # (G*G,)
32u    = u.reshape(G, G)
33
34# 5.  Read off two diagnostic values.
35i_centre = G // 2
36print(f"u at centre              = {u[i_centre, i_centre].item():.5f}")
37print(f"u at top, r ~ 0.5        = {u[G//4, i_centre].item():.5f}")
38print(f"u at bottom, r ~ 0.5     = {u[3*G//4, i_centre].item():.5f}")

The centre lands at $50$ , the top probe at about $79.5$ , the bottom probe at about $20.5$ — the same answers as both Python methods. The mean-value theorem keeps re-emerging because it is a property of the function, not of the method.

Common Pitfalls

Forgetting boundedness. Including $r^{-n}$ or $\ln r$ terms in a disk solution leads to a function that blows up at the origin — physically wrong. These terms only belong on an annulus.
Confusing $r$ with $r/R$ . On a disk of radius $R \ne 1$ the radial factor must be $(r/R)^n$ , not $r^n$ . Setting $R = 1$ hides the rescaling — always check the dimensions.
Mistaking truncation error for a real feature. A truncated Fourier sum exhibits Gibbs overshoot near boundary discontinuities. That ~9% overshoot persists no matter how many terms you keep — but vanishes if you use the Poisson integral instead. Discontinuous boundaries are best handled with the integral form.
Numerical instability near $r \to 1$ . The Poisson kernel denominator can come within ULPs of zero when $\alpha \to 0$ . Increase $M$ or stay slightly away from the boundary; a clamp(min=0) is a cheap safety net.
Forgetting periodicity. $f(\theta)$ must equal $f(\theta + 2\pi)$ . If your boundary data is given as a function on $[0, 2\pi]$ with $f(0) \ne f(2\pi)$ , the resulting Fourier series secretly fits the periodic extension — which has a jump you may not have intended.

Summary

On a disk, Laplace's equation has an exceptionally clean solution structure that we can summarise in one box:

The complete recipe

Rewrite $\nabla^2 u = 0$ in polar form: $r^2 u_{rr} + r u_r + u_{\theta\theta} = 0$ .
Separate variables. The bounded basis solutions are $1, r^n\cos n\theta, r^n \sin n\theta$ .
Expand the boundary data $f(\theta)$ in a Fourier series with coefficients $a_n, b_n$ .
Multiply each mode by $r^n$ to damp it inward; sum:
$\displaystyle u(r,\theta) = a_0 + \sum_{n \ge 1} r^n\bigl(a_n\cos n\theta + b_n\sin n\theta\bigr).$
Equivalently, integrate against the Poisson kernel:
$\displaystyle u(r,\theta) = \frac{1-r^2}{2\pi}\int_0^{2\pi}\frac{f(\varphi)\,d\varphi}{1 - 2r\cos(\theta-\varphi) + r^2}.$
Read off the mean-value theorem and the maximum principle for free.

The disk is a special domain, but the lessons are general. On every two-dimensional region, harmonic functions are boundary-determined averages: every interior value is the weighted balance of its boundary, the high frequencies are suppressed as you move inward, and no interior point can be hotter than the hottest part of the rim. The disk gave us the cleanest version of this story — and a single matrix multiply that delivers the entire steady-state temperature field in one shot.