# Solve, Don’t Step

*Task 4 of 5 · [The Heat Equation & Stability](https://gpu.rocks/learn/heat-and-stability-514063bb.md) · GPU.js Learn*

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:

```js
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.

```js
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.

## Figures

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

## Goal

**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

```js
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.

## Starter code

```js
// Backward Euler: the new field appears on BOTH sides of the equation.
// Solve it with Jacobi sweeps — every cell reads the previous iterate.
const gpu = new GPU({ mode });

const D = 8;
const dx = 4;
const dt = 2;                       // 4× the explicit limit of 0.5
const ALPHA = D * dt / (dx * dx);   // = 1
const SWEEPS = 25;

const sweep = gpu.createKernel(function (uOld, guess) {
  const x = this.thread.x;
  const y = this.thread.y;
  let xl = x - 1; if (xl < 0) xl = this.constants.size - 1;
  let xr = x + 1; if (xr > this.constants.size - 1) xr = 0;
  let yd = y - 1; if (yd < 0) yd = this.constants.size - 1;
  let yu = y + 1; if (yu > this.constants.size - 1) yu = 0;
  // TODO: one Jacobi sweep —
  //   (uOld[y][x] + alpha * (the four neighbours of GUESS)) / (1 + 4 * alpha)
  return guess[y][x];
}, {
  output: [32, 32],
  constants: { size: 32, alpha: ALPHA },
});

function hottest(u) {
  let m = 0;
  for (let y = 0; y < u.length; y++) {
    for (let x = 0; x < u[y].length; x++) {
      const a = Math.abs(u[y][x]);
      if (!(a <= m)) m = a;
    }
  }
  return m;
}

function total(u) {
  let s = 0;
  for (let y = 0; y < u.length; y++) for (let x = 0; x < u[y].length; x++) s += u[y][x];
  return s;
}

let guess = seed;
for (let k = 0; k < SWEEPS; k++) {
  // TODO: the right-hand side is the OLD field and never changes during a
  // solve. This passes the current iterate instead, which throws the time
  // step away and converges to a flat field.
  guess = await sweep(guess, guess);
}

console.log('hottest after one implicit step:', hottest(guess));
console.log('total heat (unchanged at 36):', total(guess));
```

---

Interactive version: https://gpu.rocks/learn/heat-and-stability-514063bb/4

[Previous task](https://gpu.rocks/learn/heat-and-stability-514063bb/3.md) · [Next task](https://gpu.rocks/learn/heat-and-stability-514063bb/5.md)