# Stable Is Not Accurate

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

Time to put the two schemes on the same clock. Both runs below finish at
`T = 16`: the explicit one in small steps of `0.2`, the implicit
one in steps of `2` — ten times larger, and four times past the explicit
limit. A third run takes the explicit scheme at the implicit step size, for the pleasure
of watching it fail in eight steps.

The two survivors will not agree. Backward Euler is *first-order* accurate,
just like forward Euler, so a ten-times-larger step carries a ten-times-larger error —
it damps sharp features harder than the real equation does. Expect the two fields to
differ by a few percent of the peak. That is the honest trade, and it is worth stating
plainly: **unconditional stability is not accuracy**. What implicit
stepping buys you is the right to choose `dt` for the accuracy you need,
rather than having it dictated by the smallest cell in your mesh.

Nor is it free. Each implicit step here costs 25 sweeps, so the coarse run makes
`8 × 25 = 200` kernel launches against the fine run's 80 — implicit
*loses* the launch count at this size. It wins when the explicit limit gets
brutal: refine `dx` by 10× and the explicit run needs 100× the steps, while
the implicit one needs the same eight and a slightly harder solve.

## Goal

**Goal:** finish the Jacobi sweep with `α` arriving as an
argument, work out how many steps of each size reach `T`, and run all three.

## Requirements

- The sweep is the one from the last task, but `alpha` is a kernel *argument*, not a constant
- `smallSteps` and `bigSteps` are the counts that reach `T = 16` at `dt = 0.2` and `dt = 2`
- Run all three: explicit at the small step, implicit at the big one, explicit at the big one
- Log the hottest cell of each, and the largest gap between the two survivors

## Hint 1 — alpha as an argument

Only the spelling changes: a plain `alpha` where the constant used to
be, and the caller passes it. Everything else is last task's body:

```js
return (uOld[y][x] + alpha * neighbours)
  / (1 + 4 * alpha);
```

## Hint 2 — how many steps?

Steps × step size = elapsed time, so it is `T / dt`:
`16 / 0.2 = 80` and `16 / 2 = 8`.

## Same idea elsewhere

Choosing a scheme by what limits it — accuracy or stability — is the daily work
of numerical simulation everywhere: stiff chemistry and implicit thermal solvers pay for
a linear solve per step because the explicit alternative would need millions of them,
while explicit codes dominate wave propagation and particle work, where the stability
step is close to the accuracy step anyway. On a GPU the arithmetic is nearly free, so
the calculus is really about launches and memory traffic per unit of simulated time.

## Starter code

```js
// Same physics, same finish time, two step sizes.
const gpu = new GPU({ mode });

const D = 8;
const dx = 4;
const T = 16;           // finish time
const DT_SMALL = 0.2;   // 0.4× the explicit limit (0.5)
const DT_BIG = 2;       // 4×   the explicit limit
const SWEEPS = 25;

const explicitStep = gpu.createKernel(function (u, alpha) {
  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;
  const c = u[y][x];
  return c + alpha * (u[y][xl] + u[y][xr] + u[yd][x] + u[yu][x] - 4 * c);
}, { output: [32, 32], constants: { size: 32 } });

const sweep = gpu.createKernel(function (uOld, guess, alpha) {
  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: last task's sweep, with alpha coming in as an argument
  return guess[y][x];
}, { output: [32, 32], constants: { size: 32 } });

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 gap(a, b) {
  let m = 0;
  for (let y = 0; y < a.length; y++) {
    for (let x = 0; x < a[y].length; x++) {
      const d = Math.abs(a[y][x] - b[y][x]);
      if (!(d <= m)) m = d;
    }
  }
  return m;
}

async function runExplicit(dt, steps) {
  const alpha = D * dt / (dx * dx);
  let u = seed;
  for (let i = 0; i < steps; i++) u = await explicitStep(u, alpha);
  return u;
}

async function runImplicit(dt, steps) {
  const alpha = D * dt / (dx * dx);
  let u = seed;
  for (let i = 0; i < steps; i++) {
    let guess = u;
    for (let k = 0; k < SWEEPS; k++) guess = await sweep(u, guess, alpha);
    u = guess;
  }
  return u;
}

// TODO: how many steps of each size land exactly on time T?
const smallSteps = 0;
const bigSteps = 0;

const fine = await runExplicit(DT_SMALL, smallSteps);
const coarse = await runImplicit(DT_BIG, bigSteps);
const doomed = await runExplicit(DT_BIG, bigSteps);

console.log('explicit, dt =', DT_SMALL, 'x', smallSteps, 'steps → hottest', hottest(fine));
console.log('implicit, dt =', DT_BIG, 'x', bigSteps, 'steps → hottest', hottest(coarse));
console.log('explicit at dt =', DT_BIG, '→ hottest', hottest(doomed));
console.log('largest gap between the two survivors:', gap(fine, coarse));
```

---

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

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