WI20

Distributed Systems & Fields

PDEs and vector fields in continuous media — waves, heat diffusion, and electromagnetic field structure, solved both analytically and on a grid.

Concepts learned

Tech

Python
NumPy
Matplotlib

Hero demo Charge Field Plotter

N = 2

max ∣ E ∣ = 84.63

max ∣ V ∣ = 8.25

A classic dipole: +1 on the left, -1 on the right. Field sampled on a 20×20 grid over [-2, 2] × [-2, 2].

nx (grid columns)20

ny (grid rows)20

Arrow scale1.0

Show potential (0=off, 1=on)0

charges 2 · max |E| = 84.63

Reflection

Pre-interview preview. Akwasi’s first-person reflection on this course is pending — track at issue #48.

1D wave equation (string)

Plucked-string wave equation with reflective boundary conditions.

Walk me through this step by step

Pluck a guitar string and watch the shape travel, oscillate, return, and oscillate again. One simple rule governs every plucked string, every taut wire, every acoustic-pipe mode. That rule is the wave equation, and it falls out of nothing more exotic than Newton’s second law applied to a tiny segment of the string.
Here is the idea in plain English. How fast a point on the string accelerates up or down depends on how curved the string is right there. If the string is straight at that point, it doesn’t accelerate. If it’s bent upward (like the bottom of a smile), it gets yanked up. If it’s bent downward (like the top of a frown), it gets pushed down. Curvature creates acceleration. That’s it.
Where does the curvature-acceleration rule come from? Think of a tiny segment of string under tension. Tension pulls on both ends. If the string is straight there, the two pulls point in opposite directions and cancel — no net force. If the string is curved, both pulls have a small vertical component pointing toward the “inside” of the curve. Newton’s $F = ma$ F = ma turns that vertical force into vertical acceleration.
In symbols, that becomes
$\frac{\partial ^{2} u}{\partial t ^{2}} = c^{2} \frac{\partial ^{2} u}{\partial x ^{2}}$ \frac{\partial^2 u}{\partial t^2} = c^2 \, \frac{\partial^2 u}{\partial x^2}
where u(x,t) is the string’s height at position x and time t. The left side is “how fast the string is accelerating right now.” The right side is “how curved the string is right now,” times a constant $c^{2}$ c^2. The constant c is the wave speed — bigger for tighter or lighter strings.
The boundary matters. A guitar string is clamped at both ends — $u (0, t) = u (L, t) = 0$ u(0, t) = u(L, t) = 0 for all time. That is what makes it different from a whip or a free-hanging rope. The clamps mean only certain “shapes” can vibrate cleanly — the modes. Every other initial shape is just a sum of these modes, oscillating at different rates.
The allowed shapes are $sin (mπ x / L)$ \sin(m\pi x / L) for m = 1, 2, 3, … The “Fundamental” preset shows mode 1: a single bump from one end to the other. The “Second harmonic” preset shows mode 2: a full sine wave with a node — a permanently still point — in the middle. (Press a string at the 12th fret of a guitar and you have forced mode 2.)
The demo doesn’t use that closed form to draw the animation, though — it solves the equation directly on a grid. Replace the smooth string with N points spaced dx apart. For each interior point, “value next” = 2·(value now) − (value before) + curvature term. The curvature term is just (left neighbour + right neighbour − 2·middle), which is the discrete version of $u_{xx}$ u_{xx}. This is the leapfrog scheme.
There is one inequality the numerical scheme must respect: $c \cdot d t / d x \leq 1$ c \cdot dt / dx \le 1. Physically: in one timestep dt, the true wave moves a distance c·dt. If that distance is more than one grid spacing dx, the wave has “skipped” past data the scheme needed, and the simulation blows up. This is called the CFL condition, and it is the one trap every PDE solver shares.
Try it yourself.

Load the “High resolution (200 pts)” preset. Notice that dt drops to 0.002 (versus 0.005 in the fundamental preset) precisely because dx shrank by the same factor. Doubling resolution roughly doubles the cost in time too. That’s why “just use a finer grid” is usually too easy a wish in simulation work.
The featured problem “Standing-wave frequencies on a fixed-end string by separation of variables” just below derives the mode shapes from scratch by assuming the solution factors as X(x) · T(t). For the travelling-wave (rather than standing-wave) picture, the featured problem “d’Alembert solution to the 1D wave equation” decomposes any initial pulse into a left-mover plus a right-mover. And to see what happens when two waves meet on a 2D membrane, jump to the 2D wave interference demo immediately below.

c = 1.00

r = 0.395

t = 0.00 s

Standing wave mode 1 on a 1 m string with wave speed c = 1.00 m/s, 80 grid points, timestep 0.005 s, CFL ratio r = 0.395 (stable when r ≤ 1).

c (wave speed)1.00 m/s

Speed of transverse waves on the string.

dt (timestep)0.005 s

Smaller dt keeps the CFL ratio r = c·dt/dx ≤ 1 for stability.

nGrid (grid points)80

Number of points discretising the 1 m string.

mode (harmonic)1

Initial standing-wave mode number m: u₀(x) = sin(mπx/L).

CFL r = 0.395 · t = 0.00 s · mode 1

2D wave interference

Two-source interference pattern with adjustable phase, frequency, and source separation.

c = 1.00

r = 0.020

n = 0

max ∣ u ∣ = 0.000

nGrid (square)60

Grid resolution n × n cells.

Wave speed c1.00

Propagation speed. Combined with dt and dx, sets the CFL ratio r.

Timestep dt0.020

Integrator step size. CFL stability requires c·dt/dx ≤ 1/√2.

Source: single0

Forcing pattern: single, two-sources (interference), or corner.

step n = 0 · max |u| = 0.000 · r = 0.020

1D heat equation

Finite-difference solver for heat diffusion along a rod; tune the boundary conditions and watch the profile relax.

Walk me through this step by step

You pull a metal rod out of a fire and grab the cool end. The heat doesn’t leap to your hand — it crawls along the rod, smoothing out the hot spot a little more every second. Eventually the whole rod settles to one warm temperature. The heat equation is the rule that governs that crawl, and it turns up wherever something diffuses: temperature in a wall, ink in still water, even the price of a stock option.
Forget formulas for a moment. The whole physics is this: heat flows from hot to cold, and the flow is faster when the temperature gap is bigger. Two beads sitting at 100° and 0° pull on each other much harder than two beads at 100° and 99°. So a sharp temperature spike levels out quickly at first and then more slowly as it fades. That is all “Fourier’s law” really says — no Greek letters required.
To put this on a computer, chop the rod into N evenly-spaced beads, each holding a temperature u_i. Every timestep, each interior bead nudges itself toward the average of its two neighbours. If a bead is already sitting at the neighbour average, nothing happens — it is already “in line”. If it is hotter than the average it cools; cooler, it warms. That nudge-toward-the- average is the discrete Laplacian, and smoothing — not wave-like propagation — is the only thing the heat equation ever does.
Make it concrete. Take five beads with temperatures (0, 10, 20, 30, 0) and a smoothing strength $r = 0.25$ r = 0.25. Focus on bead i = 3, currently at 30. Its neighbours are 20 and 0, so the discrete second derivative is $(20 - 2 \cdot 30 + 0) = - 40$ (20 - 2 \cdot 30 + 0) = -40 and the new value is $30 + 0.25 \cdot (- 40) = 20$ 30 + 0.25 \cdot (-40) = 20 — bead 3 dropped ten degrees in one step, precisely because it was sticking out way above the local average. Do this for every interior bead in parallel, leave the pinned endpoints alone, and you have advanced one timestep.
Written for any interior bead, the update rule is
$u_{i}^{n + 1} = u_{i}^{n} + r (u_{i + 1}^{n} - 2 u_{i}^{n} + u_{i - 1}^{n})$ u_i^{n+1} = u_i^n + r\,(u_{i+1}^n - 2 u_i^n + u_{i-1}^n)
where the superscript n labels the timestep and $r = α d t / d x^{2}$ r = \alpha\, dt / dx^2 is the smoothing strength per step. The bracketed term is the discrete second derivative — three samples reading off how “bent” the temperature profile is at that point. Concave-down hot peak: the bracket is negative and the bead cools. Concave-up valley: it warms. Curvature drives change here too, exactly as it did for the wave equation, but with a first time derivative rather than a second.
Shrink dx → 0 and dt → 0 and that update rule becomes the continuous heat equation:
$\frac{\partial u}{\partial t} = α \frac{\partial ^{2} u}{\partial x ^{2}}$ \frac{\partial u}{\partial t} = \alpha \, \frac{\partial^2 u}{\partial x^2}
The constant $α$ \alpha (thermal diffusivity) sets how fast smoothing happens — bigger for copper, smaller for wood. Sine-shaped modes are eigenfunctions: feed the rod a clean $sin (mπ x / L)$ \sin(m\pi x / L) bump and it decays as $exp (- α (mπ / L)^{2} t)$ \exp(-\alpha (m\pi/L)^2 t) while keeping its shape. Higher modes (more wiggles) decay much faster than the fundamental — the rod forgets fine detail first and broad bumps last.
There is one trap. The smoothing strength $r = α d t / d x^{2}$ r = \alpha\, dt / dx^2 must stay at or below $1/2$ 1/2 for this explicit forward Euler scheme — that is the CFL condition for diffusion. Step too big and the smoothing overshoots: a bead that ought to drift toward its neighbours’ average instead leaps past, then leaps back further next step, and the simulation explodes into ringing negative temperatures inside a few iterations. Halve dx for resolution and you must quarter dt — diffusion is unforgiving that way.
Try it yourself.

Load the Sine mode 3 preset, then the Gaussian bump preset. The mode-3 sine wave keeps its shape and just shrinks — it is an eigenfunction, so only its amplitude decays. The Gaussian bump instead spreads out and flattens, because it is a sum of many sine modes and the wigglier ones die first. By a few hundred steps the bump looks almost identical to the mode-1 sine; every initial condition relaxes toward the lowest mode and then quietly fades away.
That contrast with the 1D wave equation (string) demo just above is the whole point of putting the two side by side. The wave equation has a second time derivative — energy sloshes between potential and kinetic and is preserved, so a pluck rings forever. The heat equation has a first time derivative, so energy leaks out monotonically and the rod relaxes to silence. Same Laplacian on the right-hand side, very different behaviour on the left. The Laplace in a rectangle demo immediately below is what is left when even the time derivative is zero — the steady state the heat equation generalises toward as t → ∞.

t = 0.000 s

r = 0.240

max ∣ u ∣ = 0.999

1D heat equation with initial condition "Sine mode 1". α = 0.100, dt = 0.0010 s, stability r = 0.240.

Diffusivity α0.100

Thermal diffusivity. Larger α → faster smoothing toward zero.

Time step dt0.0010 s

FTCS step size. Stability requires r = α·dt/dx² ≤ 0.5.

Grid points nGrid50

Number of spatial samples on [0, L]. Smaller dx → tighter stability.

t = 0.000 s

Laplace in a rectangle

Iterative Laplace solver on a rectangular domain with mixed Dirichlet boundaries.

iterations: 357converged: ✓residual: 9.70e-7ω: 1.50

nGrid30

Grid resolution (used for both nx and ny).

ω (omega, SOR factor)1.50

Successive over-relaxation: 1 = Gauss-Seidel; ≈1.7 often fastest.

top (boundary)100

bottom (boundary)0

left (boundary)0

right (boundary)0

iter 357 · residual 9.70e-7 · converged ✓

Vector field visualizer

Type any 2D vector field; render its arrows, divergence, and curl.

max ∣ F ∣ = 1.044

max ∣\nabla \cdot F ∣ = 0.000

max ∣ (\nabla \times F)_{z} ∣ = 0.000

Coloring0

magnitude / divergence (red = source, blue = sink) / curl (red = CCW, blue = CW)

Field0

uniform · radial · vortex · saddle — current: uniform

nx16

Samples along x.

ny16

Samples along y.

max |F| 1.04 · max |div| 0.00 · max |curl| 0.00

RC / RL step response

First-order circuit step response paired with its frequency response.

τ = 1.00 e - 1 s

t_{90} = 2.30 e - 1 s

Mode (circuit)0

0 = RC charge · 1 = RC discharge · 2 = RL

R (resistance)100 Ω

Series resistance Ω.

C (capacitance)0.001000 F

Capacitance — only relevant for RC modes.

L (inductance)0.500 H

Inductance — only relevant for RL mode.

Vstep (step voltage)5 V

Magnitude of the step input.

τ = 1.000e-1 s · t90 = 2.303e-1 s

Biot-Savart B-field

Magnetic field from a current loop or coil, visualised as field lines and arrows.

Plan view · |Bz| heatmap (red out, blue in)

On-axis Bz(z) for circular loop · z ∈ [−3R, 3R]

I (current)1.00 A

Loop current in amperes.

R (loop radius)1.00 m

Loop radius / half-spacing in metres.

nSegments24

Number of straight segments approximating curved carriers.

grid resolution18

Heatmap grid cells per axis.

carrier geometry0

circular-loop (1/3)

max |Bz| = 0 T · on-axis peak 628.319 nT

Faraday’s law / induced EMF

Slide a magnet through a coil; live flux and induced EMF traces.

Rotating coil in uniform B

Φ(t) and ε(t) over one period

N (turns)1

Number of loops in the coil. EMF scales linearly with N.

B (field)0.50 T

Uniform field strength.

A (area)0.050 m²

Loop area enclosing the flux.

ω (omega, angular freq)60 rad/s

Rotation rate. Peak EMF = N·B·A·ω.

peak ε = 1.500 V · period T = 0.105 s

Featured problems

d'Alembert solution to the 1D wave equation

The 1D wave equation $u_{tt} = c^{2} u_{xx}$ u_{tt} = c^2 u_{xx} on the infinite line with initial displacement $u (x, 0) = φ (x)$ u(x,0) = \varphi(x) and initial velocity $u_{t} (x, 0) = ψ (x)$ u_t(x,0) = \psi(x). Solve in closed form and read off what kinds of disturbances propagate at what speed.

Change variables: $ξ = x - c t$ \xi = x - ct, $η = x + c t$ \eta = x + ct. The PDE becomes $u_{ξ η} = 0$ u_{\xi\eta} = 0, which integrates twice to give any sum of a right-moving and left-moving wave:

u (x, t) = F (x - c t) + G (x + c t)

Match the two initial conditions: $F (x) + G (x) = φ (x)$ F(x) + G(x) = \varphi(x) and $- c F^{'} (x) + c G^{'} (x) = ψ (x)$ -c F'(x) + c G'(x) = \psi(x). Solving for $F$ F and $G$ G yields the d’Alembert formula:

u (x, t) = \frac{1}{2} [φ (x - c t) + φ (x + c t)] + \frac{1}{2 c} \int_{x - c t}^{x + c t} ψ (s) d s

The geometry is striking: the value of $u$ u at $(x, t)$ (x,t) depends only on initial data inside the past light cone $[x - c t, x + c t]$ [x-ct,\, x+ct] — finite signal speed falls out of the math. An initial displacement splits into two half-amplitude copies travelling outward at speed $\pm c$ \pm c. An initial velocity profile contributes a smeared, integrated tail — which is why striking a drumhead (impulsive velocity) makes a longer-lived disturbance than plucking a string (instantaneous displacement). The Charge Field Plotter above tells the static side of the same story for electrostatics; the wave equation is where time enters and everything starts propagating.

Magnetic field on the axis of a circular current loop

A circular loop of radius $R$ R carries steady current $I$ I. Find the magnetic field $B (z)$ B(z) on the axis, a distance $z$ z from the loop’s centre.

The Biot–Savart law for a current element $I d ℓ$ I\, d\boldsymbol\ell at position $r^{'}$ \mathbf r' contributes

d B = (μ_{0} /4 π) I d ℓ \times (r - r^{'}) /∣ r - r^{'} ∣^{3}

at the field point. Each loop element is distance $R^{2} + z^{2}$ \sqrt{R^2 + z^2} from the axis point, and the cross product gives a vector of magnitude $R d ℓ$ R\, d\ell tilted off-axis by angle $θ$ \theta with $sin θ = R / R^{2} + z^{2}$ \sin\theta = R/\sqrt{R^2+z^2}. Symmetry kills the perpendicular components when integrated around the loop, leaving only the axial projection:

B_{z} (z) = \frac{μ _{0} I}{4 π} \oint \frac{R}{( R ^{2} + z ^{2} ) ^{3/2}} d ℓ = \frac{μ _{0} I R ^{2}}{2 ( R ^{2} + z ^{2} ) ^{3/2}}

Two limits worth checking. At the centre ( $z = 0$ z=0) $B = μ_{0} I / (2 R)$ B = \mu_0 I / (2R), the standard single-loop result. For $z ≫ R$ z \gg R, expanding gives $B \approx μ_{0} I R^{2} / (2 z^{3})$ B \approx \mu_0 I R^2 / (2 z^3) — the far field of any planar current loop is a magnetic dipole with moment $m = I π R^{2}$ m = I \pi R^2, and the $1/ z^{3}$ 1/z^3 fall-off is the signature of a dipole.

Stack two such loops separated by a distance equal to their radius (Helmholtz coils) and the cubic non-uniformity in the midplane cancels — the field is uniform to leading order over a surprisingly large central region. The Biot–Savart visualizer above shows this directly: drop two loops in coaxially, vary the separation, and watch the central field flatten near $d = R$ d = R.

Standing-wave frequencies on a fixed-end string by separation of variables

A string of length $L$ L with both ends pinned obeys the 1D wave equation $u_{tt} = c^{2} u_{xx}$ u_{tt} = c^2 u_{xx} with boundary conditions $u (0, t) = u (L, t) = 0$ u(0, t) = u(L, t) = 0 for all time. Show that the only frequencies that can stand on this string are $f_{n} = n c / (2 L)$ f_n = n c / (2 L) for $n = 1, 2, 3, \dots$ n = 1, 2, 3, \ldots.

Assume a separated solution $u (x, t) = X (x) T (t)$ u(x, t) = X(x) T(t). Substituting into the PDE and dividing by $X (x) T (t)$ X(x) T(t) gives

\frac{T ^{''} ( t )}{c ^{2} T ( t )} = \frac{X ^{''} ( x )}{X ( x )} = - λ

The left side depends only on $t$ t, the right only on $x$ x, so both must equal a constant — call it $- λ$ -\lambda with the sign chosen so non-decaying oscillation is possible. The spatial equation $X^{''} + λ X = 0$ X'' + \lambda X = 0 with $X (0) = X (L) = 0$ X(0) = X(L) = 0 is a classic Sturm–Liouville problem whose only non-trivial solutions are $X_{n} (x) = sin (nπ x / L)$ X_n(x) = \sin(n \pi x / L) with eigenvalues $λ_{n} = (nπ / L)^{2}$ \lambda_n = (n \pi / L)^2.

The temporal equation $T^{''} + c^{2} λ_{n} T = 0$ T'' + c^2 \lambda_n T = 0 oscillates at angular frequency $ω_{n} = c λ_{n} = nπ c / L$ \omega_n = c \sqrt{\lambda_n} = n \pi c / L, i.e. linear frequency $f_{n} = n c / (2 L)$ f_n = n c / (2 L). The fundamental is $f_{1} = c / (2 L)$ f_1 = c / (2 L); every higher mode is an integer multiple of it — the harmonic series that gives plucked-string instruments their consonant character.

Lab 3 used this exact relation in reverse. With $L$ L measured directly and $f_{1}$ f_1 read off the spectrum analyser, $c = 2 L f_{1}$ c = 2 L f_1 recovered the wave speed on the string to within a few percent of the value predicted from string density and tension. The wave-equation visualizer above lets you set fixed ends and pluck initial conditions, then watch the d’Alembert travelling waves bounce, recombine, and self- organise into the predicted standing-wave envelope.

Method of images — a point charge above a grounded conducting plane

A point charge $+ q$ +q sits a distance $d$ d above an infinite grounded conducting plane at $z = 0$ z = 0. Find the electrostatic potential in the upper half-space, the induced surface charge on the conductor, and the force pulling the real charge towards the plane. Solving Laplace’s equation directly with the boundary condition $Φ (x, y, 0) = 0$ \Phi(x, y, 0) = 0 is painful. Solving it by the method of images is two lines.

Replace the conductor with a fictitious image charge $- q$ -q at $(0, 0, - d)$ (0, 0, -d), in the empty half-space below. The superposition of real and image potentials,

Φ (x, y, z) = \frac{q}{4 π ε _{0}} [\frac{1}{x ^{2} + y ^{2} + ( z - d ) ^{2}} - \frac{1}{x ^{2} + y ^{2} + ( z + d ) ^{2}}]

is identically zero on $z = 0$ z = 0, satisfies Laplace’s equation everywhere in $z > 0$ z > 0 (the only singularity is the real charge), and matches the physical configuration at infinity. By the uniqueness theorem for Laplace’s equation with Dirichlet boundary data, this is the answer. The induced surface charge follows from $σ = - ε_{0} \partial Φ/ \partial z ∣_{z = 0} = - q d / [2 π (r^{2} + d^{2})^{3/2}]$ \sigma = -\varepsilon_0 \, \partial \Phi / \partial z|_{z=0} = -q d / [2\pi (r^2 + d^2)^{3/2}], which integrates to a total of $- q$ -q on the plane.

The force on the real charge is the force exerted by the image charge alone (the conductor is just a mathematical fiction): $F = - q^{2} / (16 π ε_{0} d^{2})$ F = -q^2 / (16 \pi \varepsilon_0 d^2), attractive, falling off as $1/ d^{2}$ 1/d^2 as if the charge were attracted by a single mirror partner. The wave- and field-visualisers on this page run the same uniqueness logic in motion — replace a boundary with the right virtual source, and the messy PDE collapses into a free-space superposition.

What’s coming

The author’s written reflection lands alongside the v7 interview. The demos and featured problems above are wired and playable today.