The walk(M, N, L) function
M walkers, N steps of length L, angle uniform on [0, 2π). Not a single Python loop: numpy arrays do all the work.
Solution
walk.pyimport numpy as np
def walk(M, N, L=1.0, rng=None):
rng = np.random.default_rng(rng)
theta = rng.uniform(0, 2*np.pi, size=(M, N)) # all angles
x = np.zeros((M, N + 1))
y = np.zeros((M, N + 1))
x[:, 1:] = np.cumsum(L*np.cos(theta), axis=1) # Σ steps
y[:, 1:] = np.cumsum(L*np.sin(theta), axis=1)
return x, y # (M, N+1)- One matrix of angles,
theta, of shape(M, N): a row is a walker, a column is a step number.default_rngproduces all M·N numbers in one call. - A step is a vector:
(L cos θ, L sin θ), computed element-wise for the whole matrix at once. - Position is the running sum of steps:
np.cumsum(..., axis=1)adds along each row. That is the "loop over steps", only it runs in C instead of Python. - Column 0 is the origin. That is why the arrays are one column wider:
(M, N+1).
What happens inside
theta.shape = (5, 10): each arrow is one random angle
Five walkers, N = 1000
Same rules, same starting point, and completely different fates.
Equal scales on both axes. ★ is the start (0, 0), circles are the end points. The dashed circle is the "typical" radius √N·L ≈ 31.6 L.
Each trajectory is one realisation of a random process: the sum of 1000 independent random vectors. The rule is the same, but every walker gets its own sequence of angles.
1. No direction is preferred. ⟨r⟩ = 0: on average a walker goes nowhere, so each one drifts off in its own random direction.
2. The distance does not self-average. The typical scale is √N·L ≈ 32 L, but for large N the distribution of r² is exponential:
The spread equals the mean itself: 100 %! About one walker in five ends up closer than 0.5√N·L, and one in ten ends up farther than 1.5√N·L. A single trajectory tells you nothing about the average.
3. The shape is a fractal. The path is self-similar (fractal dimension 2): dense "tangles" where the walker lingers alternate with long random excursions. In the plane the walk is recurrent: it crosses itself again and again, and every walker draws its own pattern.
Check: 1000 walkers, r²/(NL²) at the end
⟨r²⟩ as a function of N
M = 1000 walkers, N = 1000 steps. We average r² over the ensemble at every step and compare it with the theory, NL².
The cloud of all M walkers at step 1000. Hover over the chart on the left and the cloud jumps to that step. The solid circle is √⟨r²⟩, the dashed one is √N·L.
The position after N steps is the sum of the step vectors. Square it:
The first sum is exactly NL². In the second, every term \(L^2\cos(\theta_i-\theta_j)\) averages to zero: the steps are independent and each direction is uniform. Therefore
This holds for every N, not just asymptotically. It follows that √⟨r²⟩ ∝ √N: to get twice as far you need four times as many steps. That is diffusion.
The slope of the least-squares line comes out at about 1·L² (about 1 for L = 1). In my Python run (walk.py, seed = 1) it is 0.91 for aN+b and 0.95 for the fit through the origin.
Why not exactly 1? For a single walker σ(r²) = ⟨r²⟩, so the relative error of a mean over M walkers is ≈ 1/√M ≈ 3 %. On top of that, the points on the curve are strongly correlated (they are the same walkers at different steps), so the whole curve wanders together. A deviation of 5–8 % at M = 1000 is normal.
Move the M slider: at M = 10 000 the curve almost sits on the theory line, at M = 10 it jumps around. The shaded band shows ±1σ statistical error of the mean.
Code for part 3
walk.pyx, y = walk(1000, 1000, L=1.0, rng=1)
r2 = np.mean(x**2 + y**2, axis=0) # over walkers → (N+1,)
n = np.arange(1001)
slope, intercept = np.polyfit(n, r2, 1) # least squares: a*N + b
slope0 = np.sum(n*r2) / np.sum(n*n) # fit through the origin
plt.plot(n, r2); plt.plot(n, n*L**2, '--')The same data on log–log axes
A logarithm turns a power law into a straight line, and the exponent into its slope.
If \(\langle r^2\rangle = C\,N^{\alpha}\), then
A power law becomes a straight line, and its slope is the exponent α itself. The fit gives α ≈ 1 (0.97 in the Python run), i.e. ⟨r²⟩ ∝ N¹: normal diffusion.
Why this is a good test:
- The slope does not depend on the prefactor. Change L in part 3: the line moves up or down by log L², but the slope stays at 1. We test the law, not the constant.
- The eye judges a straight line reliably, but not curvature on a linear plot: N and N1.15 look almost the same on linear axes (see below).
- All scales are visible at once: every decade (1–10, 10–100, 100–1000) gets the same space, so small N is not squashed into a corner.
- Different regimes have different slopes: ballistic (α = 2), diffusive (α = 1) and subdiffusive (α < 1) motion are told apart at a glance.
Experiment: N versus N1.15
Both curves meet at N = 1000. On linear axes the difference is easy to mistake for noise; on log–log axes they are clearly two straight lines with different slopes.
Changing the rules
A line, a grid, 3D and a "lazy" walker. Does ⟨r²⟩ = NL² still hold?
| Rule | Step | ⟨|l|²⟩ | theory ⟨r²⟩/N | simulation ⟨r²⟩/N | α (log–log) |
|---|
The derivation in part 3 used only two facts: the steps are independent and zero on average. Then the cross terms vanish for any geometry and any dimension:
- (a) Line ±L: l² = L² always → NL². The cleanest case.
- (b) Grid N/S/E/W: each step is orthogonal or opposite to any other, ⟨l_i·l_j⟩ = 0, |l|² = L² → NL². The anisotropy of the lattice does not affect ⟨r²⟩.
- (c) 3D: |l| = L, ⟨l⟩ = 0 → NL² again. The dimension changes the diffusion coefficient \(D = L^2/(2d\,\tau)\), but ⟨r²⟩ = 2dDt = NL² stays the same.
With probability ½ the step is zero; otherwise it is an ordinary step of length L. So
The NL² law breaks, but only in the coefficient. Growth is still linear (α = 1; on log–log it is a line with the same slope, shifted down by log 2). The slope of the linear fit is ≈ 0.5.
Physically: out of N "ticks" the walker actually makes only N/2 steps on average, so its clock runs half as fast and the diffusion coefficient is halved. The cloud spreads just like the ordinary one, but takes twice as long. In general, with a probability of moving p: ⟨r²⟩ = pNL².
Why "two uniform angles" is wrong
✗ θ ∈ U[0, π], φ ∈ U[0, 2π): points bunch up at the poles.
✓ (g₁, g₂, g₃)/|g|, gᵢ ~ N(0, 1): uniform over the sphere.
The area element of a sphere is \(dA=\sin\theta\,d\theta\,d\varphi\). If θ is uniform, every latitude band gets the same number of points, but the bands near the poles have tiny area, so the density goes as 1/sin θ. The histograms under the spheres show z = cos θ: for the correct method it is flat (Archimedes' hat-box theorem); for the wrong one it has humps at ±1.
Why the Gaussian vector works: the density \(\propto e^{-(g_1^2+g_2^2+g_3^2)/2}\) depends only on the length |g|. It is spherically symmetric, so the direction is uniform. An alternative: \(\cos\theta\sim U[-1,1]\).
A curious detail: even with the wrong angles ⟨l⟩ = 0 and |l| = L, so ⟨r²⟩ = NL² still holds! What breaks is isotropy: ⟨z²⟩ = NL²/2 instead of NL²/3, and the cloud is stretched along the z axis (last row of the table).
def walk_3d(M, N, L=1.0, rng=None):
rng = np.random.default_rng(rng)
v = rng.standard_normal((M, N, 3)) # isotropic vector
v /= np.linalg.norm(v, axis=2, keepdims=True) # → onto the sphere
steps = L * v
return np.cumsum(steps, axis=1) # (+ start at 0)
def walk_lazy(M, N, L=1.0, rng=None):
rng = np.random.default_rng(rng)
theta = rng.uniform(0, 2*np.pi, size=(M, N))
move = rng.random((M, N)) < 0.5 # moves?
dx, dy = L*np.cos(theta)*move, L*np.sin(theta)*move
return np.cumsum(dx, axis=1), np.cumsum(dy, axis=1)Summary
What to take away.
The typical distance grows as the square root of the number of steps. That is why diffusion is slow over large distances.
⟨r²⟩ = N⟨|l|²⟩ for any independent zero-mean steps: in 1D, 2D, 3D and on a lattice.
The slope on log–log axes is the exponent. It depends neither on L nor on how lazy the walker is.
The spread of r² for a single walker equals the mean. The law is a property of the ensemble, not of an individual trajectory.