Notebooks · posts/make-the-heat-equation-step-stable

Make the heat-equation step stable past r = 1/2

Adam Zettel · September 25, 2026 · written with an agent

The explicit (FTCS) step for $u_{t} = u_{xx}$ multiplies the grid's most jagged mode by about $1 - 4 r$ each step, and for $r = 0.6$ that is $- 1.4$ : a checkerboard that should fade grows to 215 in 16 steps. Rewrite heat_step as backward Euler, which needs a tridiagonal solve per step and damps every mode for every $r$ .

u_{i}^{n + 1} = u_{i}^{n} + r (u_{i - 1}^{n} - 2 u_{i}^{n} + u_{i + 1}^{n}), r = \frac{Δ t}{Δ x ^{2}} = 0.6 \Rightarrow blows up

FTCS stands for forward in time, centered in space: the new value at each point is the old one plus $r$ times the second difference of its neighbors. It is stable only for $r \leq 1/2$ , so halving $Δ x$ forces a time step four times smaller. Backward Euler takes the second difference at the new time instead,

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

a tridiagonal linear system for the new values, with $- r$ , $1 + 2 r$ , $- r$ on its three diagonals. Viv's tridiag_solve(lower, diagonal, upper, rhs) solves one.

The grid has 19 interior points with $u = 0$ at both ends, and the arrays include the two ends. On this grid $sin (mπ i /20)$ is a mode of both methods: backward Euler multiplies it by exactly $1/ (1 + 4 r sin^{2} (mπ /40))$ per step, which the first two tests use. The other expected numbers come from posts/refs/make-the-heat-equation-step-stable.py, which runs backward Euler with numpy.linalg.solve on the full matrix.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

This is a problem: its tests fail until the functions in it are written. Edit a cell and run it to try; the cells below it show whether the tests pass.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 32–125

# interior points; the arrays also hold the two ends, where u stays 0
N_IN = 19
NSTEP = 16
R_MAIN = 0.6
PI = 3.141592653589793

#| One time step of $u_t = u_{xx}$ with $r = \Delta t / \Delta x^2$ on an array whose
#| two ends are held at 0. As written it is FTCS. Make it backward Euler.
def heat_step(u: Arr, r) {
    n = len(u)
    v = zeros(n)
    for i from 1 to n - 1 {
        v[i] = u[i] + r * (u[i - 1] - 2.0 * u[i] + u[i + 1])
    }
    return v
}

#| The mode $\sin(m\pi i/(N+1))$ on the grid, ends included.
def mode(m) {
    u = zeros(N_IN + 2)
    for i to N_IN + 2 {
        u[i] = sin(m * PI * i / (N_IN + 1.0))
    }
    return u
}

#| The checkerboard +1, -1, +1, ... on the interior, 0 at the ends.
def checkerboard() {
    u = zeros(N_IN + 2)
    s = 1.0
    for i from 1 to N_IN + 1 {
        u[i] = s
        s = -s
    }
    return u
}

#| u after NSTEP calls of heat_step with ratio r.
def run(u0: Arr, r) {
    u = u0 + 0.0
    for k to NSTEP {
        u = heat_step(u, r)
    }
    return u
}

#| The largest |u[i]|.
def max_abs(u: Arr) {
    m = 0.0
    for i to len(u) {
        m = max(m, abs(u[i]))
    }
    return m
}

show max_abs(run(checkerboard(), R_MAIN))
# FTCS multiplies the most jagged mode by about 1 - 4r per step; it is stable
# while that stays at or above -1: 1 - 4r = -1 at r = (1 - (-1)) / 4
ftcs_jagged_factor = 1.0 - 4.0 * R_MAIN
ftcs_r_limit = (1.0 - (-1.0)) / 4.0
four_r = 4.0 * R_MAIN
show ftcs_jagged_factor
show ftcs_r_limit
show four_r

test "the smoothest mode decays by 1 / (1 + 2.4 sin^2(pi/40))^16 = 0.790844 at r = 0.6" {
    g = 1.0 / (1.0 + 4.0 * R_MAIN * pow(sin(PI / 40.0), 2.0))
    u = run(mode(1.0), R_MAIN)
    assert_near(u[5] / mode(1.0)[5], pow(g, NSTEP), 1.0e-12)
    assert_near(pow(g, NSTEP), 0.7908438775198426, 1.0e-13)
}
test "the most jagged mode, m = 19, decays to 3.362e-9 of itself instead of growing" {
    g = 1.0 / (1.0 + 4.0 * R_MAIN * pow(sin(19.0 * PI / 40.0), 2.0))
    u = run(mode(19.0), R_MAIN)
    assert_near(u[5] / mode(19.0)[5], pow(g, NSTEP), 1.0e-15)
    assert_near(pow(g, NSTEP), 3.361996814033325e-09, 1.0e-20)
}
test "the checkerboard at r = 0.6 fades to numpy's 0.0068342 after 16 steps" {
    u = run(checkerboard(), R_MAIN)
    assert_near(max_abs(u), 0.006834176520732032, 1.0e-14)
    assert_near(u[10], 0.00318114592290907, 1.0e-14)
}
test "and at r = 5, ten times the FTCS limit, to numpy's 0.0012278" {
    u = run(checkerboard(), 5.0)
    assert_near(max_abs(u), 0.0012277856452884454, 1.0e-14)
}
test "the ends stay at 0" {
    u = run(checkerboard(), R_MAIN)
    assert_near(u[0], 0.0, 0.0)
    assert_near(u[N_IN + 1], 0.0, 0.0)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

# interior points; the arrays also hold the two ends, where u stays 0
N_IN = 19
NSTEP = 16
R_MAIN = 0.6
PI = 3.141592653589793

#| One time step of $u_t = u_{xx}$ with $r = \Delta t / \Delta x^2$ on an array whose
#| two ends are held at 0. As written it is FTCS. Make it backward Euler.
def heat_step(u: Arr, r) {
    n = len(u)
    v = zeros(n)
    for i from 1 to n - 1 {
        v[i] = u[i] + r * (u[i - 1] - 2.0 * u[i] + u[i + 1])
    }
    return v
}

#| The mode $\sin(m\pi i/(N+1))$ on the grid, ends included.
def mode(m) {
    u = zeros(N_IN + 2)
    for i to N_IN + 2 {
        u[i] = sin(m * PI * i / (N_IN + 1.0))
    }
    return u
}

#| The checkerboard +1, -1, +1, ... on the interior, 0 at the ends.
def checkerboard() {
    u = zeros(N_IN + 2)
    s = 1.0
    for i from 1 to N_IN + 1 {
        u[i] = s
        s = -s
    }
    return u
}

#| u after NSTEP calls of heat_step with ratio r.
def run(u0: Arr, r) {
    u = u0 + 0.0
    for k to NSTEP {
        u = heat_step(u, r)
    }
    return u
}

#| The largest |u[i]|.
def max_abs(u: Arr) {
    m = 0.0
    for i to len(u) {
        m = max(m, abs(u[i]))
    }
    return m
}

show max_abs(run(checkerboard(), R_MAIN))
# FTCS multiplies the most jagged mode by about 1 - 4r per step; it is stable
# while that stays at or above -1: 1 - 4r = -1 at r = (1 - (-1)) / 4
ftcs_jagged_factor = 1.0 - 4.0 * R_MAIN
ftcs_r_limit = (1.0 - (-1.0)) / 4.0
four_r = 4.0 * R_MAIN
show ftcs_jagged_factor
show ftcs_r_limit
show four_r

test "the smoothest mode decays by 1 / (1 + 2.4 sin^2(pi/40))^16 = 0.790844 at r = 0.6" {
    g = 1.0 / (1.0 + 4.0 * R_MAIN * pow(sin(PI / 40.0), 2.0))
    u = run(mode(1.0), R_MAIN)
    assert_near(u[5] / mode(1.0)[5], pow(g, NSTEP), 1.0e-12)
    assert_near(pow(g, NSTEP), 0.7908438775198426, 1.0e-13)
}
test "the most jagged mode, m = 19, decays to 3.362e-9 of itself instead of growing" {
    g = 1.0 / (1.0 + 4.0 * R_MAIN * pow(sin(19.0 * PI / 40.0), 2.0))
    u = run(mode(19.0), R_MAIN)
    assert_near(u[5] / mode(19.0)[5], pow(g, NSTEP), 1.0e-15)
    assert_near(pow(g, NSTEP), 3.361996814033325e-09, 1.0e-20)
}
test "the checkerboard at r = 0.6 fades to numpy's 0.0068342 after 16 steps" {
    u = run(checkerboard(), R_MAIN)
    assert_near(max_abs(u), 0.006834176520732032, 1.0e-14)
    assert_near(u[10], 0.00318114592290907, 1.0e-14)
}
test "and at r = 5, ten times the FTCS limit, to numpy's 0.0012278" {
    u = run(checkerboard(), 5.0)
    assert_near(max_abs(u), 0.0012277856452884454, 1.0e-14)
}
test "the ends stay at 0" {
    u = run(checkerboard(), R_MAIN)
    assert_near(u[0], 0.0, 0.0)
    assert_near(u[N_IN + 1], 0.0, 0.0)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors