The heat equation: adding time

Finite Elements

Poisson's equation has no time. Heat does. We need a new kind of matrix — the mass matrix — and we need to decide how to march the solution forward in time. Both decisions have visible consequences, and the visualizer at the end of the chapter lets you feel them.

The PDE and its weak form

The 1D heat equation:

\frac{\partial u}{\partial t} = α \frac{\partial ^{2} u}{\partial x ^{2}}, x \in (0, 1), t > 0

With Dirichlet boundary conditions $u (0, t) = u (1, t) = 0$ and an initial profile $u (x, 0) = u_{0} (x)$ . The constant $α > 0$ is the thermal diffusivity. Multiply by a test function and integrate:

\int_{0}^{1} \frac{\partial u}{\partial t} v d x + α \int_{0}^{1} u^{'} v^{'} d x = 0 for all admissible v .

Same integration-by-parts trick as Poisson on the second term; the time derivative just rides along. Now expand $u (x, t) = \sum_{j} c_{j} (t) ϕ_{j} (x)$ — coefficients that depend on time but basis functions that don't.

The mass matrix

Plug the expansion in. The first integral becomes:

j \sum \overset{c}{˙}_{j} (t) \int_{0}^{1} ϕ_{i} (x) ϕ_{j} (x) d x = j \sum M_{ij} \overset{c}{˙}_{j} (t)

where $M_{ij} = \int ϕ_{i} ϕ_{j} d x$ is the mass matrix. It plays the same role for the time derivative that $K$ played for the spatial Laplacian.

For 1D hat functions on a uniform mesh of width $h$ , the integrals are short to compute by hand. Each interior node's self-overlap integrates to $2 h /3$ ; each adjacent-node overlap integrates to $h /6$ ; everything else is zero (different basis functions don't share elements):

M = \frac{h}{6} 41141 14 ⋱ 1 ⋱

Tridiagonal, symmetric, positive-definite. Same shape as $K$ , different entries.

The semi-discrete system

Putting $M$ and $K$ together gives a system of ordinary differential equations in time:

M \overset{c}{˙} (t) + α K c (t) = 0

"Semi-discrete" because we discretized in space (FEM) but not yet in time. The next step picks a numerical method for the time derivative; the choice has a name and a stability story.

Three time-stepping schemes

Approximate $\overset{c}{˙} (t) \approx (c^{n + 1} - c^{n}) /Δ t$ and decide where to evaluate the right-hand side: at the old time, the new time, or somewhere between.

Forward Euler (explicit)

\frac{c ^{n + 1} - c ^{n}}{Δ t} = - α M^{- 1} K c^{n}

The right-hand side uses only the known $c^{n}$ , so the update is direct — solve $M c^{n + 1} = M c^{n} - α Δ t K c^{n}$ once per step. Conditionally stable: blows up unless $Δ t \leq h^{2} / (2 α)$ . Halving the mesh size requires quartering the time step. Painful for fine meshes.

def step_forward_euler(c, M, K, alpha, dt):
    """Explicit. One step: c_new = c - α dt M^-1 K c.
    Stable only when dt < h² / (2α)."""
    rhs = M @ c - alpha * dt * (K @ c)
    return solve_tri(M, rhs)

Backward Euler (implicit)

\frac{c ^{n + 1} - c ^{n}}{Δ t} = - α M^{- 1} K c^{n + 1}

The right-hand side uses the unknown $c^{n + 1}$ , so we solve $(M + α Δ t K) c^{n + 1} = M c^{n}$ at every step. Unconditionally stable. First-order accurate in time, like forward Euler.

def step_backward_euler(c, M, K, alpha, dt):
    """Implicit. (M + α dt K) c_new = M c. Unconditionally stable."""
    A = M + alpha * dt * K
    rhs = M @ c
    return solve_tri(A, rhs)

Crank-Nicolson

\frac{c ^{n + 1} - c ^{n}}{Δ t} = - \frac{α}{2} M^{- 1} K (c^{n + 1} + c^{n})

Average of the two. Solve $(M + \frac{1}{2} α Δ t K) c^{n + 1} = (M - \frac{1}{2} α Δ t K) c^{n}$ . Unconditionally stable AND second-order accurate in time. The default choice when accuracy matters and you can afford an implicit solve per step.

def step_crank_nicolson(c, M, K, alpha, dt):
    """Implicit, 2nd-order in time. (M + α dt K/2) c_new = (M - α dt K/2) c."""
    A = M + 0.5 * alpha * dt * K
    rhs = (M - 0.5 * alpha * dt * K) @ c
    return solve_tri(A, rhs)

All three solve a tridiagonal system per step in 1D — same Thomas-algorithm cost. The difference is what the matrix on the left looks like, and whether it depends on $Δ t$ .

See it

The visualizer evolves an initial profile under each scheme. The dashed gray line shows the initial condition for reference; the solid accent curve is the current solution. Switch schemes, switch initial conditions, push $Δ t$ past the CFL bound on forward Euler and watch the explicit scheme blow up while the implicit ones stay stable.

scheme

Ndt

CFL bound (FE): dt < h² / 2α = 0.0098

Things to read off the animation:

Diffusion smooths. Sharp features — the step, the bumps — round into broader, lower humps. The total integral decreases over time because heat leaks out through the Dirichlet boundaries.
The CFL bound is real. With forward Euler at $N = 32$ and the largest $Δ t$ , the curve oscillates wildly and diverges within a few steps. Switch to backward Euler or Crank-Nicolson at the same step size and the solution is fine.
Backward Euler is dissipative. Even well-resolved features lose amplitude faster under BE than under Crank-Nicolson — BE has a built-in numerical dissipation that adds to the physical one. CN preserves the right amount of energy decay.
The two-bump initial condition shows the bumps merging as their support diffuses outward. Eventually only the lowest mode survives, decaying as $e^{- α π^{2} t} sin (π x)$ .

Where this is going

We've added time. The next chapter adds vector fields: linear elasticity turns the scalar PDE into a vector PDE, with two displacement components per node and a stress tensor that needs visualizing. The mass matrix story carries over to dynamic elasticity (wave equation, basically); the Crank-Nicolson generalization is called the Newmark scheme, and it's the workhorse of time-domain structural dynamics. Static elasticity, which is what we'll build next, drops the time derivative — it's Poisson's equation again, just vectorial.