Copyright is owned by the Author of the thesis. Permission is given for a copy to be downloaded by an individual for the purpose of research and private study only. The thesis may not be reproduced elsewhere without the permission of the Author. An equation-free approach for heterogeneous networks Sidra Zafar School of Mathematical and Computational Sciences Massey University This dissertation is submitted for the degree of Doctor of Philosophy September 2023 I would like to dedicate this thesis to my sweet little girl, Bareerah Fatima. Acknowledgements First and foremost, I want to thank Almighty Allah for providing me with the strength and determination to accomplish my research. A very special thanks to my co-supervisor, Dr. Alona Ben Tal, for her selfless and consistent dedication to my intellectual and professional development, for her suggestions and the constant encouragement she gave me during my time at Massey University, and for believing in me even (and especially) when I didn’t seem to believe in myself. I would also like to thank my main supervisor Distinguished Professor Robert Mclach- lan and my other co-supervisor, Dr. David Simpson for their continuous support, expert feedback on this thesis and for giving me the flexibility to work on my own way. I am grateful to the School of Mathematical and Computational Sciences, Massey University, New Zealand for providing me with a scholarship during my PhD. The generous support has greatly helped me to advance my work. I appreciate Dr. Carlo Laing’s support during the first two years of my PhD. I would also like to thank my parents, especially my mother, who taught me the importance of education and encouraged me every step of the way to pursue my dream of getting higher education from overseas. I am very thankful to my my mother-in-law and father-in-law, for taking care of my daughter while I was away from her during my PhD. Last, but not least, I would like to thanks my children, Bareerah Fatima and Muham- mad Affan Raees, for their endless love and joy, as well as my husband, Muhammad Saad, whose dedication, affection and confidence in me, has lifted a burden off my shoulder. Together, we have accomplished it! Abstract Oscillators exist in many biological, chemical and physical systems. Often when oscilla- tors with different periods of oscillations are coupled, they synchronize and oscillate with the same period. Examples include groups of synchronously flashing fire-flies or chirping crick- ets. There are two questions of interest in this work. (1) Under what conditions will a system of coupled oscillators synchronize? and (2) Can a large system of synchronized oscillators be represented by a smaller number of variables? We study these questions for Kuramoto like models which are coupled in different ways. Examples include spatially extended and all-all coupling. We study conditions under which synchronization occurs in small and large networks by varying the coupling strength, calculating stabilities of synchronized solutions and creating bifurcation diagrams of the steady state solution as a function of coupling strength in one and two-dimensions. We use an equation free approach to approximate the coarse scale behaviour of a large coupled network, for which the equations are known, by a low dimensional description of variables, for which no governing equations are available in closed form. Our results show that a small number of variables can reproduce the behaviour of the stable solutions in the full systems. However, the equation-free approach did not work as well for the unstable solutions. Possible reasons for this are explored in the thesis. Contents List of Figures xi List of Tables xv 1 Introduction 1 1.1 Kuramoto model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Spatially extended Kuramoto network . . . . . . . . . . . . . . . . . . . . 2 1.3 Equation-free approach and its application to the Kuramoto model . . . . . 4 1.4 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Preliminary Analysis of the Kuramoto Model 7 2.1 From spatially extended network to all-all Kuramoto model . . . . . . . . . 7 2.2 Calculating the stationary solutions . . . . . . . . . . . . . . . . . . . . . . 9 2.2.1 Theorem : ""Synchronization in spatially extended and all-all Ku- ramoto networks: The equalization of angular velocities" . . . . . . 11 2.2.2 Computing the fixed points and their stabilities . . . . . . . . . . . 14 2.3 Numerical continuation of the Kuramoto models . . . . . . . . . . . . . . 16 2.3.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Variation of the coupling strength in the spatially extended network . . . . . 18 2.4.1 Variation of coupling strength in the spatially extended Kuramoto model with N = 100 and M = 5. . . . . . . . . . . . . . . . . . . . 20 2.5 Varying twisted states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.6 Varying the coupling range . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.6.1 N = 100, M = 25 . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.6.2 N = 100 and M = 50, all-all Kuramoto model . . . . . . . . . . . . 28 2.7 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 31 3 An Equation-Free Approach for the All-All Kuramoto Network 33 3.1 Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 x Contents 3.2 Methodology of applying an equation-free approach to the all-all Kuramoto network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3 Coarse-projective integration of all-all Kuramoto network . . . . . . . . . . 37 3.3.1 Varying the integration time and projection times. . . . . . . . . . . 38 3.4 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 41 4 Creating Bifurcation Diagrams Using the Equation-Free Approach in One Dimension 43 4.1 Finding the derivatives of coarse-grained coefficients using an equation-free approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2 Finding fixed points of the coarse-grained system . . . . . . . . . . . . . . 45 4.3 Bifurcation Diagrams when θi are functions of their positions . . . . . . . . 47 4.4 Bifurcation Diagrams when θi are functions of their angular frequencies . . 54 4.5 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 59 5 Creating Bifurcation Diagrams Using the Equation-Free Approach in Two Dimensions 61 5.1 Bifurcation Diagrams when θi are functions of both positions and angular frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.2 The equation-free approach with two summations for other twisted solutions 67 5.3 Varying the coupling range . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.4 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 79 6 Conclusions 81 Appendix A Appendix 87 A.1 Numerical Continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 A.2 List of Shifted Legendre Polynomials . . . . . . . . . . . . . . . . . . . . 90 A.3 Algorithm for applying Equation-free approach on heterogenous networks . 93 x List of Figures 2.1 From spatially extended to all-all network . . . . . . . . . . . . . . . . . . 8 2.2 State of oscillators for a spatially extended Kuramoto network in a stationary frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 State of oscillators with respect to time for spatially extended Kuramoto network in a rotating frame . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4 A general description of a connection within a network . . . . . . . . . . . 11 2.5 Stationary state of the spatially extended Kuramoto model using method-2 . 13 2.6 Wrapping process for solving the differential equations numerically in spa- tially extended networks. . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.7 Solutions of spatially extended Kuramoto model with N = 6 and M = 2 . . 19 2.8 Solutions of spatially extended Kuramoto model with N = 100 and M = 5 . 21 2.9 Bifurcation diagram for various twisted states as a function of K. . . . . . . 23 2.10 State of 100 oscillators for various twisted states from stable and unstable branches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.11 Solutions of spatially extended Kuramoto model with N = 100 and M = 25. 27 2.12 Solutions of spatially extended Kuramoto model with N = 100 and M = 50. 29 2.13 Bifurcation diagrams for various twisted states by varying M. . . . . . . . . 30 2.14 A bifurcation diagram created with our own software compared to the one created with CLMATCONT-L. . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1 Evolution of θi as a function of time, t and as a function of ωi. . . . . . . . 34 3.2 Snapshot of θi as a function of ωi. . . . . . . . . . . . . . . . . . . . . . . 36 3.3 Schematic diagram of the coarse-projective integration. . . . . . . . . . . . 37 3.4 The effects of variation in N1 and N2 on a coarse projective integration with three coarse grained coefficients. . . . . . . . . . . . . . . . . . . . . . . . 38 3.5 Error between the solution obtained by full integration and the solution obtained by projective integration. . . . . . . . . . . . . . . . . . . . . . . 39 xi xii List of Figures 3.6 The effects of variation in N2 on a coarse projective integration with three coarse grained coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.7 Error of the coarse-projective integration. . . . . . . . . . . . . . . . . . . 41 4.1 A flow diagram showing how to find the derivatives for the coarse-grained system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2 Schematic description of a numerical continuation technique for an unknown function g. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.3 Comparison of the solutions of twist-1 equation-free approach with the twist-1 detailed system for N = 100, M = 5. . . . . . . . . . . . . . . . . . 49 4.4 Comparison of the fitting of the stable and unstable branches for twist-1. . . 50 4.5 Comparison of the bifurcation diagrams of twist-2 obtained by the equation- free approach and the detailed system. . . . . . . . . . . . . . . . . . . . . 51 4.6 Comparison of the fitting of the stable and unstable branches for twist-2. . . 52 4.7 Bifurcation diagrams of twist-2 obtained by the equation-free approach and other twisted solutions of the detailed system. . . . . . . . . . . . . . . . . 53 4.8 Comparison of the fitting of the stable and unstable branches for twist-2. . . 54 4.9 Comparison of the solutions of twist-1 equation-free approach with the twist-1 detailed system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.10 Comparison of the fitting of the stable and unstable branches at K = 10 for twist-1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.11 Comparison of the solutions of twist-0 equation-free approach with the twist-0 detailed system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.12 Overlapping solutions of 3 and 4 coefficients for a1,a2,a3. . . . . . . . . . 58 4.13 Relationship between the state of the oscillators and the their natural frequen- cies for the all-all Kuramoto model twist-0. . . . . . . . . . . . . . . . . . 59 5.1 Comparison of the fitting of the stable and unstable branches for twist-1 with spatial function vk = (i−1)k. . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.2 Comparison of the solutions of twist-1 equation-free approach using two- summation with the twist-1 solution of the detailed system. . . . . . . . . . 64 5.3 Comparison of the two-dimensional (using two summations) fitting of the stable and unstable solutions for twist-1 with 8 combinations of coefficients. 65 5.4 Bifurcation diagrams for a11 of twist-1, N = 100, M = 5 with various coeffi- cients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.5 Comparison of the solutions of twist-2 equation-free approach using two summations with the twist-2 detailed system. . . . . . . . . . . . . . . . . 69 xii List of Figures xiii 5.6 Comparison of the solutions of twist-3 equation-free approach using two summations with the twist-3 detailed system. . . . . . . . . . . . . . . . . 70 5.7 Comparison of the solutions of twist-4 equation-free approach using two summations with the twist-4 detailed system. . . . . . . . . . . . . . . . . 71 5.8 Comparison of the fitting of the stable and unstable branches for twist-1 with the number of coefficients 3 and 4. . . . . . . . . . . . . . . . . . . . . . . 73 5.9 Comparison of the two-dimensional fitting (using two summations) of stable and unstable branches with 6 coefficients. . . . . . . . . . . . . . . . . . . 74 5.10 Comparison of the two-dimensional fitting (using two summations) of stable and unstable branches with 8 coefficients. . . . . . . . . . . . . . . . . . . 75 5.11 Comparison of the solutions of twist-1 equation-free approach using two summations with the twist-1 detailed system for N = 100,M = 10 using 6 combinations of coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.12 Comparison of the solutions of twist-1 equation-free approach using two summations with the twist-1 detailed system for N = 100, M = 10 using 8 combinations of coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.13 Bifurcation diagrams for a11 of twist-1, N = 100, M = 10 with various coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 A.1 Graphical representation for computing Eq.(A.2). . . . . . . . . . . . . . . 88 A.2 Phase portrait representation of the system described by Eq. (A.1) with respect to the variable x and parameter µ . . . . . . . . . . . . . . . . . . . 90 xiii List of Tables 2.1 Bifurcation points for various twisted states with N = 100 and M = 5 and the natural angular velocities are uniformly distributed (randomly). . . . . . 23 5.1 The absolute error differences at several points from the stable and unstable branches of bifurcation diagrams with various coefficients. . . . . . . . . . 67 5.2 The L2-norm fitting at several points from the stable and unstable branches of bifurcation diagrams. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.3 Comparison of bifurcation points for various twisted states. . . . . . . . . . 72 5.4 Error differences between true solution and approximated solution for various twisted states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.5 L2-norm of fitting at two different points for various twisted states. . . . . . 72 5.6 Error difference between the detailed system’s true solution and the equation- free’s estimated solutio for various coefficients for N = 100,M = 10. . . . . 78 5.7 L2-norm for detailed solution and the approximated solution at K = 2 and K = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 xv Chapter 1 Introduction Synchronization is a natural phenomenon that can be observed widely in nature. Examples include the rhythmic clapping of human audiences, the flashing of fireflies, the chirping of crickets, and the coordinated movements of motor neurons in the human body [7, 49, 1]. Mathematical models that capture synchronization can be traced back to 1967, when Winfree proposed a mathematical model for a network of globally coupled oscillators [62]. Inspired by Winfree’s model, a much simpler and solvable model was introduced by Kuramoto [31, 30]. Kuramoto’s remarkable model has inspired a wide spectrum of research, and various modifications have been developed and applied to a variety of fields such as neuro-science [13], power systems [14, 21], chemical engineering [30, 29], geophysics [58], and physics [56, 2, 16, 15]. This thesis contributes to the study of Kuramoto-based networks with a focus on simplifying them using the equation-free approach. We first review briefly the Kuramoto model and some of the research that grew out of it. We then review the relevant literature on the equation-free approach and its previous application to Kuramoto-based networks. 1.1 Kuramoto model Consider a population of N oscillators, each having a phase θi and an angular frequency ωi. Thus, for uncoupled oscillators, the rate of change of the phase is defined as θ̇i = ωi, (i = 1, . . . ,N). The interaction between the oscillators was included by Kuramoto as a sinusoidal function of the phase differences with a coupling strength, K, leading to the following equations: θ̇i = ωi + K N N ∑ j=1 sin(θ j−θi), i = 1, . . . ,N (1.1) 2 This model demonstrates that coupled oscillators can synchronize when the coupling strength is larger than a critical value. The original Kuramoto model was studied for the globally coupled, equally weighted oscillators. Since the development of the first globally connected Kuramoto model, numerous studies have been conducted. Several methodologies for addressing phase/frequency synchronization, partial synchronization, phase balancing, as well as the inclusion of noise and inertia in the Kuramoto model have all been addressed using a variety of analytical and numerical techniques. [1, 15, 45, 10, 11, 25, 18, 48]. In parallel with the studies on the traditional Kuramoto model, one has witnessed the emergence of an overwhelming number of works focused on effects caused by the introduction of heterogeneous connection patterns, letting the interactions to be no longer only restricted to global coupling. This generalization of the Kuramoto model in complex networks is obtained by including the connectivity in the coupling term as [2, 4]. θ̇i = ωi +K N ∑ j=1 Ai j sin(θ j−θi) i = 1, . . . ,N (1.2) where Ai j is the adjacency matrix with its elements distributed such as Ai j = 1 if there is a connection between nodes i and j otherwise Ai j = 0. After that, most of the attention has been aimed at determining how network structure influences the onset of synchronization. Initial studies combined the small world networks with the Kuramoto model [60, 24, 42, 37]. The synchronization phenomenon for a Kuramoto model has been observed on the finite-size scaling network [12, 23] and on scale-free net- works, highlighting the role played by highly connected nodes (hubs) in the synchronization phenomenon [36]. 1.2 Spatially extended Kuramoto network Studying the Kuramoto model with various complex networks that have different network architectures demonstrates the importance of coupling topologies in studying the synchronization of coupled oscillators. Hence, coupling range becomes one of the main characteristics in analyzing different kinds of complex Kuramoto networks. The cases in which oscillators interact only with their close neighbors are known as local coupling [28], whereas the case in which each oscillator interacts with all the others is known as all-all coupling or global coupling. It should be noted that in the case of an all-all network, the coupling topology loses its meaning as all the oscillators are driven by the same mean field from all the oscillators. This fact makes it simpler, and thus many analytical methods have been proposed for studying its collective properties [46, 59, 49, 3]. However, in the 3 case of non-locally coupled oscillators, many numerical techniques can be applied, and one kind of solution that has been observed is the so-called chimera state in which some oscillators synchronize while others do not [33, 47]. Another type of solution in nonlocally coupled oscillators are traveling waves. These are defined by a fixed profile and move with a constant speed. These traveling waves are studied for many neural systems [17, 27, 50, 51]. One specific type of traveling wave is the uniformly twisted wave for which the state of the oscillators changes linearly with the spatial position. One type of non-locally coupled Kuramoto model for which uniformly twisted states (details are discussed in section 2.5) appear is referred to as the spatially extended Kuramoto model. Consider oscillators θi, moving with their individual angular frequencies ωi and coupled to their M nearest neighbors on both sides [34]: θ̇i = ωi + K 2M+1 M ∑ j=−M sin(θi+ j−θi), i = 1, . . . ,N (1.3) Here, N is the number of oscillators, K is the coupling strength, and periodic boundary conditions for index i are taken. Periodic boundary conditions in the given model of oscillators mean that the system forms a ring where each oscillator is coupled to its M nearest neighbors on both sides, treating the first and last oscillators as neighbors. This circular arrangement enables the emergence of stable synchronous states and uniformly twisted traveling waves, known as q-twisted states, characterized by constant phase shifts between neighboring oscillators, resulting in spatially linear patterns. The dynamics of such a model depend on the number of oscillators N and the coupling range M. The given system Eq. (1.3) always has a stable synchronous state, θi = ωit +C. However, as the coupling range varies, new states emerge in the form of uniformly twisted traveling waves θi = ωit + 2πqi N +C, i = 1, ...,N, where 2πq N is the phase difference between the oscillators. These so-called q-twisted states (characterized by an integer q) are simple traveling waves with constant phase shifts between neighboring oscillators, resulting in spatially linear patterns. Such models have been studied previously, and analysis of the twisted states has been performed. For example, in the case of identical oscillators, [61] discovered that many twisted states exist for such types of networks, and the basins of attraction follow a Gaussian distribution in relation to a winding number q. Later, this idea has been extended to a heterogeneous unimodal frequency distribution network, and partially coherent twisted states have been observed [44]. Thus, these partially coherent twisted states are defined by the expression of a spatially dependent, complex order parameter whose argument differs from 2πq across the ring. The research conducted by [44] considered the Lorentzian frequency distribution for non-local coupling phase oscillators, and the stability of the resulting partially 4 coherent twisted states was theoretically studied using the Ott-Antonsen or (OA) approach [46]. The Ott-Antonsen ansatz provides a powerful method to simplify the analysis by projecting the dynamics of the coupled oscillator system onto a lower-dimensional manifold. Laing (2016) also discussed the same model for heterogeneous networks but included different forms of delays in their coupling. Again, the OA approach is used to derive evolution equations for a spatially-dependent order parameter. More recent work is carried out by [63] where two different frequency distributions are considered, and the stability analysis of two different types of twisted waves is done by reducing the model to a continuum limit using the OA approach. 1.3 Equation-free approach and its application to the Ku- ramoto model For decades, reducing high-dimensional models to lower dimensions has been a prominent research direction for understanding and analyzing the complicated dynamics of larger and complex networks. Many strategies for reducing different models have been used in the past, such as the perturbation technique [43], and, more recently, the data-driven reduction technique [54]. A recently developed technique, known as an equation-free approach [26, 19], was used to approximate the macroscopic behavior of a detailed system for which the equations are known, by a low-dimensional description of variables for which no governing equations are available in closed form. This can be accomplished by performing short bursts of appropriately designed computational experiments with the detailed, "fine-scale" network dynamic evolution model to obtain the numerical information of the appropriate coarse observables. Since the macroscopic variables change across a much longer time scale compared to the fine-scale variables, the separation of time scales should be considered when choosing coarse variables. The transitions between the microscopic and the macroscopic scales depend on two operators: the restriction and lifting operators. The restriction operator (many-one relation), in general, maps the detailed description of the system to the coarse- grained level, whereas the lifting operator (one-many relation) maps from the coarse-grained to the detailed system. Equation-free methods have been successfully implemented on various problems, including molecular dynamics [9], collective animal behavior dynamics [41], cell population [8], dynamical models on static networks [52], as well as several specific network models [57, 20]. However, this technique mainly depends on (a) identifying and defining a suitable set of coarse variables for which a closed, reduced description of network evolution can theoretically be obtained, and (b) the ability to switch back and forth between 5 the two levels of system description—the "fine" and the "coarse-grained" levels. Thus, one of the important tasks in coarse-graining with this approach is selecting appropriate coarse observables. Diffusion maps [55] and the polynomial chaos expansion [35] have been used for finding these coarse observables. The equation-free approach with polynomial expansion has also been used to reduce the Kuramoto model to lower dimensions [40, 41]. Moon et al. investigated an original Kuramoto model (all-all coupling) consisting of a finite population of nonidentical coupled oscillators. The coarse-grained variables were defined by realizing the correlation between the phase angle and natural frequencies, where both are considered as random variables. The phase angle is then expanded in terms of the Hermite polynomials of the heterogeneous frequencies with the expansion coefficients selected as coarse variables. Later, this work was extended where the coarse-graining approach was implemented on the Kuramoto model with some network structure (spectral graphs) [52]. More recently, Bertlan (2017) used the polynomial expansion method for coarse-graining the Kuramoto model and described the state of oscillators as a function of both heterogeneous parameters (frequencies and connectivity). 1.4 Discussion and Conclusions Synchronization has been studied in a wide range of Kuramoto-based networks. In this thesis, we investigate two types of networks, namely the all-all and the spatially-extended network, and observe the transition between them. Unlike many previous studies, our focus is on the simplification of the networks using the equation-free approach. Although some studies have used the equation-free approach, our study is unique in the sense that the equation-free approach has not been used for this specific Kuramoto model (spatially extended network) previously, which we considered. As mentioned in section 1.2, the spatially extended Kuramoto model has been studied in the past, but to represent the system in lower dimensions, the Ott ansatz [46] was used. The technique of using polynomial expansion for an equation-free approach on the all-all Kuramoto model has also been studied in [40, 41], where the state of the oscillators is expressed as a function of only one of its heterogeneous parameters, i.e., frequencies. In one study, Bertalan (2017) expressed the state of the oscillator in terms of two heterogeneous parameters: intrinsic frequencies and the degree of connectivity. However, our idea is to express the state of oscillators as a function of two of its heterogeneous parameters: frequency and position. For the quickest convergence of the expansion, although we initially used 6 shifted Legendre polynomials, later we use different types of expansions in our study (Details are discussed in Chapter 3 and Chapter 4). The outline of this thesis is as follows: In Chapter 2, we provide a preliminary analysis of the Kuramoto model, which we will use as a benchmark in the later chapters. In Chapter 3, we formally introduce the equation-free framework and demonstrate it with coarse-projective integration. In Chapter 4, we use the equation-free approach in one dimension (when θ is a function of i and ω only) to generate bifurcation diagrams. Finally, in Chapter 5, we improve the results calculated in Chapter 4 by using the equation-free approach in two dimensions (when θ is a function of both ω and i), and in Chapter 6, we present the conclusions. Chapter 2 Preliminary Analysis of the Kuramoto Model In this chapter, we analyze the all-all Kuramoto model in which each oscillator is connected with every other oscillator and the spatially extended Kuramoto model in which each oscillator is connected to its neighbouring oscillators. We first show that when the coupling range, M, is large enough, the spatially extended network is the same as the all-all network. We then discuss two different methods to calculate the stationary states. Lastly, we use the numerical continuation technique to create bifurcation diagrams when the coupling strength, K, varies. 2.1 From spatially extended network to all-all Kuramoto model Our model of all-all Kuramoto model is defined as θ̇i = ωi + K N−1 N ∑ j=1 sin(θ j−θi), i = 1, . . . ,N (2.1) Here, θi and ωi are the phase and natural angular velocity of oscillator i respectively, N is the number of oscillators, and K ≥ 0 is the coupling strength. We divide K by (N−1), rather than by N, as has been done, for example, in [1, 56]. This is because when i = j, sin(θ j−θi) = 0, hence, we divide the coupling strength by the number of actual connections. As we will illustrate later in the section, this small modification enables us to show that the all-all network is a special case of the spatially extended network. 8 We define the spatially-extended Kuramoto model as: θ̇i = ωi + K 2M j=M ∑ j=−M sin(θi+ j−θi), i = 1, . . . ,N (2.2) Here again, N represents the number of oscillators, θi is the phase of oscillator i, K ≥ 0 is the coupling strength, M is the coupling range and ωi is the natural frequency drawn from a uniform distribution. In contrast to the previous system (all-all Kuramoto model) where each oscillator was connected to every other oscillator, here, each oscillators is connected to 2M other oscillators. When M is large enough the spatially extended network becomes an all-all. As an example, consider a spatially extended network of N = 6 oscillators. When M = 1 each oscillator is connected to its very next neighbour on each side, see Figure 2.1(a). When M = 2 each oscillator is connected to its two closest neighbours, see Figure 2.1(b). However, once the value of the coupling range M gets higher than half of the total number of oscillators, (in this case, it is M = 3), the model becomes all-all Kuramoto model, as shown in Figure 2.1(c). Therefore, for a spatially extended network, if M > ( N 2 − 1 ) , then we consider it as all-all Kuramoto model with the edges counted not more than once. Figure 2.1 From spatially extended to all-all network (a) Each oscillator is connected to its closest neighbour. (b) Each oscillator is connected to its two closest neighbours. (c) Each oscillator is connected to all other oscillators. This example shows N = 6 oscillators. 9 2.2 Calculating the stationary solutions 0 2 4 6 8 10 12 0 1 2 3 4 5 6 1 2 3 4 5 6 rad/sec=0.5 t θ i Figure 2.2 State of oscillators for a spatially extended Kuramoto network in a stationary frame for N = 6 oscillators with coupling range M = 2. Figure 2.2 shows a solution of a spatially extended Kuramoto network with N = 6 and M = 2. As can be seen, the oscillators θi change with time. In order to calculate the periodic steady state, we use two methods to make our system stationary. In the first method, we make our system stationary by replacing the stationary frame with a rotating frame. In the second method, we consider the relative positions between the oscillators. The two methods are implemented with our own software and the results are then verified with CLMATCONT-L (an advanced Matlab package used for studying the dynamical systems and their bifurcation of large-scale problems). Method-1 (Analyzing spatially extended network in a rotating frame.) We demonstrate this method on a network of six oscillators. We solve Eq. (2.2) with M = 2 and natural angular velocity uniform randomly distributed in the interval [0,1]. We initially placed all of the oscillators within an interval [0,2π] distributed randomly. Then, using the ode45 subroutine in MATLAB, we solve the system of equations, as shown in Figure 2.2. 10 Figure 2.2 shows the state of each oscillator as it changes with respect to time. We can see that after a short transition time, all the oscillators start moving with the same angular speed Ω. To make our system stationary, we subtract Ω from our system of equations, (Θ̇i = θ̇i−Ω) . This creates a new set of phases, Θi on the rotating frame. Our modified model of a spatially extended network then becomes Θ̇i = ωi + K 2M j=M ∑ j=−M sin(θi+ j−θi)−Ω, i = 1, . . .N (2.3) Because we have a degree of freedom, we may choose where to attach the rotating frame. Here, "degree of freedom" refers to the flexibility to choose the reference frame for the system. In the modified model (Eq. (2.3)), a rotating frame is introduced by subtracting the common angular speed Ω from each oscillator’s phase, and the choice to attach this frame at Θ1 = 0 establishes a consistent phase reference for the oscillators. The simulations of Eq. (2.3) is shown below in Figure 2.3. 0 2 4 6 8 10 12 0 2 4 6 1 2 3 4 5 6 t Θ i Figure 2.3 State of oscillators with respect to time for spatially extended Kuramoto network in a rotating frame for N = 6 and M = 2. Each coloured symbol in Figure 2.3 indicates the state of one of the oscillators in the rotating frame. As oscillators are arranged in a ring, the two lines generated at 0 and 2π are 11 considered to be one. We can observe in Figure 2.3 that in a rotating frame, all the oscillators reached a steady state. One of the interesting results of our computation is the rate at which the pattern drifts. For both our networks, we can show that the speed of the oscillators Ω is identical to the average frequency of the oscillators. This is proved here for any generalized case. 2.2.1 Theorem : ""Synchronization in spatially extended and all-all Kuramoto net- works: The equalization of angular velocities" For the spatially extended and all-all Kuramoto networks, the synchronized angular velocity of the oscillators is equal to the average frequency of the oscillators. Figure 2.4 A general description of a connection within a network. Connection X is shown between two nodes XL and XR. Proof. X represents the connection that links the two nodes, which we have identified as the left node, XL, and the right node, XR, shown in Figure 2.4. Each of the nodes XL and XR will have a differential equation associated with it in the form of Eq. (2.1) and Eq. (2.2) that includes the terms Node XL : sin(XL−XR) (2.4) Node XR : sin(XR−XL) (2.5) Since sin(−θ) =−sin(θ), Eq. (2.5) can be written as: sin(XR−XL) =−sin(XL−XR) (2.6) This will be true for all the connections in the network. Hence, after summing up, all the differential equations involving all the sine terms cancel each other, resulting in: Θ̇1 + Θ̇2 + Θ̇3 + . . .+ Θ̇N = (ω1 +ω2 +ω3 + . . .+ωN)−NΩ (2.7) In a steady state, Θ̇i = 0 i = 1, . . .N (2.8) 12 which implies, Ω = ∑ N i=1 ωi N . (2.9) Thus, in a steady state of spatially extended and all-all Kuramoto networks, the synchronized speed with which all the oscillators rotate is the same as the average natural velocities of all the oscillators. Method-2 (relative distances in the spatially extended network.) For our rotating pattern shown in Figure 2.2, we can observe that the distance between the oscillators remains constant because they move at the same speed. Therefore, in this method, instead of looking at the rotating frame, we look at the relative distances between two oscillators. Let xi denote the relative distance between each oscillator and θ1 xi = θi+1−θ1 i = 1, ...,N−1 (2.10) The differential equations then become: ẋi = θ̇i+1− θ̇1 i = 1, . . . ,N−1 (2.11) Solving Eq. (2.11), we obtain Figure 2.5, which illustrates the stationary state of the Kuramoto Model spatially extended network. In this case, M = 2, N = 6, and the angular velocities are uniformly distributed between [0,1]. 13 0 2 4 6 8 10 0 2 4 6 x 1 x 2 x 3 x 4 x 5 t x i Figure 2.5 Stationary state of the spatially extended Kuramoto model using method-2 with N = 6 and M = 2. In the case of all-all Kuramoto model, we can easily find the generalized form of the differential equations: ẋ1 = (ω2−ω1)+ K (N−1) [−2sin(x1)+ sin(x2− x1)+ sin(x3− x1)+ . . .+ sin(xN−1− x1) −sin(x2)− sin(x3)− . . .− sin(xN−1)] ẋ2 = (ω3−ω1)+ K (N−1) [−2sin(x2)+ sin(x1− x2)+ sin(x3− x2)+ . . .+ sin(xN−1− x2) −sin(x1)− sin(x3)− . . .− sin(xN−1)] ẋ3 = (ω4−ω1)+ K (N−1) [−2sin(x3)+ sin(x1− x3)+ sin(x2− x3)+ . . .+ sin(xN−1− x3) −sin(x1)− sin(x2)− sin(x4)− . . .− sin(xN−1)] ... ẋi = (ωi+1−ω1)+ K (N−1) [ sin(−xi)+∑ N−1 j=1 j ̸=i sin(x j− xi)−∑ N−1 j=1 sin(x j) ] , i = 1, ...,N−1 (2.12) In the spatially extended case, we solve the differential equations numerically as illustrated in Figure 2.6. We start the process by initializing the values of xi. With the help of Eq. (2.10), we calculate θi, next we evaluate θ̇i from Eq. (2.2). Once, we calculate θ̇i, then the last step is to switch from θ̇i’s to ẋi’s by using Eq. (2.11). 14 ẋi(t0) Input θ1 = 0 θi+1(t0) = xi(t0)+θ1(t0) θi(t0) θ̇i(t0) Output ẋi(t0) Figure 2.6 Wrapping process for finding the differential equations numerically in spatially extended networks. 2.2.2 Computing the fixed points and their stabilities Stable fixed points in both methods 1 and 2 indicate a synchronized solution. We now describe how we calculate fixed points and their stability in both methods. Fixed points for stationary pattern in a rotating frame (method-1) To find the fixed points for our system Eq. (2.3), we put Θ̇i = 0, i = 1, ...,N (2.13) Or, Θ̇i = fi(θθθ)−Ω = 0, i = 1, ...,N (2.14) where fi(θθθ) is given in Eq. (2.2). By using the Newton-Raphson method, we solve the system of equations for a stationary pattern. θ1 θ2 θ3 θ4 ... θN  (n+1) =  θ1 θ2 θ3 θ4 ... θN  (n) − JJJ(−1)  f1(θ1,θ2,θ3, ...,θN−1,θN)−Ω f2(θ1,θ2,θ3, ...,θN−1,θN)−Ω f3(θ1,θ2,θ3, ...,θN−1,θN)−Ω f4(θ1,θ2,θ3, ...,θN−1,θN)−Ω ... fN(θ1,θ2,θ3, ...,θN−1,θN)−Ω  (n) (2.15) Or, θθθ (n+1) = θθθ (n)− (JJJ−1)( fff (θθθ (n))−ΩΩΩ) (2.16) 15 In Eq. (2.16), n indicates the iteration number, JJJ is the Jacobian of fi(θθθ) and a superscript -1 indicates the inverse of the Jacobian. We can calculate the Jacobian as: Jik = ∂ fi(θθθ) ∂ θk i = 1, ...,N k = 1, ...,N (2.17) Thus, we can find the fixed point of our stationary pattern. The stability of our system is calculated by examining the eigenvalues of the Jacobian. If all the eigenvalues have negative real parts, then the fixed point is asymptotically stable. Fixed point for relative positions (method-2) The fixed point for Eq. (2.11) can be found using the Newton-Raphson method. In case of all-all Kuramoto model, we used the generalized form defined in Eq. (2.12) to find the Jacobian of the steady state. For all xi defined in Eq. (2.12), let g function be defined as ẋi = gi(xxx), i = 1, ...,N−1 (2.18) Then the Jacobian is defined as: Jik = ∂gi ∂xk , i = 1, ...,N−1 k = 1, ...,N−1 (2.19) In a summation form, the Jacobian can be written as: For the diagonal values: Jii = K (N−1) [ −2cos(xi)− N−1 ∑ p=k p̸=i cos(xp− xi) ] , i = 1, ...,N−1 (2.20) For the off-diagonal values: Jik = K (N−1) [ (cos(xk− xi)− cos(xk)) ] , i = 1, ...,N−1 k = 1, ...,N−1 (2.21) In the case of spatially extended network, we apply the same wrapping approach as illustrated in Figure 2.6 to compute the Jacobian.Let ĴJJ be the Jacobian with respect to the state variables θi and JJJ be the Jacobian with respect to xi. We expand the trigonometric term involved in Eq. (2.2) in terms of sine and cosine and use the connectivity matrix Ai j (already defined in Eq. (1.2)) to calculate the ĴJJ. 16 For the off-diagonal values: Ĵik = (−sinθi) ( Aik d dθi cosθi ) +(cosθi) ( Aik d dθi sinθi ) , i = 1, . . .N, k = 1, . . .N (2.22) For the diagonal values: Ĵii = d dθi (−sinθi) ( Aik cosθi ) + d dθi (cosθi) ( Aik sinθi ) , i = 1, . . .N (2.23) Using Eq. (2.11), Eq. (2.18) can be expressed as ẋi = gi(θθθ(xxx)), i = 1, ...,N−1 (2.24) With ẋi = gi(xxx), θ̇i = fi(θθθ), i = 1, ...,N−1, we can express Eq. (2.19) as Jik = ∂gi ∂xk = ∂gi ∂θk+1 = [ ∂ fi+1 ∂xk+1 − ∂ f1 ∂xk+1 ] , i = 1, ...,N−1, k = 1, ...,N−1 (2.25) Or, equivalently, Jik = Ĵ(i+1,k+1)− Ĵ(1,k+1), i = 1, ...,N−1, k = 1, ...,N−1 (2.26) Hence, by using the relationship between J and Ĵ as defined in Eq. (2.26), we effectively switch from θi to xi. 2.3 Numerical continuation of the Kuramoto models The coupling strength, K, plays a crucial role in the existence and stability of fixed points in both of our Kuramoto networks. Previous work on the Kuramoto model (covered in Chapter 1) shows that by decreasing the value of one of the system’s parameters, there exists a point where the solution is no longer stable and the system’s state changes from coherent to incoherent. Everything could be synchronized and the solution converges, but altering just one parameter, K in our example, causes the entire network to become unstable. Thus, in order to fully understand the dynamics of our models, we investigate their stability by incorporating variations in the coupling strength K. This can be accomplished in two different ways. One method involves integrating forward with some initial condition in time until it reaches a steady state. Then we slightly change the K value and repeat the process. Although 17 the process is simple, there are many drawbacks. First, it takes a lot longer to reach a steady state. Secondly, this approach can only identify stable fixed points. Last but not least, for certain fixed K values, several stable fixed points could co-exist, and which one is reached depends entirely on the initial conditions we introduced. The second method, which is more efficient, is the technique of numerical continua- tion. Details about the numerical continuation method we use in this thesis is discussed in appendix A.1. Here we briefly review the numerical continuation technique when it applies to Eq. (2.14) and follow its solutions as the parameter K varies. 2.3.1 Methodology By using the numerical continuation technique, we follow the solutions of our spa- tially extended and all-all Kuramoto models as the parameter K varies. The system of equations defined in Eq. (2.2) and Eq. (2.1) can be written at steady state in a vector form as F(θθθ ,K) = 0 (2.27) As discussed in the previous section, we find the fixed point with the help of the Newton-Raphson method. Let (θθθ 000,K0) be one of these solutions, we would like to find another nearby solution (θθθ 111,K1), that satisfies Eq. (2.27). Following ideas discussed in appendix A.1 and Eq. (A.2), we get (θθθ 111−θθθ 000) T θ̇θθ +(K1−K0)K̇−∆s = 0 (2.28) Here, superscript T is the transpose. Eq. (2.28) states the condition that the new point (θθθ 111,K1) lies on the hyper-plane perpendicular to the tangent vector (θ̇θθ 000, K̇0) and simultaneously at a distance ∆s away from its closest point (θθθ 000,K0). We also assume that the (N +1) dimensional column vector, i.e.( θ̇θθ 000 K̇0 ) (2.29) is the null vector of the N × (N + 1) matrix (Fθθθ | FK) where subscripts indicate partial derivatives. Hence, Fθθθ is the N×N Jacobian of F with respect to θθθ and FK is the column vector of derivatives with respect to K that are evaluated at (θθθ 000,K0). Once the vector Eq. (2.29) has been identified and normalized, it can be used in Eq.(2.28). 18 Thus, solving Eq. (2.27) and Eq. (2.28) simultaneously,( θθθ (i) 111 K(i) 1 ) = ( θθθ (i−1) 111 K(i−1) 1 ) − JJJ−1 (i−1) ( F(θθθ (i−1) 111 ,K(i−1) 1 ) (θθθ (i−1) 111 −θθθ 000) T θ̇θθ 000 +(K(i−1) 1 −K0)K̇0−∆s ) (2.30) for i = 1,2, ...Ne, where Ne is defined as the number of iterations performed for convergence. Here, JJJ(i) = ( Fθθθ FK θ̇θθ 000 K̇0 ) (2.31) is the (N +1)× (N +1) Jacobian of the augmented system, and the partial derivatives are evaluated at (θθθ 111 (i),K(i) 1 ). As above, we take Ne Newtonian iterations (i.e. we assume that Eq. (2.30) has converged in such a way that θθθ 111 = θθθ (Ne) 111 , K1 = K(Ne) 1 for an initial condition θθθ (0) 111 = θθθ 000+ θ̇θθ 000∆s, K(0) 1 = K0+ K̇0∆s. The stability of the fixed point (θθθ 111,K1) depends on the eigenvalues of Fθθθ evaluated at this point, whereas this matrix Fθθθ has already been calculated as the top left N×N block in the Jacobian JJJ, see Eq. (2.31). We find the next solution (θθθ 222,K2) in exactly the same way. 2.4 Variation of the coupling strength in the spatially ex- tended network In this section, we consider the spatially extended Kuramoto model with different numbers of oscillators and vary the coupling strength, K, using the numerical continuation technique. The results are computed in MATLAB using software we wrote (provided in the appendix A.3). We start the computation by considering the number of oscillators as N = 6 with a coupling range M = 2. The random natural angular velocities are uniformly distributed in an interval [0,1]. Initially, all the oscillators are placed at [0,2π] interval. First, we calculate the average angular velocity of all the oscillators and subtract it from the system, Eq. (2.3). We then solve this system using the subroutine ode45 for some time t as shown previously in Figure 2.3. It should be noted that the same set of random angular velocities is reused for each K value, and the initial guess for the steady-state numerical integration is computed only once. The obtained steady-state solution from the direct integration serves as the initial guess for the numerical continuation technique. Subsequently, the numerical continuation technique is employed to generate a bifurcation diagram, depicting the stability and characteristics of all the stationary solutions, as the coupling strength K varies. 19 0 5 10 15 20 0 0.2 0.4 0.6 A 1 B 1 B 2 B 3 A 3 A 2 K Θ 2 1 2 3 4 5 6 0 2 4 6 8 (A1) i Θ i K = 15, λmax =−0.06 1 2 3 4 5 6 0 2 4 6 8 (B1) i Θ i K = 15, λmax = 6.01 1 2 3 4 5 6 0 2 4 6 8 (A2) i Θ i K = 5, λmax =−0.02 1 2 3 4 5 6 0 2 4 6 8 (B2) i Θ i K = 5, λmax = 1.74 1 2 3 4 5 6 0 2 4 6 8 (A3) i Θ i K = 2.5, λmax =−0.009 1 2 3 4 5 6 0 2 4 6 8 (B3) i Θ i K = 2.5, λmax = 0.70 Figure 2.7 Solutions of spatially extended Kuramoto model with N = 6 and M = 2 using numerical continuation scheme as discussed in section 2.3.1. The figure at the top shows a bifurcation diagram of oscillator 2 as a function of the coupling strength K. Solid magenta line shows the stable solution, blue dashed line is the unstable solution. Sample steady state solutions for all the oscillators along the bifurcation curve are also shown. 20 In Figure 2.7, the results of spatially extended Kuramoto model as K varies with a numerical continuation scheme is presented. Here, the magenta solid line shows the stable solution for one of the individual oscillator Θ2 in a steady state whereas the unstable solution for the same value of K are represented by blue dashed line. We investigate the overall behaviour of the state of oscillators for several values of K i.e. K = 15, K = 5, K = 1 (highlighted by A1, A2, A3 on the stable branch and B1, B2, B3 on the unstable branch) and calculate the maximum eigenvalue at the same values. It can be seen that the system’s solution is stable as eigenvalues are negative until K = 0.5159, but as soon as we cross this point, one of the eigen value becomes positive and the solution becomes unstable. This is known as the point of saddle-node bifurcation. 2.4.1 Variation of coupling strength in the spatially extended Kuramoto model with N = 100 and M = 5. In order to see how far the state of oscillators disperse for a stable and unstable branch, we increase the total number of oscillators to N = 100, for a spatially extended Kuramoto model with randomly uniformed distributed natural angular velocities, drawn from the interval [0,1] and a coupling range M = 5. We place the oscillators with the initial conditions such that they all lie within an interval [0,2π] and then evolve θi(t) for some time. We then apply the same approach as we did before, i.e. calculating a rotating pattern in a stationary state by subtracting the average frequency from the original system, then finding the fixed points of the stationary pattern with the Newton-Raphson method, and then varying the coupling strength K using numerical continuation technique. The results are shown in Figure 2.8. The top figure in the middle of Figure 2.8 shows the state of one of the oscillators. The green line represents the stable branch whereas the blue line represents the unstable branch. The system is stable until K = 1.313 so this is the point of bifurcation. Again, we consider several values of K at different distances from the point of bifurcation. For K = 15 (A1), at the stable branch, all the oscillators appear to be within the interval [0,2π], however for the unstable branch, at the same value of K,(B1) they appear to lie approximately between [0,π]. This is also true for K = 5 and K = 2.5. 21 0 5 10 15 20 0 0.2 0.4 B 1 A 2 A 1 B 2 B 3 A 3 K Θ 2 0 50 100 0 5 10 (A1) i Θ i K = 15, λmax =−0.06 0 50 100 0 2 4 (B1) i Θ i K = 15, λmax = 6.01 0 50 100 0 5 10 (A2) i Θ i K = 5, λmax =−0.02 0 50 100 0 2 4 (B2) i Θ i K = 5, λmax = 1.74 0 50 100 0 2 4 6 (A3) i Θ i K = 2.5, λmax =−0.009 0 50 100 0 5 10 (B3) i Θ i K = 2.5, λmax = 0.70 Figure 2.8 Solutions of spatially extended Kuramoto model with N = 100 and M = 5 using numerical continua- tion scheme as discussed in section 2.3.1. The figure at the top shows a bifurcation diagram of oscillator 2 as a function of the coupling strength K. Solid green line shows the stable solution, blue dashed line is the unstable solution. Sample steady state solutions for all the oscillators along the bifurcation curve are also shown. 22 2.5 Varying twisted states In the previous section, we made an intriguing observation regarding the behavior of states in a spatially extended Kuramoto model when we increased the value of a coupling strength. Specifically, we discovered that the patterns generated by the oscillator’s phases were distributed within the range of [0,2π] on a stable branch, but shifted to the range of [0,π] on an unstable branch. This motivates us to consider setting up various initial conditions and analyse the behaviour of these states for stable and unstable branches. Thus, in the following section, we will analyse our spatially extended network by providing different initial conditions and observe the states of the oscillators. Twisted states in spatially extended Kuramoto model For the spatially extended Kuramoto model, we consider the same number of oscilla- tors, i.e. N = 100 but vary the number of twists using different values of q. We call twist-0 to solutions in which the initial condition for the state of the oscillators remain the same, i.e. θi = 0 for all i, twist-1 to solutions in which the state of oscillator phases grow for [0,2π]. Similarly, twist-2 lies on [0,4π], twist-3 lies on [0,6π] and lastly twist-4 lies on the interval [0,8π]. Hence, the general form for varying twisted states would then become 2πq where q is an integer or winding number used for varying the twists. In Figure 2.9, we present all the twisted states we calculated (twist-0, twist-1, twist-2, twist-3, twist-4) as a function of the coupling strength K for one of the oscillator states Θ2. Although their bifurcation point is different, their shape is the same, i.e. they all undergo saddle-node bifurcations. Moreover, the dashed line represents the unstable branch whereas the solid line is for the stable branch. We can also observe that the bifurcation point of twist-1 is lower than any other twists. In addition, the bifurcation point continues to increase for twist-2 and other higher twists. We can also observe that the higher twisted states exist for larger values of K. 23 0 5 10 15 20 0 0.2 0.4 0.6 Stable for q =0 Unstable for q=0 Stable for q=1 Unstable for q=1 Stable for q=2 Unstable for q=2 Stable for q=3 Unstable for q=3 Stable for q=4 Unstable for q=4 K Θ 2 Figure 2.9 The bifurcation diagram illustrates the behavior of twist-0, twist-1, twist-2, twist-3, and twist-4 (q =0, 1, 2, 3, 4) as a function of the coupling strength K. The system parameters are set to N = 100 and M = 5, and the natural angular velocities are uniformly distributed (randomly). To compare our results numerically for all the twisted states, we present the corre- sponding point of bifurcation for each twisted states in Table 2.1 below. Number of Twists Bifurcation Point q=0 1.3571 q=1 1.3131 q=2 1.7762 q=3 2.7542 q=4 5.3154 Table 2.1 Bifurcation points for various twisted states with N = 100 and M = 5 and the natural angular velocities are uniformly distributed (randomly). 24 0 50 100 0 2 4 6 8 K=15 (Stable Branch) K=15 (Unstable Branch) Twist-0 i Θ i 0 50 100 0 2 4 6 8 K=15 (Stable Branch) K=15 (Unstable Branch) Twist-1 i Θ i 0 50 100 0 2 4 6 8 K=15 (Stable Branch) K=15 (Unstable Branch) Twist-2 i Θ i 0 50 100 0 2 4 6 8 K=15 (Stable Branch) K=15 (Unstable Branch) Twist-3 i Θ i 0 50 100 0 2 4 6 8 K=15( Stable Branch) K=15(Unstable Branch) Twist-4 i Θ i Figure 2.10 The states of 100 oscillators with different twisted states from both stable and unstable branches are observed. The system parameters are set to N = 100 and M = 5, and the coupling strength is fixed at K = 15. The natural angular velocities are uniformly distributed (randomly). Figure 2.10 shows the state of the oscillators for various twisted states at K = 15. We consider the same condition, i.e. N = 100 oscillators and M = 5 for each of the twisted states. For all the twisted states, the blue pattern shows the results of the unstable branch. The stable branch of twist-0 is shown in magenta color. There exists a discontinuity for the unstable branch of twist-0. Although it appears from the unstable branch of twist-0 that the state of oscillators varies from 0 to 2π , if we place all these oscillators in a ring, it shows that there is a discontinuity of approximately 0.9852π in the middle. This suggests that all the oscillators 25 will cover less than half of the circle and eventually become closer to 0 or equivalently to 2π . Therefore, we call it the twist-0 unstable branch. For twist-1, Θi varies from 0 to 2π for the stable branch (in green color) and the oscillators cover more than half of the circle in the unstable branch, ranging from 0 to 3.365 radians. For twist-2, with all the same conditions as above, the phases of the oscillators range from 0 to 4π for the stable branch (shown in red color), and for the unstable branch, the number of full circles covered is (10.09 radians−0.416 radians) 2π ≈ 1.5364. So, the oscillators cover approximately 1.5364 full circles when arranged in a ring between 0.416 and 10.09 radians. To be even more precise, the oscillators complete one full revolution (2π radians) and then cover an additional 0.5364×2π radians, which is approximately 3.365 radians beyond the one full revolution. This means they cover one complete revolution and a little over half of another revolution. For twist-3, the oscillators cover approximately 2.6783 full circles. In the stable branch, the phases of the oscillators range from 0 to 6π radians (shown in cyan color). For the unstable branch, the number of full circles covered is (16.86 radians−0.5529 radians) 2π ≈ 2.6783. So, the oscillators cover approximately 2.6783 full circles when arranged in a ring between 0.5529 and 16.86 radians. To be even more precise, the oscillators complete two full revolutions (2π radians each) and then cover an additional 0.6783× 2π radians, which is approximately 4.259 radians beyond the two full revolutions. This means they cover two complete revolutions and a little over half of another revolution. For twist-4, the oscillators cover approximately 3.7355 full circles. In the stable branch, the phases range from 0 to 8π radians (shown in brown color). For the unstable branch, when the oscillators are arranged in a ring and cover the region from 0.5728 to 23.67 radians, they cover approximately 3.7355 full circles, which is almost three and three-quarters of a full revolution. 2.6 Varying the coupling range Previously, we were considering N = 100 and M = 5 for our spatially extended Kuramoto model. In this section, we gradually increase the coupling range, M for a spatially extended network and analyze the stability by evaluating the point of bifurcations using the same numerical continuation scheme as discussed in section 2.3.1. First, we consider a case for which N = 100 and M = 25. Then, we increase the coupling range further by considering N = 100 and M = 50. For both cases, we discuss the twisted states and the state of oscillators at various values of K. We can observe that as the coupling range increases to a certain limit, 26 the only existing twisted states is twist-0. Moreover, when the value of the coupling range becomes equal to half of the number of oscillators, the system has the same bifurcation point as for the all-all Kuramoto model. 2.6.1 N = 100, M = 25 We consider the number of oscillators as N = 100 but increase the coupling range from M = 5 to M = 25. Here in Figure 2.11, the middle bigger figure presents the bifurcation diagram for which the coupling strength varies using numerical continuation scheme. The unstable branch is presented in blue color while the pink color shows the stable branch. The point of bifurcation is at K = 2.24. We pick three different values of K to check the state of the oscillators at these values; the stable and unstable branches are shown in pink and blue color respectively, the state of the oscillator show twist-1. 27 0 5 10 15 20 0 0.1 0.2 0.3 A 1 B 1 B 2 A 2 A 3 B 3 K Θ 2 0 50 100 0 2 4 6 8 (A1) i Θ i K = 15, λmax =−0.08 0 50 100 0 2 4 6 (B1) i Θ i K = 15, λmax = 3.43 0 50 100 0 2 4 6 8 (A2) i Θ i K = 5, λmax =−0.02 0 50 100 0 2 4 6 (B2) i Θ i K = 15, λmax = 0.96 0 50 100 0 2 4 6 8 (A3) i Θ i K = 15, λmax =−0.01 0 50 100 0 2 4 6 8 (B3) i Θ i K = 15, λmax = 0.39 Figure 2.11 Solutions of spatially extended Kuramoto model with N = 100 and M = 25 using numerical continuation scheme as discussed in section 2.3.1. The figure at the top shows a bifurcation diagram of oscillator 2 as a function of the coupling strength K. Solid pink line shows the stable solution, blue dashed line is the unstable solution. Sample steady state solutions for all the oscillators along the bifurcation curve are also shown. 28 2.6.2 N = 100 and M = 50, all-all Kuramoto model We now increase the coupling range further while keeping the other conditions the same. Figure 2.12 shows the stable branch in magenta color whereas the unstable branch is shown with the blue color. The point of bifurcation lies at K = 0.6431. We can observe from the state of oscillators at various K values that the phases are having twist-0. 29 0 5 10 15 20 0 0.05 0.1 0.15 0.2 0.25 A 2 B 3 A 3 B 2 B 1 A 1 K Θ 2 0 20 40 60 80 100 0 2 4 6 8 (A1) i Θ i K = 15, λmax =−0.15 0 20 40 60 80 100 0 2 4 6 8 (B1) i Θ i K = 5, λmax = 15.16 0 20 40 60 80 100 0 2 4 6 8 (A2) i Θ i K = 5, λmax =−0.05 0 20 40 60 80 100 0 2 4 6 8 (B2) i Θ i K = 5, λmax = 5.06 0 20 40 60 80 100 0 2 4 6 8 (A3) i Θ i K = 1, λmax =−0.009 0 20 40 60 80 100 0 2 4 6 8 (B3) i Θ i K = 1, λmax = 0.81 Figure 2.12 Solutions of spatially extended Kuramoto model with N = 100 and M = 50 using numerical continuation scheme as discussed in section 2.3.1. The figure at the top shows a bifurcation diagram of oscillator 2 as a function of the coupling strength K. Solid magenta line shows the stable solution, blue dashed line is the unstable solution. Sample steady state solutions for all the oscillators along the bifurcation curve are also shown. 30 Next, we consider the same number of oscillators, with K = 15 and varying the cou- pling range starting from M = 2 to M = 50. We are interested in finding out the relationship between the twisted states and the coupling range M. Figure 2.13, shows all the coupling range values that give twist-3 solutions. We can see that twist-3 exists for all the values from 2 to 9. However, when we slightly increase the coupling range to M = 10, the state of the oscillators change- we can no longer get twist-3 and instead get twist-2. This state of oscillators exists until M = 14. After increasing the coupling range further we get the twist-1 solutions for all the M values ranging from 15 to 32. However, as soon as the coupling range increases to M = 33, twist-1 state of the oscillators turned to twist-0. In general, we can conclude that the higher twisted states exists at the lower values of M. 0 5 10 15 0 2 4 6 0.6 0.8 1 0.1 0.2 M=2 M=33 M=44 M=50 M=23 M=15 Twist-0 K Θ 2 0 5 10 15 0 2 4 6 14 14.5 15 0.05 0.1 2 3 4 5 0.5 1 M=2 M=23 M=27 M=32 M=15 Twist-1 K Θ 2 0 5 10 15 0 2 4 6 8 10 11 12 13 14 15 0 0.2 0.4 M=2 M=14 M=12 M=10 Twist-2 K Θ 2 0 5 10 15 0.2 0.4 0.6 0.8 M=5 M=9 M=2 M=8 Twist-3 K Θ 2 Figure 2.13 Bifurcation diagrams for various twisted states. We can see that the higher twisted states exist for lower values of M. But as we increase the coupling range and the more oscillators connect to each other, the number of possible twisted states decreases until the only possible solution is twist-0. In order to check the numerical integrity of our results we compare the calculation we did by solving Eq. (2.3) where N = 100, M = 50 with our own software (discussed in appendix A.3) and the calculation done by CLMATCONT-L (an advanced numerical bifurcation software, built in MATLAB, for studying large scale dynamical problems, [6]) 31 when solving the Eq. (2.11) with N = 100, M = 50 using method-2. In Figure 2.14, we present a comparison of results from CLMATCONT-L and our own software. 0 5 10 15 20 25 0 0.1 0.2 0.3 00 LP 99 11 11.5 12 7.5 8 8.5 10-3 max=-0.05max=-0.009 max=-0.17 max=0.81 max=5.06 max=15.16 K Θ 2 Figure 2.14 A bifurcation diagram created with our own software compared to the one created with CLMATCONT-L. Our results are shown in solid red for the stable branch and the blue dashed line for the unstable branch. The CLMATCONT-L results are shown in purple for the stable branch and the cyan dashed line for the unstable branch. Here, "00" means the K value from where the bifurcation starts, "LP" means the limit point or the bifurcation point", which in this case is K = 0.6431 and "99" means the last value of K where the numerical continuation ends. As can be seen, our calculation agrees well with the calculation obtained in CLMATCONT-L. 2.7 Discussion and Conclusions In this Chapter, we numerically analyzed the transition between spatially extended Kuramoto network and all-all Kuramoto network. We discussed two methods, one in which we analyzed the system in a rotating frame and a second in which we considered the relative position of the oscillators. We used both methods to find the fixed points. We also presented a proof to show that for both networks, the synchronized angular velocity of the oscillators is equal to the average of the natural angular velocities. Finally we analyzed various twisted states in both networks and showed numerically that the twisted states exist for lower values of the coupling range, M. Thus, by gradually increasing the coupling range, more oscillators connect with each other and the only twist we get for the all-all network is twist-0. We ended the chapter by comparing our numerical results with those obtained by another software (CLMATCONT-L) and showed that we get the same results. This will become important 32 when we use our software to implement the equation free approach as later we compare the solutions obtained from an equation-free approach with the actual detailed system. Chapter 3 An Equation-Free Approach for the All-All Kuramoto Network In this chapter, we discuss how a large system of synchronized oscillators, the all-all Kuramoto network, can be represented by a smaller number of variables using an "equation- free approach". We discuss the equation-free approach and implement a coarse-projective integeration to check and improve the efficiency of our calculations. 3.1 Development We begin by making the observation (also made by [5]) that once all the oscillators synchronize, their state becomes a function of their natural frequencies ωi, as shown in Figure 3.1. The figure shows the solution of Eq. (2.1) with N = 100, coupling strength K = 0.7 whereas the frequencies are uniformly distributed between [0,1] interval. Figure 3.1 (a) shows the state of the oscillators as a function of time t and Figures 3.1 (b-c) show the state of oscillators as a function of the natural frequencies at times t = 0, t = 12 and t = 50 respectively. 34 0 10 20 30 40 50 0 2 4 6 8 t θ i( t) t = 0 t = 12 t = 50 0 0.5 1 0 2 4 6 8 0 0.5 1 1 2 3 4 5 0 0.5 1 3 4 5 6 ωi ωi ωi θ i θ i θ i Figure 3.1 Solution of the dynamical system in Eq. (2.1) with N = 100, K = 0.7, and ωi are the uniformly distributed random frequencies between [0,1]. (a) Evolution of θi as a function of time, t. (b)-(c) θ as a function of ωi at t = 0, t = 12 and t = 50 respectively. It can be seen that once the oscillators synchronized, θi becomes a function of ωi. This correlations between θi and ωi allows us to expand θi (microscopic variables) in certain polynomial classes of ωi, known as polynomial expansion, (given in the appendix A.2). Since we have taken ωi from a uniform distribution, we will use the shifted Legendre polynomials (set of functions similar to the Legendre polynomials, but orthogonal on the interval [0,1]) which provides the fastest convergence of the expansion [64]. Thus, we can write θi(t)≈ n ∑ k=1 ak(t)φk(ωi), i = 1 . . .N (3.1) 35 where φk(ωi) are the shifted Legendre polynomials, ak(t) are n low dimensional variables and θi(t) are N high dimensional variables. Given a particular detailed realization of the oscillator state, these polynomial coefficients ak’s are estimated by minimizing the quantity: N ∑ i=1 [ θi− n ∑ k=1 ak(t)φk(ωi) ]2 (3.2) This is referred to as a Restriction operator in an equation-free approach. It can be reformulated as an over-determined linear system whose least-squares solution can be calculated by the backslash operator in Matlab. We also need a Lifting operator to execute an equation-free approach, which involves considering ak with a particular ωi realization. Thus, it is defined as: θi(t) = n ∑ k=1 ak(t)φk(ωi), i = 1, . . . ,N (3.3) With these two operators, Eq. (3.2) and Eq. (3.3), we now implement an equation- free approach. Here, we present an overview of the method for applying an equation-free approach. 3.2 Methodology of applying an equation-free approach to the all-all Kuramoto network In an equation-free approach, an adequately initialised short burst of detailed simula- tions is used to estimate the appropriate quantity for the solution of low-dimensional systems. The general steps consist of: 1. Identifying a collection of macroscopic variables that sufficiently describe the coarse- grained dynamics. We assume that the dynamical system at the level of coarse ob- servables exists but we don’t have the equations explicitly. For the all-all Kuramoto network, our macroscopic variables are ak. 2. Drawing ωi randomly, and saving the selection for all future calculations. 3. Constructing a lifting operator , as defined in Eq.(3.3), mapping the coarse description to (one or more) consistent fine-scale realizations. 4. Evolving the lifted, fine-scale initial conditions for a certain time horizon. 36 5. Restricting the resulting fine-scale description to the coarse observables, i.e finding the polynomial coefficients ak of the final state variables θi(t), using the restriction operator, specified in Eq.(3.2). 6. Repeating the procedure for steps 3, 4 and 5 as needed. We follow the above steps and compare the actual solutions (without using an equation-free approach) and the approximated solution (with an equation-free approach). 0 0.5 1 3 4 5 6 ωi Θ i Figure 3.2 Snapshot of θi as a function of ωi. Blue circles represent the actual solutions of Eq. (2.1) with N = 100 and K = 0.7 whereas the red dots represent the solution of an equation-free approach using only 2 shifted Legendre polynomials (φ1 = 1,φ2 = 2ωi−1). Figure 3.2 shows a comparison between the steady state solution as calculated using Runge- Kutta method with the full network of 100 oscillators (shown in blue circles) and the steady state solution as calculated using the equation-free approach with only 2 shifted Legendre polynomials (shown in red dots). These two shifted Legendre polynomials are φ1 = 1,φ2 = 2ωi− 1. In both simulations, the value of the coupling strength is K = 0.7. As can be seen the relationship between θi and ωi is captured well with the equation-free approach. 37 3.3 Coarse-projective integration of all-all Kuramoto net- work We now briefly demonstrate the coarse-projective integration of low-dimensional variables. The method involves short simulations of the high dimensional system to derive the time projection of the low dimensional variables.Let the high dimensional variables be presented as θθθ = [θ1,θ2, ...,θN ] ∈ RN (3.4) The low dimensional polynomial coefficients are then defined as A = [a1,a2, ...,ak] ∈ Rk (3.5) Given θθθ(t1), we integrate Eq. (2.1) with N1 time steps. Using the restriction operator, specified in Eq. (3.2), we calculate A in the last two time steps. Next, we use the last two values from A to extrapolate the value of A at some time t2 in the future. We then lift from A(t2) to θθθ(t2) by using the lifting operator defined in Eq. (3.3). The procedure is repeated as shown in the Figure 3.3. It should be noted that if the cost of restricting, extrapolating and lifting is small compared to the cost of integrating the network for the entire time, this procedure will be faster than directly integrating the network. Figure 3.3 Schematic diagram of the coarse-projective integration. 38 3.3.1 Varying the integration time and projection times. In this section, we explore the effect of varying N1, (the integration time of the high dimensional system) on the coarse-projective integration when only three coarse-grained variables are used. We consider all-all Kuramoto network as discussed in section 3.1. (N = 100, K = 0.7, random uniformly distributed frequencies). 0 5 10 15 20 -5 0 5 10 N 1 N 2 (a) Time C oe ffi ci en ts 0 5 10 15 20 -5 0 5 10 15 (b) Time C oe ffi ci en ts 0 5 10 15 20 -5 0 5 10 15 (c) Time C oe ffi ci en ts 0 5 10 15 20 -5 0 5 10 15 (d) Time C oe ffi ci en ts Figure 3.4 The effects of variation in N1 and N2 on a coarse projective integration with three coarse grained coefficients a1, a2 and a3. For (a) N1 = 1 and N2 = 4, (b) N1 = 2 and N2 = 3, (c) N1 = 3, N2 = 2, (d) N1 = 4, N2 = 1. The red dots represent the results of the coarse-grained variables after applying the coarse-projective integration, while the black line shows the result of the coarse-grained variables when there is no coarse projection involved. Following the process described in the Figure 3.3, we consider n = 3, so that our choice of coarse-grained variables becomes a1, a2, and a3 to estimate θ1,θ2, ...,θ100. In order to perform the coarse-projective integration, local time derivatives of the coarse observables are estimated by using the last two values from each curve and then projecting them for some time. After the projection, we lift the coarse variables again to the consistent fine-scale realizations and use them as an initial condition for another short burst of direct detailed 39 integration. The resulting values are then restricted and the whole process repeats again. Thus, we vary N1 from 1 to 4 in each of the cases, given in Figure 3.4, with the projection size N2 varies correspondingly. In Figure 3.4 (a) N1 = 1 and N2 = 4, Figure 3.4 (b) N1 = 2 and N2 = 3, Figure 3.4 (c) N1 = 3, N2 = 2 and Figure 3.4 (d) N1 = 4, N2 = 1. Figure 3.5 shows the error between the true solution (when the full integration is carried out) and the solution we evaluated with the coarse-projective integration at the value of time t = 20. Figure 3.5(a) indicates that coarse-projection error increases as the projection length N2 increases. 1 2 3 4 -0.01 0 0.01 0.02 0.03 (a) N1 E rr or 1 2 3 4 1 2 3 4 (b) N1 N 2 Figure 3.5 (a)-Error between the solution obtained by full integration and the solution obtained by projective integration for N1 =1,2,3,4 with three coarse-grained coefficients. The smallest rectangle represents a1, the medium is for a2, while the largest rectangle is for a3. (b)-Relationship between N1 and N2 shows N2 = 5−N1. In Figure 3.6, we take the size of the projection N2 as 5,7,9,11 while keeping N1 = 3. We follow the same procedure as discussed above. The total time to be considered is 40 units. The cyan circles represent the coarse grain coefficients with the projective integration, while a full detailed simulation is given by the blue line. It can be seen in Figure 3.6 that as the length of the projection increases, the integration loses its accuracy. 40 0 10 20 30 0 5 10 15 20 25 (a) Time C oe ffi ci en ts 0 10 20 30 0 5 10 15 20 25 (b) Time C oe ffi ci en ts 0 10 20 30 0 5 10 15 20 25 (c) Time C oe ffi ci en ts 0 10 20 30 0 5 10 15 20 (d) Time C oe ffi ci en ts Figure 3.6 The effects of variation in N2 on a coarse projective integration with three coarse grained coefficients a1, a2 and a3, while considering N1 = 3. For (a) N2 = 5, (b) N2 = 7, (c) N2 = 9, and (d) N2 = 11. The cyan circles represent the results of the coarse-grained variables after applying the coarse-projective integration, while the blue line show the results of the coarse-grained variables when there is no coarse projection. In Figure 3.7, we calculate the error between the coarse variables of the true solutions and the one with the coarse projective integration at t = 40. It can be observed that the error of each coefficient a1, a2 and a3 increases as the projection duration increases. 41 N2 E rr or Figure 3.7 Error of the coarse-projective integration when N2 is taken as 5, 7, 9, 11 for a1 (smallest rectangle), a2 (medium rectangle) and a3 (the largest rectangle). 3.4 Discussion and Conclusions In this chapter, we utilized polynomial expansion to reduce the dimensionality of a large network, which contains a heterogeneous population of Kuramoto oscillators with an all-all coupling structure. The idea of estimating the polynomial coefficient when there exists a rapid correlation between the state of the oscillators and the heterogeneous frequencies was taken from Bertalan et al. [5], who used an equation-free approach on a Kuramoto model. However, the authors expanded the coarse-grained coefficients using univariate polynomials. Laing and Troy [35] also used a polynomial expansion, but they solved a different set of differential equations with the frequencies taken from a Gaussian distribution and chose Hermite polynomials to achieve faster convergence of the expansion. Xiu and Karniadakis [64] showed that for a uniform distribution of frequencies, Legendre polynomials provide faster convergence, and in the case of a gamma distribution, Laguerre polynomials are best. Hence, we considered the uniform distribution of frequencies and used the shifted Legendre polynomials with an equation-free approach to demonstrate the coarse-projective integration for a detailed system. Chapter 4 Creating Bifurcation Diagrams Using the Equation-Free Approach in One Dimension In Chapter 2, we applied a numerical continuation technique to create bifurcation diagrams of a large system of coupled oscillators. In this chapter, we use the equation-free approach to create bifurcation diagrams for the coarse-grained system. This helps us estimate the behaviour of a high dimensional system by a smaller number of variables (the coarse- grained variables) and study the overall emergent behaviour of the detailed system. First, we briefly discuss the implementation of an equation-free approach along with the numerical continuation technique. Then, we compute the bifurcation diagrams at the coarse level and compare them with the results of the detailed system- (already computed in Chapter 2). 4.1 Finding the derivatives of coarse-grained coefficients using an equation-free approach In Chapter 2, we observed that when M (coupling range) is small, θi becomes a monotonically increasing function of their positions therefore we now consider this case and express our high dimensional variables, θi as: θi(t) = n ∑ k=0 ak(t)vk(i), i = 1 . . .N (4.1) where, ak are the time dependent coarse grained coefficients, vk are the spatial functions, n is the maximum number of coarse-grained variables, i is the position of the oscillators and N is 44 the total number of oscillators. We considered vk as vk = (i−1)k, k = 0, ...,n and we choose a0 = 0. This gives us vk = 0 at i = 1 and this polynomial representation is allowed when θi goes beyond 2π (we do not use modulus of 2π as given in Figure 2.10). Eq. (4.1) is our lifting operator, denoted by L , creating the fine scale initial conditions consistent with the prescribed values of the macroscopic observable, L(aaa) = θθθ . (4.2) The restriction operator is defined as ℜ : min [ N ∑ i=1 [ θi(t)− n ∑ k=0 ak(t)vk(i) ]2] (4.3) Or, ℜ(θθθ) = aaa (4.4) Thus, given the detailed variables, θi, we find the macroscopic variables such that the relationships in Eq.(4.3) are minimized. This calculation can be done in Matlab by using a backslash operator. The two operators given in Eq. (4.1) and Eq. (4.3) must satisfy the condition of consistency, i.e ℜ(L(aaa)) = aaa (4.5) In order to find the derivatives for a given set of coarse-grained coefficients, we use the lifting operator to find the values of θi and insert them into Eq. (2.3) to find the derivatives of the detailed system in a rotating frame. We then note that θ̇i(t) = ∑ n k=0 ȧk(t)vk(i) and use the restriction operator ℜ(θ̇θθ) = ȧaa, to get the derivatives of the coarse-grained coefficients. A flow diagram is given in Figure 4.1. 45 Coarse-grained coefficients Lifting-Operator Derivatives of the detailed system in a rotating frame Restriction operator Using an equation-free approach to find the derivatives of coarse-grained coefficients Coarse-derivatives Given ሶ𝜃𝑖 = 𝑓 𝜃𝑖 , min ෍ 𝑖=1 𝑁 ሶ𝜃𝑖(𝑡) −෍ 𝑘=0 𝑛 ሶ𝑎𝑘(𝑡)𝑣𝑘(𝑖) 2 𝜃𝑖(𝑡) = ෍ 𝑘=0 𝑛 𝑎𝑘(𝑡)𝑣𝑘(𝑖) 𝑎𝑘(𝑡) ሶ𝑎𝑘(𝑡) 𝑖 = 1,… ,𝑁 Figure 4.1 A flow diagram showing how to find the derivatives for the coarse-grained system for a given set of coarse-grained variables. Here, the spatial functions are defined as vk = (i−1)k,k = 1, ...,n. 4.2 Finding fixed points of the coarse-grained system If the equations of the coarse-grained system were previously known, we would have been able to represent them as daaa dt = ggg(aaa,K) (4.6) 46 Given the first guess aaa(1), we can find the next guess aaa( j+1) using the equation: aaa( j+1) = aaa( j)− zzz (4.7) zzz = J−1ggg(aaa( j)) (4.8) where J is the Jacobian defined as J =  ∂g1 ∂a1 ∂g1 ∂a2 . . . ∂g1 ∂an ∂g2 ∂a1 ∂g2 ∂a2 . . . ∂g2 ∂an ... ... ... . . . ... ∂gn ∂a1 ∂gn ∂a2 . . . ∂gn ∂an  (4.9) We perform the iterations until ∥ggg(aaa,K)∥<ε (ε = 10−10 in our simulations). Once we find the value of aaa for which ∥ggg(aaa,K)∥ < ε at some K = K0, we apply the numerical continuation technique described in Chapter 2 (section 2.3) to create a bifurcation diagram for the coarse-grained system. This algorithm is summarized in Figure 4.2. The very first initial guess is calculated from a stable stationary solution of the detailed system. The difficulty with the algorithm described in Figure 4.2 is that the differential equations of the coarse-grained system are unknown. These can be found as described in Figure 4.1. Therefore, to approximate the Jacobian in Eq. (4.9), we use finite differences. ∂gi ∂ak = gi(aaa+ εek,K)−gi(aaa,K) ε , i = 1, ...,n (4.10) Here, ek is a column unit vector having 1 at its k-th entry and 0s elsewhere. Recall that the rotating frame of the detailed system has a degree of freedom, which implies that the coarse-grained coefficients also have a degree of freedom. To remove this degree of freedom, we choose aaa such that θ1 = 0 (to be consistent with our choice in Chapter 2). For simplicity, we choose the spatial functions as vk = (i−1)k, k = 1, . . . ,n. 47 𝐽𝒛 = 𝒈(𝒂(𝒋)) 𝒂(𝒋+𝟏) =𝒂(𝒋) − 𝐳 Initial guess of 𝒂 at 𝐾 = 𝐾0 𝒂(𝟏) No Yes Finding 𝒂∗such that 𝒈(𝒂∗) = 0 with the Newton Raphson’s method 𝒈(𝒂(𝒋+𝟏)) < 𝜖 Is 𝑗 = 1 𝒂∗ = 𝒂(𝒋+𝟏) 𝐾 = 𝐾0 + ∆𝐾 Use numerical continuation to guess new 𝒂(𝟏) Figure 4.2 Schematic description of a numerical continuation technique for an unknown function g. Note that the function g is unknown and is estimated using the equation-free approach described in Figure 4.1. 4.3 Bifurcation Diagrams when θi are functions of their positions We start with an example of twist-1 for N = 100 oscillators and M = 5 with randomly distributed frequencies, but instead of solving the full system (as we did in Chapter 2, Figure 2.8), we apply an equation-free approach as discussed in Figure 4.2 and Figure 4.1. We increase the number of coefficients from 2 to 4 in Figure 4.3. Note that we consider the 48 spatial functions as vk = (i−1)k,k = 1, ...,n. The results of an equation-free approach are shown by dark pink color whereas the cyan color shows the results of the detailed system. We can observe that as the number of coefficients is increased, we get closer to the point of bifurcation of the detailed system, which is K = 1.313. Also, in the case of higher coefficients, the stable and unstable branches both converged to the detailed system’s stable and unstable branches respectively. 49 2 coefficients 3 coefficients 4 coefficients Bifurcation point: 1.17 Bifurcation point: 1.19 Bifurcation point: 1.21 0 5 10 -0.02 0 0.02 0.04 0.06 K a1 0 5 10 0.05 0.1 0.15 K a1 2 4 6 8 10 0.1 0.2 0.3 K a1 0 5 10 0 2 4 6 10-4 K a2 0 5 10 -6 -4 -2 0 10-3 K a2 0 5 10 -0.015 -0.01 -0.005 0 K a2 0 2 4 6 8 10 0 1 2 3 4 10-5 K a3 0 2 4 6 8 10 0 1 2 3 10-4 K a3 Key: Stable solution from equation-free approach Unstable solution from equation-free approach Stable solution from detailed system(N = 100, M = 5) Unstable solution from detailed system(N = 100, M = 5) 0 5 10 -1 -0.5 0 10-6 K a4 Figure 4.3 Comparison of the solutions of twist-1 equation-free approach (pink color) with the twist-1 detailed system (cyan color) for N = 100, M = 5 and randomly distributed frequencies between 0 and 1. The spatial functions are v1 = i−1, v2 = (i−1)2,v3 = (i−1)3,v4 = (i−1)4 where i is the oscillator number. The point of bifurcation for the detailed system is K = 1.313, whereas the point of bifurcation with an equation-free approach for 2, 3 and 4 coefficients are K = 1.17, K = 1.19 and K = 1.21 respectively. It can be seen that for four coefficients, the solution converges to both the stable and unstable branches. In Figure 4.4, we show one of the solutions from Figure 4.3 along with the fitting at K = 10 of the stable and unstable solutions marked by A1 and B1 respectively. The blue circles represent the phase of the oscillators in a moving frame of the detailed system, and the red dots represent the estimated solution from the equation-free approach. Clearly, for 50 both the stable and unstable branches, the fitting is good. Note that for the stable branch, the range of states of the oscillators is [0,2π] while the range of the states of the unstable branch is [0,π]. This illustrates that the distinction between the different twists is not a clear cut. 2 4 6 8 10 0.1 0.2 0.3 B1 A1 K a 1 0 50 100 0 2 4 6(A1) i θ i K = 10, (Stable Branch) 0 50 100 0 1 2 3 4(B1) i θ i K = 10, (Unstable Branch) Key: Stable solution from detailed system (N = 100, M = 5) Unstable solution from detailed system (N = 100, M = 5) Stable solution from equation-free approach Unstable solution from equation-free approach Solution of twist-1 detailed system Fitting from twist-1 equation-free approach Figure 4.4 Comparison of the fitting of the stable and unstable branches for twist-1 at K = 10 (marked by A1 and B1) when the number of coefficients is 4. Next, we show an example of twist-2. Again, we estimate our solution for N = 100 oscillators and M = 5. The angular frequencies are equally distributed within the interval [0.3,0.7]. The point of bifurcation for the detailed system is K = 2.09, whereas the point of bifurcation with the 2, 3 and 4 coefficients is K = 1.56, K = 1.74 and K = 1.76 respectively. Figure 4.5 shows the comparison of the detailed system (turquoise color) with an equation- free approach (pink color). 51 2 coefficients 3 coefficients 4 coefficients Bifurcation point: 1.56 Bifurcation point: 1.74 Bifurcation point: 1.76 2 4 6 8 10 -0.2 -0.1 0 0.1 K a1 2 4 6 8 10 -0.3 -0.2 -0.1 0 0.1 K a1 2 4 6 8 10 -0.4 -0.3 -0.2 -0.1 0 K a1 2 4 6 8 10 0 5 10 10-4 K a2 2 4 6 8 10 K 2 4 6 8 10 12 14 10-3 K a2 2 4 6 8 10 5 10 15 10-3 K a2 2 4 6 8 10 -1 -0.5 0 10-4 K a3 2 4 6 8 10 -1 -0.5 0 10-4 K a3 Key: Stable solution from equation-free approach Unstable solution from equation-free approach Stable solution from detailed system (N = 100, M = 5) Unstable solution from detailed system (N = 100, M = 5) 2 4 6 8 10 -10 -5 0 10-8 K a4 Figure 4.5 Comparison of the bifurcation diagrams of twist-2 obtained by the equation-free approach and the detailed system. The natural frequencies are equally distributed between 0.3 and 0.7. The spatial functions are i−1, (i−1)2,(i−1)3,(i−1)4 where i is the oscillator number. The point of bifurcation for the detailed system is K = 2.09, whereas the point of bifurcation with an equation-free approach for 2, 3 and 4 coefficients is K = 1.56, K = 1.74 and K = 1.76 respectively. 52 2 4 6 8 10 -0.4 -0.3 -0.2 -0.1 0 B 1 A 1 K a 1 0 50 100 0 5 10 (A1) i θ i K = 10, (Stable Branch) 0 50 100 0 5 10 15 (B1) i θ i K = 10, (Unstable Branch) Key: Stable solution from detailed system (N = 100, M = 5) Unstable solution from detailed system (N = 100, M = 5) Stable solution from equation-free approach Unstable solution from equation-free approach Solution of twist-2 detailed system Fitting from twist-2 equation-free approach Figure 4.6 Comparison of the fitting of the stable and unstable branches for twist-2 at K = 10 (marked by A1 and B1) when the number of coefficients is 4. We can observe in Figure 4.5 that for twist-2, the stable branch of the equation-free approach converges to the stable branch of the detailed system but the unstable branch does not. Moreover, the point of bifurcation for an equation-free approach is always ahead of the detailed system in all these cases. Additionally, our analysis of the fitting depicted in Figure 4.6, indicates that the stable branch adheres to twist-2 behavior, whereas the unstable branch does not. We will now investigate the possibility that the equation-free solution for the unstable branch follows other twists. 53 2 4 6 8 10 -0.3 -0.2 -0.1 0 0.1 twist-4 twist-3 Equation-free approach solution twist-2 twist-1 twist-0 (a) K a 1 2 4 6 8 10 -0.3 -0.2 -0.1 0 0.1 twist-4 twist-3 Equation-free approach solution twist-2 twist-1 twist-0 (b) K a 2 2 4 6 8 10 -0.3 -0.2 -0.1 0 0.1 twist-4 twist-3 Equation-free approach solution twist-2 twist-1 twist-0 (c) K a 3 2 4 6 8 10 -0.3 -0.2 -0.1 0 0.1 twist-4 twist-3 Equation-free approach solution twist-2 twist-1 twist-0 (d) K a 4 Figure 4.7 Bifurcation diagrams of twist-2 obtained by the equation-free approach (with four coefficients, pink line) and other twisted solutions of the detailed system with the frequencies equally distributed between 0.3 and 0.7. The spatial functions are (i−1), (i−1)2,(i−1)3,(i−1)4 where i is the oscillator number. In Figure 4.7, we consider the case of 4 coefficients from Figure 4.5 and compare the solutions with several twisted states of the detailed system. As can be seen the solutions of all the twisted states are very close by. For a1, the stable branch of twist-2 in both the detailed and coarse-grained systems are indistinguishable but the stable solution of the coarse-grained system continues towards twist-1 while the unstable branch then converges to the stable solution of twist-3. For a2, the unstable branch of twist-2 obtained by the equation-free approach converges to the solution of twist-1. For a3, it is in the middle of twist-3 and twist-1 while a4 seems to converge to twist-3 but is negligible. The solution we see for twist-2 might capture an emergent property of the coarse-grained system. Next, we look at the fitting for one of the solutions at K = 10 (marked by A1). In Figure 4.8, we plot a solution of twist-3 obtained by the detailed system and compare it with the equation-free solution of twist-2. We can observe that the equation-free twist-2 solution is getting closer to detailed system twist-3. 54 2 4 6 8 10 -0.3 -0.2 -0.1 0 0.1 twist-2 twist-3 twist-4 twist-1 twist-0 A1 K a 1 0 50 100 0 5 10 15 20 (A1) K a 1 Key: Solution of stable twist-3 detailed system Solution of unstable twist-2 equation-free approach Figure 4.8 Comparison of the fitting of the stable and unstable branches at K = 10 (marked by A1) for twist-2 when the number of coefficients is 4. 4.4 Bifurcation Diagrams when θi are functions of their angular frequencies In this section, we return to the case where in a synchronized state, θi, become a function of their natural angular frequencies, ωi, as we considered in Chapter 3. The algorithm for creating a bifurcation diagram in this case is similar to the algorithm we described in section 4.3 but with vk(i) replaced with φk(ωi). To remove the degree of freedom and for simplicity we choose the functions φk = ωk instead of the Legendre functions we used in Chapter 3. For convenience, we choose ω1 = 0 with a0 = 0 and the other ωi distributed equally between 0 and 1. Figure 4.9 shows the bifurcation diagrams of a spatially extended network with N = 100, M = 5. 55 2 coefficients 3 coefficients 4 coefficients Bifurcation point: 4.53 Bifurcation point: 4.80 Bifurcation point: 4.804 5 6 7 8 9 10 0 5 10 K a1 5 6 7 8 9 10 -30 -25 -20 -15 -10 -5 K a1 5 6 7 8 9 10 -30 -25 -20 -15 -10 -5 K a1 4 6 8 10 K 0 5 10 15 K a2 5 6 7 8 9 10 40 60 80 100 K a2 5 6 7 8 9 10 60 80 100 K a2 5 6 7 8 9 10 -70 -60 -50 -40 -30 K a3 5 6 7 8 9 10 -80 -60 -40 K a3 Key: Stable solution from equation-free approach Unstable solution from equation-free approach Stable solution from detailed system(N = 100, M = 5) Unstable solution from detailed system(N = 100, M = 5) 5 6 7 8 9 10 -2 0 2 4 K a4 Figure 4.9 Comparison of the solutions of twist-1 equation-free approach (pink) with the twist-1 detailed system (blue) with equally distributed frequencies between 0 and 1. The functions are φk = (ω)k,k = 1, ...4, where ω is the angular frequencies of the oscillators. The actual point of bifurcation for the detailed system is K = 4.865 Here, in Figure 4.9, the solution of twist-1 equation-free approach is represented by the pink color, whereas the blue one represents the solution of twist-1 detailed system. The actual point of bifurcation for the detailed system is K = 4.865. We can observe from these figures that the stable branch of an equation-free approach converges to the stable solution of the detailed system but this is not the case for the unstable branch. This behaviour can be observed in Figure 4.10, where we see the fitting for the detailed system at K = 10. We can 56 see that for the stable branch of twist-1 (marked by A1), the fitting is good. However, this is not the case for the unstable branch (marked by B1). 5 6 7 8 9 10 -30 -25 -20 -15 -10 -5 A 1 B 1 K a 1 0 0.5 1 0 2 4 6 (A1) i θ i K = 10, (Stable Branch) 0 0.5 1 0 5 10(B1) i θ i K = 10, (Unstable Branch) Key: Stable solution from detailed system (N = 100, M = 5) Unstable solution from detailed system (N = 100, M = 5) Stable solution from equation-free approach Unstable solution from equation-free approach Solution of twist-2 detailed system Fitting from twist-2 equation-free approach Figure 4.10 Comparison of the fitting of the stable and unstable branches at K = 10 for twist-1 when the number of coefficients is 4. Now, we consider the case N = 100 and M = 50 which is the all-all Kuramoto model and form its bifurcation diagram. The bifurcation point for a detailed system is K = 0.6431. In Figure 4.11, we show the results for twist-0 with the detailed and equation-free approach. The blue line represents the solution from the detailed system whereas the bold pink line represents the solution from the equation-free approach with 3 and 4 coefficients. The point of bifurcations are K = 0.62, K = 0.63 and K = 0.632 from 2, 3 and 4 coefficients respectively. 57 2 coefficients 3 coefficients 4 coefficients Bifurcation point: 0.62 Bifurcation point: 0.63 Bifurcation point: 0.632 0 5 10 0 5 10 K a1 0 5 10 0 10 20 30 40 K a1 0 5 10 0 10 20 30 40 K a1 0 5 10 -10 -5 0 K a2 0 5 10 -80 -60 -40 -20 0 K a2 0 5 10 -100 -50 0 K a2 0 5 10 0 20 40 60 K a3 0 5 10 0 50 100 150 K a3 Key: Stable solution from equation-free approach Unstable solution from equation-free approach Stable solution from detailed system (N = 100, M = 50) Unstable solution from detailed system (N = 100, M = 50) 0 10 20 30 40 50 -60 -40 -20 0 K a4 Figure 4.11 Comparison of the solutions of twist-0 equation-free approach (pink) with the twist-0 detailed system (blue) with the equally distributed frequencies between 0 and 1. The functions are φk = (ω)k,k = 1, ...,n, where ω is the angular frequencies of the oscillators. We can observe that in all the subfigures of Figure 4.11, the stable branches from the detailed system and the equation-free approach converge but this is not the case for the unstable branch. For the unstable branch, the equation-free solution of 3 and 4 coefficients converge to each other. This is shown in Figure 4.12, where we can see the overlapping solutions of 3 and 4 coefficients for a1, a2 and a3. 58 0 5 10 0 10 20 30 40 K a 1 0 5 10 -80 -60 -40 -20 0 K a 2 0 5 10 0 20 40 60 K a 3 Key: Stable branch with 4 coefficients at K = 10 Unstable branch with 4 coefficients at K = 10 Stable branch with 3 coefficients at K = 10 Unstable branch with 3 coefficients at K = 10 Figure 4.12 Overlapping solutions of 3 and 4 coefficients for a1,a2,a3. The stable and unstable branch with 3 and 4 coefficients converges for a1,a2,a3. Now, we investigate the reason for why the unstable branch of the equation-free approach and the detailed system does not converge. By considering the same value K = 10, from the stable and unstable branch (marked by A1 and B1), we check the fitting for a detailed system. We can observe in Figure 4.13 that there exists an outlier in the all-all Kuramoto model. Because of this outlier, the unstable branches of the detailed and equation-free solutions do not converge. 59 0 5 10 0 10 20 30 40 A 1 B 1 K a 1 0 0.5 1 0 0.05 0.1(A1) i θ i K = 10, (Stable Branch) 0 0.5 1 0 2 4 6(B1) i θ i K = 10, (Unstable Branch) Key: Stable solution from detailed system (N = 100, M = 50) Unstable solution from detailed system (N = 100, M = 50) Stable solution from equation-free approach Unstable solution from equation-free approach Solution of twist-0 detailed system Fitting from twist-0 equation-free approach Figure 4.13 Relationship between the state of the oscillators and the their natural frequencies for the all-all Kuramoto model twist-0. Top figure shows the bifurcation diagram. (B1) Unstable solution at K = 10. (A1) Stable solution at K = 10. The blue circles represent the oscillators in a detailed system (Eq. (2.1)), and the red dots represent the estimated solution from the equation-free approach. 4.5 Discussion and Conclusions In this chapter, we explained the methodology of applying an equation-free approach along with the numerical continuation technique to form bifurcation diagrams. We used two representations for θi and discussed both the spatially extended network and all-all Kuramoto model. We considered an example (see Figure 4.3, and Figure 4.9) where the equation-free solution captures the behaviour of the detailed system well. And some other examples where 60 it failed to do so (see Figure 4.5). We have also seen the emergent behaviour of the coarse grained system where the solution of one of the twisted state follows the solution of other twisted states. Lastly, we also observed that in the all-all Kuramoto model, the stable branch from the equation-free approach and the detailed system converge but the unstable branch does not due to the existence of an outlier. Chapter 5 Creating Bifurcation Diagrams Using the Equation-Free Approach in Two Dimensions In this chapter we attempt to improve the results of the previous chapter by considering the equation-free approach in two dimensions using two summations. We first motivate our approach with an example. We then explain the approach and show some examples with M = 5 and M = 10 for different twisted solutions. 5.1 Bifurcation Diagrams when θi are functions of both positions and angular frequencies. Consider the example of twist-1 from Chapter 4 (section 4.3) when N = 100, M = 5 and the number of coefficients are 3 and 4 with the spatial function vk = (i− 1)k. We already observed in Figure 4.4 that for the given example the fitting is good at K = 10 for both the stable and unstable branches. Figure 5.