Koopman analysis in oscillator synchronization
KKoopman analysis in oscillator synchronization
Jing Hu and Yueheng Lan
1, 2, ∗ School of Science, Beijing University of Posts and Telecommunications, Beijing 100876, China State Key Lab of Information Photonics and Optical Communications,Beijing University of Posts and Telecommunications, Beijing 100876, China (Dated: September 22, 2020)Synchronization is an important dynamical phenomenon in coupled nonlinear systems, which has been stud-ied extensively in recent years. However, analysis focused on individual orbits seems hard to extend to complexsystems while a global statistical approach is overly cursory. Koopman operator technique seems to well balancethe two approaches. In this paper, we extend Koopman analysis to the study of synchronization of coupled oscil-lators by extracting important eigenvalues and eigenfunctions from the observed time series. A renormalizationgroup analysis is designed to derive an analytic approximation of the eigenfunction in case of weak couplingthat dominates the oscillation. For moderate or strong couplings, numerical computation further confirms theimportance of the average frequencies and the associated eigenfunctions. The synchronization transition pointscould be located with quite high accuracy by checking the correlation of neighbouring eigenfunctions at di ff erentcoupling strengths, which is readily applied to other nonlinear systems. Keywords: synchronization, Koopman operator, renormalization group, average frequencies, synchroniza-tion transition point.
PACS numbers: 02.30.-f , 02.70.-c , 05.10.Cc , 05.45.Xt
I. INTRODUCTION
With the rapid development of science and technology,huge amount of data is generated and collected each day, con-cerning with complex systems involving natural phenomena,engineering operations or human activities. Nevertheless, itis nearly impossible to get all the details of a large systemand measurement could only be made for a small numberof observables which are functions of the system state. Inthe theory of low-dimensional nonlinear dynamics, the statespace reconstruction is possible even for a scalar observa-tion [1] since information is highly mixed when the dynamicsis chaotic. However, in high-dimensional systems, the obser-vations are so scant that much information is simply missedfrom the data. On the other hand, measurement is often con-taminated with various noise that could be very disturbing.How to extract important information about the system fromthese partial observations remains a major challenge in boththeory and practice.At this point, the traditional focus on a single orbit maynot be so appealing and it should be more fruitful to view theevolution of all observables as a whole, which may suppresspossible noise and excavate stable information simultaneouslyfrom all the components as well as from their temporal corre-lation. The evolution of an observable as a function in thestate space is described by the Koopman operator proposedby Bernard Koopman and John von Neumann in 1930 [2, 3].However, in recent years, with the rise of computer power anddata science, this alternative framework of nonlinear dynam-ics analysis is extensively explored and a plethora of its appli-cations are found in di ff erent fields [4–8]. In this new formu-lation, we focus on spectral properties of the Koopman oper-ator, especially, the dominant eigenvalues and eigenfunctions. ∗ Email address : [email protected]
In high-dimensional systems, however, the Koopman operatoris usually represented by a very large matrix which has numer-ous eigenvalues and a reliable way is still lacking on how toidentify and interpret the dominant ones, not to mention ananalytic approach to the problem, concerning its complexity.Nevertheless, recent progress in the renormalization group(RG) theory may shed some light, which is capable of extract-ing global information from local expansion [9]. Renormal-ization group investigates the change of physical laws acrossdi ff erent scales and was first proposed for the perturbative cal-culation in quantum physics and statistical physics [10]. Later,it was extended to the treatment of nonlinear dynamics forboth flows and mappings [9, 11, 12]. The starting point ofthe renormalization group is the removal of divergences (orresonance terms) from the perturbation series such that sta-ble characteristics of system structure and dynamics are ex-tracted which are insensitive to details, so it can be regarded asan altermative formulation of asymptotic analysis [9, 11, 12].Usually asymptotic analysis techniques such as multiple scalemethod (MS), boundary layer method (BL) and WKB approx-imation are sophisticated and quite daunting to use because oftheir complexity and limitations [13–15]. On the contrary, therenormalization group method does naive perturbation expan-sion and requires little prior knowledge and hence is very con-venient to apply in practice. In Ref. [16], the renormalizationgroup analysis was further extended to the treatment of low-dimensional structures embedded in high-dimensional phasespace. In the current paper, we will use renormalization grouptechnique to carry out a perturbation analysis of synchroniza-tion process in the well-known Kuramoto model of coupledoscillators based on the Koopman operator.Synchronization exists ubiquitously and is intensely stud-ied, especially in coupled nonlinear systems, such as rhythmicapplause of concert audience, synchronous vibration of pace-maker cells, fluid flow patterns, outbreak of infectious dis-eases, etc. Kuramoto model is a prototyped model for the a r X i v : . [ n li n . C D ] S e p study of synchronization and has become a paradigm [17–22], though plenty of mysteries still exist concerning its dy-namical behaviour. Thus, it constitutes a good test bed of theKoopman analysis for high-d nonlinear systems. Recently, Ottand Antonsen proposed a nice technique referred to as the QAmethod [23, 24], which rigorously defines a low-dimensionalsub-manifold in the limit of large oscillator population thatcontains the synchronized state. Its success in di ff rent con-texts [25, 26] indicates possible low-dimensional behaviour inthe synchronization of networked systems though the dynam-ics could be weakly chaotic. Therefore, it seems possible topin down possible synchronization transitions with single in-dex defined in the Koopman analysis of the Kuramoto modeleven when the number of osciallators is finite.Although the renormalization group approach is good foroscillator systems with weak coupling, due to the limitation inthe expansion order or convergence, especially when the sys-tem size is large, it is hard to get the exact equation of motionfor a real system with moderate or strong coupling. A lot ofprogress has been made in numerical analysis of the synchro-nization of Kuramoto model on complex networks, arousinga wave of research [20, 21, 27, 28]. In Ref. [20], global tran-sition was found by checking the complex order parameter.But this method can not accurately get the critical couplingstrength, not to mention the local synchronization point. InRef. [21], the local and global phase transition points can befound by plotting the frequency tree of the oscillation, but westill need to know the evolution of all oscillators. In this pa-per, we will use Koopman operator to analyze the Kuramotomodel, based on the Hankel matrix and singular value decom-position (SVD) for dynamic mode decomposition (DMD). Inview of the time delay embedding theory, we reconstruct lo-cal e ff ective dynamical modes from a few oscillators. Fromthe approximation of the Koopman operator, the frequencystructure is identified and the important frequency componentis obtained. At the same time, the phase transition point canbe located by considering the correlation function of the dom-inant eigenfunctions at di ff erent coupling strengths. The reli-ability of the Koopman method will be verified in comparisonwith a direct numerical computation.This paper is organized as follows: In Section II, we will in-troduce Koopman operator, including its definition, eigenval-ues and eigenfunctions of spectral decomposition, and com-monly used algorithms for its computation. In Section III,we take the nearest-neighbor Kuramoto model as an exam-ple to carry out a perturbation analysis by a renormalizationgroup method, which yields an approximate analytical solu-tion. Then, the Koopman operator is applied to analyze theevolution data when the coupling is strong, in which eigen-values and eigenfunctions are used to unfold the frequencystructure of the system and locate the phase transition point.In section IV, the method is extended to the Kuramoto modelon a complex network. The last section is the conclusion ofthe whole paper. II. KOOPMAN OPERATOR
Typical orbits are meandering around in a complex way inthe phase space of a nonlinear system when the dynamics ischaotic, which poses considerable challenge to its description.Nevertheless, Koopman realized that the evolution of state-space observables can be described by a linear operator, whichis later called the Koopman operator [2, 3].The overall dynam-ics could be decomposed into simpler modes characterized byits eigenvalues and eigenfunctions, which laid the foundationof a statistical approach to nonlinear evolution. Comparedwith the classical approach to individual state space trajec-tories, this method is more suitable for global analysis andoptimal control. In this section, we will sketch the main in-gredient of the theory and several numerical schemes for itsapplication.
A. Definition
Consider a continuous-time dynamical system in the n -dimensional phase space M ,˙ x = F ( x ) , (1)where x ∈ M is the state vector and F : M → M gives thevector field which is usually nonlinear. S t ( x ) denotes the so-lution to Eq. (1) with the initial value x , which defines adi ff erentiable flow in the phase space. For fixed t , S t ( x ) ac-tually gives a di ff eomorphism. The Koopman operator K isa linear operator, which acts on di ff erentiable functions in thephase space of a dynamical system and gives the evolution ofthe function along orbits.More explicitly, based on a di ff erentiable function g ( x ) inthe phase space M , we may define its evolution g ( t , x ) ≡ g ( S t ( x )) . (2)All the observable functions constitute an infinite-dimensionallinear vector space. A family of Koopman operators K t ( t ∈ [0 , ∞ )) may be defined in this linear space K t g ( x ) = g ( S t ( x )) . (3)For a fixed t , K t maps the observable g ( x ) to g ( t , x ) , whichis actually more naturally defined in a discrete dynamical sys-tem. For a map y = T ( x ) , the Koopman operator K acts on g ( x ) as follows K g ( x ) = ˜ g ( x ) = g ( T ( x )) , (4)which maps a function value to the next one on an orbit. B. Eigenvalues and eigenfunctions
From the above definition, it is easy to see that Koopmanoperator is closely related to the dynamics of the system. Infact, all the orbit information could be excavated from its ac-tion on a complete set of observables. Alternatively, a linearoperator could also be characterized by its eigenfunctions andeigenvalues which may be more stable under perturbationsand thus good for numerical computation. In Ref. [2], Koop-man himself related the spectral properties of the Koopmanoperator to conservativity, integrability and ergodicity. Later,Mezi´c used eigenfunctions of Koopman operator to extract in-variant sets and ergodic blocks in the state space [7, 8, 29–32].Below, some spectral properties of the Koopman operator arediscussed.Let the eigenfunction ϕ k ( x ) of Koopman operator K corre-spond to the eigenvalue λ k , then K ϕ k ( x ) = ϕ k ( T ( x )) = λ k ϕ k ( x ), k = , , . . . . (5)If we have the eigenvalues λ k , λ k with the eigenfunctions ϕ k ( x ), ϕ k ( x ), the product λ k , λ k is also an eigenvalue withthe eigenfunction ψ ( x ) = ϕ k ( x ) ϕ k ( x ), since K ψ ( x ) = ψ ( F ( x )) = ϕ k ( F ( x )) ϕ k ( F ( x )) = λ k λ k ϕ k ( x ) ϕ k ( x ) . (6)Specifically, the n th ( n ∈ Z ) power of λ k is also an eigen-value, corresponding to the eigenfunction ϕ nk ( x ) = ( ϕ k ( x )) n .Thus for the Koopman operator, producs of eigenfunctions areeigenfunctions, which may hinder the search for the principaleigenfunctions which look relatively simple but dominate theevolution in some sense.The Koopman operator K can be uniquely decomposed intoa regular and a singular part corresponding to the continuousand the discrete spectra [8]. If the asymptotic dynamics is sim-ple (fixed points, periodic orbits, invariant tori), the spectrumof the Koopman operator would have enough eigenvalues sothat the eigenfunctions could be used as a basis to expand ob-servables g ( x ) as follows [5, 7] g ( x ) = (cid:88) k b k ϕ k ( x ) , (7)where b k is the expansion coe ffi cient. For a finite number ofthe basis functions, the expansion is approximate.On an ergodic component, the Koopman operator is uni-tary so its spectrum concentrates on a unit circle in the com-plex plane [8], while other eigenvalues mark the growth ordecay of the corresponding modes and thus signal transientstates. For λ k = ϕ k ( x ), we have ϕ k ( x p ) = K ϕ k ( x p − ) = ϕ k ( x p − ) = . . . = ϕ k ( x ) ( p ∈ Z repre-sents the discrete time). Along an orbit, the eigenvalue doesnot change with time and all the points on the orbit assumethe same value, which maps out an invariant subset of the dy-namical system [29, 30]. For a long chaotic orbit, the subsectshould be the ergodic component that contains the orbit. For λ k = e I θ ( I = √− θ is real, ifthe corresponding eigenfunction ϕ k ( x ) robustly exists, then ϕ k ( x p ) = K ϕ k ( x p − ) = e I θ ϕ k ( x p − ) = . . . = e Ip θ ϕ k ( x ) , (8)which depicts periodic motion of certain type [8]. C. Dynamic Mode Decomposition
With the development of computer and the great progressin numerical calculation, to obtain, process and analyze a large amount of dynamics data has become possible, andthis data-driven approach enables new techniques and newparadigms [33]. The traditional point by point description ofthe system focuses on a single orbit, while the operator theorydescribes all orbits as a whole which thus naturally providesa global view of system dynamics. DMD is closely relatedto the spectral analysis of the Koopman operator [34, 35]. Infact, even at the early times of its development [8], spectraldecomposition is proposed for model reduction and mode de-compositions based on data, either from experiments or sim-ulation. DMD can be regarded as a numerical approximationof the Koopman spectral analysis. Below, di ff erent numericalschemes for the DMD will be introduced based on the Koop-man operator.Most often, a time series is available: x , x τ , x τ , . . . , x ( n − τ , x n τ , which is used to extract key dynamical modes or informationabout the system under investigation. For this purpose, a dic-tionary of observables g i ( x ) , i = , , . . . , m are selected to ap-proximate the functional space, which leads to the followingdata matrix: X = ( g ( x ) , g ( x ) , . . . , g m ( x )) = g ( x ) g ( x ) . . . g m ( x ) g ( x τ ) g ( x τ ) . . . g m ( x τ ) ... ... . . . ... g ( x ( n − τ ) g ( x ( n − τ ) . . . g m ( x ( n − τ ) , Y = (˜ g ( x ) , ˜ g ( x ) , . . . , ˜ g m ( x )) = g ( x τ ) g ( x τ ) . . . g m ( x τ ) g ( x τ ) g ( x τ ) . . . g m ( x τ ) ... ... . . . ... g ( x n τ ) g ( x n τ ) . . . g m ( x n τ ) , (9)where each column of the data matrix X , Y is an n -point rep-resentation of the observables and the points are picked upuniformly in time from a long orbit. The number n of repre-sentative points should be large enough to give good spatialresolution. The projection K τ of the Koopman operator inthis finite-dimensional function space, therefore, could be ob-tained from the equation below: K τ X = Y , or K τ = Y X T ( XX T ) − , (10)by a matrix inversion (or pseudo-inversion). The selec-tion of basis functions is very important, which depends onthe features to be explored. We need enough basis func-tions to describe these features and reach reliable conclu-sions [34, 36, 37]. For a good basis, m could be small and thecomputation load is thus much reduced, especially for high-dimensional systems. But an optimal selection of the trunca-tion and the form of basis functions is a di ffi cult problem.Commonly used basis functions include Gaussian basis,Fourier basis, polynomial basis and so on, which may be goodon certain occasions but may not be so on others. Fortunately,it is possible to numerically extract a natural basis based on theHankel matrix, which directly uses the delay embedding of themeasurement values, i.e. g ( x ) = x . Delay embedding is a ge-ometric reconstruction scheme for attractors in nonlinear sys-tems and extremely useful in the modern theory of time seriesanalysis [38, 39]. In Ref. [5], it is proved that the eigenfunc-tions and eigenvalues of Koopman operator can be obtainedby applying DMD to the Hankel data matrix with limited ob-servations. By combining delay embedding with DMD, weare able to extract interesting dynamics information about theinaccessible state space from the partial data. More explicitly,we write the Hankel matrix as follows: H = x x τ . . . x m τ x τ x τ . . . x ( m + τ ... ... . . . ... x ( n − τ x n τ . . . x ( n + m − τ , (11)where each column could be viewed as an observable basedon the n points along a typical trajectory in the phase space.As is well known, a non-trivial smooth function defined in thephase space will soon become rather rugged upon the actionof the Koopman operator if the system is chaotic. Therefore,the high wave-number undulations could not be captured inthis representation and the best that we can hope for is a gooddescription of features (or averages of features) which couldbe depicted with the n points. As the column vectors of theHankel matrix are usually not even close to orthogonal, it is agood practice to first apply the SVD to the Hankel matrix toextract important directions forming an orthogonal basis [1]. As a result, the large-scale structures are contained in thefirst, say, r columns. We follow this procedure in the analysisbelow by writing the Hankel matrix as H = U Σ V T = | | | u u . . . u m | | | σ . . . σ . . . ... ... . . . ... . . . − v T −− v T − ... − v Tn − , (12)where as a common practice the positive singular values arearranged in a decreasing order, reflecting the relative impor-tance of the corresponding columns of U and V . The columnsof U could be viewed as a new orthogonal basis while eachcolumn of V T is the coordinates of the observable in this newbasis divided by the corresponding singular value. In anotherword, the original Hankel matrix becomes Σ V T in the newframe. The selection of the truncation order r depends on therelative magnitude of the singular values as well as the fea-tures under investigation. For safe, here we keep all the modeswith σ i ≥ − and thus ´ H = ´ U ´Σ ´ V T = | | | u u . . . u r | | | σ . . . σ . . . ... ... . . . ... . . . σ r − v T −− v T − ... − v Tr − . (13) A Koopman analysis based on ´Σ ´ V could be performed sim-ilar to what has been done in Eq. (10) . The eigenfunctionscould be represented in the original function space by a ma-trix multiplication. As non-essential parts are ignored, thecomputational complexity is greatly reduced. As a result, weachieve both e ffi ciency and stability in the computation. TheSVD is essentially a linear transformation, which is widelyused in numerical analysis for dimension reduction or noisefiltering [1, 40]. In the current context, the Koopman analysiswill further reduce the dimension by embedding the temporalevolution in the analysis while maintaining the robustness bykeeping only the dominant eigenfunctions. III. APPLICATION TO A KURAMOTO MODEL WITHNEAREST-NEIGHBOR COUPLING
Synchronization is closely related to our life, so it has al-ways been a hot research topic [41] . In 1975 , Kuramotoproposed his model [42] based on Winfree’s research [43] ,which is classic for the description of synchronous behaviorof a large number of coupled oscillators. It was originallyused in biological rhythm, chemical oscillations, and later ex-tended to other applications [19–22] , which describes limit-cycle oscillators with distinct natural frequencies and sinu-soidal couplings, where amplitude information is ignored andonly phase information is considered. The model is simpleenough for various analytical derivation, rich enough to dis-play a variety of synchronization modes, and flexible enoughto adapt to a variety of physical scenarios.Without lack of generality, we consider the classical andsimplest Kuramoto model with a general coupling structure:˙ θ i = ω i + K N (cid:80) j = a i j sin( θ j − θ i ), i = , · · · , N , (14)where N is the total number of oscillators. θ i is the phase ofthe i -th oscillator and ω i is the natural frequency of the i -thoscillator drawn from a common distribution ρ ( ω ). The pa-rameter K is the coupling strength (also known as couplingcoe ffi cient) between oscillators. a i j = a i j = ff erent networks. A. Kuramoto model on a circle
In this section, we consider the nearest-neighbor Kuramotomodel [44–46], which includes only nearest-neighbor cou-plings on a circle:˙ θ i = ω i + K (sin( θ i + − θ i ) + sin( θ i − − θ i )), i = , · · · , N . (15)Here, the circle topology indicates a periodic boundary condi-tion: θ = θ N , θ N + k = θ k .It can be seen from Eq. (15) that an oscillator only inter-acts with its two neighbours in an attractive manner. Whenthe coupling is weak or the natural frequency di ff erence islarge, the oscillators tend to rotate independently and the over-all dynamics is incoherent. With increasing coupling strength,the average frequency of neighboring oscillators tends to con-dense. When a threshold is reached, a small group of oscil-lators fall into synchrony and more or less acts like one ”big”oscillator. For large enough coupling strength, all oscillatorsget synchronized and rotate together with a frequency whichis the algebraic average of the natural frequencies of these N oscillators. To characterize the long-time rotation behaviour,an average frequency is defined for each oscillator i :ˆ ω i = (cid:104) ˙ θ i (cid:105) = lim T →∞ T (cid:90) T ˙ θ i ( t ) dt , (16)which is a relatively simple index and amenable to numericalcomputation but may be used to select important modes in theKoopman analysis as will be shown below.To facilitate the analytic computation, the change of thevariable z i = e I θ i ( I = √− z i = I ω i z i + K ( z i + − z i ¯ z i + + z i − − z i ¯ z i − ), i = , · · · , N , (17)where the bar over a symbol denotes complex conjugation.The boundary condition then becomes z = z N , z N + k = z k .In the following, without loss of generality, Koopman anal-ysis will be carried out for Eq. (17) with N = − − . − . . . B. Renormalization group analysis
After Goldenfeld et al extended the RG analysis to the treat-ment of di ff erential dynamical systems, great progress hasbeen made in both applications and theoretical understand-ing. In fact, many singular perturbation methods can be under-stood in terms of a renormalization process. The amplitude orphase equations could be regarded as RG equations and thusderivable in the framework of the RG theory [47]. A geo-metric interpretation is proposed based on the classical theoryof envelopes in di ff erential geometry [48–55], which leads toapplications in the computation of central manifold [56], het-eroclinic orbit in nonlinear systems [16], and so on. Here,we apply the RG technique to a perturbation analysis of thenearest-neighbor Kuramoto model (17).The RG analysis starts from a naive series expansion of thesolution to the di ff erential equation. In fact, an oscillator ro-tates with an average frequency after transient, rather than itsnatural frequency. Hence, we introduce a new parameter ˜ ω i indicating the average frequency of the i -th oscillator and getthe transformed formula:˙ z i = I ˜ ω i z i + ε K z i + − z i ¯ z i + + z i − − z i ¯ z i − ) − I ε ( (cid:52) ω i ) z i , (18)where (cid:52) ω i = ˜ ω i − ω i , being an unknown frequency correctionterm of the i -th oscillator, which may be used to eliminatepossible resonance terms.Next, we make an expansion in ε : z i = z i + ε z i + ε z i + · · · , (cid:52) ω i = (cid:52) ω i + ε (cid:52) ω i + ε (cid:52) ω i + · · · , (19)where the first subscript i in z ik and (cid:52) ω ik represents the indexof the oscillator, and the second subscript k indicates the ex-pansion order. Put the perturbation expansion into Eq. (18)and compare di ff erent powers of ε , we get a set of di ff erentialequations: ε (0) :˙ z i = I ˜ ω i z i ε (1) :˙ z i = I ˜ ω i z i + K z ( i + − z i ¯ z ( i + + z ( i − − z i ¯ z ( i − ) ε (2) :˙ z i = I ˜ ω i z i + K z ( i + − z i ¯ z ( i + − z i z i ¯ z ( i + + z ( i − − z i ¯ z ( i − − z i z i ¯ z ( i − ) − I (cid:52) ω i z i · · · (20)The first equation is a linear di ff erential equation, which caneasily be solved. z i = A i e I ˜ ω i ( t − t ) , (21) where A i is the initial value, being a function of t , to be usedas the renormalization variable later.The higher-order equations can be solved based on low-order results. For example, at the first order, to eliminate theresonance term, we set (cid:52) ω i = K ˜ ω i − − ˜ ω i + K ˜ ω i + − ˜ ω i leading to: z i = K I ˜ ω i − ˜ ω i − ( z ( i − + ¯ z ( i − z i ) + K I ˜ ω i − ˜ ω i + ( z ( i + + z i ¯ z ( i + ) . (22)Because of the limited space, higher order results will not be listed. The approximate solution is:´ z i = z i + ε K (cid:34) I ˜ ω i − ˜ ω i − ( z ( i − + ¯ z ( i − z i ) + I ˜ ω i − ˜ ω i + ( z ( i + + z i ¯ z ( i + ) (cid:35) + · · · , (23)where ´ z i = ´ z i [ t ; t , A ( t ) ] is a function of the initial time t and the initial vector A ( t ) . The higher the order of ε , the more accurate is the approximate analytic solution.We can get the polynomial expansion of A i on z i fromEq. (23) after setting t = t : A i = c i z i − + (1 + c i + ( − c i + ) z i + − c i + z i + + c i z i − z i + − c i + z i ¯ z i + − c i z i − ¯ z i − ( − c i + ) z i z i + + . . . , (24)where c i = K ε I ( ˜ ω i − ˜ ω i − ) .The RG principle tells that the physical quantity should beindependent on the starting point of the renormalization [47].Adopting the idea from Ref. [16], we focus on the dynamicsinduced by the coordinate z i on some submanifold, along theeigen-direction in the complex domain near the origin. Theinitial vector is A = (0 , . . . , A i , . . . , A i component. According tothe renormalization equation d ´ z i [ t ; t , A i ( t )] dt (cid:12)(cid:12)(cid:12)(cid:12) t = t = , (25)the equation for dA i ( t ) / dt could be obtained. To removenon-autonomous terms, we take a limit t = t . The resultingevolution equation for A i ( t ) is˙ A i = i ˜ ω i A i , (26)which looks simple and its solution is A i ( t ) = A i , e i ˜ ω i t ,where A i , is the initial value. Based on Eq. (3) in Sec. II A andEq. (5) in Sec. II B , K τ A i , = A i ( τ ) = A i , e i ˜ ω i τ , which indi-cates that e i ˜ ω i τ and A i , (or A i ) is an eigenvalue and eigenfunc-tion of the Koopman operator K τ respectively. Here, we keepthe simple oscillation dynamics by eliminating the resonanceterms order by order, and it turns out that the renormalizedfrequency ˜ ω i matches very well with the average frequencyˆ ω i computed from the time series. C. Koopman analysis
The above analytic approach is possible only for lowestorders and converging well for weak couplings. In order tostudy more general cases, numerical computation has to chipin. In this section, we will analyze the frequency structureof the coupled system, then extract dominant frequencies andlocate phase transition points according to the eigenvalues ofthe Koopman operator. TABLE I: Comparison of coe ffi cients of the principaleigenfunction obtained by the RG and the Koopmancomputation for the Kuramoto model on a circle of N = K = .
08 . The last line is the ratio of the two. z z z ¯ z ¯ z z RG 1.00000 0.14588 0.14588 0.01064Koopman 43.98444 5.84375 5.79867 0.42083RG / Koopman 0.02274 0.02496 0.02516 0.02529
We use the observable(s) Z t to construct the Hankel matrix: H = Z Z τ . . . Z m τ Z τ Z τ . . . Z ( m + τ ... ... . . . ... Z ( n − τ Z n τ . . . Z ( n + m − τ , (27)where Z t can be a row or a column vector composed of z ti = z i ( t ) , and the superscript t is the time indicator. According tothe algorithm in Sec. II C, with the evolutionary step τ = . n = m = Z t = ( z t ) , only the observations from os-cillator 1. In Fig. 1 , we draw the expansion coe ffi cient of theobservable( Z , Z τ , . . . , Z ( n − τ ) T = ( z , z τ , . . . , z ( n − τ ) T (28)with respect to eigenfunctions ordered by frequencies withdi ff erent coupling strengths. For ease of plotting, we limitthe abscissa to the interval [ − , (a) K = .
61 (b) K = . K = .
43 (d) K = . FIG. 1: Expansion coe ffi cients of the observable ( z , z τ , . . . , z ( n − τ ) T on the eigenfunctions of the Koopman operator for theKuramoto model on a circle of N = ff erent coupling strengths when Z t = ( z t ) .it can be observed that the components become more andmore complex with the increase of the coupling strength, andthe frequency proportion of nonlinear components increasessteadily. Of course, some local synchronizations in the pro-cess are not clearly observed from the figure. When reach-ing a critical coupling strength, all oscillators are synchronousand experiencing a global phase transition, where the mainfrequency is the average of individual generic frequencies asshown in Fig. 1(d). It seems that the system can be treatedas ”one” cluster rotating at a common frequency to some ex-tent after synchronization. To sum up, eigenfunctions describemodes with di ff erent frequencies, and the expansion coe ffi -cients on the eigenfunctions of the observable give a kind of”spectrum”. From the change of spectrum with the couplingstrength, we find the global synchronization point K = . K = .
08. The ratio of the correspond- TABLE II: The critical coupling strength for the Kuramotomodel on a circle of N = numerical 0.16 0.20 0.26 0.32 0.38 0.44 1.45Koopman 0.16 0.20 0.28 0.34 0.40 0.44 1.45 ing coe ffi cients is about 0.025, which indicates that the tworesults are proportional to each other and the two approacheswell match in case of small couplings. D. Tracking of bifurcations
As the coupling strength increases, the actively involvedeigen-modes increase first, at some point start to decrease andcollapse to a principal one beyond the critical point. How tofind the important ones and how to determine the bifurcationFIG. 2: The average frequency v.s. coupling strength K forthe Kuramoto model on a circle of N = Z t = ( z ti ) , i = , , . . . ,
6. The blue curve is the numericalsolution, and the red hollow circle is the Koopman solution.FIG. 3: Correlation coe ffi cient (red dashed line) and theaverage frequency (blue solid lines) v.s. the coupling strength K for the Kuramoto model on a circle of N = Z t = ( z ti ) , i = , , . . . , N .points along the way are key questions that we need to checkwithin the framework of the Koopman operator. Throughthe analysis above, we find that no matter what the couplingstrength is, the largest expansion coe ffi cient is always attachedto the average frequency which is thus regarded as the princi-pal mode of the oscillator.If we take the observable Z t = ( z ti ) , i = , , . . . , N , theeigenfunction expansion of the motion of each oscillator isdi ff erent and changes with coupling strength, while the dom-inant eigenfunctions always correspond to the average fre-quencies of individual oscillators. In virture of this observa-tion, we may extract six sets of spectra for di ff erent coupling strengths but with N = Z t = ( z ti , z ti + ) or a column vector Z t = ( z ti , z ti + ) T , thekoopman analysis gives a result similar to what is in Fig. 2 sowe will not show it again.From the above discussion we know that there is a corre-spondence between the eigenvalues and the average frequen-cies of the oscillators, while the eigenfunctions map out rela-tive importance of di ff erent modes in the phase space.Thus, itis possible to find out the phase transition point by checkingthe change of the eigenfunction with the coupling strength.With Z t = ( z tn ) , we may define the correlation coe ffi cient ofthe eigenfunction ϕ n of oscillator n: ρ n ( K j ) = Cov ( ϕ n ( K j − ) , ϕ n ( K j )) (cid:112) Var ( ϕ n ( K j − )) Var ( ϕ n ( K j )) , (29)where in this section, we take K j − K j − = .
01. The eigen-function ϕ n ( K j ) receives the largest expansion coe ffi cient ofthe observable at certain coupling strength K j , which can beused to predict some of the synchronization points. For bettercharacterization, we take the summation: ρ = N (cid:80) Nn = ρ n ( K j ), and draw the ” ρ - K ” graph with the red dashed line in Fig. 3. In most cases, the correlation coe ffi cient ρ is close to 1, be-cause the mode changes smoothly with the coupling strengthif there is phase transition. When a local or a global syn-chronization is about to occur, there is a rapid change of theprincipal eigenvalue and eigenfunction, resulting in a perceiv-able change of ρ . In the graph, henceforth, a dip in a smoothbackground is featured for ρ< IV. APPLICATION TO A KURAMOTO MODEL ON ACOMPLEX NETWORK
In the above section, Koopman operator is used to ana-lyze the nearest-neighbor Kuramoto model of N = N =
10 oscillators.the Koopman analysis to the Kuramoto model on a complexnetwork. First, we discuss the network topology.
A. Kuramoto model on a complex network
In 1998, Watts and Strogatz proposed that biological or so-cial networks are between completely regular and random,and established the WS model to generate this kind of net-works [57]. This model starts from a completely regular net-work and then randomly reconnects the nodes in the network.As a result, all the path lengths in the network are short and theclustering degree is high, which is consistent with the small-world characteristics [58] and carries on rich statistical anddynamics behaviors. So it is interesting to extend applicationof the Koopman method to systems defined on WS networks.For ease of observation, we generate a network topology of N =
10 oscillators in Fig. 4 , where the ten circled red dotssymbol ten oscillators, black numbers marking indices of os-cillators, and black edges indicating the interaction topology.The edge connection between vertex i and j gives a i j = B. Koopman analysis
With the network topology in Fig. 4, here we focus on twotypical distributions of frequencies: randomly distributed inthe interval [ − ,
1] with uniform probability (called random-uniform distribution, RUD) and chosen linearly within [ − , ρ - K ” graph for this Kuramoto model with reddashed line in Fig. 7 and the natural frequency distributionis LD in Fig. 7(a) and RUD in Fig. 7(b) respectively. Thedips in the correlation profile correspond well to phase tran-sition points in the frequency tree. Tab. III show the criticalcoupling strength obtained with direct numerical computationand the frequency jump points identified with the Koopmananalysis, where K j − K j − = . K = .
009 is not detected by the Koopmantechnique in Tab. IV with K j − K j − = . K = .
060 covers twodips in the correlation coe ffi cient at K = . , .
063 , whichmay indicate possible occurrence of very subtle bifurcation.The marching stepsize of the coupling strength is smaller thanin Sec. III C, resulting in a computation error less than 0.01.In conclusion, the results of the Koopman analysis remainsreliable for the Kuramoto model on a complex network.
V. CONCLUSIONS
With the help of spectral properties of Koopman operator,we are able to divide the phase space, extract important dy-namics patterns, and determine global properties of couplednonlinear systems. In this paper, this technique has been ex-tended to the analysis of Kuramoto oscillators on a ring ora complex network. After a Hankel matrix is built from thetime series and is preprocessed with SVD, we carry out theDMD based on an approximation of the Koopman operator.Eigenfunctions and eigenvalues are obtained to unfold thefrequency structure, and dominant modes are determined bychecking the weight of di ff erent eigenfunctions. With an RGanalysis valid in case of small coupling, the eigen-frequenciesturn out to be the average frequencies of the oscillators, whichholds true even for moderate or strong couplings justified bythe Koopman approach. The correlation of eigenfunctionsat neighbouring coupling strength is a manifestation of thechange rate of the important dynamical modes, through whichlocal or global synchronization points can be found with smallerror bars. These new techniques may provide extra tools forthe extraction of important dynamical information in complexsystems.In the above analysis, near the bifurcation point, the fre-quency derived from the Koopman analysis sometimes does0 (a) LD (b) RUD FIG. 5: Two frequency distributions for the Kuramoto model on a network. (a) a linear distribution; (b) a realization from theuniform distribution on [ − , (a) LD (b) RUD FIG. 6: The average frequency vs the coupling strength K for the Kuramoto model on a network shown in Fig. 4, for di ff erentgeneric frequency distributions: (a) linear; (b) random as displayed in Fig. 5 , and for di ff erent computations: the directnumerical average (blue solid lines) and the Koopman analysis (red hollow circles).TABLE III: The critical coupling strength obtained by the direct numerical and the Koopman computation for the Kuramotomodel on a complex nextwork shown in Fig. 4 with the natural frequency linearly distributed in the interval [ − ,
1] . numerical 0.088 0.100 0.116 0.140 0.190 0.254 0.384 0.388Koopman 0.086 0.098 0.116 0.134 0.190 0.246 0.380 0.392
TABLE IV: The critical coupling strength obtained by the direct numerical and the Koopman computation for the Kuramotomodel on a complex nextwork shown in Fig. 4 with the natural frequency from the uniform distribution in the interval [ − ,
1] . numerical 0.009 0.045 0.060 0.072 0.087 0.090 0.222 0.267 0.387 0.414Koopman - 0.045 0.054,0.063 0.072 0.084 0.090 0.222 0.267 0.387 0.417 (a) LD (b) RUD FIG. 7: The correlation coe ffi cient (red dashed line) of the eigenfunction and the average frequency (blue solid lines) vs thecoupling strength K . All the parameters are the same as in Fig. 6 .not match well with the average frequency. How to improvethe computation and pin down important frequencies shouldbe better probed. We use the correlation function as an in-dicator for the bifurcation events, which seems valid in mostcases but may not be the best ones since a majority of othereigenfunctions are ignored. Some important features of dy-namics may be lost in the current treatment and so it wouldbe very helpful if one can find a way to pick up contributionsfrom them. Also, in the current investigation, enough data isprovided for our purpose. If only a limited amount of datais available, what kind of information can be extracted fromthe Koopman analysis is an interesting question, which maydepend on the functional basis we use to construct the approx-imate Koopman operator.Looking back on the whole analysis, the Koopman tech-nique does not require governing equations of motion, and atime series su ffi ces for all the above operations. To reducethe complexity associated with high dimensionality and easethe di ffi culty of basis selection, we combined the concept of DMD with the time delay embedding and the SVD. Underthe premise of ensuring stability and reliability, the Koopmananalysis helps us reconstruct global information from partialobservation and at the same time minimizes the impact by un-avoidable noise. The importance of the average frequency un-expectedly announces again in this type of analysis of com-plex systems although rhythms in nature appear everywhereand have been encapsulating us already for many years up tothis very day. The Koopman analysis seems to provide a newframe for the analysis and understanding of all these fascinat-ing nonlinear systems. ACKNOWLEDGMENTS
This work was supported by the National Natural ScienceFoundation of China under Grants No.11775035 and also bythe Fundamental Research Funds for the Central Universitieswith Contract No.2019XD-A10. [1] S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, and J. N.Kutz, Nat. Commun. , 1 (2017).[2] B. O. Koopman, Proc. Natl. Acad. Sci. U.S.A. , 315 (1931).[3] B. Koopman and J. v. Neumann, Proc. Natl. Acad. Sci. U.S.A. , 255 (1932).[4] S. L. Brunton, B. W. Brunton, J. L. Proctor, and J. N. Kutz,PLoS One (2016).[5] H. Arbabi and I. Mezic, SIAM J. Appl. Dyn. Syst. , 2096(2017).[6] I. Mezi´c and A. Banaszuk, Phys. D , 101 (2004).[7] A. Mauroy and I. Mezi´c, IEEE Trans. Autom. Control , 3356(2016).[8] I. Mezi´c, Nonlinear Dyn. , 309 (2005).[9] N. Goldenfeld, Lectures on phase transitions and the renormal-ization group (CRC Press, 2018). [10] J. Zinn-Justin,
Quantum field theory and critical phenomena (Clarendon Press, 1996).[11] L.-Y. Chen, N. Goldenfeld, Y. Oono, and G. Paquette, PhysicaA , 111 (1994).[12] G. Paquette, L.-Y. Chen, N. Goldenfeld, and Y. Oono, Phys.Rev. Lett. , 76 (1994).[13] A. F. Vakakis and M. Azeez, Nonlinear Dyn. , 245 (1998).[14] J. Kevorkian, J. Cole, and A. H. Nayfeh, Bulletin of the Amer-ican Mathematical Society , 414 (1982).[15] C. M. Bender and S. A. Orszag, Advanced mathematical meth-ods for scientists and engineers I: Asymptotic methods and per-turbation theory (Springer Science & Business Media, 2013).[16] Y. Lan, Phys. Rev. E , 012914 (2013).[17] Z. Zheng, G. Hu, and B. Hu, Phys. Rev. Lett. , 5318 (2012). [18] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, Phys. Rev.Lett. , 1804 (1996).[19] J. A. Acebr´on, L. L. Bonilla, C. J. P. Vicente, F. Ritort, andR. Spigler, Rev. Mod. Phys. , 137 (2005).[20] H. Chiba, G. S. Medvedev, and M. S. Mizuhara, Chaos: An In-terdisciplinary Journal of Nonlinear Science , 073109 (2018).[21] Y. Wu, J. Xiao, G. Hu, and M. Zhan, EPL (Europhysics Letters) , 40005 (2012).[22] S. H. Strogatz, Phys. D , 1 (2000).[23] E. Ott and T. M. Antonsen, Chaos: An Interdisciplinary Journalof Nonlinear Science , 037113 (2008).[24] E. Ott and T. M. Antonsen, Chaos: An interdisciplinary journalof nonlinear science , 023117 (2009).[25] E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P. So, andT. M. Antonsen, Phys. Rev. E , 026204 (2009).[26] N.-P. Wu, H.-Y. Cheng, Q.-L. Dai, and H.-H. Li, Chin. Phys.Lett. , 070501 (2016).[27] C. A. Moreira and M. A. de Aguiar, Physica A , 487 (2019).[28] J. Zhang and J. Zhu, Automatica , 122 (2019).[29] M. Budiˇsi´c and I. Mezi´c, Phys. D , 1255 (2012).[30] I. Mezi´c and S. Wiggins, Chaos: An Interdisciplinary Journalof Nonlinear Science , 213 (1999).[31] M. Budiˇsi´c and I. Mezi´c, in Proceedings of the 48h IEEE Con-ference on Decision and Control (CDC) held jointly with 200928th Chinese Control Conference (IEEE, 2009) pp. 3162–3168.[32] I. Mezi´c, in (IEEE, 2015) pp. 7034–7041.[33] E. Kaiser, J. N. Kutz, and S. L. Brunton, Proceedings of theRoyal Society A Mathematical Physical, and Engineering Sci-ences (2017).[34] C. W. Rowley, I. Mezi´c, S. Bagheri, P. Schlatter, and D. S.Henningson, J. Fluid Mech. , 115 (2009).[35] Mezi and Igor, Annu. Rev. Fluid Mech. , 357 (2013). [36] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz,J. Neurosci. Methods , 1 (2016).[37] P. J. Schmid, J. Fluid Mech. , 5 (2010).[38] T. Sauer, J. A. Yorke, and M. Casdagli, J. Stat. Phys. , 579(1991).[39] F. Takens, in Dynamical systems and turbulence, Warwick 1980 (Springer, 1981) pp. 366–381.[40] C. W. Rowley, Int. J. Bifurcation Chaos , 997 (2005).[41] A. Pikovsky, J. Kurths, M. Rosenblum, and J. Kurths, Syn-chronization: a universal concept in nonlinear sciences , Vol. 12(Cambridge university press, 2003).[42] A. T. Winfree, J. Theor. Biol. , 15 (1967).[43] Y. Kuramoto, Lect. Notes Phys. , 420 (1975).[44] Z. Zheng, G. Hu, and B. Hu, Phys. Rev. Lett. , 5318 (1998).[45] O. Kogan, J. L. Rogers, M. Cross, and G. Refael, Phys. Rev. E , 036206 (2009).[46] H. F. El-Nashar, P. Muruganandam, F. F. Ferreira, and H. A.Cerdeira, Chaos: An Interdisciplinary Journal of Nonlinear Sci-ence , 013103 (2009).[47] L. Y. Chen, N. Goldenfeld, and Y. Oono, Phys. Rev. Lett. ,1311 (1994).[48] T. Kunihiro, Progress of theoretical physics , 503 (1995).[49] T. Kunihiro, Jpn. J. Ind. Appl. Math. , 51 (1997).[50] S.-I. Ei, K. Fujii, and T. Kunihiro, Ann. Phys. , 236 (2000).[51] Y. Hatta and T. Kunihiro, Ann. Phys. , 24 (2002).[52] T. Kunihiro, Progress of Theoretical Physics Supplement ,459 (1998).[53] T. Kunihiro, Phys. Rev. D , R2035 (1998).[54] T. Kunihiro and J. Matsukidaira, Phys. Rev. E , 4817 (1998).[55] T. Kunihiro and K. Tsumura, J. Phys. A: Math. Gen. , 8089(2006).[56] Chiba and Hayato, J. Math. Phys. , 102703 (2008).[57] D. J. Watts and S. H. Strogatz, Nature , 440 (1998).[58] G. Marconi, Nobel Lecture11