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
![\theta_{j}^{t+1}=\arg \left[\sum_{k \sim j} e^{i \theta_{k}^{t}}\right]+\eta \xi_{j}^{t}](https://s0.wp.com/latex.php?latex=%5Ctheta_%7Bj%7D%5E%7Bt%2B1%7D%3D%5Carg+%5Cleft%5B%5Csum_%7Bk+%5Csim+j%7D+e%5E%7Bi+%5Ctheta_%7Bk%7D%5E%7Bt%7D%7D%5Cright%5D%2B%5Ceta+%5Cxi_%7Bj%7D%5E%7Bt%7D&bg=ffffff&fg=333333&s=0&c=20201002)
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,
anim = FuncAnimation(fig,animate,np.arange(1, 200),interval=1, blit=True)
plt.show()
The result is the following animation (the colour indicates the orientation):
Francesco Turci 
Thanks, I think that the FuncAnim needs to be assigned to a variable to work like anim=FuncAnimation(fig,animate,np.arange(1, 200),interval=1, blit=True)
otherwise the animation will not be updated