What could be a good potential function for flocks of birds?

I am interested in modelling flocks of birds but I have difficulties finding a good candidate for my potential function $$V(\vec{r})$$. It would need to have the following characteristics:

• A short range repulsion (birds do not want to crash into their neighbours)
• A long range attraction (birds aim to stay in packs)
• (From wikipedia: "Alignment - steer towards average heading of neighbours")

What function could be a good candidate for this?

In general, for a $$N$$-body system:

In the case of gravity, we may use the potential $$V_g(\mathbf{x})=\sum_{i=1}^{n}-\frac{G m_{i}}{\left|\mathbf{x}-\mathbf{x}_{\mathbf{i}}\right|}$$

In the case of flocking, I tried using Lennard-Jones potential: $$V_{\mathrm{LJ}}=4 \varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]$$ since it satisfies the short range repulsion and the long range attraction, but the resulting modelling does not look like bird flocking at all...

Any ideas/advice are welcome and much appreciated. I can share screen shots showing the results.

I don't know how you did your simulation, but you should at least include some kind of damping term in the motion of the birds, otherwise every configuration starting out-of-equilibrium will keep oscillating forever. Maybe add a term like $$\overrightarrow{F_{\mathrm{viscous}}} = - \alpha \overrightarrow{v}$$. If you are not only interested in the equilibrium position but also the dynamics of the flock for small perturbations, you can also add a random Langevin force $$\overrightarrow{F_{\mathrm{Langevin}}}(t) = \eta(t)$$, with $$\left\langle \eta(t) \eta(t') \right\rangle = A \delta(t-t')$$.

For the exact form of the interaction force between the birds, I don't know if it is possible to capture all of the complexity of the phenomenon simply by adding pairwise, conservative force to the problem, but we can play around and see if the results are good enough in the end.

If you want a "homogeneous" density of birds, you can try $$\overrightarrow{F}_{\overrightarrow{x}\to\overrightarrow{x'}}= \frac{k}{N} \left(\overrightarrow{x} - \overrightarrow{x'}\right) - \beta \frac{\left(\overrightarrow{x} - \overrightarrow{x'}\right)}{\left|\overrightarrow{x} - \overrightarrow{x'} \right|^2},$$

where $$\overrightarrow{F}_{\overrightarrow{x}\to\overrightarrow{x'}}$$ is the force exerted by a bird at position $$\overrightarrow{x}$$ on a bird at position $$\overrightarrow{x'}$$, $$N$$ is the number of birds in the flock (a normalization factor), and $$k$$ and $$\beta$$ are positive constants whose relative value will change the average spacing between birds at equilibrium. You also need to fix the relative value of the Langevin force and viscous force to somewhat realistic values so that you reach an equilibrium fast enough but you can still see the small perturbations in position. I will come with a more detailed explanation for the form of the force and a simulation later.

Reason for the choice of the force:

The first term in $$\overrightarrow{F}_{\overrightarrow{x}\to\overrightarrow{x'}}$$, $$\frac{k}{N} \left(\overrightarrow{x} - \overrightarrow{x'}\right)$$, simply correspond to the force of an harmonic oscillator of strengh $$k$$ centered on the center of mass of the flock. So, assuming the center of mass is centered around $$0$$, the force is simply $$\overrightarrow{F}(\overrightarrow{x}) = -k \overrightarrow{x}$$.

The reason for the second term is actually a little bit subtle. By analogy with the 3D Coulomb force, you can show that the field $$\overrightarrow{f}(\overrightarrow{x})$$ created by a density of birds $$\rho(\overrightarrow{x}')$$ follows a 2D analog to Gauss law in electromagnetism, which is $$\overrightarrow{\nabla} \cdot \overrightarrow{f} = 2 \pi \beta \rho$$. You can then show that a constant density of birds $$\rho(r) = N/(\pi R^2)$$ for $$0 < r < R$$, and $$\rho(r>R) = 0$$ is creating a repulsive force of the form $$\overrightarrow{f}(\overrightarrow{x}) = \beta N \overrightarrow{x}/R^2$$ (for $$\left| \overrightarrow{x} \right| < R$$), which exactly compensates the harmonic oscillator force when $$R = \sqrt{\beta N/k}$$. This is of course only true for a continuous "density" of birds, which is not the case in reality but for $$N \gg 1$$, you can assimilate $$\rho$$ to the average density of birds. In short, for $$N \gg 1$$, the birds at equilibrium should arrange themselves in a somewhat homogeneous configuration within a radius $$R = \sqrt{\beta N/k}$$.

Simulation

I implemented the problem in Python using all the forces that I described earlier, with $$k = \beta = 1$$, $$A = 3$$ and $$\alpha = 2$$. I integrate in time steps of $$dt = 0.01$$ from $$t = 0$$ to $$t = T = 10$$. The mass $$m$$ of each bird is also set to $$1$$. Initially, the position of each bird is drawn randomly and uniformly in a square smaller than the expected radius, so that you can see the out-of-equilibrium dynamics.

Here is what the trajectories look like on a typical run

And here is a nice gif of the time evolution of the bird flock:

As you can see, despite some random fluctuations in position, the birds keep moving around an "equilibrium" configuration where the number of birds per surface area is approximately constant within a given radius.

If you're worried that the birds get too close to each other and it is unrealistic, try decreasing the $$A$$ parameter (controlling the amplitude of the random motion of birds). See below a gif for $$A = 0.5$$.

Code used

import numpy as np
import matplotlib.pyplot as plt

plt.rc('font', size=20)

N = 160 #number of birds

sizex = 20 #size of the image for plotting
sizey = 20

k = 1.0 #strength of the (long-range) attractive force
beta = 1.0 #strength of the (short-range) repulsive force
alpha = 2.0 #"viscosity"
A = 3.0 #intensity of the random erratic movement of each bird

dt = 0.01 #timestep
T = 10 #final time

expected_radius = (beta*N/k)**0.5 #expect radius of the bird flock by analogy with the continuous case

X_pos_0 = (np.random.random(N)-0.5)*0.5*sizex #draw the initial X position of birds randomly
X_pos_0 -= np.mean(X_pos_0) #center it on 0
Y_pos_0 = (np.random.random(N)-0.5)*0.5*sizey #same in the Y direction
Y_pos_0 -= np.mean(Y_pos_0)

X_pos = np.zeros((N, int(T/dt)+1)) #array of X position of each bird
Y_pos = np.zeros((N, int(T/dt)+1)) #array of Y position of each bird
X_speed = np.zeros((N, int(T/dt)+1)) #array of X velocity of each bird
Y_speed = np.zeros((N, int(T/dt)+1)) #array of Y velocity of each bird

X_pos[:,0], Y_pos[:,0] = X_pos_0, Y_pos_0

@np.vectorize
def repulsive_force_x(X,Y): #function used to compute the pairwise forces at each timestep of the algorithm
if X==0 and Y==0:
return 0.0
else:
return -beta*X/(X**2+Y**2)

@np.vectorize
def repulsive_force_y(X,Y):
if X==0 and Y==0:
return 0.0
else:
return -beta*Y/(X**2+Y**2)

time = np.linspace(0, T, int(T/dt)+1)
for i in range(int(T/dt)):
langevin_force_x = np.random.normal(0, A, N)/dt**0.5 #the Langevin force has to be normalized like this in the limit dt -> 0
langevin_force_y = np.random.normal(0, A, N)/dt**0.5
viscous_force_x = -alpha*X_speed[:, i]
viscous_force_y = -alpha*Y_speed[:, i]
attractive_force_x = -k*(X_pos[:,i]-np.mean(X_pos[:,i]))
attractive_force_y = -k*(Y_pos[:,i]-np.mean(Y_pos[:,i]))
X1, X2 = np.meshgrid(X_pos[:,i], X_pos[:,i])
distance_matrix_X = X1-X2 #array containing all pairwise distances x-x' along X
Y1, Y2 = np.meshgrid(Y_pos[:,i], Y_pos[:,i])
distance_matrix_Y = Y1-Y2 #array containing all pairwise distances y-y' along Y
force_X = np.sum(repulsive_force_x(distance_matrix_X, distance_matrix_Y), axis=1)
force_Y = np.sum(repulsive_force_y(distance_matrix_X, distance_matrix_Y), axis=1)
X_speed[:, i+1] = X_speed[:, i] + dt*(langevin_force_x + viscous_force_x + attractive_force_x + force_X) #v(t+dt) = v(t) + a(t)*dt (m=1)
Y_speed[:, i+1] = Y_speed[:, i] + dt*(langevin_force_y + viscous_force_y + attractive_force_y + force_Y)
X_pos[:, i+1] = X_pos[:, i] + dt*X_speed[:, i+1] #x(t+dt) = x(t) + v(t+dt)*dt (symplectic integration)
Y_pos[:, i+1] = Y_pos[:, i] + dt*Y_speed[:, i+1]
if i%(int(T/dt)//10) ==0: #progression of the calculation
print(10*i//(int(T/dt)//10), '%')

cmap_cool = plt.get_cmap('cool')

#plot initial and final positions only
fig, ax = plt.subplots(1, figsize=(15,12))
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(X_pos[:, 0], Y_pos[:, 0], 'o', label="Initial position", color=cmap_cool(0)) #initial position in blue
ax.plot(X_pos[:, -1], Y_pos[:, -1], 'o', label="Final position", color=cmap_cool(0.99)) #final position in purple
ax.axis("scaled")
ax.set_xlim(left=-24, right=24)
plt.legend(loc='best')

#plot trajectories too
fig, ax = plt.subplots(1, figsize=(15,12))
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(X_pos[:, 0], Y_pos[:, 0], 'o', label="Initial position", color=cmap_cool(0)) #initial position in blue
ax.plot(X_pos[:, -1], Y_pos[:, -1], 'o', label="Final position", color=cmap_cool(0.99)) #final position in purple
ax.axis("scaled")
ax.set_xlim(left=-24, right=24)

for i in range(20):
for bird in range(N):
ax.plot([X_pos[bird, i*(int(T/dt)//20)], X_pos[bird, (i+1)*(int(T/dt)//20)]], [Y_pos[bird, i*(int(T/dt)//20)], Y_pos[bird, (i+1)*(int(T/dt)//20)]], '--', color=cmap_cool(i/20), lw=1)
plt.legend(loc='best')



Hope this helps.

• The velocity dependent force needs something that repels $\vec v = 0$, or even small speeds, as birds don't fly well at $v=0$. With the exception of non-flocking humming birds, of course. – JEB Dec 9 '20 at 16:31
• Thank you for your comment. I forgot to mention it, but I am assuming the whole flock is already moving at a constant velocity $V_0$ and I look at the reference frame moving with the bird flock. This is just to decouple the average motion of the flock to the relative motion of birds. Whether or not this represents any physical reality is up to the reader. – QuantumApple Dec 9 '20 at 20:46