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

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

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

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$ .

Open in the workbench Download .vivOutputs below are from the site build’s run; Run all runs it again in your browser.

Look inside this program

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

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.

Your task

Write these 5 functions:

heat_step(u: Arr, r): 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.
mode(m): The mode \sin(m\pi i/(N+1)) on the grid, ends included.
checkerboard(): The checkerboard +1, -1, +1, ... on the interior, 0 at the ends.
run(u0: Arr, r): u after NSTEP calls of heat_step with ratio r.
max_abs(u: Arr): The largest |u[i]|.

Make these tests pass 4 of 5 fail

the smoothest mode decays by 1 / (1 + 2.4 sin^2(pi/40))^16 = 0.790844 at r = 0.6 assert_near: |0.788087 - 0.790844| = 0.00275737 > 1e-12 (fails)
the most jagged mode, m = 19, decays to 3.362e-9 of itself instead of growing assert_near: |183.794 - 3.362e-09| = 183.794 > 1e-15 (fails)
the checkerboard at r = 0.6 fades to numpy's 0.0068342 after 16 steps assert_near: |215.155 - 0.00683418| = 215.148 > 1e-14 (and 1 more) (fails)
and at r = 5, ten times the FTCS limit, to numpy's 0.0012278 assert_near: |2.88182e+20 - 0.00122779| = 2.88182e+20 > 1e-14 (fails)
the ends stay at 0 (passes)

Printed

215.155

ftcs_jagged_factor: -1.4

ftcs_r_limit: 0.5

four_r: 2.4

Your task

Solutions 0