耦合網路的同步研究

國 立 交 通 大 學

應 用 數 學 系

博

士

論

文

耦合網路的同步研究

The Study of Synchronization in Coupled

Networks

研究生： 梁育豪

指導教授： 莊重 教授

耦合網路的同步研究

The Study of Synchronization in Coupled

Networks

研究生：梁育豪 Student: Yu-Hao Liang

指導教授：莊重 Advisor: Jonq Juang

國 立 交 通 大 學

應 用 數 學 系

博 士 論 文

A Dissertation

Submitted to Department of Applied Mathematics College of Science

National Chiao Tung University in Partial Fulfillment of Requirements

for the Degree of Doctor of Philosophy

in

Applied Mathematics

September 2011

I

耦合網路的同步研究

研究生：梁育豪 指導教授： 莊 重 教授

國立交通大學

應用數學系

摘 要

本篇論文的目的是為研究一些常見耦合網路的同步化現象。對於

混沌耦合網路，基於矩陣測度的概念，我們發展出一些同步化平面的

全域穩定性判別定理。對比於其他研究，本篇論文所建立的方法是較

容易做驗證的：我們僅需確認單一未耦合系統的向量場結構，便能辦

別出耦合系統是否得能達成同步。此外，此處所考慮的耦合網路涵蓋

相當的廣泛，包含了時變的耦合網路。特別地，為了顯示我們發展出

的 定 理 其 可 用 性 以 及 其 和 生 物 模 型 的 相 關 性 ， 我 們 將 拿

Hindmarsh-Rose 神經元的耦合系統作為一同步化的探討例子（在這

裡，我們也將顯示一些神經元達成同步化時會出現的一些有趣動態行

為）。

對於耦合映像晶格系統(CMLs)，我們考慮同步化平面的鄰域穩定

性問題。我們將給出此穩定性的充分必要條件，並且提供了一個簡單

II

的建構

同步化曲線

演算法。利用此同步化曲線，我們可以容易地解釋

波長分歧與系統大小相依性的問題。進一步地，我們考慮小波變換方

法對 CMLs 所造成的影響。我們將在數值上以及理論上顯示出此變換

方法在 CMLs 模型中是可以提高同步化現象的發生。

III

The Study of Synchronization in Coupled

Networks

Student: Yu-Hao Liang Advisor: Prof. Jonq Juang

Department of Applied Mathematics

National Chiao Tung University

Abstract

The purpose of this thesis is to study the onset of synchronization in

some common models. For the coupled chaotic systems, based on the

concept of matrix measure, we derive some criteria for the global stability of

the synchronous manifold. Comparing with other developed criteria, our

criteria are easy to apply. By merely checking the structure of the vector

field of the single oscillator, we shall be able to determine if the coupled

system could acquire synchrony. Moreover, the considered coupled

networks could be quite general including the time-varying networks.

Specifically, to illustrate the applicability of our developed criteria and to be

much relative to the real phenomena, the coupled Hindmarsh-Rose neurons

are considered. Some interesting phenomena under synchronization are

showed. For the coupled map lattices (CMLs), ones consider the local

stability problem. We shall give the necessary and sufficient conditions for

the local stability of the synchronous manifold for the arbitrary coupled

networks. Also, a simple algorithm for the construction the so-called

synchronization curve is given. Using it, one can simply explain the

phenomena of the wavelength bifurcation and the size dependence problem.

In addition, the effect of the wavelet transform method developed by Wei

and et al is considered for the CMLs. We shall show that both numerically

and theoretically, it enhances the onset of the synchronization.

IV

誌 謝

首先感謝我的指導教授 莊重老師，在我六年的交大學生生涯中

不論是在研究、課業、生活上的種種幫助。在研究有困難時，感謝老

師您對我的鼓勵；感謝老師您經常鼓勵我去參加國內外的研討會，讓

我更能看清研究本質、培養國際觀。另外，在與老師相處的這些年中，

老師的研究態度深深的影響著我。相信這對於我未來的研究是會有很

大得助益的。真的很感謝老師您這些年來的照顧。另外，也感謝口試

委員 林文偉教授、郭忠勝教授、李明佳教授以及陳國璋教授，在論

文口試時給我的寶貴意見，使論文更甄完善。而這些建議更提供我作

為未來研究生涯中的一個研究方向。

在六年的交大應數系學習中，深深得發覺交大老師們對於學生有

著濃濃得關懷。感謝許元春老師、賴明治老師經常地鼓勵我應該到國

外看看別人的研究，也願意提供給我一些出國的管道。雖然最後我還

是待在交大應數系裡，但是老師們對我的關心，真的是由衷的感謝。

感謝陳秋媛老師，常常地問候我的研究情況，也常常的給予我鼓勵，

這些都讓我覺得十分的窩心，讓我覺得應數就像是個大家庭一樣。感

謝石至文老師，您教授的常微分方程、生物系統等課程，讓我對動態

系統的領域充滿了研究興趣。因為有了您的教導，我才會選擇這塊研

究領域。感謝林文偉老師，您在論文研究上給予我的許多寶貴方向、

意見，以及其他許多在學習、研究方面的幫助。另外，感謝李明佳老

師、林琦琨老師、李榮耀老師、梁耀文老師、鄧清政老師、陳宗麟老

師，感謝您們曾經給過我的指導與幫助。也感謝系辦的張麗君、陳盈

吟、張慧珊、宋雅鈴小姐們，感謝您們給予的許多幫忙，有很多地方

都麻煩您們了。

在這博士班四年的學習中，感謝許多朋友、同學、學長、學姊、

學弟、學妹：李金龍、張郁泉、張靖尉、黃俊銘、陳怡樺、洪千茹、

周芳竹、鄧仁益、魏士傑、黃仕翰、黃韋強、陳冠羽、邱鈺傑、廖康

伶、呂明杰、林光暉、蔡佳穎、林信佑（未來Google全球CEO =.=|||）、

V

蔡宗憲、葉淑娟、李芷萱、陳玉雯、方怡中、連威翔的互相提攜鼓勵，

感謝您們使得我這段求學過程充滿了樂趣、更加的多采多姿。其中，

特別感謝同學黃俊銘、魏士傑、黃仕翰，你們在我有困難時，對我的

幫助。另外，我也想特別感謝周芳竹學妹以及林信佑學弟（未來Google

全球CEO =.=|||），那段一起在研究小間的日子總是那麼令人難忘，

而你們畢業後還一直常來看我、陪我聊天、給我加油打氣。真的謝謝

你們，我會想念你們的。

最後，我要感謝我親愛的父母（梁有忠、王瓊珠）、姊姊（梁翠

洺）、弟弟（梁育瑋）、外婆（王吳秀桃）、阿姨、姨丈、舅舅、舅

媽、表弟、表妹們，謝謝你們的一路相挺、始終支持著我。另外，也

要感謝我已過世的奶奶（梁蘇灼），謝謝你從小到大一直疼愛著我。

最後，我想要感謝老天，謝謝您賜給我的一切，這一切的一切使

我的生命是這麼的美好。

梁育豪 2011/9/9 於交大

Contents

1 Introduction 1

1.1 Synchronization in lattices of coupled systems . . . 1

1.2 Modeling . . . 2

1.2.1 Model I: Coupled Chaotic Systems . . . 2

1.2.2 Model II: Coupled Map Lattices . . . 4

1.2.3 Model III: Pulse Coupled Systems . . . 5

1.3 Organization and results of the thesis . . . 7

2 Synchronization in Model I 9 2.1 Preliminaries . . . 9

2.2 Global synchronization with time-invariant coupling . . . 13

2.3 Global synchronization with time-varying coupling . . . 24

2.3.1 Matrices of the Coordinate Transformation . . . 25

2.3.2 Synchronization Criteria . . . 35

3 Applications for Model I 42 3.1 Synchronization in coupled Lorenz and coupled Duffing systems . . . . 42

3.2 Synchronization in Hindmarsh-Rose Neurons with Chemical and Elec-trical Synapses . . . 52

3.2.1 Introduction of the Hindmarsh-Rose Neurons . . . 52

3.2.2 Synchronization of the Hindmarsh-Rose Neurons . . . 54

4 Synchronization in Model II 73 4.1 Local Synchronization Criteria . . . 73 4.2 Application: Local Synchronization in coupled logistic maps . . . 79 4.3 Wavelet Transform Method for Coupled Map Lattices . . . 82

Chapter 1

Introduction

1.1

Synchronization in lattices of coupled systems

Synchronization is a common, significant phenomenon occurring in the natural world no matter from the micro view or the macro view. It depicts a group of different individuals displaying the same behaviors at the same time under certain interactions. The purpose of it is to get them to solve problems cooperatively. For instance, in the human brain, there are about 1010 _{neurons. To integrate separately processing} information in the brain, they have to synchronize their activity [14,31]. In the nights of the South-East Asia, such as Thailand, Malaysia or Borneo, male fireflies accumulate along the river banks flashing on and off simultaneously [13] in order to attract the female fireflies cooperatively.

Correspondingly, in the technology field, synchronization now and then is de-signed to occur in the experiment for some purpose, especially in the design of the electronic circuit system. In engineering, it is studied as a tool for transmitting infor-mation by using chaotic signals and monitoring dynamical systems [36].

So as an interesting, important topic, synchronization has drawn a great deal of attention and is intensively studied in many fields, including neuron science, biol-ogy, physics, engineering, and other fields of science [1,16,24,26,32–34,48,58,59,62,69, 71,73,77,84,90,94,96,99]. Consequently, several elegant theories and articles concerning synchronization have been rapidly constructed and published in the past few decades. The general approaches involved for driving the synchronization criteria are roughly the

Lyapunov function method (global results) and the master-stability function method (local results), the analysis of the transversal Lyapunov exponents calculated from the linearized equations for the perturbations transversal to the synchronous manifold.

Since in the real world the number of coupled units is usually large, the in-creasing interest in synchronization phenomena has led many researchers to consider synchronization in large networks of coupled systems with different coupling configu-rations [1,2,27,48,66,74,85,86,97,100]. As a result, one of the most important questions in synchronization phenomena is that how the coupling strengths and coupling config-urations of the network influence the stability of the synchronous state. Furthermore, one may ask the following controlling problem: Can one slightly modify the coupling configuration of coupled systems to dramatically reduce the coupling strength needed to acquire synchrony [51,52,98]? In this thesis, we shall focus on these issues by con-sidering various models that are frequently used to explain these phenomena.

1.2

Modeling

Based on the disparate individuals, the dynamics of interesting qualities, like move-ment, velocity, energy, potential, and so modeling are in the various fashions. Similarly, connection and commutation between themselves have several types. Depending upon the different models considered, there are different explanations as to why the units in the coupled systems synchronize. In this context, we will introduce three different kinds of models (Model I ∼ Model III) that are often used to describe the real collec-tive networks, and synchronization within. Specifically, the first two models (Model I ∼ Model II) are to be extensively studied. By it, in this section, we start with the introduction of the models.

1.2.1

Model I: Coupled Chaotic Systems

The model under consideration in this subsection is of continuous time with continuous state coupling. Specifically, we consider a unit with the interesting dynamics governed by a set of ordinary differential equations, saying dx_{dt = f (x, t). Here, x ∈ R}n_{, and f is a}

vector-valued function from Rn_{×R to R}n_{denoted by f (x, t) = (f}

1(x, t), · · · , fn(x, t))T. Moreover, when there are connections/commutations between a group of such m units, the induced whole dynamics under interaction is then described by

dxi dt = f (xi, t) + d · m X j=1 gij(t) Dxj, i = 1, 2, . . . , m, (1.1a)

where xi = (xi1, xi2, . . . , xin)T ∈ Rn, d is the coupling strength, D = (dij) ∈ Rn×n is the inner coupling matrix, and the quantity gij(t) describes the coupling weight from the unit j to the unit i. Let x = (x1, x2. . . , xm)T, and G(t) = (gij(t)) ∈ Rm×m. Then G(t) represents the (outer) coupling configuration of the network at time t. Equivalently, (1.1a) becomes ˙ x=    f(x1, t) .. . f(xm, t)    + d(G(t) ⊗ D)x =: F (x, t) + d(G(t) ⊗ D)x, (1.1b) where ⊗ denotes the Kronecker product.

Such type of model can been seen in some nervous systems or designs of the elec-tronic circuit systems [99]. For this model, general approaches to deal with the local stability of the synchronization state, include the master stability function-based crite-ria [1,72,74,75,81] and the matrix measure critecrite-ria [17,18]. For the global stability, the developing method includes the Lyapunov function-based criteria [2–4,8,78,100–104] and the matrix measure approach [17,18], and et. al. Among the Lyapunov function-based criteria, the connection graph approach proposed by Belykh [2–4,8] is sui generis since the proposed criteria combine some graph theories to avoid the direct compu-tation of the eigenvalues of the coupling matrix G(t). (Not surprisingly, quantities related to the eigenvalues of matrix G(t) should play the decisive roles in the synchro-nization phenomena realistically and theoretically.) What involved term to replace the role of eigenvalues of G(t) is the total length of all paths passing through an edge in the network connection graph inducing from G(t). Nevertheless, the method is limited to the networks with cooperative couplings, i.e., gij(t) ≥ 0, ∀i 6= j. As well, despite

the fact that the criterion based on the analysis of Lyapunov function guarantees the onset of synchronization, it is not a general method since there is no procedure for constructing the Lyapunov function for an arbitrary system. Similar problem occurs in the matrix measure approach. Furthermore, such criteria may not hold when the number of the coupled systems is large.

In the thesis, we shall present some criteria for the onset of synchronization in this model. Instead of constructing the Lyapunov function, the developed criteria are based on the completely different version of the matrix measure approach as proposed in [17,18]. The approach can overcome the drawbacks mentioned above.

1.2.2

Model II: Coupled Map Lattices

The model considered in this subsection is essentially similar to that given in Model I, except that the dynamics of units and the intersection ways are governed by some maps. It implies the equations of the motion then read:

xi(k + 1) = f (xi(k)) + d · m X

j=1

gijf(xj(k)), i = 1, . . . , m. (1.2a) Here x ∈ Rn and f is a vector-valued function from Rn_{× R to R}n. Based on the structure of the modeling, ones usually call (1.2a) as the coupled map lattices (CMLs). Set G = (gij). Then in the vector-matrix form, (1.2a) becomes

x_{(k + 1) = F (x(k)) + d(G ⊗ I)F (x(k))}, (1.2b)

where x(k) = (x1(k), . . . , xm(k))T, and F (x(k)) = (f (x1(k)), . . . , f (xm(k)))T.

This model, first introduced in 1980’s [20,56,93], has been the subject of much recent research. It is studied in the populations, chemical reactions, information pro-cessing, and biological networks and et. al [57]. Many dynamical behaviors have been observed, including spatiotemporal chaos and synchronization. In this thesis, we shall specifically focus on the issue of the synchronization.

The method developed to deal with the local stability problem for the synchro-nization state is based on the master-stability function criteria [1,19,28–30,42,106].

For the global stability problem, the methods include constructing the Lyapunov func-tion [64,65], applying the matrix measure criteria [61] and et. al. Nevertheless, compar-ing with Model I (1.1a), these methods are more limited. More precisely, the criteria limit to the special forms of the coupling matrix G.

In the thesis, we shall consider the local stability problem for the synchronization state in the model with coupling matrix G arbitrary given and study some correspond-ing problems.

1.2.3

Model III: Pulse Coupled Systems

The model introduced in this section is to explain the flashing synchrony in the fireflies [13,69], and the rhythmic activity of cells of the heart pacemaker [47,68,76,89], of cells of the pancreas [82] and of neural networks [12,13,25,76,80,88]. The most difference between this model and Models I, II is the interaction fashion. Interaction in Models I, II is “continuous in time”, while that in this model is “fleeting and intermittent in time”. Moreover, the reset mechanism occurs herein. Such a coupling fashion is called to be pulse-coupled, and the model is called the integrate-and-fire model. We start with introducing the Peskin’s model [76].

Let m be the number of the units in the coupled system, and the individual state be denoted by xi (i = 1, 2, · · · , m), where xi are subject to the dynamics

dxi

dt = −rixi+ Ii, 0 ≤ xi ≤ 1, (1.3)

with Ii > ri ≥ 0. As time progresses, suppose t−1 is the first time that some ith unit reaches the threshold 1, i.e., xi(t−) = 1 (one also says such unit fires). Then xi is reset, i.e.

xi(t−) = 1 → xi(t+) = 0, (1.4a)

and the firing effect from i to j (j 6= i) yields immediately: xj(t−) →



xj(t−) + gji if xj(t−) + gji < 1,

1 if xj(t−) + gji ≥ 1. (1.4b)

Suppose in this instant, the kth unit fires, i.e., xk(t−) + gki ≥ 1, then similar process like (1.4) occurs again at once. This process at time t is lasted until all effects of the firing units are evaluated. In addition, the whole evolution dynamics shall repeat above process once and once again. For the onset of synchronization in the integrate-and-fire model, it means that all units fire simultaneously after some fixed time T .

This model was later generalized by Mirollo and Strogatz [69]. It is assumed that the state variable xi evolves according to xi = fi(φi), where fi : [0, 1] → [0, 1] is smooth, strictly increasing, and satisfies fi(0) = 0 and fi(1) = 1, and the dynamics of phase φi is governed by dφi dt = 1 Ti , _{0 ≤ φ}i ≤ 1. (1.5)

Moreover, as time progresses, suppose t− is the first time that φi(t−) = 1 (corre-spondingly, xi(t−) = 1) for some i (one also says such unit fires). Then let Ej(t−) = f (Ti

Tj(1 − φi(t−)) + φj(t−)), j = 1, 2, · · · , m, and undergo the process of the reset and

firing effects as given in (1.4) for Ej(t−), j = 1, 2, · · · , m to get Ej(t+). Then the phase φi at time t+ is defined as φi(t+) = g(Ej(t+)). Here g is the inverse function of f . In addition, the whole evolution dynamics shall repeat above process once and once again.

As follows, we define some terminologies in the model. We say φ = (φ1, · · · , φm) is in the firing state if φi = 0 for some i. For a given initial phase φ(0), we say < ti > is its firing time series if it is the increasing sequence that records the successive time when phase φ(ti) =: φ(i) is in the firing state. Let r be the map defined in the set of firing states by r(φ(i)_{) = φ}(i+1) _{(Such definition is well-defined). Then synchronization} in the model is defined to satisfy that for any φ in the firing state, there is N ∈ N such that rn_{(φ) = (0, · · · , 0), ∀n ≥ N.}

Note that the Peskin’s model is one of the special model proposed by Mirollo and Strogatz with fi(φ) = Ii ri(1 − e −riTiφ_). and T = ln( Ii Ii−ri) r .

In the article [76], Peskin conjectured that, first, for identical units, the coupled system approaches synchronization for almost all initial conditions, and second, this remains true even when the units are not quite identical. Herein, “identical ” means fi ≡ f, Ti ≡ T , and gij ≡ g > 0. For the first conjecture, Mirollo and Strogatz [69] give a rigorous proof for the case that f′′

i < 0. The second part of Peskin’s conjecture was verified by Urbanczik and Senu [91] with flat units, i.e., f′′

i ≡ 0. However, Bottani [11] numerically showed that even concave-upward units, i.e., fi′′ > 0, can synchronize, provided that the concavity is not too large. In the article [15], Chang and Juang show that if the stability condition holds (see, Eq. (2.9) therein), then the nonidentical coupled system with f′′

i < 0 will achieve synchrony for almost all initial conditions; if both the stability condition (see, Eq. (2.9) therein) and absorption condition (see, Eq. (3.19) and (3.20) therein) hold, then the coupled system with f′′

i > 0 will achieve synchrony for almost all initial conditions. The holding of the stability condition implies a group of units reaching the threshold at the same time will remain coordinated in the future, while the holding of the absorption condition implies the number of firing units grows larger and larger and ultimately all units fire simultaneously. We comment that the stability condition requires that gij > 0.

1.3

Organization and results of the thesis

The organization of the thesis is as follows. In Chapter 2, we study the global synchro-nization in Model I (1.1), including the cases that coupling matrix G(t) therein is time independence and time dependence. Some criteria for the onset of synchronization are given. In Chapter 3, we take the coupled Lorenz equations, coupled Duffing equations, and Hindmarsh-Rose neurons as examples to see the applications of the given criteria. Specifically, we study furthermore the dynamics of the coupled Hindmarsh-Rose neu-rons under synchronization. The developing method for these criteria includes roughly the concepts of matrix measure, and coordinate transformation. In Chapter 4, we study the local synchronization in Model II with coupling matrix G generally given. Some criteria for the local stability of the synchronization state are given. In addition, the

phenomenon of wavelength bifurcations, as the terminology using in [37], and the ap-plication of wavelet transform method in Model II are discussed. Here, the wavelet transform method is an approach which is first given out concerning about Model I (see, e.g., [51,52,83,98]) in order to bring the synchronization easier to happen. It will be shown that this method also applies well in the CMLs (1.2). In Chapter 5, we summarize the results in this thesis and give some directions of the future work. We remark that most results in the thesis are adopted from [49,50,53–55].

Chapter 2

Synchronization in Model I

Before going into details of the derivation of synchronization for Model I, we first give some needed preliminaries.

2.1

Preliminaries

First, we introduce the concept of matrix measures. The following definitions and properties of matrix measures can be referenced to the book by M. Vidyasagar [92]. Definition 2.1.1. ( [92]) Let || · ||i be an induced matrix norm on Rn×n. The matrix

measure of matrix K on Rn×n _{is defined to be µ}

i(K) = lim

ǫ→0+

||I+ǫK||i−1

ǫ .

Lemma 2.1.1. ( [92]) Let || · ||k be an induced k-norm on Rn×n, where k = 1, 2, ∞.

Then each of matrix measure µk(K), k = 1, 2, ∞, of matrix K = (kij) on Rn×n is,

respectively, µ_∞(K) = max i {kii+ X j6=i |kij|}, µ1(K) = max j {kjj + X i6=j |kij|}, and µ2(K) = λmax(KT + K)/2.

Theorem 2.1.1. ( [92]) (i) µi(αA) = αµi(A), ∀ α ≥ 0 (ii) µi(A + B) ≤ µi(A) + µi(B). (iii) If λ is an eigenvalue of A, then Reλ ≤ µ(A). (iv) Consider the differential

equation ˙x_{(t) = K(t)x(t) + v(t), t ≥ 0, where x(t) ∈ R}n_{, K(t) ∈ R}n×n, and K(t), v(t)

are piecewise-continuous. Let || · ||i be a norm on Rn, and || · ||i, µi denote, respectively,

the corresponding induced norm and matrix measure on Rn×n_{. Then whenever t ≥ t}

0 ≥ 0, we have ||x(t)||i ≤ ||x(t0)||iexp Z t t0 µi(K(s))ds  + Z t t0 exp Z t s µi(K(τ ))dτ  ||v(s)||ids. (2.1) Next, we introduce a function being of type K, which generates a monotone dynamics of the system of linear differential equations. For completeness and ease of the references, we also recall the definitions of the above described concepts and their properties [41,92].

Let Rn

+ = {x = (x1, x2, . . . , xn)T ∈ Rn : xi ≥ 0, i = 1, . . . , n} be the nonnegative cone. Let a, b ∈ Rn_{. We write a ≤ b if b − a ∈ R}n

+.

Definition 2.1.2. We say that a function f = (f1, . . . , fn) : D ⊂ Rn → Rn is of type

K on D if, for each i, fi(a) ≤ fi(b) whenever a = (a1, . . . , an) and b = (b1, . . . , bn) are

in D with a ≤ b and ai = bi.

The following theorem amounts to saying that a vector field being of type K is a sufficient condition to generate a monotone dynamics.

Theorem 2.1.2. ( [41]) Let f (t, x) be of type K on Rn _{for each fixed t and let x(t)}

be a solution of ˙x(t) = f (t, x) on [a, b]. Let z(t) be continuous on [a, b] and satisfy

Dlz(t) ≤ f(t, z). Here Dlx(t) = lim

h→0−

x_(t+h)−x(t)

h . Then z(t) ≤ x(t) for a ≤ t ≤ b

provided that z(a) ≤ x(a).

Consider linear system of differential equations in the homogeneous case

y′ = A(t)y. (2.2)

Here A(t) is an n × n matrix. Then clearly if lim

t→∞

Rt

0 µ2(A(s))ds

t ≤ −r for some r > 0. Then y(t) converges to zero exponentially. The following propositions play one of the critical steps in obtaining our main results.

Proposition 2.1.1. Suppose ξ(t), η(t) and ζ(t) are nonnegative functions on [0, ∞) satisfying the following inequalities

Dlξ(t) ≤ a1(t)ξ(t) + a4(t)η(t) + a5(t)ζ(t), (2.3a) Dlη(t) ≤ a6(t)ξ(t) + a2(t)η(t) + a7(t)ζ(t), (2.3b) Dlζ(t) ≤ a8(t)ξ(t) + a9(t)η(t) + a3(t)ζ(t). (2.3c)

Here ai(t), i = 4, 5, . . . , 9, are nonnegative functions on [0, ∞). Then ξ(t), η(t), and

ζ(t) converge to zero exponentially provided that lim

t→∞ Rt 0µ2(A(s))ds t ≤ −r, for some r > 0, where A(t) =   aa16(t) a(t) a42(t) a(t) a57(t)(t) a8(t) a9(t) a3(t)   , (2.4)

(ii) Suppose, in addition, that, ai(t), i = 1, . . . , 9, are constants. Then ξ(t), η(t), and

ζ(t) converge to zero exponentially provided that all eigenvalues of A are negative. Proof. Let ξ(t), η(t) and ζ(t) satisfy the following equation.

˙ξ = a1(t)ξ + a4(t)η + a5(t)ζ, ξ(0) = ξ(0), ˙η = a6(t)ξ + a2(t)η + a7(t)ζ, η(0) = η(0), ˙ζ = a8(t)ξ + a9(t)η + a3(t)ζ, ζ(0) = ζ(0).

It is easily checked that the above system is of type K. Following from Theorem 2.1.2 that ξ(t) ≥ ξ(t), η(t) ≥ η(t), and ζ(t) ≥ ζ(t), for all t ≥ 0, we see that the first statement of the proposition holds as claimed. The second assertion of the proposition is obvious.

Proposition 2.1.2. Let ˙x_{= A(t)x. Here A(t) is an n×n matrix. Suppose lim}

t→∞A(t) =

A. Then x(t) converges to the origin exponentially provided that all real parts of

eigenvalues of A are negative.

Proof. For any ǫ > 0, there is a Pǫ such that A can be decomposed into a Jordan form

PǫAPǫ−1 = D + ǫQ,

where D is a diagonal matrix with the diagonal entries being all the eigenvalues of A, and Q is a matrix with its entries being either 0 or 1. Then PǫA(t)Pǫ−1 = Pǫ(A(t) − A)P−1

ǫ + ǫQ + D. It follows µ2(PǫA(t)Pǫ−1) ≤ µ2(Pǫ(A(t) − A)Pǫ−1) + ǫµ2(Q) + µ2(D) ≤ µ2(Pǫ(A(t) − A)Pǫ−1) + ǫn + µ2(D). Since lim

t→∞A(t) = A, we get

µ2(PǫA(t)Pǫ−1) ≤ (n + 1)ǫ + µ2(D), whenever t ≥ tǫ for some tǫ > 0. Hence,

lim

t→∞

Rt

0 µ2(PǫA(t)Pǫ−1)ds

t ≤ (n + 1)ǫ + µ2(D).

By the arbitrariness of ǫ, take ǫ = −µ2(D)

2(n+1). Then, we have lim t→∞ Rt 0 µ2(PǫA(t)Pǫ−1)ds t ≤ µ2(D) 2

Thus, all solutions of ˙y = (PǫA(t)Pǫ−1)y converges to the origin, and so are those of ˙

x= A(t)x.

Proposition 2.1.3. Let A and G be matrices of dimension m × m and n × n,

respec-tively, and Ip be the p × p identity matrix. Let λi, i = 1, · · · , k, be all the eigenvalues

of G. Then the real parts of the eigenvalues of (A ⊗ In) +  I1 0 0 0  ⊗ G 

are negative provided that all real parts of the eigenvalues of matrices Mi := A + λi  I1 0 0 0  are negative.

Proof. For any ǫ > 0, there is Pǫ such that

where D is a diagonal matrix with the diagonal entries being all the eigenvalues of G, and Q is a matrix with its entries being either 0 or 1. Then

(Im⊗ Pǫ)  (A ⊗ In) +  I1 0 0 0  ⊗ G  (Im⊗ Pǫ−1) =  (A ⊗ In) +  I1 0 0 0  ⊗ D  + ǫ  I1 0 0 0  ⊗ Q By taking ǫ sufficiently small, we get real parts of the eigenvalues of

(A ⊗ In) +  I1 0 0 0  ⊗ G  are negative iff those of

(A ⊗ In) +  I1 0 0 0  ⊗ D (2.5)

are negative. Then, the proof is completed by noting that after some permutation, matrix in (2.5) becomes diag(W1, · · · , Wn), where Wi = Mj, for some j = 1 . . . , k.

2.2

Global synchronization with time-invariant

cou-pling

In this section, we first consider synchronization of (1.1) with the time-invariant cou-pling. It means we consider the synchronization in the following system.

dxi dt = f (xi, t) + d · m X j=1 gijDxj, i = 1, 2, . . . , m, (2.6a) where xi = (xi1, xi2, . . . , xin)T ∈ Rn, and f is a vector-valued function from Rn× R to Rn _{denoted by f (x, t) = (f1(x, t), · · · , fn(x, t))}T_{. Or equivalently,} ˙x =    f(x1, t) .. . f(xm, t)    + d(G ⊗ D)x =: F (x, t) + d(G ⊗ D)x, (2.6b) where x = (x1, x2. . . , xm)T, and G = (gij).

As one usually concerns, synchronization is the phenomenon that units in a group have their dynamical behaviors get closer and closer as time progresses, and eventually they tend to be identical. So, mathematically, we define synchronization in the same sense. Mention that such definitions are also set for the time-varying coupling.

Definition 2.2.1. (i) The synchronous manifold M of Model I (1.1) is defined as the

set M = {x = (x1, · · · , xm)T : xi = xj, 1 ≤ i, j ≤ m}. (ii) Model I (1.1) is said to be

locally synchronized if the synchronous manifold M is asymptotically stable under the given coupling strength d.

Definition 2.2.2. Model I (1.1) is said to be globally synchronized if, under the given

coupling strength d, for all initial conditions xi(t0) (i = 1, 2, . . . , m) in Rn,

lim

t→∞kxi(t) − xj(t)k = 0, ∀1 ≤ i, j ≤ m.

We next give the definition of the bounded dissipation of a system.

Definition 2.2.3. Model I (1.1) is called to be bounded dissipative (with respect to α) if there is a bounded region Bmn(α) =: {x : kxk ≤ α} such that for each parameter d > 0, and each initial value x(0), there is a time t0, such that x(t) lies in Bmn(α)

whenever t ≥ t0.

To prove global synchronization of coupled chaotic systems, one needs to assume bounded dissipation, which plays the role of an a priori estimate. Without such an a priori estimate, as in the case of the R¨ossler system, global synchronization is much more difficult to obtain. Only local synchronization was reported numerically in litera-ture (see e.g., [75]). An interesting question in this direction is how bounded dissipation of the coupled system is related to the uncoupled dynamics and its connectivity topol-ogy. Not much general theorems have been provided so far. In just the case that G(t) and D are specially given, it was shown in [7] that bounded dissipation of the single oscillator implies that of the coupled oscillators. Moreover, the absorbing domain of the coupled system is a topological product of the absorbing domain of each individual system.

Now, we impose the conditions on coupling matrices G and D. We assume, throughout the section, that

(i)λ = 0 is a simple eigenvalue of G and e= _√1

m(1, 1, . . . , 1) T

1×m is its corresponding eigenvector; (2.7a) (ii) All nonzero eigenvalues of G have negative real part. (2.7b)

Such assumption above is to ensure the invariant property of the synchronous man-ifold M and make the dynamics of each unit under synchronization be the same as that without coupling. We further assume that coupling matrix D is, without loss of generality, of the form

D=  Ik 0 0 0  n×n. (2.7c) The index k, 1 ≤ k ≤ n, means that the first k components of the individual system are coupled. If k 6= n, then the system is said to be partial-state coupled. Otherwise, it is said to be full-state coupled.

To study synchronization of equation (2.6), we permute the state variables in the following way: ˜ xi =    x1i ... xmi    , and ˜x=    ˜ x1 ... ˜ xn    . (2.8)

Then (2.6b) can be written equivalently as

˙˜x =    ˜ f1( ˜x, t) .. . ˜ fn( ˜x, t)    + d(D ⊗ G) ˜x=: ˜F( ˜x_{, t) + d(D ⊗ G) ˜}x, (2.9a) where ˜ fi( ˜x, t) =    fi(x1, t) .. . fi(xm, t)    . (2.9b) The purpose of such a reformulation is two-fold. First, a transformation of coordinates of ˜x is to be applied to (2.9a) so as to isolate the synchronous manifold. Second, once the synchronous manifold is isolated, proving synchronization of (2.6), is then equiva-lent to showing that the origin is asymptotically stable with respect to reduced system (2.12). To do this, we first make a coordinate change to decompose the synchronous subspace. Let A be an m × m matrix of the form

A=         1 ₋₁ 0 _{· · ·} 0 0 . .. ... ... ... .. . . .. ... ... 0 0 _{· · ·} 0 1 ₋₁ 1 √ m · · · · 1 √ m 1 √ m         m×m =:  E eT  , (2.10a)

where e is given as in (2.7a). It is then easy to see that EET _{is invertible and that}

A−1 = ET(EET)−1, e
_. (2.10b) Setting b A= In⊗ A, (2.10c) we see that b A_{(D ⊗ G) b}A−1 = (In⊗ A)(D ⊗ G)(In⊗ A−1) = D ⊗ AGA−1 = D ⊗  EGET(EET₎−1 ₀ eT_GET_(EET₎−1 ₀  =: D ⊗ ¯_hGT 0₀  . (2.10d) We remark, via (2.10d), that σ( ¯G_{) = σ(G) − {0}, where σ(·) takes the spectrum of a} matrix. Multiplying bAto the both side of equation (2.9a), we get

˙˜y =: bA ˙x˜ = bA ˜F( ˜x, t) + d bA_{(D ⊗ G) b}A−1y˜ = bA ˜F( bA−1y˜_{, t) + d( D ⊗} ¯G 0 hT ₀  ) ˜y. (2.11) Let ˜y =    ˜ y1 ... ˜ yn    . Then ˜yi =        x1,i− x2,i ... x_{m−1,i− x}m,i 1 √_m m P j=1 xj,i        . Setting ˜yi =   ¯ yi 1 √ m m P j=1 xj,i   , and ¯ y=    ¯ y1 .. . ¯ yn    ,

we have that the dynamics of ¯y is satisfied by following equation

Here ¯F is obtained from bA ˜F( bA−1y, t) accordingly.˜

The task of obtaining global synchronization of system (2.6) is now reduced to showing that the origin is globally and asymptotically stable with respect to system (2.12). To this end, the space ¯y is broken into two parts ¯yc, the coupled space, and

¯

yu, the uncoupled space. ¯ y=  ¯ yc ¯ yu  , and ¯F( ¯y, t) = ¯F_¯c( ¯y, t) Fu( ¯y, t)  , respectively. (2.13) Here ¯yc =    ¯ y1 .. . ¯ yk    , and ¯yu =    ¯ yk+1 .. . ¯ yn  

. The dynamics on the coupled space with respect to the linear part is under the influence of ¯G, which is asymptotically stable. The dynamics of the nonlinear part on coupled space can then be controlled by choosing large coupling strength. As a matter of fact, it is easier to obtain synchronization of coupled chaotic systems with a larger coupled space. On the other hand, the uncoupled space has no stable matrix ¯Gto play with. Thus, its corresponding vector field ¯Fu( ¯y, t) must have a certain structure to make the trajectory stay closer to the origin as time progresses. As we shall explain latter.

Now, assume that ¯Fc( ¯y, t) satisfies a uniformly Lipschitz condition with a uni-formly Lipschitz constant b1. That is,

k ¯Fc( ¯y, t)k ≤ b1k ¯yk (2.14a)

whenever ¯y in the ball B(m−1)n(α), and for all time t. Since the estimate in the

right-hand side of (2.14a) depends on the whole space ¯y, condition (2.14a) is a mild as-sumption provided that the coupled system is bounded dissipative. Write ¯Fu( ¯y, t) as

¯

Fu( ¯y, t) = U (t) ¯yu+ ( ¯Fu( ¯y, t) − U(t) ¯yu)

=: U (t) ¯yu+ ¯Ru( ¯y, t). (2.14b)

(i) U (t) is a block diagonal matrix of the form U (t) = diag(U1(t), · · · , Ul(t)) where each Uj(t), j = 1, . . . , l, are matrices of size (m − 1)kj× (m − 1)kj. Here

l X

j=1 kj = n − k, and kj ∈ N. Assume that there exists a constant γ > 0 such that matrix

measures µi(Uj(t)) ≤ −γ, for all t and all j. (2.14c)

(ii) Let ¯Ru( ¯y, t) =    Ru1( ¯y, t) ... Rul( ¯y, t)  

. Here l is the number given in (i). Then Ruj( ¯y, t), j = 1, . . . , l, satisfy a strong uniformly Lipschitz condition with a strong uniformly Lipschitz constant b2. Specifically, let ¯yu =

   ¯ yu1 .. . ¯ yul  

, written in accordance with the block structure of U (t). Then we assume that

kRuj( ¯y, t)k ≤ b2k      ¯ yc ¯ yu1 ... ¯ y_{u j−1}     k (2.14d)

whenever ¯y in the ball B_(m−1)n(α), and for all j = 1, . . . , l and all time t.

Specifically, we break the vector field ¯Fuinto (time dependent) linear part U (t) ¯yu and nonlinear part ¯Ru( ¯y, t). We will further break U (t) into certain block diagonal form if necessary. Note that form (2.14b) can always be achieved since the remainder term ¯Ru still depends on the whole space ¯y. To take control of the dynamics on the linear part, we assume that the matrix measure of each diagonal block Uj(t) is negative. As to contain corresponding dynamics on the nonlinear part, we assume that (2.14d) holds. Note that though the nonlinear terms Ruj( ¯y, t) could possibly depend on the whole space, their norm estimates are required to depend only on the coupled space and uncoupled subspaces with their indices proceeding j. In this set up, the nonlinear dynamics on uncoupled space can be iteratively controlled by choosing large

coupling strength. We also remark that if (2.14c) and (2.14d) are satisfied for l, the number of diagonal blocks, being one, then we do not need to further break U (t). Such further breaking is needed only if (2.14c) and (2.14d) are not satisfied. The proof in the following theorem gives exactly how the above strategy can be realized.

Theorem 2.2.1. Let G and D be given as in (2.7). Assume that ¯F satisfies (2.14),

and system (2.12) is bounded dissipative with respect to α. Let λ1 = max{λj|λj∈ Re(σ( ¯G))}.

If d > cb1 −λ1− ǫ  1 + (b2 γ) 2 l 2 =: dmin, (2.15)

where ǫ ≥ 0 and c is some constant depending on G and ǫ, then lim

t→∞y(t) = 0.¯

Proof. Since system (2.12) is bounded dissipative with respect to α, without loss of generality, we may assume that k ¯y_{(t)k ≤ α for all time t ≥ t}0. Using (2.14b), we write (2.12) as  _˙¯y c ˙¯yu  =  d(Ik⊗ ¯G) 0 0 U(t)   ¯ yc ¯ yu  + ¯F_¯c( ¯y, t) Ru( ¯y, t)  . (2.16a) Applying the variation of constant formula to (2.16a) on ¯yc, we get

¯ yc(t) = e(t−t0)d(Ik⊗ ¯ G₎ ¯ yc(t0) + Z t t0 e(t−s)d(Ik⊗ ¯G)_F¯ c( ¯y(s), s)ds. (2.16b)

Let λ1 = max{ λj|λj ∈ Re( σ(G) − {0} ) }. Then

ketd(Ik⊗ ¯G)_{k ≤ ce}tdν _(2.16c)

for ν = λ1+ ǫ and some constant c. Here 0 < ǫ < −λ1. Thus, k ¯yc(t)k ≤ ce(t−t0)dνk ¯yc(t0)k + cb1 Z t t0 ed(t−s)ν_{k ¯}y_(s)kds ≤ ce(t−t0)dν_{α +}α d cb1 |ν| =: ce (t−t0)dν_{α +} α dc0. (2.16d)

Let δ > 1, we see that

k ¯yc(t)k ≤ α

whenever t ≥ t0,1 for some t0,1 > 0. We then apply Theorem 2.1.1 on ¯yu1 and the resulting inequality is k ¯yu1(t)k ≤ k ¯yu1(t0,1)k exp (Z t t0,1 µi(U1(s))ds ) + Z t t0,1 exp Z t s µi(U1(τ ))dτ  kRu1( ¯y(s), s)kds. It then follows from (2.14c), (2.14d), and (2.17a) that

k ¯yu1(t)k ≤ αe−γ(t−t0,1)+ α d b2 γc0δ ≤ α d b2 γc0δ 2 _=: α dc1δ 2 , (2.17b)

whenever t ≥ t1,1 for some t1,1 ≥ t0,1. Inductively, we get k ¯yuj(t)k ≤ α d  b2 γ v u u tXj−1 i=0 c2 i   δj+1 _=: α dcjδ j+1_{, j = 2, . . . , l, (2.17c)}

whenever t ≥ tj,1(≥ tj−1,1). Letting tl,1 = t1 and summing up (2.17a)-(2.17c), we get k ¯y_{(t)k =} v u u tXl j=1 k ¯yuj(t)k2+ k ¯yc(t)k2 ≤ α d  1 + (b2 γ) 2 l 2 _cb 1 |ν|δ l+1 _{=: hα,} whenever t ≥ t1. Choosing d >  1 + (b2 γ)2 l 2 _cb 1

|ν|δl+1, we see that the contraction factor h is strictly less than 1, and k ¯y_{(t)k contracts as time progresses. To complete the proof} of the theorem, we note that δ > 1 can be made arbitrary close to 1. Consequently, if d >1 + (b2 γ) 2 l 2 _cb 1

|ν|, then h can still be made to be less than 1.

Remark 2.2.1. (i) In case that ¯G is symmetric, then c and ǫ can be chosen to be 1

and 0, respectively. (ii) b1 and b2 could possibly depend on α.

Corollary 2.2.1. Suppose ¯F and G are given as in Theorem 2.2.1. Let

D =  _˜ D_k×k 0 0 0  n×n, where Re( σ( ˜D) ) > 0. (2.18)

Assume, in addition, that either σ(G) or σ( ˜D) has no complex eigenvalue. Then

assertions in Theorem 2.2.1 still hold true, except dmin needs to be replaced by

dmin= c b1 (−λ1− ǫ) · min{Re( σ( ˜D) )}  1 + (b2 γ ) 2 l 2 . (2.19)

Proof. Assumption on D is to ensure that (2.16c) is still valid. Other parts of the proof are similar to those in Theorem 2.2.1 and are thus omitted.

We next turn our attention to finding conditions on the nonlinearities fi(u, t), i = 1, . . . , n, u ∈ Rn_{, so that assumptions (2.14) are satisfied. To this end, we need} the following notations. Let ˜xi and ˜x be given as in (2.8). Define

[ ˜xi]− =    x1,i ... x_m−1,i    , and [ ˜x]−=    [ ˜x1]− ... [ ˜xn]−    . (2.20)

We then break ˜F as given in (2.9a) into two parts so that the breaking is in consistent with ¯y in (2.13). Specifically, we shall write

˜ F( ˜x, t) = ˜F_˜c( ˜x, t) Fu( ˜x, t)  . (2.21) We are now in the position to state the following propositions.

Proposition 2.2.1. Suppose that fi(x, t), i = 1, 2, . . . , k satisfy a Lipschitz condition in Bn(α₂) with a Lipschitz constant b1. That is

|fi(u, t) − fi(v, t)| ≤ b1

k ku − vk,i = 1, 2, . . . , k, (2.22)

for all u, v in Bn(α₂) and all time t. Then (2.14a) holds true.

Proof. Note that bA ˜F( ˜x, t) =

   A ˜f1( ˜x, t) ... A ˜fn( ˜x, t)    ,

where A is given as in (2.10a), and so

[A ˜fi( ˜x, t)]−=    fi(x1, t) − fi(x2, t) ... fi(xm−1, t) − fi(xm, t)    , i = 1, 2, . . . , n. (2.23) Since ¯ Fc( ¯y, t) =    [A ˜f1( ˜x, t)]− ... [A ˜fk( ˜x, t)]−    , we conclude that (2.14a) holds.

From the above proposition, we see that the nonlinearities on the corresponding coupled space are only assumed to be Lipchitz. The following proposition is very useful in the sense that by checking how each component fi of the nonlinearity f is formed, one would then be able to conclude whether (2.14c) and (2.14d) are satisfied.

Proposition 2.2.2. Let u = (u1, . . . , un)T and v = (v1, . . . , vn)T be vectors in Bn(α₂).

Let wp =

p X

i=0

ki, p = 1, . . . , l, where k0 = k, the dimension of coupled space, and

k1, . . . , kl and l are given as in (2.14c). Write fwp−1+i(u, t)−fwp−1+i(v, t), i = 1, . . . , kp,

as fwp−1+i(u, t) − fwp−1+i(v, t) = kp X j=1

qwp−1+i,wp−1+j(u, v, t)(uwp−1+j − vwp−1+j) + rwp−1+i(u, v, t).

(2.24a) We further assume that the followings are true.

(i) For p = 1, . . . , l, let Qu,v,p = (qwp−1+i,wp−1+j(u, v, t)), where 1 ≤ i, j ≤ kp.

Then µ_∗(Vp) < −γ for all p, u, v in Bn(α₂) and all time t, where ∗ = 1, 2, ∞.

(2.24b)

(ii) Let rp = rwp−1+1(u, v, t), . . . , rwp(u, v, t)

T . We have that krpk ≤ b2k    u1− v1 .. . uwp−1 − vwp−1    k (2.24c)

for all p, u, v in Bn(α₂) and all time t.

Then (2.14c) and (2.14d) hold true for ∗ = 1, 2, ∞.

Proof. Since ri(u, v, t) depend on the whole space, fi(u, t) − fi(v, t) can always be

written as the form in (2.24a). Using (2.24a) and (2.23), we have that the matrices Up(t) in the linear part of ¯Fu( ¯y, t) take the form

Up(t) =

m−1_X

w=1

where xw are given as in (2.6a), and (Dw)ij =



1 i = j = w,

0 otherwise, 1 ≤ i, j ≤ m − 1.

It then follows from Lemma 2.1.1, and (2.25) that µ_∗(Up(t)) < −γ for ∗ = 1 or ∞. For ∗ = 2, we have that m−1_[ w=1 σ{Qx_w,xw+1,p(t) + Qxw,xw+1,p(t)
T } = σ (_m−1 X w=1  Qx_w,xw+1,p(t) ⊗ D_w+ Qx_w,xw+1,p(t)
T ⊗ Dw ) = σ Up(t) + UpT(t)
,

where σ(A) is the spectrum of A. We remark that the first equality above can be verified by the definition of eigenvalues due to the structure of Up(t). It then follows from Lemma 2.1.1 that µ2(Up(t)) < −γ. The remainder of the proof is similar to that of Proposition 2.2.1, and is thus omitted.

Remark 2.2.2. The upshot of Proposition 2.2.2 is that by only checking the “structure” of the vector field f of the single oscillator, one should be able to determine if our main result can be applied. To be precise, we begin with saving notations by setting f

as f = f (x, t) = (f1(x, t), . . . , fn(x, t))T. We then check the form of the difference

of “uncoupled” part of dynamics. That is, we write fi(u, t) − fi(v, t) in the form of

(2.24a) with i = k + 1, . . . , n. If (2.24b)-(2.24c) can be satisfied, then l = 1 gets the job done. Otherwise, we further break the uncoupled states into a set of smaller pieces to see if the resulting (2.24b)-(2.24c) are satisfied.

We are now ready to state the main theorems of the paper.

Theorem 2.2.2. Assume that system (2.6) is bounded dissipative. Let coupling

matri-ces G and D satisfy (2.7) and the nonlinearities fi(x, t), i = 1, 2, . . . , n, satisfy (2.22)

and (2.24). Suppose d is greater than dmin, as given in (2.15). Then system (2.6) is

Proof. The proof is direct consequences of Propositions 2.2.1 and 2.2.2, and Theorem 2.2.1.

Theorem 2.2.3. Coupled system (2.6) with D given as in Corollary 2.2.1, is globally synchronized provided that the coupled system is bounded dissipative, the nonlinearities fi(x, t), i = 1, 2, . . . , n, satisfy (2.22) and (2.24), and d is greater than dmin. Here dmin is given in (2.19).

The actual way to apply the above Theorem 2.2.1, 2.2.2 to the real question is postponed to the Chapter 3.

2.3

Global synchronization with time-varying

cou-pling

In this section, we consider the synchronization of Model (1.1). As given in (2.9a), we can write the coupled system equivalently as

˙˜x =    ˜ f1( ˜x, t) ... ˜ fn( ˜x, t)    + d(D ⊗ G(t)) ˜x=: ˜F( ˜x_{, t) + d(D ⊗ G(t)) ˜}x, (2.26) where ˜xi, ˜x and ˜fi( ˜x, t) are defined as in (2.8), (2.9b). Our basic strategy to get the criteria of synchronization such as Theorem 2.2.2 is similar to that gotten there. What makes the most difference is that in this case Eq. (2.16b) does not hold again. To deal with it, we will apply Theorem 2.1.1. Nevertheless, in Theorem 2.1.1, what is concerned is the matrix measure of a matrix not the eigenvalue of a matrix. Thus, one needs to deal with the corresponding problem carefully. To do so, we make the usage of the concept of “coordinate transformation”.

First, instead of defining E as in (2.10a), herein we let E be an (m − 1) × m full-rank matrix with all its row sums being zero. Such a matrix is to be termed a coordinate transformation. Define

A=  E eT  . (2.27a)

Then A−1 _{= E}T_(EET₎−1_{, e
and} AG(t)A−1 =  EG(t)ET_(EET₎−1 ₀ eT_G(t)ET_(EET₎−1 ₀  =: ¯GE(t) 0 h(t)T ₀  . (2.27b) Let bA= In⊗ A and ˜y= bAx. Multiplying b˜ Ato both sides of Equation (2.26), we get

˙˜y = bA ˜F( bA−1y, t) + d˜  D_⊗ ¯GE(t) 0 h(t)T ₀  ˜ y. Let ˜y= ( ˜y1, . . . , ˜yn)T. Then ˜ yi =  Ex˜i Pm j=1xji/√m  =:  ¯ yi ei  . (2.28) Setting ¯y = ( ¯y1, . . . , ¯yn)T, we have that the dynamics of ¯y is now satisfied by the following equation

˙¯y = d(D ⊗ ¯GE(t)) ¯y+ ¯F( ¯y, t), (2.29a) where

¯

F( ¯y, t) = (In⊗ E) · ˜F( bA−1y, t).˜ (2.29b) Since the rank and the row sums of E are m − 1 and 0, respectively, we conclude that the task of obtaining global synchronization of system (1.1) is now reduced to showing that the origin is globally and asymptotically stable with respect to system (2.29a). The choice of a coordination transformation will greatly influence how negative the matrix measure of ¯GE(t) could be, which plays the important role, among others, to determine the global stability of (2.29a) with respect to the origin.

2.3.1

Matrices of the Coordinate Transformation

In what follows we shall address the question of how to choose a matrix E of the coordinate transformation, and its corresponding properties. To make the origin an asymptotically stable equilibrium of system (2.29a), one would like to have the matrix measure of ¯GE(t) as smaller a negative number as possible. In fact, such an optimal choice E can be achieved provided that the outer coupling matrix G(t) is symmetric, nonpositive definite.

Definition 2.3.1. Denote by C the set of (m − 1) × m coordinate transformations, i.e., C_{= {E ∈ R}(m−1)×m _{: E is full-rank, and all its row sums are zero}.}

Let O ⊆ C be such that

O_{= {E ∈ C : E such that matrix A = (E}T_{, e)}T _{is orthogonal}.}

Theorem 2.3.1. Assume that all eigenvalues of outer coupling matrix G(t) have non-positive real parts. Then inf

E_∈Cµ2( ¯GE(t)) ≥ Re λ2(G(t)). Here Re λ2(G(t)) is the second largest real part of eigenvalues of G(t). If, in addition, G(t) is symmetric for all t, then the above equality can be achieved by choosing any E in O.

Proof. It follows from (2.27b) that the spectrum σ( ¯GE(t)) of ¯GE(t) is equal to σ(G(t))−

{0}. Using the fact that Re λ(K) ≤ λmax(K+K

T

2 ) for any real matrix K, we have, via Lemma 2.1.1, that µ2( ¯GE(t)) ≥ Re λ₂(G(t)). In particular, if E ∈ O and G(t) is symmetric, then ¯GE(t) = EG(t)ET

is symmetric and h(t) = 0. Here h(t) is given as in (2.27b). Therefore, µ2( ¯GE(t)) = λ₂(G(t)). We have just completed the proof of the theorem.

The theorem above amounts to saying that if G(t) is symmetric, nonpositive definite, then any choice of E in O yields the smallest possible matrix measure of

¯

GE(t). This, in turn, gives one the best possible position to study the stability of equation (2.29a) with respect to the origin.

Remark 2.3.1. In those earlier papers (see, e.g., [17,53,104]), the choice of the coor-dinate transformations is either

E1 =      1 −1 0 _{· · ·} 0 1 0 ₋₁ . .. ... ... ... ... ... 0 1 0 _{· · ·} 0 ₋₁      or E2 =      1 −1 0 _{· · ·} 0 0 1 ₋₁ . .. ... ... ... ... ... 0 0 · · · 0 1 ₋₁      . (2.30)

with periodic boundary conditions, i.e., G_{(t) ≡}          −2 1 0 _{· · ·} 0 1 1 ₋₂ 1 0 _{· · ·} 0 0 . .. ... ... ... ... ... ... ... ... ... 0 0 _{· · ·} 0 1 ₋₂ 1 1 0 _{· · ·} 0 1 ₋₂          m×m , (2.31)

the corresponding matrix measure of ¯GE_i, i = 1, 2 is positive whenever m > 7 (see

Table 2.1), while µ2( ¯GE) = λ₂(G) < 0 for all E ∈ O regardless the size of G.

m 4 5 6 7 8 9

E1 −1.78 −1 −0.51 −0.19 0.05 0.23

E2 −1.78 −1 −0.51 −0.19 0.05 0.23

Table 2.1: The table gives the matrix measures of ¯GE_i, i = 1, 2, with various size of G, which is given in (2.31). Since G is a circular matrix, the matrix measures of ¯G with respect to E1 and E2 are equal. Note that the matrix measure of ¯GE is λ2(G), ∀E ∈ O, which is negative regardless the size of G.

Theorem 2.3.2. For any outer coupling matrix G(t), and any coordinate

transforma-tions Ep, Eq in O, µ2( ¯GEp(t)) = µ2( ¯GEq(t)).

Proof. Since for any x ∈ Rm−1_{, there is z = E}

qETpx such that

xTEp(G(t) + G(t)T)EpTx= zTEq(G(t) + G(t)T)EqTz. By the definition of matrix measure, we have that µ2( ¯GEp(t)) = µ2( ¯GEq(t)).

In this section, matrix E in Equation (2.29a) is assumed to lie in D unless oth-erwise stated. For ease of the notations, we shall drop the subscript E of ¯GE(t) if E _{∈ O. The remainder of the subsection is devoted to finding the matrix measure of}

¯

G(t) where its corresponding coupling matrix G(t) appears often in many applications. Proposition 2.3.1. Assume that for each t, G(t) is a node-balancing matrix, i.e., its row sums and column sums are equal. Then

µ2( ¯G(t)) = λ2(

G(t) + G(t)T

whenever all eigenvalues of G(t) + G(t)T _{are nonpositive.} Proof. If G(t) is as assumed, then it follows from (2.27b) that

AG(t)A−1 = ¯G(t) 0

0 0

 . Consequently, (2.32) holds as asserted.

In what follows, some outer coupling matrices are to be provided. Their corre-sponding matrix measures of ¯G(t) and ¯GE_i(t), i = 1 or 2, are to be compared.

Example 1. ( [2]) Consider the regular coupled network by adding to the pristine world G (the ring of 2K-nearest coupled oscillators) an additional global coupling such that the coupling p(t), 0 ≤ p(t) ≤ 1 is placed on all free spots of the matrix G (see, e.g., [2]). Specifically, the resulting coupling matrix G(t) can be represented by a circular matrix of the form

G_{(t) = circ(−g(t),} K z }| { 1, . . . , 1, m−2K−1 z }| { p(t), . . . , p(t), K z }| { 1, . . . , 1), (2.33)

where g(t) = 2K + (m − 2K − 1)p(t). Since G(t) is symmetric, we have that µ2( ¯G(t)) = λ2(G(t)) = max 1≤j≤m−1 −g(t) + K X l=1 (ωlj + ω(m−l)j) + p(t) m−K−1_X l=K+1 ωlj ! .

Here ω = exp(2πi/m). The matrix measures µ2( ¯G(t)) and µ2( ¯GE_i(t)), i = 1, 2, with m = 8, K = 2, and p(t) = t, t ∈ [0, 1] are recorded in Figure 2.1.

Example 2. ( [52,98]) Let G = G(m)_{β , 0 ≤ β ≤ 1 be the diffusive matrix of size m × m} with mixed boundary conditions. That is, if m > 2,

G(m)_β =          −1 − β 1 0 _{· · ·} 0 β 1 ₋₂ 1 0 _{· · ·} 0 0 . .. ... ... ... ... ... . .. ... ... ... 0 0 _{· · ·} 0 1 ₋₂ 1 β 0 _{· · ·} 0 1 _{−1 − β}          m×m, (2.34)

0 0.2 0.4 0.6 0.8 1 -8 -7 -6 -5 -4 -3 -2 -1 µ2( ¯GEi(t)) µ2( ¯G(t)) m a tr ix m ea su re t

Figure 2.1: The matrix measures of ¯G(t) and ¯GE_i(t), i = 1, 2, with G being given in (2.33), m = 8, K = 2, and p(t) = t, are, respectively, represented by the solid line and the dotted lines above. Lines for ¯GE_i(t), i = 1, 2 are coincided since G(t) is circular for all t. and if m = 2, G(2)_β =  −1 − β 1 + β 1 + β _{−1 − β}  .

For such G, µ2( ¯G) = λ2(G) < 0. However, λ2(G) would move closer to the origin as the number of nodes increases. As a result, synchronization of the network is more difficult to be realized as the number m of nodes increases. In [52,98], a wavelet transformation method is proposed to alter the connectivity topology of the network. In doing so, λ2(G(t)) = λ2(p(t)) becomes a quantity depending on wavelet parameter p(t). By choosing suitable p(t), which is a wavelet transformation method [52,98] applied to the coupling matrix G(m)_β , one would expect that λ2(p(t)) will move away from the origin regardless the number of the nodes. Under such a reconstruction, the resulting coupling matrix G(t) is of the following form

G(t) = G(m)_β + p(t)(G(

m k)

β ⊗ ¯ee¯T), (2.35)

where ¯e= (1, . . . , 1)T _{∈ R}k_{. Here we assume p(t) ≥ 0 and k = 2}l _{for some l ∈ N, and} m = Nk for some N ∈ N − {1}. Since the reconstructed matrix G(t) is symmetric,

µ2( ¯G(t)) = λ2(G(t)) < 0. The matrix measures µ2( ¯G(t)) and µ2( ¯GE_i(t)), i = 1, 2, with m = 8, β = 1_{2, l = 1, and p(t) = t, t ∈ [0, 1] are recorded in Figure 2.2.}

0 0.2 0.4 0.6 0.8 1 -2.5 -2 -1.5 -1 -0.5 0 µ2( ¯G_E1(t)) µ2( ¯G_E2(t)) µ2( ¯G(t)) m a tr ix m ea su re t

Figure 2.2: The matrix measures of ¯G(t) and ¯GE_i(t), i = 1, 2, with G given in (2.35), m = 8, β = 1₂, l = 1, and p(t) = t, are, respectively, represented by the solid line and the dotted lines above.

Example 3. Let G(t) = circ(

m

z }| {

−2, 2, 0, . . . , 0), a circulant matrix. Since G(t) is a node-balancing matrix, µ2( ¯G(t)) = λ2(G(t)) < 0. Note that the values of µ2( ¯GE_i), i = 1, 2, are positive provided that m > 5 (see Table 2.2).

m 4 5 6 7 8 9

E1 −0.83 −0.17 0.24 0.54 0.78 0.98

E2 −0.83 −0.17 0.24 0.54 0.78 0.98

Table 2.2: The table gives the matrix measures of ¯GE_i, i = 1, 2, with various size of G, which is given in Example 3.

Proposition 2.3.2. Let E = (e1, . . . , em−1)T ∈ O. If, in addition, {ei}m−1i=1 are

diagonal matrix. Moreover,

µ2( ¯G(t)) = λ2(G(t)), (2.36)

whenever all eigenvalues of G(t) are nonpositive.

Proof. Note that ¯G(t) = EG(t)ET _{= (e}T

i G(t)ej). Hence, ¯G(t) is a diagonal matrix. Therefore, the assertion in (2.36) holds as asserted.

Example 4. ( [17]) Let G(t) describe a star-typed coupled network of the form

G(t) =      −d1(t) d1(t) . .. ... −d1(t) d1(t) 1 _{· · ·} 1 _{−(m − 1)}      m×m. (2.37)

Here d1(t) is a real number. We next show that a set {ei}m−1i=1 of column vectors can be chosen so that E = (ei, . . . , em−1) ∈ O and that {ei}m−1i=1 are pairwise G(t)-conjugate. Define αi = (i(i + 1))−1/2, i = 1, . . . , m − 1. Let

eT_i = ( i z }| { αi, . . . , αi, −iαi, m−i−1 z }| { 0, . . . , 0)

for all i = 1, . . . m − 1. Then ei, i = 1, . . . , m − 1 are orthonormal vectors. Moreover, they are also G(t)-conjugate. To see this, we first note that d1(t) is an eigenvalue of G(t) and its associated eigenvectors are ei, i = 1, . . . , m − 2. Therefore, eTiG(t)ej = 0 for all 1 ≤ i 6= j ≤ m − 2. Some direct computation would yield that eTi G(t)em−1 = 0 for i = 1, . . . , m − 2 and that eT

m−1G(t)em−1 = −d1(t) − (m − 1). By Proposition 2.3.2, we have that

µ2( ¯GE(t)) = max{−d1(t), −d1(t) − (m − 1)} = −d1(t). (2.38) The matrix measures µ2( ¯G(t)) and µ2( ¯GE_i(t)), i = 1, 2, with m = 8 and d₁(t) = t, t ∈ [0, 1] are demonstrated in Figure 2.3.

The remainder of the subsection is to address the system with even more complex topology.

0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 1.5 2 2.5 µ2( ¯G_E1(t)) µ2( ¯G_E2(t)) µ2( ¯G(t)) m a tr ix m ea su re t

Figure 2.3: The matrix measures of ¯G(t) and ¯GE_i(t), i = 1, 2, with G being given in (2.37), m = 8, and d1(t) = t, are, respectively, represented by the solid line and the dotted lines above.

Proposition 2.3.3. Let G(t) = O(t) + P (t) with O(t) and P (t) having all its row sums zero. Suppose further that P (t) is node-balancing. Then

µ2( ¯G(t)) ≤ µ2( ¯O(t)) + λ2(

P(t) + P (t)T

2 ),

whenever all eigenvalues of P (t) + P (t)T _{are nonpositive.}

Proof. Noting that ¯G(t) = EG(t)ET _{= ¯O(t) + EP (t)E}T_{, we easily conclude that the}

above inequality holds as asserted.

Example 5. ( [2]) Consider the outer coupling matrix G(t) to be of the random type. Specifically, G(t) is of the form:

G_{(t) = circ(−2K,} K z }| { 1, . . . , 1, m−2K−1 z }| { 0, . . . , 0, K z }| { 1, . . . , 1) + P (t) =: O + P (t), (2.39) where P (t) =: (pij(t)) is a symmetric matrix with all its row sums being zero, and satisfies pij(t) ≡ 0 for (i, j) with i − j mod m ≤ K or j − i mod m ≤ K, and pij(t) =

Sij(q) for (q − 1)τ ≤ t < qτ for all remaining pairs (i, j) with i 6= j. Here each of Sij(q) is a random variable that takes the value 1 with probability p and 0 with probability1 − p.

The random variables Sij(q) are assumed to be all independent. To each realiza-tion ω of this stochastic process S(1), S(2), . . ., where S(q) = {Sij(q) : i − j mod m ≤ K or j − i mod m ≤ K}, i.e., to each switching sequence ω, there corresponds a time-varying system described by Equation (1.1b).

Since P (t) is symmetric, by Proposition 2.3.3,

µ2( ¯G(t)) ≤ µ2( ¯O) + λ2(P (t)) ≤ µ2( ¯O) = λ2(O) < 0.

Let G(t) ≡ G. Generally speaking, infE¯_∈Cµ2( ¯GE¯) 6= µ2( ¯GE) for any E ∈ O. Nevertheless, µ2(GE) produces a good upper bound of infE¯_∈Cµ2( ¯GE¯).

To support the observation, we conclude this section by providing some additional network topologies where the matrix measure of its corresponding ¯GE(t), E ∈ O is smaller than that of ¯GE_i, i = 1, 2. As a matter of fact, µ₂( ¯GE_i), i = 1, 2, switch signs as the number of nodes increases. In contrast, µ2( ¯GE) mostly remains negative as the size of the system grows.

Example 6. Consider a generalized wheel-typed coupled network of the form as illus-trated in Figure 2.4(a). The inner nodes have the strong all-to-all connections. The outer nodes are only directly connected with their nearest neighbors. The communi-cations between the inner and outer nodes are through one way going from each inside node to its nearest outside node. Specifically, such a network can be written as the following. G_{(t) ≡}  G1 G2 G3 G4  m×m, (2.40)

where G1 =      −(m2 − 1) 1 · · · 1 1 . .. ... ... .. . . .. ... 1 1 _{· · ·} 1 ₋₍m_{2 − 1)}      m 2× m 2

corresponding to the all-to-all coupling, G2 = 0, G3 = 0.1I, and G4 = G(

m 2) 1 − 0.1I. Here G( m 2)

1 is the diffusive matrix with periodic boundary conditions and of size m2 × m

2.

The numerical computation suggests that the matrix measures of ¯GE_i, i = 1, 2, are positive provided that m ≥ 4 while that of ¯GE, E ∈ O remains negative (see Table 2.3).

m 4 6 8 10 5000

E1 0.11 0.32 0.53 0.74 517.47

E2 0.23 0.56 0.96 1.44 34843.01

E _{−0.1 −0.1 −0.1 −0.1 −0.1}

Table 2.3: The table gives the matrix measures of ¯GE_i, i = 1, 2, and ¯GE, E ∈ D with various size of G, which is given in (2.40).

Example 7. Consider the prism-typed coupled network of the form as illustrated in Figure 2.4(b). The difference between the generalized wheel-typed network and the one considered here lies only on how the inner nodes communicate with each other (see Figure 2.4). Specifically, such a network can be written as the following.

G_{(t) ≡}  G1 G2 G3 G4  m×m, where G1 = G( m 2) 1 , G2 = 0, G3 = 0.1I, and G4 = G( m 2)

1 − 0.1I. The numerical computation suggests (see Table 2.4) that the matrix measures of ¯GE_i, i = 1, 2, are positive provided that m ≥ 4, while that of ¯GE, E ∈ O stays negative until m = 86. The example demonstrates that a coordinate transformation E, E ∈ O, is indeed a good candidate among all coordinate transformations.

m 4 6 8 86 88

E1 0.34 0.32 0.35 4.65 4.72

E2 0.34 0.56 0.73 4.79 4.86

E _{−0.1 −0.1 −0.1 −0.0006 0.0004}

Table 2.4: The table gives the matrix measures of ¯GE_i, i = 1, 2, and ¯GE, E ∈ D with various size of G, which is given in Example 7.

1 2 N 3 N-1 N+1 N+2 N+i 2N 2N-1 N+3 i (a) 1 2 N 3 N-1 N+1 N+2 N+i 2N 2N-1 N+3 i (b)

Figure 2.4: Coupling Topologies: (a) Generalized wheel-typed coupled network with m = 2N, and (b) Prism-typed coupled network with m = 2N. Networks (a) and (b) appear in Examples 6 and 7, respectively.

2.3.2

Synchronization Criteria

In the section, we turn our attention back to the dynamics of (2.29a), and analyze the stability of the origin of the system. Let ¯y, ¯yc, and ¯yu be defined as in (2.13). Then, like (2.16a) , Equation (2.29a) can be rewritten as in the form

 ˙¯yc ˙¯yu  =  d(Ik⊗ ¯G(t)) 0 0 U(t)   ¯ yc ¯ yu  + ¯F_¯c( ¯y, t) Ru( ¯y, t)  , (2.41) where ¯Ru( ¯y, t) =: ¯Fu( ¯y, t) − U(t) ¯yu for some matrix U (t). Now we impose the condi-tions for coupling matrices G and D.

(i) λ = 0 is a simple eigenvalue of G(t), ∀ t ≥ 0 and e= √1

m(1, 1, . . . , 1) T

1×m is its corresponding eigenvector; (2.42a)

(ii) There is some λ > 0 such that µ2( ¯G(t)) ≤ −λ, ∀t ≥ 0. (2.42b) (iii) Coupling matrix D is of the form

D=  Ik 0 0 0  n×n. (2.42c) We are now in a position to state our first main theorem in the time-varying coupled system.

Theorem 2.3.3. Let coupling matrices G(t) and D satisfy (2.42). Suppose that ¯F, given in (2.29a) or (2.41), satisfies (2.14a), (2.14c), and that (2.14d), and system (2.29a) is bounded dissipative with respect to α. Then lim

t→∞y(t) = 0 for any initial¯ value provided that the coupling strength d satisfies the following inequality

d > b1 λ  1 + b 2 2 γ2 l 2 . (2.43)

Proof. For any initial condition ¯y(0), there is t0 > 0 such that k ¯y(t)k ≤ α for all t ≥ t0.

Applying the matrix measure inequality (2.1) and hypotheses (2.14a), (2.42b) on ¯yc, for any t ≥ t0, we have that

k ¯yc(t)k ≤ k ¯yc(t0)ke−λd(t−t0)+ b1α λd ≤ (e−λd(t−t0)₊ b1 λd)α =: (e−λd(t−t0)_{+ c} 0 1 d)α. Let δ > 1. We see that

k ¯yc(t)k ≤ α

dc0δ, (2.44a)

whenever t ≥ t0,1for some t0,1> t0. Similarly, applying inequality (2.1) and hypotheses (2.14c), (2.14d) on ¯yu1, k ¯yu1(t)k ≤ α d  b2 γc0  δ2 =: α dc1δ 2_{, (2.44b)}