Diffusion Simulation

One process, two descriptions

The domain [−10, 10] is cut into n_x cells of width Δx. All of the concentration, and all of the walkers, start in the middle cell. At every time step Δt:

• Each walker hops one cell left with probability r, one cell right with probability r, and stays put otherwise, where r = DΔt/Δx².
• The concentration is updated by the explicit finite-difference (FTCS) scheme below.

The fraction of walkers expected in cell i after one step is r times the fraction in each neighbor plus (1 − 2r) times its own, which is exactly the FTCS update. So the finite-difference solution is the expected histogram of the walkers, step by step, and the only difference between the blue bars and the orange line is sampling noise. The lab reports that difference as an L1 distance, next to the value expected from binomial noise alone; quadruple the number of walkers and both halve.

The walls follow the same rule on both sides: a walker that tries to hop through a wall stays where it is, and the finite-difference scheme copies the edge cell as its missing neighbor. Both make the walls reflecting, so no mass is lost; the lab reports the total drift, which is rounding error. Press Play to watch the peak flatten, and set T = 10 with D = 2 to see the profile reach the walls and level off.

The diffusion equation and its exact solution

If particles move down concentration gradients with flux J = −D ∂u/∂x (Fick's law) and are neither created nor destroyed, the concentration obeys

the same equation that describes heat conduction. For a unit mass released at x₀ on an infinite line, the solution is a Gaussian whose variance grows linearly in time:

Reflecting walls at ±L impose ∂u/∂x = 0 there. The method of images handles them exactly: mirror the source in each wall, mirror the mirrors, and add all the Gaussians. The dashed curve is that sum, so it stays exact even after the profile reaches the walls. The variance then stops growing like 2Dt and approaches L²/3 ≈ 33.3, the variance of a uniform distribution on [−10, 10].

The link to single particles is Einstein's: a particle doing a random walk or Brownian motion has a position distribution that obeys this equation, and D sets its mean squared displacement.

Finite differences and the stability limit

Replacing the time derivative with a forward difference and the space derivative with a centered difference gives the FTCS scheme:

Its error is first order in Δt and second order in Δx. The scheme is stable only for r ≤ ½. The walker picture shows why: the probability of staying put, 1 − 2r, becomes negative beyond that, and in the PDE the same negative weight lets short-wavelength wiggles grow at every step until they swamp the solution. The lab therefore picks the number of steps so that r ≤ 0.4 and reports the value it used.

The price is that Δt must shrink with Δx²: doubling the resolution means twice the cells and four times the steps, eight times the work. Implicit schemes such as Crank–Nicolson remove the stability limit at the cost of solving a linear system per step. Two exact properties survive the discretization: total mass is conserved, and away from the walls the discrete variance grows by exactly 2rΔx² = 2DΔt per step, so the finite-difference variance equals 2Dt until the profile feels the walls. The remaining L1 difference from the exact solution comes from starting with the mass spread over one cell instead of a point and from the truncation error. It is largest for small T, when the Gaussian is only a few cells wide.

Example result

With D = 1, T = 2, 101 cells and 4,000 walkers (seed 42), the scheme takes 128 steps with r = 0.398. The finite-difference variance is 4.000 against 2DT = 4.000, and the walkers’ variance is 3.996. The finite-difference solution differs from the exact one by 0.0016 in L1, while the walker histogram differs from the finite-difference solution by 0.091, in line with the 0.089 expected from sampling noise.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

// FTCS (forward time, centered space) step of ∂u/∂t = D ∂²u/∂x²:
//   u_i ← u_i + r (u_{i−1} − 2u_i + u_{i+1}),  r = DΔt/Δx².
// Reflecting (no-flux) walls: a missing neighbour is replaced by the cell itself.
export function stepPDE(u, r) {
  const n = u.length, out = new Float64Array(n);
  for (let i = 0; i < n; i++) {
    const left = i > 0 ? u[i - 1] : u[i];
    const right = i < n - 1 ? u[i + 1] : u[i];
    out[i] = u[i] + r * (left - 2 * u[i] + right);
  }
  return out;
}

// The matching random walk on the same cells: hop left with probability r,
// right with probability r, stay otherwise; a hop into a wall is cancelled.
// The expected histogram of these walkers obeys exactly the FTCS update above.
export function stepWalkers(pos, r, n, rng) {
  for (let k = 0; k < pos.length; k++) {
    const u = rng.uniform();
    if (u < r) { if (pos[k] > 0) pos[k]--; } else if (u < 2 * r) { if (pos[k] < n - 1) pos[k]++; }
  }
}

def step_pde(u, r):
    """FTCS step u_i += r*(u_{i-1} - 2u_i + u_{i+1}); a missing neighbour at a wall is the cell itself."""
    n = len(u)
    out = [0.0] * n
    for i in range(n):
        left = u[i - 1] if i > 0 else u[i]
        right = u[i + 1] if i < n - 1 else u[i]
        out[i] = u[i] + r * (left - 2 * u[i] + right)
    return out


def step_walkers(pos, r, n, rng):
    """Hop left w.p. r, right w.p. r, else stay; hops into a wall are cancelled."""
    for k in range(len(pos)):
        u = rng.uniform()
        if u < r:
            if pos[k] > 0:
                pos[k] -= 1
        elif u < 2 * r:
            if pos[k] < n - 1:
                pos[k] += 1

Run the full model from a terminal:

node simulation.js --diffusion 1 --time 2 --seed 42
python simulation.py --diffusion 1 --time 2 --seed 42

Download the model

Assumptions & limitations

• One dimension, constant D. No drift (advection), no reactions, a uniform medium.
• Reflecting walls at ±10. Both the walkers and the finite-difference scheme conserve mass; the exact comparison uses the method of images for the same walls.
• Explicit time stepping. FTCS with r ≤ 0.4; the time step is chosen automatically from D, T and Δx.
• Lattice walkers. Walkers hop between cell centers, so their histogram is exact in expectation for the discrete scheme, not for the continuous equation. Both converge to the continuum as Δx → 0.
• Independent walkers. No crowding or interaction; concentration is the density of non-interacting particles.

Related simulations

• OS-012Random walkCoin-flip steps on a line or grid: distance grows like √n, and in 1-D and 2-D every walker comes back.
• OS-013Brownian motionParticles jiggle with Gaussian steps; their mean squared displacement grows as 4Dt, as Einstein predicted.
• OS-010Galton boardBalls bouncing left or right through rows of pegs pile up into a binomial bell curve.
• OS-011Central limit theoremAverage n draws from almost any distribution and the means turn normal — at a speed you can measure.

Background: Euler vs Runge–Kutta: Time Stepping in Simulations · How Accurate Is a Monte Carlo Estimate?

References

1. Crank, J. (1975). The Mathematics of Diffusion, 2nd ed. Oxford University Press.
2. Einstein, A. (1905). Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Annalen der Physik, 322(8), 549–560. doi.org/10.1002/andp.19053220806
3. Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing, 3rd ed. Cambridge University Press.