6.5.3 Practice set: finite differences

The explicit scheme

The explicit scheme and stability (lesson 6.5.1)

0 of 8 correct

1.
A rod of length 2 is split into 40 equal cells. Find $\Delta x$
2.
With $r = 0.5$ , a point has value 0.3 and neighbours 0.1 and 0.9. Find its value one step later
3.
Four cells, ends held at 0, start $(0, 0, 1, 0, 0)$ , $r = 0.25$ . Find the middle value after two steps
4.
The ends of a rod of length 1 are held at 1 and 5. What value does the point $x = 0.25$ settle on?
5.
For $u_t = 2u_{xx}$ with $\Delta x = 0.1$ , the explicit scheme has $r = \frac{2\Delta t}{\Delta x^2}$ . What is the largest stable $\Delta t$ ?
6.
At $r = 0.625$ , one step multiplies a zig-zag by $1 - 4r$ . A zig-zag of size 0.01 takes two steps. How big is it then?
7.
With $r$ fixed at 0.4, a grid of 50 cells has largest error $2 \times 10^{-4}$ . Estimate the largest error on 100 cells
8.
A stable grid has $r = 0.1$ . Someone triples the number of cells and keeps $\Delta t$ . Find the new $r$

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.

Implicit schemes and Crank–Nicolson

Implicit schemes and Crank–Nicolson (lesson 6.5.2)

0 of 8 correct

1.
Implicit scheme, $r = 0.5$ , one inside point between ends held at 0, old value 1. Find its new value
2.
Implicit scheme, $r = 1$ , one inside point between ends held at 1 and 3, old value 2. Find its new value
3.
One implicit step multiplies a zig-zag by $\frac{1}{1 + 4r}$ . For which $r$ does a step shrink it to a tenth?
4.
Halving both steps of the implicit scheme takes its error from 0.01 to about what?
5.
Crank–Nicolson has error $6 \times 10^{-4}$ on 20 cells with $\Delta t = \frac{\Delta x}{10}$ . Estimate its error on 40 cells, still with $\Delta t = \frac{\Delta x}{10}$
6.
Rannacher's start replaces the first Crank–Nicolson step, of length 0.02, by two implicit steps. How long is each?
7.
solve_banded((1, 1), ab, b) solves a 5 × 5 tridiagonal system. ab has 3 rows of 5 entries. How many of the 15 are unused?
8.
Crank–Nicolson, $r = 2$ , one inside point between ends held at 0, old value 1. Find its new value

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.

Coding exercise

Exercise

Three schemes in one function

Write heat_scheme(u0, r, steps, scheme="cn") for $u_t = u_{xx}$ on a grid, as in lessons 6.5.1 and 6.5.2. u0 holds the grid values at time 0, ends included; r is $\frac{\Delta t}{\Delta x^2}$ . Return the grid values after steps steps, as a NumPy array.

scheme is "explicit", "implicit" or "cn" (Crank–Nicolson). Hold the two end values fixed at every step, and include them in the boundary terms of the implicit schemes.
Solve the implicit systems with scipy.linalg.solve_banded. Never change u0 itself; steps = 0 returns a copy of it.
Raise ValueError if u0 has fewer than 3 values or if r is not positive. Raise it if steps is not a whole number of at least 0, or scheme is not one of the three. Raise it too for the explicit scheme with $r > \frac{1}{2}$ .

The tests use the spike and the sine of the lessons and a rod whose ends are held at 0 and 1. They add a straight line that must not move, and a tent with $r = 50$ .

See the tests

import math

import numpy as np


def expect_value_error(call, why):
    try:
        call()
    except ValueError:
        return
    raise AssertionError(f"expected ValueError: {why}")


def test_implicit_spike():
    """One implicit step with r = 1 spreads the spike (0, 0, 1, 0, 0) to 1/7, 3/7, 1/7 inside"""
    u = np.asarray(heat_scheme([0, 0, 1, 0, 0], 1, 1, "implicit"))
    expect_close(u, np.array([0, 1, 3, 1, 0]) / 7, what="heat_scheme(spike, 1, 1, 'implicit')", rel=1e-9, abs_tol=1e-12)


def test_cn_spike():
    """One Crank–Nicolson step with r = 1 gives 2/7, 1/7, 2/7 inside: the middle dips"""
    u = np.asarray(heat_scheme([0, 0, 1, 0, 0], 1, 1, "cn"))
    expect_close(u, np.array([0, 2, 1, 2, 0]) / 7, what="heat_scheme(spike, 1, 1, 'cn')", rel=1e-9, abs_tol=1e-12)


def test_explicit_sine():
    """One explicit step with r = 0.25 on sin(πx) with 4 cells multiplies it by 1 − sin²(π/8)"""
    x = np.linspace(0, 1, 5)
    u0 = np.sin(math.pi * x)
    u = np.asarray(heat_scheme(u0, 0.25, 1, "explicit"))
    expect_close(u, (1 - math.sin(math.pi / 8) ** 2) * u0, what="one explicit step", rel=1e-9, abs_tol=1e-12)


def test_accuracy():
    """On 40 cells to t = 0.1 with 40 steps (r = 4), Crank–Nicolson is within 2e-4 of e^(−π²t) sin(πx) and the implicit scheme within 5e-3"""
    x = np.linspace(0, 1, 41)
    exact = math.exp(-math.pi**2 * 0.1) * np.sin(math.pi * x)
    cn = np.max(np.abs(np.asarray(heat_scheme(np.sin(math.pi * x), 4, 40, "cn")) - exact))
    im = np.max(np.abs(np.asarray(heat_scheme(np.sin(math.pi * x), 4, 40, "implicit")) - exact))
    assert cn < 2e-4, f"Crank–Nicolson error should be under 2e-4, got {cn:.2e}"
    assert 1e-3 < im < 5e-3, f"implicit error should be between 1e-3 and 5e-3, got {im:.2e}"


def test_ends_held():
    """The ends keep their values: from 0 inside with ends 0 and 1, the profile settles on the line u = x"""
    u0 = np.zeros(11)
    u0[-1] = 1
    for scheme, r, steps in [("explicit", 0.4, 250), ("implicit", 5, 20), ("cn", 2, 50)]:
        u = np.asarray(heat_scheme(u0, r, steps, scheme))
        expect_close(u[[0, -1]], np.array([0.0, 1.0]), what=f"{scheme}: the end values", abs_tol=1e-12)
        expect_close(u, np.linspace(0, 1, 11), what=f"{scheme}: the profile at t = 1", abs_tol=1e-3)


def test_line_is_steady():
    """A straight line between the end values does not change, for every scheme"""
    line = np.linspace(2, 6, 9)
    for scheme in ["explicit", "implicit", "cn"]:
        expect_close(np.asarray(heat_scheme(line, 0.5, 10, scheme)), line, what=f"{scheme} on a line", abs_tol=1e-12)


def test_large_r_stays_bounded():
    """With r = 50 the implicit scheme keeps a tent between 0 and 0.5"""
    x = np.linspace(0, 1, 21)
    u = np.asarray(heat_scheme(np.minimum(x, 1 - x), 50, 5, "implicit"))
    assert u.min() >= 0 and u.max() <= 0.5, f"values should stay between 0 and 0.5, got {u.min():.3g} to {u.max():.3g}"


def test_input_untouched():
    """Zero steps return the start unchanged, and u0 itself is never modified"""
    u0 = np.array([0.0, 0.5, 1.0, 0.5, 0.0])
    keep = u0.copy()
    expect_close(np.asarray(heat_scheme(u0, 1, 0, "cn")), keep, what="heat_scheme(u0, 1, 0)", abs_tol=0)
    heat_scheme(u0, 1, 3, "implicit")
    expect_close(u0, keep, what="u0 after the call", abs_tol=0)


def test_bad_inputs():
    """An unstable explicit r, a non-positive r, a fractional step count, too few points or an unknown scheme raise ValueError"""
    u0 = [0, 0.5, 1, 0.5, 0]
    expect_value_error(lambda: heat_scheme(u0, 0.6, 5, "explicit"), "the explicit scheme is unstable for r = 0.6")
    expect_value_error(lambda: heat_scheme(u0, 0, 5, "cn"), "r = 0 is not a time step")
    expect_value_error(lambda: heat_scheme(u0, -1, 5, "implicit"), "r must be positive")
    expect_value_error(lambda: heat_scheme(u0, 1, 2.5, "cn"), "2.5 steps is not a whole number")
    expect_value_error(lambda: heat_scheme([0, 1], 1, 1, "cn"), "two values leave no inside point")
    expect_value_error(lambda: heat_scheme(u0, 1, 1, "euler"), "'euler' is not one of the three schemes")

The explicit scheme

Implicit schemes and Crank–Nicolson

Coding exercise

Free further reading