Solving PDEs

Why it matters

Heat flow, diffusion, waves, electrostatics, fluid flow and the Schrödinger equation are all partial differential equations (PDEs): equations in several variables, usually space and time. Most real geometries have no closed-form solution, so physicists discretize them on a grid. The heat equation is the best place to learn this, because it is simple and it shows the main pitfall, numerical instability, very clearly.

The heat equation

For temperature u(x, t) along a thin rod,

du/dt = alpha * d^2u/dx^2

alpha is the thermal diffusivity in m^2/s (about 9.7e-5 for aluminium, 1.1e-4 for copper, about 1.4e-7 for water). The equation says a point warms up when it is colder than the average of its neighbours, so bumps smooth out over time.

Finite differences

Put the rod on a grid with spacing dx and step time by dt. Replace each derivative with a difference:

d^2u/dx^2  ~  (u[i+1] - 2*u[i] + u[i-1]) / dx^2
du/dt      ~  (u_new[i] - u[i]) / dt

Rearranging gives the explicit FTCS update (forward in time, centred in space):

r = alpha * dt / dx^2
u_new[i] = u[i] + r * (u[i+1] - 2*u[i] + u[i-1])

With NumPy slicing the whole interior updates in one line:

import numpy as np

def heat_step(u, r):
    new = u.copy()
    new[1:-1] = u[1:-1] + r * (u[2:] - 2 * u[1:-1] + u[:-2])
    return new          # end points unchanged: fixed-temperature boundaries

The stability limit

The explicit scheme is only stable if

r = alpha * dt / dx^2 <= 1/2

Break this and small round-off ripples grow exponentially until the numbers explode, even though each step looks harmless. The practical consequence: halving dx forces you to cut dt by four. For a fine grid, implicit methods (backward Euler, Crank-Nicolson) are stable for any dt at the cost of solving a linear system each step; SciPy's sparse solvers handle that in real projects.

Worked example: a rod between two baths

An aluminium rod 0.5 m long starts at 20 °C. At t = 0 one end is held at 100 °C and the other at 0 °C. With 51 grid points, dx = 0.01 m. The largest stable time step is

dt_max = 0.5 * dx^2 / alpha = 0.5 * 1e-4 / 9.7e-5 = 0.52 s

Using dt = 0.4 s (r = 0.39), the profile bends toward a straight line from 100 to 0 °C, which is the steady state (where d^2u/dx^2 = 0). The slowest-decaying part of the starting error has a time constant of L^2 / (pi^2 alpha) = 0.25 / (9.87 * 9.7e-5), about 260 s, so after about 20 minutes the rod is within roughly a degree of the straight line.

Testing against an exact solution

Always test a PDE solver on a case you can solve by hand. With both ends held at 0 and a starting profile u = sin(pi x / L), the exact solution keeps its shape and decays:

u(x, t) = sin(pi x / L) * exp(-alpha * pi^2 * t / L^2)

Run the solver to some time t, compare the peak with exp(-alpha pi^2 t / L^2), and you have a quantitative check. Refining the grid should shrink the difference.

Common mistakes

• Ignoring the stability condition and blaming the physics when values blow up.
• Updating in place so new values feed into neighbours in the same step. Compute from a copy (or the slicing form above).
• Overwriting boundary values each step. Decide the boundary condition (fixed temperature, insulated, ...) and apply it explicitly.
• Mixing units, for example dx in cm with alpha in m^2/s.
• Trusting one grid. Repeat with half the spacing; if the answer moves, refine further.

Recap

• The heat equation u_t = alpha u_xx smooths bumps over time.
• FTCS: u_new = u + r (u[i+1] - 2u[i] + u[i-1]) with r = alpha dt / dx^2.
• Stable only for r <= 1/2; implicit methods remove the limit.
• Check solvers against exact solutions and steady states.

Try it

Open the Do tab to build an explicit heat-equation solver with NumPy, including a stability guard and an exact-solution test.