Electrostatic Potential
Charges in a grounded box, solved honestly.
Advanced / PDE / field / linear-solver / relaxation
Placing charges stamps a charge density onto a grid, and the potential comes from solving the Poisson equation with the box walls held at zero. Successive over-relaxation - Gauss-Seidel with momentum - runs a few sweeps per tick, so you watch the field propagate outward from each new charge, and the residual tells you how close to solved the field is. The banded rendering is the equipotential structure.
Tagged: Poisson equation, SOR, equipotential, Gauss-Seidel
Live simulation
Click to place charges with the selected tool; the field re-relaxes live. Keyboard: Space runs or pauses, . steps, r resets.
Speed
Click to place charges with the selected tool; the field re-relaxes live.
Mathematics
Charges on a grounded box solve the Poisson equation: the potential's Laplacian is the (negative) charge density, with the walls pinned at zero volts. Discretized, every interior cell must satisfy a weighted average of its neighbors plus its own source term.
∇²V = −ρ/ε, V|_boundary = 0
discrete: 4 V_ij − V_ij−1 − V_ij+1 − V_i−1j − V_i+1j = ρ_ij
SOR: V ← V + ω (ρ + Σ V_neighbors − 4V) / 4, ω ≈ 1.75
The solver is Gauss-Seidel with relaxation (ω > 1 accelerates it), red-black ordered so it needs no second buffer. The sweep loop stops when the largest per-cell update drops below a residual threshold tied to the charge scale - and that residual is displayed, not assumed. The banded color mapping quantizes the potential into 16 levels: the band edges are equipotential contours.
How it works
Clicking stamps a charge into the density array as a smooth Gaussian blob (sharp point charges make the discrete system unhappy), then resets the sweep counter; the field you see is the solver physically relaxing toward the solution a few sweeps per tick.
The tool selector chooses positive, negative, or erase; erase finds the nearest charge within radius. The same linear system is what Conjugate Gradient solves - there three solvers race on it; here SOR animates it.
Implementation
src/engine/experiments/electrostatic-potential.ts · 232 lines · the module running above
TypeScript · fixed-timestep engine
| 1 | import { BaseExperiment } from "../core/base-experiment"; |
| 2 | import { sampleInto } from "../render/colormaps"; |
| 3 | import type { ParameterDef, PointerState, ThemeColors } from "../core/types"; |
| 4 | import { divergingFor } from "../render/diverging-cache"; |
| 5 | |
| 6 | interface Charge { |
| 7 | x: number; |
| 8 | y: number; |
| 9 | q: number; // +1 / -1 |
| 10 | } |
| 11 | |
| 12 | /** |
| 13 | * Electrostatics as a Poisson problem: 4V - sum(neighbors) = rho on a grid |
| 14 | * with the box walls grounded (V = 0). Charges stamp rho; successive |
| 15 | * over-relaxation (Gauss-Seidel + momentum) solves the same linear system |
| 16 | * the Conjugate Gradient experiment races on, a few sweeps per tick so the |
| 17 | * field is seen propagating, and the residual is reported, not hidden. |
| 18 | */ |
| 19 | export class ElectrostaticPotentialExperiment extends BaseExperiment { |
| 20 | readonly id = "electrostatic-potential"; |
| 21 | |
| 22 | private n = 80; |
| 23 | private drawSize = 0; |
| 24 | private offsetX = 0; |
| 25 | private offsetY = 0; |
| 26 | private v = new Float32Array(0); |
| 27 | private rho = new Float32Array(0); |
| 28 | private image: ImageData | null = null; |
| 29 | private scratch: HTMLCanvasElement | null = null; |
| 30 | private charges: Charge[] = []; |
| 31 | private sweeps = 0; |
| 32 | private residual = Infinity; |
| 33 | private pointer = { x: 0, y: 0, down: false, inside: false }; |
| 34 | private lastCell = { x: -1, y: -1 }; |
| 35 | |
| 36 | protected params(): ParameterDef[] { |
| 37 | return [ |
| 38 | { |
| 39 | key: "tool", |
| 40 | label: "Pointer tool", |
| 41 | options: [ |
| 42 | { value: "positive", label: "Place + charge" }, |
| 43 | { value: "negative", label: "Place - charge" }, |
| 44 | { value: "erase", label: "Erase nearest" }, |
| 45 | ], |
| 46 | defaultValue: "positive", |
| 47 | }, |
| 48 | { |
| 49 | key: "resolution", |
| 50 | label: "Grid", |
| 51 | options: [ |
| 52 | { value: "64", label: "64 × 64" }, |
| 53 | { value: "80", label: "80 × 80" }, |
| 54 | { value: "96", label: "96 × 96" }, |
| 55 | ], |
| 56 | defaultValue: "80", |
| 57 | }, |
| 58 | { key: "omega", label: "SOR relaxation w", min: 1.2, max: 1.9, step: 0.05, defaultValue: 1.75 }, |
| 59 | ]; |
| 60 | } |
| 61 | |
| 62 | protected onReset(): void { |
| 63 | this.n = Number(this.str("resolution")); |
| 64 | const cell = Math.min(this.width, this.height) / this.n; |
| 65 | this.drawSize = cell * this.n; |
| 66 | this.offsetX = Math.floor((this.width - this.drawSize) / 2); |
| 67 | this.offsetY = Math.floor((this.height - this.drawSize) / 2); |
| 68 | this.v = new Float32Array(this.n * this.n); |
| 69 | this.rho = new Float32Array(this.n * this.n); |
| 70 | // Default scene: a dipole, so the first frame already teaches something. |
| 71 | this.charges = [ |
| 72 | { x: 0.35, y: 0.5, q: 1 }, |
| 73 | { x: 0.65, y: 0.5, q: -1 }, |
| 74 | ]; |
| 75 | this.sweeps = 0; |
| 76 | this.residual = Infinity; |
| 77 | this.stampRho(); |
| 78 | } |
| 79 | |
| 80 | protected onParameterChange(key: string): void { |
| 81 | if (key === "resolution") this.reset(); |
| 82 | } |
| 83 | |
| 84 | resize(width: number, height: number): void { |
| 85 | this.width = width; |
| 86 | this.height = height; |
| 87 | this.reset(); |
| 88 | } |
| 89 | |
| 90 | private stampRho(): void { |
| 91 | this.rho.fill(0); |
| 92 | for (const charge of this.charges) { |
| 93 | const cx = charge.x * this.n; |
| 94 | const cy = charge.y * this.n; |
| 95 | const radius = Math.max(1.5, this.n / 32); |
| 96 | for (let y = Math.max(1, Math.floor(cy - 2 * radius)); y < Math.min(this.n - 1, Math.ceil(cy + 2 * radius)); y++) { |
| 97 | for (let x = Math.max(1, Math.floor(cx - 2 * radius)); x < Math.min(this.n - 1, Math.ceil(cx + 2 * radius)); x++) { |
| 98 | const d = Math.hypot(x - cx, y - cy); |
| 99 | if (d <= 2 * radius) { |
| 100 | this.rho[y * this.n + x] += (charge.q * 8) * Math.exp(-(d * d) / (radius * radius)); |
| 101 | } |
| 102 | } |
| 103 | } |
| 104 | } |
| 105 | this.sweeps = 0; |
| 106 | this.residual = Infinity; |
| 107 | } |
| 108 | |
| 109 | onPointer(state: PointerState): void { |
| 110 | this.pointer = state; |
| 111 | if (!state.inside || !state.down) { |
| 112 | this.lastCell = { x: -1, y: -1 }; |
| 113 | return; |
| 114 | } |
| 115 | const gx = ((state.x - this.offsetX) / this.drawSize) * this.n; |
| 116 | const gy = ((state.y - this.offsetY) / this.drawSize) * this.n; |
| 117 | const cell = { x: Math.round(gx), y: Math.round(gy) }; |
| 118 | if (cell.x === this.lastCell.x && cell.y === this.lastCell.y) return; |
| 119 | this.lastCell = cell; |
| 120 | const tool = this.str("tool"); |
| 121 | if (tool === "erase") { |
| 122 | let best = -1; |
| 123 | let bestDist = 0.08; |
| 124 | for (let i = 0; i < this.charges.length; i++) { |
| 125 | const d = Math.hypot(this.charges[i].x - gx / this.n, this.charges[i].y - gy / this.n); |
| 126 | if (d < bestDist) { |
| 127 | bestDist = d; |
| 128 | best = i; |
| 129 | } |
| 130 | } |
| 131 | if (best >= 0) this.charges.splice(best, 1); |
| 132 | } else if (cell.x > 0 && cell.y > 0 && cell.x < this.n - 1 && cell.y < this.n - 1) { |
| 133 | this.charges.push({ x: gx / this.n, y: gy / this.n, q: tool === "positive" ? 1 : -1 }); |
| 134 | if (this.charges.length > 40) this.charges.shift(); |
| 135 | } |
| 136 | this.stampRho(); |
| 137 | } |
| 138 | |
| 139 | protected onUpdate(): void { |
| 140 | const n = this.n; |
| 141 | const v = this.v; |
| 142 | const omega = this.num("omega"); |
| 143 | let maxDelta = 0; |
| 144 | const sweepsPerTick = 8; |
| 145 | // Red-black SOR keeps the update in place without a second buffer. |
| 146 | for (let sweep = 0; sweep < sweepsPerTick; sweep++) { |
| 147 | if (this.residual < 1e-5 * Math.max(1, this.totalAbsRho() / (n * n))) break; |
| 148 | maxDelta = 0; |
| 149 | for (let color = 0; color < 2; color++) { |
| 150 | for (let y = 1; y < n - 1; y++) { |
| 151 | for (let x = 1 + ((y + color) % 2); x < n - 1; x += 2) { |
| 152 | const i = y * n + x; |
| 153 | const sigma = v[i - 1] + v[i + 1] + v[i - n] + v[i + n]; |
| 154 | const update = (omega * (this.rho[i] + sigma - 4 * v[i])) / 4; |
| 155 | v[i] += update; |
| 156 | maxDelta = Math.max(maxDelta, Math.abs(update)); |
| 157 | } |
| 158 | } |
| 159 | } |
| 160 | this.sweeps += 1; |
| 161 | this.residual = maxDelta; |
| 162 | } |
| 163 | } |
| 164 | |
| 165 | private totalAbsRho(): number { |
| 166 | let sum = 0; |
| 167 | for (let i = 0; i < this.rho.length; i++) sum += Math.abs(this.rho[i]); |
| 168 | return sum; |
| 169 | } |
| 170 | |
| 171 | render(ctx: CanvasRenderingContext2D, theme: ThemeColors): void { |
| 172 | ctx.fillStyle = theme.bg; |
| 173 | ctx.fillRect(0, 0, this.width, this.height); |
| 174 | if (!this.image || this.image.width !== this.n) { |
| 175 | this.image = new ImageData(this.n, this.n); |
| 176 | this.scratch = null; |
| 177 | } |
| 178 | const diverging = divergingFor(theme); |
| 179 | // Quantized into 16 levels: the band edges are the equipotential contours. |
| 180 | let maxAbs = 1e-6; |
| 181 | for (let i = 0; i < this.v.length; i++) maxAbs = Math.max(maxAbs, Math.abs(this.v[i])); |
| 182 | const data = this.image.data; |
| 183 | for (let i = 0; i < this.v.length; i++) { |
| 184 | const normalized = this.v[i] / maxAbs; |
| 185 | const banded = Math.round(normalized * 8) / 8; |
| 186 | sampleInto(diverging, Math.min(1, Math.max(0, 0.5 + banded * 0.5)), data, i * 4); |
| 187 | } |
| 188 | if (!this.scratch) { |
| 189 | this.scratch = document.createElement("canvas"); |
| 190 | this.scratch.width = this.n; |
| 191 | this.scratch.height = this.n; |
| 192 | } |
| 193 | const offCtx = this.scratch.getContext("2d")!; |
| 194 | offCtx.putImageData(this.image, 0, 0); |
| 195 | ctx.imageSmoothingEnabled = true; |
| 196 | ctx.drawImage(this.scratch, this.offsetX, this.offsetY, this.drawSize, this.drawSize); |
| 197 | // Charges |
| 198 | for (const charge of this.charges) { |
| 199 | const px = this.offsetX + charge.x * this.drawSize; |
| 200 | const py = this.offsetY + charge.y * this.drawSize; |
| 201 | ctx.fillStyle = charge.q > 0 ? theme.accent : theme.viz[1]; |
| 202 | ctx.beginPath(); |
| 203 | ctx.arc(px, py, 5, 0, Math.PI * 2); |
| 204 | ctx.fill(); |
| 205 | ctx.fillStyle = theme.bg; |
| 206 | ctx.font = "bold 9px ui-monospace, monospace"; |
| 207 | ctx.textAlign = "center"; |
| 208 | ctx.textBaseline = "middle"; |
| 209 | ctx.fillText(charge.q > 0 ? "+" : "\u2212", px, py + 0.5); |
| 210 | } |
| 211 | } |
| 212 | |
| 213 | getMetrics() { |
| 214 | return { |
| 215 | charges: this.charges.length, |
| 216 | netQ: this.charges.reduce((sum, c) => sum + c.q, 0), |
| 217 | sweeps: this.sweeps.toLocaleString("en-US"), |
| 218 | residual: this.residual === Infinity ? "n/a" : this.residual.toExponential(2), |
| 219 | status: this.residual < 1e-5 ? "Converged" : "Solving", |
| 220 | }; |
| 221 | } |
| 222 | |
| 223 | describe(): string { |
| 224 | const net = this.charges.reduce((sum, c) => sum + c.q, 0); |
| 225 | return `${this.charges.length} charges (net ${net >= 0 ? "+" : ""}${net}) on a grounded box; potential solved with ${this.sweeps.toLocaleString("en-US")} SOR sweeps, residual ${this.residual === Infinity ? "n/a" : this.residual.toExponential(1)}.`; |
| 226 | } |
| 227 | |
| 228 | entities(): number { |
| 229 | return this.n * this.n; |
| 230 | } |
| 231 | } |
| 232 | |
Performance
- Time per tick
- O(cells) per sweep, 8 sweeps per tick until converged
- Space
- Two float arrays: potential and charge density
- Notes
- Grounded walls are the boundary condition: V is never updated on the outer ring, so field lines have somewhere real to terminate.
Live numbers (FPS, frame time, steps per second, microseconds per step, entity count) are in the performance overlay on the simulation above: the pulse icon in its top-right corner.