Embedding dynamical networks into distributed models

Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Embedding dynamical networks into distributed models Giacomo Innocenti a, Paolo Paoletti b,⇑ a b

Dipartimento di Ingegneria dell’Informazione, Università degli Studi di Firenze, via S. Marta 3, I-50139 Firenze (FI), Italy Centre for Engineering Dynamics, School of Engineering, University of Liverpool, L69 3GH Liverpool, UK

a r t i c l e

i n f o

Article history: Received 9 June 2014 Received in revised form 25 November 2014 Accepted 17 December 2014 Available online 29 December 2014 Keywords: Nonlinear networks Embedding space Travelling waves

a b s t r a c t Large networks of interacting dynamical systems are well-known for the complex behaviours they are able to display, even when each node features a quite simple dynamics. Despite examples of such networks being widespread both in nature and in technological applications, the interplay between the local and the macroscopic behaviour, through the interconnection topology, is still not completely understood. Moreover, traditional analytical methods for dynamical response analysis fail because of the intrinsically large dimension of the phase space of the network which makes the general problem intractable. Therefore, in this paper we develop an approach aiming to condense all the information in a compact description based on partial differential equations. By focusing on propagative phenomena, rigorous conditions under which the original network dynamical properties can be successfully analysed within the proposed framework are derived as well. A network of Fitzhugh–Nagumo systems is finally used to illustrate the effectiveness of the proposed method. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction The study of networks composed by the interconnection of similar nonlinear dynamical systems is an intriguing field of investigation, motivated by its wide range of applications, ranging from the understanding of natural systems up to the control of artificial ones. For example, in the last decades cellular neural networks and analogue processor arrays have been introduced in engineering applications as a paradigm for analogue computation, image analysis and pattern generation [1,2]. Similarly, dynamical networks provide suitable models for the description of biological systems and their phenomena, such as, for example, the generation of spatio-temporal patterns in nervous systems of various organisms, see e.g. [3,4]. Recently, opinion dynamics [5] and infectious diseases spreading [6] have also been investigated by exploiting a dynamical network formulation of the problem. However, despite the interplay between the individual node dynamics and the interconnections turns out to be relevant in many scientific fields, its role in affecting the network behaviour is still not completely understood. In this regard, the main obstacle is represented by the intrinsically large dimension of the complete system state space, which scales linearly with the number of nodes and that, therefore, can easily reach the point where traditional analytical methods and numerical techniques fails or cannot be efficiently applied anymore. A dramatic reduction of the model complexity can be achieved by assuming that the dynamics of each node is actually a sampling of the dynamics of a continuous system governed by an underlying partial differential equation (PDE). In such a

⇑ Corresponding author. E-mail addresses: giacomo.innocenti@unifi.it (G. Innocenti), [email protected] (P. Paoletti). http://dx.doi.org/10.1016/j.cnsns.2014.12.009 1007-5704/Ó 2014 Elsevier B.V. All rights reserved.

22

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

framework each node is required to evolve close (in some norm) to the behaviour of the related PDE observed where that node is collocated, see e.g. [7] for a linear case of study and [8,9] for some preliminary results regarding chains of nonlinear systems. Notably, the continuous spatial coordinate system used to describe the nodes of the network does not necessarily represent a physical space, meaning that any interconnection structure can be represented within the proposed framework. It is also worth underlining that there exist phenomena, such as propagation failure, which may occur on networks but not on PDEs [10]. However, such degenerate situations seem to play a negligible role for all the applications mentioned before and this justifies further investigations on the relationships between networks and the corresponding PDE models. Remarkably, this problem has only recently received some attention in the literature and the results are still limited to special cases, such as the linear framework or simple topologies (see, e.g. [7–9,11–13]). For example, in [11] the authors derive a PDE from the original linear network using a finite difference method. Then, the study of the whole system stability is carried out by estimating the least stable eigenvalue. However, no guarantee of actual dynamic equivalence between the original network and the corresponding PDE is given. A different PDE approach is developed in [12], where the evolution of a continuum of mobile agents in the plane is described through an a priori fixed dynamics featuring a nonlinear time-varying Burgers-like PDE. Then, the class of the actual discrete network equations which fit in this representation is obtained by formal power series descriptions, and control strategies derived from Burgers PDE theory are used. In [13] some limitations of spatially invariant topologies are overcome by using distance-dependent coupling functions among the nodes, assuming that they act as spatially decaying operators. Such an approach, that takes into account only linear systems, leads to modelling the network as a continuum of agents, whose stability is then studied via functional operator equations. Also nonlinear networks have only recently attracted some attention, mainly thanks to the advances in consensus in nonlinear spaces. However, due to the complexity of the problem, there are only fragmented, although encouraging, results. In [14] a time-varying reference state, evolving according to a nonlinear model, is introduced. Then, the resulting problem is solved by assuming a globally defined reference frame and bidirectional information flows among all the agents. On the other hand, in [15] the nonlinear consensus problem is derived from a geometric interpretation of the linear case. In particular, the Euclidean mean is substituted with the so-called Induced Arithmetic Mean. Global consensus is then proved to strongly depend on the graph topology. In this regard, also in [16] the consensus is obtained under specific conditions on the graph. Moreover, here the agents are assumed to have low order dynamics subjected to global Lipschitz conditions on the nonlinearity (thus excluding, for example, all the class of polynomial functions). Similarly, in [17] the low order dynamics assumption, as well as the global Lipschitz condition, is used to ensure the consensus on strongly connected graphs. In this paper we propose an analytical method for the analysis of networks composed by individual nodes governed by similar (possibly, but not necessarily, the same) nonlinear ordinary differential equations (ODEs) and we discuss numerical implementation aspects. In Section 2 the network is embedded in a continuous spatial coordinate system in order to define a PDE whose dynamics interpolates the behaviour of each node. Conditions for the PDE to correctly represent the network dynamics are stated in Section 3. There we focus mainly on propagating phenomena where solutions of the PDE can be effectively expressed in a moving coordinate frame, thus obtaining an ODE whose solution provides the profile of the propagating wave. Several notions for the stability of this kind of systems are discussed in Section 4. Finally, an illustrative example is presented in Section 5 to show the effectiveness of the proposed approach and how it can be successfully used to analyse the behaviour of a non trivial nonlinear dynamical network. Some final remarks and future outlooks are discussed in Section 6. Some technical results are included in Appendix A for the sake of completeness. 1.1. Notation and preliminaries Throughout the paper we denote the sets of integer, real and complex numbers as, respectively, Z; R and C, whereas N0 represents the set of natural numbers including zero. We denote as Lp ðRm ; Rn Þ the Banach space of the measurable functions f ðxÞ : Rm ! Rn such that

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Z p p kf kLp :¼ kf ðxÞk dx < 1;

1 6 p < 1:

Rm

Given a function f ðxÞ 2 Lp ðRm ; Rn Þ, its derivative is defined, when it exists, as the linear operator

2 @f

1

@x1

6 df ðxÞ 6 :¼ 6 ... 4 dx

@f n @x1

...

@f 1 @xm

...

@f n @xm

3 7 7 7; 5

such that

lim

dx!0Rm

    df ðxÞ f ðx þ dxÞ  f ðxÞ  dx dx kdxk

¼ 0;

8x 2 R m :

The weak derivative of f ðxÞ, instead, is defined, when it exists, as the function @ x f ðxÞ 2 Lp ðRm ; Rn Þ such that

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

Z Rm

f ðxÞ

d/ðxÞ dx ¼  dx

Z Rm

23

@ x f ðxÞ/ðxÞdx;

for each test function /ðxÞ : Rm ! Rn with continuous derivatives of all orders and compact support in Rm . Higher order weak derivatives are addressed by the multi-index notation, that is, given the vectors T

x ¼ ½x1 ; . . . ; xm T 2 Rm ;

h ¼ ½ h 1 ; . . . ; h m  2 Nm 0;

we refer to the following quantities as follows

h! ¼ h1 ! . . . hm !;

xh ¼ xh11 . . . xhmm ;

@ hx ¼

@ h1 @xh11

...

@ hm @xhmm

¼

@ jhj @xh11

. . . @xhmm

:

We denote by Hk;p ðRm ; Rn Þ the Sobolev space defined as

n o Hk;p ðRm ; Rn Þ :¼ f ðxÞ 2 Lp ðRm ; Rn Þ : @ hx f ðxÞ 2 Lp ðRm ; Rn Þ; 8h 2 Nm 0 : jhj 6 k ; where the derivatives are intended in the weak sense. In particular, H0;p ðRm ; Rn Þ ¼ Lp ðRm ; Rn Þ. We also point out that Hk;p ðRm ; Rn Þ equipped with the norm

kf kHk;p :¼

 X  h  @ x f 

jhj6k

Lp

is a Banach space. It is also worth mentioning that Hk;2 ðRm ; Rn Þ is a Hilbert space equipped with the internal product

< f ; g >:¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi XZ @ hx f ðxÞ@ hx gðxÞdx; jhj6k

Rm

where gðxÞ is the complex conjugate of gðxÞ. We denote as Hk;p ðRm ; Rn Þ  Rq the Banach space given by the space Hk;p ðRm ; Rn Þ  Rq equipped with the norm

kðf ; fÞk ¼ kf kLp þ kfk;

f 2 Hk;p ðRm ; Rn Þ; f 2 Rq :

For a functional operator F : Lp ðRm ; Rn Þ ! Lq ðRr ; Rs Þ the induced ðLp ; Lq Þ-norm is defined as

kF kðLp ;Lq Þ :¼ sup

kf kLp <1

kFðf ÞkLq : kf kLp

We use the little-o notation oðf Þ to refer to any function g, whose Lp -norm grows at least polynomially faster than kf kLp ,   i.e. g is oðf Þ if limkxk!1 kgðxÞkLp =kf ðxÞkLp ¼ 0. A (functional) operator H, which is continuously (Frétchet) differentiable on a certain set U, is said to belong to the class C 1 ðUÞ and its derivative operator is denoted by DH. We also denote as vecðAÞ : Rv w ! Rv w the linear operator such that

vecðAÞ ¼ ½a11 ; . . . ; av 1 ; a12 ; . . . ; av 2 ; . . . ; a1w ; . . . ; av w T ;   being A ¼ aij for i ¼ 1; . . . ; v and j ¼ 1; . . . ; w. Eventually, for the sake of conciseness we will alternatively use f  g to denote the composite function f ðgðxÞÞ : Rn ! Rm , given f : Rh ! Rm and g : Rn ! Rh . 2. Problem set up In this paper we consider networks composed by different interconnected agents governed by the dynamical equations

X d ui ðtÞ ¼ F i  ui ðtÞ þ C i;j  uj ðtÞ; dt j2M

i 2 N;

ð1Þ

i

where N # Zw is a set of multi-indexes labelling the nodes, ui : R ! Rn is the ith node’s state, F i : Rn ! Rn are static nonlinearities describing the agents’ autonomous dynamics, C i;j : Rn ! Rn is the function linking the ith node to the jth one, and Mi # N is the set of the neighbours of the ith node. Our goal is to embed the network (1) into a proper space Rm ; m being its dimension, acting as a system of (potentially virtual) spatial coordinates. To this aim, let us denote as x 2 Rm the set of coordinates representing a point in the embedding space, and more specifically let xi denote the position associated to the ith node. We point out that in some applications xi might actually represent the physical location of the ith agent, but in general such assumption is not required, and so xi simply represents a label to identify the node in the abstract embedding space Rm . Such a distinction among actual and abstract coordinates is in fact the crucial feature for the developing of the sought embedding theorem.

24

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

The first step consists in introducing a displacement map E # Zw such that 8i; j 2 N there exists satisfied. Such a map always exists and it allows one to rewrite (1) as

X d u ðtÞ ¼ F i  ui ðtÞ þ Bi;ji  uj ðtÞ; dt i j2Eþi

j 2 E so that j ¼ i þ j is

i 2 N;

ð2Þ

where the operators Bi;ji are defined as

 Bi;ji ¼

C i;j ; j 2 Mi ; 0;

otherwise;

and E þ i denotes the set

E þ i ¼ fj 2 Zw : j ¼ i þ j;

j 2 E g:

Then, let us also introduce the arbitrary functions, / : R  Rm ! Rn ; F : Rn  Rm ! Rn and Bj : Rn  Rm ! Rn , which spatially interpolate respectively the states of the nodes, the functions F i and the links Bi;ji , i.e.

/ðt; xi Þ ¼ ui ðtÞ;

/ðt; xÞ 2 Hlþ1;p ðR  Rm ; Rn Þ;

FðuðtÞ; xi Þ ¼ F i  uðtÞ;

ð3Þ

Fðu; xÞ 2 Hlþ1;p ðRn  Rm ; Rn Þ;

Bji ðuðtÞ; xi Þ ¼ Bi;ji  uðtÞ;

ð4Þ

Bji ðu; xÞ 2 Hlþ1;p ðRn  Rm ; Rn Þ;

ð5Þ

where l P 0. Since / is thought to interpolate the solutions at values of x different from the nodes locations, one can conveniently assume that in these regions it shows a regular behaviour. For instance, the study could be reduced to functions /, which are continuous with respect to x, or which are bounded on compact sets. Similar assumptions can be enforced onto F i and Bi;ji as well. Then, denoting for the sake of compactness the operator Fð ; xÞ as FðxÞ and the operators Bji ð ; xÞ as Bji ðxÞ, and rearranging the indexes, Eq. (2) becomes

X d /ðt; xi Þ ¼ Fðxi Þ  /ðt; xi Þ þ Bj ðxi Þ  /ðt; xiþj Þ; dt j2E

i 2 N:

Remembering that / is regular but arbitrary, just assume that each coupling operator Bj ðxi Þ; j 2 E, is sufficiently smooth that Bj ðxi Þ  /ðt; xÞ can be developed for every x within a ball of radius r according to the Taylor expansion formula:

Bj ðxi Þ  /ðt; x þ dxÞ ¼

Z 1 X lþ1     @ hx Bj ðxi Þ  /ðt; xÞ dxh þ dxh ð1  aÞl @ hx Bj ðxi Þ  /ðt; x þ adxÞ da; 8kdxk < r: h! h! 0 jhj6l jhj¼lþ1

X1

ð6Þ

Then, assuming that the radius r is sufficiently large with respect to the embedding coordinates, let us denote the spatial displacement between two neighbouring nodes as

xiþj ¼ xi þ eij mij ; where eji 2 Rþ : eji < r represents  the distance between the nodes in the embedding space, and tion of the displacement, i.e. mij  ¼ 1. Therefore, the dynamics of the ith node is governed by

XX 1 jhj   X d /ðt; xi Þ ¼ Fðxi Þ  /ðt; xi Þ þ eij mhij @ hx Bj ðxi Þ  /ðt; xi Þ þ Rlþ1 ð/ðt; xi Þ; eij ; mij Þ; dt h! j2E jhj6l j2E

mij 2 Rm identifies the direc-

ð7Þ

where we have defined the residual function

Rlþ1 ð/ðt; xi Þ; eij ; mij Þ :¼

X l þ 1 jhj Z 1   eij mhij ð1  aÞl @ hx Bj ðxi Þ  /ðt; xi þ aeij mij Þ da: h! 0 jhj¼lþ1

It is worth observing that the displacement map E undergoes substantial simplifications when the network topology features some sort of regularity. In such a case, indeed, the connections show a repeated pattern that reduces the number of operators Bj ðxi Þ which are not null, and therefore the necessary displacement indexes boil down to those needed to represent the pattern itself. For instance, if the embedding space describes a nonuniform lattice with tetrahedral pattern, only four connecting functions per node are necessary to describe the whole link topology, i.e. E counts only four indexes. A special case emerges when the connecting operators Bj ðxi Þ do not depend on xi , i.e.

Bj ðxi Þ ¼ Bj ;

8i 2 N :

ð8Þ

In such a condition every system has a structurally similar neighbourhood and it is affected by those neighbours in the very same way. This regularity ensures that the behaviour of each agent is dynamically equivalent to that of all the others, no matters where it is located in the embedding space. The resulting network is then a ‘‘translation (or shift) invariant set’’ with

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

25

respect to the embedding coordinate x. This class of systems is particularly suited for the application of the following PDE approach, and it covers important models such as lattices, cellular automata, chemical reaction networks and transportation systems, which are widely used in several fields counting Biology, Chemistry, Engineering and Quantum Physics (see, e.g. [18–22] and references therein). Observe that the above dynamics (7) is defined point-wise in Rm , but /; F and the Bj ’s are interpolating functions defined everywhere in the embedding domain. Therefore, we may switch to the study of the PDE obtained from (7) after removing the dependency from the index i:

XX 1 jhj   X d /ðt; xÞ ¼ FðxÞ  /ðt; xÞ þ ej mhj @ hx Bj ðxÞ  /ðt; xÞ þ Rlþ1 ð/ðt; xÞ; ej ; mj Þ: dt h! j2E jhj6l j2E

ð9Þ

We anticipate here that for the development of the final result in the following section, each operator in (9) will have first to be proven continuous and continuously differentiable with respect to /. In this respect, it is worth underlining that the partial derivatives of the composite functions Bj ðxÞ  /ðt; xÞ are polynomials in Bj ; /ðt; xÞ and their partial derivatives up to the order jhj, whose coefficients can be computed, for instance, according to the Faà di Bruno formula [23]. Therefore, by the regularity property of the derivative operator (see Lemma Appendix A.6), it is straightforward to verify that, if Bj is con  tinuous and continuously differentiable up to the order jhj þ 1 with respect to x, then @ hx Bj ðxÞ  /ðt; xÞ is a (potentially nonlinear) continuous and continuously differentiable operator with respect to /. The potentially nonlinear nature of such operator, however, is not crucial for the following development, since it only implies the need of collecting the nonlinear part with FðxÞ in order to isolate a single overall nonlinear operator. Therefore, for the sake of simplicity, in the following we limit the analysis to the case of shift invariant F and linear shift invariant connecting functions Bj , i.e

FðxÞ ¼ F; 8x 2 Rm ; Bj ðxÞ  /ðt; xÞ :¼ Aj /ðt; xÞ; 8x 2 Rm ; Observe that the role of the unique nonlinear operator is left to F alone and that the translation invariance property (8) is satisfied as well. We stress once more that the general case (9) can be approached exactly as it will be presented in the following section for (10) under the assumption that the connecting functions Bj ; j 2 E, are sufficiently regular, i.e. that Bj belongs to the class C lþ2 ðRn Þ. Also the absence of the shift invariant property in the general case only implies more complicated formulas, but it does not invalidate the proposed approach. Hence, Eq. (9) conveniently reduces to

XX 1 jhj X d /ðt; xÞ ¼ F ð/ðt; xÞÞ þ ej mhj Aj @ hx /ðt; xi Þ þ Rlþ1 ð/ðt; xÞ; ej ; mj Þ; dt h! j2E jhj6l j2E where

Rlþ1 ð/ðt; xÞ; ej ; mj Þ ¼

ð10Þ

Z 1 X l þ 1 jhj ej mhj Aj ð1  aÞl @ hx /ðt; x þ aej mj Þda: h! 0 jhj¼lþ1

Again, for the sake of simplicity and without loss of generality we will restrict the analysis to l ¼ 1, since in such a case the formulas admit compact explicit forms, being the extension to l > 1 straightforward, though complicated by the necessity of using higher order chain rule formulas, such as the Faà di Bruno formula [23]. In order to study the propagative and synchronisation phenomena that are the main focus of this paper, we introduce the moving coordinate

n :¼ kx  ct 2 Rm ;

k 2 Rmm ; c 2 Rm ;

where detðkÞ – 0 so that the coordinate change ðt; xÞ # ðt; nÞ is invertible. Moreover, we also introduce a new function w such that

/ðt; xÞ ¼: wðt; kx  ctÞ ¼ wðt; nÞ 2 Rn : Observe that the chain rule for derivatives implies

X h d @ wðt; nÞ ¼ wðt; nÞ  @ n wðt; nÞch ; dt @t jhj¼1 d d wðt; nÞ ¼ wðt; nÞk; dx dn X h h @ ix wðt; nÞ ¼ @ n wðt; nÞki ;

8i 2 Nm 0 : jij ¼ 1;

jhj¼1

where ki denotes the column of k corresponding to the only non zero element of i. Hence, if l ¼ 1, Eq. (10) can be rewritten in this new coordinate system as

@ t wðt; nÞ 

X jhj¼1

ch @ hn wðt; nÞ ¼ Fðwðt; nÞÞ þ

X XX X X   Aj wðt; nÞ þ ej mij khi Aj @ hn wðt; nÞ þ R2 wðt; kx  ctÞ; ej ; mj ; j2E

j2E jij¼1jhj¼1

j2E

26

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

and we eventually obtain

@ t wðt; nÞ ¼ Fðwðt; nÞÞ þ G0 wðt; nÞ þ G1 ð@ n wðt; nÞ; e; m; k; cÞ þ C2 ðwðt; nÞ; e; m; k; cÞ;

ð11Þ

where for convenience we have introduced the following operators and sets of coefficients:

 

e ¼ ej

j2E

;

X Aj ; G0 ¼

 

m ¼ mj

j2E

;

G1 ð@ n w; e; mÞ ¼

j2E

C2 ðwðt; nÞ; e; m; k; cÞ ¼

X XX

! i h j j k i Aj

em

jhj¼1

þ c In @ hn wðt; nÞ; h

j2E jij¼1

Z 1 X   X X 2 jhj h   R2 wðt; kx  ctÞ; ej ; mj ¼ ej mj Aj ð1  aÞ2 @ hx w t; kx  ct þ aej mj da; h! 0 j2E j2E jhj¼2

being In 2 Rnn the identity matrix. Eq. (11) will be the basis for the main results about the equivalence between the original network behaviour and the dynamics of the PDE solutions in the embedded domain. Furthermore, special relevance will be given to the stability properties of both such solutions.

3. Main results The main focus of this paper are patterns and propagative phenomena on networks, so they can be represented as stationary solutions w0 2 H2;2 ðR  Rm ; Rn Þ of (11), i.e.

@ t w0 ðt; nÞ ¼ Fðw0 ðt; nÞÞ þ G0 w0 ðt; nÞ þ G1 ð@ n w0 ðt; nÞ; e; m; k; cÞ þ C2 ðw0 ðt; nÞ; e; m; k; cÞ ¼ 0H0;p : It is worth stressing that problem (11) depends on the topology of the network, which is summarised by e and m. Therefore, also its stationary solutions depend on the embedding space, i.e. w0 ðt; nÞ ¼ w0 ðt; n; e; mÞ. We point out that here we are restricting the solution space to H2;2 (instead of a generic H2;p ) as this space has the structure of a Hilbert space and therefore in Section 4 we will be able to decompose any neighbourhood of a given solution in orthogonal subspaces corresponding to unstable and stable directions. However, the results presented in this section are valid for any Hlþ1;p with l P 1. Since w0 is stationary, it does not depend explicitly on t, i.e. w0 ðt; nÞ ¼ w0 ðnÞ. Moreover, for each . 2 Rm let us define the ntranslation operator sð ; .Þ as

sðwðt; nÞ; .Þ ¼ wðt; n þ .Þ: Observe that sðw0 ; .Þ is a solution of (11) for each

ð12Þ m

. 2 R as well. Indeed,

Fðsðw0 ðnÞ; .ÞÞ þ G0 sðw0 ðnÞ; .Þ þ G1 ð@ n sðw0 ðnÞ; .Þ; e; m; k; cÞ þ C2 ðsðw0 ðnÞ; .Þ; e; m; k; cÞ ¼ Fðw0 ðn þ .ÞÞ þ G0 w0 ðn þ .Þ þ G1 ð@ n w0 ðn þ .Þ; e; m; k; cÞ þ C2 ðw0 ðn þ .Þ; e; m; k; cÞ ¼ 0H0;p ¼ @ t w0 ðn þ .Þ; and so, by changing variable into n0 ¼ n þ ., the equivalence is proven. Notice also that such a result derives from the spatial invariance (8) initially assumed for the original network. Hence, we define the set of the solutions related to w0 as





W0 ¼ w 2 H2;2 ðR  Rm ; Rn Þ : wðt; nÞ ¼ sðw0 ðnÞ; .Þ; . 2 Rm :

ð13Þ

Observe that also the above set depends on the choice of the embedding space, i.e. W0 ¼ W0 ðe; mÞ. By Proposition Appendix A.8 it follows that W0 is a manifold, since it is the image of Rm through the .-operator sðw0 ; .Þ. Theorem 3.1. Let us derive from the original problem (11) the reduced problem

@ t wðt; nÞ ¼ Fðwðt; nÞÞ þ G0 wðt; nÞ þ G1 ð@ n wðt; nÞ; e; m; k; cÞ;

ð14Þ

consisting in a finite order PDE. Let F 2 C 1 ðRm Þ almost everywhere, JF 2 Rm  Rm being its Jacobian matrix, and define the linear operator T 1 : H2;2 ðR  Rm ; Rn Þ ! L2 ðR  Rm ; Rn Þ as

T 1 ðwÞdw :¼ @ t dw  JFðwÞdw  G0 dw  G1 ð@ n dw; e; m; k; cÞ:

ð15Þ

^ 2 W0 of (11) there exists a solution w ~ of (14) that is diffeomorphic to w, ^ if there exists a sufficiently small Then, for each solution w ^ such that neighbour J of w



sup kT 1 ðwÞkðH2;2 ;L2 Þ þ kC2 kðH2;2 ;L2 Þ < r; w2J

for some

r 2 ð0; 1Þ. In particular, w~ 2 J .

ð16Þ

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

27

Proof. Since this proof is quite long and technical, let us first describe just its outline before entering into the details. First, we characterise the functions w, which are solutions of the original problem (11) and of the reduced problem (14), in terms of the kernel of a special functional operator Q interpolating Eqs. (11) and (14). For each solution of the original problem (11) we will then build a new operator H, that will be sufficiently regular to apply the Implicit Function Theorem. By investigating the invertibility range of the derivative of such operator, we will finally find conditions on the interpolation of the original network nodes for the operator itself to be invertible on a domain that is sufficiently large to cover the solutions of the reduced problem as well. Therefore, by the Implicit Function Theorem the original and reduced problem solutions will be related by a diffeomorphism, concluding the proof. Let us then introduce the operator Q ðw; fÞ : H2;2 ðR  Rm ; Rn Þ  R ! L2 ðR  Rm ; Rn Þ defined as

Q ðw; fÞ :¼ @ t w  FðwÞ  G0 w  G1 ð@ n w; e; m; k; cÞ  fC2 ðw; e; m; k; cÞ: By Lemmas Appendix A.4, Appendix A.6 and Appendix A.5, it follows that @ t ; G0 ; G1 and C2 are linear, bounded and continuously differentiable operators with respect to w. Moreover, F is continuous and continuously differentiable almost everywhere in Rm by assumption and thus, by applying the chain rule to each w, we also have

Fðw þ dwÞ ¼ FðwÞ þ JFðwÞdw þ oðdwÞ: Observe that Q is also continuous with respect to f, since C2 is bounded. Moreover, by the definition it follows that the (Frétchet) derivative of Q with respect to f is just C2 . Therefore, Q is continuous and continuously differentiable in each ðw; fÞ and, by the chain rule, its derivative is given by the linear operator

DQðw; fÞðdw; dfÞ ¼ @ t dw  JFðwÞdw  G0 dw  G1 ð@ n dw; e; m; k; cÞ  fC2 ðdw; e; m; k; cÞ  dfC2 ðw; e; m; k; cÞ: In particular, we have

Dw Q ðw; fÞdw ¼ @ t dw  JFðwÞdw  G0 dw  G1 ð@ n dw; e; m; k; cÞ  fC2 ðdw; e; m; k; cÞ: Let us now define for each

. 2 Rm the parametric operator H. : H2;2 ðR  Rm ; Rn Þ  R ! L2 ðR  Rm ; Rn Þ as

H. ðw; fÞ :¼ Q ðw; fÞ þ ðw  sðw0 ; .ÞÞ:

ð17Þ

Because of the properties of Q ; H. turns out to be continuous and continuously differentiable in each ðw; fÞ and its derivative is given by

DH. ðw; fÞðdw; dfÞ ¼ DQðw; fÞðdw; dfÞ þ dw: In particular, we have

Dw H. ðw; fÞdw ¼ Dw Q ðw; fÞdw þ dw;

ð18Þ

that does not depends on . any more. A number of observations are in order. First, notice that Q ðw; 1Þ ¼ 0L2 for each w 2 W0 , i.e. Q ðsðw0 ; .Þ; 1Þ ¼ 0L2 , 8. 2 Rm , and that, similarly, H. ðsðw0 ; .Þ; 1Þ ¼ 0L2 . Then, observe that for all . þ d. 2 Rm

 d 0L2 ¼ Q ðsðw0 ; . þ d.Þ; 1Þ ¼ Q sðw0 ; .Þ þ sðw0 ; .Þd. þ oðd.Þ; 1 dn ¼ Q ðsðw0 ; .Þ; 1Þ þ Dw Q ðsðw0 ; .Þ; 1Þ

d d sðw0 ; .Þd. þ oðd.Þ ¼ Dw Q ðsðw0 ; .Þ; 1Þ sðw0 ; .Þd. þ oðd.Þ: dn dn

It follows that each dw satisfying dw ¼ ðd=dnÞsðw0 ; .Þd. for some d. 2 Rm belongs to ker Dw Q ðsðw0 ; .Þ; 1Þ, i.e.

Dw Q ðsðw0 ; .Þ; 1Þdw ¼ Dw Q ðsðw0 ; .Þ; 1Þ

d sðw0 ; .Þd. ¼ 0L2 ; dn

8d. 2 Rm :

Hence, Dw Q ðsðw0 ; .Þ; 1Þ is not injective and so not even invertible. However, observe that

Dw Q ðw; fÞdw ¼ T 1 ðwÞdw  fC2 ðdw; e; m; k; cÞ 2 H2;2 ðR  Rm ; Rn Þ;

8ðw; fÞ;

ð19Þ

^ 2 W0 the operator with T 1 as defined in (15), and so by (16) it follows that in the neighbourhood J of the nominal solution w (19) satisfies



  sup Dw Q ðw; fÞðH2;2 R;L2 Þ < supðjfjÞ kT 1 ðwÞkðH2;2 ;L2 Þ þ kC2 kðH2;2 ;L2 Þ < 1;

ðw;fÞ2J I

ð20Þ

f2I

  where I is the connected open set r1 ; r1 containing both f ¼ 1 and f ¼ 0. Hence, the application of Corollary Appendix A.2 implies that the operator Dw H. ðw; fÞ is invertible in any ðw; fÞ 2 J  I independently from .. ^ of To conclude the proof, we then only need to show that the invertibility domain of such operator includes the solution w ^ and define P. ðw; fÞ : J  I ! J  I as the reduced problem. To this end, let . 2 Rm such that sðw0 ; .Þ ¼ w,

28

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39





P. ðw; fÞ ¼ D1 w H. ðsðw0 ; .Þ; 1ÞH . ðw; fÞ; f : Denoting as Nu and Id respectively the null and the identity operators, since

" DP. ðw; fÞ ¼

D1 D1 w H . ðsðw0 ; .Þ; 1ÞDw H. ðw; fÞ w H . ðsðw0 ; .Þ; 1ÞDf H . ðw; fÞ Nu

Id

# ;

we have that DP. ðw; fÞ is invertible for all the ðw; fÞ such that Dw H. ðw; fÞ is invertible as well. Therefore, DP. ðw; fÞ is invertible for each ðw; fÞ 2 J  I . Then, since it results that

"

DP1 . ðw; fÞ ¼

D1 Dw H1 . ðw; fÞDf H. ðw; fÞ w H . ðw; fÞDw H . ðsðw0 ; .Þ; 1Þ Nu

Id

#

;

we also have " DP1 . ðsðw0 ; .Þ;1ÞDP. ðw; fÞ  Id ¼

  1  # D1 Dw H. ðsðw0 ; .Þ;1Þ Df H. ðw; fÞ  Df H. ðsðw0 ; .Þ; 1Þ w H . ðsðw0 ; .Þ;1Þ Dw H. ðw;fÞ  Dw H. ðsðw0 ; .Þ; 1Þ : Nu Nu

Because DHðw; fÞ is continuous in J  I , we can always build a small open neighbourhood of ðsðw0 ; .Þ; 1Þ, namely V  W  J  I , such that

    sup DP1 < 1: . ðsðw0 ; .Þ; 1ÞDP. ðw; fÞ  Id 2;2 ðH R;L2 Þ ðw;fÞ2VW

ð21Þ

^ ¼ sðw ; .Þ, and observing that Recalling that w 0

      ^ 1Þ ^  C2  limDw H. ðw; 0Þ  Dw H. ðw;  2;2 2 ¼ limT 1 ðwÞ  T 1 ðwÞ  2;2 2 ¼ kC2 kðH2;2 R;L2 Þ < r; ^ ðH R;L Þ w!w^ ðH R;L Þ w!w         ^ ^ w;  C ðw; e ; m Þ þ C ð e ; m Þ ¼ lim limDf H. ðw; 0Þ  Df H. ðw; 1Þ 2;2   2;2 2 ¼ 0; 2 2 ^ ðH R;L2 Þ w!w^ ðH R;L Þ w!w

one can always choose V sufficiently small to satisfy (21) for W containing zero, i.e. the invertibility of Dw H. ðw; fÞ can be granted on a V  W such that 0 2 W. Therefore, by Theorem Appendix A.3 it follows that also H. is invertible in that V  W. ^ 2 V; H. ðw; ^ 1Þ ¼ 0 p and Dw H. ðw; ^ 1Þ is invertible, the hypothesis of the Implicit Function Since H. is C 1 ðV  WÞ; w L Theorem (on functional Banach spaces, see also [24, chapt. 22]) are satisfied. Therefore, there exists an open neighbourhood ^ ðu. ðfÞ; fÞ 2 V  U and H. ðu. ðfÞ; fÞ ¼ 0 p ; 8f 2 U. U # W of f ¼ 1 and a unique C 1 ðWÞ function u. : W ! V such that u. ð1Þ ¼ w; L In order to investigate the extent of U, let p1 ðw; fÞ ¼ w 8ðw; fÞ 2 V  U and observe that u. ðfÞ p1  P1 . ð0Lp ; fÞ in W. Indeed, ^ 1Þ ¼ ð0 p ; 1Þ. Moreover, P. ðw; fÞ is invertible in the set of the couples ð0Lp ; fÞ belonging to V  W is not empty, since Pðw; L V  W and in particular in each ð0Lp ; fÞ 2 V  W, i.e. for each of these ð0Lp ; fÞ there exists a unique couple 1 ðf . ðfÞ; g . ðfÞÞ 2 V  W such that P. ð0Lp ; fÞ ¼ ðf . ðfÞ; g . ðfÞÞ. Thus, we have

ð0Lp ; fÞ ¼ P. ðf . ðfÞ; g . ðfÞÞ ¼ ðH. ðf . ðfÞ; g . ðfÞÞ; g . ðfÞÞ;

ð22Þ

which is equivalent to

g . ðfÞ ¼ f; H. ðf . ðfÞ; fÞ ¼ 0Lp : ^ Hence, u. f . ¼ p1  P1 . ð0Lp ; fÞ 2 V 8f 2 W, i.e. U W. Because of the latter equivalence ðw; 1Þ, representing a solution of the original problem (11), is diffehomorphic to any ðu. ðfÞ; fÞ and in particular to ðu. ð0Þ; 0Þ, which corresponds to the solution ~ :¼ u. ð0Þ of the reduced problem (14). h w ^ 2 W0 , for possibly different values of r. Then, Proposition 3.2. Assume that the conditions of Theorem 3.1 are satisfied for each w the reduced order problem (14) admits a set of solutions that is diffeomorphic to W0 . Proof. Observe that by Proposition Appendix A.8 the set W0 turns out to be a manifold parametrized by .. Moreover, by the same proposition it follows that the operator s in (12) is continuous and continuously differentiable and, thus, that H. as defined in (17) is continuous and continuously differentiable with respect to . as well. Hence, the same holds for the solutions of problem (22) and in particular for ðu. ð0Þ; 0Þ, that represents the solution of the reduced problem at .. Therefore, the set of the solutions of problem (14) is continuously parametrized by . and each of them is diffeomorphic to a distinct solution of W0 related to that .. h Remark 1. Condition (16) can be relaxed by considering the induced norms applied to the restrictions of the operators T 1 ðwÞ and C2 to the only set J .

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

29

^ w ~ be soluLet us define the operator Q 1 ðwÞ :¼ Q ðw; 1Þ þ C2 ðw; e; m; k; cÞ, and observe that DQ 1 ðwÞ ¼ T 1 ðwÞ. Moreover, let w; ^ 1Þ ¼ 0Rn ; Q ðw; ~ 0Þ ¼ 0Rn , then tions of problems (11) and (14), respectively, as in Theorem 3.1. Since, Q ðw; ^ ¼ C2 ðw; ^ e; m; k; cÞ; Q 1 ðwÞ ~ ¼ 0. Therefore, thanks to continuity and continuous differentiability of Q 1 , it holds that Q 1 ðwÞ ^ e; m; k; cÞ ¼ Q 1 ðwÞ ^  Q 1 ðwÞ ~ T 1 ðwÞð ^ w ^  wÞ. ~ C2 ðw; ^ of the original problem measures the distance between w ^ and its diffeomorphic Remark 2. The value of C2 at the solution w     ~ solution of the reduced problem. The smaller C2 ðw; ^ e; m; k; cÞ ^ and w ~ are. The same effect could be w, is, the closer w 2;2 H

obtained by bigger values of the induced norm of T 1 ðwÞ (possibly restricted to J , as observed in Remark 1), but condition (16)    ^ e; m; k; cÞ prevents such a possibility in order to make operator Q a contraction as pointed out in (20). Since C2 ðw;  H2;2

^ through its Taylor expansion, it follows that using w ~ instead of w ^ will provide better intrinsically depends on the nature of w results for those solution of the original problem which are better approximated by their truncated Taylor series.

Remark 3. Theorem 3.1 also applies to higher order expansions of the connecting links among the systems, i.e. for l > 1. In such a case, G1 and C2 are substituted by two higher order operators, namely Gl and Clþ1 , which share with the first ones similar properties of continuity and continuous differentiability, but whose explicit forms are too complicated to be studied    ^ e; mÞ remains unchanged from Remark 2. Therefore, the higher in the general case. However, the meaning of Clþ1 ðw;  H2;2

order version of problem (14) is suitable to describe waves behaviours which are more complex with respect to the interpolating network. Remark 4. For the sake of the simplicity in this paper only Banach spaces of the form used in (3) have been considered. However, Theorem 3.1 holds true for any functional Banach space provided that the network satisfies equivalent assumptions. For instance, the extension to Banach space of periodic functions over a compact interval is immediate, and it will be illustrated in the example of Section 5. Remark 5. A complete analysis of the propagating phenomena in the network requires to study problems (11) and (14) with respect to different values of the parameters k and c characterising the nature of the sought wave.

4. On the stability properties of travelling wave solutions Theorem 3.1 shows that, under the conditions commented in Remarks 1–4, stationary solutions of (14) correspond to as much rest states, travelling waves, pulses and fronts in the original dynamical network governed by (1). The extent of that result is not limited to the existence condition alone. Indeed, the theorem states that the solutions in the embedding space are diffeomorphic to the corresponding ones in the original space, i.e. that the dynamical properties are preserved as well. Therefore, the stability analysis can be equivalently carried out in the embedding space without any loss. Corollary 4.1. Under the conditions of Theorem 3.1 the stability properties of W0 and those of the corresponding solutions in the original discrete network (1) are the same.

Proof. The result directly follows by the existence of a diffeomorphism between them. h Then, the stability of the class of solutions of (1) represented by W0 can be investigated directly via the PDE representation (14) by recurring to any fitting method or technique. However, it is worth stressing that with respect to a set W0 of non solitary solutions such as (13) there is no unique definition of stability. Therefore, following the most common approaches in the literature (see, e.g. [25–27]), we introduce two different definitions of stability, also providing some useful insights about the conditions which have to be met in order to enforce the corresponding behaviour. Definition 4.1. Let w0 ðt; nÞ ¼ w0 ðnÞ be a stationary solution of (14). Then, w0 ðnÞ is (strictly) spectrally (or linearly) stable if for each sufficiently small perturbation Dwðt; nÞ the corresponding trajectory wðt; nÞ ¼ w0 ðnÞ þ Dwðt; nÞ converges exponentially to w0 ðnÞ, i.e. if there exists k 2 C, with RðkÞ < 0, and lðnÞ 2 Rn such that

  limDwðt; nÞ  lðnÞekt  ¼ 0; 8n 2 Rm :

t!1

ð23Þ

Observe that considering condition (23) is equivalent to investigating the linear behaviour of the system dynamics around w0 ðnÞ. Therefore, let us consider the linearisation of the operators in (14). With respect to the nonlinear term F we have

Fðw0 ðnÞ þ Dwðt; nÞÞ ¼ Fðw0 ðnÞÞ þ @ w0 Fðw0 ðnÞÞDwðt; nÞ þ    Moreover, notice that G1 is a linear operator with respect to @ n and therefore it can be more conveniently rewritten as

30

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

^ 1 ðe; m; k; cÞ vecð@Þ; G1 ð@; e; m; k; cÞ ¼ G ^ 1 ðe; m; k; cÞ. Therefore, the linearised dynamics around w ðnÞ reads for some suitable matrix G 0

^ 1 ðe; m; k; cÞ vecð@ n DwÞ þ G0 Dw: @ t Dw ¼ @ w0 Fðw0 ÞDw þ G

ð24Þ

This notation will be useful for the numerical implementation of this theory, as shown in the next section. Then, by substituting Dwðt; nÞ ¼ lðnÞekt into (24) one obtains

^ 1 ðe; m; k; cÞ vecð@ n DwÞ ¼ kDw; Ln ðDwÞ ¼ @ w0 Fðw0 ÞDw þ G0 Dw þ G

ð25Þ

where the linear parametric operator Ln as been introduced for the sake of conciseness. Notice that (25) is a spectral problem over Ln . Definition 4.2. The spectrum of the linear operator Ln : H2;p ðR  Rm ; Rn Þ ! H2;p ðR  Rm ; Rn Þ, denoted as spðLn Þ, is the set of all the values k 2 C, referred to as eigenvalues, such that the operator Ln  kId is not invertible (Id being the identity operator). If kerðLn  kIdÞ does not contain only the null function Nu 2 H2;p ðR  Rm ; Rn Þ, then its other functions are denoted as the eigenfunctions associated to k. The set of all the linear combinations of the eigenfunctions of k is called eigenspace and it is denoted as spanðkÞ. Remark 6. When p ¼ 2 the Sobolev space H2;2 ðR  Rm ; Rn Þ is also a Hilbert space and so the spectral problem provides a complete characterisation of the local behaviour around the nominal solution. Remark 7. The spectrum of the parametric linear operator Ln defined on a Banach space of periodic functions over a compact interval is such that spðLn Þ is constant with respect to n, while for each k 2 spðLn Þ the associated eigenspace spanðkÞ changes accordingly and periodically. In the literature the spectrum is usually parted into several sets according to the reason Ln  kId cannot be inverted for (see [28,27] and the references within for details), but for the purpose of this paper we will only refer to the spectrum as a whole. The spectral stability of the stationary solutions set W0 can then be formulated as follows. Definition 4.3. The set W0 defined as in (13) is spectrally stable if the spectrum spðLn Þ of the operator Ln , defined as in (25), for each n comprises only values kðnÞ with negative real part but a single kðnÞ ¼ 0. The above definition is motivated by the observation that for a set such as W0 in every arbitrary small neighbourhood of any w0 ðnÞ 2 W0 there are infinite other stationary solutions belonging to W0 itself. Therefore, by choosing the initial perturbation value as lðnÞ ¼ w0 ðn þ .Þ  w0 ðnÞ for an arbitrary small ., the corresponding solution w0 ðnÞ þ Dwðt; nÞ ¼ w0 ðn þ .Þ turns out to be stationary, i.e. it is naturally associated to k ¼ 0. The analytic computation of spðLn Þ as n varies is in general a formidable task. Nevertheless, in the literature some numerical approaches have been presented, see e.g. [29–31]. One of the most efficient methods recurs to numerical continuations tool in order to build the spectra from the knowledge of spðLn Þ for a certain n (see [30,32] and references therein for detailed explanations). It is worth observing that spectral stability requires the perturbed trajectory to be absorbed exactly on the initial nominal solution. In many situations such a condition turns out to be too restrictive and, then, one allows the perturbed solution starting from w0 ðnÞ to land on a certain w0 ðn þ .Þ or to continuously move along and close to W0 . In order to take into account for the phase drift, i.e. the chance of the perturbed solution to get close to W0 regardless the original starting nominal solution, a different definition of stability is needed (see [27] for a detailed review of the topic). Definition 4.4. The set W0 defined as in (13) is orbitaly (or nonlinearly) stable if for each w0 ðnÞ 2 W0 and for any small initial perturbation Dwð0; nÞ the corresponding perturbed trajectory wðt; nÞ ¼ w0 ðnÞ þ Dwðt; nÞ converges to the set W0 itself, i.e.

lim inf kwðt; nÞ  w0 ðn þ .Þk ¼ 0; 8n 2 Rm :

t!1.2Rm

ð26Þ

Remark 8. Spectral stability directly implies orbital stability, but not viceversa. A particular situation where orbital stability occurs is when the phase drift tends to a constant shift for each perturbed solution. In such a case, starting from close to w0 ðnÞ 2 W0 the trajectory moves diverging to an other stationary solution w0 ðn þ .Þ 2 W0 . This is equivalent to the existence of a front interconnecting them, i.e. the presence of an heteroclinic solution. Definition 4.5. Let mðt; nÞ be a non stationary solution of (14). Then, mðt; nÞ is said a (exponential) front starting from w0 ðnÞ and landing to w0 ðn þ .Þ, if it exists lu ðnÞ; ls ðn þ .Þ 2 Rn and ku ; ks 2 C along with the additional condition

Rðks Þ < 0 < Rðku Þ; such that

ð27Þ

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

  lim mðt; nÞ  w0 ðnÞ  lu ðnÞeku t  ¼ 0;

t!1

  limmðt; nÞ  w0 ðn þ .Þ  ls ðn þ .Þeks t  ¼ 0;

t!1

8n 2 Rm :

31

ð28Þ

In the above definition ku and ks are eigenvalues of, respectively, Ln and Lnþ. , since they characterise the linear dynamics around w0 ðnÞ and w0 ðn þ .Þ, i.e.

ku 2 spðLn Þ;

ks 2 spðLnþ. Þ:

ð29Þ

Then, condition (27) means that, in order for such a front to exist, the spectrum at w0 ðnÞ has to comprise at least a positive real part eigenvalue, i.e. w0 ðnÞ cannot be spectrally stable. This property reflects the need for a set of motions around W0 which are more complex than those featuring spectral stability. In this latter case, indeed, the perturbed trajectory is prevented from moving away from the nominal solution, and even a displacement along W0 is freezed, neither diverging nor converging. On the other hand, orbital stability allows a trajectory to move away from the local region where the nominal solution is settled and to reach a different part of W0 . Hence, around w0 ðnÞ a diverging path has necessarily to be recorded. Remark 9. In presence of fronts defined as in Definition 4.5, the spectrum of the operator Ln must comprise positive real part eigenvalues for some values of n. Thus, spectral instability does not imply orbital instability, but orbital stability is consistent with spectral instability. Orbital stability is a wider concept than spectral stability, and, despite being a weaker kind of stability, it turns out being pretty common in real processes. For instance, it is a common experience that in many media, such as water, sand, clouds or fluttering fabrics, the perturbation of a travelling wave may result in a phase shift of the original motion. Unfortunately, orbital stability is a very challenging problem, that requires to figure out any possible mode a perturbed trajectory may approach or diverge with respect to a set W0 of non solitary solutions. Such a study is beyond the scope of this paper and, besides, no general theory has been developed yet. Anyway, in the literature there exist some tailored results, developed for specific problems; for instance, examples of orbital stability analysis for travelling waves can be found in [33–35]. Hereafter, we illustrate through very simple insights a number of issues, which in general prevent the investigation of the orbital stability of W0 even with respect to the heteroclinic phenomenon of the connecting fronts. First of all, orbital stability is not a local problem. In fact, the perturbation is expected to move along W0 , and then the dynamics around the nominal solution w0 ðnÞ is not sufficient to characterise its behaviour. Therefore, the linearisation around W0 is not an efficient approach, since it is unable to completely describe the whole phenomenon, even in the case of a small front. Possibly, a study based on the Poincaré section technique might manage to grasp the orbital stability property, but the related linear map design would be unfeasible. Recently, a novel harmonic balance approach has been proposed in [36], but the results are only meant to provide an approximation of the period travelling wave. The linearisation approach fails also because the nominal solutions are non isolated, then preventing the local spectrum to provide information about the dynamics tangent to W0 , which in turn is the most interesting from the orbital stability point of view. This is a crucial complication in the problem, even though orbital stability of isolated waves can still be carried out only for specific systems and through tailored techniques, see e.g. [37–39]. Considering the Hilbert case p ¼ 2, a local positive real part eigenvalue at w0 ðnÞ denotes the existence of a diverging direction from this nominal solution. Indeed, remember that each eigenfunction is related to a specific local linear behaviour. However, since around w0 ðnÞ there are infinite other nominal solutions, a perturbation directed this way could actually gain a component parallel to W0 but with respect to some other w0 ðn þ .Þ nearby. In general, a diverging behaviour from the nominal solution w0 ðnÞ does not imply that the perturbation is distancing itself from the set W0 as a whole, but, rather, that its position is simply varying. Therefore, not only an unstable spectrum Ln does not denote orbital instability, but also it turns out to be necessary in order for the system to exhibit phenomena such as connecting fronts. The spectra of W0 are just not sufficient to study orbital stability. Finally, it is worth reporting that energy-like methods for analysing the stability of a specific class of systems have been proposed in the literature, see for example [40]. 5. Ring of coupled Fitzhugh–Nagumo systems: a practical tutorial The approach presented so far can be indifferently applied to a wide variety of behaviours exhibited by dynamical networks, ranging from convergence to homogeneous states and stationary fronts, to the propagation of solitary waves, pulses and other waves. In order to show how Theorem 3.1 can be exploited in a real example, we focus on the propagation of travelling waves in a homogeneous network of nonlinear nodes coupled to the nearest neighbours via linear connections and subject to periodic boundary conditions, i.e. a ring of linearly interconnected nodes. Efficient analytic methods and numerical techniques for the investigation of the properties of problem (14) will be presented, and direct numerical simulations of the discrete network dynamics will be used to quantitatively compare the real solution with the prediction obtained via the proposed PDE-based approach. Let us then consider a network composed by nodes governed by the Fitzhugh–Nagumo dynamics [41]

1 3

v_ j ¼ v j  wj  v 3j þ c þ dv jþ1 þ ev j1 ;

ð30Þ

_ j ¼ v j  awj þ g; w

ð31Þ

32

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

where c and g are external constant excitation currents. By defining two interpolating functions v ðx; tÞ; wðx; tÞ such that v j ðtÞ ¼ v ðxj ; tÞ and wj ðtÞ ¼ wðxj ; tÞ, and expanding v jþ1 and v j1 via Taylor series as in (6), we obtain v jþ1 ¼ v j þ e@ x v þ   , v j1 ¼ v j  e@ x v þ   , and similarly for w. We can then write the associated PDEs (10) for the interpolating functions as

@t v ¼ v  w 

1 3 v þ c þ ðd þ eÞv þ ðd  eÞe@ x v ; 3

ð32Þ

@ t w ¼ v  aw þ g;

ð33Þ

where we truncated the Taylor expansion to the first order term according to Theorem 3.1. By changing the coordinates ðx; tÞ # ðn ¼ kx  ct; tÞ we can rewrite (32) and (33) in the moving frame as

@t v ¼ v  w 

1 3 v þ c@ n v þ c þ ðd þ eÞv þ keðd  eÞ@ n v ; 3

ð34Þ

@ t w ¼ v  aw þ c@ n w þ g:

ð35Þ

A travelling wave propagating along the ring will then be a stationary solution of (34) and (35) satisfying

@nv ¼ 

@nw ¼

1 c þ keðd  eÞ



1 3

v  w  v 3 þ c þ ðd þ eÞv

ð36Þ

;

1 ðv þ aw  gÞ: c

ð37Þ

Eqs. (36) and (37) do not admit non-trivial bounded solutions for c > 0, whereas for c < 0 the dynamics converges to a stable limit cycle (on the moving coordinate n) corresponding to the propagation of a travelling wave along the ring. In Fig. 1 we show an example of such periodic solution, obtained by setting a ¼ 0:8; c ¼ g ¼ 0; d ¼ 0:27; e ¼ 0:1; N ¼ 50; e ¼ 0:05; c ¼ 2. We point out that the presence of periodic boundary conditions imposes additional constraints on the admissible values for the wave number k that must be tuned by imposing v ð0Þ ¼ v ðLÞ; wð0Þ ¼ wðLÞ where L is the length of the ring. For the results shown in Fig. 1 we used k ¼ 10:77, corresponding to having a single period of the travelling waves fitted in the whole length L ¼ 2:5 of the ring. In this section we will only focus on this situation for simplicity, but by choosing different values for k all the analysis can be easily performed for the generic case where multiple profiles fit within the length L. In Fig. 2 we show a comparison between the solution predicted by the reference ODEs (36) and (37) (shaded surface) and the direct simulation of the original network (30) and (31) (solid lines). We note that the reference ODE prediction is close to the real behaviour exhibited by the discrete network. To quantitatively assess the quality of the ODE predictions we also compared them with direct numerical simulations of (30) and (31) for various values of N, ranging from N ¼ 25 to N ¼ 100. Given that the network length is kept constant, this corresponds to gradually decreasing kC2 k in Theorem 3.1. As shown in Fig. 3, larger values of N induces smaller prediction errors, in agreement with the fact that the Taylor expansion used to define the PDE becomes more accurate as the state difference between adjacent nodes decreases.

1.5

v w

1

profile

0.5

0

−0.5

−1

−1.5 0

5

10

ξ

15

20

25

Fig. 1. Travelling wave profile obtained by solving the reference ODEs (36) and (37) with a ¼ 0:8; c ¼ g ¼ 0; d ¼ 0:27; e ¼ 0:1; N ¼ 50; L ¼ 2:5; c ¼ 2 and imposing periodic boundary conditions.

33

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

1.5 1 0.5 v

0

−0.5 −1 −1.5 50 40 30 20

t

10 0

0.5

0

1

2.5

2

1.5 x

Fig. 2. Comparison between predictions from reference ODEs (36) and (37) (shaded surface) and direct simulation of the network dynamics (30) and (31) with the same set of parameters used for Fig. 1.

30

amplitude wave speed c

Prediction error [%]

25

20

15

10

5

0

30

40

50

60

N

70

80

90

100

Fig. 3. Prediction error on wave amplitude and travelling speed c as a function of the number N of nodes placed in a ring of fixed length L ¼ 2:5.

Finally, we exploited the reference ODEs (36) and (37) to obtain the dispersion curve for the travelling waves propagating on the original network. In Fig. 4 we report a comparison between the estimated dispersion curve and the real dispersion curve calculated by simulating the discrete network (30) and (31). We note that in both cases the dispersion curve is linear. Moreover, the wave number k is always predicted with high accuracy, whereas the error on the travelling wave speed c translates to a difference in slope between the predicted and true dispersion curve. ^ WÞ ^ of (36) and (37) is known, its stability with respect to (at least infinitesimal) perturbations needs to Once a solution ðV; be analysed to check if such travelling wave can actually be sustained by the original network of dynamical systems. As described in Section 3, the first step in such direction is studying the so-called spectral stability of the solution, i.e. the computation of the spectrum of the operator Ln defined in Eq. (25). For this purpose, we consider a perturbation ðv ; wÞ around the ^ WÞ ^ and we make the ansatz v ðn; tÞ ¼ v ðnÞekt and wðn; tÞ ¼ wðnÞekt . By plugging such ansatz in (34) and nominal solution ðV; (35) and linearising the resulting equations, we obtain that the functions v ðnÞ and wðnÞ must satisfy the linear ODEs

^ 2 v þ ðd þ eÞv þ keðd  eÞ@ n v ;  c@ n v þ kv ¼ v  w  V  c@ n w þ kw ¼ v  aw; or, in a more compact form,



 " bðnÞþk v d v a ¼ Aðn; kÞ ¼ dn w w  1c bðnÞ ¼ 1 þ d þ e  V^ 2 ðnÞ;

1

a aþk c

#

v w

;

a ¼ keðd  eÞ þ c:

ð38Þ

ð39Þ

34

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

3.5

ODE network

3 2.5

c

2 1.5 1 0.5 0

0

2

4

6

8 k

10

12

14

16

Fig. 4. Comparison between predicted (solid) and true (dashed) dispersion curve for a travelling wave propagating on the ring.

Eq. (38), together with the periodic boundary conditions, represents a time-varying (the role of time being played by n here) boundary value problem whose analytical solutions cannot be easily obtained. We therefore discretize all the terms involving the derivative with respect to n using central finite difference over a grid of N mesh evenly spaced points on n 2 ½0; L, thus obtaining





ad v iþ1  v i1 ¼ kv i  bi v i þ wi ;

ð40Þ

  cd wiþ1  wi1 ¼ kwi  v i þ awi ;

ð41Þ

where d ¼ 1=2h; 

v

i

¼ v ðihÞ; w ¼ wðihÞ and h is the distance between two successive points on the grid [30,31]. By defining  i

v ¼ v 1 ; v 2 ; . . . ; v Nmesh ; w1 ; w2 ; . . . ; wNmesh we can write (40) and (41) in matrix form as Aeig v ¼ kv;

ð42Þ

where we have defined

2 " Aeig ¼

diagðbi Þ þ TðadÞ

I

I

diagðaÞ þ TðcdÞ

# ;

0

y

6 y 0 6 6 6 TðyÞ ¼ 6 0 y 6 . 6 . 4 . y

0

0

0

...

y

0

...

0

y

...

. . . 0 y

y

3

0 7 7 7 0 7 7: 7 7 5 0

Eq. (42) is the discrete counterpart of the general form (25) and we can obtain an approximation of the desired spectrum by evaluating the eigenvalues of Aeig . In Fig. 5 we report the approximated spectrum obtained using N mesh ¼ 400 points for the solution shown in Fig. 1. We note that the system has both stable (green) and unstable (red) eigenvalues and therefore the travelling wave is not spectrally stable, that is there exists perturbations that are not homogeneously reabsorbed on the nominal profile. However, as pointed out in Section 3, spectral stability is only a sufficient condition for the more general orbital stability, because, for example, a perturbation may not decay exponentially toward the nominal profile, but the perturbed solution may still converge to a shifted profile. To support this claim, in Fig. 6 we show the evolution of the network dynamics (30) and (31) when the initial conditions are initialized by using the nominal profile plus one of the unstable eigenvectors (in the specific case the eigenvector corresponding to k ¼ 0:02 þ 1:99i). Although the linearised dynamics would predict the response to be unstable, the perturbed solution actually converges to a shifted version of the nominal one. In fact, the functional distance

distðtÞ ¼

1

e

Z

L

½Vðx; tÞ  v ðx; tÞ2 dx;

ð43Þ

0

converges to a constant value. We conclude this section by pointing out that both finding the nominal solution and performing its stability analysis on the original network (having 2N ¼ 100 states) would have required a much larger computational effort than the one required by the proposed techniques (that effectively works on two-dimensional ODEs).

35

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

40 30 20

Im(λ)

10 0

−10 −20 −30 −40 −1

−0.8

−0.6

−0.4

−0.2

0 Re( λ)

0.2

0.4

0.6

0.8

1

Fig. 5. Spectrum of the linearised dynamics obtained by solving (42) with N mesh ¼ 400.

Y nominal start Y nominal end Y perturbed start Y perturbed end

2

1 v7

v

1 0

0

−1

−1 −2 0

nominal perturbed

2

0.5

1

x

1.5

2

−2 0

2.5

50

100 t

150

200

40

dist

30 20 10 0 0

50

100

150

200

250 t

300

350

400

450

500

Fig. 6. Response of the network to a perturbation to the initial profile oriented along a linearly unstable eigenvalue. Top left: initial and final profile of nominal and perturbed solution. Top right: temporal evolution of the seventh node in the network, showing that the perturbation get reabsorbed to the (shifted) nominal trajectory. Bottom panel: evolution of the functional distance distðtÞ defined in (43) showing that the distance between nominal and perturbed trajectory remains bounded.

6. Final remarks In this paper we have presented a technique for analysing the behaviour of networks of dynamical systems. An embedding in a (possibly) virtual spatial coordinate system allowed us to recast the study of several phenomena occurring on the network to the analysis of a related PDE whose solutions interpolate the dynamical behaviour exhibited by the nodes. Propagating solutions like fronts and travelling waves can be then obtained as stationary solutions of the PDE expressed in a moving coordinate n, thus effectively reducing the problem to the analysis of a nonlinear ODE. It is worth pointing out that the computational effort required to get such solution by directly simulating the whole network scales linearly with the number of nodes, whereas the burden required for solving the simplified ODE is not affected by the number of nodes. Moreover, the definition of the interpolating PDE allowed us to exploit some of the known techniques (such as spectral stability theory) for PDE stability analysis to study the stability of the corresponding solution on the original network. We also showed that spectral stability is not the most appropriate notion for studying the stability of travelling waves propagating on networks. An extensive example has been included to better illustrate the steps necessary to analyse the behaviour of a network composed by nodes governed by a non trivial nonlinear dynamics. The same example also highlights some of the benefits, such as reduced computational burden, that a simplified description of the network dynamics can provide, and it suggests that the proposed techniques represent a very effective tool for analysing the behaviour of networks composed by a large

36

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

number of nodes. Such promising results justify the exploration of extensions of the technique to more general frameworks and to the design of distributed feedback control systems, potentially opening new opportunities in such field. Appendix A. Theorems, lemmas and propositions Lemma Appendix A.1. Let X; Y be Banach spaces, U # X open and connected, a 2 ð0; 1Þ, and f : X ! Y a continuous operator. Moreover, define the operator F : U ! Y as

FðxÞ ¼ x þ f ðxÞ; and assume that f satisfies

kf ðx1 Þ  f ðx2 ÞkY 6 akx1  x2 kX ;

8x1 ; x2 2 U:

ðA:1Þ

Then, the image of U through F, namely FðUÞ, is open and F : U ! HðUÞ is a homeomorphism, i.e. F is continuous and invertible on U and it inverse F 1 : FðUÞ ! U is continuous on FðUÞ. Proof. See [24, chapt. 22]. h Corollary Appendix A.2. If f is linear and its operator norm kf k < 1, then condition (A.1) is satisfied for U ¼ X. Proof. By the linearity property, it is sufficient to observe that (A.1) is equivalent to

kf ðx1 Þ  f ðx2 ÞkY kf ðx1  x2 ÞkY kf ðxÞkY ¼ sup ¼ sup 6 kf k 6 a; x  x x  x k k k k 1 2 X 1 2 X x1 ;x2 2U x1 ;x2 2U x2Ux2 kxkX sup

8U # X:



Theorem Appendix A.3 (Inverse Function Theorem). Consider three Banach spaces, namely X; Y and Z, and a connected open set U  X  Y. Let H : X  Y ! Z be an operator and assume that H 2 C 1 ðUÞ and also that its derivative DHðx; yÞ is invertible for all ðx; yÞ 2 U. Moreover, let V 0 2 U be a connected open subset, ðx0 ; y0 Þ 2 V 0 and

    sup DH1 ðx0 ; y0 ÞDHðx; yÞ  Id < 1;

ðx;yÞ2V 0

where V 0 is the closure of V 0 and Id the identity operator. Then, H is invertible in V 0 ; H maps open subsets of V 0 in open subsets of  Z, and the restriction of H to V 0 , i.e. HjV 0 : V 0 ! HðV 0 Þ, is invertible and its inverse operator H1  : HðV 0 Þ ! V 0 is C 1 ðHðV 0 ÞÞ. V0

Proof. See [24, Lemma 22.28 and Theorem 22.29]. h Lemma Appendix A.4. Let f 2 Hk;p ðRm ; Rn Þ, with k P 0, and H : Hk;p ðRm ; Rn Þ ! Hh;p ðRm ; Rn Þ, with h 6 k. If H is a linear operator, it is (Lipschitz) continuous if and only if it is bounded. Proof. Assume that H is continuous. Then, for each

 > 0 it exists d 2 Rþ

such that

  Hðf  0 k;p Þ h;p ¼ kHf k h;p 6 ; H H H   for all f such that f  0Hk;p Hk;p ¼ kf kHk;p 6 d. Then, if we choose

kHf kHh;p

 ¼ 1, it follows that for each f we also have

      kf kHk;p   kf kHk;p  kf kHk;p d 1  H d f  ¼ H f ¼ 6  1 ¼ kf kHk;p ;  d   d d d kf kHk;p kf kHk;p Hh;p Hh;p

i.e. 1=d is an upper bound of the operator norm of H. Now, assume that H is bounded and that q > 0 is its norm. Then, for each df 2 Hk;p ðRm ; Rn Þ. such that df – 0Hk;p , we have

kHðf þ df Þ  Hf kHh;p ¼ kHdf kHh;p 6 qkdf kHk;p ;

8f :

Therefore, letting df ! 0Hk;p it follows that

lim kHðf þ df Þ  Hf kHh;p ¼ 0Hh;p ;

df !0Hk;p

and, since M does not depend on f, we also have that H is Lipschitz continuous. h

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

37

Lemma Appendix A.5. Consider an operator H : Hk;p ðRm ; Rn Þ ! Hh;p ðRm ; Rn Þ, with 0 6 h 6 k. If H is a linear bounded operator, then it is also continuously (Frétchet) differentiable. Proof. By the definition, H is differentiable at a certain f 2 Hk;p ðRm ; Rn Þ, if it exists an open set U  Hk;p ðRm ; Rn Þ containing the origin 0Hk;p , and a linear bounded operator DH : Hk;p ðRm ; Rn Þ ! Hh;p ðRm ; Rn Þ, such that

Hðf þ df Þ ¼ Hf þ DHdf þ oðdf Þ; for each f 2 U. Equivalently, the differentiability property can be rewritten as

lim

df !0Hk;p

kHðf þ df Þ  Hf  DHdf kHh;p ¼ 0: kdf kHk;p

Then, choose DH ¼ H and observe that this latter is bounded by Lemma Appendix A.4. Moreover, it follows that

kHðf þ df Þ  Hf  Hdf kHh;p ¼ 0 8df ; kdf kHk;p and for each f. Since DH does not depend on the chosen f where we are evaluating the differentiability property, it follows that the derivative of H is constant, and then continuous, for each f. h Lemma Appendix A.6. Consider f 2 Hkþi;p ðRm ; Rn Þ, with k > 0 and i P 0, and let h 2 Nm 0 be a multi-index such that jhj ¼ k. Then, the operator Dh : Hkþi;p ðRm ; Rn Þ ! Hi;p ðRm ; Rn Þ defined as

Dh ðf Þ ¼ @ hx f ; is linear and bounded. Proof. The linearity property directly derives from the definition of (weak) derivative. Then, since f belongs to the Sobolev space Hkþi;p ðRm ; Rn Þ, we have

   j  @ x f 

H0;p

< 1;

m for each multi-index j 2 Nm 0 such that jjj 6 k þ i. Moreover, for each j 2 N0 such that jjj 6 i, observe that

@ jx Dh ðf Þ ¼ @ xhþj f ; with jh þ jj ¼ k þ jjj 6 k þ i. Hence, it follows that

   h  D f 

Hi;p

¼

 X  j h  @ x D ðf Þ jjj6i

H0;p

¼

 X  hþj  @ x f  jjj6i

H0;p

¼

X jjj 6 k þ i j  h 2 Nm 0

   j  @ x f 

H0;p

< 1:



Lemma Appendix A.7. Let f 2 Hkþi;p ðRm ; Rn Þ, with k P 2 and i P 0, and u 2 Rm , such that u – 0Rm . Then, the operator Hku , defined as

Hku ðf ÞðxÞ ¼

X uh Z h! jhj¼k

1 0

ð1  aÞ@ hx f ðx þ auÞda;

h 2 Nm 0;

^ k : Hkþi;p ðRm ; Rn Þ ! Hi;p ðRm ; Rn Þ is linear and bounded. is linear in Hkþi;p ðRm ; Rn Þ. Moreover, the restriction H u Proof. The linearity property directly follows from the properties of the sum, the (weak) derivative and the integral. First, observe that k X uh h 1 d f ðx þ auÞ ¼ @ f ðx þ auÞ: k! dak h! x jhj¼k

Then, by integrating by parts, we have

38

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

^ k ðf ÞðxÞ ¼ H ^ k f ðxÞ ¼ H u u

Z

1

ð1  aÞ 0

X uh

jhj¼k

h!

"

@ hx f ðx þ auÞda ¼

Z

k

1

ð1  aÞ

0

1 d f ðx þ auÞda k! dak

#1

Z 1 k1 k1 1 d 1 d f ðx þ auÞ þ f ðx þ auÞda k1 k1 k! da 0 k! da 0 " #1 k2 1 X uh h 1 d ¼ @ x f ðxÞ þ f ðx þ a uÞ k jhj¼k1 h! k! dak2 ¼ ð1  aÞ

0

X uh h 1 X uh h 1 ¼ @ x f ðxÞ þ @ x f ðx þ uÞ  @ hx f ðxÞ : k jhj¼k1 h! kðk  1Þ jhj¼k2 h!

ðA:2Þ

^ k and Lemma Appendix A.6 the statement is proved. h Thus, by exploiting the linearity of H u Proposition Appendix A.8. The operator s : Hk;p ðRm ; Rn Þ  Rm ! Hk;p ðRm ; Rn Þ, with k > 0, defined as in (12) is continuous and continuously differentiable with respect to both its arguments. Proof. For each given

. 2 Rm the operator s is linear by definition, indeed

sðwðt; nÞ þ dwðt; nÞ; .Þ ¼ wðt; n þ .Þ þ dwðt; n þ .Þ: Hence, it follows that

lim ksðw þ dw; .Þ  sðw; .ÞkHk;p ¼

dw!0Hk;p

lim ksðdw; .ÞkHk;p ¼ lim kdwkHk;p ¼ 0;

dw!0Hk;p

dw!0Hk;p

i.e. sðw; .Þ is continuous and then bounded in H ðR ; R Þ at each w. Thus, by Lemma Appendix A.5 it follows that for each fixed . the operator s is linear, bounded and continuously differentiable with respect to w. Then, since w 2 Hk;p ðRm ; Rn Þ we have k;p

m

sðwðt; nÞ; . þ d.Þ ¼ wðt; n þ . þ d.Þ ¼ wðt; n þ .Þ þ ¼ sðwðt; nÞ; .Þ þ

n

d wðt; n þ .Þd. þ oðd.Þ dn

d sðwðt; nÞ; .Þd. þ oðd.Þ: dn

ðA:3Þ

Hence, for each given w we have

   d lim ksðw; . þ d.Þ  sðw; .ÞkHk;p ¼ lim  sðw; .Þd. þ oðd.Þ   d.!0 m d.!0 m dn R

¼ 0;

Hk;p

R

since for k > 0 the derivative is a bounded operator, as a direct consequence of Lemma Appendix A.6. This implies that s is continuous at each .. Moreover, observe that

lim

d.!0Rm

    d sðw; .Þd. sðw; . þ d.Þ  sðw; .Þ  dn kd.kRm

Hk;p

¼ 0:

Thus, s is differentiable at each ., and, since ðd=dnÞsðw; .Þ is the combination of continuous operators, by the previous lemmas we also have that s is continuously differentiable. Finally, by applying the chain rule and the above observations, we have

sðw þ dw; . þ d.Þ ¼ sðw; .Þ þ

d sðw; .Þd. þ sðdw; .Þ þ oðdw; d.Þ ¼ sðw; .Þ þ Dsðw; .Þðdw; d.Þ þ oðdw; d.Þ; dn

ðA:4Þ

which implies, by the previous results on the continuity property of the single partial operators, that s is continuously differentiable with respect to both its arguments. h References [1] Roska T, Chua L. The CNN universal machine: an analogic array computer. IEEE Trans Circuits Syst II: Analogue and Digital Signal Process 1993;40:163–73. [2] Nishio T, Nishio Y. Periodic pattern formation and its applications in cellular neural networks. IEEE Trans Circuits Syst I: Regul Pap 2008;55:2736–42. [3] Chen E, Stiefel K, Sejnowski T, Bullock T. Model of travelling waves in a coral nerve network. J Comp Physiol A: Neuroethol Sensory Neural Behav Physiol 2008;194:195–200. [4] Alexander RM. Principles of animal locomotion. Princeton University Press; 2006. [5] Han S-G, Um J, Kim B. Voter model on a directed network: role of bidirectional opinion exchanges. Phys Rev E 2010;81:057103. [6] Christensen C, Albert I, Grenfell B, Albert R. Disease dynamics in a dynamic social network. Physica A 2010;389:2663–74.

G. Innocenti, P. Paoletti / Commun Nonlinear Sci Numer Simulat 24 (2015) 21–39

39

[7] Sarlette A, Sepulchre R. A PDE viewpoint on basic properties of coordination algorithms with symmetries. Proceedings of 48th IEEE conference on decision and control. p. 5139–44. [8] Innocenti G, Paoletti P. Traveling waves propagation on networks of dynamical systems,. Proceedings of 19th IFAC world congress. p. 8445–50. [9] Innocenti G, Paoletti P. A virtual space embedding for the analysis of dynamical networks. Proceedings of 52nd IEEE conference on decision and control. p. 1331–6. [10] Gilli M, Roska T, Chua L, Civalleri P. On the relationship between CNNs and PDEs, CNNA 2002. Proceedings of the 2002 7th IEEE international workshop on cellular neural networks and their applications, 2002.. p. 16–24. [11] Barooah P, Mehta P, Hespanha J. Mistuning-based control design to improve closed-loop stability margin of vehicular platoons. IEEE Trans Autom Control 2009;54:2100–13. [12] Meurer T, Krstic M. Finite-time multi-agent deployment: a nonlinear PDE motion planning approach. Automatica 2011;47:2534–42. [13] Motee N, Jadbabaie A. Optimal control of spatially distributed systems. IEEE Trans Autom Control 2008;53:1616–29. [14] Ren W. Multi-vehicle consensus with a time-varying reference state. Syst Control Lett 2007;56:474–83. [15] Sepulchre R. Consensus on nonlinear spaces. Annu Rev Control 2011;35:56–64. [16] Liu K, Xie G, Ren W, Wang L. Consensus for multi-agent systems with inherent nonlinear dynamics under directed topologies. Syst Control Lett 2013;62:152–62. [17] Wen G, Duan Z, Yu W, Chen G. Consensus of multi-agent systems with nonlinear dynamics and sampled-data information: a delayed-input approach. Int J Robust Nonlinear Control 2013;23:602–19. [18] Angeli D, Sontag ED. Translation-invariant monotone systems, and a global convergence result for enzymatic futile cycles. Nonlinear Anal: Real World Appl 2008;9:128–40. [19] Baskar L, De Schutter B, Hellendoorn J, Papp Z. Traffic control and intelligent vehicle highway systems: a survey. IET Intel Transport Syst 2011;5:38–52. [20] Kari J. Theory of cellular automata: a survey. Theor Comput Sci 2005:3–33. [21] Ostoja-Starzewski M. Lattice models in micromechanics. Appl Mech Rev 2002;55:35–60. [22] Giamarchi T. Quantum physics in one dimension. Oxford University Press; 2004. [23] Johnson WP. The curious history of Faà di Bruno’s formula. Am Math Monthly 2002;109:217–34. [24] Driver BK. Analysis tools with applications. Lecture notes of the course in Partial Differential Equations; 2002. [25] Grillakis M, Shatah J, Strauss W. Stability theory of solitary waves in the presence of symmetry I. J Funct Anal 1987;74:160–97. [26] Grillakis M, Shatah J, Strauss W. Stability theory of solitary waves in the presence of symmetry II. J Funct Anal 1990;94:308–48. [27] Sandstede B. Stability of travelling waves. Handbook of dynamical systems, vol. 2. Elsevier; 2002. p. 983–1055. [28] Edmunds D, Evans W. Spectral theory and differential operators. Oxford University Press; 1987. [29] Deconinck B, Kutz J. Computing spectra of linear operators using the Floquet–Fourier–Hill method. J Comput Phys 2006;219:296–321. [30] Rademacher J, Sandstede B, Scheel A. Computing absolute and essential spectra using continuation. Physica D 2007;229:166–83. [31] Guba O, Lorenz J. Continuous spectra and numerical eigenvalues. Math Comput Modell 2011;54:2616–22. [32] Sherratt JA. Numerical continuation methods for studying periodic travelling wave (wavetrain) solutions of partial differential equations. Appl Math Comput 2012;218:4684–94. [33] Jones CKRT. Stability of the travelling wave solution of the Fitzhugh–Nagumo system. Trans Am Math Soc 1984;286:431–69. [34] Gallay T, Ha˘ra˘gus M. Orbital stability of periodic waves for the nonlinear Schrödinger equation. J Dyn Differ Equ 2007;19:825–65. [35] Pava JA, Natali FMA. Stability and instability of periodic travelling wave solutions for the critical Korteweg–de Vries and nonlinear Schrödinger equations. Physica D 2009;238:603–21. [36] Leung A, Yang H, Guo Z. Periodic wave solutions of coupled integrable dispersionless equations by residue harmonic balance. Commun Nonlinear Sci Numer Simul 2012;17:4508–14. [37] Leto J, Choudhury S. Solitary wave families of a generalized microstructure PDE. Commun Nonlinear Sci Numer Simul 2009;14:1999–2005. [38] Zhang X, Xu R. Global stability and travelling waves of a predator-prey model with diffusion and nonlocal maturation delay. Commun Nonlinear Sci Numer Simul 2010;15:3390–401. [39] Zhao H-Q, Weng P, Liu S-Y. Uniqueness and stability of travelling wave fronts for an age-structured population model in a 2d lattice strip. Commun Nonlinear Sci Numer Simul 2014;19:507–19. [40] Iooss G, Kirchgässner K. Travelling waves in a chain of coupled nonlinear oscillators. Commun Math Phys 2000;211:439–64. [41] FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane. Biophys J 1961;1:445–66.