6.5.1 The explicit scheme and stability

A problem with a known answer

The heat equation of lesson 6.4.1 is $u_t = u_{xx}$ . Here $u(x, t)$ is the temperature at position $x$ and time $t$ . $u_t$ is its rate of change in time, and $u_{xx}$ its second partial derivative in $x$ .

Take a rod from $x = 0$ to $x = 1$ whose ends are held at temperature 0, and start it with the profile $\sin(\pi x)$ . These are one initial condition and two boundary conditions, as in lesson 6.3.2. The solution is known:

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

It is lesson 6.4.1's decaying wave $e^{-k^2t}\sin kx$ with $k = \pi$ . Check it by computing partials, as lesson 6.3.1 does: differentiating in $t$ gives $u_t = -\pi^2 e^{-\pi^2 t}\sin(\pi x)$ , and differentiating twice in $x$ gives the same. The profile keeps its shape and shrinks: at $t = 0.1$ the middle value is $e^{-0.1\pi^2} \approx 0.372708$ .

Most problems have no formula like this. A finite-difference scheme finds the solution approximately, from numbers alone, and this exact answer lets you measure how well it does.

The grid

Split the rod into $M$ equal cells of width $\Delta x = \frac{1}{M}$ , and time into steps of length $\Delta t$ . The grid points are $x_i = i\,\Delta x$ for $i = 0, 1, \ldots, M$ , and the times are $t_n = n\,\Delta t$ for $n = 0, 1, 2, \ldots$

Rule: A grid for u(x, t)

\begin{gathered} x_i = i\,\Delta x, \quad t_n = n\,\Delta t \\[4pt] u_i^n \approx u(x_i, t_n) \end{gathered}

$u_i^n$ is the scheme's value at point $i$ and time step $n$ . The upper index $n$ counts time steps; it is not a power. The bottom row $u_i^0$ is the initial condition, and the two end columns $u_0^n$ and $u_M^n$ are the boundary values.

The grid for $u_t = u_{xx}$ on a rod: the initial condition fills the bottom row and the boundary conditions fill the end columns. One step of the explicit scheme takes three neighbouring values at one time and gives the middle value one step later.

The explicit update

Replace each derivative in $u_t = u_{xx}$ by a difference of grid values, both at time $t_n$ .

Time: the forward difference of lesson 2.2.2, $u_t \approx \frac{u_i^{n+1} - u_i^n}{\Delta t}$ . Lesson 2.8.1 gave its error: about $\frac{1}{2}u_{tt}\,\Delta t$ .
Space: the second difference of lesson 2.4.2, $u_{xx} \approx \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}$ , with error about $\frac{1}{12}u_{xxxx}\,\Delta x^2$ (lesson 2.8.1).

Set the two equal and solve for the one value at the new time:

Rule: The explicit scheme

\begin{gathered} u_i^{n+1} = u_i^n + r\left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right) \\[6pt] r = \frac{\Delta t}{\Delta x^2} \end{gathered}

Each new value is the old one plus $r$ times its second difference. So it rises where the profile bends up and falls where it bends down, as heat does (lesson 6.4.1). The scheme is explicit: each new value comes straight from old ones.

This is the step lesson 6.4.1's cell took, with $\Delta x = 0.05$ and $\Delta t = 0.0005$ , so $r = 0.2$ . This lesson asks how large $r$ may be, and how the error depends on the grid.

$r$ is a pure number here, because $u_t = u_{xx}$ measures time in units of distance squared. For an equation $u_t = a\,u_{xx}$ with a constant $a > 0$ , the same steps give $r = \frac{a\,\Delta t}{\Delta x^2}$ .

One step by hand

Take $M = 4$ cells, so $\Delta x = 0.25$ , and $r = 0.25$ , so $\Delta t = r\,\Delta x^2 = 0.015625$ . Start from $\sin(\pi x)$ on the grid and apply the update at the three inside points:

$x_i$	0	0.25	0.5	0.75	1
start, $u_i^0$	0	0.7071	1	0.7071	0
second difference		$-0.4142$	$-0.5858$	$-0.4142$
one step, $u_i^1$	0	0.6036	0.8536	0.6036	0
exact at $t = 0.015625$	0	0.6061	0.8571	0.6061	0

At the middle, $u_2^1 = 1 + 0.25 \times (-0.5858) = 0.8536$ . The ends stay at 0: they are boundary values, not updated.

Every inside value shrank by the same factor, 0.8536, so the scheme keeps the sine's shape. The exact solution shrinks by $e^{-\pi^2\Delta t} \approx 0.8571$ in the same time, so the scheme cools the rod slightly too fast.

That common factor is no accident. The second difference of a sine on the grid is a fixed multiple of the sine, so each step multiplies it by $g = 1 - 4r\sin^2(\pi\Delta x/2)$ . Here that is $1 - \sin^2(\pi/8) \approx 0.8536$ .

g

import math import numpy as np M, r = 4, 0.25 dx = 1 / M dt = r * dx**2 x = np.linspace(0, 1, M + 1) u = np.sin(math.pi * x) second = u[2:] - 2 * u[1:-1] + u[:-2] new = u.copy() new[1:-1] = u[1:-1] + r * second print(dt) print(np.round(u, 4), np.round(second, 4)) print(np.round(new, 4)) print(np.round(math.exp(-math.pi**2 * dt) * np.sin(math.pi * x), 4)) print(np.round(new[1:-1] / u[1:-1], 4), round(1 - 4 * r * math.sin(math.pi * dx / 2)**2, 4))

⌘/Ctrl + Enter runs the cell

Boundary values

The update needs a neighbour on each side, so it runs only at the inside points $i = 1, \ldots, M - 1$ . At every step the two ends are set by the boundary conditions instead. In code, updating u[1:-1] leaves u[0] and u[-1] alone.

Hold the left end at 0 and the right end at 1, and start with 0 everywhere inside. Take $M = 10$ and $r = 0.4$ , so $\Delta t = 0.004$ .

After one step only the point next to the hot end has changed: $u_9 = 0 + 0.4 \times (1 - 0 + 0) = 0.4$ .
After two steps, $u_9 = 0.48$ and $u_8 = 0.16$ . Heat enters one cell per step.
After 25 steps ( $t = 0.1$ ) the middle value is 0.2674; after 125 steps ( $t = 0.5$ ) it is 0.4957.

The profile settles on the straight line $u = x$ between the two end values. There $u_{xx} = 0$ , so nothing changes any more: a steady state. At $t = 1$ the scheme is within $2.91 \times 10^{-5}$ of the line everywhere.

u = x

import numpy as np def explicit(u0, r, steps): u = u0.copy() for _ in range(steps): u[1:-1] = u[1:-1] + r * (u[2:] - 2 * u[1:-1] + u[:-2]) return u M, r = 10, 0.4 x = np.linspace(0, 1, M + 1) u0 = np.zeros(M + 1) u0[-1] = 1 # right end held at 1, left end at 0 print(np.round(explicit(u0, r, 1), 4)) for t in [0.01, 0.1, 0.5, 1]: steps = round(t / (r / M**2)) u = explicit(u0, r, steps) print(t, steps, np.round(u, 4), f"{np.max(np.abs(u - x)):.2e}")

⌘/Ctrl + Enter runs the cell

Stability: r at most one half

A larger $r$ means a longer time step and fewer steps to reach a given time. But $r$ cannot be as large as you like. Gather the terms of the update by grid value:

u_i^{n+1} = r\,u_{i-1}^n + (1 - 2r)\,u_i^n + r\,u_{i+1}^n

The three weights add up to 1. When $r \le \frac{1}{2}$ , none is negative, so each new value is an average of old ones and lies between them: nothing can grow. When $r > \frac{1}{2}$ , the middle weight is negative.

The zig-zag

A zig-zag is a profile whose grid values alternate: $+1, -1, +1, -1$ , and so on. At a point with value $-1$ , the second difference is $1 - 2(-1) + 1 = 4$ ; at a point with value $+1$ it is $-4$ . So it is $-4$ times the value everywhere, and one step multiplies the whole zig-zag by $1 - 4r$ .

$r$	0.25	0.4	0.5	0.55	0.6
factor $1 - 4r$	0	$-0.6$	$-1$	$-1.2$	$-1.4$

For $r < \frac{1}{2}$ the factor lies between $-1$ and 1, so the zig-zag dies away. At $r = \frac{1}{2}$ it flips sign each step but keeps its size. Above $\frac{1}{2}$ it flips sign and grows, by 1.4 times per step at $r = 0.6$ .

A corner in the profile contains a little of every wave the grid can hold, the zig-zag included, and so does rounding error. So above $\frac{1}{2}$ there is always a zig-zag to grow.

Rule: Stability of the explicit scheme

\begin{gathered} r = \frac{\Delta t}{\Delta x^2} \le \frac{1}{2} \\[6pt] \text{that is, } \Delta t \le \tfrac{1}{2}\Delta x^2 \end{gathered}

A scheme is stable if errors do not grow from step to step. The zig-zag is the fastest wave a grid can hold and the first to grow, so this limit on $r$ is the whole condition. It is stated here, not proved.

Start from a tent with a corner at the middle, $\min(x, 1 - x)$ , on $M = 20$ cells. Here is the largest value $|u_i|$ after each number of steps:

steps	10	20	30	40	50
$r = 0.4$	0.3875	0.3407	0.3048	0.2745	0.2481
$r = 0.5$	0.377	0.3238	0.2833	0.2494	0.2201
$r = 0.55$	0.3879	0.4022	0.7223	2.641	13.31
$r = 0.6$	0.4605	2.325	47.93	1189	30,370

The tent starts at 0.5 and should only cool. For $r \le \frac{1}{2}$ it does; above, the values explode. After 20 steps at $r = 0.6$ , the seven values around the corner are $-1.2909$ , 2.0903, $-1.6636$ , 2.3253, $-1.6636$ , 2.0903 and $-1.2909$ : a zig-zag.

r = 0.6

import numpy as np def explicit(u0, r, steps): u = u0.copy() for _ in range(steps): u[1:-1] = u[1:-1] + r * (u[2:] - 2 * u[1:-1] + u[:-2]) return u zig = np.array([1.0, -1, 1, -1, 1, -1, 1]) print(zig[2:] - 2 * zig[1:-1] + zig[:-2]) print(explicit(zig, 0.6, 1)) x = np.linspace(0, 1, 21) tent = np.minimum(x, 1 - x) for r in [0.4, 0.5, 0.55, 0.6]: print(r, [f"{np.max(np.abs(explicit(tent, r, s))):.4g}" for s in [10, 20, 30, 40, 50]]) print(np.round(explicit(tent, 0.6, 20)[7:14], 4))

⌘/Ctrl + Enter runs the cell

The explorer runs the scheme on a kink, $\max(x, 0)$ , with its ends held, beside the exact solution. Set the time to 0.5 and pick the $r$ presets in turn. At 0.25 and 0.5 the dots sit on the curve; at 0.55 a zig-zag appears at the corner.

At 0.6 the zig-zag swamps everything. Finally, at $r = 0.55$ , move the cells from 30 to 60 and watch it blow up sooner.

Explore

The scheme updates each grid value from its two neighbours, with r = Δt/Δx². Keep r at or below 1/2 and it follows the exact curve; push r past 1/2 and the error grows from step to step until it swamps the answer.

Dashed: the profile at time 0. Dark: the exact solution. Dots: the scheme on its grid.

Time t = 0.1

r = Δt/Δx² = 0.4

Grid cells = 30

Grid: Δx = 0.2, Δt = r Δx² = 0.016, 6 steps to t = 0.096
Largest error: 0.002

The price of stability

The limit ties the time step to the square of the space step. Halve $\Delta x$ and the largest stable $\Delta t$ falls 4 times. Reaching the same time then takes 4 times as many steps, each with twice as many points.

cells $M$	10	100	1,000
largest $\Delta t$	0.005	$5 \times 10^{-5}$	$5 \times 10^{-7}$
steps to $t = 0.1$	20	2,000	200,000

Lesson 6.5.2 removes this limit with schemes that are stable for every $r$ .

The error as the grid gets finer

Back to the sine, whose exact solution is known. Keep $r = 0.4$ , run to $t = 0.1$ , and measure the largest error over the grid:

cells $M$	$\Delta x$	$\Delta t$	steps	largest error	error ratio
10	0.1	0.004	25	$4.294 \times 10^{-3}$
20	0.05	0.001	100	$1.063 \times 10^{-3}$	4.04
40	0.025	0.00025	400	$2.649 \times 10^{-4}$	4.01
80	0.0125	0.0000625	1,600	$6.620 \times 10^{-5}$	4.00

Each halving of $\Delta x$ divides the error by 4, so the error shrinks like $\Delta x^2$ . The two parts of lesson 2.8.1's error explain it. The time error, about $\frac{1}{2}u_{tt}\,\Delta t$ , equals $\frac{r}{2}u_{tt}\,\Delta x^2$ once $\Delta t = r\,\Delta x^2$ ; the space error is about $\frac{1}{12}u_{xxxx}\,\Delta x^2$ .

Rule: Error of the explicit scheme

\begin{gathered} \text{error} \approx C_1\,\Delta t + C_2\,\Delta x^2 \\[4pt] \text{with } r \text{ fixed: } \propto \Delta x^2 \end{gathered}

$C_1$ and $C_2$ depend on the solution, not on the grid. The scheme is first order in time and second order in space. With $r$ fixed, halving $\Delta x$ cuts the error about 4 times and costs about 8 times the work.

Stability and accuracy are different questions. Stability says the errors do not grow; the step sizes say how large they are to begin with.

t = 0.1

import math import numpy as np def explicit(u0, r, steps): u = u0.copy() for _ in range(steps): u[1:-1] = u[1:-1] + r * (u[2:] - 2 * u[1:-1] + u[:-2]) return u r, t, previous = 0.4, 0.1, None for M in [10, 20, 40, 80]: x = np.linspace(0, 1, M + 1) dt = r / M**2 steps = round(t / dt) u = explicit(np.sin(math.pi * x), r, steps) error = np.max(np.abs(u - math.exp(-math.pi**2 * t) * np.sin(math.pi * x))) print(M, steps, f"{error:.3e}", round(previous / error, 2) if previous else "") previous = error print(round(math.exp(-math.pi**2 * t), 6))

⌘/Ctrl + Enter runs the cell

Case study 6.5.4 uses these ideas on the Black–Scholes equation, with the improved schemes of lesson 6.5.2.

Common mistakes

Halving $\Delta x$ but keeping $\Delta t$ . That multiplies $r$ by 4: a stable $r = 0.4$ becomes $1.6$ , and the answer explodes. Cut $\Delta t$ by 4 too.
Using $\Delta x$ instead of $\Delta x^2$ in $r$ . With $\Delta x = 0.1$ and $\Delta t = 0.004$ , $r$ is 0.4, not 0.04.
Updating in place in a loop. Computing u[i] from a u[i - 1] already changed this step gives a different scheme. Keep the old values until every new one is made; the slice line does this.
Updating the ends. The end values come from the boundary conditions, not from the update.
Reading $u_i^n$ as a power. The upper index $n$ counts time steps.
Taking stable for accurate. At $r \le \frac{1}{2}$ errors do not grow, but a coarse grid still has large ones.

The explicit scheme

0 of 12 correct

1.
A grid has $\Delta x = 0.1$ and $\Delta t = 0.002$ . Find $r = \frac{\Delta t}{\Delta x^2}$
2.
With $r = 0.25$ , a point has value 0.6 and neighbours 0.5 and 0.9. Find its value one step later
3.
With $r = 0.4$ , a point has value 1 and both neighbours 0. Find its value one step later
4.
With $r = 0.3$ , a point has value 0 and neighbours 0 and 1. Find its value one step later
5.
An explicit scheme has $\Delta x = 0.05$ . What is the largest stable time step $\Delta t$ ?
6.
One step multiplies a zig-zag by $1 - 4r$ . Find that factor for $r = 0.7$
7.
At $r = 0.75$ , one step multiplies a zig-zag by $-2$ . Starting from size $10^{-6}$ , how big is it after 10 steps?
8.
A stable grid has $\Delta x = 0.1$ and $\Delta t = 0.004$ . Someone halves $\Delta x$ and keeps $\Delta t$ . Find the new $r$
9.
Run to $t = 0.1$ with $r = 0.5$ on $M = 50$ cells of a rod of length 1. How many steps does it take?
10.
With $r$ fixed, the largest error is $8 \times 10^{-4}$ on a grid. About what is it after halving $\Delta x$ ?
11.
With $r$ fixed, halving $\Delta x$ doubles the points and multiplies the steps by how much?
12.
The ends of a rod are held at 2 and 6. What value does the middle of the rod settle on?

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.

x_i = i\,\Delta x

The explicit scheme and stability

A problem with a known answer

The grid

Rule: A grid for u(x, t)

The explicit update

Rule: The explicit scheme

One step by hand

Boundary values

Stability: r at most one half

The zig-zag

Rule: Stability of the explicit scheme

The price of stability

The error as the grid gets finer

Rule: Error of the explicit scheme

Common mistakes

The explicit scheme

Summary card

Free further reading