The Vicsek model is one of the simplest models for active matter. It displays interesting features, such as swarming.

Large scale simulations are often needed in order to provide firm statements on the statistical properties of this kind of models. However, for a pedagogical and illustrative purpose it may be useful to have an elementary code to play with. For this purpose, I have written a relatively simple Python code which implements the model with a few clever tricks to make simulations of a few thousands of agents possible on a standard laptop.

We follow Gregoire and Chaté in the formalism: point-wise particles move synchronously at constant speed v in discretised time of steps Δt=1. The particles have an orientation described by an angle θ which evolves taking into account all particles k within a given radius of interaction

For the neighbourhood calculations, cell-lists would be ideal, but they are too complex for the kind of elementary code that we want to write. What we are going to use is the `kd-tree`

quick nearest neighbour lookup algorithm as implemented in `Scipy`

and some clever sparse matrix manipulation. For visualisation, we employ `matplotlib`

, so that the resulting code is just 60 lines with only very popular libraries.

import numpy as np import scipy as sp from scipy import sparse from scipy.spatial import cKDTree import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation L = 32.0 rho = 3.0 N = int(rho*L**2) print(" N",N) r0 = 1.0 deltat = 1.0 factor =0.5 v0 = r0/deltat*factor iterations = 10000 eta = 0.15 pos = np.random.uniform(0,L,size=(N,2)) orient = np.random.uniform(-np.pi, np.pi,size=N) fig, ax= plt.subplots(figsize=(6,6)) qv = ax.quiver(pos[:,0], pos[:,1], np.cos(orient[0]), np.sin(orient), orient, clim=[-np.pi, np.pi]) def animate(i): print(i) global orient tree = cKDTree(pos,boxsize=[L,L]) dist = tree.sparse_distance_matrix(tree, max_distance=r0,output_type='coo_matrix') #important 3 lines: we evaluate a quantity for every column j data = np.exp(orient[dist.col]*1j) # construct a new sparse marix with entries in the same places ij of the dist matrix neigh = sparse.coo_matrix((data,(dist.row,dist.col)), shape=dist.get_shape()) # and sum along the columns (sum over j) S = np.squeeze(np.asarray(neigh.tocsr().sum(axis=1))) orient = np.angle(S)+eta*np.random.uniform(-np.pi, np.pi, size=N) cos, sin= np.cos(orient), np.sin(orient) pos[:,0] += cos*v0 pos[:,1] += sin*v0 pos[pos>L] -= L pos[pos<0] += L qv.set_offsets(pos) qv.set_UVC(cos, sin,orient) return qv, FuncAnimation(fig,animate,np.arange(1, 200),interval=1, blit=True) plt.show()

The result is the following animation (the colour indicates the orientation):