# Semi-linear Poisson-mediated Flocking in a Cucker-Smale Model

Christos N. Mavridis  Amoolya Tirumalai  John S. Baras  Ion Matei  Electrical
and Computer Engineering Department and the Institute for Systems Research,
University of Maryland, College Park, MD 20742 USA, (emails: $\\{$mavridis,
ast256<EMAIL_ADDRESS> Palo Alto Research Center (PARC), Palo Alto, CA
94304 USA
(email<EMAIL_ADDRESS>
###### Abstract

We propose a family of compactly supported parametric interaction functions in
the general Cucker-Smale flocking dynamics such that the mean-field
macroscopic system of mass and momentum balance equations with non-local
damping terms can be converted from a system of partial integro-differential
equations to an augmented system of partial differential equations in a
compact set. We treat the interaction functions as Green’s functions for an
operator corresponding to a semi-linear Poisson equation and compute the
density and momentum in a translating reference frame, i.e. one that is taken
in reference to the flock’s centroid. This allows us to consider the dynamics
in a fixed, flock-centered compact set without loss of generality. We approach
the computation of the non-local damping using the standard finite difference
treatment of the chosen differential operator, resulting in a tridiagonal
system which can be solved quickly.

###### keywords:

Control of Distributed Parameter Systems, Networked Control Systems, Large
Scale Systems

theoremstyle IFACtheoremstyle

††thanks: This material is based upon work supported by the Defense Advanced
Research Projects Agency (DARPA) under Agreement No. HR00111990027. This paper
has been accepted for publication in the 24th International Symposium on
Mathematical Theory of Networks and Systems (MTNS 2020), Cambridge, UK.

##  1 INTRODUCTION

Collective motion of autonomous agents is a widespread phenomenon appearing in
numerous applications ranging from animal herding to complex networks and
social dynamics (Okubo, 1986; Cucker and Smale, 2007; Giardina, 2008).

In general, there are two broad approaches when investigating the underlying
dynamics for flocks or swarms: the microscopic, particle models described by
ordinary differential equations (ODEs) or stochastic differential equations,
and the macroscopic continuum models, described by partial differential
equations (PDEs). Agent-based models assume behavioral rules at the individual
level, such as velocity alignment, attraction, and repulsion (Cucker and
Smale, 2007; Giardina, 2008; Ballerini et al., 2008) and are often used in
numerical simulations and in learning schemes where the interaction rules are
inferred (Matei et al., 2019). As the number of interacting agents gets large,
the agent-based models become computationally expensive (Carrillo et al.,
2010). Considering pairwise interactions, the growth is $O(N^{2})$, where $N$
is the number of agents. As we approach the mean-field limit, it is useful to
consider the probability density of the agents. Using Vlasov-like arguments
(Carrillo et al., 2010), we can construct an equation analogous to the Fokker-
Planck-Kolmogorov equation. We can then define momentum and density and
construct a system of compressible hydrodynamic PDEs (Carrillo et al., 2010;
Shvydkoy and Tadmor, 2017).

In flocking dynamics (Cucker and Smale, 2007; Carrillo et al., 2010), the
velocity alignment term is not only nonlocal but can also be nonlinear
(Shvydkoy and Tadmor, 2017; Mao et al., 2018). The computation of the
corresponding hydrodynamic equations with nonlocal forces becomes quite costly
due to the approximation of the convolution integrals or integral transforms
using the various quadrature methods. The simplest ‘quadrature’ method is the
Riemann sum, whose complexity is $O(n^{2})$, where $n$ is the number of grid
points, when estimating a convolution integral as a convolution sum in one
dimension. On the other hand, an equivalent solution may be obtained using
finite differences if the interaction kernel is associated with a differential
operator. If that operator can be put into a sparse form, ideally a
tridiagonal form, a solution can be obtained efficiently.

In this work, we modify the classical Cucker-Smale model of nonlocal particle
interaction for velocity consensus (Cucker and Smale, 2007; Ha et al., 2009).
We propose a family of parametric interaction functions in $\mathbb{R}^{d}$,
$d\in\\{1,2,3\\}$, that are Green’s functions for appropriately defined linear
partial differential operators, which allow us to speed-up computation of the
nonlocal interaction terms. We investigate the conditions under which time-
asymptotic flocking is achieved in the microscopic formulation in a centroid-
fixed frame. We solve the macroscopic formulation using the Kurganov-Tadmor
MUSCL finite volume method (Kurganov and Tadmor, 2000) and a second-order
finite difference discretization of our chosen differential operator. The
method is compared to bulk variables computed from the microscopic formulation
for validation.

The rest of the manuscript is organized as follows: Section 2 introduces the
agent-based Cucker-Smale flocking dynamics and the macroscopic Euler
equations. Section 3 describes the conversion of the Euler equations to an
augmented system of PDEs, and the formulation of the boundary value problem.
In Section 4 a family of interaction functions is proposed and the computation
process is explained. Finally, Section 5 compares the numerical results and
Section 6 concludes the paper.

##  2 Mathematical Models

In this section we introduce the Cucker-Smale dynamics under general
interaction functions, define time-asymptotic flocking, and present the mean-
field macroscopic equations.

###  2.1 The Cucker-Smale Model

Consider an interacting system of $N$ identical autonomous agents with unit
mass in $\mathbb{R}^{d}$, $d\in\\{1,2,3\\}$. Let $x_{i}(t),\
v_{i}(t)\in\mathbb{R}^{d}$ represent the position and velocity of the
$i^{th}$-particle at each time $t\geq 0$, respectively, for $1\leq i\leq N$.
Then the general Cucker-Smale dynamical system (Cucker and Smale, 2007) of
$(2N)$ ODEs reads as:

$\begin{cases}\frac{dx_{i}}{dt}&=v_{i}\\\
\frac{dv_{i}}{dt}&=\frac{1}{N}\sum_{j=1}^{N}\psi(x_{j},x_{i})(v_{j}-v_{i})\end{cases}$
(1)

where $x_{i}(0)$, are $v_{i}(0)$ are given for all $i=1,\ldots,N$, and
${\psi:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}}$ represents
the interaction function between each pair of particles.

The center of mass system $(x_{c},v_{c})$ of
$\left\\{(x_{i},v_{i})\right\\}_{i=1}^{N}$ is defined as

$x_{c}=\frac{1}{N}\sum_{i=1}^{N}x_{i},\quad
v_{c}=\frac{1}{N}\sum_{i=1}^{N}v_{i}$ (2)

When $\psi$ is symmetric, i.e., $\psi(x,s)=\psi(s,x)$, system (1) implies

$\frac{dx_{c}}{dt}=v_{c},\quad\frac{dv_{c}}{dt}=0$ (3)

which gives the explicit solution

$x_{c}(t)=x_{c}(0)+tv_{c}(0),\ t\geq 0$ (4)

###  2.2 Asymptotic Flocking

We investigate the additional assumptions on the initial conditions and the
interaction function $\psi$, such that system (1) converges to a velocity
consensus, a phenomenon known in the literature as _time-asymptotic flocking_
, defined in terms of the center of mass system as

######  Definition 1 (Asymptotic Flocking)

An $N-$body interacting system
$\mathcal{G}=\left\\{(x_{i},v_{i})\right\\}_{i=1}^{N}$ exhibits time-
asymptotic flocking if and only if the following two relations hold:

  * •

(Velocity alignment):
$\lim_{t\rightarrow\infty}\sum_{i=1}^{N}\left\|v_{i}(t)-v_{c}(t)\right\|^{2}=0$
,

  * •

(Spatial coherence): $\sup_{0\leq
t\leq\infty}\sum_{i=1}^{N}\left\|x_{i}(t)-x_{c}(t)\right\|^{2}<\infty$ .

We consider the new variables

$(\hat{x}_{i},\hat{v}_{i}):=(x_{i}-x_{c},v_{i}-v_{c})$ (5)

which correspond to the fluctuations around the center of mass system, and
define $\hat{x}:=(\hat{x}_{1},\ldots,\hat{x}_{N})$,
$\hat{v}:=(\hat{v}_{1},\ldots,\hat{v}_{N})$,
$|\hat{x}|=\left(\sum_{i=1}^{N}\|\hat{x}_{i}\|^{2}\right)^{1/2}$, and
$|\hat{v}|=\left(\sum_{i=1}^{N}\|\hat{v}_{i}\|^{2}\right)^{1/2}$, where
$\|\cdot\|$ represents the standard $l_{2}$-norm in $\mathbb{R}^{d}$. Based on
Definition 1, asymptotic flocking is achieved if

$|\hat{x}(t)|<\infty,t\geq 0,\ \text{and
}\lim_{t\rightarrow\infty}|\hat{v}(t)|=0$ (6)

We first notice that

$\frac{d|\hat{x}|^{2}}{dt}=2\left<\frac{d\hat{x}}{dt},\hat{x}\right>\leq
2|\hat{x}||\hat{v}|$ (7)

which implies

$\frac{d|\hat{x}|}{dt}\leq|\hat{v}|$ (8)

Suppose the interaction function $\psi$ is chosen such that
$\psi(x,s)=\tilde{\psi}(\|x-s\|)$, with
$\tilde{\psi}:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ being a non-negative
and non-increasing function. Then $(\hat{x}_{i},\hat{v}_{i})$ are governed by
the dynamical system (1), and

$\displaystyle\frac{d|\hat{v}|^{2}}{dt}$
$\displaystyle=-\frac{1}{N}\sum_{1\leq i,j\leq
N}\tilde{\psi}(\|\hat{x}_{j}-\hat{x}_{i}\|)\|\hat{v}_{j}-\hat{v}_{i}\|^{2}$
(9) $\displaystyle\leq-\frac{1}{N}\tilde{\psi}(2|\hat{x}|)\sum_{1\leq i,j\leq
N}\|\hat{v}_{j}-\hat{v}_{i}\|^{2}$
$\displaystyle=-\frac{2}{N}\tilde{\psi}(2|\hat{x}|)|\hat{v}|^{2}$

which implies

$\frac{d|\hat{v}|}{dt}\leq-\frac{2}{N}\tilde{\psi}(2|\hat{x}|)|\hat{v}|:=-\phi(|\hat{x}|)|\hat{v}|$
(10)

where we have used the fact that $\sum_{i=1}^{N}\hat{v}_{i}(t)=0$, $t\geq 0$,
and

$\max_{1\leq i,j\leq N}\|\hat{x}_{i}-\hat{x}_{j}\|\leq 2|\hat{x}|$ (11)

The following Theorem by (Ha et al., 2009) provides sufficient conditions for
time-asymptotic flocking.

###### Theorem 1

Suppose $(|x|,|v|)$ satisfy the system of dissipative differential
inequalities (8), (10) with $\phi\geq 0$. Then if
$|v(0)|<\int_{|x(0)|}^{\infty}\phi(s)ds$, there is a $x_{M}\geq 0$ such that
$|v(0)|=\int_{|x(0)|}^{x_{M}}\phi(s)ds$, and for every $t\geq 0$, $|x(t)|\leq
x_{M}$, and $|v(t)|\leq|v(0)|e^{-\phi(x_{M})t}$.

The following is an immediate consequence of Theorem 1.

###### Proposition 1

Let $\mathcal{G}=\left\\{(x_{i},v_{i})\right\\}_{i=1}^{N}$ be an $N-$body
interacting system with dynamics given by (1). Suppose
$\psi(x,s)=\tilde{\psi}(\|x-s\|)$, with
$\tilde{\psi}:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ being a non-negative
and non-increasing function. Then if
$|v(0)-v_{c}(0)|<\int_{|x(0)-x_{c}(0)|}^{\infty}\frac{2}{N}\tilde{\psi}(2s)ds$,
$\mathcal{G}$ exhibits time-asymptotic flocking.

###  2.3 The Mean-Field Limit

Consider the empirical joint probability distribution of the particle
positions and velocities $\\{x_{i},v_{i}\\}_{i=1}^{N}$

$F_{xv}^{N}(t,x,v):=\frac{1}{N}\sum_{i=1}^{N}\delta(x_{i},v_{i})$ (12)

where $\delta(\cdot,\cdot)$ is the Dirac measure on $\mathbb{R}^{2d}$. As the
number of particles $N\rightarrow\infty$, we can use McKean-Vlasov arguments
to show that the empirical distribution converges weakly to a distribution
whose density $f_{xv}$ evolves according to the forward Kolmogorov equation
(Carrillo et al., 2010)

$\begin{split}&\partial_{t}f_{xv}+\nabla_{x}\cdot(vf_{xv})+\nabla_{v}\cdot(Af_{xv})=0\\\
&A:=\int_{\mathbb{R}^{2d}}\psi(x,s)(w-v)f_{xv}(t,s,w)dsdw.\end{split}$ (13)

We define

$\begin{split}\rho(t,x)&:=\int_{\mathbb{R}^{d}}f_{xv}(t,x,v)dv\\\
m(t,x)&:=\rho(t,x)u(t,x):=\int_{\mathbb{R}^{d}}vf_{xv}(t,x,v)dv.\end{split}$
(14)

which are the marginal probability and momentum density functions.
Substituting these into (13) yields the following $(d+1)$ compressible Euler
equations with non-local forcing:

$\begin{cases}\partial_{t}{\rho}+\nabla_{x}\cdot(\rho u)=0\\\
\partial_{t}{(\rho u)}+\nabla_{x}\cdot(\rho u\otimes
u)=\rho\mathcal{L}_{\psi}(\rho u)-\rho u\mathcal{L}_{\psi}\rho\end{cases}$
(15)

where $u$ is the mean velocity, $\rho(0,x)$ and $u(0,x)$ are given and

$\mathcal{L}_{\psi}f(t,x):=\int_{\mathbb{R}^{d}}\psi(x,s)f(t,s)ds.$ (16)

##  3 Semi-linear Poisson Mediated Flocking

###  3.1 Conversion to a system of PDEs

We think of the function $\psi$ as a Green’s function, i.e., as the impulse
response of a linear differential equation, represented by the operator
$\mathcal{L}_{x}$, such that

$\mathcal{L}_{x}y(t,x)=g(t,x)$ (17)

implies

$y(t,x)=\int_{\mathbb{R}^{d}}\psi(x,s)g(t,s)ds$ (18)

which results in

$\mathcal{L}^{-1}_{\psi}=\mathcal{L}_{x}$ (19)

for all $t\geq 0$, where

$\displaystyle\mathcal{L}_{x}\psi(x,s)=\delta(x-s),\ x,s\in\mathbb{R}^{d}.$
(20)

Then the following proposition holds:

###### Proposition 2

Suppose $\psi$ is a Green’s function with respect to a linear differential
operator $\mathcal{L}_{x}$. Then system (15) is equivalent to the augmented
system of ($2d+2$) partial differential equations:

$\begin{cases}\partial_{t}{\rho}+\nabla_{x}\cdot(\rho u)=0\\\
\mathcal{L}_{x}y=\begin{bmatrix}\rho u&\rho\end{bmatrix}^{T}\\\
\partial_{t}{(\rho u)}+\nabla_{x}\cdot(\rho u\otimes u)=\sum_{i=1}^{d}(\rho
y_{i}-\rho u_{i}y_{d+1})\cdot\hat{e}_{i}\end{cases}$ (21)

where $\left\\{\hat{e}_{i}\right\\}_{i=1}^{d}$ is the standard basis in
$\mathbb{R}^{d}$.

###  3.2 The Boundary Value Problem

Due to the time-dependence of the center of mass (4), $x_{i}$, $i=1,\ldots,N$,
will escape any fixed and open bounded domain $\Omega\subset\mathbb{R}^{d}$,
unless in the trivial case where $v_{c}(0)=0$. Because of the flocking
behavior (Definition 1), the position fluctuations with respect to the center
of mass are uniformly bounded, i.e.,

$\sup_{0\leq
t\leq\infty}\sum_{i=1}^{N}\left\|x_{i}(t)-x_{c}(t)\right\|^{2}<\infty$ (22)

and, therefore we can define a Boundary Value Problem (BVP) in the moving
domain

$\Omega_{c}(t)=\left\\{x+x_{c}(t):x\in\Omega\right\\}$ (23)

where it is assumed that $0_{d}\in\Omega$, $0_{d}$ being the origin of
$\mathbb{R}^{d}$.

We notice that solving system (21) for $(x,u)$, $x\in\Omega_{c}$ is equivalent
to solving it for the fluctuation variables $(\hat{x},\hat{u})$ (5), with
$\hat{x}\in\Omega$.

We note that the boundedness of the domain has an effect on both the Green’s
function and the flocking behavior of the system of interacting particles,
which should satisfy

$x_{i}(t)-x_{c}(t)\in\Omega,\ i=1,\ldots,N,\ t\geq 0.$ (24)

##  4 One-Dimensional Case

The BVP of the augmented system of PDEs (21) for $d=1$, on
$\Omega=\\{\hat{x}\in[-\frac{L}{2},\frac{L}{2}]\\}$ reads as:

$\begin{cases}\partial_{t}{\rho}+\partial_{\hat{x}}(\rho u)=0\\\
\mathcal{L}_{\hat{x}}y=\begin{bmatrix}\rho u&\rho\end{bmatrix}^{T}\\\
\partial_{t}(\rho u)+\partial_{\hat{x}}(\rho u^{2})=\rho y_{1}-\rho
uy_{2}\end{cases}$ (25)

with homogeneous Dirichlet boundary conditions and initial conditions

$\rho(0,\hat{x})=\rho_{0}(\hat{x}),\quad u(0,\hat{x})=u_{0}(\hat{x})$ (26)

which are smooth functions.

We select the linear partial differential operator

$\mathcal{L}_{x}=-\frac{1}{2k}(\frac{\partial^{2}}{\partial
x^{2}}-\lambda^{2})$ (27)

with $k\neq 0$ and $\lambda\neq 0$, for which the associated parametric family
of Green’s functions with homogeneous Dirichlet boundary conditions on $[0,L]$
reads as:

$\hat{\psi}(x,s)=\begin{cases}c_{1}(s)(e^{\lambda x}-e^{-\lambda x})&s\leq
x\\\ c_{2}(s)(e^{\lambda(x-2L)}-e^{-\lambda x})&s>x\end{cases}$ (28)
$\displaystyle c_{1}(s)$
$\displaystyle=\frac{k}{\lambda(e^{-2L\lambda}-1)}(e^{\lambda(s-2L)}-e^{-\lambda
s})$ (29) $\displaystyle c_{2}(s)$
$\displaystyle=\frac{k}{\lambda(e^{-2L\lambda}-1)}(e^{\lambda s}-e^{-\lambda
s})$

The solution over any interval of length $L$ can be obtained by a simple
translation of coordinates.

The profile of the Green’s function $\hat{\psi}$ and the effect of the bounded
domain on on it is illustrated in Fig. 1, where, for different fixed values of
$x$, $\hat{\psi}(x,s)$ is compared to the function

$\psi(x,s)=\frac{k}{\lambda}e^{-\lambda\|x-s\|}$ (30)

which is the Green’s function corresponding to $\mathcal{L}_{x}$ in an
infinite domain. We note that the parameters $(k,\lambda,L)$ generate a family
of interaction functions (see also (Mavridis et al., 2020)) that can simulate
widely used interaction functions as the one found in the original Cucker-
Smale model (Cucker and Smale, 2007):

$G(x,s)=\frac{K}{(1+\|x-s\|^{2})^{\gamma}}$ (31)

for given parameters $(K,\gamma)$.

Figure 1: Illustration of $\hat{\psi}(x,\cdot)$ (28) for different values of
$x$, and for $\lambda=1$, $k=4$ on $[-\pi,\pi]$. The function
$\psi(x,s)=\frac{k}{\lambda}e^{-\lambda\|x-s\|}$, which is the Green’s
function for $\mathcal{L}_{x}$ in infinite domain, is depicted in the dashed-
dotted lines.

###  4.1 Asymptotic Flocking

Next we provide sufficient conditions such that the solution
$\left\\{(x_{i}(t),v_{i}(t))\right\\}_{i=1}^{N}$, $t\geq 0$, of system (1)
with interaction function $\hat{\psi}$ as defined in (28), (29), satisfy the
flocking conditions in Definition 1, with $\hat{x}_{i}(t)\in\Omega$, for all
$t\geq 0$.

Similar to Section 2.2, we notice that

$\frac{d|\hat{x}|}{dt}\leq|\hat{v}|$ (32)

From (11) and the fact that $\|\hat{x}_{i}\|\leq\max_{1\leq i,j\leq
N}\|\hat{x}_{i}-\hat{x}_{j}\|$, we get

$|\hat{x}|\leq\frac{\hat{x}_{M}}{2}\implies\|\hat{x}_{i}\|\leq\hat{x}_{M},\
i=1,\ldots,N$ (33)

Therefore, we are interested in showing asymptotic flocking with
$|\hat{x}(t)|\in[0,\frac{\hat{x}_{M}}{2}]$, for all $t\geq 0$.

For any given initial conditions $|\hat{x}_{i}(0)|$, there is a large enough
value of $L$ such that there exist an $\hat{x}_{M}\in[0,\frac{L}{2})$ for
which

$\hat{x}_{M}>2|\hat{x}(0)|$ (34)

From (28), (29) it follows that for $|\hat{x}|\leq\frac{\hat{x}_{M}}{2}$,

$\displaystyle\hat{\psi}(x_{j},x_{i})$
$\displaystyle\geq\hat{\psi}(-\hat{x}_{M},\|\hat{x}_{j}-\hat{x}_{i}\|)$ (35)
$\displaystyle\geq\hat{\psi}(-\hat{x}_{M},2|\hat{x}|)$

which implies that

$\displaystyle\frac{d|\hat{v}|^{2}}{dt}$
$\displaystyle=-\frac{1}{N}\sum_{1\leq i,j\leq
N}\hat{\psi}(\hat{x}_{j},\hat{x}_{i})\|\hat{v}_{j}-\hat{v}_{i}\|^{2}$ (36)
$\displaystyle\leq-\frac{2}{N}\hat{\psi}(-\hat{x}_{M},2|\hat{x}|)|\hat{v}|^{2}$

and

$\frac{d|\hat{v}|}{dt}\leq-\frac{2}{N}\hat{\psi}(-\hat{x}_{M},2|\hat{x}|)|\hat{v}|:=-\phi(|\hat{x}|)|\hat{v}|$
(37)

Next we notice that the Lyapunov function

$V(|x|,|v|):=|\hat{v}|+\int_{\alpha}^{|\hat{x}|}\phi(s)ds,\ \alpha\geq 0$ (38)

is non-increasing along the solutions of $(|\hat{x}(t)|,|\hat{v}(t)|)$ of the
system of dissipative differential inequalities (8) and (10), for
$|\hat{x}(t)|\leq\frac{\hat{x}_{M}}{2}$, since

$\displaystyle\frac{d}{dt}V(|\hat{x}|,|\hat{v}|)$
$\displaystyle=\frac{d|\hat{v}|}{dt}+\phi(|\hat{x}|)\frac{d|\hat{x}|}{dt}$
(39) $\displaystyle\leq\phi(|\hat{x}|)\left(-|v|+\frac{d|\hat{x}|}{dt}\right)$
$\displaystyle\leq 0$

which implies that

$|\hat{v}(t)|+\int_{|\hat{x}_{0}|}^{|\hat{x}|}\phi(s)ds\leq|\hat{v}(0)|,\
|\hat{x}|\leq\frac{\hat{x}_{M}}{2}$ (40)

Choosing the initial velocity $|\hat{v}(0)|$ such that
$|\hat{v}(0)|<\int_{|\hat{x}(0)|}^{\hat{x}_{M}/2}\phi(s)ds$, and, since $\phi$
is non-negative for $|\hat{x}(t)|\leq\frac{\hat{x}_{M}}{2}$, there exists a
$\bar{x}\in[|\hat{x}(0)|,\frac{\hat{x}_{M}}{2}]$ for which

$|\hat{v}(0)|=\int_{|\hat{x}(0)|}^{\bar{x}}\phi(s)ds$ (41)

Suppose there exists a $t^{*}\geq 0$, such that
$\hat{x}^{*}:=|\hat{x}(t^{*})|\in(\bar{x},\frac{\hat{x}_{M}}{2}]$. Then

$\int_{|\hat{x}(0)|}^{\hat{x}^{*}}\phi(s)ds>|v(0)|$ (42)

which contradicts (40). Therefore

$|\hat{x}(t)|\leq\bar{x}\leq\frac{\hat{x}_{M}}{2},\ t\geq 0$ (43)

and from (10) and the Grönwall-Bellman inequality

$|\hat{v}(t)|\leq|\hat{v}(0)|e^{-\phi(\bar{x})t},\ t\geq 0.$ (44)

###  4.2 Conservation of Mass and Momentum

###### Lemma 1

The operator (27) $\mathcal{L}_{x}$ on $C^{\infty}_{\mathbb{R},C}(\Omega)$,
the space of compactly supported test functions, is self-adjoint and
invertible, and therefore has a self-adjoint inverse $\mathcal{L}_{x}^{-1}$ on
$C^{\infty}_{\mathbb{R},C}(\Omega)$.

###### Proof.

Self-adjointness of the inverse follows immediately from self-adjointness of
$\mathcal{L}_{x}$ and the existence of the inverse (Taylor, 2010). It is clear
that $\mathcal{L}_{x}$ has an inverse since the Green’s function is
nontrivial.

We shall now show that the operator $\mathcal{L}_{x}$ is self-adjoint on
$C^{\infty}_{\mathbb{R},C}(\Omega)$. Consider two functions $u,w\in
C^{\infty}_{\mathbb{R},C}(\Omega)$, $u\neq w$, the space of test functions,
and associated $f_{u},f_{w}\in C^{\infty}_{\mathbb{R},C}$,
$f_{u}:=\mathcal{L}_{x}u,f_{w}:=\mathcal{L}_{x}w$. Let
$\Omega:=[-\frac{L}{2},\frac{L}{2}]$. We have

$\int_{\Omega}(w\mathcal{L}_{x}u-u\mathcal{L}_{x}w)dx=-\frac{1}{2k}\int_{\Omega}(w\partial_{x}^{2}u-u\partial_{x}^{2}w)dx.$
(45)

since the semi-linear term drops out. Using Green’s second identity, and the
compact support of $u,w$, we have that

$\int_{\Omega}(w\partial_{x}^{2}u-u\partial_{x}^{2}w)dx=\int_{\partial\Omega}(w\partial_{\mathbf{n}}u-u\partial_{\mathbf{n}}w)dx=0.$
(46)

Thus, $\mathcal{L}_{x}$ is self-adjoint and has a self-adjoint inverse, i.e.

$\int_{\Omega}(f_{w}\mathcal{L}_{x}^{-1}f_{u}-f_{u}\mathcal{L}_{x}^{-1}f_{w})dx=\int_{\Omega}(f_{w}u-f_{u}w)dx=0.$
(47)

∎

###### Proposition 3

If $y$ is compactly supported, and $psi$ is as given, then mass and momentum
are conserved, i.e.

$\frac{d}{dt}\int_{\Omega}\begin{bmatrix}\rho&\rho
u\end{bmatrix}^{T}d\hat{x}=\int_{\Omega}\begin{bmatrix}0&\rho y_{1}-\rho
uy_{2}\end{bmatrix}^{T}d\hat{x}=0.$ (48)

###### Proof.

We obtain (48) by simply integrating the conservation laws in (25) over the
entire space and apply the Leibniz rule. The conclusion follows directly from
the self adjointness of the inverse in (47). The proposition holds for any
self-adjoint alignment operator. ∎

###  4.3 Computational Methods

For compactness, we re-write the PDEs (25) as

$\begin{cases}\partial_{t}U+\partial_{\hat{x}}F(U)=S(U,Y)\\\
\mathcal{L}_{\hat{x}}Y=U\\\ \end{cases}$ (49)

with $U=\begin{bmatrix}\rho,\rho u\end{bmatrix}^{T}$,
$Y=\begin{bmatrix}y_{2},y_{1}\end{bmatrix}^{T}$, $F=\begin{bmatrix}\rho u,\rho
u^{2}\end{bmatrix}^{T}$, and $S=\begin{bmatrix}0,\rho y_{1}-\rho
uy_{2}\end{bmatrix}^{T}$. Recall the transformation $m=\rho u$. From this, the
flux Jacobian is given by

$\mathbf{D}_{U}F:=\begin{bmatrix}0&1\\\ -u^{2}&2u\end{bmatrix}$ (50)

which is not diagonalizable, and thus the system is only weakly hyperbolic.
Its eigenvalues are $\pm u$. With these notations established, we now detail
the numerical solution of the PDEs.

####  4.3.1 Hyperbolic Solver.

To solve the hyperbolic system, we apply the finite volume method (LeVeque,
2002). To begin, we define the sequence of points
$\\{\hat{x}_{0},...,\hat{x}_{i},...,\hat{x}_{N}\\}$ which are the centers of
the cells $I_{i}:=[\hat{x}_{i-\frac{1}{2}},\hat{x}_{i+\frac{1}{2}})$. Then, we
average the PDE over these cells, which gives

$\frac{1}{\lambda(I_{i})}\frac{d}{dt}\int_{I_{i}}Ud\hat{x}=-\frac{1}{\lambda(I_{i})}\int_{I_{i}}\partial_{\hat{x}}Fd\hat{x}+\frac{1}{\lambda(I_{i})}\int_{I_{i}}Sd\hat{x}$
(51)

where $\lambda(\cdot)$ denotes the length of an interval. Suppose these are
identical, so $\Delta\hat{x}:=\lambda(I_{i})\forall i$. Then, using the
divergence theorem, and replacing the integrals of $U,F,S$ with their cell-
averages, i.e. their midpoint values $\bar{U},\bar{F},\bar{S}$, we obtain

$\frac{d}{dt}\bar{U}_{i}=-\frac{1}{\Delta\hat{x}}(\bar{F}_{i+\frac{1}{2}}-\bar{F}_{i-\frac{1}{2}})+\bar{S}_{i}$
(52)

where
$\bar{U}_{i}:=\bar{Y}(\hat{x}_{i}),\bar{F}_{i}:=\bar{F}(\hat{x}_{i}),\bar{S}:=\bar{S}(\hat{x}_{i})$.
In this work, we employ the second-order strong stability preserving Runge-
Kutta scheme (Kurganov and Tadmor, 2000) for time integration. For the fluxes,
we assume piecewise linearity and use the Kurganov-Tadmor flux (Kurganov and
Tadmor, 2000). The fluxes are given by

$\begin{split}\bar{F}_{i+\frac{1}{2}}&:=\frac{1}{2}[F^{*}_{i}+F^{*}_{i+1}-\max\\{|u^{*}_{i}|,|u^{*}_{i+1}|\\}(U^{*}_{i+1}-U^{*}_{i})]\\\
U^{*}_{i+1}&:=U_{i+1}-\frac{\Delta\hat{x}}{2}minmod(\frac{U_{i+2}-U_{i+1}}{\Delta\hat{x}},\frac{U_{i+1}-U_{i}}{\Delta\hat{x}})\\\
U^{*}_{i}&:=U_{i}+\frac{\Delta\hat{x}}{2}minmod(\frac{U_{i+1}-U_{i}}{\Delta\hat{x}},\frac{U_{i}-U_{i-1}}{\Delta\hat{x}})\\\
\end{split}$ (53)

where $minmod(a,b):=\frac{1}{2}(sign(a)+sign(b))\min(|a|,|b|)$.

####  4.3.2 Elliptic Solver.

To solve the elliptic equations, we apply the classical second-order finite
difference method, which is

$\frac{y^{j}_{i+1}-2y^{j}_{i}+y^{j}_{i-1}}{\Delta\hat{x}^{2}}-\lambda^{2}y_{i}^{j}=-2kU^{j}_{i}$
(54)

Over the interior points, this yields linear equations

$(\frac{1}{\Delta\hat{x}^{2}}\mathbf{A}-\lambda^{2}\mathbf{I})y^{j}_{int}=-2kU^{j}_{int}-\frac{1}{\Delta\hat{x}^{2}}\begin{bmatrix}y^{j}_{0}&0&\ldots&0&y^{j}_{N}\end{bmatrix}^{T},$
(55) $\mathbf{A}=\begin{bmatrix}-2&1&0&\ldots&\ldots&0\\\ 1&-2&1&0&\ldots&0\\\
0&1&-2&1&\ldots&0\\\ \vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&0&0&0&1&-2\end{bmatrix}$ (56)

The matrix in (55) is tridiagonal, so banded matrix algorithms (Golub and Van
Loan, 2013) can be used to solve the corresponding system of equations. As
shown in Fig. 2, using finite differences is much faster than a convolution
(Riemann) sum, even when the embarrassing parallelism of the sum is exploited.

####  4.3.3 Particle Solver.

We solve the system of particle equations using the velocity Verlet algorithm
(Mao et al., 2018). Given a system of ODEs of the form

$\begin{cases}\frac{dx}{dt}&=v\\\ \frac{dv}{dt}&=a(x,v,t),\end{cases}$ (57)

with appropriate initial conditions and a time-discretization at steps
$\\{0,1,...,i,...\\}$ with increment $\Delta t$, the discretization is

$\begin{split}v_{i+\frac{1}{2}}&=v_{i}+\frac{1}{2}a(x_{i},v_{i},t_{i})\Delta
t\\\ x_{i+1}&=x_{i}+\Delta tv_{i+\frac{1}{2}}\\\ v_{i+1}&=v_{i}+\frac{\Delta
t}{2}[a(x_{i},v_{i},t_{i})+a(x_{i+1},v_{i+\frac{1}{2}},t_{i+1})].\end{split}$
(58)

Figure 2: Computation Times for Nonlocal Terms using Finite Differences and
Riemann Sum.

##  5 Numerical Results and Higher Dimensions

In this section we present numerical simulations of one-dimensional nonlocal
flocking dynamics, by solving $(a)$ the agent-based Cucker-Smale model using
the velocity Verlet method, and $(b)$ the macroscopic model with initial
conditions whose support is the interval $[-\pi,\pi]$. Our aim is to verify
that the agent based and continuum based approaches to the flocking problem
produce similar results.

In the following, the initial density and velocity are given by

$\displaystyle\rho_{0}(\hat{x})$
$\displaystyle=\frac{\pi}{2L}\cos(\frac{\pi\hat{x}}{L}),$ (59) $\displaystyle
u_{0}(\hat{x})$ $\displaystyle=-c\sin(\frac{\pi\hat{x})}{L}),\
\hat{x}\in[-\frac{L}{2},\frac{L}{2}],$ (60)

i.e. it is assumed that $\rho_{0}(\hat{x})=u_{0}(\hat{x})=0,\
\forall\hat{x}\notin[-\frac{L}{2},\frac{L}{2}]$, where we have used $L=2\pi$.

###  5.1 Cucker-Smale Model Simulation

In all simulations, we take $\lambda=1$, $k=4$. For the particle simulation,
we use $N=10^{4}$ particles. For the macro-scale simulation, we use
$\Delta\hat{x}=\frac{2\pi}{600}$ as the spatial increment. In both
simulations, we take $\Delta t=.001$ as the time increment.

In both cases, the support of the initial profile shrinks as the bulk comes
together. The semi-linear Poisson-forced Euler system is highly dissipative,
and the momentum profile is damped until it flattens (although it is conserved
over the domain), and the system attains an equilibrium distribution. Fig. 3
shows the agreement between the particle model and the macro-scale model.

Figure 3: Evolution of the Probability Densities $\rho(t,\hat{x})$ and
Momentum Densities $m(t,\hat{x})$ as computed by solving the macro-scale model
and the particle model (dashed-line).

###  5.2 Higher Dimensions

In higher dimensions, the radial symmetry of the interaction function $\psi$
suggests the use of a _singular kernel_. Singular kernels have been
extensively studied in the literature and, under mild assumptions in the
initial conditions, have been shown to result in flocking behavior while, at
the same time, avoiding collisions (Ahn et al., 2012).

In the BVP of the augmented system of PDEs (21) with the initial and boundary
conditions (26), we select the linear differential operator (see also
(Mavridis et al., 2020)):

$\mathcal{L}_{x}=-k^{-d/2}(\nabla_{x}^{2}-\lambda^{2})$ (61)

and $\Omega=B_{d}(0,r):=\left\\{x\in\mathbb{R}^{d}:\|x\|<r\right\\}$, which
results in a Green’s function of the form

$\displaystyle\hat{\psi}(x,s)=\psi(x-s)+\phi(x,s)$ (62)

where $\psi$ is given by

$\displaystyle\psi(x,s)$ $\displaystyle=\tilde{\psi}(\|x-s\|)$ (63)
$\displaystyle=\left(\frac{k}{2\pi}\right)^{d/2}\left(\frac{\lambda}{\|x-s\|}\right)^{d/2-1}K_{d/2-1}(\lambda\|x-s\|)$

with $K_{\alpha}(\cdot)$ being the modified Bessel function of the second kind
of order $\alpha$, and $\phi$ is a function such that

$\displaystyle\mathcal{L}_{s}\phi(x,s)=0,\ s\in B_{d}(0,r)$ (64)
$\displaystyle\phi(x,s)=-\psi(x,s),\ s\in\partial B_{d}(0,r)$

For $s\in\partial B_{d}(0,r)$ we have

$\displaystyle\|x-s\|^{2}$
$\displaystyle=\|x\|^{2}-2\left<x,s\right>+\|s\|^{2}$ (65)
$\displaystyle=\|x\|^{2}\|\frac{s}{r}-\frac{rx}{\|x\|^{2}}\|^{2}$

and it can be shown that

$\phi(x,s)=-\tilde{\psi}(\frac{1}{r}\|x\|\|s-r^{2}\frac{x}{\|x\|^{2}}\|).$
(66)

The interaction function $\hat{\psi}$ is affected by the bounded domain in the
same way as in the one-dimensional case, and depends on the parameter values
$k$ and $\lambda$ as illustrated in Fig.4 for the $2$-dimensional case.

Figure 4: The effect of the parameters $k$, $\lambda$ on the profile of the
interaction function $\hat{\psi}((0,0.5),s)$, $s\in B_{2}(0,1)$. Left:
$(k,\lambda)=(1,0.5)$. Right: $(k,\lambda)=(2,10)$.

##  6 Conclusion

A family of compactly supported parametric interaction functions in the
general Cucker-Smale flocking dynamics was proposed such that the macroscopic
system of mass and momentum balance equations with non-local damping terms can
be converted to an augmented system of coupled PDEs in a compact set. We
approached the computation of the non-local damping using the standard finite
difference treatment of the chosen differential operator, which was solved
using banded matrix algorithms. The expressiveness of the proposed interaction
functions may be utilized for parametric learning from trajectory data.

## References

  * Ahn et al. (2012) Ahn, S., Choi, H., Ha, S.Y., and Lee, H. (2012).  On collision-avoiding initial configurations to cucker-smale type flocking models.  _Communications in Mathematical Sciences_ , 10.  10.4310/CMS.2012.v10.n2.a10.
  * Ballerini et al. (2008) Ballerini, M., Cabibbo, N., Candelier, R., Cavagna, A., Cisbani, E., Giardina, I., Lecomte, V., Orlandi, A., Parisi, G., Procaccini, A., et al. (2008).  Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study.  _Proceedings of the national academy of sciences_ , 105(4), 1232–1237.
  * Carrillo et al. (2010) Carrillo, J.A., Fornasier, M., Toscani, G., and Vecil, F. (2010).  Particle, kinetic, and hydrodynamic models of swarming.  In _Mathematical modeling of collective behavior in socio-economic and life sciences_ , 297–336. Springer.
  * Cucker and Smale (2007) Cucker, F. and Smale, S. (2007).  Emergent behavior in flocks.  _IEEE Transactions on automatic control_ , 52(5), 852–862.
  * Giardina (2008) Giardina, I. (2008).  Collective behavior in animal groups: theoretical models and empirical studies.  _HFSP journal_.
  * Golub and Van Loan (2013) Golub, G.H. and Van Loan, C.F. (2013).  _Matrix Computations_.  The Johns Hopkins University Press, fourth edition.
  * Ha et al. (2009) Ha, S.Y., Liu, J.G., et al. (2009).  A simple proof of the cucker-smale flocking dynamics and mean-field limit.  _Communications in Mathematical Sciences_ , 7(2), 297–325.
  * Kurganov and Tadmor (2000) Kurganov, A. and Tadmor, E. (2000).  New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations.  _Journal of Computational Physics_ , 160(1), 241 – 282.
  * LeVeque (2002) LeVeque, R.J. (2002).  _Finite Volume Methods for Hyperbolic Problems_.  Cambridge Texts in Applied Mathematics. Cambridge University Press.  10.1017/CBO9780511791253.
  * Mao et al. (2018) Mao, Z., Li, Z., and Karniadakis, G. (2018).  Nonlocal flocking dynamics: Learning the fractional order of pdes from particle simulations.  _arXiv preprint arXiv:1810.11596_.
  * Matei et al. (2019) Matei, I., Mavridis, C., Baras, J.S., and Zhenirovskyy, M. (2019).  Inferring particle interaction physical models and their dynamical properties.  In _2019 IEEE Conference on Decision and Control (CDC)_ , 4615–4621. IEEE.
  * Mavridis et al. (2020) Mavridis, C.N., Tirumalai, A., and Baras, J.S. (2020).  Learning interaction dynamics from particle trajectories and density evolution.  In _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE.
  * Okubo (1986) Okubo, A. (1986).  Dynamical aspects of animal grouping: swarms, schools, flocks, and herds.  _Advances in biophysics_.
  * Shvydkoy and Tadmor (2017) Shvydkoy, R. and Tadmor, E. (2017).  Eulerian dynamics with a commutator forcing ii: Flocking.  _arXiv_.
  * Taylor (2010) Taylor, M. (2010).  _Partial Differential Equations I: Basic Theory_.  Applied Mathematical Sciences. Springer New York.