How heat forgets: solving diffusion in a single call

The problem

Drop ink in still water and it spreads. Warm one corner of a metal sheet and the heat creeps outward. Blur a photo and the fine detail goes first while the broad shapes hang on.

All three are the same process, and it is one of the most useful things to be able to predict — it is how you work out where a pollutant will be tomorrow, how a chip sheds heat, how an image filter behaves. The usual way to find out is to simulate every moment in between: chop time into a thousand slices and step through them one at a time. Each step is slightly wrong, and each wrong answer becomes the input to the next. Want a different moment? Run the whole thing again.

This walkthrough skips the stepping. You ask for a time and get the answer in two requests — and then, for fun, you ask for a time in the past, and watch it fail in exactly the way the physics says it should.

You need a running HeatherDB and numpy. No GPU, no solver library.

Why this works at all

Any pattern can be built out of ripples: one big slow wave, plus a finer one on top, plus a finer one still, all the way down to the tiniest wiggle. Add enough ripples together and you can make any shape you like. Splitting a pattern into those ripples is what a Fourier transform does, and numpy does it in one line.

Here is the useful fact about spreading. Each ripple fades on its own. The tiny wiggles vanish quickly, the big slow waves hang around, and no ripple interferes with any other on the way down. That is why a blurred photo keeps its broad shapes and loses the fine ones — the fine ones simply faded first.

So predicting the future of a spreading pattern is not hard once you look at it this way. You do not simulate anything. You take the list of ripples, multiply each one by how much it has faded, and add them back up.

Two things make this cheap on a database. First, multiplying two ripple-lists together is exactly what the engine's bind operation does — bind blends two patterns, and blending in the everyday world is multiplication in ripple-world. Second, "how much has each ripple faded after t seconds" is just "how much it fades in one second, raised to the power t" — and the engine has an operation that raises a pattern to a power. So the fading is a pattern, and the time is its exponent.

That is the whole trick. One pattern that means one second of spreading, one exponent that means how long, and one blend to apply it.

Step 1 — build the clock

N, ALPHA = 64, 0.005
x  = np.arange(N) * 2 * np.pi / N
u0 = np.sin(x) + 0.5 * np.sin(3 * x)        # two modes; the finer one dies first
w  = np.fft.rfftfreq(N, d=1.0 / N)

spec = np.exp(-ALPHA * w**2); spec[0] = 1.0; spec[-1] = 1.0
clock = np.fft.irfft(spec, n=N)

One detail will bite you otherwise. Two entries in that ripple-list — the flattest one and the sharpest one — have to be set to 1 by hand, which is what those last two lines do. Raising a negative number to a fractional power (the square root of −1, say) does not give you a normal number back, and the engine would have to hand you something imaginary. Pinning those two entries keeps every answer real.

Step 2 — evolve to any time

u = bind(u0, power(clock, t))     # one pow, one bind, any t

Against the closed-form solution:

  t      max|err| vs closed form
  0.5      6.661e-16
  4.0      9.992e-16
 20.0      3.220e-15

Read the column, not the rows. The error does not grow with t. There is no step for error to accumulate in — t = 20 costs exactly what t = 0.5 costs, because both are one exponent. A marching solver cannot give you this; its error is a function of how many steps it took to arrive.

Step 3 — run it backwards

Negate the exponent. There is no separate inverse routine.

inverse = power(clock, -2.0)
recovered = bind(forward, inverse)

  reversing from t = 2
  clean          2.130e-12
  1e-06 noise    5.275e-03
  1e-03 noise    7.095e+00
  t = 8, 1e-6    2.255e+10

On clean data you get the original field back to 2e-12. Add a part per million of noise and accuracy drops four orders of magnitude. Reverse from t = 8 instead of t = 2 with that same tiny noise and the answer has no information in it at all.

This is not the method failing — it is the honest answer. Spreading destroys the fine detail. Going backwards means amplifying whatever is left of it, and if what is left is mostly measurement noise, you amplify the noise. Un-mixing milk from coffee has the same problem: nothing is stopping you in principle, but the information you would need is gone.

The practical rule: how far back you can go depends on how clean your data is. Cleaner measurements buy you more history, and there is a point past which no amount of cleverness recovers anything.

Worth stating plainly, because it is the practical takeaway: a deblurring routine that hands you a clean picture here has chosen a prior on your behalf and not mentioned it. This one shows you the blow-up instead.

The whole thing

import json, math, urllib.request, base64
import numpy as np

BASE = "http://localhost:6380"
AUTH = base64.b64encode(b"admin:quickstart-password").decode()

def vec(op, body):
    """Call one of the stateless /vec routes. No collection, no writes —
    this is arithmetic the engine performs on vectors you hand it."""
    req = urllib.request.Request(
        f"{BASE}/vec/{op}", method="POST", data=json.dumps(body).encode(),
        headers={"Content-Type": "application/json", "Authorization": "Basic " + AUTH})
    return np.array(json.load(urllib.request.urlopen(req))["result"])

# bind blends two patterns together. Left unnormalised, blending is exactly
# ripple-by-ripple multiplication, which is what makes it a physics operator.
def bind(a, b): return vec("bind", {"a": list(map(float, a)), "b": list(map(float, b)),
                                    "normalize": False})
def power(a, t): return vec("pow", {"a": list(map(float, a)), "t": float(t)})

N, ALPHA = 64, 0.005
x  = np.arange(N) * 2 * np.pi / N
u0 = np.sin(x) + 0.5 * np.sin(3 * x)        # two modes; the finer one dies first
w  = np.fft.rfftfreq(N, d=1.0 / N)

# The heat kernel for one unit of time. DC and Nyquist are held at 1 so that
# fractional powers of this vector stay real — the constraint /vec/pow documents.
spec = np.exp(-ALPHA * w**2); spec[0] = 1.0; spec[-1] = 1.0
clock = np.fft.irfft(spec, n=N)

def closed_form(t):
    s = np.exp(-ALPHA * w**2 * t); s[0] = 1.0; s[-1] = 1.0
    return np.fft.irfft(np.fft.rfft(u0) * s, n=N)

print("  t      max|err| vs closed form")
for t in (0.5, 4.0, 20.0):
    u = bind(u0, power(clock, t))            # any timestamp: one pow, one bind
    print(f"{t:5.1f}      {np.abs(u - closed_form(t)).max():.3e}")

print("\n  reversing from t = 2")
forward = bind(u0, power(clock, 2.0))
inverse = power(clock, -2.0)                 # negative exponent runs it backwards
rng = np.random.default_rng(0)
print(f"  clean          {np.abs(bind(forward, inverse) - u0).max():.3e}")
for eps in (1e-6, 1e-3):
    noisy = forward + eps * rng.standard_normal(N)
    print(f"  {eps:.0e} noise    {np.abs(bind(noisy, inverse) - u0).max():.3e}")

far = bind(u0, power(clock, 8.0)) + 1e-6 * rng.standard_normal(N)
print(f"  t = 8, 1e-6    {np.abs(bind(far, power(clock, -8.0)) - u0).max():.3e}")

What to change for your own field

An image instead of a line. Run the same clock across the rows, then across the columns. Twice the calls, nothing new to learn.

A different process. Spreading is one pattern in that slot; others go in the same slot. Drift — everything sliding steadily in one direction, like pollution carried by wind — is a different pattern. So is a camera blur. Build the pattern once, raise it to a power, blend it.

Edges wrap around. The blend treats your data as a loop, so heat leaving the right-hand edge reappears on the left. If your problem has real walls — a metal plate that ends, a room with sides — you either pad the data out with empty space so nothing reaches the edge, or you use a different method. Check this before anything else; it is the thing most likely to rule your problem out.

The same two calls, other processes

Only the pattern in the slot changes.

A wave that steepens into a shock. Traffic bunching up, or water piling into a breaker, is Burgers' equation — and it feeds back on itself, which the trick does not allow. There is a change of variables (Cole–Hopf) that converts it into plain spreading, does the one blend, and converts back. Against a 15,000-step conventional simulation, the error is 1.67e-10 early and 2.26e-12 later, and the wave still steepens exactly as it should — the steepening is real, it has just been moved out of the way of the clock.

Option prices. Black–Scholes is spreading plus drift plus a discount, and each of those three is recognisable in the pattern: how fast prices smear out is volatility, the sliding is the drift, and the shrinking is the discount rate. Price is one blend of the payoff with the clock; the option's delta is one more blend. Across four maturities and five spot prices, the worst price error is 1.61e-4 on options worth $0.06 to $32, and the worst delta error 7.15e-6.

A swinging pendulum with friction. The degenerate case: one ripple instead of a whole field. Stepping it 25 times with the clock holds 8.29e-16 of error where the textbook step-by-step method reaches 59.2. Collapsing those 25 steps into a single exponent costs 6.31e-16 — the same answer, one call.

Where the method stops

Three conditions, and all three are real:

The process must not feed back on itself. Spreading does not care how much has already spread. A wave that steepens because it is already steep does — that is why Burgers needs the change of variables above, which converts the feedback away before the clock touches it.

It must behave the same everywhere. If one half of your plate conducts heat faster than the other, the ripples stop fading independently and the trick breaks.

The edges must wrap, or not matter. As above.

Miss any one and you want an ordinary solver. Meet all three — and a surprising amount of everyday physics does — and the answer is two calls away, at whatever moment you ask for.

More depth: Solve a PDE by binding and Physics in the frequency domain.