Vicsek Model

Basic Animation (Particles Moving)

In this chapter I will guide you to build a basic animation of the Vicsek model using Matplotlib. This will be a foundation for more complex visualizations and interactive explorations in later chapters.

You already have the Vicsek functions implemented (e.g. initialize_particles, vicsek_equations, vicsek_order_parameter).
Your job here is to assemble them into a 2D particle animation (with a short tail) that shows the recent trajectory of each particle.

Imports and Model Functions

Create a new script named vicsek_animation.py. At the end of this page you will have a complete version.

Fill in the missing imports:

import sys

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np
from matplotlib.axes import Axes

# make sure Python can find your package/module
sys.path.append(".")

# import Vicsek utilities from your project
from amlab.collective_motion.vicsek import initialize_particles, vicsek_equations, vicsek_order_parameter

Common import issues

If you get ModuleNotFoundError, you likely need sys.path.append(".") (or adjust it to your repo root).
Prefer importing from your package (e.g. amlab...) rather than relative imports.

The Main Function

We will put everything inside a function so you can run it with different dt.

Create a function skeleton and choose parameter defaults.

def run_simulation(dt: float = 1.0) -> None:
    """Run an interactive Vicsek animation (Matplotlib)."""
    # TODO: choose defaults
    num_boids = 
    noise_eta = 
    box_size = 
    radius_interaction = 
    v0 = 

    # Tail length for trajectory visualization
    TAIL_LEN = 20

    # TODO: initialize the state (xy, theta)

    # We stop here for now

Create the Figure Layout

We want a large space for particles, and a smaller area for sliders + extra diagnostics.

Create the subplots and name the axes.

def make_figure():
    # Create a single grid
    fig, axs = plt.subplots(
        1, 1,
        figsize=(12, 8),
        height_ratios=[4, 1],  # top row larger than bottom row
    )
    ax_plane = axs[0]

    return fig, ax_plane

Draw Particles

In order to convey the movement of the particles, we want them to look like “tadpoles”: a big head with a smaller tail. In practice, we will be plotting their current position as the “head”, using a bigger marker, and their recent trajectory as the “tail”, using many small points. This will show the direction of movement and give a better sense of the dynamics. We will use two artists:

one for the tail (many tiny points),
one for the current positions (bigger markers).

Translating this to code, we will - create xy_tail with shape (2, N, TAIL_LEN), and - create the two Line2D artists with ax_plane.plot(...).

def init_particles_plot(ax_plane: Axes, xy: np.ndarray, box_size: float, tail_len: int):
    # build a tail tensor by repeating xy across a third axis
    xy_tail = np.repeat(xy[:, :, np.newaxis], tail_len, axis=2)

    # tail (small grey points)
    (plt_particles,) = ax_plane.plot(
        xy_tail[0].flatten(), xy_tail[1].flatten(),
        linestyle="", marker=".", markersize=2, color="grey"
    )

    # current positions (bigger black circles)
    (plt_current,) = ax_plane.plot(
        xy[0], xy[1],
        linestyle="", marker="o", markersize=3, color="black"
    )

    # set plot limits
    ax_plane.set_xlim(0, box_size)
    ax_plane.set_ylim(0, box_size)
    ax_plane.set_aspect("equal")

    return xy_tail, plt_particles, plt_current

Tail trick (fast and simple)

Store the last TAIL_LEN positions:

shift the tensor with np.roll(...),
write the newest positions into the last slice.

This avoids creating new arrays each frame.

Write the Animation Update Function

This is the core loop:

advance the Vicsek dynamics by one step (vicsek_equations).
update xy_tail.
update the particle artists.

Use this checklist:

call vicsek_equations(...)
update xy_tail with np.roll(...)
plt_particles.set_data(...) and plt_current.set_data(...)

You can follow the structure of the update_animation function below, which includes comments for each step.

def update_animation(frame: int):
    nonlocal xy, xy_tail, theta, noise_eta, v0, radius_interaction, box_size

    # TODO: advance the Vicsek equations

    # TODO: update the tail tensor
    # TODO: update the particle artists with set_data

    return (plt_particles, plt_current)

Tip: blitting requires you to return artists If you set blit=True, your update function must return an iterable of the artists that changed, as shown above. If you forget to return them, the animation will run but the plot will not update (it will look frozen).

Move the function update_animation inside run_simulation so it can access the state variables with nonlocal.

Solution: Simple update_animation()

def update_animation(frame: int):
    nonlocal \
        xy, \
        xy_tail, \
        theta, \
        noise_eta, \
        v0, \
        radius_interaction, \
        box_size
    xy, theta = vicsek_equations(
        xy,
        theta,
        v0=v0,
        dt=dt,
        radius_interaction=radius_interaction,
        box_size=box_size,
        noise=noise_eta,
    )

    # Update tails
    xy_tail = np.roll(xy_tail, shift=-1, axis=2)
    xy_tail[:, :, -1] = xy
    plt_particles.set_data(xy_tail[0].flatten(), xy_tail[1].flatten())
    plt_current.set_data(xy[0], xy[1])

    return (plt_particles, plt_current)

Create the Animation Object

Matplotlib uses FuncAnimation. Inside run_simulation, after defining update_animation, create the animation object:

ani = animation.FuncAnimation(fig, update_animation, interval=0, blit=True)

interval=0 makes it as fast as your computer can handle.
You can set a larger interval (e.g. 20 ms) to reduce CPU usage.

Full Script

At this point you have built every piece.
Below is a sample of how your run_simulation function might look with all the pieces assembled.

Solution: Simplified run_simulation()

This is a simplified version that assumes you have not implemented the order plot or noise plot, neither the sliders, yet. It focuses on the particle animation.

def run_simulation(dt: float = 1.0) -> None:
    # Parameters
    num_boids = 100
    noise_eta = 0.5
    box_size = 10.0
    radius_interaction = 1.0
    v0 = 0.03

    # Initialize state
    xy, theta = initialize_particles(num_boids, box_size)

    # Order parameter history and noise diagnostic
    # will be initialized here

    # Create figure and axes
    fig, ax_plane, ax_order, ax_sliders, ax_noise = make_figure()

    # Initialize plots
    xy_tail, plt_particles, plt_current = init_particles_plot(ax_plane, xy, box_size, tail_len=20)
    # The order plot and noise plot will go here

    # Sliders will be added here

    # Define the animation update function
    def update_animation(frame: int):
        # See the previous step for the full implementation
        ...

    # Create the animation object
    ani = animation.FuncAnimation(fig, update_animation, interval=0, blit=True)

    # The slider callbacks will come here

    plt.show()

Important: For your script to work, remember to call the run_simulation() function at the end of your script:

if __name__ == "__main__":
    run_simulation()

Checklist (Debugging)

If your animation is blank:

Did you return the artists in update_animation?
Did you set axis limits (set_xlim, set_ylim)?
Are xy and theta updated each frame?
Are you calling plt.show() at the end of run_simulation()?
Are you calling run_simulation() at the end of your script?

Extension

Add the order parameter plot and noise diagnostic plot.
Add sliders for noise_eta and v0.
Add a slider for radius_interaction and see how it affects the dynamics.

Go to the next chapter to see how to add the order parameter plot and noise diagnostic plot.

Next: Adding Diagnostics and Sliders

--- title: "Vicsek Model" subtitle: "Basic Animation (Particles Moving)" format: html --- In this chapter I will guide you to build a basic animation of the Vicsek model using Matplotlib. This will be a foundation for more complex visualizations and interactive explorations in later chapters. You already have the Vicsek functions implemented (e.g. `initialize_particles`, `vicsek_equations`, `vicsek_order_parameter`). Your job here is to **assemble** them into a 2D particle animation (with a short *tail*) that shows the recent trajectory of each particle. --- ## Imports and Model Functions Create a new script named `vicsek_animation.py`. At the end of this page you will have a complete version. Fill in the missing imports: ```python import sys import matplotlib.pyplot as plt import matplotlib.animation as animation import numpy as np from matplotlib.axes import Axes # make sure Python can find your package/module sys.path.append(".") # import Vicsek utilities from your project from amlab.collective_motion.vicsek import initialize_particles, vicsek_equations, vicsek_order_parameter ``` ::: {.callout-tip collapse="true"} ## Common import issues - If you get `ModuleNotFoundError`, you likely need `sys.path.append(".")` (or adjust it to your repo root). - Prefer importing from your package (e.g. `amlab...`) rather than relative imports. ::: --- ## The Main Function We will put everything inside a function so you can run it with different `dt`. Create a function skeleton and choose parameter defaults. ```python def run_simulation(dt: float = 1.0) -> None: """Run an interactive Vicsek animation (Matplotlib).""" # TODO: choose defaults num_boids = noise_eta = box_size = radius_interaction = v0 = # Tail length for trajectory visualization TAIL_LEN = 20 # TODO: initialize the state (xy, theta) # We stop here for now ``` --- ## Create the Figure Layout We want a large space for particles, and a smaller area for sliders + extra diagnostics. Create the subplots and name the axes. ```python def make_figure(): # Create a single grid fig, axs = plt.subplots( 1, 1, figsize=(12, 8), height_ratios=[4, 1], # top row larger than bottom row ) ax_plane = axs[0] return fig, ax_plane ``` --- ## Draw Particles In order to convey the movement of the particles, we want them to look like "tadpoles": a big head with a smaller tail. In practice, we will be plotting their current position as the "head", using a bigger marker, and their recent trajectory as the "tail", using many small points. This will show the direction of movement and give a better sense of the dynamics. We will use **two** artists: - one for the **tail** (many tiny points), - one for the **current** positions (bigger markers). Translating this to code, we will - create `xy_tail` with shape `(2, N, TAIL_LEN)`, and - create the two `Line2D` artists with `ax_plane.plot(...)`. ```python def init_particles_plot(ax_plane: Axes, xy: np.ndarray, box_size: float, tail_len: int): # build a tail tensor by repeating xy across a third axis xy_tail = np.repeat(xy[:, :, np.newaxis], tail_len, axis=2) # tail (small grey points) (plt_particles,) = ax_plane.plot( xy_tail[0].flatten(), xy_tail[1].flatten(), linestyle="", marker=".", markersize=2, color="grey" ) # current positions (bigger black circles) (plt_current,) = ax_plane.plot( xy[0], xy[1], linestyle="", marker="o", markersize=3, color="black" ) # set plot limits ax_plane.set_xlim(0, box_size) ax_plane.set_ylim(0, box_size) ax_plane.set_aspect("equal") return xy_tail, plt_particles, plt_current ``` ::: {.callout-tip collapse="true"} ## Tail trick (fast and simple) Store the last `TAIL_LEN` positions: - shift the tensor with `np.roll(...)`, - write the newest positions into the last slice. This avoids creating new arrays each frame. ::: --- ## Write the Animation Update Function This is the core loop: 1. advance the Vicsek dynamics by one step (`vicsek_equations`). 2. update `xy_tail`. 3. update the particle artists. Use this checklist: - [ ] call `vicsek_equations(...)` - [ ] update `xy_tail` with `np.roll(...)` - [ ] `plt_particles.set_data(...)` and `plt_current.set_data(...)` You can follow the structure of the `update_animation` function below, which includes comments for each step. ```python def update_animation(frame: int): nonlocal xy, xy_tail, theta, noise_eta, v0, radius_interaction, box_size # TODO: advance the Vicsek equations # TODO: update the tail tensor # TODO: update the particle artists with set_data return (plt_particles, plt_current) ``` **Tip:** blitting requires you to return artists If you set `blit=True`, your update function **must** return an iterable of the artists that changed, as shown above. If you forget to return them, the animation will run but the plot will not update (it will look frozen). Move the function `update_animation` inside `run_simulation` so it can access the state variables with `nonlocal`. ::: {.callout-tip collapse="true"} ## Solution: Simple `update_animation()` ```python def update_animation(frame: int): nonlocal \ xy, \ xy_tail, \ theta, \ noise_eta, \ v0, \ radius_interaction, \ box_size xy, theta = vicsek_equations( xy, theta, v0=v0, dt=dt, radius_interaction=radius_interaction, box_size=box_size, noise=noise_eta, ) # Update tails xy_tail = np.roll(xy_tail, shift=-1, axis=2) xy_tail[:, :, -1] = xy plt_particles.set_data(xy_tail[0].flatten(), xy_tail[1].flatten()) plt_current.set_data(xy[0], xy[1]) return (plt_particles, plt_current) ``` ::: --- ## Create the Animation Object Matplotlib uses `FuncAnimation`. Inside `run_simulation`, after defining `update_animation`, create the animation object: ```python ani = animation.FuncAnimation(fig, update_animation, interval=0, blit=True) ``` - `interval=0` makes it as fast as your computer can handle. - You can set a larger interval (e.g. 20 ms) to reduce CPU usage. --- ## Full Script At this point you have built every piece. Below is a sample of how your `run_simulation` function might look with all the pieces assembled. ::: {.callout-tip collapse="true"} ## Solution: Simplified `run_simulation()` This is a simplified version that assumes you have not implemented the order plot or noise plot, neither the sliders, yet. It focuses on the particle animation. ```python def run_simulation(dt: float = 1.0) -> None: # Parameters num_boids = 100 noise_eta = 0.5 box_size = 10.0 radius_interaction = 1.0 v0 = 0.03 # Initialize state xy, theta = initialize_particles(num_boids, box_size) # Order parameter history and noise diagnostic # will be initialized here # Create figure and axes fig, ax_plane, ax_order, ax_sliders, ax_noise = make_figure() # Initialize plots xy_tail, plt_particles, plt_current = init_particles_plot(ax_plane, xy, box_size, tail_len=20) # The order plot and noise plot will go here # Sliders will be added here # Define the animation update function def update_animation(frame: int): # See the previous step for the full implementation ... # Create the animation object ani = animation.FuncAnimation(fig, update_animation, interval=0, blit=True) # The slider callbacks will come here plt.show() ``` ::: **Important:** For your script to work, remember to call the `run_simulation()` function at the end of your script: ```python if __name__ == "__main__": run_simulation() ``` --- ## Checklist (Debugging) If your animation is blank: - Did you return the artists in `update_animation`? - Did you set axis limits (`set_xlim`, `set_ylim`)? - Are `xy` and `theta` updated each frame? - Are you calling `plt.show()` at the end of `run_simulation()`? - Are you calling `run_simulation()` at the end of your script? --- ## Extension - Add the order parameter plot and noise diagnostic plot. - Add sliders for `noise_eta` and `v0`. - Add a slider for `radius_interaction` and see how it affects the dynamics. Go to the next chapter to see how to add the order parameter plot and noise diagnostic plot. [Next: Adding Diagnostics and Sliders](vicsek-animation-interactive.qmd){.btn .btn-primary}