Task 4 of 5

Solve, Don’t Step

The whole problem is that the explicit step evaluates the Laplacian on the field it is leaving. Evaluate it on the field it is arriving at instead — that is backward Euler — and the step size limit vanishes entirely. Unconditionally stable, at any dt, forever.

The catch is visible the moment you write it down. The unknown is on both sides:

u' = u + α·∇²u'
⟺  (1 + 4α)·u' − α·(neighbours of u') = u

That is not a formula you evaluate, it is a linear system — one equation per cell, 1,024 of them on this grid, all coupled. Solve it the way GPUs like: rearrange each equation for its own cell, then iterate.

u'[c] = ( u[c] + α·(neighbours of u') )
        / (1 + 4α)

Every cell reads only the previous iterate, so all 1,024 can be computed at once — that is a Jacobi sweep, and the Iterative Solvers module takes it much further (red-black ordering, residuals, why it beats Gauss–Seidel on a GPU). Here 25 sweeps is plenty, because the 1 in 1 + 4α makes this system diagonally dominant and easy. One thing must not slip: u on the right is the old time level and never changes during the solve. Only the guess moves.

one word changes — 'new' — and the step becomes a system of equations

Goal: write the Jacobi sweep, then iterate it 25 times to take a single implicit step at dt = 2 — four times the explicit limit.

Requirements

The sweep returns (uOld[y][x] + α · (four neighbours of guess)) / (1 + 4α)
The centre term comes from uOld; only the four neighbours come from guess
Iterate SWEEPS times, passing the same seed as uOld every time

Hint 1 — where the division comes from

Collect the unknown cell on the left of u' = u + α·(l + r + up + dn − 4u'): the −4α·u' moves over as +4α·u', giving (1 + 4α)·u' = u + α·(l + r + up + dn). Divide.

Hint 2 — the sweep body

const neighbours = guess[y][xl] + guess[y][xr]
  + guess[yd][x] + guess[yu][x];
return (uOld[y][x] + this.constants.alpha * neighbours)
  / (1 + 4 * this.constants.alpha);

Hint 3 — why the starter’s loop is wrong

sweep(guess, guess) replaces the right-hand side with the current iterate every sweep, which throws away the one piece of information that makes this a time step. It still converges — to ∇²u = 0, a flat field. The fix is one word: guess = await sweep(seed, guess);

Same idea elsewhere

"The implicit step is a linear solve" is the fork in the road for every production simulator: implicit thermal and structural codes hand (I − αL) to a Krylov solver with a preconditioner, and GPU fluid solvers run exactly this Jacobi (or a multigrid V-cycle) for the pressure projection every frame. cuSPARSE, rocSPARSE and every WebGPU fluid demo you have seen are all standing on this one rearrangement.

All tasks in The Heat Equation & Stability

This page is an interactive exercise — the editor, the GPU runner and your saved progress need JavaScript. The text above is the full brief.