Opinion Dynamics on a General Compact Riemannian Manifold

  • Received: 01 March 2017 Revised: 01 July 2017
  • Primary: 34C40; Secondary: 37N40, 91B14

  • This work formulates the problem of defining a model for opinion dynamics on a general compact Riemannian manifold. Two approaches to modeling opinions on a manifold are explored. The first defines the distance between two points using the projection in the ambient Euclidean space. The second approach defines the distance as the length of the geodesic between two agents. Our analysis focuses on features such as equilibria, the long term behavior, and the energy of the system, as well as the interactions between agents that lead to these features. Simulations for specific manifolds, S1,S2, and T2 , accompany the analysis. Trajectories given by opinion dynamics may resemble n body Choreography and are called "social choreography". Conditions leading to various types of social choreography are investigated in R2 .

    Citation: Aylin Aydoğdu, Sean T. McQuade, Nastassia Pouradier Duteil. Opinion Dynamics on a General Compact Riemannian Manifold[J]. Networks and Heterogeneous Media, 2017, 12(3): 489-523. doi: 10.3934/nhm.2017021

    Related Papers:

    [1] Aylin Aydoğdu, Sean T. McQuade, Nastassia Pouradier Duteil . Opinion Dynamics on a General Compact Riemannian Manifold. Networks and Heterogeneous Media, 2017, 12(3): 489-523. doi: 10.3934/nhm.2017021
    [2] Yuntian Zhang, Xiaoliang Chen, Zexia Huang, Xianyong Li, Yajun Du . Managing consensus based on community classification in opinion dynamics. Networks and Heterogeneous Media, 2023, 18(2): 813-841. doi: 10.3934/nhm.2023035
    [3] Mayte Pérez-Llanos, Juan Pablo Pinasco, Nicolas Saintier . Opinion fitness and convergence to consensus in homogeneous and heterogeneous populations. Networks and Heterogeneous Media, 2021, 16(2): 257-281. doi: 10.3934/nhm.2021006
    [4] Sergei Yu. Pilyugin, Maria S. Tarasova, Aleksandr S. Tarasov, Grigorii V. Monakov . A model of voting dynamics under bounded confidence with nonstandard norming. Networks and Heterogeneous Media, 2022, 17(6): 917-931. doi: 10.3934/nhm.2022032
    [5] Clinton Innes, Razvan C. Fetecau, Ralf W. Wittenberg . Modelling heterogeneity and an open-mindedness social norm in opinion dynamics. Networks and Heterogeneous Media, 2017, 12(1): 59-92. doi: 10.3934/nhm.2017003
    [6] Marina Dolfin, Mirosław Lachowicz . Modeling opinion dynamics: How the network enhances consensus. Networks and Heterogeneous Media, 2015, 10(4): 877-896. doi: 10.3934/nhm.2015.10.877
    [7] Sabrina Bonandin, Mattia Zanella . Effects of heterogeneous opinion interactions in many-agent systems for epidemic dynamics. Networks and Heterogeneous Media, 2024, 19(1): 235-261. doi: 10.3934/nhm.2024011
    [8] Sergei Yu. Pilyugin, M. C. Campi . Opinion formation in voting processes under bounded confidence. Networks and Heterogeneous Media, 2019, 14(3): 617-632. doi: 10.3934/nhm.2019024
    [9] Sharayu Moharir, Ananya S. Omanwar, Neeraja Sahasrabudhe . Diffusion of binary opinions in a growing population with heterogeneous behaviour and external influence. Networks and Heterogeneous Media, 2023, 18(3): 1288-1312. doi: 10.3934/nhm.2023056
    [10] Rainer Hegselmann, Ulrich Krause . Opinion dynamics under the influence of radical groups, charismatic leaders, and other constant signals: A simple unifying model. Networks and Heterogeneous Media, 2015, 10(3): 477-509. doi: 10.3934/nhm.2015.10.477
  • This work formulates the problem of defining a model for opinion dynamics on a general compact Riemannian manifold. Two approaches to modeling opinions on a manifold are explored. The first defines the distance between two points using the projection in the ambient Euclidean space. The second approach defines the distance as the length of the geodesic between two agents. Our analysis focuses on features such as equilibria, the long term behavior, and the energy of the system, as well as the interactions between agents that lead to these features. Simulations for specific manifolds, S1,S2, and T2 , accompany the analysis. Trajectories given by opinion dynamics may resemble n body Choreography and are called "social choreography". Conditions leading to various types of social choreography are investigated in R2 .



    The emergence of a group's global behavior from local interactions among individual agents is a fascinating feature of opinion dynamics. When local rules imply global patterns in a population, we are observing a phenomenon called self-organization. Traditionally, interest focuses on understanding the complex rules of interacting opinions which lead to certain global configurations, such as classic consensus, alignment, clustering, or the less studied dancing equilibrium [3]. For instance, in bounded-confidence models such as the one proposed by Hegselmann and Krause, the radius of interaction determines the clustering of the system [11]. Motsch and Tadmor studied the influence of the shape of the interaction potential on the convergence to consensus of the Hegselmann-Krause system [15]. Ha, Ha and Kim looked at the Cucker-Smale second-order alignment model and provided a condition on the interaction potential ensuring convergence of the system to alignment [9]. Cristiani, Frasca and Piccoli studied the effect of anisotropic interactions on the behavior of the group [6].

    The dynamics of an opinion formation system depend on the state-space [1] and interaction network [11]. Models on the Euclidean space in one dimension (for opinion dynamics) or in two or three dimensions (with applications to groups of animals or robots) have been extensively studied and are well understood. However, such models are locally linear, which may be a limitation when one strives to capture more complex phenomena and better represent reality [21]. In this line of thought, the Kuramoto model on the sphere S1 addresses the problem of synchronizing a large number of oscillators [12,22]. There exist numerous applications to this model [7,18,19,20]. Similarly, applications to satellite or ground vehicle coordination have motivated the development of models on special orthogonal groups [16,17]: satellites evolve on SO(3) while ground vehicles evolve on SE(2) or SE(3). A nonlinear model of opinion formation on the sphere was also developed in [3]. A discussion of the so called ''flocking realizability problem'' for a sphere is given in [5], which focuses on the particular equilibrium that we refer to as consensus for a holonomic dynamical system on a sphere.

    The present work defines a general model of opinion dynamics on a Riemannian manifold. We investigate how the manifold on which the model is defined affects the global configurations resulting from opinion dynamics. These are the first steps to build a robust theory of opinion dynamics on general Riemannian manifolds.

    There is an inherent difficulty in defining opinion dynamics on a general Riemannian manifold. Using the Riemannian distance, an agent will move towards a point by following the manifold's geodesics, which are well defined only locally. On a larger scale, there might not exist a unique geodesic. Another challenge is the extreme complexity of computing geodesics, even on a relatively simple manifold such as the torus [8]. One way around this issue is to consider the embedding of the manifold into a Euclidean space. Each agent's velocity is defined by projection of the other agents' influence onto the tangent space at that point. This is the choice made in [3].

    Other than the mentioned practical aspect, there is an intrinsic rationale for choosing one approach over the other. When evolving along the geodesics of the manifold, one assumes that each agent has a global understanding of the manifold's geometry and is able to choose the shortest path among all possible ones. On the other hand, the approach based on the projection of the desired destination onto the tangent space implies that each agent only holds local information about the space in which it evolves. It chooses to move in the direction which locally seems to bring it closer to the target.

    We explore these two specific approaches for our generalized model. The first method, Approach A, uses projections in the Euclidean space in which the manifold is embedded. The second method, Approach B, uses only geodesics defined on the manifold to define strength and direction of interaction. We exhibit properties of the interaction matrix that lead to specific kinds of equilibria. Simulations and examples compare the two methods. Dancing equilibria for Approach B are shown (dancing equilibria were studied for Approach A in [3]).

    We use the sphere and torus as sample manifolds to evaluate these approaches. Specifically, we simulate dynamics on the following manifolds: S1,S2 and T2. These examples allow us to directly compare the two approaches, and see if one is more appropriate for a given manifold. We show the influence of the manifold's geometry on the dynamics by examining the dynamics resulting from the same interaction matrix in S2, T2 and R2.

    Opinion dynamics trajectories can resemble nbody choreography, that is, solutions to the well known nbody problem. These dynamics drive agents along orbits which either are periodic, or have a periodic feature, and that may be shared by multiple agents. We refer to opinion dynamics trajectories along such orbits as ''Social Choreography''. We show that a simple example of Social Choreography in R2 does not hold on S2 or T2, see Figures 10 and Figures 11. We then focus on the case of R2 and investigate initial conditions and properties of the interaction matrix which give rise to Social Choreography. Future investigations can consist of exploring similar properties of trajectories on general manifolds.

    This work will primarily discuss two approaches to define opinion dynamics on a Riemannian manifold. Let M be a Riemannian manifold. Let NN represent the number of agents with opinions evolving on M. We denote by x:=(xi)i{1,...,N}MN the set of opinions. For each i{1,...,N}, ˙xiTxiM. The opinions xi evolve according to the following general dynamics:

    ˙xi=Nj=1aijΨ(d(xi,xj))νij (1)

    where

    aijR is the interaction coefficient of the pair of agents i and j,

    Ψ:RR is the interaction potential,

    d(,):M×MR+ represents the distance between opinions,

    νijTxiM is a unit vector giving the direction of the influence of j over i.

    Each of these terms is further specified in the following.

    The evolution of each agent's opinion depends on the opinions of all other agents, with influences weighted by the interaction coefficients aij. More specifically, an agent xj's influence on xi is determined by two elements: the direction of influence νijTxiM and the magnitude of influence Ψ(d(xi,xj))R+. We propose and study two different approaches for the choices of d and νij. Approach A uses the embedding of M in Rn to define d(xi,xj), whereas Approach B is intrinsic to M, with distance and direction of influence based on geodesics.

    Approach A. Assume that M of dimension m is embedded in a Euclidean space Rn, with nm. Agent xj acts on agent xi via a projection onto TxiMRm. Now considering points (xi,xj)M2 as points of Rn, the difference xjxi is a vector of Rn. Given a vector subspace Y of Rn, we denote by ΠYy the projection of yRn onto YRn and define dP(,) as follows:

    dP(xi,xj)=ΠTxiM(xjxi) (2)

    where denotes the Euclidean norm on Rn. The same projection also defines the direction of influence of xj on xi:

    νij={ΠTxiM(xjxi)ΠTxiM(xjxi) if ΠTxiM(xjxi)00   otherwise. (3)

    With the specific choice ΨId, system (1)-(2)-(3) becomes:

    ˙xi=Nj=1aijΠTxiM(xjxi). (4)

    This is the approach used in [3], applied to the sphere S2.

    Notice that the magnitude of influence, dP(xi,xj), is symmetric for the sphere in the sense that dP(xi,xj)=dP(xj,xi), but not symmetric for a general Riemannian manifold (see Figure 1). However, it is a continuous function defined for all pairs of points (xi,xj)M2. The originality of this approach is that the influence of xj on xi is not related to a notion of distance between the points. The use of the projection of xjxi onto TxiM reflects the concept of ''local visibility.'' For the situation of two agents evolving on a one dimensional manifold, if xjxiTxiM, then a local displacement of xi does not affect the distance between the points xixj. Indeed, a first order Taylor expansion gives: xi(ε)=xi(0)+ε˙xi(0)+o(ε).

    Figure 1.  An example of a manifold M such that dp(xi,xj)dp(xj,xi), Using system (4), an agent is subject to ''local visibility'', and movement of xi along TxiM (dashed line through xi) will not bring xi closer to xj in this local sense.

    Supposing that xj is fixed, we have:

    xi(ε)xj2=xi(ε)xj,xi(ε)xj=xi(0)xj,xi(0)xj+2ε˙xi(0),xi(0)xj+o(ε) (5)

    so if xjxi(0)Txi(0)M, then xi(ε)xj2=xi(0)xj2+o(ε). Hence if xi only has local visibility, all directions of displacement seem equivalent (at first order), which justifies the influence of xj over xi to be zero if their difference is orthogonal to the tangent space of M at xi. This is illustrated in Figure 1.

    Approach B. This second approach defines d and νij using the manifold M itself, and does not require any reference to the space in which M is immersed. This would make Approach B a natural way to define system dynamics, however the complete knowledge of the geodesics between any two points on the manifold may be unrealistic. Furthermore, the geometry of the manifold may introduce difficulties to the uniqueness of νij, particularly at the cut-locus of a point.

    Definition 2.1. The cut locus of a point qM is the set of points CL(q)M for which there are multiple geodesics between q and pCL(q) (see also [4]).

    Let γij:[0,1]M denote a geodesic connecting xi to xj,γij(0)=xi and γij(1)=xj. We then define the distance between xj and xi as the length of a geodesic, i.e. denoting by gy:TyM×TyMR+ the Riemannian metric at point yM,

    dG(xi,xj)=10gγij(s)(˙γij(s),˙γij(s))ds. (6)

    The direction of influence is determined by the same geodesic:

    νij={0 if xj=xi or if xjCL(xi)˙γij(0)gxi(˙γij(0),˙γij(0)) otherwise. (7)

    Unlike in Approach A, the magnitude of influence is a symmetric function: dG(xi,xj)=dG(xj,xi). Furthermore, this approach ensures that the magnitude of influence of one agent on another is a function of the exact Riemannian distance between the agents.

    Interaction networks. In finite-dimensional systems such as system (1), the set of interacting agents can be described by vertices of a graph. A directed edge exists from a vertex i to a vertex j if and only if aij0. The system depends on the interaction network, and likewise, if the coefficients aij are chosen to be functions of the state, the interaction network may change as a result of the dynamics. Two main types of interaction networks have been proposed in the literature: metric interactions and topological interactions. If interactions between agents occur only locally, only the neighbors of agent i influence agent i. Metric interactions define the set of neighbors of agent i, given a radius r>0, as

    Sri(x)={j{1,,N},d(xi,xj)r}, (8)

    where d(,) can represent either the projection or the geodesic distance, as specified in each of the two approaches described above (see equations (2) and (6)). The other main type of interactions specifies that an agent is influenced by only its k closest neighbors. We call these topological interactions [1]. We define the relative separation between two agents as αij=card{k:d(xi,xk)d(xi,xj)}, The set of neighbors of agent i is then defined as the set of its k closest neighbors, i.e. for a given kN,

    Ski(x)={j{1,,N},αijk}. (9)

    Figures 2 and 3 illustrate differences between the metric and topological networks for the specific example of S1, with each of the approaches A and B.

    Figure 2.  The set of agents that influence x1 depends on how the interaction network is defined. In (a) and (b) the dashed lines show the projection of agents onto the tangent space of xi, (TxiS1). The agents depicted in red with larger dots influence x1. With the same configuration on S1, four combinations are possible (approach {A, B} type {Metric, Topological}). Each combination implies x1 interacts with a different set of agents.
    Figure 3.  The agent x1 is influenced by different agents depending on how the interaction network is defined. These networks may change as the dynamics move the agents on S1. Each agent xj,j{1,,6} will have a network describing which other agents influence xj. The interaction networks corresponding to systems from Figure 2.

    Resolution of discontinuities. The definitions of νij for approaches A and B (given by equations (3) and (7)) allow discontinuities of νij at certain points. Thus, one must impose conditions on the interaction potential ΨC0(R+,R+), in order to ensure the continuity of the right-hand side of the system (1), and hence the existence and uniqueness of a solution. Table 1 lists the discontinuities of νij and gives necessary conditions on Ψ to ensure the continuity of Ψ(d(xi,xj))νij.

    Table 1.  Possible discontinuities of the right-hand side of (1). The bottom row of the table show conditions for Ψ so that the system is continuous.
    ApproachABA and B
    Critical pointsxjNi xjCL(xi) xj=xi
    DiscontinuitieslimxjNiνij=1
    νij=0 for xjNi
    limxjCL(xi)νij=1
    νij=0 for xjCL(xi)
    limxjxiνij=1
    νii=0
    Condition on ΨΨ(0)=0Ψ(d)=0 for all dϵΨ(0)=0

     | Show Table
    DownLoad: CSV

    Firstly, notice that in both approaches, νij is discontinuous at the point xi=xj. Indeed, if xi=xj, νij=0, whereas almost everywhere else, νij=1. To ensure the continuity of Ψ(d(xi,xj))νij at this point, we impose the following condition:

    Ψ(0)=0. (10)

    In Approach A, we created a discontinuity of νij at the points xjN(xi), where we denote by N(q) the set N(q):={qM|ΠTpM(qp)=0}. For convenience of notation, we will use interchangeably the notations N(xi) and Ni. More specifically, we have limxjNiνij=1 but νij=0 if xjNi (see also Table 1). However, from the definition of dP (see equation (2)), we have limxjNidP(xi,xj)=0 and d(xi,xj)=0 for xjNi. Hence a sufficient condition for Ψ(d(xi,xj))νij to be continuous is again:

    Ψ(0)=0. (11)

    In Approach B, there is a discontinuity for xjCL(xi). Denoting by Bgeo(p,ρ) the geodesic ball of center p and radius ρ, we require the following condition on the influence function Ψ:

    Ψ(d)=0 for all dϵ (12)

    where ϵ:=inf{ρ>0|pM,Bgeo(p,ρ)CL(p)=}. This distance ϵ, also known as injectivity radius, is known to exist and be greater than 0 for any compact Riemannian manifold (see [4]).

    Notice that in the case of the geodesics approach (Approach B), the condition Ψ(d)=0 for all dϵ is incompatible with the use of the topological network (9). Indeed, if agent j is among the k closest neighbors of agent i, the topological network would require: aij0. However, the interaction between i and j would be canceled if dG(xi,xj)>ϵ. On the other hand, the metric interaction network as defined by (8) is compatible with Approach A, and with Approach B if the interaction radius is smaller than the injectivity radius: rϵ. For simplicity purposes, in the rest of this paper, we will consider that the interaction coefficients aij are constant, thus not requiring the need to differentiate between metric and topological networks. While models with constant interaction coefficients are our focus here, these models are quite restrictive, and exclude all models with dynamic interactions.

    Definition 2.2. The configuration x1=...=xN is called consensus. On the sphere, Sn, A configuration such that, for every j{2,,N},eitherxj=x1orxj=x1, which is not a concensus is called antipodal equilibrium.

    Proposition 1. The consensus configuration is an equilibrium for system (1).

    Proof. In both approaches A and B, if xi=xj, then νij=0. Hence if x1=...=xN, then for all i{1,,N}, ˙xi=0.

    Proposition 2. Let N>d+1. Then for every ˉx=(ˉx1,,ˉxN)MN, there exists a square matrix A=(aij)i,j{1,...,N} such that ˉx is an equilibrium for system (1).

    Proof. The configuration ˉx=(ˉx1,,ˉxN) is an equilibrium if and only if

    ddtˉxi=Nj=1aijΨ(d(ˉxi,ˉxj))νij=0.

    This is a system of at most Nd equations in the N2N unknowns aij,ij, notice that Ψ(d(xi,xi)νii=0, and diagonal values of A do not change the system. So if N>d+1 there exists a nontrivial choice of the interaction coefficients for which ˉx is an equilibrium.

    Definition 2.3. The kinetic energy of System (1)-(2)-(3) is the quantity

    EP(t):=12Ni=1˙xi(t)2. (13)

    The kinetic energy of System (1)-(6)-(7) is the quantity

    EG(t):=12Ni=1gxi(˙xi(t),˙xi(t)). (14)

    Proposition 3. Let M be a general Riemannian compact manifold. Consider the dynamics given by projection onto the tangent space (Approach A) given by (4). If the interaction matrix A=(aij)i,j{1,...,N}2 is symmetric, then

    limtEP(t)=0. (15)

    Proof. Let F(t)=12Ni=1Nj=1aijxixj2. Using the symmetry of A, we prove that

    ddtF(t)=4EP(t). (16)

    Indeed, notice that

    xi(Nj=1aijxixj2)=2ΠTxiMNj=1aij(xjxi)=2˙xi.

    Then we compute

    ddtF(t)=Nk=1xk12Ni,j=1aij(xixj2),˙xk=Nk=1xk[12Ni=1aik(xixk2)+12Nj=1akj(xkxj2)],˙xk=Nk=12xk12Ni=1aik(xixk2),˙xk=Nk=12ΠTxkMNj=1akj(xjxk),˙xk=2Nk=1˙xk2=4EP(t). (17)

    where the third equality uses the property: aij=aji for all i,j.

    Since EP(t)0, F(t) is a non-decreasing function. Moreover F(t) and d2dt2F(t) are bounded, since M is a compact manifold. Hence F(t) converges as t.

    Now, +0EP< and there exists c>0 such that suptddtEP(t)<c. By contradiction, assume lim supEP=α>0, then there exists a sequence (tn) such that tn+1>tn+α2c and EP(tn)α2. Then 0EP>nαcα212=+. Hence lim suptEP(t)=0.

    This shows that ddtF(t)0 when t, which implies that limtEP(t)=0.

    Remark 1. Propositions 2 and 3 are generalizations of results proven for the case M=S2 in [3].

    Remark 2. Proposition 3 assumes that ΨId which creates a discontinuity for Approach B (see Table 1). A result is shown for the more restricted case of M=S2 and Ψ()sin(), (Corollary 1).

    Definition 2.4. Let x solve the differential equation (1). A dancing equilibrium is a configuration in which for all pairs of agents (i,j), the distance dP(xi,xj) (in Approach A) or dG(xi,xj) (in Approach B) is constant. In the context of a system of oscillators, this equilibrium is also known as a phase-locked state or entrainment state [10].

    Remark 3. This definition is a generalization of the concept of dancing equilibrium described in [3].

    Remark 4. It follows immediately from definition 2.4 that the kinetic energy of a system in dancing equilibrium is constant.

    We study both approaches A and B in the case M=S1, i.e. for the one-dimensional sphere embedded in R2. Let (θi)i{1,...,N}[0,2π]N such that for all i{1,...,N}, xi=(cosθi,sinθi)T.

    Approach A. The projection onto an agent's tangent space can be rewritten as:

    ΠTxiNj=1aij(xjxi)=Nj=1aij(cosθjsinθj)(cosθisinθi),(sinθicosθi)(sinθicosθi)=Nj=1aij(sinθicosθj+sinθjcosθi)(sinθicosθi)=Nj=1aijsin(θjθi)(sinθicosθi). (18)

    So System (1)-(2)-(3) becomes:

    for alli{1,,N},˙θi(sinθicosθi)=Nj=1aijΨ(|sin(θjθi)|)sgn(sin(θjθi))(sinθicosθi) (19)

    where sgn() is the sign function defined by:

    for allxR,sgn(x)={1 if x>01 if x<00 if x=0. (20)

    We can then specify:

    for all(i,j){1,,N}2,dP(xi,xj)=|sin(θjθi)|,νPij=sgn(sin(θjθi)). (21)

    This gives the system of scalar equations:

    for alli{1,,N},˙θi=Nj=1aijΨ(|sin(θjθi)|)sgn(sin(θjθi)). (22)

    In particular, in the case ΨId, the system becomes the Kuramoto model [12].

    for alli{1,,N},˙θi=Nj=1aijsin(θjθi). (23)

    Approach B. For M=S1, the geodesics distance dG and the vector νGij are given by:

    dG(xi,xj)=arccos(cos(θjθi)),    νGij=sgn(sin(θjθi)). (24)

    System (1)-(6)-(7) is written:

    for alli{1,,N},    ˙θi=Nj=1aijΨ(arccos(cos(θjθi)))sgn(sin(θjθi)). (25)

    In order for the system to be well defined, the interaction function Ψ must satisfy the conditions given in Table 1. Notice that the injectivity radius is constant over S1, with ϵ=π. Possible choices involve choosing Ψ from a family of function defined as follows:

    Ψa(d)={1adfordadπaπford>a (26)

    where a(0,π) (see Figure 8).

    Figure 4.  Initial (empty circles) and final positions (filled circles) of 4 agents initially on the vertices of a rectangle with Approach B (left) and Approach A (right), with A=1, ΨId (Approach A) and Ψ=Ψ3π/4(Approach B) (see equation (26)). Notice that with Approach A, initial and final positions are identical since any rectangle configuration is an equilibrium. However, with Approach B, the system reaches a square configuration, the only possible equilibrium with pairwise distinct positions.
    Figure 5.  Evolution of the system (27) with Approach A (left) Approach B (center) when the interaction matrix satisfies condition (28) for the projection distance. Right: Kinetic energy.
    Figure 6.  Evolution of the system (27) with Approach A (left) Approach B (center) when the interaction matrix satisfies condition (28) for the geodesic distance. Right: Kinetic energy.
    Figure 7.  A dancing equilibrium for Approach A. The energy becomes constant in time after initial fluctuations.
    Figure 8.  The left side shows candidates for the choice of function Ψ. The right side shows how choice of function determines the energy of the system, for the case of a=3π4 the system forms an antipodal equilibrium.

    Another possible choice is: Ψ:xsin(x). Notice that for the specific choices Ψ=Id for Approach A and Ψ:xsin(x) for Approach B, the two approaches A and B are equivalent.

    We first examine the different equilibria for both approaches.

    Theorem 3.1. Consider Approach A, System (22). Let NN be even. Suppose that for all i{1,,N} for all j{1,N2}, aij=ai(j+N2). Then any configuration that is centrally symmetric, i.e.

    forallj{1,,N2},θj+N2=θj+π

    is an equilibrium.

    Proof. Using the hypotheses from Theorem 3.1, we can easily compute:

    ˙θi=Nj=1aijΨ(sin(θjθi))sgn(sin(θjθi))=N/2j=1[aijΨ(sin(θjθi))sgn(sin(θjθi))+ai(j+N2)Ψ(sin(θj+N2θi))sgn(sin(θj+N2θi))]=N/2j=1[aijΨ(sin(θjθi))sgn(sin(θjθi))+aijΨ(sin(θj+πθi))sgn(sin(θj+πθi))]=0.

    Interestingly, Theorem 3.1 is not applicable to Approach B. We illustrate the different behaviors of the two systems by studying the specific example of four agents initially in a rectangular configuration. According to Theorem 3.1, this configuration is an equilibrium for Approach A, independently of the choice of interaction function Ψ. However, one can easily prove that in the geodesics-based Approach B, with N=4 and the choice Ψ:=Ψa with a=3π4, the only equilibrium for which all agents have pairwise distinct positions is obtained by a regular polygon, i.e. all agents are evenly spaced out on the circle. This is illustrated by numerical simulations shown in Figure 4.

    This highlights the fundamentally different behaviors of the systems (1)-(2)-(3) and (1)-(6)-(7) in the case M=S1.

    In both approaches A and B, conditions on the interaction matrix A can be found such that the system forms a dancing equilibrium (see Definiton 2.4).

    Theorem 3.2. Consider the dynamics on S1 given by:

    foralli{1,,N},˙θi=Nj=1aijΨ(d(xi,xj))νij (27)

    where d(,) and ν are given either by Approach A (21) or Approach B (24). LetCR and suppose that for all i{1,,N},

    aij={CΨ(d(xi(0),xj(0)))νijifΨ(d(xi(0),xj(0)))00otherwise. (28)

    Then the system is in a dancing equilibrium.

    Proof. If the interaction matrix satisfies (28), then at t=0,

    for alli{1,,N},    ˙θi(0)=Nj=1C=CN

    so for all (i,j){1,,N}2, ˙θi(0)˙θj(0)=0. Then d(xi,xj) does not change in time, and (28) holds for all time.

    Numerical simulations show the evolution of the system (27) with condition (28) for the projection or the geodesic distance, see Figures 5 and 6.

    We study both approaches A and B for M=S2, i.e. for a two dimensional sphere embedded in R3. We use spherical coordinates: let (θi)i{1,...,N}[0,2π]N, and (ϕi)i{1,...,N}[0,π]N such that for all i{1,...,N},

    xi=(cosθsinϕ,sinθsinϕ,cosϕ)T.

    Choice of influence function. We choose an influence function Ψ(d) between two agents xi and xj so that the right-hand side of the system is continuous, the discontinuities are shown in Table 1. For Approach B, the only point in CL(xi) for a given xi is the antipodal point (this is an end point of a diameter for which xi is the other end point.) As in the case of S1, for Approach B, we choose a function Ψ from a family of functions of the form Ψa, see equation (26) (choices of Ψ are shown in Figure 8).

    Approach A. On S2, the derivative for system (1)-(2)-(3) with ΨId reduces to the sum of all projections onto the tangent space of agent xi, weighted by the corresponding interaction term aij. This is rewritten as:

    ΠTxiNj=1aij(xjxi)=ΠTxiNj=1aij(xj)=Nj=1aij(xjxj,xixi)=Nj=1aij((cosθjsinϕjsinθjsinϕjcosϕj)(cosθjsinϕjsinθjsinϕjcosϕj),(cosθisinϕisinθisinϕicosϕi)(cosθisinϕisinθisinϕicosϕi)).

    Approach B. The geodesic distance dG(xi,xj) from (6) between two points xi, and xj on S2 is given by:

    dG(xi,xj)=2arcsin(xixj2),

    and the direction toward xj from xi is

    νij=xjxj,xixixjxj,xixi,

    where is the standard norm in R3.

    Noticing that for S2, Approach A with ΨId is equivalent to Approach B with Ψsin, we extend the results of Proposition 3.

    Corollary 1. Consider the dynamics given by geodesic distance (Approach B) on S2, system (1)-(6)-(7), and let Ψsin. If the interaction matrix A=(aij)i,j{1,...,N}2 is symmetric, then

    limtEG(t)=0. (29)

    Proof. System (1)-(6)-(7), with Ψ()sin() reads:

    ˙xi=Nj=1aijsin(dG(xi,xj))νij. (30)

    Considering the embedding of the system in R3, we notice that for all i,j{1,,N},

    aijsin(dG(xi,xj))νij=aijΠTxiM(xjxi), (31)

    thus the system is (4), and by proposition 3, limtEP(t)=0. Finally,

    limtEP(t)=0limt˙xi=0  for  all  ilimtEG(t)=0. (32)

    We use a fourth order Runge-Kutta scheme to approximate the trajectories. The derivative is calculated as a vector in R3, and then we express this vector in spherical coordinates, ˙θtθ and ˙ϕtϕ where tθ and tϕ are the unit vectors in the direction of the azimuth angle(˙θ) and the polar angle(˙ϕ) respectively for the ith agent. {tθ,tϕ} form an orthonormal basis for the tangent space of xi. Using the angular derivatives avoids having to calculate the agent's trajectory in R3 and then project onto the sphere for every iteration which would cause significant numerical errors.

    For an agent xi=(θi,ϕi) we can write tθi and tϕi as

    tθ=(sinθcosθ0),tϕ=(cosθcosϕsinθcosϕsinϕ).

    We express the derivative of an agent as

    ˙xi=xiθi˙θ+xiϕi˙ϕ. (33)

    By direct computation, we get:

    xiθi=(sinθisinϕicosθisinϕi0),xiϕi=(cosθicosϕisinθicosϕisinϕi).

    We can also express the derivative of xi as the projection of the derivative in R3 with tθi and tϕi

    ˙xi=˙xi,tθitθi+˙xi,tϕitϕi. (34)

    It follows from (33) and (34) that

    ˙θi=1sinϕi˙xi,tθitϕi     and     ˙ϕi=˙xi,tϕitϕi. (35)

    Singularities: In (35), the factor 1sinϕ causes a singularity around ϕ=kπ for a non-negative integer k. To avoid this practical problem, before each iteration of the RK4 scheme, we identify critical agents that have a polar angle close to 0 or π (ϕ=kπ), for non-negative integer k; for these critical agents, we rotate all agents π2 around the x-axis to calculate the derivative. This is a concern for both Approach A and Approach B.

    An additional concern is singularities for agents forming consensus. In equations (3) and (7), νij is normalized, and when agent i and agent j are very close together, dividing by dP and dG causes a singularity. We avoid this problem in our simulations by defining a minimum distance dmin between agents for the sake of the normalization term in the denominator. When d(xi,xj)<dmin then

    νij=xjxj,xixidmin

    We ran simulations using different choices of Ψ to see how this choice can impact the system. We show two simulations, the first uses Approach A (Figure 7); the second uses Approach B (Figure 9). In both pictures the same interaction matrix is used (given in the appendix), and we see that our choice of Ψ (Figure 8) may dramatically change the system behavior. The second example shows the effect of curvature on the system (Figure 10) for comparison to T2 and R2 (Figure 11). Another example in the appendix shows unexpected behavior using Approach A.

    Figure 9.  A comparison of the effect of the choice of influence function for Approach B. For a=3π4 an antipodal equilibrium occurs (see Definition 2.2).
    Figure 10.  Dynamics with Approach A on S2, using the interactions matrix (36). If the agents' initial positions are close enough to each other, the agents with will form trajectories that remain in a neighborhood of their initial position.
    Figure 11.  Trajectories of three agents interacting according to the matrix A given in (36). Left: Dynamics in R2, with periodic trajectories on a unique orbit. Center: Dynamics on M=T2 with small initial mutual distances. Right: Dynamics on M=T2 with large initial distances.

    Example 4.1. Five agents with a general interaction matrix A and Ψ as defined in (26) with a{π4,π2,3π4}. The behavior of the system can change dramatically from our choice of Ψ (Figure 8), and dancing equilibria may arise from A with certain configurations, (Figure 9).

    Example 4.2. To assess the influence of the curvature of S2 on the dynamics, observe a simple case involving 3 agents evolving according to the interaction matrix:

    A=(011101110) (36)

    In Section 6.2, we prove that those dynamics in R2 lead to periodic trajectories on a single orbit shared by all three agents, the orbit's parameters being fully determined by the initial conditions (see Theorem 6.4). However, the same dynamics on the sphere do not give rise to periodic trajectories. In sections 5.3, we also discuss the dynamics with this interactions matrix on T2, to assess the effect of curvature of the manifold.

    We now study how the general dynamics given by equation (1) apply to the specific case of the torus T2R3. Let (ex,ey,ez) denote the Euclidean basis of R3. Let (R,r)(R+)2, with R>r. We define the manifold T2 as the torus obtained by rotating the circle (xR)2+z2=r2 around the z-axis. Hence T2 is defined by the equation (Rx2+y2)2+z2=r2. The parametric equations for such a torus are:

    {x=(R+rcosθ)cosϕy=(R+rcosθ)sinϕz=rsinθ for (ϕ,θ)[0,2π)2.

    The angles ϕ and θ are respectively referred to as the toroidal and poloidal angles. A set of points with the same toroidal angle is called a meridian.

    We first investigate the behavior of system (1) with Approach B (using the geodesic distance) in the case of T2. Unlike in the cases of S1 and S2 presented in the sections 3 and 4, there exists no simple expression for the geodesic distance between two points on the torus. In 1903, Bliss studied and classified the different kinds of geodesic lines on the standard torus [2], using elliptic functions. Gravesen et al. determined the structure of the cut loci of a torus of revolution [8].

    Several challenges arise when defining Approach B on T2. Firstly, computing the Riemannian distance between two points is highly non-trivial. One could consider approximating it numerically, but in the numerical discretization of equations (1)-(6)-(7), N(N1)/2 geodesics would have to be computed per time-step. That would require tremendous computing power.

    Secondly, assuming that one is able to efficiently compute the geodesics on T2, one must take into account the cut-loci of each point to ensure that the dynamics (1)-(6)-(7) are well-defined. A method to guarantee well-defined dynamics would be to use a bounded confidence model [11], where the neighborhood of influence for an agent xi at point p is of smaller radius than the closest element in the cut locus of p. See section 2 for conditions on Ψ to make the right hand side of equation (1) continuous.

    For simplicity, we thus focus on Approach A, where the dynamics are a function of the projection of each vector xjxi onto the tangent space at xi. We will show that some restrictions still apply to the interaction function Ψ, but they are less restrictive and more easily determined than in Approach B.

    Equations (1)-(2) reads:

    ˙xi=Nj=1aijΨ(ΠTxiT2(xjxi))νij,i{1,,N}. (37)

    The vector νij depends on the influence that xj has over xi. It is zero if ΠTxiT2(xjxi)=0, and it is a unit vector otherwise. Let Ni be the set of points that have no influence on xi (see Table 1). Then, given i,j{1,,N}, νij has the following expression:

    νij={ΠTxiT2(xjxi)ΠTxiM(xjxi) if xjNi0 if xjNi. (38)

    Let xiT2. We start by determining the set Ni. For all i, we define the vectors uϕi=cosϕiex+sinϕiey and uθi=cosθiuϕi+sinθiez, so that each agent's position vector reads: xi=Ruϕi+ruθi. With these notations, uθi is the normal to the tangent space at the point xi. A basis for the tangent space at a point xi(ϕi,θi) is given by the two tangent vectors tϕi=(sinϕi,cosϕi,0) and tθi=(sinθicosϕi,sinθisinϕi,cosθi). Notice that xi,tϕi=0. Hence the condition ΠTxiT2(xjxi)=0 reads:

    {xj,tϕi=0xjxi,tθi=0.

    After computations, we get:

    xj,tϕi=0sin(ϕjϕi)=0ϕj=ϕi+kπ,kZ.

    If ϕj=ϕi, the second condition becomes:

    xjxi,tθi=0sin(θjθi)=0θj=θi+kπ,kZ.

    If ϕj=ϕi±π, the second condition becomes:

    sin(θi+θj)=2Rrsinθi.

    Notice that this last equation only has a solution if |sinθi|r2R. The set of positions that have no influence on xi thus comprises up to four points on the torus, depending on the values of r, R and sinθi. We then have: Ni={(ϕi,θi), (ϕi,θi), (ϕi,θisgn(sinθi)arcsin(|2Rrsinθi|)),(ϕi,πθi+sgn(sinθi)arcsin(|2Rrsinθi|))}. To ensure the continuity of the right-hand side of equation (37), one must impose the conditions of table 1.

    We now go back to equation (37). We study the specific case where ΨId, which indeed satisfies (1). Then the system becomes:

    ˙xi=ΠTxiT2(Nj=1aij(xjxi)). (39)

    Hence the velocity reads:

    ˙xi=Nj=1aij(xjxi)Nj=1aij(xjxi),uθiuθi=αiαi,uθiuθi(Nj=1aij)xi,tθitθi

    where αi:=Nj=1aijxj is the sum of the influences of all agents on agent i. Notice that with the same notation, the system does not reduce to the simple form ˙xi=αiαi,xixi for the same dynamics on the sphere (see [3]). This is due to the fact that on the torus, the position vector xi does not define the normal to the tangent space at xi, unlike in the cases of S1 and S2.

    The velocity of each agent is given by:

    ˙xi=(˙ϕisinϕi(R+rcosθi)r˙θisinθicosϕi˙ϕicosϕi(R+rcosθi)r˙θisinθisinϕir˙θicosθi)=˙ϕi(R+rcosθi)tϕi+r˙θitθi. (40)

    From (39) and (40) we get the angular velocities:

    {˙ϕi=1(R+rcosθi)Nj=1aij(xjxi),tϕi˙θi=1rNj=1aij(xjxi),tθi. (41)

    Notice that unlike in the case of S2, see equation (35) here the derivatives ˙ϕi and ˙θi are not singular. This makes numerical simulations straightforward, not requiring the approximations described in Section 4.2.

    We now analyze the dynamics (1)-(2)-(3) on T2. We identify families of initial conditions that trivialize the dynamics.

    Proposition 4. Consider the dynamics (1)-(2)-(3) on M=T2. Let Pz:={(x,y,z)R3|z=0}. Let xi(t) be the position of the ithe agent at time t. If for all i{1,,N}, xi(0)T2Pz, then for all t0, for all i{1,,N}, xi(t)T2Pz.

    Proof. Suppose that for all i{1,,N}, xi(0)T2Pz. Then for all i{1,,N}, θi(0)=0 or θi(0)=π. Hence, for all i,j{1,,N},

    tθi(0)=(00±π)andxj(0)xi(0)=((R+rcosθj)cosϕj(R+rcosθi)cosϕi(R+rcosθj)sinϕj(R+rcosθi)sinϕi0)

    From equation (41) we get: for alli{1,,N}, ˙θi=0. By uniqueness of solution, for all i{1,,N}, θi(t)=θi(0). All the initial velocities belong to the plane Pz. Hence all agents remain on Pz at all time.

    Remark 5. As a consequence of Proposition 4, if all agents are initially in T2Pz, all agents initially on the bigger circle θ=0 remain on the major circle at all time and all agents on the minor circle θ=π remain on the minor circle at all time. In particular, if all agents are initially all on the same circle (i.e. for all i{1,,N}, θi=0 or for all i{1,,N}, θi=π), then the torus dynamics simplify to the dynamics on S1 given by (22) or (23).

    Proposition 5. Consider the dynamics (1)-(2)-(3) on M=T2. Let ˜ϕ[0,2π] and let P˜ϕ:={(x,y,z)R3|y=tan(˜ϕ)x}. If for all i{1,,N}, xi(0)T2P˜ϕ, then for all t0, for all i{1,,N}, xi(t)T2P˜ϕ.

    Proof. Suppose without loss of generality that ˉϕ=0. Similarly to the proof for Proposition 4, we can show that for all i{1,,N}, ˙ϕi(0)=0. By uniqueness of solution, for all i{1,,N}, ϕi(t)=ϕi(0). Hence all agents remain in P˜ϕ at all time.

    Remark 6. As a consequence of Proposition 5, if all agents are initially in T2Pˉϕ, all agents initially on the circle ϕ=ˉϕ remain on that circle at all time and all agents on the circle ϕ=ˉϕ remain on that circle at all time. In particular, if all agents are initially all on the same circle (i.e. for alli{1,,N}, ϕi=ˉϕ or for alli{1,,N}, ϕi=ˉϕ), then the torus dynamics simplify to the dynamics on S1 given by (22) or (23).

    To assess the influence of the curvature of the manifold on the dynamics, we compare a simple case involving 3 agents evolving according to the interaction matrix given in equation (36). As in the case of S2, the dynamics on the torus do not give rise to periodic trajectories (as opposed to the dynamics in R2, see Theorem 6.4). Instead, since T2 can locally be identified with R2, if the initial mutual distances are small enough, the dynamics resemble those in R2. More specifically, the trajectories are quasi-periodic with a gradual shift of the center of mass (see Figure 11). However, if the initial distances between agents are large, the geometry and curvature of the torus changes radically the behavior of the system.

    As seen in Sections 3 and 5.3, when the interaction matrix A satisfies certain properties, for instance given by (28) on S1 or by (36) in R2, then the trajectories exhibit special properties of symmetry or periodicity. In [3], configurations on S2 in which all mutual distances between agents remain constant were named dancing equilibrium.

    In this section, we investigate systems with similar properties of periodicity or symmetry. We use the term social choreography, drawing a parallel with the well-known "n-body choreographies" discovered by Moore [13,14] in the context of point masses subject to gravitational forces. In the n-body problem, the interaction potentials between masses are predetermined, as they depend exclusively on the masses and distances between agents. Hence the conditions for a n-body choreography to occur only depend on the initial state of the system. In the case of social choreography, there are more degrees of freedom, as we design the interaction matrix as well as to set the initial conditions.

    We study sufficient conditions on the interaction matrices for the trajectories of the system to be periodic or symmetric by focusing on the Euclidean space R2, with the specific choice of interaction potential ΨId. A future direction of this paper can consist in extending these results to general Riemannian manifolds. In R2 and with ΨId, both approaches A and B are equivalent and the system simply reads as:

    for alli{1,,N},˙xi=Nj=1aij(xjxi). (42)

    We define the kinetic energy E=EG=EP as in Definition 2.3:

    E(t):=12Ni=1˙xi(t)2. (43)

    A simple case of social choreography is that of a system with periodic trajectories, which we define as follows:

    Definition 6.1. Let (xi)i=1...N be a solution of (42). We refer to the system as having periodic trajectories if there exists τ>0 such that

    for alli{1,...,N},for allt>0,     xi(t+τ)=xi(t).

    We will examine possible periodic behaviors of the system in sections 6.2, 6.3 and 6.4.

    We now give sufficient conditions on the interaction matrix and on the initial conditions for the system to be invariant by rotation.

    Theorem 6.2. Let kN such that k divides N. Let Pk=(0INkIk0) be the matrix of change of basis from (e1,,eN) to (ek,,eN,e1,,ek1). Let R(θ) denote the rotation matrix in R2 for the angle θ[0,2π). Suppose that initially, the system is invariant by rotation of angle 2kπN, that is:

    foralli{1,,N},R(2kπN)xi(0)={xi+k(0)ifi+kNxi+kN(0)ifi+k>N.

    Suppose that the interaction matrix A is invariant by change of basis, i.e. P1kAPk=A. Then the system remains invariant by rotation of angle 2kπN at all time:

    forallt>0,foralli{1,,N},R(2kπN)xi(t)={xi+k(t)ifi+kNxi+kN(t)ifi+k>N.

    Proof. Let AMN(R) be the interaction matrix, i.e. A=(aij)i,j=1,N, and define D=(jaij). Let x=(x1,,xN) denote the set of all xi's. It is a vector of length N with entries in R2. Let XMN×2(R) denote the corresponding matrix of RN×2 such that for all i{1,,N}, for all j{1,2},Xij is the j-th coordinate of xi. With these notations, ˙X=˜AX, where ˜A=AD. We denote by (e1,,eN) the canonical orthonormal basis of (R)N such that X=Ni=1eixTi.

    From the definition of the matrix X, the condition

    for alli{1,,N},R(2kπN)xi(0)={xi+k(0)ifi+kNxi+kN(0)ifi+k>N

    can be rewritten as: PkX(0)=(R(2kπN)X(0)T)T. Let Y:=PkX and Z:=(R(2kπN)XT)T. From the theorem's hypotheses, Y(0)=Z(0). Let us show that Y and Z have the same evolution. One can easily prove that P1k˜APk if and only if P1kAPk. Then notice that

    ˙X=˜AX=P1k˜APkX.

    From that we compute:

    ˙Y=Pk˙X=Pk(P1k˜APkX)=˜APkX=˜AY.

    Similarly,

    ˙Z=(R(2kπN)˙XT)T=(R(2kπN)(˜AX)T)T=(R(2kπN)XT˜AT)T=˜AZ.

    Since Y and Z satisfy the same differential equation and Y(0)=Z(0), then Y(t)=Z(t) for all t0. This implies that at all time,

    for alli{1,,N},R(2kπN)xi(t)={xi+k(t)ifi+kNxi+kN(t)ifi+k>N.

    Another example of social choreography is that of a system in which all agents share one unique orbit. Such choreographies have been discovered in the context of the n-body problem, for instance the "figure 8" orbit for three equal masses [13].

    Definition 6.3. Let (xi)i=1...N be a solution of (42). We say that the system has a unique orbit if the orbits of all points are identical, i.e.

    for alli,j{1,...,N},{zM|t>0,xi(t)=z}={zM|t>0,xj(t)=z}.

    To illustrate Theorem 6.2, we study the evolution of N agents initially positioned at regular intervals on a circle, with an interaction matrix and initial conditions given by:

    A=(01001100000110010)and for alli{1,,N},xi(0)=(cos(2iπN)sin(2iπN)). (44)

    Notice that ˜A=A, and the system satisfies the conditions of Theorem 6.2 with k=1. Hence for all i{1,,N1}, R(2πN)xi(t)=xi+1(t) and R(2πN)xN(t)=x1(t). The 2N-dimensional system then reduces to a 2-dimensional one for the two coordinates x11 and x12 of x1, and all the other variables can be recovered by rotation of x1:

    ˙x1=x2xN=R(2πN)x1R(2πN)x1.

    This can be written as:

    (˙x11˙x12)=(02sin(2πN)2sin(2πN)0)(x11x12).

    Solving this linear system yields:

    {x11(t)=x11(0)cos(2sin(2πN)t)x12(0)sin(2sin(2πN)t)=cos(2sin(2πN)t)x12(t)=x11(0)sin(2sin(2πN)t)+x12(0)cos(2sin(2πN)t)=sin(2sin(2πN)t).

    This proves that all agents share one common circular orbit, and their trajectories are periodic of period 2π(2sin(2πN))1. Figure 15 provides a numerical illustration of this behavior, with 10 agents initially positioned at regular intervals on the unit circle.

    Figure 12.  Evolutions of the coordinates of the three agents evolving on T2 with interaction matrix A from equation (36), with small initial mutual distances. Left: Evolution of ϕ. Center: Evolution of θ. Right: Evolution of the kinetic energy.
    Figure 13.  Left: Evolution of 12 agents with the conditions of Theorem 6.2, with k=3, resulting in diverging trajectories. Dark to light color scale indicates earlier to later time. Right: corresponding exploding kinetic energy. The interaction matrix A and the initial positions were generated according to a random algorithm, with the conditions of Theorem 6.2.
    Figure 14.  Left: Evolution of 12 agents with the conditions of Theorem 6.2, with k=3, resulting in convergence to consensus. Dark to light color scale indicates earlier to later time. Right: corresponding kinetic energy converging to zero. The interaction matrix A and the initial positions were generated according to a random algorithm, with the conditions of Theorem 6.2.
    Figure 15.  Evolution of 10 agents with initial conditions and interaction matrix given in (44). The agents have periodic trajectories along one shared circular orbit.

    Another interesting example is that of 3 agents interacting according to the interaction matrix given previously, which, reduced to N=3, gives:

    A=(011101110). (45)

    Theorem 6.4. Let N=3. Consider the system (42) with interaction matrix given by (45). Then there exists a unique orbit shared by all agents, and all three trajectories are periodic.

    Proof. The x and y-coordinates of the systems are decoupled, so that the 6-dimensi-onal system can be reduced to two 3-dimensional ones. Notice that ˜A=A. Then for each coordinate j{1,2}, the system reads:

    (x1jx2jx3j)(t)=exp(tA)(x01jx02jx03j)

    with

    etA=13(1+2cos(3t)1cos(3t)+3sin(3t)1cos(3t)3sin(3t)1cos(3t)3sin(3t)1+2cos(3t)1cos(3t)+3sin(3t)1cos(3t)+3sin(3t)1cos(3t)3sin(3t)1+2cos(3t)).

    Due to the special structure of etA, this can be rewritten as:

    (x1jx2jx3j)(t)=13(x01jx02jx03jx02jx03jx01jx03jx01jx02j)(1+2cos(3t)1cos(3t)+3sin(3t)1cos(3t)3sin(3t)).

    This shows that all three trajectories are periodic, or period 2π3. One can compute the positions of each agent after a third of a period and notice that:

    (x1jx2jx3j)(t+2π33)
    =13(x01jx02jx03jx02jx03jx01jx03jx01jx02j)(1cos(3t)3sin(3t)1+2cos(3t)1cos(3t)+3sin(3t))=(x2jx3jx1j)(t).

    This shows that there is one unique shared orbit.

    Other conditions on the interaction matrix A give rise to different kinds of periodic behaviors. Here we provide sufficient conditions for the system to exhibit periodic trajectories, such that each orbit is shared by two agents.

    Theorem 6.5 (Coupled periodic trajectories). Let N be even. Suppose that initially, the system is invariant by rotation of angle 4πN, that is:

    foralli{1,,N},R(4πN)xi(0)={xi+2(0)ifi+2Nxi+2N(0)ifi+2>N.

    Let a,b>0 and let

    A=(0a00ba0b00ba00ab00a0). (46)

    Then the system is periodic of period τ=πabsin(2π/N). Furthermore, if N is divisible by 4, opposite agents share orbits two by two, i.e.:

    forallt>0,foralli{1,,N2},xi(t+τ)=xi+N2(t),

    and the kinetic energy is periodic with period τ/2.

    Proof. First remark that the system satisfies the hypotheses of Theorem 6.2, so

    for allt>0,for alli{1,,N},R(4πN)xi(t)={xi+2(t)ifi+2Nxi+2N(t)ifi+2>N.

    Hence the system is entirely known from the positions of the first two agents, since all others can be obtained by simple rotations. We show that this 2N-dimensional problem can be rewritten as a 4-dimensional one. Indeed, using the fact that xN=R(4π/N)x2 and x3=R(4π/N)x1, the system

    {˙x1=a(x2x1)b(xNx1)˙x2=b(x3x2)a(x1x2)

    becomes:

    {˙x1=(˙x11˙x12)=a[(x21x22)(x11x12)]b[(cos(4πN)sin(4πN)sin(4πN)cos(4πN))(x21x22)(x11x12)]˙x2=(˙x21˙x22)=b[(cos(4πN)sin(4πN)sin(4πN)cos(4πN))(x11x12)(x21x22)]a[(x11x12)(x21x22)].

    This can be rewritten in matrix form as:

    (˙x11˙x12˙x21˙x22)=A4(x11x12x21x22) (47)

    where

    A4:=(a+b0abcos(4πN)bsin(4πN)0a+bbsin(4πN)abcos(4πN)a+bcos(4πN)bsin(4πN)ab0bsin(4πN)a+bcos(4πN)0ab)(x11x12x21x22).

    One can easily show that this reduced interaction matrix A4 has two purely imaginary conjugate eigenvalues, iλ and iλ, each of multiplicity 2, where λ=2absin(2πN). Hence the solution of the system (47) can be written as a weighted sum of the functions tcos(λt) and tsin(λt). This implies that the system is periodic, of period

    τ=2πλ=πabsin(2πN).

    Furthermore, if N is divisible by 4, according to Theorem 6.2, xN2+1=x1 and xN2+2=x2. This implies that for all t>0, x1(t+τ)=x1(t)=xN2+1(t) and x2(t+τ)=x2(t)=xN2+2(t), so the agents x1 and xN2+1 share an orbit, as well as all pairs of agents xi and xN2+i for i{1,,N2}.

    As a consequence, the kinetic energy is periodic, of period τE=τ=π/(absin(2πN)). If N is divisible by 4, every half period, the system is rotated by an angle π, so the kinetic energy is periodic with period τE=τ/2.

    Remark 7. Notice that the agents sharing orbits do not interact with one another, as shown in Figure 16.

    Figure 16.  Left: Directed graph corresponding to the matrix A given in (44). Full arrows represent positive coefficients (aij>0) while dashed ones represent negative coefficients (aij<0). Right: Weighted directed graph corresponding to the matrix A given in (46). Thin arrows represent the weighted edges |aij|=a while bold ones represent the weight |aij|=b. Nodes with the same color and symbol share orbits but are not directly connected in the graph.

    An example of such a choreography is given in Figure 17.

    Figure 17.  Left: Periodic trajectories of 8 agents sharing orbits two by two, in the situation of Theorem 6.5. Matrix A from (46) was constructed with (a,b)=(1,3). The initial positions x1(0) and x2(0) were randomly generated and the other 6 were obtained by rotation. The period is τ=2π/6. Right: Corresponding kinetic energy, of period τ/2.

    Remark 8. As a slight generalization, we provide numerical simulations illustrating a similar behavior, but with slightly different conditions: the periodic evolution of 9 agents on three distinct orbits shared three by three, see figures 18 and 19.

    Figure 18.  Left: evolution of 9 agents with periodic trajectories, each orbit shared by 3 agents. Right: periodic kinetic energy.
    Figure 19.  Isolated orbits of the evolution shown in Figure 18. Left: trajectories of agents 3, 6, 9. Middle: trajectories of agents 1, 4, 7. Right: trajectories of agents 2, 5, 8).

    In sections 6.2 and 6.3, we provided conditions for the trajectories of the system to be periodic. Here, we explore further the notion of periodicity by studying systems with drift, displaying helical trajectories but periodic kinetic energy.

    Definition 6.6. Let (xi)i=1...N be a solution of (42). We call the corresponding trajectories helical trajectories if there exists vR2 and τR such that

    for alli{1,...,N},for allt>0,xi(t+τ)=xi(t)+τv.

    Notice that this definition generalizes the notion of periodic trajectories recalled in Definition 6.1, which corresponds to the case v=0. When v0, the system has a drift term, meaning that the relative positions between agents remain periodic but their absolute positions evolve in space.

    Theorem 6.7. Sufficient conditions for helical trajectories. Let N=4. Let (a,b,c,d)(R+)4 such that the interaction matrix reads

    A=(0a0da0b00b0cd0c0). (48)

    Then the system exhibits helical trajectories.

    Proof. First notice that the first and second components xi1 and xi2 of the i-th agent's position are decoupled, so that the system in matrix form reads

    ˙xj=(˙x1j˙x2j˙x3j˙x4j)=(daa0daabb00bbccd0ccd)(x1jx2jx3jx4j):=˜A(x1jx2jx3jx4j),for j{1,2}. (49)

    Hence the projections of x on the first and second axes solve the same differential equation. The matrix ˜A has three distinct eigenvalues:

    λ1=0,    λ2=i(a+c)(b+d)    and     λ3=i(a+c)(b+d).

    There is one eigenvector associated with λ1: v1:=(1,1,1,1)T. One can show that the vectors x(t)=v1 and x(t)=v1t+ν are both solutions of System (49), where, denoting Δ:=bcdabc+abdacd,

    ν:=1Δ(ab+bc+Δ,abcd+Δ,ab+ad+Δ,Δ)T.

    Let v2 denote the eigenvector associated with λ2 and let vR2 and vI2 denote respectively its real and imaginary components, i.e. v2:=vR2+ivI2. Then the solution of System (49) can be written as:

    xj(t)=Cj1v1+Cj2(v1t+ν)+Cj3[vR2cos(λ2t)vI2sin(λ2t)]           +Cj4[vR2sin(λ2t)+vI2cos(λ2t)]

    where (C1,C2,C3,C4)R4 are constants depending on the initial conditions. Let τ=2πλ2. Then for all t>0, for all i{1,,4}, for all j{1,2}, xij(t+τ)=xij(t)+Cj2τ. This can be rewritten as:

    for alli{1,,4},for allt>0,xi(t+τ)=xi(t)+(C12C22)τ.

    Theorem 6.8. A system with helical trajectories has periodic kinetic energy.

    Proof. Supose that (xi)i=1N has helical trajectories, i.e. there exists τR, vR2 such that for all i{1,,N}, for all t0, xi(t+τ)=xi(t)+τv. Then ˙xi(t+τ)=˙xi(t) and so E(t+τ)=E(t).

    An example using Approach A which shows unexpected behavior in the first example (A.1), as well as the interactions matrix and initial positions used for simulations shown in Figure 7 and 9.

    Example 7.1. 15 agents with a general interaction matrix A. We use the same notion of a general interactions matrix as used in [3]. The interactions matrix A is composed of integers aij that are uniformly chosen between -5 and 5 inclusive. ΨId. Generally, a system with this kind of interaction matrix will exhibit simple oscillating kinetic energy, as in [3]. We show a rare simulation using this general interactions matrix A in the appendix (Figure 22 and 23).

    A=(5110415110205124153441221)X=(1.57551.73991.65230.56195.30262.70082.49710.72880.65710.52811.21800.08402.28121.01295.49491.74413.76852.59031.62182.52662.25210.07675.57661.16715.65821.54532.81461.46411.68920.1310)
    Figure 20.  Left: Trajectories of 4 agents with helical trajectories. Parameters for matrix A (48) chosen to be (a,b,c,d)=(1,2,3,4). Dark to light color indicates earlier to later time. Right: Corresponding kinetic energy. The period is τ=2π((a+c)(b+d))1/2=π/6 (see proof of Theorem 6.7).
    Figure 21.  Evolution of the first and second coordinates of 4 agents with helical trajectories.
    Figure 22.  Energy of the system using Approach A, 15 agents, and a general interaction matrix (left). A snapshot of the energy oscillations to match with trajectories in Figure 23 (right).
    Figure 23.  An agent's trajectory simulated with Approach A, 15 agents, and a general interaction matrix. The trajectory in shown the top right is of a second agent. The agents oscillate with amplitudes that increase with time, eventually the trajectory approximates a great circle, after which the oscillations resume with smaller amplitudes.

    Example 7.2. Five agents with a general interaction matrix A and Ψ as defined in (26). Simulations are shown in Figure 7 and 9 with a{π4,π2,3π4}.

    A=(5110415110205124153441221)X=(6.17432.84734.58832.76352.16062.56913.66980.81910.6771138672)
    [1] A. Aydoğdu, M. Caponigro, S. McQuade, B. Piccoli, N. Pouradier Duteil, F. Rossi and E. Trélat, Interaction network, state space and control in social dynamics, in Active Particles Volume 1, Theory, Methods, and Applications (eds. N. Bellomo, P. Degond and E. Tadmor), Birkhäuser-Springer, 2017.
    [2] The geodesic lines on the anchor ring. Annals of Mathematics (1902) 4: 1-21.
    [3] A nonlinear model of opinion formation on the sphere. Discrete and Continuous Dynamical Systems -Series A (2014) 35: 4241-4268.
    [4] J. Cheeger and D. G. Ebin, Comparison Theorems in Riemannian Geometry, Vol. 365, AMS Chelsea Publishing, 1975.
    [5] D. Chi, S. -H. Choi and S. -Y. Ha, Emergent behaviors of a holonomic particle system on a sphere, Journal of Mathematical Physics, 55 (2014), 052703, 18pp. doi: 10.1063/1.4878117
    [6] Effects of anisotropic interactions on the structure of animal groups. Journal of Mathematical Biology (2011) 62: 569-588.
    [7] Synchronization in complex oscillator networks and smart grids. Proceedings of the National Academy of Sciences (2013) 110: 2005-2010.
    [8] J. Gravesen, S. Markvorsen, R. Sinclair and M. Tanaka, The Cut Locus of a Torus of Revolution, Technical University of Denmark, Department of Mathematics, 2003.
    [9] Emergent behavior of a Cucker-Smale type particle model with nonlinear velocity couplings. IEEE Transactions on Automatic Control (2010) 55: 1679-1683.
    [10] Collective synchronization of classical and quantum oscillators. EMS Surveys in Mathematical Sciences (2016) 3: 209-267.
    [11] R. Hegselmann and U. Krause, Opinion dynamics and bounded confidence: Models, analysis and simulation, Journal of Artificial Societies and Social Simulation 5 (2002).
    [12] Cooperative dynamics of oscillator community a study based on lattice of rings. Progress of Theoretical Physics Supplement (1984) 79: 223-240.
    [13] Braids in classical dynamics. Physical Review Letters (1993) 70: 3675-3679.
    [14] New periodic orbits for the n-body problem. ASME. J. Comput. Nonlinear Dynam. (2006) 1: 307-311.
    [15] Heterophilious dynamics enhances consensus. SIAM Review (2014) 56: 577-621.
    [16] Coordinated motion design on lie groups. Automatic Control, IEEE Transactions on (2010) 55: 1047-1058.
    [17] Consensus optimization on manifolds. SIAM Journal on Control and Optimization (2009) 48: 56-76.
    [18] Synchronization and balancing on the N-torus. Sys. Cont. Let. (2007) 56: 335-341.
    [19] Stabilization of planar collective motion: All-to-all communication. Automatic Control, IEEE Transactions on (2007) 52: 811-824.
    [20] Stabilization of planar collective motion with limited communication. Automatic Control, IEEE Transactions on (2008) 53: 706-719.
    [21] P. Sobkowicz, Modelling opinion formation with physics tools: Call for closer link with reality, Journal of Artificial Societies and Social Simulation, 12 (2009), p11.
    [22] From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators. Physica D: Nonlinear Phenomena (2000) 143: 1-20.
  • This article has been cited by:

    1. Johan Markdahl, 2020, Consensus seeking gradient descent flows on boundaries of convex sets, 978-1-5386-8266-1, 830, 10.23919/ACC45564.2020.9147648
    2. Seung-Yeal Ha, Doheon Kim, Franz Wilhelm Schloder, Emergent Behaviors of Cucker–Smale Flocks on Riemannian Manifolds, 2021, 66, 0018-9286, 3020, 10.1109/TAC.2020.3014096
    3. Johan Markdahl, Counterexamples in synchronization: Pathologies of consensus seeking gradient descent flows on surfaces, 2021, 134, 00051098, 109945, 10.1016/j.automatica.2021.109945
    4. M A Lohe, Higher-dimensional generalizations of the Watanabe–Strogatz transform for vector models of synchronization, 2018, 51, 1751-8113, 225101, 10.1088/1751-8121/aac030
    5. Seung-Yeal Ha, Seungsu Hwang, Dohyun Kim, Sun-Chul Kim, Chanho Min, Emergent behaviors of a first-order particle swarm model on the hyperboloid, 2020, 61, 0022-2488, 042701, 10.1063/1.5066255
    6. M.A. Lohe, WITHDRAWN: Solitons in complex systems of chiral fields with Kuramoto interactions, 2020, 25900544, 100048, 10.1016/j.csfx.2020.100048
    7. Johan Markdahl, Synchronization on Riemannian Manifolds: Multiply Connected Implies Multistable, 2021, 66, 0018-9286, 4311, 10.1109/TAC.2020.3030849
    8. Benedetto Piccoli, Nastassia Pouradier Duteil, 2021, Chapter 12, 978-3-030-82945-2, 289, 10.1007/978-3-030-82946-9_12
    9. Razvan C. Fetecau, Seung-Yeal Ha, Hansol Park, Emergent Behaviors of Rotation Matrix Flocks, 2022, 21, 1536-0040, 1382, 10.1137/21M1404569
    10. Mohamad Amin Sharifi Kolarijani, Anton V. Proskurnikov, Peyman Mohajerin Esfahani, Macroscopic Noisy Bounded Confidence Models With Distributed Radical Opinions, 2021, 66, 0018-9286, 1174, 10.1109/TAC.2020.2994284
    11. Seung-Yeal Ha, Dohyun Kim, Hansol Park, Existence and Emergent Dynamics of Quadratically Separable States to the Lohe Tensor Model, 2022, 21, 1536-0040, 1166, 10.1137/21M1409664
    12. Nathalie Ayi, Nastassia Pouradier Duteil, Mean-field and graph limits for collective dynamics models with time-varying weights, 2021, 299, 00220396, 65, 10.1016/j.jde.2021.07.010
    13. Johan Markdahl, Daniele Proverbio, La Mi, Jorge Goncalves, Almost global convergence to practical synchronization in the generalized Kuramoto model on networks over the n-sphere, 2021, 4, 2399-3650, 10.1038/s42005-021-00689-y
    14. M. A. Lohe, Solitons in complex systems of chiral fields with Kuramoto interactions, 2021, 31, 1054-1500, 023138, 10.1063/5.0039991
    15. Hyunjin Ahn, Seung-Yeal Ha, Woojoo Shim, Emergent dynamics of a thermodynamic Cucker-Smale ensemble on complete Riemannian manifolds, 2021, 14, 1937-5077, 323, 10.3934/krm.2021007
    16. Dohyun Kim, Jeongho Kim, Aggregation and disaggregation of active particles on the unit sphere with time-dependent frequencies, 2022, 27, 1531-3492, 2247, 10.3934/dcdsb.2021131
    17. Amic Frouvelle, Jian-Guo Liu, 2019, Chapter 16, 978-3-030-15095-2, 457, 10.1007/978-3-030-15096-9_16
    18. Hyungjin Huh, Dohyun Kim, Asymptotic behavior of gradient flows on the unit sphere with various potentials, 2021, 270, 00220396, 47, 10.1016/j.jde.2020.07.016
    19. Johan Markdahl, Johan Thunberg, Jorge Goncalves, High-dimensional Kuramoto models on Stiefel manifolds synchronize complex networks almost globally, 2020, 113, 00051098, 108736, 10.1016/j.automatica.2019.108736
    20. Sean McQuade, Benedetto Piccoli, Nastassia Pouradier Duteil, Social dynamics models with time-varying influence, 2019, 29, 0218-2025, 681, 10.1142/S0218202519400037
    21. Christian Düll, Piotr Gwiazda, Anna Marciniak-Czochra, Jakub Skrzeczkowski, Measure differential equation with a nonlinear growth/decay term, 2023, 73, 14681218, 103917, 10.1016/j.nonrwa.2023.103917
    22. Sui Tang, Malik Tuerkoen, Hanming Zhou, On the Identifiability of Nonlocal Interaction Kernels in First-Order Systems of Interacting Particles on Riemannian Manifolds, 2024, 84, 0036-1399, 2067, 10.1137/23M1599835
    23. Lunxiao Tang, Tao Yu, Maokang Luo, Aggregation of Multi-Agent System on Manifolds With Riemannian Geometry, 2023, 11, 2169-3536, 55456, 10.1109/ACCESS.2023.3280070
    24. La Mi, Jorge Gonçalves, Johan Markdahl, Polarization of Multi-agent Gradient Flows Over Manifolds With Application to Opinion Dynamics, 2024, 69, 0018-9286, 1288, 10.1109/TAC.2023.3281980
    25. Piotr Lisowski, Roman Urban, On Some Dynamics in Conceptual Spaces, 2024, 0925-8531, 10.1007/s10849-024-09422-8
    26. Wouter Jongeneel, Emmanuel Moulay, 2023, Chapter 6, 978-3-031-30132-2, 77, 10.1007/978-3-031-30133-9_6
    27. Benedetto Piccoli, Control of multi-agent systems: Results, open problems, and applications, 2023, 21, 2391-5455, 10.1515/math-2022-0585
    28. Shanshan Peng, Jianquan Lu, Almost global synchronization of amplitude-dependent high-dimensional Kuramoto model, 2025, 471, 01672789, 134448, 10.1016/j.physd.2024.134448
    29. Zhengyang Qiao, Yicheng Liu, Xiao Wang, Consensus and bipartite consensus in graphon models for opinion dynamics on the sphere, 2025, 472, 01672789, 134503, 10.1016/j.physd.2024.134503
    30. Animesh Chakravarthy, Debasish Ghose, 2024, Steering Opinions over n-Dimensional Spherical Manifolds, 979-8-3503-1633-9, 3057, 10.1109/CDC56724.2024.10886840
  • Reader Comments
  • © 2017 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5678) PDF downloads(252) Cited by(30)

Figures and Tables

Figures(23)  /  Tables(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog