Boids flocking simulation

Three local rules, one flock

Craig Reynolds introduced boids in 1987 to animate flocks for computer graphics without scripting each bird. His insight was that a flock needs no leader and no global plan: each individual steering by what its nearby flockmates do is enough. The lab starts every boid at a random position with a random heading and lets the rules run.

• Order appears. Polarization, the length of the average heading vector, starts near zero because random headings cancel out, and climbs toward 1 as neighborhoods align and merge.
• Groups merge. Small aligned groups form first; when two meet, alignment and cohesion fuse them. The green curve tracks the share of boids in the largest group. Because groups are chains of boids within r of each other, even the random start contains a few large, loose groups, and the curve jumps as groups touch and separate before one flock holds together.
• Spacing is regulated. Separation keeps a typical distance between neighbors; cohesion keeps groups from dispersing.

Each rule can be switched off by setting its weight to zero. Without alignment the boids clump but mill about with disordered headings. Without separation they collapse onto each other, and the nearest-neighbor distance drops almost to zero. Without cohesion they still align but spread out more loosely. A small perception radius leaves many separate flocks.

The steering rules

Boid i has position x_i and velocity v_i. Its neighbors N_i are the boids within the perception radius r, measured across the wrapped edges. With d_ij = x_j − x_i, each step every boid updates its velocity by

where v̄_i is the mean velocity and d̄_i the mean offset of the neighbors (the vector to their center). The w's are the weights in the lab; the fixed gains k_s = 0.05, k_a = 0.05 and k_c = 0.005 put the three terms on comparable scales. Separation is weighted by 1/distance, so it dominates at close range.

1. For every pair within r, accumulate both boids' neighbor sums (each pair is visited once).
2. For every boid, add Δv, then clamp the speed to [v_max/2, v_max] by rescaling.
3. Move every boid by its new velocity and wrap positions into the 100 × 100 world.

All boids update synchronously from the same snapshot. Reynolds' boids also had a limited field of view and steered around obstacles; this version keeps only the three flocking rules.

Measuring a flock

• Polarization = |(1/n) Σ v_i/|v_i||, between 0 and 1. It is the standard order parameter in physics models of collective motion. For independent random headings it is not zero but about √(π/4n), which is why the starting value is small but positive.
• Flocks are connected groups: two boids are linked if they are within r of each other, and links chain. A flock can be a long, thin string of boids.
• Mean nearest-neighbor distance shows the spacing that separation and cohesion negotiate.

Flocking dynamics are sensitive to tiny differences: a last-bit rounding change grows until the flock takes a different path. So that a seed gives the same flight in the browser and in Python, the model uses only addition, subtraction, multiplication, division and square roots — operations that IEEE 754 arithmetic rounds identically in both languages — and no trigonometry. The neighbor search compares every pair, so the cost grows with n²; that is why the number of boids is capped at 300.

Example result

With 120 boids, weights 1.5 / 1 / 1, radius 10 and seed 42, polarization rises from 0.058 at the random start to 0.999 after 600 steps. The boids end in 1 flock of 120, with a mean nearest-neighbor distance of 3.09 units.

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.

// One synchronous step. Every boid looks at the others within radius r:
//   separation: push away from boids closer than r/2, weighted by 1/distance
//   alignment:  steer toward the neighbors' mean velocity
//   cohesion:   steer toward the neighbors' mean position
// Then clamp the speed to [max/2, max] and move; the world wraps around.
export function step(b, p) {
  const { n, x, y, vx, vy } = b;
  const r2 = p.radius * p.radius, s2 = r2 / 4;
  const cnt = new Float64Array(n), ax = new Float64Array(n), ay = new Float64Array(n);  // mean velocity
  const cx = new Float64Array(n), cy = new Float64Array(n);                             // mean offset
  const sx = new Float64Array(n), sy = new Float64Array(n);                             // separation
  for (let i = 0; i < n; i++) {
    for (let j = i + 1; j < n; j++) {
      const dx = wrapDelta(x[j] - x[i]), dy = wrapDelta(y[j] - y[i]);
      const d2 = dx * dx + dy * dy;
      if (d2 >= r2) continue;
      cnt[i]++; cnt[j]++;
      ax[i] += vx[j]; ay[i] += vy[j]; ax[j] += vx[i]; ay[j] += vy[i];
      cx[i] += dx; cy[i] += dy; cx[j] -= dx; cy[j] -= dy;
      if (d2 < s2 && d2 > 0) { sx[i] -= dx / d2; sy[i] -= dy / d2; sx[j] += dx / d2; sy[j] += dy / d2; }
    }
  }
  const vmax = p.max_speed, vmin = vmax / 2;
  for (let i = 0; i < n; i++) {
    if (cnt[i] > 0) {
      vx[i] += p.separation * K_SEP * sx[i] + p.alignment * K_ALI * (ax[i] / cnt[i] - vx[i]) + p.cohesion * K_COH * (cx[i] / cnt[i]);
      vy[i] += p.separation * K_SEP * sy[i] + p.alignment * K_ALI * (ay[i] / cnt[i] - vy[i]) + p.cohesion * K_COH * (cy[i] / cnt[i]);
    }
    const speed = Math.sqrt(vx[i] * vx[i] + vy[i] * vy[i]);
    if (speed > vmax) { vx[i] = (vx[i] / speed) * vmax; vy[i] = (vy[i] / speed) * vmax; }
    else if (speed < vmin && speed > 0) { vx[i] = (vx[i] / speed) * vmin; vy[i] = (vy[i] / speed) * vmin; }
    x[i] += vx[i]; y[i] += vy[i];
    if (x[i] < 0) x[i] += WORLD; else if (x[i] >= WORLD) x[i] -= WORLD;
    if (y[i] < 0) y[i] += WORLD; else if (y[i] >= WORLD) y[i] -= WORLD;
  }
  b.step++;
}

def step(b, p):
    """One synchronous step: separation (closer than r/2, weighted 1/distance),
    alignment (mean neighbor velocity) and cohesion (mean neighbor offset) within r;
    speed clamped to [max/2, max]; the world wraps around."""
    n, x, y, vx, vy = b["n"], b["x"], b["y"], b["vx"], b["vy"]
    r2 = p["radius"] * p["radius"]
    s2 = r2 / 4
    cnt, ax, ay = [0] * n, [0.0] * n, [0.0] * n
    cx, cy, sx, sy = [0.0] * n, [0.0] * n, [0.0] * n, [0.0] * n
    for i in range(n):
        xi, yi, vxi, vyi = x[i], y[i], vx[i], vy[i]
        for j in range(i + 1, n):
            dx = x[j] - xi
            if dx > HALF:
                dx -= WORLD
            elif dx < -HALF:
                dx += WORLD
            dy = y[j] - yi
            if dy > HALF:
                dy -= WORLD
            elif dy < -HALF:
                dy += WORLD
            d2 = dx * dx + dy * dy
            if d2 >= r2:
                continue
            cnt[i] += 1
            cnt[j] += 1
            ax[i] += vx[j]
            ay[i] += vy[j]
            ax[j] += vxi
            ay[j] += vyi
            cx[i] += dx
            cy[i] += dy
            cx[j] -= dx
            cy[j] -= dy
            if d2 < s2 and d2 > 0:
                sx[i] -= dx / d2
                sy[i] -= dy / d2
                sx[j] += dx / d2
                sy[j] += dy / d2
    vmax = p["max_speed"]
    vmin = vmax / 2
    for i in range(n):
        if cnt[i] > 0:
            vx[i] += p["separation"] * K_SEP * sx[i] + p["alignment"] * K_ALI * (ax[i] / cnt[i] - vx[i]) + p["cohesion"] * K_COH * (cx[i] / cnt[i])
            vy[i] += p["separation"] * K_SEP * sy[i] + p["alignment"] * K_ALI * (ay[i] / cnt[i] - vy[i]) + p["cohesion"] * K_COH * (cy[i] / cnt[i])
        speed = math.sqrt(vx[i] * vx[i] + vy[i] * vy[i])
        if speed > vmax:
            vx[i] = vx[i] / speed * vmax
            vy[i] = vy[i] / speed * vmax
        elif speed < vmin and speed > 0:
            vx[i] = vx[i] / speed * vmin
            vy[i] = vy[i] / speed * vmin
        x[i] += vx[i]
        y[i] += vy[i]
        if x[i] < 0:
            x[i] += WORLD
        elif x[i] >= WORLD:
            x[i] -= WORLD
        if y[i] < 0:
            y[i] += WORLD
        elif y[i] >= WORLD:
            y[i] -= WORLD
    b["step"] += 1

Run the full model from a terminal:

node simulation.js --count 120 --steps 600 --seed 42
python simulation.py --count 120 --steps 600 --seed 42

Download the model

Assumptions & limitations

• Two dimensions, wrapped edges. The world is a 100 × 100 torus: no walls, no obstacles, no ground.
• All-around perception. Boids see every neighbor within r in every direction; there is no blind spot behind them and no limit on how many neighbors they track.
• No noise, no predators, no goals. Once the random start is drawn the dynamics are deterministic.
• Steering gains are modeling choices. The fixed gains and the separation radius r/2 are one reasonable setting; other boids implementations use different forms and give different flocks.
• A behavioral sketch, not animal data. Real flocks and schools show similar patterns, but these rules are not fitted to any species.

Related simulations

• OS-024Schelling segregationAgents with mild preferences for similar neighbors sort themselves into strongly segregated clusters.
• OS-022Game of LifeTwo rules on a grid of cells: gliders, guns, oscillators and soups that settle or keep going.
• OS-013Brownian motionParticles jiggle with Gaussian steps; their mean squared displacement grows as 4Dt, as Einstein predicted.
• 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.

Background: Pseudo-Random Numbers and Seeds in Simulation

References

1. Reynolds, C. W. (1987). Flocks, herds and schools: A distributed behavioral model. ACM SIGGRAPH Computer Graphics, 21(4), 25–34. doi.org/10.1145/37402.37406
2. Wilensky, U. & Rand, W. (2015). An Introduction to Agent-Based Modeling. MIT Press.