Time-Dependent Heat Conduction

Method of lines and implicit time stepping

Run AllReset

The heat equation describes how temperature diffuses through a material over time. In this tutorial we solve it with the method of lines: discretize space first (yielding a large ODE system), then integrate in time with an implicit scheme that is stable for any step size. Everything runs in your browser.

Problem: Solve $\dfrac{\partial u}{\partial t} = \Delta u$ on the unit square with $u = 0$ on the boundary, starting from a hot Gaussian blob in the center.

Step 1: Spatial Discretization

As in the Poisson tutorial, we use a structured grid with $n$ interior nodes per direction, spacing $h = 1/(n+1)$ , and the 5-point Laplacian stencil. Semi-discretizing in space turns the PDE into a system of ordinary differential equations — one per node:

$\frac{du}{dt} = -\frac{1}{h^2} L \, u$

where $L$ is the (positive definite) stiffness matrix.

Code [1]Run

# Grid parameters
n = 16                 # interior nodes per direction
h = 1.0 / (n + 1)
N = n * n

println("=== Spatial grid ===")
println("Nodes: ${n} x ${n} = ${N} unknowns")
println("Spacing h = ${h}")

Step 2: Assemble the Laplacian

Matrix entries are collected as (row, column, value) triplets; sparse sums duplicates, matching FEM assembly semantics.

Code [2]Run

# 5-point stencil Laplacian L (the 1/h^2 factor is applied in the time stepper)
rows = zeros(5 * N)
cols = zeros(5 * N)
vals = zeros(5 * N)

idx = 1
for j in 1:n {
  for i in 1:n {
    k = (j - 1) * n + i
    rows[idx] = k
    cols[idx] = k
    vals[idx] = 4.0
    idx = idx + 1
    if i > 1 {
      rows[idx] = k
      cols[idx] = k - 1
      vals[idx] = -1.0
      idx = idx + 1
    }
    if i < n {
      rows[idx] = k
      cols[idx] = k + 1
      vals[idx] = -1.0
      idx = idx + 1
    }
    if j > 1 {
      rows[idx] = k
      cols[idx] = k - n
      vals[idx] = -1.0
      idx = idx + 1
    }
    if j < n {
      rows[idx] = k
      cols[idx] = k + n
      vals[idx] = -1.0
      idx = idx + 1
    }
  }
}

used = idx - 1
r2 = zeros(used)
c2 = zeros(used)
v2 = zeros(used)
for t in 1:used {
  r2[t] = rows[t]
  c2[t] = cols[t]
  v2[t] = vals[t]
}

L = sparse(r2, c2, v2, N, N)
println("Laplacian assembled: ${N} x ${N}, ${nnz(L)} nonzeros")

Step 3: Initial Condition

A Gaussian temperature blob centered at $(0.5, 0.5)$ :

Code [3]Run

# u0(x, y) = exp(-50 ((x-0.5)^2 + (y-0.5)^2))
u = zeros(N)
for j in 1:n {
  for i in 1:n {
    k = (j - 1) * n + i
    x = i * h
    y = j * h
    d2 = (x - 0.5) * (x - 0.5) + (y - 0.5) * (y - 0.5)
    u[k] = exp(-50.0 * d2)
  }
}

heat0 = sum(u) * h * h
println("Initial temperature field set")
println("Peak: ${max(u)}")
println("Total heat: ${heat0}")

Step 4: Why Implicit? Stability

The simplest time stepper, explicit (forward) Euler, is only stable for tiny steps: $\Delta t \le h^2/4$ — for our grid about $0.0009$ . Backward Euler instead solves

$\left(I + \frac{\Delta t}{h^2} L\right) u^{k+1} = u^k$

at each step. It costs a sparse solve per step, but is stable for any $\Delta t$ — this is why production codes use implicit integrators for diffusion.

Code [4]Run

dt = 0.002                     # ~2x above the explicit stability limit
alpha = dt / (h * h)

# System matrix M = I + alpha * L (triplets again; duplicates sum)
mr = zeros(used + N)
mc = zeros(used + N)
mv = zeros(used + N)
for t in 1:used {
  mr[t] = r2[t]
  mc[t] = c2[t]
  mv[t] = alpha * v2[t]
}
for k in 1:N {
  mr[used + k] = k
  mc[used + k] = k
  mv[used + k] = 1.0
}
M = sparse(mr, mc, mv, N, N)

println("Time step dt = ${dt} (explicit limit: ${h * h / 4.0})")
println("System matrix ready: ${nnz(M)} nonzeros")

Step 5: March in Time

Each step is one conjugate-gradient solve. Watch the temperature peak decay as heat drains through the cold boundary.

Code [5]Run

nsteps = 10

println("  step    time      peak")
for s in 1:nsteps {
  (u, flag, relres, iters) = pcg(M, u, 1e-10, 500)
  println("   ${s}     ${s * dt}    ${max(u)}")
}

heat = sum(u) * h * h
println("Total heat: ${heat} (started at ${heat0})")

The peak decreases monotonically — diffusion smooths the blob out, with no oscillations even though our step exceeds the explicit stability limit.

Step 6: Visualize the Final Temperature Field

Code [6]Run

# Embed into the padded (n+2) x (n+2) grid including the cold boundary
m = n + 2
U = zeros(m * m)
for j in 1:n {
  for i in 1:n {
    k = (j - 1) * n + i
    g = j * m + (i + 1)
    U[g] = u[k]
  }
}
Ugrid = reshape(U, m, m)

xs = linspace(0.0, 1.0, m)
(X, Y) = meshgrid(xs, xs)

surf(X, Y, Ugrid)
title("Temperature after ${nsteps} implicit steps")

Step 7: Fast-Forward to Steady State

As $t \to \infty$ the blob diffuses away completely: with no heat source and cold boundaries, the steady state is $u \equiv 0$ . Because backward Euler is unconditionally stable, we can fast-forward with 10× larger steps:

Code [7]Run

dt2 = 0.02
alpha2 = dt2 / (h * h)

mr2 = zeros(used + N)
mc2 = zeros(used + N)
mv2 = zeros(used + N)
for t in 1:used {
  mr2[t] = r2[t]
  mc2[t] = c2[t]
  mv2[t] = alpha2 * v2[t]
}
for k in 1:N {
  mr2[used + k] = k
  mc2[used + k] = k
  mv2[used + k] = 1.0
}
M2 = sparse(mr2, mc2, mv2, N, N)

for s in 1:5 {
  (u, flag, relres, iters) = pcg(M2, u, 1e-10, 500)
}

println("Peak after fast-forward: ${max(u)}")

Summary

The method-of-lines recipe for time-dependent PDEs:

Discretize space — sparse stiffness matrix from triplets, exactly as in the stationary case
Pick a time integrator — implicit (backward Euler) for diffusion: one sparse solve per step, unconditionally stable
March — pcg solves each step's SPD system in a few dozen iterations
Monitor physics — peak decay and total heat give immediate sanity checks

Swap the Laplacian for a different operator (variable conductivity, convection) and only the assembly cell changes — the time-stepping machinery stays identical.

Workbench

Clear

No variables in workbench

Next Steps

Finite Elements 1: Poisson Problem

The stationary problem behind this tutorial

Sparse Matrices

The sparse toolbox used for assembly and solving