Vicsek Lab: Update Rule

Module 2: Synchronization and Collective Behaviour

Implementation: the update rule

All boids have the same speed module $v_0$, but their headings change based on neighbors. Every time step $\Delta t$, each boid averages the headings of neighbors within radius $r$ and adds noise $\eta$. Then it moves in the new heading direction. This happens simultaneously to all boids. Follow the equations in Vicsek et al. (1995) to implement the update step. A summary follows.

The positions of the $i$th particle is updated as:

\[ \mathbf{x}_i(t + \Delta t) = \mathbf{x}_i(t) + v_0 \Delta t \begin{bmatrix} \cos(\theta_i(t + \Delta t)) \\ \sin(\theta_i(t + \Delta t)) \end{bmatrix}. \tag{1}\]

where $v_0$ is the speed. Vicsek et al. (1995) used $v_0 = 0.03$.

It is important to remember that particles stay within the box. You apply periodic boundary conditions, meaning that a particle that goes out of the box on one side reappears on the opposite side. You can implement this with the modulo operator, if you assume that the positions are in the range $[0, L]$:

\[ \mathbf{x}_i(t + \Delta t) = \mathbf{x}_i(t + \Delta t) \mod L, \tag{2}\]

To use Equation 1 you first need the new heading $\theta_i(t + \Delta t)$:

\[ \theta_i(t + \Delta t) = \langle \theta_j(t) \rangle_{j \in \mathcal{N}_i} + \frac{\eta}{2} \xi_i(t), \tag{3}\]

where $\mathcal{N}_i$ is the set of neighbors within radius $r$ of particle $i$, including $i$ itself, $\langle \cdot \rangle$ denotes the average direction, computed as $\operatorname{atan2}(\langle\sin\theta\rangle, \langle\cos\theta\rangle)$, and $\xi_i(t)$ is a random variable uniformly distributed in $[-1, 1]$. Vicsek et al. (1995) used $r = 1$ and noise values $\eta$ between $0$ and $5$ in their Figure 2.

Steps

Outline the steps for one time step before you code:

For each particle $i$, find its neighbors within radius $r$.
Compute the average heading of these neighbors.
Add noise to the average heading to get the new heading $\theta_i(t + \Delta t)$.
Update the position of each particle based on its new heading.
Apply periodic boundary conditions to ensure particles stay within the box.

For the first step you compare every particle with every other particle. You can use nested loops or a distance matrix. Start with the loop version. Once you know the neighbors, average their headings as unit vectors, then add noise following Equation 3.

Do not average angles directly

The plain mean of $0.1$ and $2\pi - 0.1$ is $\pi$, the opposite direction. Average $\cos\theta_j$ and $\sin\theta_j$ separately and recover the angle with np.arctan2.

Exercise

Complete the loop that computes the new heading of every particle. Use the variable names in the template: radius_interaction is $r$ and noise is $\eta$.

# What we have:
xy # the positions of the particles, shape (2, N)
theta # the headings of the particles, shape (N)
noise # the noise level eta
# ----

num_boids = xy.shape[1]  # N

# Loop over each particle
for i in range(num_boids):
    # Initialize list to store neighbor headings
    ls_theta_neighbors = []
    # Loop over all other particles to find neighbors
    for j in range(num_boids):
        # TODO: Compute distance between particle i and j
        # (j == i gives distance 0, so i counts as its own neighbor)

        # TODO: If distance is less than radius,
        # add the heading of particle j to the list of neighbors
    # After finding neighbors...
    # TODO: Compute the average direction with np.arctan2 and add noise

# Result: updated theta for all particles

Hint

np.linalg.norm(xy[:, i] - xy[:, j]) is the distance between particles i and j. After the inner loop, convert the list to an array and use np.arctan2(np.mean(np.sin(...)), np.mean(np.cos(...))). For the noise, noise / 2 * np.random.uniform(-1, 1) is uniform on $[-\eta/2, \eta/2]$.

Solution

num_boids = xy.shape[1]  # N
theta_new = theta.copy()  # new headings, filled in below

# Loop over each particle
for i in range(num_boids):
    # Initialize list to store neighbor headings
    ls_theta_neighbors = []
    # Loop over all other particles to find neighbors
    for j in range(num_boids):
        # Compute distance between particle i and j (i itself is included)
        dist = np.linalg.norm(xy[:, i] - xy[:, j])
        # If distance is less than radius,
        # add the heading of particle j to the list of neighbors
        if dist <= radius_interaction:
            ls_theta_neighbors.append(theta[j])
    # After finding neighbors
    # Compute the average direction...
    ls_theta_neighbors = np.array(ls_theta_neighbors)
    theta_new[i] = np.arctan2(
        np.mean(np.sin(ls_theta_neighbors)), np.mean(np.cos(ls_theta_neighbors))
    )
    # ...and add noise, uniform in [-eta/2, eta/2]
    theta_new[i] += noise / 2 * np.random.uniform(-1, 1)

This covers the neighbor search, the average and the noise. Now compute the new positions from the new headings, following Equation 1, and apply the periodic boundary conditions of Equation 2.

Exercise

Complete the two templates below. v0 is the speed $v_0$, dt the time step and box_size the side $L$.

# After updating theta for all particles
for i in range(num_boids):
    # TODO: Compute velocity vector based on new heading

    # TODO: Update position based on velocity and time step

# Result: updated xy for all particles

# After updating positions for all particles
# TODO: Apply periodic boundary conditions to xy

# Result: updated xy with periodic boundaries

Solution

# After updating theta for all particles
xy_new = xy.copy()
for i in range(num_boids):
    # Compute velocity vector based on new heading
    v = v0 * np.array([np.cos(theta_new[i]), np.sin(theta_new[i])])
    # Update position based on velocity and time step
    xy_new[:, i] = xy_new[:, i] + dt * v

# Apply periodic boundary conditions to xy
xy = np.mod(xy_new, box_size)
# We overwrite xy with the new positions after applying periodic boundaries

Exercise

Put everything together in a function that takes the current positions and headings and returns them after one time step.

def vicsek_equations(xy, theta, dt=1.0, eta=0.1, box_size=10, radius_interaction=1, v0=0.03):
    xy_new = xy.copy()
    theta_new = theta.copy()
    # TODO: new headings (loop over particles and neighbors)

    # TODO: new positions

    # TODO: periodic boundary conditions
    return xy_new, theta_new

Solution

def vicsek_equations(xy, theta, dt=1.0, eta=0.1, box_size=10, radius_interaction=1, v0=0.03):
    xy_new = xy.copy()
    theta_new = theta.copy()
    num_boids = len(theta)
    for i in range(num_boids):
        ls_theta_neighbors = []
        for j in range(num_boids):
            dist = np.linalg.norm(xy[:, i] - xy[:, j])
            if dist <= radius_interaction:
                ls_theta_neighbors.append(theta[j])
        ls_theta_neighbors = np.array(ls_theta_neighbors)
        theta_new[i] = np.arctan2(
            np.mean(np.sin(ls_theta_neighbors)), np.mean(np.cos(ls_theta_neighbors))
        )
        theta_new[i] += eta / 2 * np.random.uniform(-1, 1)

        v = v0 * np.array([np.cos(theta_new[i]), np.sin(theta_new[i])])
        xy_new[:, i] = xy_new[:, i] + dt * v

    xy_new = np.mod(xy_new, box_size)
    return xy_new, theta_new

The optional page Vectorized Update removes the loops and discusses a limitation of the distance computation.

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.

--- title: "Vicsek Lab: Update Rule" subtitle: "Module 2: Synchronization and Collective Behaviour" format: html --- ## Implementation: the update rule All boids have the same speed module $v_0$, but their headings change based on neighbors. Every time step $\Delta t$, each boid averages the headings of neighbors within radius $r$ and adds noise $\eta$. Then it moves in the new heading direction. This happens **simultaneously** to all boids. Follow the equations in @vicsek1995novel to implement the update step. A summary follows. The positions of the $i$th particle is updated as: $$ \mathbf{x}_i(t + \Delta t) = \mathbf{x}_i(t) + v_0 \Delta t \begin{bmatrix} \cos(\theta_i(t + \Delta t)) \\ \sin(\theta_i(t + \Delta t)) \end{bmatrix}. $$ {#eq-vicsek-xi} where $v_0$ is the speed. @vicsek1995novel used $v_0 = 0.03$. It is important to remember that particles stay within the box. You apply **periodic boundary conditions**, meaning that a particle that goes out of the box on one side reappears on the opposite side. You can implement this with the modulo operator, if you assume that the positions are in the range $[0, L]$: $$ \mathbf{x}_i(t + \Delta t) = \mathbf{x}_i(t + \Delta t) \mod L, $$ {#eq-vicsek-periodic} To use @eq-vicsek-xi you first need the new heading $\theta_i(t + \Delta t)$: $$ \theta_i(t + \Delta t) = \langle \theta_j(t) \rangle_{j \in \mathcal{N}_i} + \frac{\eta}{2} \xi_i(t), $$ {#eq-vicsek-theta} where $\mathcal{N}_i$ is the set of neighbors within radius $r$ of particle $i$, including $i$ itself, $\langle \cdot \rangle$ denotes the average direction, computed as $\operatorname{atan2}(\langle\sin\theta\rangle, \langle\cos\theta\rangle)$, and $\xi_i(t)$ is a random variable uniformly distributed in $[-1, 1]$. @vicsek1995novel used $r = 1$ and noise values $\eta$ between $0$ and $5$ in their Figure 2. ### Steps Outline the steps for one time step before you code: 1. For each particle $i$, find its neighbors within radius $r$. 2. Compute the average heading of these neighbors. 3. Add noise to the average heading to get the new heading $\theta_i(t + \Delta t)$. 4. Update the position of each particle based on its new heading. 5. Apply periodic boundary conditions to ensure particles stay within the box. For the first step you compare every particle with every other particle. You can use nested loops or a distance matrix. Start with the loop version. Once you know the neighbors, average their headings as unit vectors, then add noise following @eq-vicsek-theta. ::: {.callout-warning} ## Do not average angles directly The plain mean of $0.1$ and $2\pi - 0.1$ is $\pi$, the opposite direction. Average $\cos\theta_j$ and $\sin\theta_j$ separately and recover the angle with `np.arctan2`. ::: ::: {.callout-note} ## Exercise Complete the loop that computes the new heading of every particle. Use the variable names in the template: `radius_interaction` is $r$ and `noise` is $\eta$. ```python # What we have: xy # the positions of the particles, shape (2, N) theta # the headings of the particles, shape (N) noise # the noise level eta # ---- num_boids = xy.shape[1] # N # Loop over each particle for i in range(num_boids): # Initialize list to store neighbor headings ls_theta_neighbors = [] # Loop over all other particles to find neighbors for j in range(num_boids): # TODO: Compute distance between particle i and j # (j == i gives distance 0, so i counts as its own neighbor) # TODO: If distance is less than radius, # add the heading of particle j to the list of neighbors # After finding neighbors... # TODO: Compute the average direction with np.arctan2 and add noise # Result: updated theta for all particles ``` ::: ::: {.callout-tip collapse="true"} ## Hint `np.linalg.norm(xy[:, i] - xy[:, j])` is the distance between particles `i` and `j`. After the inner loop, convert the list to an array and use `np.arctan2(np.mean(np.sin(...)), np.mean(np.cos(...)))`. For the noise, `noise / 2 * np.random.uniform(-1, 1)` is uniform on $[-\eta/2, \eta/2]$. ::: ::: {.callout-tip collapse="true"} ## Solution ```python num_boids = xy.shape[1] # N theta_new = theta.copy() # new headings, filled in below # Loop over each particle for i in range(num_boids): # Initialize list to store neighbor headings ls_theta_neighbors = [] # Loop over all other particles to find neighbors for j in range(num_boids): # Compute distance between particle i and j (i itself is included) dist = np.linalg.norm(xy[:, i] - xy[:, j]) # If distance is less than radius, # add the heading of particle j to the list of neighbors if dist <= radius_interaction: ls_theta_neighbors.append(theta[j]) # After finding neighbors # Compute the average direction... ls_theta_neighbors = np.array(ls_theta_neighbors) theta_new[i] = np.arctan2( np.mean(np.sin(ls_theta_neighbors)), np.mean(np.cos(ls_theta_neighbors)) ) # ...and add noise, uniform in [-eta/2, eta/2] theta_new[i] += noise / 2 * np.random.uniform(-1, 1) ``` ::: This covers the neighbor search, the average and the noise. Now compute the new positions from the new headings, following @eq-vicsek-xi, and apply the periodic boundary conditions of @eq-vicsek-periodic. ::: {.callout-note} ## Exercise Complete the two templates below. `v0` is the speed $v_0$, `dt` the time step and `box_size` the side $L$. ```python # After updating theta for all particles for i in range(num_boids): # TODO: Compute velocity vector based on new heading # TODO: Update position based on velocity and time step # Result: updated xy for all particles ``` ```python # After updating positions for all particles # TODO: Apply periodic boundary conditions to xy # Result: updated xy with periodic boundaries ``` ::: ::: {.callout-tip collapse="true"} ## Solution ```python # After updating theta for all particles xy_new = xy.copy() for i in range(num_boids): # Compute velocity vector based on new heading v = v0 * np.array([np.cos(theta_new[i]), np.sin(theta_new[i])]) # Update position based on velocity and time step xy_new[:, i] = xy_new[:, i] + dt * v # Apply periodic boundary conditions to xy xy = np.mod(xy_new, box_size) # We overwrite xy with the new positions after applying periodic boundaries ``` ::: ::: {.callout-note} ## Exercise Put everything together in a function that takes the current positions and headings and returns them after one time step. ```python def vicsek_equations(xy, theta, dt=1.0, eta=0.1, box_size=10, radius_interaction=1, v0=0.03): xy_new = xy.copy() theta_new = theta.copy() # TODO: new headings (loop over particles and neighbors) # TODO: new positions # TODO: periodic boundary conditions return xy_new, theta_new ``` ::: ::: {.callout-tip collapse="true"} ## Solution ```python def vicsek_equations(xy, theta, dt=1.0, eta=0.1, box_size=10, radius_interaction=1, v0=0.03): xy_new = xy.copy() theta_new = theta.copy() num_boids = len(theta) for i in range(num_boids): ls_theta_neighbors = [] for j in range(num_boids): dist = np.linalg.norm(xy[:, i] - xy[:, j]) if dist <= radius_interaction: ls_theta_neighbors.append(theta[j]) ls_theta_neighbors = np.array(ls_theta_neighbors) theta_new[i] = np.arctan2( np.mean(np.sin(ls_theta_neighbors)), np.mean(np.cos(ls_theta_neighbors)) ) theta_new[i] += eta / 2 * np.random.uniform(-1, 1) v = v0 * np.array([np.cos(theta_new[i]), np.sin(theta_new[i])]) xy_new[:, i] = xy_new[:, i] + dt * v xy_new = np.mod(xy_new, box_size) return xy_new, theta_new ``` ::: The optional page [Vectorized Update](vicsek-vectorized.qmd) removes the loops and discusses a limitation of the distance computation.