Skip to content
Computational Lab.

Heat Diffusion

How heat forgets its shape.

Intermediate / PDE / finite-difference / field

The heat equation says that a point warms or cools in proportion to how different it is from its neighbors, and the simulation applies exactly that rule on a grid. Sharp initials smooth into diffusion, and the field relaxes toward the steady state that the boundary conditions demand.

Tagged: parabolic PDE, finite differences, stability, steady state

Live simulation

Drag inside the field to inject heat.

Speed

Drag inside the field to inject heat.

Mathematics

The heat equation says each point warms or cools in proportion to how much its temperature differs from the local average. Discretized on a grid, it becomes a weighted average applied over and over.

∂u/∂t = α ∇²u

explicit FTCS: uⁿ⁺¹ = uⁿ + α (u_up + u_down + u_left + u_right − 4 uⁿ)

stability: α ≤ ¼ (2-D, unit cells)

That ¼ limit is not negotiable: beyond it the explicit scheme amplifies rounding noise into visible oscillation. The slider refuses to cross it and the metrics panel reports whether the clamp engaged. With insulated (Neumann) boundaries the total energy is conserved, which the total-energy metric verifies live; with edges held at zero, the field must relax to the empty state.

How it works

The field is a Float32Array updated in place through a second buffer, one FTCS sweep per solver step. Boundary mode mirrors neighbors (insulated) or pins them to zero. The render path writes the field through a 256-entry perceptual colormap into an ImageData buffer scaled up to the canvas.

Dragging inside the field injects heat under the pointer. Set the boundary to insulated, inject a hot spot, and watch total energy stay constant while the peak decays: conservation and diffusion, simultaneously visible.

Implementation

src/engine/experiments/heat-diffusion.ts · 205 lines · the module running above

TypeScript · fixed-timestep engine

1import { BaseExperiment } from "../core/base-experiment";
2import { MAGMA, sampleInto } from "../render/colormaps";
3import type { ParameterDef, PointerState, ThemeColors } from "../core/types";
4
5/**
6 * 2D heat equation via explicit FTCS finite differences. The stability limit
7 * (alpha <= 0.24 in cell units) is enforced by clamping and surfaced as a
8 * metric, never silently ignored.
9 */
10export class HeatDiffusionExperiment extends BaseExperiment {
11 readonly id = "heat-diffusion";
12
13 private n = 96;
14 private offsetX = 0;
15 private offsetY = 0;
16 private drawSize = 0;
17 private u = new Float32Array(0);
18 private next = new Float32Array(0);
19 private image: ImageData | null = null;
20 private scratch: HTMLCanvasElement | null = null;
21 private steps = 0;
22 private peakScale = 1;
23 private clamped = false;
24 private pointer = { x: 0, y: 0, down: false, inside: false };
25
26 protected params(): ParameterDef[] {
27 return [
28 {
29 key: "resolution",
30 label: "Grid",
31 options: [
32 { value: "64", label: "64 × 64" },
33 { value: "96", label: "96 × 96" },
34 { value: "128", label: "128 × 128" },
35 ],
36 defaultValue: "96",
37 },
38 { key: "alpha", label: "Diffusivity α", min: 0.02, max: 0.24, step: 0.01, defaultValue: 0.2 },
39 { key: "iterations", label: "Steps per tick", min: 1, max: 20, step: 1, defaultValue: 4 },
40 {
41 key: "boundary",
42 label: "Boundary",
43 options: [
44 { value: "insulated", label: "Insulated" },
45 { value: "fixed", label: "Held at zero" },
46 ],
47 defaultValue: "insulated",
48 },
49 {
50 key: "initial",
51 label: "Initial state",
52 options: [
53 { value: "blob", label: "Center blob" },
54 { value: "bar", label: "Left bar" },
55 { value: "random", label: "Random patches" },
56 ],
57 defaultValue: "blob",
58 },
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.u = new Float32Array(this.n * this.n);
69 this.next = new Float32Array(this.n * this.n);
70 const initial = this.str("initial");
71 if (initial === "blob") {
72 this.stampHeat(this.n / 2, this.n / 2, Math.max(3, this.n / 8), 1);
73 } else if (initial === "bar") {
74 for (let y = 0; y < this.n; y++) {
75 for (let x = 0; x < Math.max(2, this.n / 10); x++) {
76 this.u[y * this.n + x] = 1;
77 }
78 }
79 } else {
80 const patches = 8;
81 for (let p = 0; p < patches; p++) {
82 this.stampHeat(this.rng() * this.n, this.rng() * this.n, 2 + this.rng() * (this.n / 10), 1);
83 }
84 }
85 this.steps = 0;
86 this.peakScale = 1;
87 this.image = null;
88 }
89
90 private stampHeat(cx: number, cy: number, radius: number, value: number): void {
91 for (let y = Math.max(0, Math.floor(cy - radius)); y < Math.min(this.n, Math.ceil(cy + radius)); y++) {
92 for (let x = Math.max(0, Math.floor(cx - radius)); x < Math.min(this.n, Math.ceil(cx + radius)); x++) {
93 const d = Math.hypot(x - cx, y - cy);
94 if (d <= radius) {
95 const i = y * this.n + x;
96 this.u[i] = Math.max(this.u[i], value * (1 - d / radius));
97 }
98 }
99 }
100 }
101
102 protected onParameterChange(key: string): void {
103 if (key !== "iterations") this.reset();
104 }
105
106 resize(width: number, height: number): void {
107 this.width = width;
108 this.height = height;
109 this.reset();
110 }
111
112 onPointer(state: PointerState): void {
113 this.pointer = state;
114 if (state.inside && state.down) {
115 const gx = ((state.x - this.offsetX) / this.drawSize) * this.n;
116 const gy = ((state.y - this.offsetY) / this.drawSize) * this.n;
117 this.stampHeat(gx, gy, Math.max(2, this.n / 24), 1);
118 this.peakScale = 1;
119 }
120 }
121
122 protected onUpdate(): void {
123 const alpha = Math.min(this.num("alpha"), 0.24);
124 this.clamped = alpha < this.num("alpha") - 1e-9;
125 const n = this.n;
126 const u = this.u;
127 const next = this.next;
128 const insulated = this.str("boundary") === "insulated";
129 const iterations = this.num("iterations");
130 for (let it = 0; it < iterations; it++) {
131 for (let y = 0; y < n; y++) {
132 for (let x = 0; x < n; x++) {
133 const i = y * n + x;
134 // Insulated (Neumann) boundary: the ghost cell mirrors the cell's
135 // own value, which makes the stencil exactly conservative.
136 const up = y > 0 ? u[i - n] : insulated ? u[i] : 0;
137 const down = y < n - 1 ? u[i + n] : insulated ? u[i] : 0;
138 const left = x > 0 ? u[i - 1] : insulated ? u[i] : 0;
139 const right = x < n - 1 ? u[i + 1] : insulated ? u[i] : 0;
140 const laplacian = up + down + left + right - 4 * u[i];
141 next[i] = u[i] + alpha * laplacian;
142 }
143 }
144 u.set(next);
145 }
146 this.steps += iterations;
147 }
148
149 render(ctx: CanvasRenderingContext2D, theme: ThemeColors): void {
150 ctx.fillStyle = theme.bg;
151 ctx.fillRect(0, 0, this.width, this.height);
152 if (!this.image) this.image = new ImageData(this.n, this.n);
153 const n = this.n;
154 let max = 0;
155 for (let i = 0; i < this.u.length; i++) {
156 if (this.u[i] > max) max = this.u[i];
157 }
158 if (max > this.peakScale || max < this.peakScale * 0.4) this.peakScale = Math.max(max, 1e-4);
159 const data = this.image.data;
160 for (let i = 0; i < this.u.length; i++) {
161 sampleInto(MAGMA, this.u[i] / this.peakScale, data, i * 4);
162 }
163 if (!this.scratch) {
164 this.scratch = document.createElement("canvas");
165 this.scratch.width = n;
166 this.scratch.height = n;
167 }
168 const offCtx = this.scratch.getContext("2d")!;
169 offCtx.putImageData(this.image, 0, 0);
170 ctx.imageSmoothingEnabled = true;
171 ctx.drawImage(this.scratch, this.offsetX, this.offsetY, this.drawSize, this.drawSize);
172 }
173
174 getMetrics() {
175 let sum = 0;
176 let max = 0;
177 for (let i = 0; i < this.u.length; i++) {
178 sum += this.u[i];
179 if (this.u[i] > max) max = this.u[i];
180 }
181 const mean = sum / this.u.length;
182 return {
183 steps: this.steps,
184 peak: max.toFixed(3),
185 mean: mean.toFixed(4),
186 totalEnergy: sum.toFixed(1),
187 stability: this.clamped ? "α clamped to 0.24" : "unconditional OK",
188 };
189 }
190
191 describe(): string {
192 let sum = 0;
193 let max = 0;
194 for (let i = 0; i < this.u.length; i++) {
195 sum += this.u[i];
196 if (this.u[i] > max) max = this.u[i];
197 }
198 return `After ${this.steps.toLocaleString("en-US")} solver steps the field peaks at ${max.toFixed(2)} with mean ${ (sum / this.u.length).toFixed(3)}.`;
199 }
200
201 entities(): number {
202 return this.n * this.n;
203 }
204}
205

Performance

Time per tick
O(cells × steps-per-tick) with a 5-point stencil
Space
Two Float32Arrays plus one ImageData buffer
Notes
The solver applies the exact equation from the mathematics section: a point changes by α times its Laplacian, nothing more.

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.