The Heat Equation: Solving PDEs with Finite Differences

A conservation law in disguise

The heat equation describes how temperature spreads through a material once you stop asking about individual molecules and start asking about a continuous field. It follows from two ideas glued together. Fourier’s law says heat flux is proportional to the temperature gradient: heat flows from hot to cold, faster where the gradient is steeper. Energy conservation says the temperature at a point only changes if more heat is flowing in than out. Combine the two and you get a single partial differential equation, with thermal diffusivity α = k/(ρc) bundling conductivity k, density ρ and specific heat c into one number.

∂u/∂t = α ( ∂²u/∂x² + ∂²u/∂y² )    ← 2D heat equation
u(x,y,t)  = temperature field
α = k / (ρ c)  = thermal diffusivity

This is the same equation that governs diffusion of ink in water and the spread of a pollutant in the atmosphere — diffusion and conduction are the same mathematics with different labels on the axes.

Turning derivatives into neighbours

A computer cannot evaluate a derivative directly, so the simulation lays a grid over the plate and replaces each derivative with a difference between neighbouring cells. The standard central-difference approximation for the second derivative in one direction is (u[i+1] − 2u[i] + u[i−1]) / Δx²; add the same thing in the other direction and you get an update rule that moves every cell toward the average of its four neighbours, weighted by how far out of step it currently is.

next[i][j] = cur[i][j] + α * dt * (
  (cur[i+1][j] - 2*cur[i][j] + cur[i-1][j]) / dx2 +
  (cur[i][j+1] - 2*cur[i][j] + cur[i][j-1]) / dy2
)   // explicit FTCS scheme, one full grid pass per frame

Why the step size cannot be arbitrary

This explicit scheme (forward-time, centred-space, or FTCS) is cheap but not unconditionally safe. Von Neumann stability analysis — checking how each Fourier mode of the error grows from one step to the next — shows the scheme only stays bounded if α·dt·(1/dx² + 1/dy²) stays at or below one half. Push dt past that line and the grid does not just get slightly wrong, it oscillates and diverges within a handful of frames, because dx already appears squared in the denominator: halving the cell size forces the time step down by a factor of four. Implicit schemes such as Crank–Nicolson solve a linear system every step instead of reading neighbours directly, which costs more per step but removes the stability limit entirely, so large time steps stay well behaved.

Transient versus steady state

Running the update rule forward in time gives the transient solution — the plate cooling and warming as heat sources switch on. Left alone long enough with fixed boundary temperatures, the field stops changing and satisfies Laplace’s equation, ∇²u = 0: the steady state. You can reach it two ways — step the transient equation until it stops moving, which can take thousands of frames, or solve the steady equation directly with a relaxation method such as Gauss–Seidel or successive over-relaxation, which iterates purely in space and converges far faster when only the final picture matters.

Boundary conditions decide the picture

Painting a fixed hot or cold spot pins a cell’s temperature every step — a Dirichlet condition. An insulated edge instead forces the gradient across it to zero — a Neumann condition — so heat reflects back inward rather than leaking out. Real walls do a mix of both, losing heat in proportion to the temperature difference with the outside air, which is the Robin condition, the one actual radiators and windows obey.

Frequently asked questions

Try it live

Everything above runs in your browser — open Heat Equation and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.