Vicsek Lab: Initialization

Module 2: Synchronization and Collective Behaviour

In this lab you implement the Vicsek model (Vicsek et al. 1995), validate it, and animate it. The rules are in The Vicsek Model.

Core activity

The core of this lab is five pages: this one, Update Rule, Order Parameter, Validation and Animation. Vectorized Update, Interactive Animation and Predator Avoidance are optional extensions. Practice: Collective Motion closes the session.

Why not solve_ivp?

In Module 1 you used ODE solvers like scipy.integrate.solve_ivp() for continuous-time systems. The Vicsek model is a discrete-time, agent-based model. At each step all particles update their heading and position from their neighbours and some noise. There is no differential equation to integrate, so you implement it with a loop over time steps.

Prediction

Exercise

Before you code, write down your expectations.

  • With noise \(\eta = 0\) and a large interaction radius, what is the order parameter after many steps?
  • With \(\eta = 2\pi\), what is it for \(N = 100\)?
  • At fixed noise, does increasing the density help or hurt alignment?

With no noise and every particle seeing every other, all headings converge to one direction and \(v_a \to 1\). With \(\eta = 2\pi\) every heading is random at every step, so \(v_a\) is of order \(1/\sqrt{N} = 0.1\). Higher density means more neighbours inside the radius, which averages out the noise and helps alignment (Vicsek et al. 1995).

Initial State

Exercise

Create the initial positions and headings. Positions are uniform on the box \([0, L]^2\) and headings are uniform on \([0, 2\pi)\).

def initialize_particles(num_boids, box_size):
    # Random initial theta, shape (N)
    theta = # TODO

    # Random initial x, y - shape (2, N)
    xy = # TODO
    return xy, theta

np.random.uniform(low, high, size) accepts a tuple as size. Use size=(2, num_boids) for the positions, so that xy[0] holds the x-coordinates and xy[1] the y-coordinates.

def initialize_particles(num_boids, box_size):
    theta = np.random.uniform(0, 2 * np.pi, num_boids)
    xy = np.random.uniform(0, box_size, (2, num_boids))
    return xy, theta

Try out your function now. Vicsek et al. (1995) used between \(N = 40\) and \(N = 10000\) particles and box sizes between \(L = 3\) and \(L = 50\). The more particles you have, the slower the simulation will be, so start with a small number like \(N=100\) and a box size of \(L=10\).

Figure 1 shows what you are building: after a few steps, nearby particles already point in similar directions.

Figure 1: Vicsek model after 20 steps: \(N = 100\) particles in a box of side \(L = 10\), noise \(\eta = 0.5\). Arrows show headings and grey dots the recent trajectory.

References

Vicsek, Tamás, András Czirók, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet. 1995. “Novel Type of Phase Transition in a System of Self-Driven Particles.” Physical Review Letters 75 (6): 1226.