Time-Dependent Heat Conduction
Method of lines and implicit time stepping
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 on the unit square with 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 interior nodes per direction, spacing , and the 5-point Laplacian stencil. Semi-discretizing in space turns the PDE into a system of ordinary differential equations — one per node:
where is the (positive definite) stiffness matrix.
# 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.
# 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
end
if i < n
rows[idx] = k
cols[idx] = k + 1
vals[idx] = -1.0
idx = idx + 1
end
if j > 1
rows[idx] = k
cols[idx] = k - n
vals[idx] = -1.0
idx = idx + 1
end
if j < n
rows[idx] = k
cols[idx] = k + n
vals[idx] = -1.0
idx = idx + 1
end
end
end
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]
end
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 :
# 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)
end
end
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: — for our grid about . Backward Euler instead solves
at each step. It costs a sparse solve per step, but is stable for any — this is why production codes use implicit integrators for diffusion.
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]
end
for k in 1:N
mr[used + k] = k
mc[used + k] = k
mv[used + k] = 1.0
end
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.
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)}")
end
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
# 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]
end
end
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 the blob diffuses away completely: with no heat source and cold boundaries, the steady state is . Because backward Euler is unconditionally stable, we can fast-forward with 10× larger steps:
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]
end
for k in 1:N
mr2[used + k] = k
mc2[used + k] = k
mv2[used + k] = 1.0
end
M2 = sparse(mr2, mc2, mv2, N, N)
for s in 1:5
(u, flag, relres, iters) = pcg(M2, u, 1e-10, 500)
end
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 —
pcgsolves 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.