Synchronization and resonance on networks are some of the most remarkable collective dynamical phenomena. The network topology, or the nature and distribution of the connections within an ensemble of coupled oscillators, plays a crucial role in shaping the local and global evolution of the two phenomena. This article further explores this relationship within a compact mathematical framework and provides new contributions on certain pivotal issues, including a closed bound for the average synchronization time in arbitrary topologies; new evidences of the effect of the coupling strength on this time; exact closed expressions for the resonance frequencies in terms of the eigenvalues of the Laplacian matrix; a measure of the effectiveness of an influencer node’s impact on the network; and, finally, a discussion on the existence of a resonant synchronized state. Some properties of the solution of the linear swing equation are also discussed within the same setting. Numerical experiments conducted on two distinct real networks—a social network and a power grid—illustrate the significance of these results and shed light on intriguing aspects of how these processes can be interpreted within networks of this kind.

Why do pink flamingos perform a pre-breeding dance in perfect unison? How does the entire Australian coral reef manage to synchronize to release millions of gametes at the same time each year? How does it happen that a number of individuals in a social network take part with increasing involvement in collective ritual events? These questions concern the problem of how multiple agents within a network can achieve an identical synchronized dynamic behavior, possibly with an amplification of the intensity of the phenomenon over time. Prompted by the ubiquity and transversality of these processes, the paper offers some novel insights into the already immense research field on the dynamics of networks of coupled identical harmonic oscillators.

Up to several thousand great flamingos, during the pre-breeding period, gather in southern France to perform in synchrony a variety of stereotyped movements. Such group performances, involving elaborate courtship rituals and synchronized dances, are believed to play a role in enabling reproductive synchronization.1 In the Australia’s Great Barrier Reef, many coral species tune their reproductive cycles with incredible precision, synchronizing each year all on the same day, just after the full moon in November, to collectively release millions of gametes in one of the largest scale mass reproduction events on the planet.2 In the epidemiological context, synchronization processes govern the dynamics of disease spread through complex networks in anthropic environments. Examples are the dynamics of childhood diseases, such as measles, mainly determined by the annual seasonal cycle, or the pulse vaccination-induced synchronization in simple epidemic models. It is a peculiar feature that ecological communities tend to adjust to periodic environmental or climatic driving variables. Although ecologists and epidemiologists have essentially opposite goals, the mathematical laws underlying the different population dynamics are very similar and frequently involve internal synchronization or resonance with external driving factors.3 We find countless examples in different settings. In the biomechanical context, Brumley et al.4 study a mechanism of self-propulsion in ancestral eukaryotes, unicellular organisms equipped with tiny mobile flagella that allow them to move in a liquid medium. Correlated changes in the directions of the flagellar beating can be explained by the raise of some kind of large scale coordination and synchronization, even in a weakly coupled system as is the case for the green alga Volvox. In the human social setting, the emergence of synchronization processes among individuals, on a physical or mental level, can be induced by collective dynamics and participation in repetitive rituals or dances, as occurs in Zikr, a Sufi mystical dance performed by Chechen Sufis and involving hundreds of men at once. However, it is well known that synchronization phenomena do not always have positive implications. Strong synchronization may indeed be related to pathological activities as is the case for the epileptic seizures and Parkinson’s disease. Hammond et al.5 have observed an abnormally synchronized oscillatory activity at multiple levels of the basal ganglia-cortical loop in humans with Parkinson’s disease. This excessive synchronization correlates with motor deficit, and its suppression by specific therapies could suggest an effective approach to improve the living conditions of patients with severe motility disorders.

These examples show that synchronization emerges, in general, from the collaboration, competition, and interaction of multiple agents when they undergo periodic and cyclic dynamics. A first interpretive model, albeit oversimplified, reduces such systems to a set of coupled identical harmonic oscillators all synchronizing on the same frequency.

The two pioneering papers on synchronization in coupled systems6 and synchronization in chaotic systems7 have stimulated a great deal of interest in the study of complete synchronization of coupled nonlinear dynamical processes. The first paradigmatic model of diffusively coupled identical oscillators on a complex network is given in Pecora et al.,8 where the master stability function is introduced for the first time. Such a function accomplishes two goals: decrease the computational load and provide a guidance in designing coupling configurations that conform to the required stability. This approach is currently the most widely applied one to determine the stability of the fully synchronized state in terms of the eigenstructure of the connectivity matrix.

The interest in the study of synchronization in coupled systems grew considerably after the discovery of the small-world and scale-free properties of many natural and artificial complex networks. Barahona and Pecora9 analyze the implications of the small-world phenomenon on the synchronization of oscillator networks of arbitrary topology. Subsequent analysis provided evidence that well-connected networks are easier to globally synchronize than regular networks. For instance, Chen et al.10 found that synchronization processes typically start first from a small fraction of hub nodes and then spread to nodes with smaller degrees. This is replicated at the cluster level, where a typical synchronization process starts from partial synchronization within clusters to evolve to global complete synchronization.

Ren11 analyzes coupled second-order linear harmonic oscillators, described as point masses, under rather mild network connectivity assumptions, showing that their positions and velocities can be synchronized in the presence or absence of a leader. A similar approach is used in Zhan et al.12 where the vibration spectra of a classical mass spring model on complex networks are investigated; and, in Shang,13 where synchronization of a leader–follower system of coupled harmonic oscillators connected by dampers is studied in the presence of random noises and time delays with an interaction topology modeled by a weighted directed graph.

In general, the presence of dissipative couplings has the effect to equate the state of the nodes so that the system enters a positively invariant set, the absorbing domain, in finite time. It is well-known that if G is a constant matrix whose eigenvalues have only negative real parts, then e t G acts as a uniform contraction; that is, | | e t G | | K e η t for some positive real constants, K and η. This fact guarantees that the solution of the equation is uniformly asymptotically stable, and the synchronization manifold remains invariant under the flow of the equation that governs the system.

For numerical estimates and model prediction purposes, Pietras and Daffertshofer14 stress the importance of reduction of complex oscillatory systems. A number of surveys are devoted to this topic.15 Arola-Fernández et al.16 study the synchronized state in a population of network-coupled heterogeneous oscillators, showing that the steady-state solution of the linearized dynamics may be written as a geometric series whose terms represent the different spatial scales of the network. They prove that this expansion converges for arbitrary frequency distributions and for both undirected and directed networks, provided that the adjacency matrix is primitive. Also, Eroglu et al.17 propose a survey on the main theory behind complete, generalized, and phase synchronization phenomena in complex networks with an approach close to the one followed in the present paper.

Investigations of the transient and steady-state dynamics of synchronous machines have significantly impacted the description of power grid networks. Bhatta et al.18 study the swing equation and its linearization on transmission networks, highlighting the role of symmetries in reducing the computational complexity of modal decomposition.19 Bhatta et al.20 extend this approach to multilayer networks.

A crucial issue is the time and speed of synchronization. Grabow et al.21 focus on the speed of convergence as a collective time scale for synchronizing systems on different topologies, ranging from completely ordered and grid-like to completely disordered and random topologies. They find, for instance, that, at fixed in-degree, topological randomness induces shorter synchronization times, whereas intermediate randomness in the small-world regime induces longer synchronization times.

Many other contributions have been given on specific issues, such as in Skardal et al.22 and Skardal et al.,23 where authors introduce a synchrony alignment function that encodes the interplay between the network structure and oscillator frequencies as an objective measure for the synchronization properties of a network of heterogeneous oscillators with natural frequencies or in Su et al.,24 where authors revisit the synchronization problem for coupled oscillators in a dynamic proximity network. Dörfler et al.25 prove a closed-form condition for synchronization that can be stated in terms of the Moore–Penrose pseudoinverse of the Laplacian matrix. Liu and Barabási,26 in the wider context of control theory, address the problem of how networks organize themselves to balance control with functionality. In particular, they review the process of local pinning synchronization, when controllers, applied to a subset of nodes, help synchronize the network. Also, in Lu et al.,27 pinning synchronization in the new setting of networks of networks is investigated. Slotine et al.28 study the role, in both nature and system design, of co-existing power leaders, to which the networks synchronize, and knowledge leaders, to whose parameters the networks adapt.

One of the most important contributions to the study of synchronization in coupled systems is undoubtedly the Kuramoto model.29 In this model, the underlying topology is fully connected, even if it was later adapted to the scenario of nearest neighbor interactions. Despite the simplicity and effectiveness of the Kuramoto model and all its variants, several works have highlighted some of its shortcomings. For example, Shahal et al.30 address the issue of the importance of synchronization in human networks for our civilization. They study, in an accurate experimental setting, the synchronization between violin players, arranged in a network framework with full control in terms of network connectivity, coupling strength, and delay. Their results show that players can adjust their playing period and cancel connections by ignoring frustrating signals to find a stable solution. Their observations reveal that the usual models for coupled networks, such as the Kuramoto model, cannot always be applied to human networks, as additional degrees of freedom enable new strategies and produce better solutions than are possible within such a model. Corroborating these observations, it is argued that, during social interactions, participants are continuously active, each modifying their actions in response to the changing actions of their partners. This continuous mutual adaptation results in an interactional synchrony to which all members contribute. Dumas et al.31 study the brain activity of the participants in this process and found that states of interactional synchrony correlate with the emergence of an interbrain synchronization network between right centroparietal regions. These regions have been found to play a central role in social interaction within an interindividual brainweb.

Another limiting aspect of classical models is the presence of strictly local interactions. For this reason, Estrada et al.32 study how long-range interactions can affect synchronization dynamics, proposing a model of coupled oscillators based on the d-path Laplacian matrices. Their analysis reveals that long-range interactions improve network synchronizability with an impact that depends on the original structure.

The dynamics of opinion formation within a social network is also a widely studied field. In this kind of network, each node can be assigned, for instance, a discrete or continuous score according to his/her agreement or disagreement on a given topic. By introducing a coupling between members in order to describe the dragging effect in their opinion formation process, each member in the network can change his/her opinion according to a synchronization mechanism or according to some kind of resonant phenomenon with some source/leader status.

Related to this issue, Jenkins33 presents a very interesting analysis of the so-called self-oscillation phenomenon, which is distinct from the more familiar phenomenon of forced resonance, whereas Tönjes et al.34 focus on networks with influencers, a group of hubs that couple strongly and connect different parts of the network. They show that subjecting influencers to an optimal noise can result in enhanced network synchronization. This effect, called coherence resonance, emerges from a synergy between a network structure and stochasticity and is highly nonlinear, vanishing as the noise is too weak or too strong.

In the light of the previous literature review, this paper aims to analyze, in a sense from first principles, some aspects of the synchronization and resonance phenomena and, in general, of the transient dynamics in networks of oscillators with a heterogeneous topological structure. A unified approach is maintained in the description of all phenomena. The main new contributions and results of the present work can be summarized as follows:

  1. a synthetic framework in which to inscribe the description of synchronization and resonance phenomena on networks of coupled identical harmonic oscillators with arbitrary topology;

  2. a closed expression of the average synchronization time that encodes information about the topological structure of the network;

  3. new evidences that, for a weakly coupled system, complete and dense networks achieve the state of perfect synchronization more rapidly than, for instance, the path network and that the complete network shows a synchronization time, which is shorter for weak coupling than for strong coupling;

  4. closed expressions for the resonance frequencies of the network as a function of the eigenvalues of the Laplacian matrix;

  5. a detailed account of how a single node in the network can serve as a leader (influencer) or a follower (receiver) depending on its topological location when a resonance phenomenon arises inside the network;

  6. a proof that the presence of dissipative terms pulls down the possible resonance frequencies to a single one, corresponding to the synchronization frequency;

  7. a new perspective on the linearization of the swing equation and its solution in the proposed setting; and

  8. some new mathematical properties of the polar decomposition of the symplectic matrix governing the synchronization phenomenon.

This paper is organized as follows. Section II introduces the basic notations and the general model. Section III has a didactic purpose: it shortly illustrates the case of networked coupled identical linear oscillators without damping. In Sec. IV, existing results are framed into the broader context described, and new ones about the mathematical properties of the involved operators and the rate of synchronization are presented. Section V is devoted to resonance phenomena and, after some recalls on the case of the single driven oscillator, provides new results related to resonant frequencies on networks when an external driving force is placed at a specific node. Section VI deals with a special case of the general equation introduced, specifically the linear swing equation and its solution. Finally, Sec. VII applies the previous findings to two real world networks, the Zachary club social network and the Syrian power grid, by illustrating a number of interesting remarks. All the proofs of the propositions in this article are given in  Appendixes A and  B.

Our purpose is basically to describe the behavior of a system of n coupled harmonic oscillators within a network structure. To model the local interactions between the n harmonic oscillators, we use an undirected graph G = ( V , E ) without loops, where V is the set of vertices or nodes and E V × V is the set of edges or links. Throughout the paper, we will assume that the network is connected, i.e., that a path exists between each pair of distinct vertices. The actual configuration is completely encoded in the adjacency matrix A = { a i j } , i , j = 1 , , n, where a i j = 1 when ( i , j ) E, that is, when nodes i and j show an actual local interaction, and 0 otherwise. We will assume a i j = a j i, which is a perfectly symmetric interaction within a pair of nodes. The degree of a node i is defined as k i = j = 1 n a i j = ( A 2 ) i i, and K = diag { k i } is the diagonal matrix of the degrees of all the nodes. We will denote by λ i and ψ i , i = 1 , , n, respectively, the eigenvalues and eigenvectors of the adjacency matrix A, where the eigenvalues are listed in decreasing order: λ 1 λ n. The Laplacian matrix is then defined as L = K A. We will denote by μ i and ϕ i , i = 1 , , n, respectively, the eigenvalues and eigenvectors of the Laplacian matrix L, where again, the eigenvalues are listed in decreasing order: μ 1 > μ 2 > > μ n = 0. In particular, ϕ n = 1 n [ 1 , 1 , , 1 ] T = 1 n u T, being u R n the vector of all 1’s. Similarly, we will denote by 1 = u u T the n square matrix of all 1’s and by I the identical matrix.

To have a concrete picture of the system of coupled oscillators, let us look at the network in Fig. 1. This example will be used throughout the paper as a toy model. We have n = 4 objects of mass m i , i = 1 , , 4, attached to a rigid support by identical springs. These objects are also connected to each other according to a network scheme with adjacency matrix A by other identical springs. All vertical springs have the same elastic constant c 1, and all the springs connecting different objects have the same elastic constant c 2. Each mass is allowed to oscillate only in the vertical direction around the equilibrium position. Let x = [ x 1 , , x n ] T be the vector of the displacements of each mass from its equilibrium position so that x i R represents the one-dimensional (vertical, in the figure) displacement of the ith object. Let F i be the force acting on mass i due to the vertical spring and F i j be the force acting on mass i due to its interaction with mass j. Then, F i = c 1 x i and F i j = c 2 ( x i x j ), and the total force on node i is given by
F i tot = c 1 x i c 2 j a i j ( x i x j ) = ( c 1 + c 2 k i ) x i + c 2 j a i j x j .
(1)
FIG. 1.

A simple example of a network of coupled harmonic oscillators.

FIG. 1.

A simple example of a network of coupled harmonic oscillators.

Close modal
The forces acting on all the nodes in the network are then collected in the force vector:
F elastic tot = c 1 x c 2 ( K A ) x = [ c 1 I + c 2 L ] x .
(2)
Let us add now the damping forces, which are proportional to the velocities x ˙ = d x ( t ) / d t with coefficients c 1 for the vertical damping and c 2 for the coupling between nodes. As above, in a similar way, they can be collected into the vector:
F damping tot = [ c 1 I + c 2 L ] x ˙ .
(3)
The most general equation of motion of the network of coupled and damped harmonic oscillators with masses M = diag { m i } is then given by
M x ¨ = [ c 1 I + c 2 L ] x [ c 1 I + c 2 L ] x ˙
(4)
which is equivalent to the system of 2 n differential equations given by
{ x ˙ = v , M v ˙ = [ c 1 I + c 2 L ] x [ c 1 I + c 2 L ] v ,
(5)
with initial conditions x ( 0 ) = x 0 and v ( 0 ) = v 0. Let us define the vector y R 2 n as
y := ( x v )
(6)
so that, in a matrix form,
y ˙ = ( x ˙ v ˙ ) = [ 0 I M 1 ( c 1 I + c 2 L ) M 1 ( c 1 I + c 2 L ) ] ( x v ) .
(7)
We will focus on the case of coupled identical harmonic oscillators so that M = I. Hence, the system (7) becomes
{ y ˙ = G y , y ( 0 ) = y 0 ,
(8)
where the matrix G is
G := [ 0 I ( c 1 I + c 2 L ) ( c 1 I + c 2 L ) ] .
(9)
We stress that coefficients c 2 and c 2 determine the strength of the coupling between the oscillators and that in the following, the coupling is meant to be weak when c 2 ( ) / c 1 ( ) < 1.
Since we are interested in resonant phenomena, we add to the system an external periodic force applied to one node. Let us imagine that such a driving force f ( t ), f : [ 0 , + ) R, acts on node h , h = 1 , , n; node h will be called the leader. The total force on node i is now
F i tot = c 1 x i c 2 j a i j ( x i x j ) + f ( t ) δ h i ,
(10)
where δ h i = 1 only if i = h, 0 otherwise. The forces acting on the nodes in the network are then described by the vector
F tot = c 1 x c 2 ( K A ) x c 1 v c 2 ( K A ) v + f ( t ) e h ,
(11)
where e h = [ 0 , 0 , , 1 , , 0 ] T, with 1 in position h, is in the standard basis of R n. System (8) is then generalized by
{ y ˙ = G y + b ( t ) , y ( 0 ) = y 0 ,
(12)
where
b ( t ) = f ( t ) [ 0 e h ] .
(13)

Finally, in addition to the cases described above, we will study a further variant of the system (12) in order to show how this formalism can accommodate the linearization of an important nonlinear problem. In particular, by setting c 1 = 0, c 2 = 1, c 1 = 1, c 2 = 0, and by giving the factor b ( t ) in Eq. (12), the form b ( t ) = [ 0 , p ] T, where p is a constant vector in R n, we will gain the linear swing equation describing the dynamic response of a network of synchronous machines to small disturbances.

Shortly, we can summarize in Table I the cases we are interested in, where c 1, c 2, c 1 , c 2 are the coupling terms in the matrix G in Eq. (9) and the last column specifies the form of the driving term b ( t ) in Eq. (12).

TABLE I.

The five cases discussed in the paper: coupled harmonic oscillators; coupled and damped harmonic oscillators; coupled and forced harmonic oscillators; coupled, damped, and forced harmonic oscillators; and the linear swing equation.

Network system c1 c2 c 1 c 2 b(t)
Coupled 
Coupled and damped 
Coupled and forced  f(t)[0, eh]T 
Coupled, damped, and forced  f(t)[0, eh]T 
Linear swing equation  [0, p]T 
Network system c1 c2 c 1 c 2 b(t)
Coupled 
Coupled and damped 
Coupled and forced  f(t)[0, eh]T 
Coupled, damped, and forced  f(t)[0, eh]T 
Linear swing equation  [0, p]T 

The first case in Sec. III is for illustrative purposes only, and it is used to properly introduce notations. The second case in Sec. IV will lead to synchronization phenomena on the network. The third and fourth cases in Sec. V are the core of the paper and lead to the resonance phenomena on networks. Finally, the fifth case in Sec. VI will allow us to discuss some properties of the linear swing equation and its solution.

The analysis of this first case is worthwhile since it allows the exact solution of system (8) to be written in a closed form useful for later discussion. Therefore, let us go back to system (8) and set c 1 = 0 and c 2 = 0,
( x ˙ v ˙ ) = [ 0 I ( c 1 I + c 2 L ) 0 ] ( x v ) .
(14)

Let us denote H = c 1 I + c 2 L and B the square root matrix of H: B 2 = H. Since, for any value of the positive coefficients c 1 and c 2, the matrix H is positive definite, B can be easily obtained as B = Φ Λ 1 / 2 Φ T, where Φ is the matrix whose columns are the eigenvectors of H and L and Λ is the diagonal matrix of the positive eigenvalues of H.

It is well known that the general solution of the equation x ¨ + B 2 x = 0, equivalent to system (14), is given by
x ( t ) = e i B t x 0 ( 1 ) + e i B t x 0 ( 2 ) ,
(15)
where x 0 ( 1 ) and x 0 ( 2 ) are complex conjugate constant vectors satisfying the following conditions:
{ x ( 0 ) =: x 0 = x 0 ( 1 ) + x 0 ( 2 ) , x ˙ ( 0 ) =: v 0 = i B ( x 0 ( 1 ) x 0 ( 2 ) ) .
(16)
Therefore,
{ x 0 ( 1 ) = 1 2 [ x 0 i B 1 v 0 ] , x 0 ( 2 ) = 1 2 [ x 0 + i B 1 v 0 ] ,
(17)
and the general solution (15) takes the usual form
x ( t ) = e i B t + e i B t 2 x 0 e i B t e i B t 2 i B 1 v 0 = cos ( B t ) x 0 + sin ( B t ) B 1 v 0 .
(18)
Let us notice that solution (18) of the system (8) can be deduced equivalently by observing that
y ( t ) = e G t y 0 ,
(19)
where now matrix G in Eq. (9) is G = [ 0 I H 0 ] .
A straightforward computation shows that, in this case,
e G t = [ I 0 0 I ] + [ 0 I H 0 ] t + 1 2 ! [ H 0 0 H ] t 2 + 1 3 ! [ 0 H H 2 0 ] t 3 + 1 4 ! [ H 2 0 0 H 2 ] t 4 + = [ cos B t B 1 sin B t B sin B t cos B t ] ,
(20)
where, as above, B 2 = H. Hence,
y ( t ) = [ cos B t B 1 sin B t B sin B t cos B t ] ( x 0 v 0 ) ,
(21)
which is equivalent to
{ x ( t ) = cos ( B t ) x 0 + B 1 sin ( B t ) v 0 , v ( t ) = B sin ( B t ) x 0 + cos ( B t ) v 0 .
(22)

Figure 2 illustrates the harmonic motion of the four nodes of the network in Fig. 1 given by solution x i ( t ) in Eq. (22). The values for the two coefficients are c 1 = c 2 = 1, and the initial conditions are x 0 = [ 1 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T. On the horizontal axis, time ranges in [ 0 , 2 π ].

FIG. 2.

Motion of the harmonic oscillators in the example network. Blue square dots refer to node 1, red circle dots to node 2, yellow triangular dots to node 3, and green kite dots to node 4.

FIG. 2.

Motion of the harmonic oscillators in the example network. Blue square dots refer to node 1, red circle dots to node 2, yellow triangular dots to node 3, and green kite dots to node 4.

Close modal

Let us remark that, if we set c 1 = c 1 = c 2 = 0 and c 2 = 1 in system (8), we have a network of free harmonic oscillators connected by identical springs.

In this case, H = L, and the system is described by x ¨ + L x = 0. Since det B = det L = 0, B 1 does not exist, and if rank ( B ) = rank ( B | v 0 ) = n 1, system (16) admits 1 solutions. In fact, if we set x 0 ( 1 ) = a i b and x 0 ( 2 ) = a + i b, system (16) is equivalent to
{ x 0 = 2 a , v 0 = 2 B b .
(23)

Therefore, there are 1 values for the imaginary part b of the initial coefficients. Let us notice that B 1 does not exist because of the null eigenvalue μ n = 0 of the Laplacian L and then of the matrix B. This eigenvalue is associated with the global translation of all the nodes with the same velocity, i.e., to the stationary state of the whole network, which, in this case, is free and not bounded to any fixed support. Let us also observe that since j B i j = i B i j = 0, i.e., the sum of the elements along a row or a column in the matrix B is null, the initial velocity vector v 0 has to satisfy i v 0 i = 2 i j B i j b j = 2 j i B i j b j = 0. Then, for the system (16) to be consistent, the initial velocity vector must satisfy the condition i v 0 i = 0, equivalent to the classical linear momentum conservation.

Let us set in system (8): c 1 = c 2 = 1 and c 2 = c 1 = 0. This choice for the values of the parameters corresponds to a system of harmonic oscillators attached to a fixed support, free to oscillate around their rest point with constant c 1 = 1 but coupled in a damping network with coupling coefficient c 2 = 1. Basically, it represents a system of n networked point-mass agents connected and interacting only by virtual dampers and can be viewed as a general distributed algorithm for the synchronization of their positions and velocities. Matrix G in Eq. (9) now becomes
G = [ 0 I I L ] .
(24)
Let us first observe that G is a symplectic matrix, that is, G T J G = J, where J is the 2 n-square antisymmetric matrix (symplectic identity),
J = [ 0 I I 0 ] .
(25)

As a consequence, we have that if λ is an eigenvalue of G, so are its complex conjugate λ , 1 λ, and, hence, 1 λ . Moreover, if λ has multiplicity m, then 1 λ has multiplicity m. Finally, both G and G 1 = J G T J have the same set of eigenvalues. In particular, denoting by λ i G ± , i = 1 , , n, the eigenvalues of G and ψ i G ± the corresponding eigenvectors, Ren11 proves the following proposition:

Proposition 1
The eigenvalues of G are the n couples of reciprocal values given by
λ i G ± = 1 2 [ μ i ± μ i 2 4 ] ,
(26)
and the corresponding normalized eigenvectors are given by
ψ i G ± = 1 1 + | λ i G ± | 2 [ ϕ i λ i G ± ϕ i ] ,
(27)
where μ i and ϕ i are the eigenvalues and the corresponding eigenvectors of the Laplacian matrix L. In  Appendix B, we propose a straightforward alternative proof for the undirected case in line with the one proposed by Ren11 for the directed case. Here, we focus on the following remarks.

Remark

Proposition 1 has the following straightforward consequence: G is a Hurwitz matrix; i.e., all its eigenvalues are negative or have a negative real part: Re ( λ i G ± ) < 0 , i = 1 , , n. Indeed, we may have the following four possibilities as in Table II.

TABLE II.

Eigenvalues of the matrix G.

If μi > 2  λ i G R λ i G ± < 0 
If μi = 2  λ i G R λ i G ± = 1 
If 0 < μi < 2  λ i G C 1 < Re ( λ i G ± ) = μ i 2 < 0 
For μn = 0  λ n G I λ n G ± = ± i 
If μi > 2  λ i G R λ i G ± < 0 
If μi = 2  λ i G R λ i G ± = 1 
If 0 < μi < 2  λ i G C 1 < Re ( λ i G ± ) = μ i 2 < 0 
For μn = 0  λ n G I λ n G ± = ± i 

Let us also notice that, since μ 1 > k max 2 for any nontrivial network, there is always a real negative eigenvalue. Moreover, as expected, the trace of the matrix G is equal to
tr G = j = 1 n ( λ j G + + λ G - j ) = j = 1 n μ j = j = 1 n k j = k tot ,
where k tot is the total degree of the network.
A general solution to system (8) is again of the usual form
y ( t ) = e G t y 0 ,
(28)
with matrix G as in Eq. (24).
Remark
Let us immediately notice that if we apply the general time evolution in Eq. (28) to an initial state equal to an eigenvector ψ i G ± as in Eq. (27), such that μ i > 0, we get that, for t + ,
y ( t ) = e G t ψ i G ± = e λ i G ± t ψ i G ± 0
since Re ( λ i G ± ) < 0. We can then identify the eigenstates of the matrix G as the initial states such that the synchronization leads to the extinction of the oscillation for all the nodes of the network.
In  Appendix A, we propose some theoretical results about the polar decomposition of the matrix G in Eq. (24), which have relevance in the synchronization process. In particular, we will show that the asymptotic behavior is entirely encoded in the unitary component U of the matrix G. Since G is a real non-singular matrix, its polar decomposition can be directly obtained as G = G ( G T G ) 1 / 2 ( G T G ) 1 / 2. Therefore, we set
G = U P ,
(29)
where
U = G ( G T G ) 1 / 2 and P = ( G T G ) 1 / 2 .
(30)

The next proposition states that, in terms of asymptotic behaviors, the evolution operators e G t and e U t are equivalent.

Proposition 2
For any initial condition y ( 0 ) = y 0, both e G t y 0 and e U t y 0 show the same asymptotic behavior y ~ ( t ), for t + , where
y ~ ( t ) := 1 n [ cos t 1 sin t 1 sin t 1 cos t 1 ] ( x 0 v 0 ) = 1 n ( cos t 1 x 0 + sin t 1 v 0 sin t 1 x 0 + cos t 1 v 0 )
(31)
and 1 is the n-square all 1’s matrix.

Let us consider again the toy model in Fig. 1. The four eigenvalues of L are μ 1 = 4 , μ 2 = 3 , μ 3 = 1 , μ 4 = 0. The eight eigenvalues of operators G , P , U and the four angles defined in Eq. (A7) in  Appendix A are listed in Table III.

TABLE III.

Eigenvalues of the matrices G, P, and U and angles θ i = arg ( λ i U + ) = 2 arctan λ i P + for the toy model.

λ G P U Angles
λ 1 +  −0.268  4.236  −0.894 + 0.447i  θ1 = 153° 
λ 2 +  −0.382  3.303  −0.832 + 0.555i  θ2 = 146° 
λ 3 +  −0.500 + 0.866i  1.618  −0.447 + 0.894i  θ3 = 117° 
λ 4 +  +i  +i  θ4 = 90° 
λ 4   i  i  … 
λ 3   −0.500 − 0.866i  0.618  −0.447 − 0.894i  … 
λ 2   −2.618  0.303  −0.832 − 0.555i  … 
λ 1   −3.732  0.236  −0.894 − 0.447i  … 
λ G P U Angles
λ 1 +  −0.268  4.236  −0.894 + 0.447i  θ1 = 153° 
λ 2 +  −0.382  3.303  −0.832 + 0.555i  θ2 = 146° 
λ 3 +  −0.500 + 0.866i  1.618  −0.447 + 0.894i  θ3 = 117° 
λ 4 +  +i  +i  θ4 = 90° 
λ 4   i  i  … 
λ 3   −0.500 − 0.866i  0.618  −0.447 − 0.894i  … 
λ 2   −2.618  0.303  −0.832 − 0.555i  … 
λ 1   −3.732  0.236  −0.894 − 0.447i  … 

In Figures 3–6, we illustrate some different synchronization processes on the same network, with different initial conditions, as specified in the caption of the figures. For each initial condition, plots represent the position and velocity of the four nodes and the phase portrait with the corresponding limiting cycle.

FIG. 3.

Synchronization of the harmonic oscillators in the example network with initial conditions x 0 = [ 1 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T: panel (a) position, panel (b) velocity, and panel (c) phase portrait.

FIG. 3.

Synchronization of the harmonic oscillators in the example network with initial conditions x 0 = [ 1 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T: panel (a) position, panel (b) velocity, and panel (c) phase portrait.

Close modal
FIG. 4.

Synchronization of the harmonic oscillators in the example network with initial conditions x 0 = [ 1 , 0.5 , 1 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T: panel (a) position, panel (b) velocity, and panel (c) phase portrait.

FIG. 4.

Synchronization of the harmonic oscillators in the example network with initial conditions x 0 = [ 1 , 0.5 , 1 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T: panel (a) position, panel (b) velocity, and panel (c) phase portrait.

Close modal
FIG. 5.

Synchronization of the harmonic oscillators in the example network with initial conditions x 0 = [ 1 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0.4 , 0.2 , 0 ] T: panel (a) position, panel (b) velocity, and panel (c) phase portrait.

FIG. 5.

Synchronization of the harmonic oscillators in the example network with initial conditions x 0 = [ 1 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0.4 , 0.2 , 0 ] T: panel (a) position, panel (b) velocity, and panel (c) phase portrait.

Close modal
FIG. 6.

Synchronization of the harmonic oscillators in the example network with initial conditions x 0 = [ 1 , 0.4 , 1 , 0.2 ] T and v 0 = [ 0 , 0.4 , 0.3 , 0 ] T: panel (a) position, panel (b) velocity, and panel (c) phase portrait.

FIG. 6.

Synchronization of the harmonic oscillators in the example network with initial conditions x 0 = [ 1 , 0.4 , 1 , 0.2 ] T and v 0 = [ 0 , 0.4 , 0.3 , 0 ] T: panel (a) position, panel (b) velocity, and panel (c) phase portrait.

Close modal

The preceding analysis gives us the opportunity to add some insights into the widely discussed topic of synchronization time and speed. What follows contributes to this topic by highlighting the role of network topology in the synchronization process in order to compare the behavior of different systems and to discriminate which systems converge more or less rapidly to the synchronized state.

In the literature, a measure of synchronization is defined as C ( t ) = 1 n i = 1 n Θ ( ε d i ( t ) ), where d i ( t ) = | | y i ( t ) y ~ i ( t ) | |, y ~ ( t ) is the synchronization solution in Eq. (31), ε is a threshold, and Θ is the step function (see, for instance, Grabow et al.21). This index provides the percentage of nodes whose distance from the synchronization solution is less than a given threshold. As the network becomes completely synchronized on this solution, C ( t ) 1.

There is actually no univocal definition of an index that measures the speed of synchronization. For example, it could also be taken as an indicator of the time at which the maximum of the absolute difference between couples of the oscillator positions, that is, max i j | x i ( t ) x j ( t ) |, or, alternatively, the mean of all the differences between couples of positions, drops below a given threshold.

In light of the closed expression of the asymptotic solution in Eq. (31), we will refer directly to the mean absolute difference between positions of nodes and such a limiting solution.

We aim here at defining a proper synchronization time as a function of the topological properties of the underlying network. Since, in the singular value decomposition, the more negative are the eigenvalues, the faster the corresponding exponentials go to zero, we have to select exponentials with less negative exponents if we are to identify the dominant behavior that slows synchronization and to set a bound on the synchronizability of the global process. Let us consider an initial state given by y 0 = ( x 0 , v 0 ) T, and let us write the expansion of the solution to problem (8) on the basis of the eigenvectors of matrix G,
y ( t ) = e G t y 0 = e t λ n G + ( ψ n G + y 0 ) ψ n G + + e t λ n G ( ψ n G y 0 ) ψ n G + e t λ n 1 G + ( ψ n 1 G + y 0 ) ψ n 1 G + + e t λ n 1 G ( ψ n 1 G y 0 ) ψ n 1 G + = e i t 1 2 n ( u i u ) ( x 0 v 0 ) 1 2 n ( u i u ) + e i t 1 2 n ( u i u ) ( x 0 v 0 ) 1 2 n ( u i u ) + e t λ n 1 G + ( ψ n 1 G + y 0 ) ψ n 1 G + + e t λ n 1 G ( ψ n 1 G y 0 ) ψ n 1 G + = 1 2 n e i t [ u T x 0 i u T v 0 ] ( u i u ) + 1 2 n e i t [ u T x 0 + i u T v 0 ] ( u i u ) + e t λ n 1 G + ( ψ n 1 G + y 0 ) ψ n 1 G + + e t λ n 1 G ( ψ n 1 G y 0 ) ψ n 1 G + = 1 n [ ( + cos t u T x 0 + sin t u T v 0 ) u ( sin t u T x 0 + cos t u T v 0 ) u ] + e t λ n 1 G + ( ψ n 1 G + y 0 ) ψ n 1 G + + e t λ n 1 G ( ψ n 1 G y 0 ) ψ n 1 G + .
The first term in the last step is equal to y ~ ( t ) in Eq. (31), and it is the synchronization solution so that
y ( t ) y ~ ( t ) = e t λ n 1 G + ( ψ n 1 G + y 0 ) ψ n 1 G + + e t λ n 1 G ( ψ n 1 G y 0 ) ψ n 1 G + .
(32)
All the 2 n 2 terms on the right hand side contain an exponential whose coefficient is λ i G ± , i = n 1 , , 2 , 1. These coefficients can be complex or real, but in both cases, their real part is always negative, according to Table II. Therefore the process is governed by the term that contains the least negative eigenvalue real part. The eigenvalues of G are expressed by Eq. (26), λ i G ± = 1 2 [ μ i ± μ i 2 4 ], as an increasing function of μ i. Since μ 1 > k max 2 for any nontrivial network, there is always a real negative eigenvalue and the maximum real value is attained for i = 1 so that the less negative real eigenvalue is λ 1 G + = 1 2 [ μ 1 + μ 1 2 4 ]. On the other hand, when λ i G ± C, the less negative real part is gained for the minimum non zero μ i, i.e., for i = n 1, corresponding to the algebraic connectivity, and it is given by Re ( λ n 1 G ± ) = μ n 1 2. Therefore, the less negative exponential coefficient is given by
λ S = max ( 1 2 [ μ 1 + μ 1 2 4 ] ; μ n 1 2 ) .
(33)
If we look at Eq. (32) by components
y i ( t ) y ~ i ( t ) = e t λ n 1 G + ( ψ n 1 G + y 0 ) ψ n 1 G + ( i ) + e t λ n 1 G ( ψ n 1 G y 0 ) ψ n 1 G ( i ) + ,
(34)
then we can say that, as t + ,
| y i ( t ) y ~ i ( t ) | e t λ S ( ψ S y 0 ) ψ S ( i ) ,
(35)
where λ S is the eigenvalue in Eq. (33) and ψ S ( i ) is the component i of the corresponding eigenvector (as far as this component is not zero). Now, we want to find out the time after which, for any node i, the difference between the actual solution y i ( t ) and the synchronization solution y ~ i ( t ) becomes less than a given threshold ε > 0: e t λ S ( ψ S y 0 ) ψ S ( i ) < ε. By solving for t, we get
t i > 1 | λ S | log | ( ψ S y 0 ) ψ S ( i ) ε | .
(36)
Finally, the average synchronization time can be bounded by
t mean > 1 n | λ S | i = 1 n log | ( ψ S y 0 ) ψ S ( i ) ε | .
(37)
Remark
When | λ S | = μ n 1 2, the greater the algebraic connectivity, the smaller is the synchronization time. Moreover, the algebraic connectivity is bounded by the inequality μ n 1 n n 1 k min with k min minimum degree in the network, and since the density of the network is given by δ = k mean n 1, we have that
μ n 1 n n 1 k min n n 1 k mean = n δ
(38)
so that, in this case, we have
t > 2 n | μ n 1 | i = 1 n log | ( ψ S y 0 ) ψ S ( i ) ε | 2 n 2 δ i = 1 n log | ( ψ S y 0 ) ψ S ( i ) ε | .
(39)

For n fixed then, as δ grows, then the synchronization time decreases.

Let us test the synchronization time in Eq. (37) on the toy networks with n = 4 and n = 5 nodes in Fig. 7.

FIG. 7.

Small graphs with different synchronization times analyzed in the text.

FIG. 7.

Small graphs with different synchronization times analyzed in the text.

Close modal

For the networks A, B, and C with n = 4, the values of the parameter λ S are λ S A = 0.292 893 2, λ S B = 0.267 949 2 ( 0.381 966 0 ), and λ S C = 0.267 949 2 ( 0.267 949 2 ).35 Then, we expect that network A would synchronize faster than B, which in turn would synchronize faster than C.

In Fig. 8, panel (a), we plot the mean of the absolute values of differences between node position at time t and the synchronization solution, as in Eq. (32), as a function of time t: we observe that the red line (network A) decays on average more rapidly than the blue line (network B), which decays more rapidly than the green line (network C). If we compute the mean synchronization times given in formula (37) with ε = 10 3, we get t mean A = 8.421, t mean B = 17.273, and t mean C = 18.264. The empirical values, computed as the mean of the minimum times t such that the absolute value of the difference between the solution for a given node and the asymptotic solution goes below ϵ, are t emp A = 8.275, t emp B = 17.525, and t emp C = 21.925.

FIG. 8.

Different speed of synchronization for the two sets of toy networks: panel (a) A–C and panel (b) D–F.

FIG. 8.

Different speed of synchronization for the two sets of toy networks: panel (a) A–C and panel (b) D–F.

Close modal

Similarly, for networks D–F with n = 5, the values of the parameter λ S are λ S A = 0.190 983 0, λ S B = 0.208 712 2 ( 0.267 949 2 ), and λ S C = 0.208 712 2 ( 0.208 712 2 ). Then, we expect that network D would synchronize faster than E, which in turn would synchronize faster than F. In Fig. 8, panel (b), we plot again the mean deviation of the node positions at time t and the synchronization solution, as in Eq. (32), as a function of time t: we observe that the red line (network D) decays on average more rapidly than the blue line (network E), which decays more rapidly than the green line (network F). Again, if we compute the mean synchronization times in formula (37) with ε = 10 3, we get t mean D = 11.849, t mean E = 19.868, and t mean F = 26.510, whereas the empirical values are t emp D = 11.140, t emp E = 21.760, and t emp F = 26.940.

This observation could be taken as a hint that a path will synchronize faster than a more integrated or even complete network. Things are actually more elaborate as argued in the following important remark.

Remark

So far, we have assumed a constant value equal to 1 for the coupling coefficient c 2 between nodes. In this remark, we want to show that different networks may exhibit different rankings in the speed and time of synchronization according to the values of the coupling coefficient. Let us consider again the triplet of networks D, E, and F in Fig. 7, and let us assume as synchronization time the lower bound in formula (37) with a threshold equal to ε = 10 6.

When coefficients in system (7) are c 1 = c 2 = 1 and c 2 = c 1 = 0, formula (37) yields
t mean D = 48.02 t mean E = 52.97 t mean F = 59.61 ,
(40)
and the empirical values are
t emp D = 38.10 t emp E = 49.74 t emp F = 60.04.
(41)
These values show again that the path network gains synchronization more rapidly than the complete network. However, if we lower the coupling coefficient to c 2 = 0.1, the same times become
t mean A = 278.57 t mean B = 82.20 t mean C = 43.07
(42)
with empirical values
t emp A = 294.50 t emp B = 105.70 t emp C = 38.24.
(43)
The different behavior is illustrated in Fig. 9, which is the analogous of Fig. 8, panel (b), for the weak coupled networks D, E, and F. This example shows that, for a weakly coupled system, the complete network achieves the state of perfect synchronization far more rapidly than the path network. This is an outcome opposite to the previous one. More surprisingly, the complete network shows a synchronization time, which is shorter for weak coupling than for strong coupling. This fact makes it interesting to analyze how this time can depend on the coupling coefficients through the eigenvalues of the involved matrices.
FIG. 9.

Synchronization for weak coupled networks.

FIG. 9.

Synchronization for weak coupled networks.

Close modal
If we set c 2 = α, with α > 0, the eigenvalues of G become
λ i G ± = 1 2 [ α μ i ± α 2 μ i 2 4 ] .
(44)

As α decreases, more eigenvalues become complex, and if α < 2 μ 1, they are all complex. Therefore, what matters is the real part, which is just α μ i 2. The maximum non-null value is then for all α μ n 1 2. For the path graph, μ n 1 = 2 [ 1 cos ( π / n ) ], and for the complete graph, μ n 1 = n. Then, we have to compare α ( 1 cos ( π / n ) ) and α n 2. Since n > 2 ( 1 cos ( π / n ) ) for every n > 2, we have that α n 2 < α ( 1 cos ( π / n ) ), which shows that, for α < 2 μ 1, the high connectivity of the complete network enables it to synchronize faster than other topologies. Similar results for networks A–C in Fig. 7 confirm these conclusions.

We move now to the third case in Table I by introducing a driving periodic force acting on or even produced by a single node. We will analyze the effect of a standard periodic sinusoidal force with frequency ω. Our purpose is to illustrate some new results about resonant phenomena on networks. In order to introduce method and notations, we quickly revise some ideas about the single driven resonant oscillator.

Let us consider a harmonic oscillator subjected to an external periodic force F ( t ): m x ¨ = F ( t ) k x. Let us set y as in Eq. (6). Then, y ˙ = A y + b, where A = ( 0 1 ω 0 2 0 ), with ω 0 = k m and b = F ( t ) m e := f ( t ) e, with e = [ 0 , 1 ] T. The general solution, with the initial condition y 0 = [ x 0 , v 0 ] T, is
y ( t ) = e A t [ 0 t e A u b ( u ) d u + c ] ,
(45)
where
e A t = [ cos ω 0 t 1 ω 0 sin ω 0 t ω 0 sin ω 0 t cos ω 0 t ] .
(46)
Let us consider a driving force F ( t ) = F 0 sin ( ω t ) so that b ( t ) = f ( t ) e with f ( t ) = F 0 m sin ( ω t ). From Eq. (46), we have
e A u b ( u ) = F 0 m [ 1 ω 0 sin ω 0 u sin ω u cos ω 0 u sin ω u ] ,
(47)
and, by a straightforward computation of the integral in Eq. (45), we get
y ( t ) = F 0 m ω ω 0 2 ω 2 [ sin ω t ω sin ω 0 t ω 0 cos ω t cos ω 0 t ] .
(48)
Let us note that, for ω ω 0, the asymptotic behaviors
x ( t ) F 0 t 2 m ω 0 cos ω 0 t v ( t ) F 0 t 2 m sin ω 0 t
(49)
show that the amplitude of both position and velocity grows linearly in time.
Incidentally, it might be useful to mention that in the case of a periodical delta-type force F ( t ) = F 0 k = 1 n δ ( t k T ), where δ is now a Dirac delta function and T = t n = 2 π ω, the previous solution turns out to be very different. Indeed, from Eq. (46), we have
e A u b ( u ) = F 0 m k = 1 n δ ( u 2 k π ω ) [ 1 ω 0 sin ω 0 u cos ω 0 u ] ,
(50)
and by computing the integral in Eq. (45), we get
0 t e A u b ( u ) d u = F 0 m [ 1 ω 0 k = 1 n sin ( 2 π ω 0 ω k ) k = 1 n cos ( 2 π ω 0 ω k ) ] .
(51)
By the Dirichlet kernels and replacing 2 n + 1 = π + ω t π, we finally obtain
y ( t ) = F 0 m ω 0 sin ω 0 2 t sin ω 0 ω π [ sin ( ω 0 2 t ω 0 ω π ) ω 0 cos ( ω 0 2 t ω 0 ω π ) ] .
(52)

Let us observe that the limit for ω ω 0 in Eq. (52) shows different behavior than in Eq. (48): when ω approaches ω 0, only the term sin ω 0 ω π in the denominator approaches 0, while other terms represent an oscillation of frequency ω 0. Therefore, the resonance produced by this kind of force preserves the proper frequency ω 0, while the amplitude goes to infinity uniformly for all t.

Let us now move to the network case. The general solution of system (12) is given by
y ( t ) = e G t [ 0 t e G u b ( u ) d u + c ] .
(53)
In the absence of dissipative terms, by setting c 1 = c 2 = 1, matrix G in Eq. (9) becomes
G = [ 0 I H 0 ] ,
(54)
where now, H = I + L. In Eq. (20), we proved that
e G t = [ cos H t ( H ) 1 sin H t H sin H t cos H t ] .
(55)

Let us assume a driving force f ( t ) = F 0 sin ( ω t ) applied on node h, h = 1 , , n. The general solution and the corresponding resonant states are described by the next proposition.

Proposition 3
The solution in Eq. (53) to system (12), with y ( 0 ) = y 0 = 0, for a sinusoidal driving force f ( t ) = F 0 sin ω t acting on node h, h = 1 , , n, is
y ( t ) = F 0 i = 1 n ω ω i 2 ω 2 ϕ i ( h ) [ ( sin ω t ω sin ω i t ω i ) ϕ i ( cos ω t cos ω i t ) ϕ i ] ,
(56)
where ω i := 1 + μ i and ϕ i are the eigenvectors of the Laplacian matrix.
The reader can find the proof of proposition (3) in  Appendix B. Here, we want to focus on its consequences. Equation (56) can be expressed by components as
x k ( t ) = F 0 ω i = 1 n ϕ i ( h ) ϕ i ( k ) ω 2 ω i 2 ( sin ω i t ω i sin ω t ω ) ,
(57)
v k ( t ) = F 0 ω i = 1 n ϕ i ( h ) ϕ i ( k ) ω 2 ω i 2 ( cos ω i t cos ω t ) .
(58)
In particular,
x k ( t ) = F 0 ω i = 1 n 1 ϕ i ( h ) ϕ i ( k ) ω 2 ω i 2 ( sin ω i t ω i sin ω t ω ) + F 0 n ω ω 2 1 ( sin t sin ω t ω )
(59)
since μ 1 > μ 2 > > μ n = 0, ω 1 > ω 2 > > ω n = 1, and ϕ n ( k ) = 1 n, k.
Let us analyze first what happens when ω ω ı ~ for some specific ı ~ = 1 , , n, that is, when the frequency of the external force approaches one of the proper frequencies of the network. The diverging term in the position components in Eq. (57) is
x k ( t ) F 0 ω ϕ ı ~ ( h ) ϕ ı ~ ( k ) ω 2 ω ı ~ 2 ( sin ω ı ~ t ω ı ~ sin ω t ω ) = F 0 ϕ ı ~ ( h ) ϕ ı ~ ( k ) ω ı ~ ( ω + ω ı ~ ) ( ω ω ı ~ ) [ ω sin ω ı ~ t ω ı ~ sin ω t ] .
(60)
Since, in general, x sin t x 0 x 0 sin t x x x 0 sin t x 0 t x 0 cos t x 0 for x x 0, we have
x k ( t ) F 0 ϕ ı ~ ( h ) ϕ ı ~ ( k ) ω ı ~ ( ω + ω ı ~ ) [ sin ω ı ~ t ω ı ~ t cos ω ı ~ t ] = F 0 ϕ ı ~ ( h ) ϕ ı ~ ( k ) 2 ω ı ~ 2 [ sin ω ı ~ t ω ı ~ t cos ω ı ~ t ] .
(61)
In a similar vein, the diverging term in the velocity components in Eq. (58), for ω ω ı ~, is
v k ( t ) F 0 ω ϕ ı ~ ( h ) ϕ ı ~ ( k ) ω 2 ω ı ~ 2 ( cos ω ı ~ t cos ω t ) = F 0 ω ϕ ı ~ ( h ) ϕ ı ~ ( k ) sin ω + ω ı ~ 2 t ω + ω ı ~ sin ω ω ı ~ 2 t ω ω ı ~ 2 F 0 2 ϕ ı ~ ( h ) ϕ ı ~ ( k ) t sin ω ı ~ t .
(62)

The first conclusion is that the amplitudes of both x k ( t ) and v k ( t ), in general, grow linearly with t as ω approaches one of the proper frequencies of the network. However, things are slightly more involved than this, as discussed in the next remark.

Remark

Let us imagine that, for i = ı ~, the product ϕ ı ~ ( h ) ϕ ı ~ ( k ) is equal to 0.

Let us begin by examining the case in which ϕ ı ~ ( k ) = 0. Then, even if ω ω ı ~, the node k cannot resonate with the source node h at such a frequency ω ı ~. In other terms, if we make the source node h oscillate with frequency ω = ω ı ~ and ϕ ı ~ ( k ) = 0, then node k is not affected by the oscillation in h. It is as if node k is transparent to the propagation of the oscillation of frequency ω ı ~ from the source node h. The network as a whole cannot transmit such a frequency to node k, and this node does not participate to the global resonant process.

On the other hand, if it is ϕ ı ~ ( h ) = 0, then an oscillation of frequency ω = ω ı ~ in the source node h cannot be spread throughout the network to any node, and no resonance phenomenon can be induced from such a node with that frequency.

The only frequency that is able to always put in resonance the entire network is ω ı ~ = ω n = 1. In fact, if we look at Eq. (59), being ϕ n ( h ) 0 and ϕ n ( k ) 0, h , k, when ω 1, any receiver node k resonates with any source node h oscillating at such a ground frequency.

More, specifically, let us observe that, for ω 1, we have
x k ( t ) F 0 i = 1 n 1 ϕ i ( h ) ϕ i ( k ) μ i ( sin t sin ω i t ω i ) + F 0 2 n ( sin t t cos t )
(63)
since 1 ω i 2 = μ i. The quantity L h k + := i = 1 n 1 ϕ i ( h ) ϕ i ( k ) μ i represents the h k-component of the Moore–Penrose pseudo-inverse of the Laplacian operator. Sometimes, it is also called vibrational communicability and denoted by G h k ( v ) := L h k +. It can be interpreted as the correlation function between the displacements x h ( t ) and x k ( t ). Formula (63) contains a term that is exactly the product of the vibrational communicability between source and receiver nodes and the driving force F 0 sin t. This term tunes the extent of the interaction between the source and receiver node when ω tends to the ground resonance frequency. More generally, the ratio
ϕ i ( h ) ϕ i ( k ) ω 2 ω i 2
(64)
in Eqs. (57) and (58) contains the information about the local communicability between nodes, and it conveys the role of the network topology in the resonant process from the source node h to the receiver node k at the ground frequency.
To exemplify, let us return to the initial toy network in Fig. 1. Let us take a look at the explicit matrix of eigenvectors,
M = [ 0.289 0.707 0.408 0.500 0.289 0.707 0.408 0.500 0.866 0 0 0.500 0.289 0 0.816 0.500 ] ,
(65)
with eigenvalues μ 1 = 4, μ 2 = 3, μ 3 = 1, μ 4 = 0. The eigen-frequencies are then ω 1 = 5, ω 2 = 2, ω 3 = 2, ω 4 = 1. Let us suppose that the driving force is applied in the first node; that is, h = 1. Then, the motion of the four nodes is described by the following functions:
x 1 ( t ) = F 0 ω [ 0.084 ω 2 5 ( sin 5 t 5 sin ω t ω ) + 0.500 ω 2 4 × ( sin 2 t 2 sin ω t ω ) + 0.166 ω 2 2 ( sin 2 t 2 sin ω t ω ) + 0.250 ω 2 1 ( sin t sin ω t ω ) ( sin 5 t 5 sin ω t ω ) ] ,
x 2 ( t ) = F 0 ω [ 0.084 ω 2 5 ( sin 5 t 5 sin ω t ω ) 0.500 ω 2 4 ( sin 2 t 2 sin ω t ω ) + 0.166 ω 2 2 ( sin 2 t 2 sin ω t ω ) + 0.250 ω 2 1 ( sin t sin ω t ω ) ( sin 5 t 5 sin ω t ω ) ] ,
x 3 ( t ) = F 0 ω [ 0.250 ω 2 5 ( sin 5 t 5 sin ω t ω ) + 0.250 ω 2 1 ( sin t sin ω t ω ) ( sin 5 t 5 sin ω t ω ) ] ,
x 4 ( t ) = F 0 ω [ 0.084 ω 2 5 ( sin 5 t 5 sin ω t ω ) 0.333 ω 2 2 ( sin 2 t 2 sin ω t ω ) + 0.250 ω 2 1 ( sin t sin ω t ω ) ( sin 5 t 5 sin ω t ω ) ] .

The four possible configurations of the phases of nodes, for the different resonant driving frequencies applied in node 1, are depicted in Fig. 10.

FIG. 10.

Phase configuration of the whole network for different resonant driving frequencies: (a) ω 1 = 5, (b) ω 2 = 2, (c) ω 3 = 2, and (d) ω 4 = 1. Red nodes have phase +, green nodes phase 0, and blue nodes phase .

FIG. 10.

Phase configuration of the whole network for different resonant driving frequencies: (a) ω 1 = 5, (b) ω 2 = 2, (c) ω 3 = 2, and (d) ω 4 = 1. Red nodes have phase +, green nodes phase 0, and blue nodes phase .

Close modal

In particular, if the driving force in node number 1 has frequency ω 2 = 2, node numbers 3 and 4 do not resonate since ϕ 2 ( 3 ) = ϕ 2 ( 4 ) = 0; similarly, when node number 1 oscillates with frequency ω 3 = 2, only node number 3 is prevented from entering the resonant state.

We also said that if the driving force is placed in a node h such that ϕ i ( h ) = 0, then if ω ω i, we have x k ( t ) = 0 , k = 1 , , n. For instance, if the force is placed in node 3, being ϕ 2 ( 3 ) = 0, this implies that for ω ω 2 = 2, all the nodes stay at rest. The network is not able to transmit such a frequency from such a node. Other frequencies may have an effect, but there will be no resonance with frequency ω 2 = 2 placed in node number 3.

In Fig. 11, we illustrate the position plots of the four nodes in the toy network when an external force f ( t ) = sin ( ω t ) is applied to node number 1 starting from initial conditions x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T. As ω approaches ω 4 = 1 from below, all the nodes are activated and enter a global resonant state, as shown in Fig. 9 panel (d). In Fig. 12, we illustrate the position plots of the four nodes under the same conditions as before but as ω approaches ω 3 = 2 from below; let us notice that all nodes are activated, except node number 3 (green line in the plot), which does not respond to node 1 stimuli at this frequency, as shown in Fig. 9 panel (c).

FIG. 11.

Approaching the resonant state with a driving force in node number 1 of the toy network and x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T; the frequencies are panel (a) ω = 0.80, panel (b) ω = 0.95, and panel (c) ω = 0.99.

FIG. 11.

Approaching the resonant state with a driving force in node number 1 of the toy network and x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T; the frequencies are panel (a) ω = 0.80, panel (b) ω = 0.95, and panel (c) ω = 0.99.

Close modal
FIG. 12.

Approaching the resonant state with a driving force in node number 1 of the toy network and x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T; the frequencies are panel (a) ω = 1.20, panel (b) ω = 1.35, and panel (c) ω = 1.40.

FIG. 12.

Approaching the resonant state with a driving force in node number 1 of the toy network and x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T; the frequencies are panel (a) ω = 1.20, panel (b) ω = 1.35, and panel (c) ω = 1.40.

Close modal
We conclude with some remarks on the case in which the coupling between nodes is due to dissipative terms in the presence of a periodic external force of the same type discussed in Subsection V B. In this case, the fourth in Table I, matrix G in Eq. (9) takes the form
G = [ 0 I I L ] .
(66)
A straightforward computation shows that now,
e G t = [ 1 2 B 1 [ ( B + L ) e B L 2 t + ( B L ) e B + L 2 t ] B 1 [ e B L 2 t e B + L 2 t ] B 1 [ e B + L 2 t e B L 2 t ] 1 2 B 1 [ ( B L ) e B L 2 t + ( B + L ) e B + L 2 t ] ] ,
(67)

where B = L 2 4 I.36 Let us now suppose again a driving force acting on node h of the form f ( t ) = F 0 sin ω t. The general solution and the corresponding resonant states are described by the next proposition.

Proposition 4
The solution y ( t ) = [ x ( t ) , v ( t ) ] T = e G t [ 0 t e G u b ( u ) d u + c ] with b ( u ) = f ( u ) [ 0 , e h ] T, y ( 0 ) = y 0 = 0, and G as in Eq. (66) for a sinusoidal driving force f ( t ) = F 0 sin ω t acting on node h, h = 1 , , n, is given by the following components:
x ( t ) = i = 1 n F 0 ϕ i ( h ) ξ i [ ( 1 ( λ i ) 2 ξ i ( ω 2 + ( λ i ) 2 ) + 1 ( λ i + ) 2 ξ i ( ω 2 + ( λ i + ) 2 ) sin ω t + ( 1 ω 2 + ( λ i ) 2 1 ω 2 + ( λ i + ) 2 ) ω cos ω t + ( ω ω 2 + ( λ i + ) 2 e λ i + t ω ω 2 + ( λ i ) 2 e λ i t ) ] ϕ i
and
v ( t ) = i = 1 n F 0 ϕ i ( h ) ξ i [ ( λ i ( 1 ( λ i ) 2 ) ξ i ( ω 2 + ( λ i ) 2 ) + λ i + ( 1 ( λ i + ) 2 ) ξ i ( ω 2 + ( λ i + ) 2 ) sin ω t + ( λ i ω 2 + ( λ i ) 2 λ i + ω 2 + ( λ i + ) 2 ) ω cos ω t + ( ω λ i ω 2 + ( λ i + ) 2 e λ i + t ω λ i + ω 2 + ( λ i ) 2 e λ i t ) ( λ i ( 1 ( λ i ) 2 ) ξ i ( ω 2 + ( λ i ) 2 ) + λ i + ( 1 ( λ i + ) 2 ) ξ i ( ω 2 + ( λ i + ) 2 ) ] ϕ i ,
where ξ i = μ i 2 4, λ i + = 1 2 ( ξ i μ i ), and λ i = 1 2 ( ξ i + μ i ) as in Eq. (26).

Notice that, with respect to Eq. (26), we dropped the index G in the eigenvalues to make the expressions lighter.

Leaving the proof of this proposition to  Appendix B, let us now focus on the x component and its behavior for t large. After some further steps, it can be written as
x ( t ) = i = 1 n F 0 ϕ i ( h ) ( ω 2 + ( λ i ) 2 ) ( ω 2 + ( λ i + ) 2 ) [ ( 1 ω 2 ) sin ω t μ i ω cos ω t ω ξ i + ω ξ i ( ( ω 2 + ( λ i ) 2 ) e λ i + t ( ω 2 + ( λ i + ) 2 ) e λ i t ) ] ϕ i ,
and, by components,
x k ( t ) = i = 1 n 1 F 0 ϕ i ( h ) ϕ i ( k ) ( ω 2 + ( λ i ) 2 ) ( ω 2 + ( λ i + ) 2 ) [ ( 1 ω 2 ) sin ω t ω ξ i μ i ω cos ω t + ω ξ i ( ( ω 2 + ( λ i ) 2 ) e λ i + t ( ω 2 + ( λ i + ) 2 ) e λ i t ) ] + F 0 n ω sin t sin ω t ω 2 1
(68)
formula in which we made the final term for i = n explicit.

The question now is what frequency or frequencies are able to induce resonance phenomena in the network as t grows to infinity, bearing in mind the fact that the presence of dissipative terms still drives synchronization between nodes. This amounts to asking what are the differences between formula (59) and formula (68).

First, let us note that the ratio before the square bracket in formula (68) is always a real value. This is due to the fact that λ i and λ i + are reciprocal and to the equality ( λ i ) 2 + ( λ i + ) 2 = μ i 2 2. Moreover, as previously discussed, the exponential terms e λ i ± t inside the brackets go to zero, as t grows. Therefore, in general, for t + , we have coefficients in the sum depending on k, which prevents nodes from synchronizing.

However, when ω 1, the last term in the sum becomes dominant, and this term is common to all the nodes. Indeed, the only frequency ω that can prompt resonance phenomena is the one such that ω 2 + ( λ i ) 2 = 0 and ω 2 + ( λ i + ) 2 = 0 and the only value that can make these quantities vanishing is the pure imaginary value λ n + = λ n = i.

We conclude that, under these particular conditions, the only resonant state is exactly the one involving complete synchronization between the nodes. Therefore, in general, a driving force does not favor the synchronization process and, vice versa, a resonant state is, in general, prevented by the presence of this kind of coupling for frequencies other than ω = 1. There is, thus, only one fully synchronized and resonant state for ω 1 and t + . This state is described by the following asymptotic behavior, still linearly growing in time:
x k ( t ) F 0 n ω sin t sin ω t ω 2 1 F 0 2 n ( sin t t cos t ) .

In Fig. 13, we illustrate the position plots of the four nodes in the toy network when an external force f ( t ) = sin ω t is applied to node number 1 starting from initial conditions x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T. As ω approaches 1 from below, all the nodes are activated and enter a global resonant state, as proven above. In Fig. 14, we illustrate the position plots of the four nodes under the same conditions as before, but as ω approaches, for instance, 2 from below; let us notice that, as expected, no resonant state arises within the network, in contrast with what shown in Fig. 12 under the same conditions.

FIG. 13.

Approaching the resonant state with a driving force in node number 1 of the toy network with dissipative couplings and x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T; the frequencies are panel (a) ω = 0.80, panel (b) ω = 0.95, and panel (c) ω = 0.99.

FIG. 13.

Approaching the resonant state with a driving force in node number 1 of the toy network with dissipative couplings and x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T; the frequencies are panel (a) ω = 0.80, panel (b) ω = 0.95, and panel (c) ω = 0.99.

Close modal
FIG. 14.

Approaching the non-resonant state with a driving force in node number 1 of the toy network with dissipative couplings and x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T; the frequencies are panel (a) ω = 1.20, panel (b) ω = 1.35, and panel (c) ω = 1.40.

FIG. 14.

Approaching the non-resonant state with a driving force in node number 1 of the toy network with dissipative couplings and x 0 = [ 0 , 0 , 0 , 0 ] T and v 0 = [ 0 , 0 , 0 , 0 ] T; the frequencies are panel (a) ω = 1.20, panel (b) ω = 1.35, and panel (c) ω = 1.40.

Close modal
Let us now return to system (12) and consider the case in which we set c 1 = 0, c 2 = 1, c 1 = γ R +, c 2 = 0 in the matrix G. Instead of considering a driving force f ( t ) acting on a single node as before, we now give the factor b ( t ) the form of a constant vector. Therefore, let us consider the system
{ y ˙ = G y + b , y ( 0 ) = y 0 ,
(69)
G being explicitly the matrix
G := [ 0 I L γ I ] ,
(70)
and b ( t ) = [ 0 , p ] T, where p is a constant vector in R n. Now, system (69) is equivalent to the equation
x ¨ = γ x ˙ L x + p ,
(71)
which is called the linearized swing equation.18,20 The swing equation plays a central role in the analysis of power systems dynamics and in a variety of different contexts, such as the description of superconducting Josephson junctions or a nonlinear second-order phase-locking loop in electronics. It governs the rotational dynamics of synchronous machines in stability studies and it is nonlinear in nature, and a closed formal solution cannot be explicitly given. Equation (71), in particular, models the dynamics of power grid networks, where nodes can be either generators (generating power) or loads (consuming power), and edges are transmission lines connecting them. The requirement for a synchronous generator is that it does not lose phase synchronization when small temporary changes occur in the output of the machine. Typical examples are small perturbations in the scheduled generation of the machine, which result in a change in its rotor angle, or small loads added to the network. The response of a power system to impacts is oscillatory, and if the oscillations are damped, the system can reach in time a steady state. Transient stability studies determine whether or not synchronization is maintained after the machine has been subjected to a transient disturbance. Each variable x i in Eq. (71) represents the corresponding node displacement and v i = x ˙ i its nodal velocity relative to a rotational reference value; γ, which we assume common to all the nodes, is the damping coefficient and the vector p has components equal (or, more, in general, proportional) to the small power deviation generated (if positive) or consumed (if negative) by nodes. The general solution of Eq. (69) can be written as
y ( t ) = y ~ + e G t ( y 0 y ~ ) ,
(72)
where y ( 0 ) = y 0 (and typically, in this context, y 0 = 0) and y ~ is such that b = G y ~ or, equivalently, such that y ~ = [ x ~ , v ~ ] T with p = L x ~ and v ~ = 0.

We now limit the discussion to balanced power grids; that is, we assume that i = 1 n p i = 0. Furthermore, we suppose that the power grid is connected so that the multiplicity of the null eigenvalue of L is equal to 1. Under these conditions, the system p = L x ~ is consistent since rank ( L | p ) = rank ( L ). The family of solutions is characterized by the fact that if x ~ 1 and x ~ 2 are two of them, then x ~ 1 x ~ 2 = w u, where w R and u R n is the vector of all 1’s. Since G [ u , 0 ] T = [ 0 , 0 ] T and e G t [ u , 0 ] T = [ u , 0 ] T, then the general solution y ( t ) = e G t y 0 + y ~ e G t y ~ is the same for any w R. This means that we can use any solution of the linear system p = L x ~ to build the general solution in Eq. (72).

Let us illustrate these ideas on the toy network in Fig. 1, postponing the study of a real power grid to Sec. VII. Let us suppose that node 3 is a generator and nodes 1, 2, and 4 are loads and that the corresponding powers are p = [ 0.50 , 0.20 , 1.05 , 0.35 ] T. Then, a solution of p = L x ~ is, for instance, x ~ = [ 0.25 , 0.15 , 0.15 , 0.20 ] T. By substituting y ~ = [ x ~ , 0 ] T and y 0 = [ 0 , 0 ] T in Eq. (71), we get the results depicted in Fig. 15 for two different values of the damping parameter γ.

FIG. 15.

Evolution of the displacement deviation of the toy network nodes for two different damping coefficients: panel (a) γ = 1 and panel (b) γ = 0.4.

FIG. 15.

Evolution of the displacement deviation of the toy network nodes for two different damping coefficients: panel (a) γ = 1 and panel (b) γ = 0.4.

Close modal

Let us imagine that in a social network, there is a debate around a certain topic and different members have different opinions on that topic. Take, for example, the case in which there are different opinions about the qualities of a particular sportsman or of a given member within the same network.

Each member can shape a positive or negative assessment with all possible gradations from full agreement to full disagreement, and he or she can share this opinion with his or her neighbors. We can model this static situation at a fixed time, for instance, by assigning + 1 to perfect agreement and 1 to total disagreement and providing an intermediate graduation of the level of disagreement or agreement continuously from 1 to + 1. The value 0 would represent the totally neutral position.

However, each member is allowed to change his or her opinion over time, adjusting the assessment made at a given instant. Our goal here is to figure out what entrainment action some nodes may play with respect to others in this opinion formation process.

Following, for instance, Cartwright and Harary37 and Altafini,38 it is reasonable to assume that if a member has a close friend, he is much more likely to be influenced by the latter’s opinion than by others’ and therefore, he will move in accordance with his close neighbors’ opinion dynamics. In other words, he will be induced to synchronize the direction and the intensity of the velocity at which he changes his opinion according to that of his friends.

The mapping between the previous theoretical model and the social problem under examination is based on the interpretation of the continuous variables x i as a measure of the level of agreement or disagreement of each member with the subject. The opinion entrainment effect induced by a node on its neighbors—and, thus, on the network—is then modeled by the couplings described above. What we intend to focus on here is not the numerical value of each individual variable at a given instant. Rather, we are interested in the collective behavior of the set of variables x i over time. Such behavior makes clear, as will be seen, the effect of the topological structure of the network in shaping processes of synchronization of members’ ideas and opinions.

A celebrated example of a social network in network theory is provided by the Zachary Club network, see Fig. 16, in which 34 members found themselves having to decide whether to agree with the position of the instructor ( + 1) or with the position of the president ( 1). We can start by assigning an initial “position” vector in such a way that node number 1 in red, the instructor, has position + 1 and node number 34 in blue, the president, has position 1. At the beginning, at time t = 0, all the other nodes are supposed to be totally neutral since they still have to build an opinion, and therefore, they are assigned a “position” equal to 0.

FIG. 16.

Zachary club network partition.

FIG. 16.

Zachary club network partition.

Close modal

We start by assuming that the social system evolves according to Eq. (28) with a matrix G as in Eq. (24). The plot in Fig. 17, panel (a), illustrates the position vs time of two different groups of nodes with initial conditions equal to + 1 for node number 1 (the instructor), 1 for node number 34 (the president), and 0 for all the other nodes. Specifically, we choose four nodes for each one of the two communities illustrated in Fig. 16 and colored in red and blue, respectively, for the group of the instructor and the group of the president. The four nodes in red are node numbers 4 , 11 , 13, and 18, and the four nodes in blue are node numbers 23 , 24 , 30, and 33. As can be seen, nodes in the two different groups tend to synchronize each other from the very beginning. As time goes on, both converge to zero due to the balanced initial conditions that see neither leader prevail over the other. What is interesting and unexpected is the fact that each group seems to assume a position, which is opposite to that of their respective leader. However, it could be easily explained by saying that if node 1 is initially in position + 1, then it can move only in the negative verse toward 0, dragging in the negative verse also nodes that are under its main influence and similarly for node 34.

FIG. 17.

Synchronization in the two different communities: panel (a) initial positions are equal to + 1 for node number 1, 1 for node number 34, and 0 for all the other nodes; panel (b) initial velocities are equal to + 1 for node number 1, 1 for node number 34, and 0 for all the other nodes.

FIG. 17.

Synchronization in the two different communities: panel (a) initial positions are equal to + 1 for node number 1, 1 for node number 34, and 0 for all the other nodes; panel (b) initial velocities are equal to + 1 for node number 1, 1 for node number 34, and 0 for all the other nodes.

Close modal

The plot in Fig. 17, panel (b), illustrates the position vs time of two groups, including the two respective leaders. In this case, initial conditions see all the nodes having the same initial position equal to 0 (including nodes 1 and 34), but we imagine now that at time t = 0, node 1 starts to move with velocity + 1 (in the positive verse) and node 34 starts to move with velocity 1 (in the negative verse). In the plot, nodes in red are numbers 1 , 4 , 11, and 18 and nodes in blue are numbers 34 , 30 , 32, and 33; dashed lines represent the two leader. As can be seen, nodes in the two different groups tend to synchronize each other, to synchronize with their leader and to stay out-of-phase with members of the other group for all the time, until both converge to 0 due to the balanced initial conditions.

Two synchronization processes with unbalanced initial conditions are illustrated in Fig. 18 for the same set of nodes, as in Fig. 17: in panel (a), node number 1 and number 34 have initial positions equal to 0.8 and 0.5, respectively; in panel (b), node number 1 and number 34 have initial velocities equal to 0.8 and 0.5, respectively.

FIG. 18.

Synchronization in the two different communities: panel (a) node number 1 and number 34 initial positions are equal to 0.8 and 0.5, respectively; panel (b) node number 1 and number 34 initial velocities are equal to 0.8 and 0.5, respectively.

FIG. 18.

Synchronization in the two different communities: panel (a) node number 1 and number 34 initial positions are equal to 0.8 and 0.5, respectively; panel (b) node number 1 and number 34 initial velocities are equal to 0.8 and 0.5, respectively.

Close modal

We compare now the synchronization times given by Eq. (37) under different conditions. In detail, we are interested in understanding which of the two leaders is more effective in terms of speed of synchronization. Keeping in mind the discussion in the last remark in Subsec. IV C, we start with strong coupling between nodes c 2 = 1. If only the instructor (node number 1) is endowed with an initial velocity (in our simulation equal to 4) and all other nodes (including number 34) are initially held at zero, then the time at which all nodes synchronize with the asymptotic solution in Eq. (31) below a threshold ε = 0.001 is equal to t = 91.70 with a mean time equal to t = 5.29. If we now assign the same initial velocity to the president (node 34) while keeping all others (including node number 1) at zero, then the same time becomes t = 96.20 with a mean time equal to t = 5.47. The two times are basically comparable with slightly higher effectiveness for the instructor than for the president, confirming the fact that the network is divided into two balanced communities. The same conclusion is validated using other initial conditions on both the velocity and the position of the two leaders.

Let us now move to the resonant phenomena. We assume now that the social system evolves according to Eq. (56) with a matrix G as in Eq. (54). The resonant frequencies given in Proposition 3 range from 1 to 4.37, and they are listed in Table IV.

TABLE IV.

Eigenvalues of the matrix L and corresponding resonance frequencies for the Zachary club network. Note that eigenvalue 2.000 and the corresponding frequency have multiplicity equal to 5.

i μi ωi i μi ωi i μi ωi
18.137  4.374 55  11  4.581  2.362 37  21  2.000  1.732 05 
17.055  4.249 14  12  4.480  2.340 94  22  1.955  1.719 03 
13.306  3.782 34  13  4.276  2.296 93  23  1.826  1.681 09 
10.921  3.452 69  14  3.472  2.114 76  24  1.762  1.661 90 
9.777  3.282 87  15  3.382  2.093 32  25  1.599  1.612 23 
6.996  2.827 76  16  3.376  2.091 93  26  1.259  1.503 13 
6.516  2.741 45  17  3.242  2.059 63  27  1.125  1.457 74 
6.332  2.707 69  18  3.014  2.003 49  28  0.909  1.381 76 
5.618  2.572 55  19  2.749  1.936 27  29  0.468  1.211 83 
10  5.379  2.525 59  20  2.487  1.867 38  30  1.000 00 
i μi ωi i μi ωi i μi ωi
18.137  4.374 55  11  4.581  2.362 37  21  2.000  1.732 05 
17.055  4.249 14  12  4.480  2.340 94  22  1.955  1.719 03 
13.306  3.782 34  13  4.276  2.296 93  23  1.826  1.681 09 
10.921  3.452 69  14  3.472  2.114 76  24  1.762  1.661 90 
9.777  3.282 87  15  3.382  2.093 32  25  1.599  1.612 23 
6.996  2.827 76  16  3.376  2.091 93  26  1.259  1.503 13 
6.516  2.741 45  17  3.242  2.059 63  27  1.125  1.457 74 
6.332  2.707 69  18  3.014  2.003 49  28  0.909  1.381 76 
5.618  2.572 55  19  2.749  1.936 27  29  0.468  1.211 83 
10  5.379  2.525 59  20  2.487  1.867 38  30  1.000 00 

Let us imagine now that an external force with frequency ω is applied to node number 1 or, in other words, that the first leader acts as an influencer by applying the same driving force in that point of the network. Let us consider two nodes belonging to the two different clusters previously considered: node number 11 in the first cluster (red colored) and node number 30 in the second cluster (blue colored). Let us consider their behavior as ω approaches two resonant frequencies, specifically ω 30 = 1 and ω 22 = 1.719 026. Results are illustrated in Figs. 19 and 20.

FIG. 19.

Approaching the resonant state with a driving force in node number 1 of the Zachary network and with x 0 = 0 and v 0 = 0; the frequencies are panel (a) ω = 0.80, panel (b) ω = 0.95, and panel (c) ω = 0.99.

FIG. 19.

Approaching the resonant state with a driving force in node number 1 of the Zachary network and with x 0 = 0 and v 0 = 0; the frequencies are panel (a) ω = 0.80, panel (b) ω = 0.95, and panel (c) ω = 0.99.

Close modal
FIG. 20.

Approaching the resonant state with a driving force in node number 1 of the Zachary network and with x 0 = 0 and v 0 = 0; the frequencies are panel (a) ω = 1.50, panel (b) ω = 1.65, and panel (c) ω = 1.71.

FIG. 20.

Approaching the resonant state with a driving force in node number 1 of the Zachary network and with x 0 = 0 and v 0 = 0; the frequencies are panel (a) ω = 1.50, panel (b) ω = 1.65, and panel (c) ω = 1.71.

Close modal

As discussed in the remark after Proposition 3, all nodes resonate at the ground frequency ω = 1. Node number 11, which belongs to the sphere of influence of the first leader, certainly resonates at the ground frequency and so does at the second frequency under consideration; on the contrary, node number 30, which falls within the sphere of influence of the second leader, while behaving similarly to node 11 at the ground frequency, does not respond to the second frequency from node 1. This is easily explained if we look at the component of the corresponding eigenvector ϕ 22. Indeed, for node number 11, it is ϕ 22 ( 11 ) = 0.540 899 038, whereas for node number 30, it is ϕ 22 ( 30 ) = 0.002 692 861. The first component is about 200 times larger than the second one, which justifies the very low influence of node 1 on node 30 at that frequency.

Let us consider now the frequency ω 9 = 2.572 554. Since ϕ 9 ( 1 ) = 0, we expect that the oscillation of node 1 with that frequency does not affect the network. This is confirmed by Fig. 21, where we selected two arbitrary nodes (numbers 13 and 29) in the two different clusters, which show no amplification of their own oscillation when the frequency of node 1 approaches the indicated value.

FIG. 21.

Approaching the resonant state with a driving force in node number 1 of the Zachary network and with x 0 = 0 and v 0 = 0; the frequencies are panel (a) ω = 2.560, panel (b) ω = 2.570, and panel (c) ω = 2.572.

FIG. 21.

Approaching the resonant state with a driving force in node number 1 of the Zachary network and with x 0 = 0 and v 0 = 0; the frequencies are panel (a) ω = 2.560, panel (b) ω = 2.570, and panel (c) ω = 2.572.

Close modal

Further confirmation comes from the frequency ω 6 = 2.827 755. Since ϕ 6 ( 1 ) = 0.0027808652 and ϕ 6 ( 4 ) = 0.823172 4260, we expect that the network globally responds much more if that oscillation starts from node 4 than from node 1. Figure 22 shows the different behavior of two randomly selected nodes (again numbers 13 and 29) when such a frequency is applied on the network starting from different source nodes, specifically nodes 1, 4, and 34. The amplitude of the oscillation of nodes 13 and 29 grows much more quickly when the source is placed in node 4 than when it is placed in node 34 or, worse, in node 1. The two leaders prove to be not so effective at this specific frequency.

FIG. 22.

Approaching the resonant state with a driving force of fixed frequency ω = 2.827 75 and with x 0 = 0 and v 0 = 0; the source nodes are node number 1 in panel (a), node number 4 in panel (b), and node number 34 in panel (c).

FIG. 22.

Approaching the resonant state with a driving force of fixed frequency ω = 2.827 75 and with x 0 = 0 and v 0 = 0; the source nodes are node number 1 in panel (a), node number 4 in panel (b), and node number 34 in panel (c).

Close modal

Finally, let us consider the case of the resonant synchronization described by proposition 4 by assuming now that the social system evolves according to Eq. (68) with a matrix G as in Eq. (66). Figure 23 shows the behavior of the two leaders when the first one approaches frequency ω = 1 from below. All the other nodes display similar behavior and enter a synchronized resonant state with the two leaders. No other frequency is able to produce a state of synchronized resonance on all nodes in the network.

FIG. 23.

Approaching the synchronized resonant state with a driving force in node number 1 of the Zachary network and x 0 = 0 and v 0 = 0; the frequencies are panel (a) ω = 0.90, panel (b) ω = 0.98, and panel (c) ω = 0.99.

FIG. 23.

Approaching the synchronized resonant state with a driving force in node number 1 of the Zachary network and x 0 = 0 and v 0 = 0; the frequencies are panel (a) ω = 0.90, panel (b) ω = 0.98, and panel (c) ω = 0.99.

Close modal
We apply now results discussed in Sec. VI to a small power grid, namely, to the main electricity transmission network of Syria, represented in Fig. 24, panel (a).39 A power grid consists of a network whose nodes are rotating machines and edges are transmission lines. Nodes can produce power (generating plants) or consume power (loads). The network topology of the Syrian grid is depicted in Fig. 24, panel (b). It consists of 19 nodes with 6 generators (nodes 1–6, in red), 13 loads (nodes 7–19, in blue), and 24 edges.
FIG. 24.

(a) The main structure of the Syrian power grid [H. Beltaifa, Arab Union of Electricity, Syrian Arab Republic; see https://energydata.info/dataset/syria-electricity-transmission-network-2017 for Syria—Electricity Transmission Network (last accessed April 13, 2020). Copyright 2017 Author(s), licensed under a Creative Commons Attribution 4.0 License.] and (b) the topology of the same network: nodes in red are generators, and nodes in blue are loads.

FIG. 24.

(a) The main structure of the Syrian power grid [H. Beltaifa, Arab Union of Electricity, Syrian Arab Republic; see https://energydata.info/dataset/syria-electricity-transmission-network-2017 for Syria—Electricity Transmission Network (last accessed April 13, 2020). Copyright 2017 Author(s), licensed under a Creative Commons Attribution 4.0 License.] and (b) the topology of the same network: nodes in red are generators, and nodes in blue are loads.

Close modal

We aim primarily to show the qualitative behavior of the solution described in Eq. (72) for specific nodes and for different values of the involved parameters. In Fig. 25, we represent the dynamic response of the selected nodes in the Syrian power grid: nodes 3, 4, and 5 correspond to the main power plants in terms of expected capacity, namely, Baniyas, Zayzoun, and Al-Thawrah plants, whereas nodes 9, 15, and 18 are located in three major Syrian cities, namely, Damascus, Palmyra, and Al-Hasakah. In panel (a), we consider a distribution of nominal balanced powers such that p i = 0.13 for generators and p i = 0.06 for loads with damping coefficient γ = 1; in panel (b), we keep the same powers but with damping coefficient γ = 0.4; and in panel (c), we assume p = [ 0.004 , 0.181 , 0.252 , 0.202 , 0.306 , 0.056 ] for the six generators and p i = 0.077 for all the motors. Figure 25 qualitatively illustrates some important features of the transient dynamics of x i ( t ), such as the first peak value, the first peak time, and the value at a steady state under different conditions (see Bhatta et al.18 for a comprehensive discussion about these values). Table V collects the values of these parameters for two representative nodes in the examined simulation and under conditions corresponding to panels in Fig. 25.

For instance, as can be observed, the plant corresponding to node 4 typically has the highest peak value, while node 3 has the shortest peak time among the examined nodes. For the sake of brevity, we do not consider a power distribution that respects the symmetries of the network, which is such that all nodes within the same symmetry cluster have the same power. In this case, it has been shown that if the vector