6.4.1 Diffusion and the heat equation

The heat equation

Picture a thin metal rod along the $x$ -axis, wrapped so that heat can only flow along it. Write $u(x, t)$ for its temperature at position $x$ and time $t$ . The heat equation says how that temperature changes.

Rule: The heat equation

\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}, \quad \text{or} \quad u_t = u_{xx}

The left side is how fast the temperature at a fixed point changes in time. The right side is how the profile bends at that point at that moment: the second partial in $x$ , with $t$ held fixed (lesson 4.1.2).

Lesson 6.3.1 met this equation and checked solutions such as $x^2 + 2t$ and $e^{-t}\sin x$ by computing their partials. This lesson explains what the equation does to a profile, and why.

Units, and the diffusivity

For a real rod the equation reads $u_t = D\,u_{xx}$ . The diffusivity $D$ says how fast heat spreads in the material. It is in units of length² per unit of time, such as m² per second, so both sides are in degrees per unit of time.

A change of clock removes $D$ . Measure time as $s = Dt$ ; by the chain rule (lesson 4.2.2), $u_t = D\,u_s$ , so the equation becomes $u_s = u_{xx}$ . With $D = 0.5$ , two seconds make one unit of $s$ .

From here on $D = 1$ . Case study 6.4.4 removes a $D$ in exactly this way.

Hot to cold: the bend and the neighbours

The second central difference (lessons 2.4.2 and 2.8.1) estimates the bend from three values a step $h$ apart. Write $m = \frac{u(x + h) + u(x - h)}{2}$ for the average of the two neighbours. Then

\begin{aligned} u_{xx} &\approx \frac{u(x + h) - 2u(x) + u(x - h)}{h^2} \\ &= \frac{2}{h^2}\big(m - u(x)\big) \end{aligned}

So $u_{xx} > 0$ means the point is colder than its neighbours on average. Heat flows in from them, and the heat equation says the point warms. A point hotter than its neighbours cools.

Rule: Which way a point moves

\begin{gathered} u_{xx} > 0: \ u \text{ rises (a dip fills in)} \\[4pt] u_{xx} < 0: \ u \text{ falls (a peak sinks)} \\[4pt] u_{xx} = 0: \ u \text{ is still, for now} \end{gathered}

The speed is the size of the bend. Only the bend counts, not the height or the slope: a straight profile such as $3x + 1$ never changes, however steep it is.

A worked example: a dip between two peaks

Take a rod of length $\pi$ whose ends are held at temperature 0, a boundary condition as in lesson 6.3.2. Start it at $u(x, 0) = \sin x + 0.5\sin 3x$ : two peaks with a dip between them. The solution is

u(x, t) = e^{-t}\sin x + 0.5\,e^{-9t}\sin 3x

Each term solves the equation, and the equation is linear, so their sum does too (lesson 6.3.1). For the second term, $u_t = -4.5\,e^{-9t}\sin 3x$ , and differentiating twice in $x$ gives $0.5\,e^{-9t}(-9\sin 3x)$ , the same.

At $t = 0$ the bend is $u_{xx} = -\sin x - 4.5\sin 3x$ . The table checks it against the neighbours with $h = 0.1$ . It looks at the dip, the left peak, and the point between them where the bend changes sign.

point	$x$	$u$	neighbours' average $m$	second difference	$u_{xx}$
dip	1.5708	0.5	0.5173	3.4672	3.5
left peak	0.7017	1.0758	1.0534	$-4.4888$	$-4.5185$
between	1.1141	0.7978	0.7978	$-0.0063$	0

The dip is colder than its neighbours, so it warms at 3.5 degrees per unit of time. The peaks cool at about 4.5. Between them the point is, for a moment, still.

time $t$	0	0.05	0.1	0.2	0.5	1
dip, $x = \frac{\pi}{2}$	0.5	0.6324	0.7016	0.7361	0.601	0.3678
peak, $x = 0.7017$	1.0758	0.8884	0.759	0.5996	0.3963	0.2375

The dip fills in until $t = 0.188$ , when it reaches 0.7365. After that the whole profile sinks, because heat leaks out through the cold ends.

The profile bends down at the peaks, so they sink; it bends up at the dip, so the dip rises. By $t = 0.2$ the dip has gone, and the single hump that is left sinks as heat leaves through the ends.

u_t - u_{xx}

import math import sympy as sp x, t = sp.symbols("x t") expr = sp.exp(-t) * sp.sin(x) + sp.Rational(1, 2) * sp.exp(-9 * t) * sp.sin(3 * x) print(sp.simplify(sp.diff(expr, t) - sp.diff(expr, x, 2))) def u(x, t): return math.exp(-t) * math.sin(x) + 0.5 * math.exp(-9 * t) * math.sin(3 * x) h = 0.1 for name, p in [("dip", math.pi / 2), ("peak", 0.7017), ("between", 1.1141)]: left, mid, right = u(p - h, 0), u(p, 0), u(p + h, 0) print(name, round(mid, 4), round((left + right) / 2, 4), round((right - 2 * mid + left) / h**2, 4)) for when in [0, 0.05, 0.1, 0.2, 0.5, 1]: print(when, round(u(math.pi / 2, when), 4), round(u(0.7017, when), 4))

⌘/Ctrl + Enter runs the cell

Fine wiggles fade fastest

The example has two waves, and they fade at different rates. For any number $k$ , the wave $e^{-k^2 t}\sin kx$ solves the heat equation, because $u_t$ and $u_{xx}$ are both $-k^2 u$ .

Rule: Decaying waves

u(x, t) = e^{-k^2 t}\sin kx

A wave with $k$ times as many wiggles fades $k^2$ times as fast: its height halves after time $\frac{\ln 2}{k^2}$ . The same holds with $\cos kx$ . Fine detail dies first, so the heat equation smooths a profile.

wiggles $k$	1	2	3	10
height left at $t = 0.1$	0.9048	0.6703	0.4066	$4.54 \times 10^{-5}$
time to halve	0.6931	0.1733	0.077	0.0069

In the example the $\sin 3x$ wave fades 9 times as fast as $\sin x$ . By $t = 0.5$ it is down to $e^{-4.5} \approx 0.0111$ of its start, which is why one hump is left.

Lesson 6.3.2 found that running the heat equation backwards makes small ripples explode. The waves show why: backwards in time each one grows like $e^{k^2 t}$ , and the finest grow fastest.

Corners and jumps

A corner is smoothed at once. Take the kink $\max(x, 0)$ : its slope jumps from 0 to 1 at $x = 0$ , so $u_{xx}$ does not exist there. Its second difference at 0 is $\frac{h - 0 + 0}{h^2} = \frac{1}{h}$ , which is 10, 100 and 1000 for $h$ = 0.1, 0.01 and 0.001.

That grows without limit as $h$ shrinks, so the corner starts to rise infinitely fast and is rounded off immediately. A jump, such as a step from 0 to 1, is smoothed at once in the same way. Away from the corner, where the kink is straight, nothing happens at first.

t = 0.1

import math for k in [1, 2, 3, 10]: print(k, f"{math.exp(-k**2 * 0.1):.4g}", round(math.log(2) / k**2, 4)) print(round(math.exp(-9 * 0.5), 4)) def kink(x): return max(x, 0) for h in [0.1, 0.01, 0.001]: print(h, round((kink(h) - 2 * kink(0) + kink(-h)) / h**2, 6))

⌘/Ctrl + Enter runs the cell

The explorer shows the heat equation smoothing three profiles on a long rod. Choose the kink and set $t = 0.01$ : the corner is already round, while at $x = 0.5$ the profile has not moved. Then choose the spike and move $t$ up to 1. The height at 0 falls from 3.3333 to about 0.28, while the spread $\sqrt{2t}$ grows to 1.4142.

Explore

Heat flows from hot to cold, so peaks sink and dips fill in. Move time forward and watch the profile spread out; a corner is smoothed at once.

Dashed: the profile at time 0. Dark: the exact solution.

Time t = 0.05

Height at x = 0: 0.1262
Spread: the bell curve that smooths it has standard deviation √(2t) = 0.3162

Total heat stays the same

The total heat in the rod between $a$ and $b$ is the area under the profile, $\int_a^b u\,dx$ . Its rate of change is $\int_a^b u_t\,dx = \int_a^b u_{xx}\,dx$ . By the fundamental theorem (lesson 2.5.3) that is the difference of the slopes $u_x$ at the two ends.

Rule: Conservation of heat

\frac{d}{dt}\int_a^b u\,dx = u_x(b, t) - u_x(a, t)

The slope at an end measures heat flowing in or out there. With no flow through the ends, the total stays the same: heat moves about, but none is made or lost.

That holds for an insulated rod, whose ends have $u_x = 0$ , and on an endless line when the profile is flat and zero far out. It fails when heat can cross the ends.

The rod of the example has its ends held at 0, so heat leaks out. At the start $u_x = 2.5$ at $x = 0$ and $-2.5$ at $x = \pi$ , so the total falls at $-2.5 - 2.5 = -5$ per unit of time. Integrating the solution gives $2e^{-t} + \frac{1}{3}e^{-9t}$ : 2.3333 at the start, 1.9452 at $t = 0.1$ , 1.2168 at 0.5 and 0.7358 at 1.

2e^{-t} + \frac{1}{3}e^{-9t}

import math from scipy.integrate import quad def u(x, t): return math.exp(-t) * math.sin(x) + 0.5 * math.exp(-9 * t) * math.sin(3 * x) def ux(x, t): return math.exp(-t) * math.cos(x) + 1.5 * math.exp(-9 * t) * math.cos(3 * x) for t in [0, 0.1, 0.5, 1]: total = quad(lambda x: u(x, t), 0, math.pi)[0] print(t, round(total, 4), round(2 * math.exp(-t) + math.exp(-9 * t) / 3, 4)) print(round(ux(0, 0), 6), round(ux(math.pi, 0), 6), round(ux(math.pi, 0) - ux(0, 0), 6))

⌘/Ctrl + Enter runs the cell

Spreading like $\sqrt{t}$

Now put one unit of heat in a narrow tent around $x = 0$ on an endless rod. As time passes the heat spreads out, but its centre stays at 0 and its total stays 1. Read the profile as a density (lesson 5.4.1): its variance, $\int x^2 u\,dx$ , grows by exactly 2 per unit of time.

Rule: Spreading

\begin{gathered} \operatorname{Var}(t) = \operatorname{Var}(0) + 2t \\[4pt] \text{width of heat from a point: } \sqrt{2t} \end{gathered}

For a total heat of 1 centred at 0; the width is the standard deviation. Stated here and checked below; lesson 6.4.2 meets it in an exact solution. Since the area stays 1 while the width grows, the peak falls like $\frac{1}{\sqrt{t}}$ .

time $t$	0.01	0.04	0.16	0.64
width $\sqrt{2t}$	0.1414	0.2828	0.5657	1.1314

Four times the time gives only twice the width. A fair random walk spreads in the same way: after $n$ steps its standard deviation is $\sqrt{n}$ (lesson 5.11.1). Heat is the same fact in continuous time, and lesson 6.4.2 makes the link exact.

Stepping a profile forward

A computer can follow the heat equation directly. Keep the profile's values on a grid of points $\Delta x$ apart, and estimate the bend at each by the second difference. Then take a short time step $\Delta t$ : each value moves by $\Delta t$ times its bend, towards its neighbours' average.

Lesson 6.5.1 studies this method, its error and how short $\Delta t$ must be. Here $\Delta x = 0.05$ and $\Delta t = 0.0005$ , a fifth of $\Delta x^2$ , which is safely small.

-8

import math import numpy as np dx, dt = 0.05, 0.0005 x = np.linspace(-8, 8, 321) u = np.maximum(0, 1 - np.abs(x) / 0.2) / 0.2 var0 = np.sum(x**2 * u) * dx t = 0 for stop in [0.25, 0.5, 1.0]: while t < stop - 1e-9: bend = (u[2:] - 2 * u[1:-1] + u[:-2]) / dx**2 u[1:-1] = u[1:-1] + dt * bend t += dt total = np.sum(u) * dx var = np.sum(x**2 * u) * dx / total print(stop, round(total, 6), round(var - var0, 4), round(u.max(), 4), round(1 / math.sqrt(4 * math.pi * stop), 4))

⌘/Ctrl + Enter runs the cell

Common mistakes

The bend drives the change, not the height or the slope. A steep but straight profile does not change at all.
Check the sign. $u$ rises where $u_{xx} > 0$ , where the profile bends up: a dip. It falls at a peak.
Corners and jumps have no $u_{xx}$ . The second difference there grows like $\frac{1}{h}$ , and the corner is rounded at once.
Heat is conserved only when none crosses the ends. Ends held at a fixed temperature let heat out or in, at the rate given by the slopes there.
The width grows like $\sqrt{t}$ , not like $t$ . From a single point it is $\sqrt{2t}$ : four times the time, twice the width.
Mind the clock. For $u_t = D\,u_{xx}$ , use the results of this lesson with $s = Dt$ in place of $t$ .

Diffusion and the heat equation

0 of 13 correct

1.
For which number $c$ does $u = 5x^2 + ct$ solve $u_t = u_{xx}$ ?
2.
A heat equation reads $u_t = 0.25\,u_{xx}$ . Measuring time as $s = 0.25t$ turns it into $u_s = u_{xx}$ . What is $s$ when $t = 8$ ?
3.
Three neighbouring temperatures, 0.1 apart, are 1.2, 1.0 and 1.0. Estimate $u_{xx}$ at the middle point by the second difference
4.
A point is hotter than the average of its two neighbours. Type 1 if it warms, $-1$ if it cools, 0 if it stays the same
5.
A rod's profile is $u(x, 0) = 4 - 3x$ . Find $u_t$ at $x = 1$ , $t = 0$
6.
At $t = 0$ , $u = \sin x + 0.5\sin 3x$ . Find $u_t$ at $x = \frac{\pi}{2}$
7.
What fraction of the wave $e^{-9t}\sin 3x$ is left at $t = 0.1$ ? Give an exact answer
8.
How long does the wave $e^{-4t}\sin 2x$ take to halve? Give an exact answer
9.
Find the second difference of $\max(x, 0)$ at $x = 0$ with $h = 0.005$
10.
An insulated rod holds a total heat of 2.5 at $t = 0$ . What is its total heat at $t = 3$ ?
11.
At one instant a rod on $[0, 1]$ has $u_x(0, t) = 1$ and $u_x(1, t) = -1$ . At what rate is its total heat changing?
12.
Heat of total 1 on an endless rod has variance 0.3 at $t = 0$ . What is its variance at $t = 0.6$ ?
13.
Heat from a single point has width $\sqrt{2t}$ . At what time is the width 1?

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

u_t = u_{xx}

Diffusion and the heat equation

The heat equation

Rule: The heat equation

Units, and the diffusivity

Hot to cold: the bend and the neighbours

Rule: Which way a point moves

A worked example: a dip between two peaks

Fine wiggles fade fastest

Rule: Decaying waves

Corners and jumps

Total heat stays the same

Rule: Conservation of heat

Spreading like $\sqrt{t}$

Rule: Spreading

Stepping a profile forward

Common mistakes

Diffusion and the heat equation

Summary card

Free further reading

The heat equation

Rule: The heat equation

Units, and the diffusivity

Hot to cold: the bend and the neighbours

Rule: Which way a point moves

A worked example: a dip between two peaks

Fine wiggles fade fastest

Rule: Decaying waves

Corners and jumps

Total heat stays the same

Rule: Conservation of heat

Spreading like t\sqrt{t}t​

Rule: Spreading

Stepping a profile forward

Common mistakes

Summary card

Free further reading

Spreading like $\sqrt{t}$