On the reliable and efficient numerical integration of the Kuramoto model and related dynamical systems on graphs

Abstract

In this work, a novel approach for the reliable and efficient numerical integration of the Kuramoto model on graphs is studied. For this purpose, the notion of order parameters is revisited for the classical Kuramoto model describing all-to-all interactions of a set of oscillators. First numerical experiments confirm that the precomputation of certain sums significantly reduces the computational cost for the evaluation of the right-hand side and hence enables the simulation of high-dimensional systems. In order to design numerical integration methods that are favourable in the context of related dynamical systems on network graphs, the concept of localized order parameters is proposed. In addition, the detection of communities for a complex graph and the transformation of the underlying adjacency matrix to block structure is an essential component for further improvement. It is demonstrated that for a submatrix comprising relatively few coefficients equal to zero, the precomputation of sums is advantageous, whereas straightforward summation is appropriate in the complementary case. Concluding theoretical considerations and numerical comparisons show that the strategy of combining effective community detection algorithms with the localization of order parameters potentially reduces the computation time by several orders of magnitude.

Keywords:

• Differential equations
• dynamical systems
• network dynamics
• Kuramoto model
• Kuramoto model on graphs
• numerical integration
• geometric integration

Mathematics Subject Classifications:

• 37M05
• 34A34
• 65L05
• 65P10

1. Introduction

In the present work, we propose a novel approach for the reliable and efficient numerical integration of nonlinear dynamical systems on network graphs and provide various numerical comparisons confirming its potential. In essence, our objective is to combine effective algorithms for the detection of communities with the concept of localized order parameters based on the precomputation of certain sums. For the sake of concreteness, we focus on Kuramoto-type models for a large set of individual oscillators and graphs that are determined by adjacency matrices with coefficients equal to one and zero, respectively. Generally, the interactions between the oscillators are described by a system of nonlinear ordinary differential equations involving the sines of the associated phases. The classical Kuramoto–Daido model [9,23] reflects the case of all-to-all coupling and thus corresponds to a complete graph. Realistic extended models incorporate complex networks [32]. Yet, the efficient numerical integration of a Kuramoto-type model on a graph comprising a high number of nodes is a relevant issue, and even the evaluation of the vector field defining the right-hand side of the system poses a major challenge. This problem actually permeates the simulation of complex dynamics on networks [28]. In order to exemplify our strategy and to illustrate its capability, we initially consider the classical Kuramoto–Daido model and revisit the well-known notion of order parameters [2,9,23,26,37,38]. An essential component of our procedure for Kuramoto-type models on graphs is the detection of communities [13,29], since this permits a transformation of the associated adjacency matrix to block structure. We demonstrate that for a submatrix comprising relatively few coefficients equal to zero the precomputation of sums is indeed advantageous and that straightforward summation is suitable in the complementary case. We conclude with theoretical considerations and numerical comparisons, which show that our approach ensures a significant reduction of the required memory capacity as well as the computational cost for the evaluation of the right-hand side. As a consequence, long-term simulations of systems involving a large number of oscillators by geometric integrators [4,19,34,35] are within reach.

Scope of the model. The starting point of our investigations is the Kuramoto–Daido model [9,23]. We henceforth refer to it as (classical) Kuramoto model. This system of coupled nonlinear ordinary differential equations is a fundamental mathematical model for the dynamical behaviour of a set of weakly coupled, nearly identical oscillators and specifies the time evolution of the associated phases. Despite its simple structure, the Kuramoto model exhibits fascinating phenomena such as synchronization and phase locking [37]. Originally introduced to describe processes in chemistry and biology [23,40], it was found to have various applications in other fields such as physics, neuroscience, and engineering [1,10,14]. From a theoretical perspective, the Kuramoto model has deep connections to effects present in Hamiltonian systems, particularly to Landau damping [11,12] and bifurcations from essential spectra [6]. A variety of Kuramoto-type models have a gradient flow structure entailing further interesting mathematical cross-connections [16].

Kuramoto model. The main terms defining the right-hand side of the classical Kuramoto model for M individual oscillators have the form (1) $\frac{1}{M}\sum _{\ell =1}^M\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),t\in [0,T],m\in \{1,2,\dots ,M\}.$(1) Here, ${\vartheta }_m$ denotes the time-dependent phase of the m-th oscillator, which takes values in the circle ${\mathbb{S}}_1=\mathbb{R}/(2\pi \mathbb{Z})$. Henceforth, we employ the convenient vector notation $\vartheta =({\vartheta }_1,\dots ,{\vartheta }_M{)}^T:[0,T]⟶{\mathbb{S}}_1^M.$With regard to the numerical simulation of a high number of oscillators, an essential requirement is the efficient evaluation of these sums at certain time grid points. Thereby, one issue is the limited memory capacity. As an example, we mention the currently widely used software Matlab, which has a maximum array size preference of about ${10}^{10}$ (74.5 GB) corresponding to a square matrix of dimension ${10}^5$. Consequently, it is desirable to avoid the creation of the matrix $({\vartheta }_{\ell }(t)-{\vartheta }_m(t){)}_{\ell ,m\in \{1,2,\dots ,M\}},t\in [0,T],$since this would restrict the dimension of the system considerably. An important alternative for the efficient evaluation of (1(1) $\frac{1}{M}\sum _{\ell =1}^M\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),t\in [0,T],m\in \{1,2,\dots ,M\}.$(1) ) relies on the macroscopic order parameters, which are given through $r(\vartheta (t)){\mathrm{e}}^{\mathrm{i}\psi (\vartheta (t))}=\frac{1}{M}\sum _{m=1}^M{\mathrm{e}}^{\mathrm{i}{\vartheta }_m(t)},t\in [0,T].$More precisely, applying the addition theorem for the sine function to (1(1) $\frac{1}{M}\sum _{\ell =1}^M\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),t\in [0,T],m\in \{1,2,\dots ,M\}.$(1) ), we obtain the following reformulation (2) $\begin{array}{rl}& S_M(\vartheta (t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t)),C_M(\vartheta (t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t)),\\ & \frac{1}{M}\sum _{\ell =1}^M\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t))=S_M(\vartheta (t))\cos ({\vartheta }_m(t))-C_M(\vartheta (t))\sin ({\vartheta }_m(t)),\\ & t\in [0,T],m\in \{1,2,\dots ,M\}.\end{array}$(2) An evident though crucial observation is that precomputing the sums $S_M(\vartheta (t))$ and $C_M(\vartheta (t))$ permits to evaluate $S_M(\vartheta (t))$ $\cos ({\vartheta }_m(t))-C_M(\vartheta (t))\sin ({\vartheta }_m(t))$ for $m\in \{1,2,\dots ,M\}$ in an efficient manner. Numerical comparisons described in Section 2 confirm that this approach reduces the number of function evaluations and thus the computation time considerably.

Kuramoto model on graphs. As indicated above, it is of high relevance to study extensions of the classical Kuramoto model in the context of dynamical networks [32]. We focus on situations, where the sum over all phases is replaced by a sum over certain phases. That is, the right-hand side of the system involves terms of the form $\sum _{\ell =1}^M{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),A_{m\ell }\in \{0,1\},t\in [0,T],\ell ,m\in \{1,2,\dots ,M\}.$The associated matrix $A=(A_{m\ell }{)}_{\ell ,m\in \{1,\dots ,M\}}\in {\mathbb{R}}^{M\times M}$has the natural interpretation as the adjacency matrix of a graph. In Section 3, we consider the complementary cases of sparse and dense adjacency matrices. Numerical tests for randomly generated matrices show that the precomputation of sums is advantageous whenever the number of coefficients equal to one is larger than the number of coefficients equal to zero, whereas straightforward summation is appropriate otherwise. Subsequently, we examine algorithms detecting communities in a graph, which yield as outputs partitions of the nodes, and extend our approach to block matrices. In our considerations, we do not presume a symmetric adjacency matrix, which would lead to a gradient system.

Generalisations. We point out that our approach applies to the more general case, where the coefficients of the adjacency matrix take values in a finite set, for instance $A_{m\ell }\in \{0,1,\dots ,J\},\ell ,m\in \{1,2,\dots ,M\}.$Provided that a suitable permutation permits the transformation to a block matrix such that a single value is prevalent in each block and the total number of blocks is relatively small, we may expect a significant gain in efficiency by the precomputation of sums. Moreover, it is straightforward to generalize our approach to systems that in addition involve multiple sums such as $\sum _{j,k,\ell =1}^M{A}_{mjk\ell }\sin ({\vartheta }_j(t)+{\vartheta }_k(t)-{\vartheta }_{\ell }(t)-{\vartheta }_m(t)),t\in [0,T],m\in \{1,2,\dots ,M\}.$Due to the fact that higher-order Kuramoto models describe interactions beyond usual graph structures, they have recently raised remarkable interest, see [3,36] and references given therein. However, as their incorporation requires laborious notation and would obstruct a comprehensible presentation of the key idea, we do not include detailed calculations here.

Remarks. We note that suitable reformulations of the classical Kuramoto model and Kuramoto-type models on graphs can also be based on complex exponentials. In view of expedient practical implementations, however, it is preferable to use real-valued quantities. On account of cross-references in captions, figures are included at the end of this manuscript.

2. Kuramoto model

In this section, we state the classical Kuramoto model and study different viewpoints on its reliable and efficient numerical integration.

Original formulation. We consider a set of M limit-cycle oscillators with time-dependent phases (3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) In the absence of external driving or damping forces, respectively, the oscillators have the intrinsic frequencies (3b) ${\omega }_m\in \mathbb{R},m\in \{1,2,\dots ,M\}.$(3b) The pairwise interactions between the oscillators are described by the following system of nonlinear ordinary differential equations (3c) $\{\begin{array}{l}{\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m+\frac{K}{M}\sum _{\ell =1}^M\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),}\\ {\displaystyle {\vartheta }_m(0)\mathrm{given},t\in (0,T),m\in \{1,2,\dots ,M\},}\end{array}$(3c) where K>0 denotes the coupling constant.

Special choices. In our numerical tests, we consider intrinsic frequencies defined by a single real number ${\omega }_0\in \mathbb{R}$ and initial phases of the form (3d) ${\omega }_m=1+{\omega }_0\frac{(2m-M-1)}{M-1},{\vartheta }_m(0)=\frac{2\pi m}{M},m\in \{1,2,\dots ,M\}.$(3d)

Reformulation. Regarding the efficient numerical integration of the classical Kuramoto model, it is to the best advantage to employ a reformulation that relies on elementary addition theorems for trigonometric functions $\begin{array}{rl}\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t))& =\sin ({\vartheta }_{\ell }(t))\cos ({\vartheta }_m(t))-\cos ({\vartheta }_{\ell }(t))\sin ({\vartheta }_m(t)),\\ & \ell ,m\in \{1,2,\dots ,M\}.\end{array}$Introducing the abbreviations (4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) the governing equations (3a(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) ) read as (4b) $\{\begin{array}{l}{\vartheta }_m^{\prime }(t)={\omega }_m+K(S_M(\vartheta (t))\cos ({\vartheta }_m(t))-C_M(\vartheta (t))\sin ({\vartheta }_m(t))),\\ {\vartheta }_m(0)\mathrm{given},t\in (0,T),m\in \{1,2,\dots ,M\}.\end{array}$(4b) Employing the standard vector notation (4c) $\omega =\left(\begin{array}{c}{\omega }_1\\ {\omega }_2\\ ⋮\\ {\omega }_M\end{array}\right)\in {\mathbb{R}}^M,\vartheta (t)=\left(\begin{array}{c}{\vartheta }_1(t)\\ {\vartheta }_2(t)\\ ⋮\\ {\vartheta }_M(t)\end{array}\right)\in {\mathbb{S}}_1^M,t\in [0,T],$(4c) and setting accordingly (4d) $\begin{array}{rl}\sin (\vartheta (t))& =\left(\begin{array}{c}\sin ({\vartheta }_1(t))\\ \sin ({\vartheta }_2(t))\\ ⋮\\ \sin ({\vartheta }_M(t))\end{array}\right),\cos (\vartheta (t))=\left(\begin{array}{c}\cos ({\vartheta }_1(t))\\ \cos ({\vartheta }_2(t))\\ ⋮\\ \cos ({\vartheta }_M(t))\end{array}\right),\\ F(\vartheta (t))& =\omega +K(S_M(\vartheta (t))\cos (\vartheta (t))-C_M(\vartheta (t))\sin (\vartheta (t))),t\in [0,T],\end{array}$(4d) the system takes the compact form (4e) $\{\begin{array}{l}{\vartheta }^{\prime }(t)=F(\vartheta (t)),t\in (0,T),\\ \vartheta (0)\mathrm{given}.\end{array}$(4e)

Potential. The classical Kuramoto model (3(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) )–(4(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ) has the intrinsic structure of a gradient system. That is, the right-hand side is given by the gradient of a real-valued potential function $\begin{array}{rl}-\mathrm{\nabla }V(\vartheta )& =({\omega }_m+\frac{K}{M}\sum _{\ell =1}^M\sin ({\vartheta }_{\ell }-{\vartheta }_m){)}_{m\in \{1,2,\dots ,M\}}\\ & =\omega +K(S_M(\vartheta )\cos (\vartheta )-C_M(\vartheta )\sin (\vartheta )),\vartheta =({\vartheta }_1,\dots ,{\vartheta }_M{)}^T\in {\mathbb{S}}_1^M,\end{array}$such that the governing equations rewrite as $\{\begin{array}{l}{\vartheta }^{\prime }(t)=-\mathrm{\nabla }V(\vartheta (t)),t\in (0,T),\\ \vartheta (0)\mathrm{given}.\end{array}$This, in particular, implies that the values of the potential decrease when time evolves $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}V(\vartheta (t))& =(\mathrm{\nabla }V(\vartheta (t)){)}^T{\vartheta }^{\prime }(t)=-{\Vert \mathrm{\nabla }V(\vartheta (t))\Vert }^2\le 0,\\ V(\vartheta (t))& \le V(\vartheta (0)),t\in [0,T].\end{array}$For our purposes, in view of its efficient evaluation, it is advantageous to reformulate the canonical potential as follows (5a) $\begin{array}{rl}V:{\mathbb{R}}^M& ⟶\mathbb{R}:\\ \vartheta & ⟼V(\vartheta )=V({\vartheta }_1,\dots ,{\vartheta }_M)\\ & =-\sum _{m=1}^M{\omega }_m{\vartheta }_m+\frac{K}{2M}\sum _{\ell ,m=1}^M(1-\cos ({\vartheta }_{\ell }-{\vartheta }_m))\\ & =-{\omega }^T\vartheta +\frac{KM}{2}(1-(C_M(\vartheta ){)}^2-(S_M(\vartheta ){)}^2).\end{array}$(5a) As mentioned below, the choice of the constant of integration is linked to the special case of synchronization.

Order parameter. The modulus $r:{\mathbb{S}}_1^M\to \mathbb{R}$ and the angle $\psi :{\mathbb{S}}_1^M\to \mathbb{R}$ of the complex order parameter are given by (5b) $r(\vartheta ){\mathrm{e}}^{\mathrm{i}\psi (\vartheta )}=\frac{1}{M}\sum _{m=1}^M{\mathrm{e}}^{\mathrm{i}{\vartheta }_m}=C_M(\vartheta )+\mathrm{i}{S}_M(\vartheta ),\vartheta \in {\mathbb{S}}_1^M.$(5b) Multiplying by ${\mathrm{e}}^{-\mathrm{i}{\vartheta }_{\ell }}$ for $\ell \in \{1,2,\dots ,M\}$ and considering the imaginary part of the resulting relation $\begin{array}{rl}& r(\vartheta ){\mathrm{e}}^{\mathrm{i}(\psi (\vartheta )-{\vartheta }_{\ell })}=\frac{1}{M}\sum _{m=1}^M{\mathrm{e}}^{\mathrm{i}({\vartheta }_m-{\vartheta }_{\ell })},r(\vartheta )\sin (\psi (\vartheta )-{\vartheta }_{\ell })=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m-{\vartheta }_{\ell }),\\ & \vartheta \in {\mathbb{S}}_1^M,\ell \in \{1,2,\dots ,M\},\end{array}$the Kuramoto model (3(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) )–(4(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ) rewrites as $\{\begin{array}{l}{\vartheta }^{\prime }(t)=\omega +Kr(\vartheta (t))\sin (\psi (\vartheta (t))-\vartheta (t)),\\ \vartheta (0)\mathrm{given},t\in (0,T).\end{array}$We point out that the order parameters contain the entire information about the interactions of the oscillators. Although this reformulation looks like a mean-field equation for a single oscillator, it completely represents the original system.

Indicator for synchronization. For configurations, where all cosine and sine values are close-by, the modulus of the complex order parameter has values nearly one $\begin{array}{rl}& \cos ({\vartheta }_m)\approx \cos ({\vartheta }_1),\sin ({\vartheta }_m)\approx \sin ({\vartheta }_1),m\in \{2,3,\dots ,M\},\\ & C_M(\vartheta )\approx \cos ({\vartheta }_1),S_M(\vartheta )\approx \sin ({\vartheta }_1),\\ & r(\vartheta )=\sqrt{(C_M(\vartheta ){)}^2+(S_M(\vartheta ){)}^2}\approx 1,\vartheta \in {\mathbb{S}}_1^M.\end{array}$Hence, this quantity indicates synchronization. Furthermore, in this situation, the above-stated choice of the constant of integration in the potential implies $V(\vartheta )\approx -{\omega }^T\vartheta ,\vartheta \in {\mathbb{S}}_1^M.$

Conserved quantity. A straightforward calculation shows that summation over all governing equations yields the identity $\frac{1}{M}\sum _{m=1}^M{\vartheta }_m^{\prime }(t)=\frac{1}{M}\sum _{m=1}^M{\omega }_m,t\in [0,T],$see (3(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) )–(4(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ). Performing integration, this implies that the mean values of the intrinsic frequencies and the initial phases determine the mean values of the phases at later times $\frac{1}{M}\sum _{m=1}^M{\vartheta }_m(t)=\frac{1}{M}\sum _{m=1}^M{\vartheta }_m(0)+t\frac{1}{M}\sum _{m=1}^M{\omega }_m,t\in [0,T].$In other words, the dynamics of the classical Kuramoto model is restricted to a time-dependent submanifold defined by the constraint (5c) $\frac{1}{M}\sum _{m=1}^M({\vartheta }_m(t)-{\vartheta }_m(0)-t{\omega }_m)=0,t\in [0,T].$(5c) We anticipate that this characteristic result extends to Kuramoto models on graphs under a certain symmetry condition, see Section 3. However, in general, this conservation property does not hold and then the quantity (5d) $\frac{1}{M}\sum _{m=1}^M({\vartheta }_m(t)-{\vartheta }_m(0)-t{\omega }_m)$(5d) reflects the deviation.

In the limit $M\to \mathrm{\infty }$, the sums in (5c(5c) $\frac{1}{M}\sum _{m=1}^M({\vartheta }_m(t)-{\vartheta }_m(0)-t{\omega }_m)=0,t\in [0,T].$(5c) )–(5d(5d) $\frac{1}{M}\sum _{m=1}^M({\vartheta }_m(t)-{\vartheta }_m(0)-t{\omega }_m)$(5d) ) are replaced by integrals. In the case of randomly chosen problem data, e.g. they are interpreted as expected phases and intrinsic frequency, respectively.

Implementation and computational cost. For systems involving a high number of oscillators $M>>1$, the above-stated approach leading to a reformulation of the classical Kuramoto model is beneficial in several respects. It permits to reduce the required memory capacity as well as the computational cost for the evaluation of the right-hand side significantly, see Table 1. Moreover, parallelization techniques can be used.

Table 1. Classical Kuramoto model.

Original formulation	$M(M-1)$
Reformulation	$4M$

Computational cost for the evaluation of the right-hand side based on the original formulation (3a(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) ) and the reformulation (4a(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ), respectively. Number of function evaluations (sine, cosine) in dependence of the total number of oscillators.

Exemplification (Euler). In order to exemplify our procedure for the efficient numerical integration of the Kuramoto model, we consider the simplest first-order one-step method, the explicit Euler method. For a time grid with associated stepsizes $0=t_0<\cdots <t_n<\cdots <t_N=T,{\tau }_n=t_{n+1}-t_n,n\in \{0,1,\dots ,N-1\},$and a prescribed initial approximation ${\vartheta }^{(0)}\approx \vartheta (0)$, the explicit Euler solution is given by the recurrence ${\vartheta }^{(n+1)}={\vartheta }^{(n)}+{\tau }_{n}F({\vartheta }^{(n)}),n\in \{0,1,\dots ,N-1\},$see (4(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ). In each time step, the evaluation of the defining function relies on the precomputation of the sums (6a) $S_M^{(n)}=S_M({\vartheta }^{(n)}),C_M^{(n)}=C_M({\vartheta }^{(n)}),$(6a) and subsequently on the computation of (6b) $F({\vartheta }^{(n)})=\omega +K(S_M^{(n)}\cos ({\vartheta }^{(n)})-C_M^{(n)}\sin ({\vartheta }^{(n)})).$(6b) Altogether, this requires $4M$ evaluations of sine and cosine functions compared to $M(M-1)$ evaluations of the sine function needed for the original formulation (3(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) ). The generalization to higher-order explicit or implicit time integration methods is straightforward. Long-term simulations are ideally based on geometric integrators. Their benefits over standard methods are demonstrated in [4,19,34,35], e.g.

Remark 2.1

In our study, the focus is on the numerical simulation of a high number of oscillators. In this situation, the reduction from $M^2$ to $M^2-M$ evaluations of the sine function by taking into account the evident identity for coinciding indices $\ell =m:\sin ({\vartheta }_{\ell }-{\vartheta }_m)=0$is of minor relevance and thus will be neglected.

Numerical comparisons. Numerical comparisons of different approaches for the evaluation of the right-hand side of the classical Kuramoto model (3(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) )–(4(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ) and the associated potential (5(5a) $\begin{array}{rl}V:{\mathbb{R}}^M& ⟶\mathbb{R}:\\ \vartheta & ⟼V(\vartheta )=V({\vartheta }_1,\dots ,{\vartheta }_M)\\ & =-\sum _{m=1}^M{\omega }_m{\vartheta }_m+\frac{K}{2M}\sum _{\ell ,m=1}^M(1-\cos ({\vartheta }_{\ell }-{\vartheta }_m))\\ & =-{\omega }^T\vartheta +\frac{KM}{2}(1-(C_M(\vartheta ){)}^2-(S_M(\vartheta ){)}^2).\end{array}$(5a) ) are displayed in Figure 1. We focus on an implementation in Matlab and expect that analogous conclusions hold for other software packages. We vary the total number of oscillators from ${10}^2$ to ${10}^8$, taking into account maximum array size preferences as mentioned in the introduction. A randomly chosen number ${\omega }_0\in (0,1)$ defines the intrinsic frequencies. The phases ${\vartheta }_1,{\vartheta }_2,\dots ,{\vartheta }_M$ are uniformly distributed in $[0,2\pi ]$. We compare

1. straightforward summation of $\sin ({\vartheta }_{\ell }-{\vartheta }_m)$ for $\ell ,m\in \{1,2,\dots ,M\}$ realized by a double loop over all rows,

2. a simple script using $\mathrm{sum}(\sin (\vartheta -{\vartheta }_m))$ for $m\in \{1,2,\dots ,M\}$ and thus involving a loop over all rows,

3. a script by Cleve Moler that generates the matrix $\sin (\vartheta -{\vartheta }^{\prime })$ and then adds up each row, ¹

4. the above-stated procedure using the precomputation of sums, see (4(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ) and (6(6a) $S_M^{(n)}=S_M({\vartheta }^{(n)}),C_M^{(n)}=C_M({\vartheta }^{(n)}),$(6a) ).

Figure 1. Classical Kuramoto model. Computational cost for the evaluation of the right-hand side and the potential. Left: Number of sine and cosine evaluations versus the total number of oscillators when using straightforward summation and the precomputation of sums, respectively. Middle: Numerical comparison of the computation time for different implementations in Matlab. Right: Corresponding results for the potential.

For a higher number of oscillators, it is probable that the evaluation of functions and the computation of sums are the most time consuming components. The numerical results confirm that our approach (iv) is favourable in this regime and permits the efficient simulation of a high number of oscillators. Similar conclusions hold for the evaluation of the potential.

Numerical integration. The numerical integration of the classical Kuramoto model (3(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) )–(4(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ) based on the precomputation of sums is illustrated in Figures 2–6. Movies showing the time evolution are found at http://techmath.uibk.ac.at/mecht/MyHomepage/Research/MovieKuramotoClassicalK1.m4v http://techmath.uibk.ac.at/mecht/MyHomepage/Research/MovieKuramotoClassicalK3.m4v http://techmath.uibk.ac.at/mecht/MyHomepage/Research/MovieKuramotoClassicalK5.m4v

Figure 2. Numerical integration of the classical Kuramoto model involving $M={10}^2$ oscillators. Coupling constant K = 1 (no synchronization). Top: Visualisation of the phases at the final time. Middle: The time series confirms decreasing potential values. Bottom: A conserved quantity is numerically preserved with high accuracy.

Figure 3. Numerical integration of the classical Kuramoto model involving $M={10}^2$ oscillators. Coupling constant $K\in \{3,5\}$ (gradual synchronization).

Figure 4. Numerical integration of the classical Kuramoto model involving $M={10}^4$ oscillators. Coupling constant K = 1 (no synchronization).

Figure 5. Numerical integration of the classical Kuramoto model involving $M={10}^4$ oscillators. Coupling constant $K\in \{3,5\}$ (gradual synchronization).

Figure 6. Long-term integration of the classical Kuramoto model based on a second-order explicit Runge–Kutta method (top) and a second-order implicit Runge–Kutta method with improved numerical preservation of a conserved quantity (bottom). To avoid a further diminishment of the relevant vertical axis, ticks along the horizontal axis are omitted.

We consider $M={10}^2$ as well as $M={10}^4$ oscillators and set T = 200. The interplay between the coupling constant $K\in \{1,3,5\}$ and the constant ${\omega }_0=2$ defining the intrinsic frequencies determines the strength of synchronization. The points and the arrow reflect the values of the phases and the complex order parameter at the final time $(\cos ({\vartheta }_m(t)-\psi (t)),\sin ({\vartheta }_m(t)-\psi (t))),m\in \{1,2,\dots ,M\},r(\vartheta (t)),t=T,$see (5b(5b) $r(\vartheta ){\mathrm{e}}^{\mathrm{i}\psi (\vartheta )}=\frac{1}{M}\sum _{m=1}^M{\mathrm{e}}^{\mathrm{i}{\vartheta }_m}=C_M(\vartheta )+\mathrm{i}{S}_M(\vartheta ),\vartheta \in {\mathbb{S}}_1^M.$(5b) ). Moreover, the graphs of the associated quantities $V(\vartheta (t)),\frac{1}{M}\sum _{m=1}^M({\vartheta }_m(t)-{\vartheta }_m(0)-t{\omega }_m),t\in [0,T],$are shown, see (5a(5a) $\begin{array}{rl}V:{\mathbb{R}}^M& ⟶\mathbb{R}:\\ \vartheta & ⟼V(\vartheta )=V({\vartheta }_1,\dots ,{\vartheta }_M)\\ & =-\sum _{m=1}^M{\omega }_m{\vartheta }_m+\frac{K}{2M}\sum _{\ell ,m=1}^M(1-\cos ({\vartheta }_{\ell }-{\vartheta }_m))\\ & =-{\omega }^T\vartheta +\frac{KM}{2}(1-(C_M(\vartheta ){)}^2-(S_M(\vartheta ){)}^2).\end{array}$(5a) ) and (5d(5d) $\frac{1}{M}\sum _{m=1}^M({\vartheta }_m(t)-{\vartheta }_m(0)-t{\omega }_m)$(5d) ). In addition, the computation time CT, measured in seconds, is displayed. For the considered time interval, a variable stepsize fourth-order explicit Runge–Kutta method leads to a reliable result. In particular, it is seen that the values of the potential decrease when time evolves and that the modulus of the complex order parameter is close to one for K = 5. When applying instead a second-order geometric integrator, we observed an improved behaviour regarding the conserved quantity over long times, see Figure 6. For comprehensive information about the benefits of geometric numerical integrators in comparison with standard integration methods, we refer to [4,19,34,35].

3. Kuramoto models on graphs

In this section, we investigate extensions of the classical Kuramoto model on graph topologies and propose suitable modifications of our approach for their efficient numerical integration.

General formulation. Henceforth, we study the following system of ordinary differential equations (7a) $\{\begin{array}{l}{\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m+\frac{K}{{\mathcal{M}}_m}\sum _{\ell =1}^M{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),}\\ {\vartheta }_m(0)\mathrm{given},t\in (0,T),m\in \{1,2,\dots ,M\}.\end{array}$(7a) In contrast to the classical Kuramoto model (3(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) ), the all-to-all coupling is replaced by the interactions of certain communities of oscillators, which are described by the associated adjacency matrix (7b) $A=(A_{m\ell }{)}_{\ell ,m\in \{1,\dots ,M\}},A_{m\ell }\in \{0,1\},\ell ,m\in \{1,2,\dots ,M\}.$(7b) We distinguish two kinds of scalings affecting in particular the strength of synchronization. In case the mth row of A is zero, the corresponding equation reduces to ${\vartheta }_m^{\prime }(t)={\omega }_m,t\in (0,T),$and hence the scaling is dispensible.

1. Uniform scaling.   With regard to the literature, a common choice is (7c) ${\mathcal{M}}_m=M,m\in \{1,2,\dots ,M\}.$(7c)

2. Non-uniform scaling.   As exemplified below, an alternative is (7d) ${\mathcal{M}}_m=\sum _{\ell =1}^M{A}_{m\ell },m\in \{1,2,\dots ,M\}.$(7d)

For our considerations, it is convenient to introduce the scaled adjacency matrix (7e) $\mathcal{A}={\left(\frac{1}{{\mathcal{M}}_m}{A}_{m\ell }\right)}_{\ell ,m\in \{1,\dots ,M\}}.$(7e)

Auxiliary identities. The decisive terms defining the right-hand side of the extended Kuramoto model (7a(7a) $\{\begin{array}{l}{\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m+\frac{K}{{\mathcal{M}}_m}\sum _{\ell =1}^M{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),}\\ {\vartheta }_m(0)\mathrm{given},t\in (0,T),m\in \{1,2,\dots ,M\}.\end{array}$(7a) ) rewrite as follows $\begin{array}{rl}& \sum _{\ell =1}^M{\mathcal{A}}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t))\\ & =\sum _{\ell =1}^M{\mathcal{A}}_{m\ell }\sin ({\vartheta }_{\ell }(t))\cos ({\vartheta }_m(t))\\ & -\sum _{\ell =1}^M{\mathcal{A}}_{m\ell }\cos ({\vartheta }_{\ell }(t))\sin ({\vartheta }_m(t)),t\in [0,T],m\in \{1,2,\dots ,M\}.\end{array}$Denoting the componentwise product of two columns by $v.\ast w=\left(\begin{array}{c}{v}_1\\ v_2\\ ⋮\\ v_M\end{array}\right).\ast \left(\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_M\end{array}\right)=\left(\begin{array}{c}{v}_1{w}_1\\ v_2{w}_2\\ ⋮\\ v_M{w}_M\end{array}\right),v,w\in {\mathbb{R}}^M,$and employing the compact vector notation (8a) $\begin{array}{r}G(\vartheta (t))=\omega +K((\mathcal{A}\sin (\vartheta (t))).\ast \cos (\vartheta (t))-(\mathcal{A}\cos (\vartheta (t))).\ast \sin (\vartheta (t)),\\ t\in [0,T],\end{array}$(8a) we obtain the following reformulation of the governing equations (8b) $\{\begin{array}{l}{\vartheta }^{\prime }(t)=G(\vartheta (t)),t\in (0,T),\\ \vartheta (0)\mathrm{given},\end{array}$(8b) see also (4(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ). Similar considerations yield the identity $\begin{array}{rl}& \sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }\cos ({\vartheta }_{\ell }(t)-{\vartheta }_m(t))\\ & =\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }\cos ({\vartheta }_{\ell }(t))\cos ({\vartheta }_m(t))+\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }\sin ({\vartheta }_{\ell }(t))\sin ({\vartheta }_m(t))\\ & =(\cos (\vartheta (t)){)}^T(\mathcal{A}\cos (\vartheta (t)))+(\sin (\vartheta (t)){)}^T(\mathcal{A}\sin (\vartheta (t))),t\in [0,T].\end{array}$Uniform scaling. In situations, where the adjacency matrix is symmetric and the scaling is uniform $A=A^T,{\mathcal{M}}_m=M,m\in \{1,2,\dots ,M\},\mathcal{A}=\frac{1}{M}A,{\mathcal{A}}^T=\mathcal{A},$the existence of a potential function is ensured and the sum over all phases leads to a conserved quantity, see also (5(5a) $\begin{array}{rl}V:{\mathbb{R}}^M& ⟶\mathbb{R}:\\ \vartheta & ⟼V(\vartheta )=V({\vartheta }_1,\dots ,{\vartheta }_M)\\ & =-\sum _{m=1}^M{\omega }_m{\vartheta }_m+\frac{K}{2M}\sum _{\ell ,m=1}^M(1-\cos ({\vartheta }_{\ell }-{\vartheta }_m))\\ & =-{\omega }^T\vartheta +\frac{KM}{2}(1-(C_M(\vartheta ){)}^2-(S_M(\vartheta ){)}^2).\end{array}$(5a) ). With regard to the classical Kuramoto model and the above stated auxiliary identity, we define the potential through (9a) $\begin{array}{rl}V:{\mathbb{R}}^M& ⟶\mathbb{R}:\\ \vartheta & ⟼V(\vartheta )=V({\vartheta }_1,\dots ,{\vartheta }_M)\\ & =-\sum _{m=1}^M{\omega }_m{\vartheta }_m+\frac{K}{2}\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }(1-\cos ({\vartheta }_{\ell }-{\vartheta }_m))\\ & =-{\omega }^T\vartheta +\frac{K}{2}(\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }-(\cos (\vartheta ){)}^T(\mathcal{A}\cos (\vartheta ))\\ & -(\sin (\vartheta ){)}^T(\mathcal{A}\sin (\vartheta ))).\end{array}$(9a) Due to the fact that the relation $\begin{array}{r}{\mathcal{A}}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t))+{\mathcal{A}}_{\ell m}\sin ({\vartheta }_m(t)-{\vartheta }_{\ell }(t))=0,t\in [0,T],\ell ,m\in \{1,2,\dots ,M\},\end{array}$holds, summation and integration with respect to time implies (9b) $\frac{1}{M}\sum _{m=1}^M({\vartheta }_m(t)-{\vartheta }_m(0)-t{\omega }_m)=0,t\in [0,T].$(9b)

Non-uniform scaling. In order to illustrate our motive for the non-uniform scaling, we consider the special case, where the adjacency matrix comprises a square submatrix with coefficients equal to one and is zero otherwise $\begin{array}{rl}A& =\left(\begin{array}{cc}{B}^{(11)}& B^{(12)}\\ B^{(21)}& B^{(22)}\end{array}\right)\in {\mathbb{R}}^{M\times M},\\ B^{(11)}& =(1{)}_{\ell ,m\in \{1,\dots ,M_0\}},B^{(12)}=B^{(21)}=B^{(22)}=0\in {\mathbb{R}}^{M_0\times M_0},\\ M_0& \in \{2,\dots ,M-1\}.\end{array}$This would correspond to a configuration with pairwise interaction of the first part of the oscillators and a decoupling of the second part $\{\begin{array}{l}{\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m+\frac{K}{{\mathcal{M}}_m}\sum _{\ell =1}^{M_0}\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),t\in (0,T),m\in \{1,2,\dots ,M_0\},}\\ {\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m,t\in (0,T),m\in \{M_0+1,M_0+2,\dots ,M\}.}\end{array}$Here, the non-uniform scaling seems to be more natural ${\mathcal{M}}_m=\sum _{\ell =1}^M{A}_{m\ell }=M_0\ne M,m\in \{1,2,\dots ,M_0\}.$We point out that a potential function of the form (9(9a) $\begin{array}{rl}V:{\mathbb{R}}^M& ⟶\mathbb{R}:\\ \vartheta & ⟼V(\vartheta )=V({\vartheta }_1,\dots ,{\vartheta }_M)\\ & =-\sum _{m=1}^M{\omega }_m{\vartheta }_m+\frac{K}{2}\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }(1-\cos ({\vartheta }_{\ell }-{\vartheta }_m))\\ & =-{\omega }^T\vartheta +\frac{K}{2}(\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }-(\cos (\vartheta ){)}^T(\mathcal{A}\cos (\vartheta ))\\ & -(\sin (\vartheta ){)}^T(\mathcal{A}\sin (\vartheta ))).\end{array}$(9a) ) remains valid for non-symmetric scaled adjacency matrices, whereas summation over all governing equations does not lead to a conserved quantity, in general.

Extension of our approach. The extension of our approach for the classical Kuramoto model based on the precomputation of sums to Kuramoto models on graphs (7(7a) $\{\begin{array}{l}{\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m+\frac{K}{{\mathcal{M}}_m}\sum _{\ell =1}^M{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),}\\ {\vartheta }_m(0)\mathrm{given},t\in (0,T),m\in \{1,2,\dots ,M\}.\end{array}$(7a) ) requires a careful incorporation of the structure of the associated adjacency matrix. The following considerations prove to be particularly expedient for situations, where the adjacency matrix has a block structure and can be divided into relatively sparse and relatively dense submatrices, respectively. In view of an efficient evaluation of the decisive terms in (7a(7a) $\{\begin{array}{l}{\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m+\frac{K}{{\mathcal{M}}_m}\sum _{\ell =1}^M{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),}\\ {\vartheta }_m(0)\mathrm{given},t\in (0,T),m\in \{1,2,\dots ,M\}.\end{array}$(7a) ), we consider a submatrix of the form $(A_{m\ell }{)}_{m\in \{M_1,\dots ,M_2\},\ell \in \{M_3,\dots ,M_4\}},$defined by positive integers $M_1,M_2,M_3,M_4\in \mathbb{N}$ such that $1\le M_1<M_2\le M$ as well as $1\le M_3<M_4\le M$. Our basic concept is to optimize the required memory capacity and the number of functions evaluations. More precisely, in order to avoid the storage of each submatrix and to accelerate the computation of the sums $\sum _{\ell =M_3}^{M_4}{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),m\in \{M_1,\dots ,M_2\},$we distinguish two complementary cases. Whenever the number of coefficients equal to zero is relatively high, we store all indices corresponding to non-zero coefficients and sum over these coefficients (10a) $\underset{A_{m\ell }=1}{\sum _{\ell \in \{M_3,\dots ,M_4\}}}{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),m\in \{M_1,\dots ,M_2\}.$(10a) Whenever the number of non-zero coefficients is relatively high, we instead store all indices corresponding to coefficients equal to zero and make use of the precomputation of sums. That is, where applicable, we employ the reformulation $\begin{array}{rl}& \sum _{\ell =M_3}^{M_4}{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t))\\ & =\sum _{\ell =M_3}^{M_4}{A}_{m\ell }\sin ({\vartheta }_{\ell }(t))\cos ({\vartheta }_m(t))-\sum _{\ell =M_3}^{M_4}{A}_{m\ell }\cos ({\vartheta }_{\ell }(t))\sin ({\vartheta }_m(t))\\ & =(\sum _{\ell =M_3}^{M_4}\sin ({\vartheta }_{\ell }(t))-\sum _{\ell =M_3}^{M_4}(1-A_{m\ell })\sin ({\vartheta }_{\ell }(t)))\cos ({\vartheta }_m(t))\\ & -(\sum _{\ell =M_3}^{M_4}\cos ({\vartheta }_{\ell }(t))-\sum _{\ell =M_3}^{M_4}(1-A_{m\ell })\cos ({\vartheta }_{\ell }(t)))\sin ({\vartheta }_m(t)),\end{array}$precompute the sums (10b) $S_{M_3,M_4}(\vartheta (t))=\sum _{\ell =M_3}^{M_4}\sin ({\vartheta }_{\ell }(t)),C_{M_3,M_4}(\vartheta (t))=\sum _{\ell =M_3}^{M_4}\cos ({\vartheta }_{\ell }(t)),$(10b) and then substract the terms that correspond to coefficients equal to zero (10c) $\begin{array}{rl}& \sum _{\ell =M_3}^{M_4}{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t))\\ & =(S_{M_3,M_4}(\vartheta (t))-\underset{A_{m\ell }=0}{\sum _{\ell \in \{M_3,\dots ,M_4\}}}\sin ({\vartheta }_{\ell }(t)))\cos ({\vartheta }_m(t))\\ & -(C_{M_3,M_4}(\vartheta (t))-\underset{A_{m\ell }=0}{\sum _{\ell \in \{M_3,\dots ,M_4\}}}\cos ({\vartheta }_{\ell }(t)))\sin ({\vartheta }_m(t)),m\in \{M_1,\dots ,M_2\}.\end{array}$(10c)

Remark 3.1

In situations, where it is desirable to evaluate the right-hand side and the potential in parallel, it is advantageous to modify the above stated approach in view of an efficient computation of $\begin{array}{rl}& \sum _{\ell =M_3}^{M_4}{A}_{m\ell }\sin ({\vartheta }_{\ell }(t))=\sum _{\ell =M_3}^{M_4}\sin ({\vartheta }_{\ell }(t))-\sum _{\ell =M_3}^{M_4}(1-A_{m\ell })\sin ({\vartheta }_{\ell }(t))\\ & \sum _{\ell =M_3}^{M_4}{A}_{m\ell }\cos ({\vartheta }_{\ell }(t))=\sum _{\ell =M_3}^{M_4}\cos ({\vartheta }_{\ell }(t))-\sum _{\ell =M_3}^{M_4}(1-A_{m\ell })\cos ({\vartheta }_{\ell }(t)),\\ & m\in \{M_1,\dots ,M_2\},\end{array}$either by straightforward summation or based on the precomputation of sums, see also (8(8a) $\begin{array}{r}G(\vartheta (t))=\omega +K((\mathcal{A}\sin (\vartheta (t))).\ast \cos (\vartheta (t))-(\mathcal{A}\cos (\vartheta (t))).\ast \sin (\vartheta (t)),\\ t\in [0,T],\end{array}$(8a) )–(9(9a) $\begin{array}{rl}V:{\mathbb{R}}^M& ⟶\mathbb{R}:\\ \vartheta & ⟼V(\vartheta )=V({\vartheta }_1,\dots ,{\vartheta }_M)\\ & =-\sum _{m=1}^M{\omega }_m{\vartheta }_m+\frac{K}{2}\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }(1-\cos ({\vartheta }_{\ell }-{\vartheta }_m))\\ & =-{\omega }^T\vartheta +\frac{K}{2}(\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }-(\cos (\vartheta ){)}^T(\mathcal{A}\cos (\vartheta ))\\ & -(\sin (\vartheta ){)}^T(\mathcal{A}\sin (\vartheta ))).\end{array}$(9a) ).

Numerical comparisons. In order to test the above described strategy (10(10a) $\underset{A_{m\ell }=1}{\sum _{\ell \in \{M_3,\dots ,M_4\}}}{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),m\in \{M_1,\dots ,M_2\}.$(10a) ) for the efficient computation of the sums $\sum _{\ell =1}^M{A}_{m\ell }\sin ({\vartheta }_{\ell }-{\vartheta }_m),{\vartheta }_m=\frac{2\pi m}{M},m\in \{1,2,\dots ,M\},$we study two situations, which we consider to be of practical relevance in view of the numerical integration of Kuramoto models on graphs (3(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) ), see also (7(7a) $\{\begin{array}{l}{\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m+\frac{K}{{\mathcal{M}}_m}\sum _{\ell =1}^M{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),}\\ {\vartheta }_m(0)\mathrm{given},t\in (0,T),m\in \{1,2,\dots ,M\}.\end{array}$(7a) ).

1. Matrices without block structure.   On the one hand, we define well-balanced adjacency matrices and compare the number of function evaluations as well as the computation time of an implementation in Matlab. For this purpose, we prescribe a threshold $p\in [0,1]$ and generate row-by-row a sequence of uniformly distributed random numbers $z\in [0,1]$. Whenever z>p, the corresponding coefficient is set to one, otherwise, it is set to zero. For p = 0.99, e.g. the resulting adjacency matrix is sparse and straightforward summation is advantageous, see Figure 7. For the complementary case p = 0.01, e.g. the adjacency matrix comprises few non-zero coefficients and the precomputation of sums is favourable, see Figure 9. For the threshold p = 0.5 both approaches lead to essentially the same counts, see Figure 8.

2. Matrices with block structure.   On the other hand, we generate adjacency matrices that are composed of well-balanced submatrices. Each block is connected with a different threshold. We use

   1. a double loop over all rows and straightforward summation of $\sin ({\vartheta }_{\ell }-{\vartheta }_m)$ provided that $A_{m\ell }\ne 0$ for $\ell ,m\in \{1,2,\dots ,M\}$,

   2. a single loop to sum over all indices that correspond to non-zero coefficients, realized by a script of the form $\mathrm{sum}(\sin (\vartheta (\mathrm{NonZero})-{\vartheta }_m))$ for $m\in \{1,2,\dots ,M\}$,

   3. the precomputation of sums and a single loop to subtract terms that correspond to coefficients equal to zero, applying the script $\mathrm{sum}(\sin (\vartheta (\mathrm{Zero})-{\vartheta }_m))$ for $m\in \{1,2,\dots ,M\}$,

   4. the procedure in (b) and (c) adapted to each submatrix.

The obtained results, displayed in Figures 10–12, confirm that the latter approach is beneficial for a higher number of oscillators, where the evaluation of functions and the computation of sums are expected to be the most time consuming components.

Figure 7. Extended Kuramoto models involving randomly generated adjacency matrices defined through the thresholds 0.99, 0.9. Evaluation of the right-hand side by means of approaches adapted to sparse matrices (straightforward summation) and dense matrices (precomputation of sums), respectively. Left: Illustration of the adjacency matrix for 100 oscillators. Middle: Number of sine and cosine evaluations versus the total number of oscillators. Right: Comparison of the computation time.

Figure 8. Corresponding results for the thresholds 0.7, 0.5, 0.3.

Figure 9. Corresponding results for the thresholds 0.1, 0.01.

Figure 10. Extended Kuramoto models involving randomly generated adjacency matrices with block structure defined through a threshold per block. Evaluation of the right-hand side by means of different approaches adapted to sparse matrices, dense matrices, and block matrices, respectively. Left: Illustration of the adjacency matrix for 100 oscillators. Middle: Number of sine and cosine evaluations versus the total number of oscillators. Right: Comparison of the computation time.

Figure 11. Corresponding results.

Figure 12. Corresponding results for a more realistic adjacency matrix describing the interactions of four communities of oscillators.

Community detection in graphs. Our numerical comparisons show that it is to the best advantage to take the underlying structure of the considered Kuramoto model on a graph (7(7a) $\{\begin{array}{l}{\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m+\frac{K}{{\mathcal{M}}_m}\sum _{\ell =1}^M{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),}\\ {\vartheta }_m(0)\mathrm{given},t\in (0,T),m\in \{1,2,\dots ,M\}.\end{array}$(7a) ) into account. A suitable reordering of the governing equations accordingly to the separation of the oscillators into communities corresponds to the transformation of the associated adjacency matrix to a well-balanced block matrix and is a fundamental means in view of their efficient numerical integration, see Figures 13 and 14. In the following, we study the performance of various algorithms for the detection of communities in graphs. We restrict ourselves to algorithms from the Python packages networkx and cdlib, which were suitable for our purposes, did not require additional parameters, and took a reasonable computation time, see Table 2. We focus on the relevant case, where interactions primarily take place within certain communities of oscillators. As illustrated in Figure 15, our starting point is a matrix (of suitably chosen dimension) with fully occupied blocks along the diagonal $\begin{array}{rl}B& =(B_{m\ell }{)}_{m,\ell \in \{1,2,\dots ,M\}}=\left(\begin{array}{cccc}{B}^{(11)}& & & 0\\ & B^{(22)}& & \\ & & B^{(33)}& \\ 0& & & B^{(44)}\end{array}\right)\in {\mathbb{R}}^{M\times M},\\ B^{(kk)}& =(1{)}_{i,j\in \{1,2,\dots ,(5-k)M/10\}}\in {\mathbb{R}}^{(5-k)M/10\times (5-k)M/10},k\in \{1,2,3,4\}.\end{array}$

Figure 13. Adjacency matrices A and $PA{P}^T$ with suitably chosen permutation matrix P. Left: The structure of the underlying graph is not evident. Right: A separation into communities of oscillators is recognizable. Matrices of this form are used for numerical tests of community detection algorithms.

Figure 14. Adjacency matrices A and $PA{P}^T$ with suitably chosen permutation matrix P. Left: The structure of the underlying graph is not evident. Right: A separation into communities of oscillators is recognizable. Symmetric adjacency matrices of this form are considered in connection with the numerical integration of Kuramoto models on graphs. See also Figure 12.

Figure 15. Left: Matrix comprising four blocks (p = 0, top) and related adjacency matrices for increasing thresholds $p\in \{0.1,0.2\}$ (middle to bottom). Right: Average of the quantity (11(11) $\sum _{\ell ,m=1}^M\left|{\tilde{C}}_{m\ell }\right|-\sum _{\ell ,m=1}^M\left|C_{m\ell }\right|,$(11) ) over eight runs, which is chosen as a measure for the performance of the considered community detection algorithms, see also Table 2. Obtained results for M = 100 (top), M = 200 (middle), M = 400 (bottom).

Table 2. Algorithms from Python packages applied for community detection in graphs.

networkx [17]	greedy_modularity [7]
cdlib [33]	louvain [5]
	rber_pots [30,31]
	rb_pots [24,31]
	significance_communities [39]
	walktrap [27]

Similarly to before, we prescribe a threshold $p\in [0,1]$ that determines the probability that a coefficient of the matrix B is changed from one to zero or from zero to one, respectively. That is, we generate uniformly distributed random numbers $z_{m\ell }\in [0,1]$ for $m,\ell \in \{1,2,\dots ,M\}$ and define $\begin{array}{r}A=(A_{m\ell }{)}_{m,\ell \in \{1,2,\dots ,M\}}\in {\mathbb{R}}^{M\times M},A_{m\ell }=\{\begin{array}{ll}1-B_{m\ell },& z_{m\ell }<p,\\ B_{m\ell },& z_{m\ell }\ge p,\end{array}m,\ell \in \{1,2,\dots ,M\}.\end{array}$The higher the threshold, the stronger the deviation of the adjacency matrix from the related block matrix, see Figures 15 and 16. The community detection algorithms were applied to the randomly permuted adjacency matrix and yield as outputs partitions of the nodes into communities. We associate them with permutations of the sequence $(1,2,\dots ,M)$ such that nodes in one community are arranged one after the other. These permutations correspond to transformations of the original matrices that can be cast into the form $\tilde{A}=PA{P}^T$with a permutation matrix $P\in {\mathbb{R}}^{M\times M}$ such that in particular the identity $P^{-1}=P^T$ is valid. For each of these matrices, we identify a matrix with fully occupied blocks along the diagonal reflecting the detected communities $\tilde{B}=\left(\begin{array}{ccc}{\tilde{B}}^{(11)}& & 0\\ & {\tilde{B}}^{(22)}& \\ 0& & \ddots \end{array}\right)\in {\mathbb{R}}^{M\times M}.$We recall that the computation of pointwise products such as $(B\sin (\vartheta )).\ast \cos (\vartheta ),(\tilde{B}\sin (\vartheta )).\ast \cos (\vartheta ),\vartheta \in {\mathbb{R}}^M,$based on the precomputation of sums requires in total 2M evaluations of cosine and sine functions and that the number of non-zero coefficients of the matrices $C=A-B,\tilde{C}=\tilde{A}-\tilde{B},$match the additional computational costs. For this reason, we determine the quantity (11) $\sum _{\ell ,m=1}^M\left|{\tilde{C}}_{m\ell }\right|-\sum _{\ell ,m=1}^M\left|C_{m\ell }\right|,$(11) in order to assess the performance of the algorithms, see [29] for detailed explanations. The numerical results obtained for $M\in \{100,200,400,800,1600\}$ and thresholds in the range $[0,0.4]$ are displayed in Figures 15–17. In case p = 0, i.e. for a random permutation of the fully occupied block diagonal matrix, each of the tested algorithms detected the four communities. For larger deviations of the adjacency matrices from the underlying block diagonal matrix, however, some algorithms failed. Overall, the community detection algorithm rber_pots, available through the Python package cdlib, provided the most reliable results for a higher number of oscillators.

Figure 16. Left: Adjacency matrices for increasing thresholds $p\in \{0.3,0.4\}$ (top to bottom). Right: Corresponding results for M = 800 (top), M = 1600 (bottom).

Figure 17. Average computation time of the considered community detection algorithms over eight runs. Left: The displayed numerical results, obtained for the threshold p = 0.2, reflect a quadratic scaling with respect to the total number of oscillators. Right: For larger deviations of the adjacency matrices from the underlying block diagonal matrix, related to increasing values of the threshold p, the computation time increases, in general.

Favourable community detection algorithm. It is notable that there exist several variants of the algorithm rber_pots [30,31] with the common objective to detect communities in a graph with associated adjacency matrix A. The standard implementation in cdlib [33] minimizes the quantity $Q=-\sum _{\ell ,m=1}^M(A_{m\ell }-p)\delta ({\sigma }_m,{\sigma }_{\ell }),$where $p\in (0,1)$ represents the mean density of edges in a graph, that is, the ratio between the numbers of actually existing and potential edges. Whenever two nodes $m,\ell \in \{1,2,\dots ,M\}$ belong to the same community, the quantity $\delta ({\sigma }_m,{\sigma }_{\ell })$ takes the value one, and it is zero otherwise. Provided that the mean edge density within a community is higher than the mean edge density of the complete network, a node is included in a community. The favourable performance of rber_pots in the considered example is explained by the fact that the four communities are well recognizable by the mean edge densities, even for p = 0.4.

Numerical integration. Following up the numerical integration of the classical Kuramoto model (3a(3a) ${\vartheta }_m:[0,T]⟶{\mathbb{S}}_1,m\in \{1,2,\dots ,M\}.$(3a) )–(4a(4a) $\begin{array}{r}{S}_M(\vartheta (t))=S_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\sin ({\vartheta }_m(t))\in \mathbb{R},\\ C_M(\vartheta (t))=C_M({\vartheta }_1(t),{\vartheta }_2(t),\dots ,{\vartheta }_M(t))=\frac{1}{M}\sum _{m=1}^M\cos ({\vartheta }_m(t))\in \mathbb{R},\end{array}$(4a) ), we finally study extended Kuramoto models on graphs (7a(7a) $\{\begin{array}{l}{\displaystyle {\vartheta }_m^{\prime }(t)={\omega }_m+\frac{K}{{\mathcal{M}}_m}\sum _{\ell =1}^M{A}_{m\ell }\sin ({\vartheta }_{\ell }(t)-{\vartheta }_m(t)),}\\ {\vartheta }_m(0)\mathrm{given},t\in (0,T),m\in \{1,2,\dots ,M\}.\end{array}$(7a) ). On the one hand, we consider the situation, where a separation of the oscillators into four communities of the same size is evident. The structure of the associated symmetric adjacency matrix is illustrated in Figure 14 (right). Here, our approach based on the precomputation of sums applies. The numerical results obtained for the common uniform scaling, a total number of $M=8^4=4096$ oscillators, coupling constant K = 3, final time T = 200, and intrinsic frequencies as well as initial phases of the form ${\omega }_0=2,{\omega }_m=1+{\omega }_0\frac{(2m-M-1)}{M-1},{\vartheta }_m(0)=\frac{2\pi m}{M},m\in \{1,2,\dots ,M\},$are displayed in Figure 18 (right). On the other hand, we consider the equivalent system without recognizable block structure, see Figure 14 (left), using straightforward summation for the evaluation of the right-hand side. In order to achieve consistency with the previous case, we permute internal frequencies and initial phases accordingly, perform the time integration, and reorder the solution values subsequently. By comparison of the numerical results shown in Figure 18, it is evident that the efficient evaluation of the decisive sums based on the block structure of the adjacency matrix is beneficial and permits a significant reduction of the computation time from approximately 890 seconds to about 54 seconds. For the sake of comparability with the classical Kuramoto, we display the values of the corresponding potentials, which decrease when time evolves, and verify the conservation property, see also (9(9a) $\begin{array}{rl}V:{\mathbb{R}}^M& ⟶\mathbb{R}:\\ \vartheta & ⟼V(\vartheta )=V({\vartheta }_1,\dots ,{\vartheta }_M)\\ & =-\sum _{m=1}^M{\omega }_m{\vartheta }_m+\frac{K}{2}\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }(1-\cos ({\vartheta }_{\ell }-{\vartheta }_m))\\ & =-{\omega }^T\vartheta +\frac{K}{2}(\sum _{\ell ,m=1}^M{\mathcal{A}}_{m\ell }-(\cos (\vartheta ){)}^T(\mathcal{A}\cos (\vartheta ))\\ & -(\sin (\vartheta ){)}^T(\mathcal{A}\sin (\vartheta ))).\end{array}$(9a) ). The analogous results for the non-uniform scaling are given in Figure 19. As expected, due to the lack of symmetry, the constraint (9b(9b) $\frac{1}{M}\sum _{m=1}^M({\vartheta }_m(t)-{\vartheta }_m(0)-t{\omega }_m)=0,t\in [0,T].$(9b) ) is notfulfilled.

Figure 18. Numerical integration of a Kuramoto model on a graph with the common uniform scaling. The associated adjacency matrix has a structure as illustrated in Figure 14 (left: matrix without recognizable block structure, right: permuted matrix with recognizable block structure). Left: Evaluation of the right-hand side by straightforward summation. Right: Employing the underlying block structure and using the precomputation of sums permits a significant reduction of the computation time CT.

Figure 19. Corresponding results for the non-uniform scaling. Synchronisation within four communities is observed. Due to the lack of symmetry of the system, the conservation property does not hold.

4. Conclusions and outlook

In summary, we have presented results which confirm that the combination of detecting communities and using localisations of order parameters permits significant improvements regarding the reliable and efficient numerical integration of dynamical systems on graphs. Our approach provides a natural way to exploit the underlying graph structure and local mean-field variables.

Popular models, where the employed concepts apply at once, are given by consensus problems [25]. However, as these systems are linear with known analytical solution representations, we have considered their numerical simulation to be of less importance and have focused on Kuramoto models as relevant instances. Generalisations to other nonlinear dynamical systems on graphs occurring in applications will be the objective of future investigations.

As mentioned in the introduction, it is straightforward to extend our approach to oscillator networks with higher-order interactions included [3,36]. Suitable modifications have to be developed for Cucker–Smale models describing flocking behaviour [8]. Other types of models, commonly encountered in neuroscience applications, incorporate individual excitable oscillators beyond a phase reduction [20,22]. A further example is the strategy of combining our approach with micro-macro numerical integration schemes such as projective integration methods [21], where one makes use in the time integration of a macroscopic evolution, e.g. for the order parameter or the probability density of a typical oscillator, in combination with direct microscopic simulation to improve the numerical simulation. In such a context, our improvements directly accelerate the microscopic integrator. Yet, some natural-looking steps are bound to be far more involved. For example, for time-dependent or adaptive network dynamics [15,18], the recomputation of community structures at every iteration is computationally inefficient and requires novel perspectives. Besides, providing a rigorous numerical analysis to assess the quality of community detection algorithms remains an open question. Furthermore, it is of relevance to study dynamical systems on graphs that incorporate stochastic perturbations.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1 The Matlab script kuramoto.m by Cleve Moler is available at https://de.mathworks.com/matlabcentral/fileexchange/72534-kuramoto-s-model-of-synchronizing-oscillators.