Marine Coralie Basile

Gwezheneg, Benoît, Yvan, Vincent, Damien Hélène Tristan

Cyril

Valentin Olivier

Pierig Le Meur

Clément. Federico

i Je tiens à remercier Philippe et Franc ¸ois, mes directeurs de thèse, pour l'aide qu'ils m'ont apportée tout au long de ces trois années. Ils ont été présents au quotidien pour me pousser à avancer, réfléchir avec moi et corriger mes nombreuses erreurs. Ils ont su m'aider à résoudre les problèmes qui se posaient tout en me laissant le choix sur la direction à donner à la thèse. Je les remercie d'avoir supporté mon pessimisme persistant et mon humeur fluctuante. Je remercie Franc ¸ois pour les conversations chiffons et potins, Philippe pour sa disponibilité, les thés offerts et surtout pour les conférences IPSO ! Je les remercie tous les deux pour les conseils cinéma, les digressions toujours inattendues, et généralement pour avoir su me redonner le sourire quand j'étais déprimée ! Je les remercie aussi d'avoir été compréhensifs face à mon choix professionnel (même si Philippe ronchonne pour le postdoc) et d'avoir tout fait pour que je puisse soutenir dans les temps. Je remercie les rapporteurs Christophe Besse et Elena Celledoni d'avoir pris de temps de relire ma thèse en détail. Je remercie également Nicolas Crouseilles, Brynjulf Owren et Gilles Vilmart d'avoir accepté de faire partie de mon jury. Je remercie les relecteurs bénévoles de mon introduction de thèse : Hélène Hivert, Marine Malo et Loïc Le Treust. Merci d'y avoir consacré du temps ! Puisque j'en suis à la phase professionnelle, je remercie Éric Darrigrand, Fabrice Mahe et Max Bauer pour ou avec lesquels j'ai encadré des TD pour leur disponibilité et leur empressement à répondre à mes questions. Je remercie aussi Carole, Hélène (Rousseaux), Marie-Aude, Chantale, Marie-Annick et toutes les secrétaires pour leur aide administrative au quotidien. Je remercie Cécile et Stéphanie de

Introduction

INTRODUCTION

L'objet de cette thèse est l'étude d'un système différentiel du type :

   ẋ = F x, z, t ε , x(0) = x 0 ∈ R n ż = 1 ε B t ε z + G x, z, t ε , z(0) = z 0 ∈ R m
, (0.0.1) avec ε un petit paramètre, B une matrice dont la résolvante est exponentiellement décroissante, F et G des fonctions régulières1 et toutes les fonctions de t ε périodiques en cette variable. On veut décrire le comportement de la solution dans la limite ε tend vers 0.

Deux dynamiques apparaissent dans ce système : une dynamique de convergence vers une variété invariante (due au terme raide 1 ε ) et une dynamique rapidement oscillante (due aux dépendances en t ε ). L'étude de la première suppose l'utilisation de techniques de type variété centrale, alors que la seconde fait appel à des méthodes de moyennisation. Notons que ces deux caractéristiques du système le rendent difficile à résoudre numériquement : la dynamique temporelle oscillant à l'échelle de temps ε, de même que la dérivée de z, il est nécessaire de choisir un pas de discrétisation ∆t de l'ordre de ε. Ici, on cherche à étudier la dynamique limite ε → 0, ce qui rend cette approche très coûteuse. On va ramener l'étude de (0.0.1) à celle d'un système se prêtant mieux à une résolution numérique.

Une dynamique de variété centrale 0.1.1 Théorie des variétés centrales

Dans un premier temps, nous nous sommes focalisés sur l'aspect "variété centrale" de ce système et nous avons démontré l'existence d'une variété centrale rapidement oscillante périodique. Dans cet objectif, nous avons adapté certains résultats du livre de J. Carr [START_REF] Carr | Applications of Centre Manifold Theory[END_REF] qui pose les bases de la théorie des variétés centrales. Pour mettre en avant les enjeux de ce domaine, nous présentons ici le contexte et les résultats de [START_REF] Carr | Applications of Centre Manifold Theory[END_REF], ainsi que différentes généralisations ( [START_REF] Mielke | A reduction principle for non autonomous systems in infinite dimensional spaces[END_REF], [START_REF] Sakamoto | Invariant manifolds in singular perturbation problems for ordinary differential equations[END_REF], [START_REF] Aulbrach | Integral manifold for Caratheodory type differential equations in Banach spaces[END_REF]). Nous expliciterons ensuite les résultats que nous avons obtenus et la direction dans laquelle ils nous ont menés.

Une variété centrale au voisinage d'un point d'équilibre en dimension finie

Dans [START_REF] Carr | Applications of Centre Manifold Theory[END_REF], les systèmes différentiels de la forme ẋ = Ax + f (x, z) , x(0) = x 0 ∈ R n ż = Bz + g(x, z), z(0) = z 0 ∈ R m , (0.1.1) sont étudiés au voisinage du point d'équilibre (0, 0). Les matrices A et B sont des matrices constantes possédant les propriétés spectrales suivantes : les valeurs propres de A sont toutes de partie réelle nulle alors que celles de B sont toutes de partie réelle strictement négative. Les fonctions f et g sont de classe C 2 , avec f (0, 0) = 0, f ′ (0, 0) = 0, g(0, 0) = 0 et g ′ (0, 0) = 0. Le système (0.0.1) pourrait s'y ramener en prenant A = 0, mais on ne peut pas choisir B = -1 ε B t ε , puisque B ne peut pas dépendre du temps. La présence du paramètre ε dans le système (0.0.1) introduit une autre différence avec (0.1.1), puisqu'elle permet de ne pas faire d'hypothèses sur F (0, 0) et G(0, 0). Définition 0.1.1 Une variété (x, h(x)) est invariante pour le système (0.1.1) lorsque z 0 = h(x 0 ) entraîne ∀t ∈ R, z(t) = h(x(t)). Définition 0.1.2 Si (x, h(x)) est une variété invariante pour le système (0.1.1), avec h une fonction régulière telle que h(0) = 0 et h ′ (0) = 0, alors h est une variété centrale associée à (0.1.1).

Grâce à un théorème de point fixe, on montre l'existence d'une telle variété au voisinage du point d'équilibre 0. Théorème 0.1.3 Il existe une variété centrale pour (0.1.1) : z = h(x) pour |x| < δ, avec h ∈ C 2 (R n ).

xi Idée de preuve: La preuve se décompose en plusieurs étapes :

1. Soit ε > 0, on utilise une fonction de troncature pour localiser l'étude sur B(0, ε). Concrètement, des fonctions F et G régulières qui coïncident sur B(0, ε) avec f et g, nulles en dehors de B(0, 2ε) sont introduites. On travaille sur le système :

ẋ = Ax + F (x, z) , x(0) = x 0 ∈ R n ż = Bz + G(x, z), z(0) = z 0 ∈ R m . (0.1.2)
Étant donné que le résultat cherché est local en x, ce n'est pas une restriction.

2. Soient p > 0 et p 1 > 0, on définit l'espace fonctionnel X = {h : R n → R m lipschitzienne, de constante de lipschitz p 1 , bornée par p, h(0) = 0}.

Muni de la norme sup adaptée (sup sur h et sur D x h), X est un espace complet.

3. Soient x 0 ∈ R n et h ∈ X fixés. On définit x(s, x 0 , h) comme la solution du système différentiel :

ẋ = Ax + F (x, h(x)) , x(0, x 0 , h) = x 0 ∈ R n .
4. On cherche à résoudre ż(t) = B z(t) + G(x(t, x 0 , h), h(x(t, x 0 , h))).

La formule de Duhamel donne z(t) = z(t 0 )e B(t-t0) + t t0 e B(t-s) G(x(s, x 0 , h), h(x(s, x 0 , h)))ds.

(0.1.3)

On cherche une solution z bornée sur R et B a toutes ses valeurs propres de partie réelle strictement négative, donc la limite t 0 → -∞ dans (0.1.3) donne

z(t) = t -∞
e B(t-s) G(x(s, x 0 , h), h(x(s, x 0 , h)))ds.

Pour t = 0, on obtient

z(0) = 0 -∞
e -Bs G(x(s, x 0 , h), h(x(s, x 0 , h)))ds = h(x 0 ). (0.1.4) Ainsi, si h est une variété centrale pour (0.1.2), l'équation (0.1.4) montre que h doit être un point fixe de l'opérateur

T : h → T h avec T h(x 0 ) = 0 -∞
e -Bs G(x(s, x 0 , h), h(x(s, x 0 , h)))ds.

On montre qu'inversement tout point fixe de T est une variété centrale pour (0.1.2).

5. On montre que pour p, p 1 et ε assez petits T : X → X.

6. On montre que pour p, p 1 et ε assez petits T est une contraction. On obtient alors l'existence d'un unique point fixe h ∈ X.

7. On montre la régularité de h.

Remarque 0.1.4 On n'a pas unicité des variétés centrales. Considérons l'exemple suivant, avec n = m = 1:

ẋ = x 2 , x(0) = x 0 ∈ R ż = -z, z(0) = z 0 ∈ R .

INTRODUCTION

Ce système admet comme solution explicite:

∀t ∈ 0, 1 x 0 , x(t) = x 0 1 -tx 0 , z(t) = z 0 e -t .
Pour toute constante C ∈ R, la fonction h définie par h(x) = Ce 1

x pour x = 0 et h(0) = 0 est une variété invariante. En effet, h(x(t)) = Ce x 0 e -t , donc si on part de z 0 = h(x 0 ) = Ce 1 x 0 , on a z(t) = h(x(t)).

Ainsi, la solution de (0.1.1) avec condition initiale (x 0 , h(x 0 )) est (x h (t), h(x h (t))), où x h est régi par une équation différentielle ordinaire: ẋh = Ax h + f (x h , h(x h )), (0.1.5) qui n'est autre que le flot sur la variété centrale. Que peut-on dire dans le cas où la condition initiale (x 0 , z 0 ) n'est pas sur la variété centrale (i.e. (x 0 , z 0 ) = (x 0 , h(x 0 ))) ?

Théorème 0.1.5 1. On suppose que la solution nulle de (0.1.5) est stable (resp. asymptotiquement stable, resp. instable). Alors la solution nulle de (0.1.1) est stable (resp. asymptotiquement stable, resp. instable).

2. On suppose que la solution nulle de (0.1.5) est stable. Soit (x(t), z(t)) solution de (0.1.1) avec x 0 et z 0 suffisamment petits. Alors il existe une solution x h (t) de (0.1.5), il existe µ > 0 tels que pour t → +∞, on ait : Le Théorème (0.1.5) se prouve à partir du Lemme (0.1.6) en utilisant un théorème de point fixe.

x(t) = x h (t) + O (e -µt ) , z (t)
L'égalité (0.1.6) montre que le comportement asymptotique de la solution de (0.1.1) est la dynamique sur la variété, d'où l'intérêt de se ramener à l'étude d'une variété centrale.

La preuve d'existence d'une variété centrale ne donne pas d'expression explicite. On cherche alors à approcher h. On dispose d'une information : h vérifie l'équation aux dérivées partielles h ′ (x)[Ax + f (x, h(x))] = Bh(x) + g(x, h(x)).

(0.1.7) où h ′ (x) = Dh(x) la jacobienne de h. On le montre en dérivant z(t) = h(x(t)) et en utilisant la deuxième équation de (0.1.1). Une fonction solution h de l'équation aux dérivées partielles (0.1.7) et vérifiant les conditions h(0) = 0 et h ′ (0) = 0 est une variété centrale pour le système (0.1.1). La résolution de (0.1.7) n'est pas plus simple que celle du système initial, mais c'est elle qui permet d'approcher h.

Définition 0.1.7 Soit V 0 un voisinage de l'origine dans R n . Pour toute application Φ : V 0 → R m de classe C 1 , on définit M Φ(x) = Φ ′ (x)[Ax + f (x, Φ(x))] -BΦ(x)g(x, Φ(x)).

Théorème 0.1.8 Soit Φ une application de classe C 1 d'un voisinage de l'origine de R n dans R m , avec Φ(0) = 0 et Φ ′ (0) = 0. On suppose que M Φ(x) = |x|→0 O (|x| q ) avec q > 1. Alors

|h(x) -Φ(x)| = |x|→0 O (|x| q ) .

UNE DYNAMIQUE DE VARI ÉT É CENTRALE xiii

Idée de preuve:

1. Comme dans la preuve du Théorème 0.1.3, on travaille avec des fonctions localisées sur une boule au voisinage de |x| = 0.

2. On considère l'opérateur T défini dans la preuve du Théorème 0.1.3 et on introduit :

S : z → T (z + Φ) -Φ.
Comme T , l'application S est une contraction.

3. Soit K > 0, on pose Y = {z ∈ X, |z(x)| ≤ K|x| q } et on montre qu'il existe un K > 0 tel que :

S : Y → Y.
4. Pour z = 0 ∈ Y , on obtient alors :

|Sz(x)| = |T (Φ)(x) -Φ(x)| ≤ K|x| q .
Ainsi, on a successivement :

|h(x) -Φ(x)| = |T h(x) -Φ(x)| car h est le point fixe de T ≤ |T h(x) -T Φ(x)| + |T Φ(x) -Φ(x)| ≤ c|h(x) -Φ(x)| + |T Φ(x) -Φ(x)| où 0 < c < 1 est la constante de contraction de T ≤ 1 1 -c |T Φ(x) -Φ(x)| ≤ K 1 -c |x| q .
Une variété centrale au voisinage d'un point d'équilibre en dimension infinie

Toujours dans [START_REF] Carr | Applications of Centre Manifold Theory[END_REF], J. Carr démontre les mêmes résultats pour un système de dimension infinie. Soit X un espace de Banach, on considère le système ẇ = Cw + N (w), w(0) ∈ X. (0.1.8)

On suppose que N : X → X est C 2 et que sa dérivée seconde est uniformément continue, avec N (0) = 0 et N ′ (0) = 0. De plus, C est le générateur d'un semi-groupe S(t) fortement continu sur X, et on suppose qu'il a les propriétés spectrales suivantes :

-X = V ⊕ Y où V est de dimension finie et Y est fermé, -V est C-invariant et si on note A la restriction de C à V , alors toutes les valeurs propres de A sont de partie réelle nulle, -En notant U (t) la restriction de S(t) à Y , Y est U (t)-invariant et ∃c 1 , µ > 0, ∀t ≥ 0, U (t) ≤ c 1 e -µt .

Dans ce contexte, une variété centrale est définie comme une variété invariante pour (0.1.8) qui est tangente à V en l'origine.

De très nombreux articles prolongent et étendent ces premiers résultats dans des contextes variés. Nous citons quelques uns de ces résultats, sans volonté d'exhaustivité. Le livre de Tony Roberts [START_REF] Roberts | Model Emergent Dynamics in Complex Systems[END_REF] présente un large panel d'applications des techniques de variétés centrales pour les sciences appliquées.

INTRODUCTION

Une variété centrale pour une équation non autonome

Dans [START_REF] Mielke | A reduction principle for non autonomous systems in infinite dimensional spaces[END_REF], A. Mielke étudie les solutions bornées d'équations différentielles non autonomes du type ẋ -Lx = f (t, λ, x) dans un espace de Banach infini X, avec t ∈ R et λ ∈ Λ un ouvert de R n . L'opérateur L est linéaire non borné et possède les propriétés suivantes :

-X = X 1 × X 2 avec X 1 de dimension finie et la restriction de L à X 1 ne possède que des valeurs propres imaginaires pures.

-La restriction de L à X 2 a une résolvante exponentiellement décroissante.

De plus, f est régulière telle que (λ 0 ∈ Λ étant fixé) pour tout t ∈ R, f (t, λ 0 , 0) = 0 et ∂ x f (t, λ 0 , 0) = 0. Dans un premier temps, l'existence d'une variété centrale de dimension finie est démontrée. Ensuite, la transmission de la périodicité (ou presque-périodicité) de f par rapport à t à la variété centrale est prouvée.

Une variété centrale au voisinage d'une courbe d'équilibre en dimension finie

Dans [START_REF] Sakamoto | Invariant manifolds in singular perturbation problems for ordinary differential equations[END_REF], la dynamique d'une variété centrale au voisinage d'une courbe d'équilibre est étudiée par K. Sakamoto.

Le système différentiel est le suivant :

ẋ = εf (x, z, ε) , x(0) = x 0 ∈ R n ż = g(x, z, ε), z(0) = z 0 ∈ R m , (0.1.9)
avec les hypothèses :

-Il existe r ∈ N * tel que f et g soient C r -bornées en tant que fonctions de (x, y, ε) et il existe h(x) une fonction aux dérivées bornées jusqu'à l'ordre r (sauf la fonction elle-même) telle que ∃ε 0 , ∀x ∈ R n , ∀0 < ε < ε 0 , g(x, h(x), ε) = 0.

(x, h(x), 0) est alors une courbe d'équilibre pour (0.1.9), puisque pour ε = 0, z(t) = h(x(t)) correspond à ż(t) = ẋ(t) 0 h(x(t)) = 0 = g(x(t), h(x(t)), 0)

-Soit µ ∈ R * + fixé. ∃k ∈ N, 0 ≤ k ≤ n tel que ∀(x, z) ∈ R n × R m
les matrices D z g(x, z, 0) ont k valeurs propres de partie réelle plus petite que -2µ, et nk valeurs propres de partie réelle plus grande que 2µ.

Alors, au voisinage de (x, h(x), ε), le système (0.1.9) admet une variété centrale (x, h ε (x), ε), telle que

h ε -h ∞ = ε→0 O(ε).
De plus, cette fonction h ε (x) peut être approchée à tout ordre en ε.

Une variété centrale pour des fonctions discontinues

Dans [START_REF] Aulbrach | Integral manifold for Caratheodory type differential equations in Banach spaces[END_REF], B. Aulbrach et T. Wanner ont généralisé l'étude de variétés centrales dans un espace de Banach infini X pour des systèmes différentiels non autonomes de la forme ẋ = A(t)x + f (t, x), au cas où les fonctions sont seulement mesurables en la variable t. Ces hypothèses de faible régularité permettent de traiter des discontinuités.

Une variété centrale périodique rapidement oscillante

Dans le cadre de cette thèse, nous avons généralisé les résultats de [START_REF] Carr | Applications of Centre Manifold Theory[END_REF] au cas d'une dynamique rapidement oscillante périodique. Dans le chapitre 1, nous adaptons la preuve d'existence de J. Carr pour démontrer l'existence d'une variété centrale périodique en t ε de la forme z(t) = h ε x(t), t ε pour le système différentiel (0.0.1).

De plus, pour R > 0 fixé, il existe ε 0 > 0 tel que pour tout ε < ε 0 , pour tout (x 0 , z 0 ) ∈ B(0, R) ⊂ R n × R m , on ait le résultat de convergence

∀t ≥ 0, z(t) -h ε x(t), t ε = O e -µ t ε .
La variété centrale nous permet donc de décrire la dynamique asymptotique de z(t). En revanche, si (x 0 , z 0 ) est quelconque, x(t) ne converge pas vers la solution de

ẋ = F x, h ε x, t ε , t ε , x(0) = x 0 ∈ R n .
En effet, comme le montre la Figure 1, la fonction x(t) ne représente pas le comportement asymptotique de x(t). Ce comportement n'est pas spécifique au cas d'une variété centrale rapidement oscillante, ce problème est déjà présent dans les travaux de J.Carr [START_REF] Carr | Applications of Centre Manifold Theory[END_REF].

INTRODUCTION

Pour capter cette dynamique asymptotique, il faut considérer x h (t) solution de ẋh = F x h , h ε x h , t ε , t ε , x h (0) = x ε 0 ∈ R n , (0.1.10) avec x ε 0 choisi comme le représente la Figure 2, en résolvant de T ∞ à 0 un système du type

ẋ = F x, h ε x, t ε , t ε , x(T ∞ ) = x(T ∞ ),
avec T ∞ > 0 un temps quelconque, puis en prenant x ε 0 = x(0). Cette technique est celle utilisée dans [START_REF] Sakamoto | Invariant manifolds in singular perturbation problems for ordinary differential equations[END_REF] par K.Sakamoto pour résoudre ce problème dans son contexte. Remarque 0.1.9 Notons que T ∞ = +∞ n'est pas nécessaire. En effet, on a l'approximation suivante :

∀t ∈ [0, T ∞ ], x(t) -x h (t) = O e -µ t ε ,
qui donne des informations pour t grand uniquement. Si on choisit T ∞ fini assez grand, mais fini, x(t) n'a pas encore parfaitement convergé vers sa dynamique limite et n'est donc pas exactement sur la variété. Cependant, l'erreur dûe à cet écart est incluse dans le O e -µ t ε de cette estimation.

La solution (x(t), z(t)) de (0.0.1) avec condition initiale (x 0 , z 0 ) converge exponentiellement rapidement vers x h (t), h ε x h (t), t ε avec condition initiale (x ε 0 , h ε (x ε 0 , 0)). La dynamique est donc celle illustrée par la Figure 3.

La preuve de notre version rapidement oscillante du théorème de variété centrale ne donne pas d'expression explicite d'une variété centrale. Cependant, h ε est solution de l'équation aux dérivées partielles :

1 ε ∂ θ h ε (x, θ) -B(θ)h ε (x, θ) = G(x, h ε (x, θ), θ) -∂ x h ε (x, θ)F (x, h ε (x, θ), θ) , (0.1.11)
et cela donne accès à une approximation de h ε . La fonction h

[r]

ε solution de (0.1.11) à l'ordre ε r est construite sous la forme d'un développement en puissance de ε : h

[r] ε = h 0 + εh 1 + ε 2 h 2 + • • • + ε r h r . C'est alors une approximation de h ε : h ε -h [r] ε ∞ = ε→0 O ε r+1 . 0.1. UNE DYNAMIQUE DE VARI ÉT É CENTRALE xvii z x z 0 x 0 x ε 0 h ε (x ε 0 , 0) (x(t), z(t)) x h (t), h ε x h (t), t ε Figure 3: La convergence de (x(t), z(t)) vers une variété centrale x h (t) , h ε x h (t), t ε . Notation 0.1.10 L'approximation d'ordre 0 en ε vaut h 0 = 0, donc h ε (x, θ) = O (ε) .
À partir de maintenant, on note la variété centrale εh(x, θ) pour mettre en évidence cette caractéristique. Il s'agit d'un abus de notations, puisque la fonction h(x, θ) dépend de ε.

La dynamique sur la variété est régie par l'équation différentielle (0.1.10), qui n'a pas le caractère raide que présentait le système (0.1.1), mais qui conserve le caractère hautement oscillant. Pour étudier (0.1.10), il faut gérer cette dynamique rapidement oscillante de la variété centrale. Des méthodes de moyennisation sont alors utilisées, qui permettent de caractériser la dynamique limite et ses perturbations d'ordre supérieur.

Moyennisation d'équations rapidement oscillantes

On cherche à approcher numériquement la solution d'une équation différentielle rapidement oscillante du type :

ẋ(t) = F ε x(t), t ε , x(0) = x 0 ∈ R n , (0.1.12) avec F ε (x, θ) une fonction régulière en (ε, x, θ) et T -périodique en θ.
Puisque la dynamique oscille à l'échelle de temps ε, il est nécessaire de choisir un pas de discrétisation ∆t de l'ordre de ε. Si c'est la dynamique limite ε → 0 que l'on cherche à approcher, cette technique devient très coûteuse, voire irréalisable. L'idée des méthodes de moyennisation est de remplacer l'étude de (0.1.12) par celle d'une équation autonome via un changement de variable. L'esprit d'une méthode de moyennisation est donné dans la proposition suivante, issue de [Cha13].

Proposition 0.1.11 Pour tout T ∞ > 0, il existe ε 0 > 0 tel que pour tout ε < ε 0 , il existe un changement de variables x 0 x(t)

Φε t = Id + O(ε) et une fonction F ε définie sur R n satisfaisant la relation ∀t ∈ [0, T ∞ ] , x (t) -Φε t ε • Ψε t (x 0 ) ≤ Cε r
x(t)

F 0 (x) = 1 T T 0 F 0 (x, θ) dθ, F 1 (x) = 1 T T 0 F 1 (x, θ) dθ - 1 2T T 0 θ 0 [F 0 (x, θ ′ ), F 0 (x, θ)] dθ ′ dθ où le crochet de Lie signifie [F 0 (x, θ ′ ), F 0 (x, θ)] := F ′ 0 (x, θ ′ )F 0 (x, θ) -F ′ 0 (x, θ)F 0 (x, θ ′ ).
Avec les notations de la Remarque 0.1.12, le système différentiel autonome associé à (0.1.12) est

ẋ = F 0 (x) , à l'ordre 0 et ẋ = F 0 (x) + ε F 1 (x), à l'ordre 1 en ε.
Dans la deuxième partie du chapitre 1, ces résultats classiques de moyennisation sont utilisés pour terminer l'étude du système différentiel (0.0.1). Ainsi, le caractère raide et le caractère hautement oscillant du système initial ont disparus. Les premiers ordres de l'équation différentielle autonome sont calculés explicitement.

Pour résumer, on peut ramener l'étude de (0.0.1), qui est un système couplé, à celle (0.1.10), puis résoudre ce dernier en utilisant des techniques de moyennisation. Mais cette approche n'est implémentable qu'à condition de savoir comment trouver x ε 0 à partir de x 0 , z 0 et ε. Cette question est donc fondamentale pour trouver des schémas numériques adaptés à l'étude de (0.0.1). La théorie des B-séries permet d'obtenir des développements formels à tout ordre de solutions d'équations différentielles et d'étudier les schémas qui les approchent. Nous l'avons adaptée pour obtenir un développement en série formelle de x ε 0 en fonction de ε, x 0 et z 0 . Mais avant d'expliquer ces prolongements, commenc ¸ons par appliquer les résultats précédents à un problème de dynamique des populations.

Application à un problème de dynamique des populations (Résumé du chapitre 2)

Avant d'expliquer ce qu'est la théorie des B-séries et les résultats auxquels elle nous a menés, il faut préciser le cheminement qu'a suivi cette thèse. En effet, le problème posé initialement n'était pas (0.0.1), mais un problème de dynamique des populations avec migrations.

UNE DYNAMIQUE DE VARI ÉT É CENTRALE xix

Le modèle considéré prend en compte à la fois les interactions entre les espèces et leurs migrations. Il s'agit donc d'une complexification de modèles écologiques plus stantards, dépendants uniquement du temps. Dans ce chapitre, on suit l'exemple des équations de Lotka-Volterra concernant l'interaction démographique proie-prédateur (mais n'importe quel autre modèle non linéaire d'interaction démograhique entre populations conviendrait). On suppose une différence d'échelle de temps entre les deux phénomènes : l'évolution démographique se déroule à l'échelle de la semaine ou du mois alors que les migrations spatiales ont lieu à l'échelle de l'heure ou de la journée. Le ratio ε entre ces deux échelles de temps est introduit. La question sur le comportement qualitatif est la suivante : comment des migrations spatiales rapides perturbent la dynamique lente de type Lotka-Volterra (et en particulier les cycles liés). On verra que cette séparation des échelles de temps fait que la répartition des populations tend, à une échelle de temps rapide, vers un "presque-équilibre" spatial, qui a son tour dépend du temps, mais cette fois à une échelle de temps lente.

L'espace est discrétisé en N sites, on introduit p ε le vecteur dont la coordonnée i représente le nombre de proies sur le site i et q ε le vecteur représentant les prédateurs.

Hypothèses 0.1.13 Les hypothèses concernant les migrations sont les suivantes :

-Les espèces peuvent se déplacer de n'importe quel site vers n'importe quel autre à tout instant.

-Les coefficients de migration sont supposés périodiques en la variable rapide t ε . C'est par exemple le cas du plancton, dont les déplacements dépendent de la luminosité et donc de l'heure.

-Le nombre d'individus est préservé lors des migrations (ils ne meurent pas pendant qu'ils migrent).

Quand la dépendance temporelle des opérateurs de migration est gelée, l'étude du système a été menée dans [START_REF] Poggiale | Lotka-Volterra's model and migrations: Breaking of the well-known center[END_REF] pour le cas de deux sites et dans [START_REF] Castella | A reduced model for spatially structured predator-prey systems with fast spatial migrations and slow demographic evolutions[END_REF] pour le cas continu en espace.

Hypothèses 0.1.14 Concernant l'interaction proie-prédateur, le point crucial est l'hétérogénéité spatiale induite par des coefficients qui diffèrent d'un site à l'autre. Ces différences témoignent par exemple de la présence de plus de nourriture sur un site, ou de plus d'endroits où se cacher pour les proies, deux situations qui entraînent une pression de prédation plus faible.

La dynamique et les hypothèses sont donc celles illustrées par la Figure 5. Le système est modélisé par

       dp ε (t) dt = 1 ε K p t ε p ε (t) + f (p ε (t), q ε (t)) , p ε (0) = p 0 dq ε (t) dt = 1 ε K q t ε q ε (t) + g (p ε (t), q ε (t)) , q ε (0) = q 0 , (0.1.14)
où p ε désigne les proies, q ε les prédateurs, f et g traduisent la dynamique d'interaction locale, ici une interaction de type Lotka-Volterra. Les migrations sont décrites par des opérateurs de Blotzmann linéaires K p et K q , avec

(K p (θ)) i,j = σ p i,j (θ) for i = j, (K p (θ)) i,i = - INTRODUCTION où σ p i,j (
θ) désigne le taux de proies migrant du site j vers le site i au temps θ et les définitions équivalentes pour K q . C'est en cherchant une variété centrale pour (0.1.14) que nous en sommes venus à considérer (0.0.1).

Dans le chapitre 2, les propriétés spectrales des opérateurs K p et K q sont établies. Elles donnent le changement de variable naturel pour séparer le terme raide (rapidement oscillant) linéaire (à spectre négatif) des termes non raides et non linéaires. En effectuant ce premier changement de variables

(p ε , q ε ) → (x ε , y ε ),
on obtient un système de la forme:

       dx ε (t) dt = F x ε (t), y ε (t), t ε dy ε (t) dt = 1 ε B t ε y ε (t) -1 ε ϕ x ε (t), t ε + G x ε (t), y ε (t), t ε . avec x ε ∈ R 2 et y ε ∈ R 2N -2 .
Il est presque du type (0.0.1), mais un terme indésirable -1 ε ϕ x ε (t), t ε est apparu. Pour régler ce problème, un deuxième changement de variables est effectué : une fonction h 0 x, t ε régulière et périodique en la variable de temps rapide est introduite. On pose

z ε (t) = y ε (t) -h 0 x ε (t), t ε , tel que (x ε (t), z ε (t)
) soit solution d'un système du type (0.0.1).

Les résultats du chapitre 1 sont alors appliqués à ce système. On montre l'existence d'une variété centrale telle que tout se ramène à l'étude de :

   ẋh = F x h , εh x h , t ε , t ε , x h (0) = x ε 0 ∈ R n z h (t) = εh u(t), t ε . (0.1.15)
Un premier avantage de cette méthode apparaît ici : on se ramène à l'étude de x h ∈ R2 , d'où une réduction des dimensions. Des développements explicites du système (0.1.15) sont calculés aux premiers ordres. De même, les premiers ordres de l'étape de moyennisation sont explicités. Tout cela donne des informations sur le comportement qualitatif de la solution de (0.1.14) : on montre que le système limite à l'ordre 0 en ε est de type Lotka-Volterra et que ses coefficients sont des moyennes en espace et en temps des coefficients originaux. L'étude des ordres supérieurs montre que les migrations peuvent entraîner une déstabilisation de ce cycle limite. De plus, on montre que la méthode "naive" consistant à moyenner les opérateurs K p et K q puis à appliquer la théorie classique de variété centrale peut mener à la prédiction d'un comportement qualitatif faux, contrairement à notre approche.

Pour finir, tous ces calculs sont menés sur un exemple de manière à illustrer les ordres de convergence obtenus théoriquement.

Une approche utilisant les séries formelles

Nous avons adapté l'étude de systèmes différentiels via les B-séries à l'étude d'une version simplifiée de notre système :

   ẋ = F (x, z) , x(0) = x 0 ∈ R n ż = -1 ε z + G (x, z) , z(0) = z 0 ∈ R m . (0.2.1)
Nous allons présenter dans un premier temps les bases de la théorie de B-séries 2 , avant d'en lister différentes extensions. Nous expliquerons ensuite comment nous l'avons adaptée à notre problème pour trouver des expressions explicites de x ε 0 et de la variété centrale.

La théorie des B-séries

Les arbres à une couleur J.C. Butcher a développé une théorie algébrique sur des objets appelés arbres (introduits dans [START_REF] Cayley | On the theory of the analytic forms called trees[END_REF], [START_REF] Merson | An operational method for the study of integration processes[END_REF]), particulièrement adaptée à l'étude des schémas numériques ([But72], [START_REF] Butcher | The numerical analysis of ordinary differential equations[END_REF]). Il s'agissait en effet de simplifier l'écriture des conditions de compatibilité des coefficients de schémas de Runge-Kutta à l'aide du groupe de Butcher. De ces premiers travaux, une théorie sur des séries formelles a été tirée dans [START_REF] Hairer | On the Butcher group and general multi-value methods[END_REF]. Ces séries sont appelées Butcher-series, ou encore B-series. Nous nous basons sur la présentation proposée dans [START_REF] Hairer | On the Butcher group and general multi-value methods[END_REF] pour la suite de cette partie.

Pour un petit paramètre ε, on considère le système différentiel

d t y = εf (y), y(t 0 ) = y 0 , (0.2.2) où d t = d dt .
Les premières dérivées de y s'écrivent:

d t y = εf (y), d 2 t y = ε 2 f ′ (y)d t y, d 3 t y = ε 3 [f ′′ (y) (d t y, d t y) + f ′ (y)d 2 t y], d 4 t y = ε 4 f (3) (y) (d t y, d t y, d t y) + 3f ′′ (y) d 2 t y, d t y + f ′ (y)d 3 t y , (0.2.3)
et en remplac ¸ant successivement y et ses dérivées par leurs expressions dans le membre de droite de (0.2.3), les dérivées de y s'écrivent uniquement à partir des dérivées de f :

d t y = εf (y), d 2 t y = ε 2 f ′ (y)f (y), d 3 t y = ε 3 (f ′′ (y) (f (y), f (y)) + f ′ (y)f ′ (y)f (y)), d 4 t y = ε 4 (f (3) (y) (f (y), f (y), f (y)) + 3f ′′ (y) (f ′ (y)f (y), f (y)) +f ′ (y)[f ′′ (y) (f (y), f (y)) + f ′ (y)f ′ (y)f (y)]).
Pour plus de clarté, on réécrit ces expressions sans expliciter la dépendance en y :

d t y = εf, d 2 t y = ε 2 f ′ f, d 3 t y = ε 3 [f ′′ f f + f ′ f ′ f ], d 4 t y = ε 4 [f (3) f f f + 3f ′′ f ′ f f + f ′ f ′′ f f + f ′ f ′ f ′ f ]. (0.2.4)
Pour automatiser ces développements, on associe un arbre à chacun de ces éléments différentiels.

Définition 0.2.1 L'ensemble T des arbres à racines est défini de manière récursive :

1. Le graphe

= u 2 = • ∈ T , l'élément noté [•, •] correspond à l'arbre , qui appartient donc à T . INTRODUCTION Définition 0.2.4 Soit f : R → R une fonction régulière. Pour tout u ∈ T , on note F u l'élément différentiel défini par : 1. F • (y) = f (y) 2. Si u = [u 1 , . . . , u n ] ∈ T , F u (y) = f (n) (y)[F u1 (y), . . . , F un (y)].
Les arbres permettent d'établir des résultats sur des méthodes numériques en comparant le développement de Taylor de la solution exacte de (0.2.2) au voisinage de t 0 et la solution d'un schéma numérique. En utilisant (0.2.4), le développement de Taylor de y(t) au voisinage de t 0 s'écrit :

y(t 0 + ∆t) = y(t 0 ) + ∆t d t y(t 0 ) + ∆t 2 2 d 2 t y(t 0 ) + ∆t 3 6 d 3 t y(t 0 ) + O(∆t 4 ) = y(t 0 ) + ε∆tf (y(t 0 )) + ε 2 ∆t 2 2 f ′ (y(t 0 ))f (y(t 0 )) + ε 3 ∆t 3 6 [f ′′ (y(t 0 ))(f (y(t 0 )), f (y(t 0 ))) + f ′ (y(t 0 ))f ′ (y(t 0 ))f (y(t 0 ))] + O(∆t 4 ) = y(t 0 ) + ε∆tF • (y(t 0 )) + ε 2 ∆t 2 2 F (y(t 0 )) + ε 3 ∆t 3 6 F (y(t 0 )) + F (y(t 0 )) + O(∆t 4 ), (0.2.5)
et l'étude de y(t) se ramène donc à l'étude de cette série. Les B-séries ont été introduites pour étudier de tels objets.

Pour les définir correctement, il nous manque encore une notion : le coefficient de symétrie. En effet, au cours des calculs, des termes du type f ′′ (x)(f (x), f (x)) apparaissent (représenté par l'arbre ). Ils présentent une symétrie, dont il faut tenir compte pour dénombrer correctement les arbres. Définition 0.2.5 Le coefficient de symétrie σ d'un arbre est défini par récurrence

1. σ • = 1, 2. Pour u = [u µ1
1 , . . . , u µp p ] ∈ T où les branches u i sont supposées distinctes et où µ i désigne le nombre d'occurences de la branche u i :

σ u = p i=1 µ i ! (σ ui ) µi .
Exemple 0.2.6 Pour calculer le coefficient de symétrie de l'arbre , on remarque que

= [•, •] = [• 2 ], puis on obtient σ = 2! (σ • ) 2 = 2.
Définition 0.2.7 Une B-série est une série formelle de la forme En revenant à (0.2.5), on remarque que les premiers termes du développement en B-séries de y(t) y ont été calculés. En effet, en prenant t 0 = 0, on a y(t) = B(a(t), y 0 ), (0.2.6) et par identification, les premiers coefficients valent :

B(a, y) = a ∅ y + u∈T ε |u| σ u a u F u (y),
a ∅ (t) = 1, a • (t) = t, a (t) = t 2 2 , a (t) = t 3 3 , a (t) = t 3 6 . ( 0
y(0) = B(a(0), y 0 ) y 0 = a ∅ (0)y 0 + u∈T ε |u| σ u a u (0)F u (y 0 ) d'où a ∅ (0) = 1 et ∀u ∈ T , a u (0) = 0.
Définition 0.2.9 La fonction 1 1 : T ∪ {∅} → R est définie par

1 1 ∅ = 1, ∀u ∈ T , 1 1 u = 0.
Ainsi, on a a(0) = 1 1. (0.2.8)

La composition des B-séries permet de calculer directement a u pour tout u ∈ T . L'ensemble des arbres muni de la loi de composition des B-séries forme un groupe, appelé Butcher-group. Ainsi, la loi de composition (3.3.3) s'écrit

(a * b) = b + 2b a • + b • (a • ) 2 + b ∅ a .
La loi de composition permet de trouver le développement en B-séries (0.2.6) de y(t) solution de (0.2.2). On commence par chercher le développement en B-série de εf (y), qui apparaît dans (0.2.2). Comme εf (y) = εF • (y), on a εf (y) = B(b, y) avec b • = 1 et pour tout arbre u ∈ T ∪ {∅}\{•}, b u = 0. Ainsi, en dérivant (0.2.6), on obtient successivement

d t y = εf (y) B(d t a(t), y 0 ) = εf (B(a(t), y 0 )) B(d t a(t), y 0 ) = B(b, B(a(t), y 0 )) B(d t a(t), y 0 ) = B(a(t) * b, y 0 ) d t a(t) = a(t
(t) = a ∅ (0) = 1, -ȧ• (t) = b • + b ∅ 0 a • (t) = 1 donc a • (t) = t, -ȧ (t) = b + b • a • (t) + b ∅ a (t) = t donc a (t) = t 2 2 , -ȧ (t) = b + 2b a • (t) + b • a • (t) 2 + b ∅ a (t) = t 2 donc a (t) = t 3 3 , ȧ (t) = b + b a • (t) + b • a (t) + b ∅ a (t) = t 2 2 donc a (t) =
-Partition à 1 élément : { }, -Partitions à 2 éléments : {•, } et { , •}, -Partition à 3 éléments : {•, •, •}.
Ainsi, par (0.2.11), on obtient

a -1 = -a + 2a a • -a 3 • .
Les B-séries permettent d'expliciter les solutions de systèmes bien plus compliqués que (0.2.2) et de les comparer avec leur approximation par un schéma numérique. Nous donnons ici quelques exemples de prolongements.

Les arbres à deux couleurs

Systèmes partitionnés: Les B-séries ont été adaptées par E. Hairer dans [START_REF] Hairer | Order conditions for numerical methods for partitioned ordinary differential equation[END_REF] en une théorie des P-séries, conc ¸ue pour des systèmes partitionnés du type :

ẋ = f 1 (x, y) ẏ = f 2 (x, y)
.

xxv On suit ici la présentation de [START_REF] Chartier | Algebraic structure of B-series[END_REF]. On définit deux sortes de noeuds : soient T • l'ensemble des arbres à racine en • = f 1 , T • l'ensemble des arbres à racine en • = f 2 et ∅ • , ∅ • les arbres vides correspondants. Les éléments différentiels sont définis pour u 1 , . . . , u p ∈ T • et v 1 , . . . , v q ∈ T • par F [u1,...,up,v1,...,vq]• (x, y) = ∂ p x ∂ q y f 1 (x, y)[F u1 (x, y), . . . , F vq (x, y)], F [u1,...,up,v1,...,vq]• (x, y) = ∂ p x ∂ q y f 2 (x, y)[F u1 (x, y), . . . , F vq (x, y)].

Les séries associées à ces arbres à deux couleurs, appelées P-séries, sont définies par

P (a, (x, y)) =    a ∅• x + u∈T• ε |u| σu a u F u (x, y) a ∅• y + v∈T• ε |v| σv a v F v (x, y)    .
En reprenant la définition des partitions d'un arbre présentée dans le cas d'arbres à une couleur pour la transposer aux arbres à deux couleurs, on montre que la règle de composition est inchangée.

Systèmes splittés: La théorie des B-séries a été adaptée dans [AMSS] à l'étude de systèmes différentiels du type:

ẏ = f 1 (t, y) + f 2 (t, y), y(0) = y 0 ∈ R d .
On définit à nouveau deux sortes de noeuds : • pour f 1 et • pour f 2 . On obtient ainsi des arbres de la forme

• correspond à F • (y) = f 1 (y), • correspond à F • (y) = f 2 (y), correspond à F (y) = f ′ 2 (y)f 1 (y).
(0.2.12) Soit T l'ensemble des arbres à deux couleurs. Les éléments différentiels associés sont

F [u1,...,un]• (y) = f (n) 1 (y)[F u1 (y), . . . , F un (y)], F [u1,...,un]• (y) = f (n) 2 (y)[F u1 (y), . . . , F un (y)].
À nouveau, la loi de composition est inchangée.

On peut définir des B-séries pour ces arbres à deux couleurs :

B(a, y) = a ∅ y + u∈T ε |u| σ u a u F u (y).
Elles sont par exemple utilisées pour démontrer la formule de Baker-Campbell-Haussdorf, liée aux méthodes numériques de splitting.

Les arbres indicés

Systèmes rapidement oscillants : L'un des objets de cette thèse est l'adaptation de la théorie des B-séries à l'étude des techniques d'approximation d'ordre élevé pour le système (0.2.1). Ce travail s'inspire des résultats obtenus dans [START_REF] Chartier | Higher order averaging, formal series and numerical integration i: B-series[END_REF] et [START_REF] Chartier | Higher order averaging, formal series and numerical integration ii: the quasi-periodic case[END_REF] pour l'étude de méthodes de moyennisation sur des systèmes rapidement oscillants du type

∀t ∈ [0, L/ε], ẏ = εf (y, tω), y(0) = y 0 ∈ R d . ( 0
(y) = f (n) k (y)[F u1 (y), . . . , F un (y)],
et la loi de composition reste inchangée. L'équation (0.2.14) peut être réécrite en termes de B-séries : On montre alors que γ vérifie l'équation de transport

εf (y, θ) = ε k∈Z e ikθ F • k (y) = k∈Z εβ • k (θ)F • k (y) = B(β(θ), y), avec ∀k ∈ Z, β • k (θ) =
∀t ≥ 0, ∀θ ∈ [0, 2π] d , ∂ t γ(t, θ) + ω.∇ θ γ(t, θ) = γ(t, θ) * β(θ), γ(0, 0) = 1 1. (0.2.17)
En utilisant l'unicité de la solution de (0.2.17), une loi de groupe sur γ est démontrée

∀t, t ′ ≥ 0, γ(t ′ , 0) * γ(t, 0) = γ(t + t ′ , 0), ∀t ≥ 0, ∀θ ∈ [0, 2π] d , γ(t, 0) * γ(0, θ) = γ(t, θ). ( 0
µ(t) = γ(t, 0), λ(θ) = γ(0, θ).
On connaît donc explicitement les développements en B-séries (à tout ordre) de Φε θ et Ψε t . Ces développements permettent d'étudier les propriétés des méthodes d'approximation numérique utilisant la moyennisation.

Remarque 0.2.17 Le cas de fonctions f presque-périodiques est traité dans [START_REF] Chartier | Higher order averaging, formal series and numerical integration ii: the quasi-periodic case[END_REF].

Remarque 0.2.18 Une étude d'erreur peut être réalisée sur ces B-séries, voir [START_REF] Chartier | A formal series approach to averaging: exponentially small error estimates[END_REF] pour le cas périodique et [START_REF] Chartier | Higher order averaging, formal series and numerical integration iii: error bounds[END_REF] pour le cas presque-périodique.

Les arbres indicés à deux couleurs (Résumé du chapitre 3)

La dernière partie de cette thèse, présentée dans le chapitre 3, a consisté à adapter toutes ces extensions des B-séries à l'étude de (0.2.1) :

   ẋ = F (x, z) , x(0) = x 0 ∈ R n ż = -1 ε z + G (x, z) , z(0) = z 0 ∈ R m .
On ne travaille pas directement avec z(t), mais plutôt avec y(t) = e t ε z(t). On remarque que y(0) = y 0 = z 0 . Le système différentiel (0.2.1) devient alors :

       ẋ = F x, e -t ε y , x(0) = x 0 ∈ R n ẏ = e t ε G x, e -t ε y , y(0) = y 0 ∈ R m . (0.2.22)
Si on renormalise (0.2.22) en posant t ← t ε , on obtient finalement le système différentiel :

   ẋ = ε F t (x, y) , x(0) = x 0 ∈ R n ẏ = ε G t (x, y) , y(0) = y 0 ∈ R m , (0.2.23) où F t (x, y) = F (x, e -t y) et G t (x, y) = e t G (
x, e -t y). On conserve cette normalisation pour toute la suite du chapitre.

Des arbres à deux couleurs et à indice sont introduits. En ce qui concerne les couleurs, il s'agit comme dans le cas des systèmes partitionnés d'associer les noeuds de type

• à F et ceux de type • à G.
Pour ce qui est des indices, ils apparaissent dans le cas des méthodes de moyennisation du fait d'un développement en série de Fourier et l'indice k fait référence à l'exponentielle e ikθ . Ici, c'est un développement de Taylor autour de z = 0 joue ce rôle et l'indice k fait référence à e -kt . Nous voulons approcher la dynamique limite de la solution de (0.2.1). Nous avons vu dans l'introduction du chapitre 1 que pour toute condition initiale

z(t) = εh (x(t), t) + O e -µt .
Ainsi, obtenir un développement de z(t) en puissances de e -t nous donne une information sur la variété centrale εh. Nous ne donnons pas plus de détails dans l'introduction sur la définition des éléments différentiels associés aux arbres. La démarche est présentée au début du chapitre 3.

Pour étudier (0.2.23), on définit α : R × T → R telle que :

x(t) y(t) = B(α(t), (x 0 , y 0 )) =    α ∅• (t) x 0 + u∈T• ε |u| σu α u (t)F u (x 0 , y 0 ) α ∅• (t) y 0 + v∈T• ε |v| σv α v (t)F v (x 0 , y 0 )    .

INTRODUCTION

On montre alors que α vérifie une équation différentielle ordinaire du type (0.2.16) avec α(0) = 1 1. On en déduit l'existence de γ(t, τ ) fonction polynomiale4 en (t, e -τ ) telle que ∀t ≥ 0, α(t) = γ (t, t) .

(0.2.24)

Alors γ est solution d'une équation de transport de type (0.2.17) avec condition initiale γ(0, 0) = 1 1. Nous nous intéressons au comportement asymptotique de la solution et donc à une condition initiale en (0, +∞), dans laquelle la solution a déjà convergé vers la variété centrale. Cependant, il n'y pas unicité des solutions polynomiales de cette équation de transport avec condition initiale en (0, +∞). Le défaut d'unicité vient du fait que la condition initiale est prise en (0, +∞) pour tous les arbres. En effet, la même équation de transport avec condition initiale en (0, +∞) sur T • et en (0, 0) sur T • admet une unique solution polynomiale. On note δ l'unique solution polynomiale de cette équation de transport avec condition initiale :

∀u ∈ T • , δ u (0, +∞) = 1 1 u , ∀v ∈ T • , δ v (0, 0) = 1 1 v .
Pour comprendre le lien qui unit γ et δ, on définit γ : R × T → R comme :

γ| T• (t) = γ| T• (t, +∞) et γ| T• (t) = γ| T• (t, 0). (0.2.25)
On note Φ 0 l'application (x, y) → Φ 0 (x, y) := B(γ -1 (0), (x, y)). Alors l'application inverse vérifie

Φ -1 0 : (x 0 , y 0 ) → (x ε 0 , y 0 ).
Le passage de x 0 à x ε 0 est donc explicite, ce qui est le premier pas pour approcher numériquement la solution de (0.2.1) par une dynamique sur une variété centrale.

Remarque 0.2.19 Cette question avait déjà été étudiée dans des articles comme [START_REF] Roberts | Appropriate initial conditions for asymptotic descriptions of the long term evolution of dynamical systems[END_REF], qui donne une méthode de calcul pour les premiers ordres sur différents exemples.

On montre alors que B(δ(t, τ ), (x ε 0 , y 0 )) correspond à B(γ(t, τ ), (x 0 , y 0 )) :

δ(t, τ ) = γ-1 (0) * γ(t, τ ),
et δ représente donc (x(t), y(t)) solution de (0.2.23) à un changement de condition initiale près.

Ensuite, on utilise l'unicité dans l'équation de transport pour démontrer une loi de groupe sur δ :

∀t, t ′ ∈ R, ∀τ ∈ R, δ(t + t ′ , τ ) = δ(t ′ ) * δ(t, τ ),
où δ est défini à partir de δ comme γ à partir de γ dans (0.2.25).

Nous étudions le comportement asymptotique de la solution, une fois que les e -kτ ont convergé vers 0 et la solution vers la variété centrale. Pour c ¸a, on définit :

δ ∞ (t) = lim τ →+∞ δ(t, τ ) 5 .
On montre alors que la variété centrale (x, εh(x)) est décrite par Π(x) = B(δ ∞ (0), x). Sur la variété, la dynamique est régie par :

ẋ∞ (t) = εF x ∞ (t), εh(x ∞ (t)) z ∞ (t) = εh(x ∞ (t))
, (0.2.26) et toutes ces fonctions possèdent des développements en B-séries explicites. Cette étude est résumée dans la Figure 6.

xxix

Φ -1 0 Π (x 0 , z 0 ) -→ (x ε 0 , z 0 ) -→ (x ε 0 , εh(x ε 0 ))
Flot de (0.2.1)

    Flot de (0.2.26) (x(t), z(t)) (x ∞ (t), z ∞ (t))
Figure 6: Variété centrale et équation exacte.

La dynamique sur la variété centrale ne correspond à celle de la solution que pour t grand, puisque l'erreur est en O (e -µt ). Comme pour les méthodes de moyennisation, on aimerait approcher la solution de (0.2.1) pour tout temps.

On définit alors δ(t) = γ-1 (0) * γ(t), qui vérifie la loi de groupe :

δ(t + t ′ ) = δ(t ′ ) * δ(t).
L'équivalent de l'équation (0.2.18) est alors :

δ(t, τ ) = δ(t) * δ(0, τ ). ( 0
(x(t), z(t)) = Φ t (x ∞ (t), z(t)).
Les développements explicites en B-séries de toutes ces fonctions sont connus. Ces résultats sont représentés sur la Figure 7.

Φ -1 0 (x 0 , z 0 ) -→ (x ε 0 , z 0 ) Flot   de (0.2.1) Flot   de (0.2.28) (x(t), z(t)) ←- (x ∞ (t), z(t)) Φ t

INTRODUCTION

Ensuite, nous expliquons comment implémenter le calcul de coefficients de B-séries en Maple, c'est-à-dire comment calculer une composée ou un inverse de B-séries, de manière à évaluer systématiquement les coefficients qui apparaissent dans l'étude de (0.2.1) menée ci-dessus.

Enfin, ces résultats sont illustrés par deux exemples : l'un avec n = m = 1, l'autre avec n = 2 et m = 1. Les B-séries sont tronquées à l'ordre 3. Le changement de condition initiale Φ -1 0 est calculé, puis (x ∞ (t), z ∞ (t)) la solution sur la variété centrale est implémentée et comparée à (x(t), z(t)). L'erreur est de la forme O e -µt + ε 4 . Ensuite, la forme normale (x ∞ (t), z(t)) est calculée, ainsi que le changement de variable Φ t et on montre que Φ t (x ∞ (t), z(t)) est une approximation en O(ε 4 ) de (x(t), z(t)).

A fast time dependent center manifold theorem: existence and approximation results

We study a differential system of the form:

   ẋ = F x, z, t ε , x(0) = x 0 ∈ R n ż = 1 ε B t ε z + G x, z, t ε , z(0) = z 0 ∈ R m , (1.1.1)
with F and G smooth functions1 , with a smooth dependence in ε. The transport operator B is related to an exponentially decreasing resolvent. The functions are T -periodic in θ = t ε , and we denote T = [0, T ]. In subsection 1.1.1 we prove (Theorem 1.1.2) that, associated with any system of the form (1.1.1), there is a fast time dependent center manifold. This means that for ε small enough there exists a function

h ε = h ε (x, θ), periodic in θ = t ε , such that if z(0) = h ε (x(0), 0) (i.
e. z(0) belongs to the center manifold initially) then for all subsequent times t we have z(t) = h ε x(t), t ε (i.e. z(t) belongs to the t ε -dependent center manifold for all times). The θ-dependent set {(x, z) ∈ R n × R m such that z = h ε (x, θ)} is called the fast time dependent central manifold, owing to the fact that it is a center manifold which varies with the fast time θ = t ε . In subsection 1.1.2 (Theorem 1.1.4) we prove that the central manifold in fact attracts all solutions to the above system, with an exponentially small convergence rate of size exp(-c t ε ) for some c > 0. In other words, even if the initial data do not satisfy the relation z(0) = h ε (x(0), 0), the solution (x, z) will nevertheless be exponentially close, to within exp(-c t ε ), to a solution (x h , z h ) to the reduced system

   ẋh = F x h , h ε x h , t ε , t ε , x h (0) = x ε 0 z h = h ε x h , t ε ,
with x ε 0 an altered initial data. We have thus removed the stiffness and reduced dimensionality of the problem, the only remaining unknown being x h .

Remark 1.1.1 The obtention of an explicit expansion for x ε 0 will be one of the motivations of the work of Chap. 3. In the first two chapters, we do not deal with this issue.

In subsection 1.1.3 we prove (Theorem 1.1.6) that h ε can be approximated as

h ε (x, θ) = h 0 (x, θ) + ε h 1 (x, θ) + • • • + ε r h r (x, θ) + O(ε r+1 ) := h [r] ε (x, θ) + O(ε r+1 )
to within any order r+1, and the formulae can be made explicit. Naturally, the solution

(x, z) is O exp(-c t ε ) + ε r+1 - close to a solution (x h [r] , z h [r] ) to the reduced system      ẋh [r] = F x h [r] , h [r] ε x h [r] , t ε , t ε , x h [r] (0) = x ε 0 z h [r] = h [r] ε x h [r] , t ε .
Subsection 1.1.4 is devoted to deriving the explicit values of the first few terms h 0 , h 1 , h 2 in the above expansion.

Existence of a fast time dependent center manifold

Theorem 1.1.2 (Existence of a fast time dependent center manifold) Consider a differential system of the form

     ẋ = F (x, z, θ) , x(0) = x 0 θ = 1 ε , θ(0) = θ 0 ż = 1 ε B(θ)z + G (x, z, θ) , z(0) = z 0 , (1.1.2)
where the functions F and G are defined and C 1 on R n × R m × T, the operator B(θ) on R m is C 1 and periodic, and its resolvent R(t, s) satisfies

∀t ≥ s, R(t, s) ≤ C 0 e -µ0(t-s) ,
for some constants C 0 ≥ 1 and µ 0 > 0. In addition, assume that F and G are bounded, with bounded first derivatives in (x, z, θ). Last, assume that F and G belong to C 0 (T; C r (R n × R m ) for some given integer r ≥ 1, and that F , G have bounded derivatives w.r.t. (x, z) up to order r.

Then, there exists ε 0 > 0, and a function

h ε (x, θ) ∈ C 1 (R n × T) ∩ C 0 (T; C r (R n ))
, defined for all 0 < ε < ε 0 , with the following property. For all x 0 ∈ R n and θ 0 ∈ T, the solution (x(t), θ(t), z(t)) of (1.1.1) with initial conditions

x(0) = x 0 , θ(0) = θ 0 , z(0) = h ε (x 0 , θ 0 ), satisfies the relation, for all t, z(t) = h ε x(t), θ 0 + t ε .
Proof:[of Theorem 1.1.2] Our proof closely follows that of [START_REF] Carr | Applications of Centre Manifold Theory[END_REF] in the case where the fast time dependence is not present. It also is inspired by [START_REF] Sakamoto | Invariant manifolds in singular perturbation problems for ordinary differential equations[END_REF]. Throughout this proof we denote by M and L real numbers such that

F (x, z, θ) ≤ M, G(x, z, θ) ≤ M, ∂ (x,z) F (x, z, θ) ≤ L, ∂ (x,z) G(x, z, θ) ≤ L,
for any (x, z, θ), where ∂ (x,z) denotes the Jacobian matrix with respect to (x, z).

First step. Center manifold as the fixed-point of an operator T . Fix a smooth function

(x, θ) ∈ R n × T → h(x, θ).
Take initial values (x 0 , z 0 , θ 0 ) ∈ R n × R m × T and denote by θ ≡ θ(t) ≡ θ(t, θ 0 ), x h ≡ x h (t) ≡ x h (t, x 0 , θ 0 ), and z h ≡ z h (t) ≡ z h (t, x 0 , θ 0 ), the solution components of the differential system

       ẋh = F (x h , h(x h , θ), θ), x h (0, x 0 , θ 0 ) = x 0 , θ = 1 ε , θ(0, θ 0 ) = θ 0 , żh = 1 ε B(θ)z h + G(x h , h(x h , θ), θ), z h (0, x 0 , θ 0 ) = z 0 .
(1.1.3)

Given that x h can be obtained independently of z h , z h can in turn be obtained as follows. Using R(t, s) the resolvent of the operator B(θ), and defining for brevity

G h (s, x 0 , θ 0 ) := G(x h (s, x 0 , θ 0 ), h(x h (s, x 0 , θ 0 ), θ 0 + s/ε), θ 0 + s/ε) (1.1.4)
the Duhamel formula provides

z h (t, x 0 , θ 0 ) = R θ 0 + t ε , θ 0 z 0 + ε t ε 0 R(θ 0 , θ 0 + u) G h (εu, x 0 , θ 0 ) du . (1.1.5)
Now, if we look for a center manifold, then the function h is such that z h (t, x 0 , θ 0 ) coincides with h(x h (t, x 0 , θ 0 ), θ 0 + t ε ) for all values of t, x 0 and θ 0 . Therefore, the function z h (t, x 0 , θ 0 ), whenever h is the seeked center manifold, should be in particular bounded for all times provided x h (t, x 0 , θ 0 ) is, and provided h(x, θ) is smooth in x and periodic w.r.t. θ. As a consequence, the initial datum z 0 in (1.1.5) can not be chosen freely but should rather be an initial value that makes z h (t, x 0 , θ 0 ) bounded for all times. The boundedness when t → +∞ is automatically ensured (thanks to the exponential decay of the resolvent R), yet the boundedness as t → -∞ selects a unique initial datum, namely2 ,

z 0 = ε 0 -∞ R(θ 0 , θ 0 + u) G h (εu, x 0 , θ 0 ) du (1.1.6)
and accordingly we recover

z h (t, x 0 , θ 0 ) = ε t ε -∞ R θ 0 + t ε , θ 0 + u G h (εu, x 0 , θ 0 ) du. (1.1.7)
All in all we finally define T h as the function which maps

(x 0 , θ 0 ) ∈ R n × T to (T h)(x 0 , θ 0 ) = ε 0 -∞ R(θ 0 , θ 0 + u) G h (εu, x 0 , θ 0 ) du, (1.1.8)
where G h is given by formula (1.1.4), and T is now defined as an operator on functions h. We have shown up to now that if h defines a center manifold, then h necessarily is a fixed point of T .

Let us show conversely that if h is a fixed point of T , then the initial relation

z h (0, x 0 , θ 0 ) = h(x 0 , θ 0 ) implies z h (t, x 0 , θ 0 ) = h(x h (t, x 0 , θ 0 ), θ 0 + t ε )
for all subsequent times t. To this aim, we fix x 0 and θ 0 and use the definition of T to rewrite the identity h = T h, namely

h x h (t, x 0 , θ 0 ), θ 0 + t ε = ε 0 -∞ R θ 0 + t ε , θ 0 + t ε + u G h εu, x h (t, x 0 , θ 0 ), θ 0 + t ε du.
Owing to the group law

∀(t, t ′ ), x h t ′ , x h (t, x 0 , θ 0 ), θ 0 + t ε = x h (t + t ′ , x 0 , θ 0 ),
we have, using the definition (1.1.4) of G h , the relation

G h εu, x h (t, x 0 , θ 0 ), θ 0 + t ε = G h (εu + t, x 0 , θ 0 ), which leads to h x h (t, x 0 , θ 0 ), θ 0 + t ε = ε 0 -∞ R θ 0 + t ε , θ 0 + t ε + u G h (εu + t, x 0 , θ 0 ) du = ε t ε -∞ R θ 0 + t ε , θ 0 + u G h (εu, x 0 , θ 0 ) du = z h (t, x 0 , θ 0 ),
and the last equality follows from (1.1.7).

Second step. The operator T maps F to F . Define F as the functional space

F = h ∈ C 1 (R n × T; R m ), such that h ∞ ≤ 1 and ∂ x h ∞ ≤ 1 where ∂ x h ∞ = ∂ x h L ∞ (R n ×T,L(R n ,R m ))
. We want to show that T maps F to itself. Given h ∈ F and the definition (1.1.8) of T h, we have for all

(x 0 , θ 0 ) ∈ R n × T T h(x 0 , θ 0 ) ≤ ε 0 -∞ R(θ 0 , θ 0 + u) G h (εu, x 0 , θ 0 ) du ≤ ε 0 -∞ C 0 e µ0u G(x h (εu, x 0 , θ 0 ), h(x h (εu, x 0 , θ 0 ), θ 0 + u), θ 0 + u) du ≤ ε 0 -∞ C 0 e µ0u M du = C 0 ε M µ 0 ≤ 1, provided ε ≤ µ 0 /(C 0 M ).
On the other hand, since h, F , and G are periodic w.r.t. θ, it is clear that x h and G h are periodic w.r.t.

θ 0 . Since R(θ 0 , θ 0 + u) is periodic w.r.t. θ 0 , it comes that T h is periodic w.r.t. θ 0 . Similarly, the C 1 smoothness of T h, provided h ∈ C 1 , is clear. It remains to prove that ∂ x (T h) ∞ ≤ 1.
To this aim, we first estimate the Jacobian of the flow ∂ x0 x h (t, x 0 , θ 0 ). From the variational equation we have

(∂ x0 x h )(t, x 0 , θ 0 ) = Id + t 0 ∂ x F (x h , h(x h , θ 0 + s/ε), θ 0 + s/ε)) • ∂ x0 x h + ∂ z F (x h , h(x h , θ 0 + s/ε), θ 0 + t/ε) • (∂ x h)(x h , θ 0 + s/ε) • ∂ x0 x h ds.
Therefore, it follows

(∂ x0 x h ) ≤ 1 + 2L [0,t] ∂ x0 x h ds, (here [0, t] denotes either [0, t] if t ≥ 0 or [t, 0] if t ≤ 0) which, owing to Gronwall's Lemma, leads to ∀t ∈ R, ∂ x0 x h (t, ., .) ≤ e 2L|t| .
(1.1.9)

Substituting this estimate into the equation obtained by differentiating T h we get

∂ x0 T h(x 0 , θ 0 ) ≤ C 0 ε 0 -∞ e µ0u ∂ x G • ∂ x0 x h + ∂ z G • ∂ x h • ∂ x0 x h du,
where all arguments of h, x h , F and G are as in (1.1.8) and have been omitted for the sake of clarity. Using (1.1.9) it then follows that

∂ x0 T h(x 0 , θ 0 ) ≤ 2 C 0 εL 0 -∞ e µ0u e 2εL|u| du = 2 C 0 ε L µ 0 -2 ε L ≤ 1, provided µ 0 -2 ε L > 0, that is ε < µ 0 /(2L), and provided ε ≤ µ 0 /[(2 C 0 + 1) L].
Third step 3. Operator T is a contraction on F . Consider h 1 and h 2 two functions belonging to F . The corresponding functions x h1 and x h2 satisfy the estimate

(x h1 -x h2 )(t, x 0 , θ 0 ) ≤ L [0,t] x h1 -x h2 + h 1 (x h1 , θ 0 + u) -h 2 (x h2 , θ 0 + u) du
where, once again, the arguments (u, x 0 , θ 0 ) of x h1 and x h2 on the r.h.s have been omitted for brevity. It is straightforward to write, say with θ = θ 0 + u, that

h 1 (x h1 , θ) -h 2 (x h2 , θ) ≤ h 1 (x h1 , θ) -h 1 (x h2 , θ) + h 1 (x h2 , θ) -h 2 (x h2 , θ) ≤ x h1 -x h2 + h 1 -h 2 ∞ .
Hence, we recover

(x h1 -x h2 )(t, x 0 , θ 0 ) ≤ L |t| h 1 -h 2 ∞ + 2 L [0,t] x h1 -x h2 du,
and by Gronwall's Lemma we obtain

x h1 (t, x 0 , θ 0 ) -x h2 (t, x 0 , θ 0 ) ∞ ≤ L |t| e 2L|t| h 1 -h 2 ∞ .
Consequently, we have

T h 1 (x 0 , θ 0 ) -T h 2 (x 0 , θ 0 ) ≤ C 0 ε L 0 -∞ e µ0u x h1 (εu, x 0 , θ 0 ) -x h2 (εu, x 0 , θ 0 ) + h 1 (x h1 (εu, x 0 , θ 0 ), θ 0 + u) -h 2 (x h2 (εu, x 0 , θ 0 ), θ 0 + u) du ≤ C 0 ε L h 1 -h 2 ∞ 0 -∞ e µ0u 2 L ε |u|e 2εL|u| + 1 du ≤ C 0 ε L 1 µ 0 + 2 L ε (µ 0 -2εL) 2 h 1 -h 2 ∞ ,
(1.1.10) so that T : F → F is a contraction for small enough values of ε.

Hence T possesses a unique fixed point h ε = h ε (x, θ) ∈ C 1 , as asserted in Theorem 1.1.2.

Fourth step. Smoothness of h ε .

The idea is to repeat the proof of the third step. We introduce the set

F r = {h ∈ C 0 (T; C r (R n )), such that h ∞ ≤ 1 and ∂ k x h ∞ ≤ 1 for all k = 1, . . . , r}
Since all derivatives up to order r of F and G are bounded, inequality (1.1.10) is simply replaced by

∂ k x (T h 1 (x 0 , θ 0 ) -T h 2 (x 0 , θ 0 )) ≤ C ε h 1 -h 2 ∞ + ∂ x (h 1 -h 2 ) ∞ + . . . + ∂ k x (h 1 -h 2 ) ∞ ,
where C is a constant that does ont depend on ε, and the norm used is the induced norm on k-linear functions. By choosing ε small enough, we again obtain a contraction map. The smoothness of the fixed point h ε in x thus follows.

Remark 1.1.3 The function h ε also depends smoothly on ε, as it is obtained as the limit of the convergent iteration

h ε = lim k→∞ T k h 0 from a ε-independent h 0 . It is thus C ∞ w.r.t. ε.

Exponential convergence of all solutions towards the fast time dependent center manifold

In this section we prove that when the initial conditions do not fulfill the requirement z(0) = h ε (x(0), 0), then the exact solution of the differential system (1.1.2) nevertheless goes exponentially fast to the center manifold. The precise statement is the following.

Theorem 1.1.4 (Error relative to the center manifold) Consider the system (1.1.2) in Theorem 1.1.2, and assume all assumptions in this Theorem are fulfilled. For any 0 < ε ≤ ε 0 , where ε 0 is as in Theorem 1.1.2, denote by x(t) and z(t) the solutions of system (1.1.2) with prescribed initial values (x 0 , z 0 , θ 0 ). Fix some µ > 0 such that µ < µ 0 . Then, there is an ε 1 satisfying 0 < ε 1 ≤ ε 0 , there exists a constant C > 0 (independent of ε and T ∞ below), such that for any 0 < ε ≤ ε 1 , the following holds.

1. Exponential convergence towards the center manifold.

∀t ≥ 0, z(t) -h ε x(t), θ 0 + t ε ≤ C e -µt/ε .

Shadowing principle.

Pick up an arbitrary T ∞ > 0. There exists an altered initial data x ε 0 (implicitely depending on T ∞ ), such that the solution components of the reduced system

             dx h dt = F x h , h ε x h , θ 0 + t ε , θ 0 + t ε x h (0) = x ε 0 z h (t) = h ε x h (t), θ 0 + t ε satisfy the following error estimate ∀t ∈ [0, T ∞ ], z(t) -z h (t) + x(t) -x h (t) ≤ Ce -µt/ε .
Moreover, if the original solution x(t) to (1.1.2) turns out to be bounded on R + , we can take T ∞ = +∞ in the above estimate.

Remark 1.1.5 The stiff, fast-oscillating, n + m dimensional system (1.1.2), with unknowns (x(t), z(t)), coincides up to exponentially small terms with the above reduced system with unknown x h , a non-stiff, n dimensional system, where the fast oscillations are kept (eventually they will be taken care of through standard averaging techniques).

Proof:[of Theorem 1.1.4] Proof of the first statement.

By construction of the function h ε , we know that any solution (

x(t), z(t), θ(t)) to (1.1.2) with z(0) = h ε ( x(0), 0) satisfies z(t) = h ε x(t), θ(0) + t
ε for all times. Hence, differentiating this last expression w.r.t. time immediately gives the relation

∂ t h ε x(t), θ(0) + t ε = 1 ε B θ(0) + t ε z(t) + G x(t), z(t), θ(0) + t ε
and, inserting the differential equation satisfied by x(t), next letting the initial data (x 0 , θ 0 ) take any value in R n ×T, eventually provides that for any x ∈ R n and any θ ∈ T, we have

1 ε ∂ θ h ε (x, θ) + ∂ x h ε (x, θ)F (x, h ε (x, θ), θ) = 1 ε B(θ)h ε (x, θ) + G(x, h ε (x, θ), θ). (1.1.11)
Coming back to the selected solution (x(t), z(t)) to (1.1.2) that is referred to in Theorem 1.1.4, we compute the time derivative of the function h ε (x(t), θ 0 + t ε ). Using the above relation immediately gives

d dt h ε x(t), θ 0 + t ε = 1 ε B θ 0 + t ε h ε x(t), θ 0 + t ε + G x(t), h ε x(t), θ 0 + t ε , θ 0 + t ε .
The Duhamel formula then leads to

z(t) -h ε x(t), θ 0 + t ε = R θ 0 + t ε , θ 0 (z 0 -h ε (x 0 , θ 0 )) + ε t ε 0 R θ 0 + t ε , θ 0 + u ∆G(εu) du
where we define

(∆G)(εu) := G(x(εu), z(εu), θ 0 + u) -G(x(εu), h ε (x(εu), θ 0 + u) , θ 0 + u).
Using the fact that G is globaly Lipschitz in z (with Lipschitz constant L), together with the exponential decay of the resolvent R (with constants C 0 and µ 0 ), we then obtain the following inequality

z(t) -h ε x(t), θ 0 + t ε ≤ C 0 e -µ0 t ε z 0 -h ε (x 0 , θ 0 )) + C 0 ε t ε 0 e -µ0( t ε -u) (∆G)(εu) du ≤ C 0 e -µ0 t ε z 0 -h ε (x 0 , θ 0 )) + C 0 L ε t ε 0 e -µ0( t ε -u) z(εu) -h ε (x(εu), θ 0 + u) du.
The Gronwall Lemma applied to the function

e +µ0 t ε z(t) -h ε (x(t), θ 0 + t ε ) provides z(t) -h x(t), t ε ≤ C 0 exp -(µ 0 -C 0 ε L) t ε z 0 -h ε (x 0 , θ 0 ) .
The first statement of the Theorem follows, upon taking 0 < µ < µ 0 -εCL with 0 < ε < µ/(C 0 L).

Proof of the second statement. We fix T ∞ > 0. We denote

x ∞ := x(T ∞ ),
where (x(t), z(t)) is, as before, the selected solution to (1.1.2) that is referred to in Theorem 1.1.4. The function x → F (x, h ε (x, θ), θ) being globally Lipschitz w.r.t. x, the backwards system, with initial datum at t = T ∞ ,

ẋh = F x h , h ε x h , θ 0 + t ε , θ 0 + t ε , x h (T ∞ ) = x ∞ , (1.1.12)
possesses a unique (global) solution x h (t). We set

x ε 0 := x h (0), with x h (t) as above. In this perspective, x h (0) is thus deduced from a scattering solution, whose initial datum is defined from its value at the final time T ∞ . With these definitions at hand, proving the second part of the Theorem reduces to proving that x(t)x h (t) is exponentially small. To do so, we rewrite the differential equation satisfied by x(t) as

ẋ(t) = F x(t), h ε x(t), θ 0 + t ε , θ 0 + t ε + δ(t)
where

δ(t) := F x(t), z(t), θ 0 + t ε -F x(t), h ε x(t), θ 0 + t ε , θ 0 + t ε .
According to the -already proved -first statement of Theorem 1.1.4, we have

δ(t) ≤ L z(t) -h ε (x(t), θ 0 + t ε ) ≤ C L e -µt/ε ,
for some C > 0, where 0 < µ < µ 0 is fixed. Thus, writing

ẋ -ẋh = F x, h ε x, θ 0 + t ε , θ 0 + t ε -F x h , h ε x h , θ 0 + t ε , θ 0 + t ε + δ(t), using the initial datum x(T ∞ ) = x h (T ∞ ) at t = T ∞ , and integrating backwards from T ∞ to t whenever 0 ≤ t ≤ T ∞ , provides x(t) -x h (t) ≤ L T∞ t (1 + ∂ x h ε ∞ ) x(u) -x h (u) + C L e -µu/ε du ≤ C 1 T∞ t x(u) -x h (u) du + C 2 ε e -µt/ε ,
for some constants C 1 , C 2 , independent of ε and T ∞ . The Gronwall Lemma provides

x(t) -x h (t) ≤ C 2 ε e -µt/ε + C 1 T∞ t e C1(t-s) C 2 ε e -µs/ε ds ≤ C 2 ε e -µt/ε + C 1 ε µ + C 1 ε e C1t e -(µ+C1 ε) t ε .
It follows, taking some 0 < ν < µ, that there is a constant C > 0, independent of ε and T ∞ , such that, for any t ∈ [0, T ∞ ], we have as desired

x(t) -x h (t) ≤ C e -ν t ε .
Now, in the special case where the original solution x(t) to (1.1.2) is bounded on [0, +∞[, we can somewhat pursue the above analysis, as follows. Taking a sequence of values of T ∞ that go to +∞, we may assume that the sequence of values x(T ∞ ) converges towards some x ∞ . This produces a sequence of functions x h,T∞ (t) as above, whose dependence on T ∞ is stressed explicitely, and we have for each T ∞ the relation

∀t ∈ [0, T ∞ ], x(t) -x h,T∞ (t) ≤ C e -ν t ε ,
where 0 < ν < µ 0 is fixed but arbitrary. This produces a sequence of values x h,T∞ (0) that is automatically bounded, thanks to the above estimate. We may therefore assume that x h,T∞ (0) goes to some x ε 0 as T ∞ → ∞. We now set x h (t) as the solution to the differential system

ẋh = F x h , h ε x h , θ 0 + t ε , θ 0 + t ε , x h (0) = x ε 0 . (1.1.13)
It is clear that x h,T∞ (t) goes to x h (t) on compact sets. We deduce (upon fixing the value of t, e.g.) that

∀t ≥ 0, x(t) -x h (t) ≤ C e -ν t ε ,
where the constants C and ν are unchanged. The result follows in the case T ∞ = +∞.

Approximation of the center manifold

In this section, we show that the function h ε can be expanded in powers of ε up to every order k ≤ r, where each coefficient-function can be computed recursively.

Theorem 1.1.6 (Approximation of h ε ) Under the assumptions and notations of Theorem 1.1.2 and 1.1.4, the following statements hold true.

1. The function h ε satisfies the following partial differential equation for all x ∈ R n and all θ ∈ T

1 ε ∂ θ h ε (x, θ) -B(θ)h ε (x, θ) = G(x, h ε (x, θ), θ) -∂ x h ε (x, θ)F (x, h ε (x, θ), θ) . (1.1.14)

The terms of the formal expansion

h ε = εh 1 + ε 2 h 2 + • • • of h ε ,
where the h n 's are independent of ε, are defined in a unique way by an equation of the form

∀n ∈ N, (∂ θ -B(θ)) h n+1 = J n (x, θ)
where J n depends only on derivatives of F and G up to order n, whenever n ≤ r. Furthermore, the function

h [r] ε := εh 1 + ε 2 h 2 + • • • + ε r h r
satisfies equation (1.1.14) up to an error term of size ε r and one has the following estimate for some positive constant

C r h ε -h [r] ε ∞ ≤ C r ε r+1 . (1.1.15) Proof:[of Theorem 1.1.6]
The validity of PDE (1.1.14) has been already established in the course of the proof of Theorem 1.1.4, see (1.1.11). Hence only the second statement is to be proved.

We now look for an expansion of h ε in powers of ε of the form

h ε (x, θ) = h 0 (x, θ) + εh 1 (x, θ) + • • • + ε n h n (x, θ) + • • •
To obtain candidate functions h n 's, we insert the previous expression into equation (1.1.14) and equate like powers of ε. At order ε -1 , this gives the homologic equation

B(θ)h 0 (x, θ) = ∂ θ h 0 (x, θ).
To solve this equation, and forthcoming analogous equations satisfied by the h n 's, we argue as follows. This is an homogeneous linear differential equation in θ, whose solution can be expressed as h 0 (•, θ) = R(θ, 0)h 0 (•, 0). The initial condition h 0 (•, 0) is a priori not prescribed, however h 0 is assumed periodic in θ. Therefore we necessarily have (Id -R(T, 0))h 0 (., 0) = 0. The estimate R(θ, 0) ≤ C 0 e -µ0θ provides the invertibility of the operator (Id -R(T, 0)). We deduce

h 0 ≡ 0. (1.1.16)
We next proceed to derive the equation satisfied by h 1 . It reads

B(θ)h 1 (x, θ) + G(x, 0, θ) = ∂ θ h 1 (x, θ), which implies h 1 (•, θ) = R(θ, 0)h 1 (•, 0) + θ 0 R(θ, ϕ) G(x, 0, ϕ) dϕ.
For h 1 to be periodic with period T , we need

(Id -R(T, 0))h 1 (•, 0) = T 0 R(T, ϕ) G(x, 0, ϕ) dϕ.
This selects the only possible choice for h 1 (.0). After easy computations we recover

h 1 (•, θ) = R(θ, 0)(Id -R(T, 0)) -1 θ θ-T R(0, ϕ) G(x, 0, ϕ) dϕ.
(1.1.17)

In other words, introducing the nonlocal operator

S : v = v(θ) → S(v)(θ) = R(θ, 0)(Id -R(T, 0)) -1 θ θ-T R(0, ϕ) v(ϕ) dϕ, (1.1.18)
we have

h 1 (x, θ) = S (G(x, 0, •)) (θ).
We may informally write

S = (∂ θ -B(θ)) -1 ,
in that for each periodic function v, operator S selects the unique periodic solution w to

(∂ θ -B(θ)) w = v.
More generally, assuming h 1 , . . ., h n are known, the function h n+1 satisfies an equation of the form

B(θ)h n+1 (x, θ) -∂ θ h n+1 (x, θ) = J n (x, θ)
where J n contains various derivatives of F and G up to order n w.r.t. x only, and is periodic w.r.t. θ. The same arguments as above allow to conclude that

h n+1 (•, θ) = R(θ, 0)(Id -R(T, 0)) -1 T θ-T R(0, ϕ) J n (x, ϕ) dϕ,
or, in other words,

h n+1 (•, θ) = S (J n (x, •)) (θ),
which, given the assumptions on F and G, provides a function h n+1 which is bounded and has bounded derivatives w.r.t. x up to order rn.

Consider now the truncated expansion h

[r] ε = εh 1 + . . .+ ε r h r of h ε , where r is fixed. The function h

[r]

ε satisfies the partial differential equation (1.1.14)

1 ε B(θ)h [r] ε (x, θ) + G(x, h [r] ε (x, θ), θ) = ∂ x h [r] ε (x, θ)F x, h [r] ε (x, θ), θ + 1 ε ∂ θ h [r] ε (x, θ) + δ ε (x, θ), (1.1.19)
up to a defect δ ε (x, θ) which is a continuous function from R n × T into R m thanks to the previous observations, and which is by construction bounded by C ε r for some constant C > 0 independent of ε, in the norm of continuous and bounded functions of (x, θ). Let us now introduce x h (t) ≡ x h (t, x 0 , θ 0 ) as the solution on the central manifold, i.e. the solution of equation (1.1.3), and consider

z h [r] (t) := h [r] ε (x h (t), θ 0 + t ε ), seen as an approximation of z h (t) := h ε (x h (t), θ 0 + t ε ) defined in (1.1.
3). We have, using (1.1.19) and the chain rule,

dz h [r] (t) dt = 1 ε B θ 0 + t ε z h [r] (t) + G (x h (t), z h [r] (t)) -δ ε x h (t), θ 0 + t ε , (1.1.20)
where we have omitted the arguments (t, x 0 , θ 0 ) of x h for brevity. Proceeding as in Theorem 1.1.2 (both h and h are bounded by construction), i.e. using the Duhamel formula in (1.1.20) to express z h [r] (t) as a function of z h [r] (0), and expressing the necessary boundedness of z h [r] (t) as t → -∞, eventually provides, as in (1.1.6), the relation

z h [r] (0) = ε 0 -∞ R(θ 0 , θ 0 + u) [G (x h (εu), z h [r] (εu)) -δ ε (x h (εu), θ 0 + u)] du.
Remembering the analogous relation

z h (0) = ε 0 -∞ R(θ 0 , θ 0 + u) [G (x h (εu), z h (εu))] du, using z h [r] (0) = h [r]
ε (x 0 , θ 0 ) and z h (0) = h ε (x 0 , θ 0 ), and taking differences, eventually provides

∆h ε (x 0 , θ 0 ) := h ε (x 0 , θ 0 ) -h [r] ε (x 0 , θ 0 ) ≤ C 0 ε 0 -∞ e µ0s [L z h (εu) -z h [r] (εu) + C ε r ] du ≤ C C 0 µ -1 0 ε r+1 + C 0 L ε 0 -∞ e µ0s h ε (x h (εu)) -h [r] ε (x h (εu)) du ≤ C C 0 µ -1 0 ε r+1 + C 0 L µ -1 0 ε ∆h ε ∞ . Therefore, provided C 0 L µ -1 0 ε < 1 we recover h ε (x 0 , θ 0 ) -h [r] ε (x 0 , θ 0 ) ≤ C C 0 µ -1 0 1 -C 0 L µ -1 0 ε ε r+1 .
Estimate (1.1.15) of the Theorem follows. Associated with the previous Theorem, the following convergence result is natural.

Theorem 1.1.7 (Shadowing principle for the truncated system) Under the asumptions and notation of Theorem 1.1.4, pick up T ∞ > 0 and take the associated modified initial data x ε 0 . Let (x(t), z(t)) be the solution of (1.1.2). Define x h [r] as the solution of following differential system, called "truncated system" in the sequel,

   dx h [r] dt = F x h [r] , h [r] ε x h [r] , θ 0 + t ε , θ 0 + t ε , x h [r] (0) = x ε 0 ,
and define

z h [r] (t) = h [r] ε x h [r] (t), t ε , where h [r] ε = εh 1 + • • • + ε r h r is as in Theorem 1.1.6. Then, we have the following estimate ∀t ∈ [0, T ∞ ], z(t) -z h [r] (t) + x(t) -x h [r] (t) ≤ C ε r+1 + e -µt/ε ,
where C > 0 and 0 < µ < µ 0 are constants independent of t and ε. Moreover, if the solution x is bounded on R + , we can take T ∞ = +∞ in the above estimate.

Proof: The result follows directly from Theorem 1.1.4 and from estimate (1.1.15).

Derivation of the first few terms of the expansion in the general case

In this subsection, we derive the explicit expressions of the first two terms h 1 and h 2 of the expansion of

h ε = εh 1 + ε 2 h 2 + • • • .
We have already shown in the course of Theorem 1.1.6 that h 0 ≡ 0 and that h 1 is given, see (1.1.17), by the equation

h 1 (x, θ) = R(θ, 0)(Id -R(T, 0)) -1 θ θ-T R(0, ϕ)G(x, 0, ϕ) dϕ,
or, in other words,

h 1 (x, θ) = S (G(x, 0, •)) (θ),
where the nonlocal operator S = (∂ θ -B(θ))

-1 is defined in (1.1.18). Now, the equation at order 1 in ε gives

∂ θ h 2 = B(θ)h 2 + ∂ z G(x, 0, θ) • h 1 -∂ x h 1 • F (x, 0, θ) =: B(θ)h 2 + J 1 (x, θ),
where we use the notation of the proof of Theorem 1.1.6, and we define

J 1 (x, θ) = ∂ z G(x, 0, θ) • h 1 (x, θ) -∂ x h 1 (x, θ) • F (x, 0, θ). (1.1.21)
Eventually we arrive at

h 2 (x, θ) = R(θ, 0)(Id -R(T, 0)) -1 θ θ-T R(0, ϕ)J 1 (x, ϕ) dϕ,
or, in other words,

h 2 (x, θ) = S (J 1 (x, •)) (θ),
where J 1 is given by (1.1.21), or equivalently by

J 1 (x, θ) = ∂ z G(x, 0, θ) • S (G(x, 0, •)) (θ) -S (∂ x G(x, 0, •)) (θ) • F (x, 0, θ).
Associated with these values of h 0 , h 1 , and h 2 , we are now in position to write the few first reduced systems associated with the original (1.1.2), i.e. associated with

ẋ = F (x, z, t ε ), ż = 1 ε B( t ε )z + G(x, z, t ε ).
They read as follows.

The zeroth order reduced system is

   ẋ(t) = F x(t), 0, t ε , with z(t) = 0.
It provides a dynamics which approximates that of (1.1.2) to within O(ε + exp(-µ t ε )). The first order reduced system is

                   ẋ(t) = F       x(t), =ε h 1 (x(t), t ε ) ε S (G(x(t), 0, •)) t ε , t ε       , with z(t) = ε S (G(x(t), 0, •)) t ε = ε h 1 (x(t), t ε ) , with S = (∂ θ -B(θ)) -1 , see (1.1.18
). It provides a dynamics which approximates that of (1.1.2) to within

O(ε 2 + exp(-µ t ε )). The second order reduced system is        ẋ(t) = F x(t), ε h 1 x(t), t ε + ε 2 h 2 x(t), t ε , t ε , with z(t) = ε h 1 x(t), t ε + ε 2 h 2 x(t), t ε ,
where

h 1 (x, θ) = S (G(x, 0, •)) (θ) and h 2 (x, θ) = ∂ z G(x, 0, θ)•S (G(x, 0, •)) (θ)-S (∂ x G(x, 0, •)) (θ)•F (x, 0, θ), and S = (∂ θ -B(θ)) -1 , see (1.1.18
). It provides a dynamics which approximates that of (1.1.2) to within O(ε 3 + exp(-µ t ε )).

Averaging out the remaining oscillations

Up to now we have successively reduced our original model (1.1.1)

dx dt = F x, z, t ε , dz dt = 1 ε B t ε + G x, z, t ε ,
to a system of the reduced form

dx h dt = F x h , h ε x h , t ε , t ε .
We are thus left with a differential system with periodic, fastly oscillating coefficients. This is amenable to a (standard) averaging procedure. Note in passing that, since h ε can be expanded as

h 0 + εh 1 + • • •+ ε r h r + O(ε r+1 ),
we may as well expand the right-hand-side of the above equation

F (x, h ε (x, t ε ), t ε ) = F 0 (x, t ε ) + εF 1 (x, t ε ) + • • • + ε r-1 F r-1 (x, t ε ) + O(ε r
), with the assumed C r smoothness of the coefficients. In this section we prove that our system can be approximated by an autonomous system of the form

dx av dt = F ε (x av ),
where

F ε = F 0 + ε F 1 + ε 2 F 2 + • • • can be computed explicitely.
Though this is an easy result now, it nevertheless is the final result of our analysis: we have eventually removed both stiffness and fast oscillations in our original model (1.1.1). Besides, we provide explicit formulae for the first few terms of this expansion.

The averaging theorem

The following averaging Theorem is standard.

Theorem 1.2.1 (Averaging Theorem -borrowed from [Cha13]) Consider a highly oscillatory differential system of the form of the form

dx dt = F ε x, t ε , x(0) = x 0 , (1.2.1)
where F ε (x, θ) is smooth in the sense that for each ε we have

F ε (x, θ) ∈ C 0 (θ ∈ T; C r (x ∈ R n ))
, for some given r ≥ 1, with globally bounded derivatives up to order r. Lastly, assume F ε can be expanded in powers of ε as

F ε (x, θ) = F 0 (x, θ) + εF 1 (x, θ) + . . . + ε r-1 F r-1 (x, θ) + O(ε r ),
where the remainder term holds in the topology of C 0 (T;

C 1 (R n )).
Then, the following holds. For all T ∞ > 0, there exists ε 0 > 0 such that for all ε < ε 0 , there exists a change of variables

Φ ε θ = Id + O(ε) ∈ C 0 (T; C r (R n ))
, and a function F ε (x) = F 0 (x) + ε F 1 (x) + . . ., belonging to C r (R n ), with bounded derivatives up to order r, such that the following holds. Introduce x av (t) as the solution to

dx av dt = F ε (x av ), x av (0) = x 0 .
We have

∀t ∈ [0, T ∞ ] , x (t) -Φ ε t ε (x av (t)) ≤ C ε r , (1.2.2)
for some C > 0 independent of ε. The expansions

Φ ε θ = Id + ε Φ 1 θ + • • • and F ε = F 0 + ε F 1 + • • • can be explicitely performed up to order O(ε r ) in the topology of C 0 (T; C 1 (R n )) resp. C 1 (R n ).
Remark 1.2.2 There are variants of this statement which lead to exponentially small errors terms of the size O(exp(-c/ε)) for some c > 0. We do not dwell on that aspect later. We simply refer to [Cha13] for the corresponding averaging statements.

Remark 1.2.3 The following formulae are proved in [Cha13]

F 0 (x) = 1 T T 0 F 0 (x, θ) dθ, F 1 (x) = 1 T T 0 F 1 (x, θ) dθ - 1 2T T 0 θ 0 [F 0 (x, θ ′ ), F 0 (x, θ)] dθ ′ dθ
where the Lie-bracket stands for

[F 0 (x, θ ′ ), F 0 (x, θ)] := ∂ x F 0 (x, θ ′ )F 0 (x, θ) -∂ x F 0 (x, θ)F 0 (x, θ ′ ).
With this notation, the approximate differential system naturally associated with

dx dt = F ε x, t ε , is thus dx av dt = F 0 (x av )
, to zero order, and dx av dt = F 0 (x av ) + ε F 1 (x av ), to first order in ε.

Application to (1.1.1)

Combining Theorems 1.2.1, 1.1.7, 1.1.6, 1.1.4, and 1.1.2 eventually gives the following.

Theorem 1.2.4 (Averaging out the fast time dependent central manfold) Under the assumptions and notation of Theorem 1.1.2, consider the system

     ẋ = F (x, z, θ) , x(0) = x 0 θ = 1 ε , θ(0) = θ 0 ż = 1 ε B(θ)z + G (x, z, θ) , z(0) = z 0 , Pick up 0 < ε ≤ ε 0 , 0 < µ < µ 0 , 0 < T ∞ < +∞,
and define x ε 0 as in Theorem 1.1.4. Then, there is an ε 1 satisfying 0 < ε 1 ≤ ε 0 , there exists a constant C > 0, independent of ε. such that for any 0 < ε ≤ ε 1 , the following holds.

There exists a change of variables

Φ ε θ = Id + O(ε) on R n , and a function F ε (x) = F 0 (x) + ε F 1 (x) + . . . + ε r-1 F r-1 (x)
, such that the solution of the reduced system

dx av dt = F ε (x av ), x av (0) = x ε 0 satisfies the following error estimate ∀t ∈ [0, T ∞ ], x(t) -Φ ε t ε (x av (t)) ≤ C e -µt/ε + ε r .
As a consequence of all the previously derived formulae we have, with the above notation

F 0 (x) = 1 T T 0 F (x, 0, θ) dθ, (1.2.3) F 1 (x) = 1 T T 0 ∂ z F (x, 0, θ) [S (G(x, 0, .))] (θ) dθ - 1 2T T 0 θ 0 [F (x, 0, θ ′ ), F (x, 0, θ)] dθ ′ dθ (1.2.4)
where we recall that S = (∂ θ -B(θ)) -1 and the Lie-bracket is

[F (x, 0, θ ′ ), F (x, 0, θ)] := ∂ x F (x, 0, θ ′ ) F (x, 0, θ) -∂ x F (x, 0, θ) F (x, 0, θ ′ ).
Here we have used the relation h ε (x, θ) = εh 1 (x, θ) + O(ε 2 ) with the identity h 1 (x, θ) = S (G(x, 0, .)) (θ), and we have written the first ordrer expansion

F (x, h ε (x), θ) = F (x, 0, θ) + ε∂ z F (x, 0, θ) h 1 (x, θ) + O(ε 2 ).

Introduction

This chapter is concerned with the analytical study of the dynamics of a prey-predator model with spatial migrations and is a follow-up of [START_REF] Castella | A reduced model for spatially structured predator-prey systems with fast spatial migrations and slow demographic evolutions[END_REF], an article which itself had been motivated by the work [START_REF] Poggiale | Lotka-Volterra's model and migrations: Breaking of the well-known center[END_REF]. The model under consideration takes into account both interactions between species and spatial migrations, and as such it is a strict elaboration of more standard, purely time dependent ecological models, such as the Lotka-Volterra equations, which we take as a paradigm throughout the present paper (but any nonlinear model for the basic demography on each site would do). The key point is, we superimpose spatial migrations to the mere temporal demographic evolution. In particular, a fundamental feature of the operating dynamics that we wish to mention right away is the occurrence of two time-scales, accounting for the fact that spatial evolutions are vastly faster than demographical ones. The qualitative question is: how do the fast spatial migrations perturb the overall underlying slow Lotka-Volterra dynamics (in particular the underlying cycles).

In this still simplified version, the space is discretized into N distinct sites amongst which species move rapidly (i.e. change from one site to another within, say, a few hours). These migrations are described by two linear operators -of the linear Boltzmann type -corresponding to migrations of preys, on the one hand, and of predators, on the other hand. We also assume, and this is a key point here, that the coefficients involved in the migrations operators depend periodically on time in a highly-oscillatory way. One may think of preys and predators migrating on the time scale of an hour, say, with migrations rates which vary on the same time scale, due to dayly variations of the environment. This is typically the case for plancton, whose motion in the vertical direction depends on the light brought by the sun during the day. In addition, and lastly, these operators are assumed to preserve the numbers of individuals, so that migration and demographic terms remain independent in the equations (individuals can not die while migrating). Note that the case when the time dependence of these coefficients is frozen is dealt with in [START_REF] Poggiale | Lotka-Volterra's model and migrations: Breaking of the well-known center[END_REF] in the case of 2 sites, as well as [START_REF] Castella | A reduced model for spatially structured predator-prey systems with fast spatial migrations and slow demographic evolutions[END_REF] in the case of continuously many sites.

Concerning the predator-prey interactions, we choose to describe the basic demographic evolution on each site through a term of Lotka-Volterra type. The point is, the coefficients involved in the Lotka-Volterra system are assumed to differ from one site to another, that is to say, spatial characteristics may vary (for instance owing to more abundant food or more spots to hide implying a lower predation pressure, or so). Naturally, the overall effect of demography only becomes apparent over intervals of time that can be gauged in weeks or months, say, rather than the previously mentionned time-scale of several hours. For this reason, we introduce at once the small parameter ε defined as the ratio between the two time-scales (time-scale of migrations vs. time-scale of demography).

The complete model shall be presented with full details in section 2.2. However, in this introductory section, it is enlightening to describe it in an abstract and concise form as follows

d dt (populations)(t) = 1 ε (migration term) t ε • (populations)(t)+(predator -prey interactions)(t), (2.1.1)
where the migration term is periodic in t/ε. The aim of this work is to conduct an analysis of the dynamics in the limit ε → 0 and to draw conclusions from the resulting asymptotic model.

To do so, the qualititative difficulty is the following. The migration term tends to bring the population to an "almost-equilibrium" that is reached within the fast time-scale t/ε, corresponding to vanishing migratory fluxes, yet the underlying equilibrium in turn depends on the fast time-scale in a periodic way, and keeps on oscillating while never being actually reached.

Our results are as follows. Mathematically speaking, and to deal with this specific situation, we first reduce the original equations through several changes of variables to a form which is amenable to center manifold techniques. This apparently harmless step turns out to be crucial in order to identify the dominant term in a kind of center manifold approach (see section 2.3). More precisely, the point here is to separate a stiff, fast-oscillating, linear term, with negative spectrum, and various non-stiff, nonlinear terms (possibly fast oscillating as well), which act at higher order from a central manifold perspective. Once this is done, and in a second step (see section 2.4), we use the Center Manifold Theorem proved in chapter 1 (see section 1.1) to deduce that the overall dynamics is roughly dictated by a reduced, fast time-dependent and periodic differential system, governing the mere evolution of the total number of preys and predators (irrespective of the sites).

As a third step, we use the previous results to study the qualitative behaviour of the solution (see section 2.5). The first terms of the expansion in ε are explicitely computed. The averaging techniques of section 1.2 are applied to this example. The limiting dynamics is then computed, as well as its correctors at first orders. Eventually, a qualitative analysis of the zero and first order averaged systems is performed. In the last section, the full methodology is worked through for an example implying two sites (see section 2.6) . The limit equations are shown to be of Lotka-Volterra type, with coefficients that are appropriate averages in space and time of the original Lotka-Volterra coefficients and of the transfer operators. The lower order dynamics therefore possesses cycles, and we show using the first order corrector that the cycles may be naturally destabilized in the original fast-oscillating system. This eventually shows that fast spatial migrations may destabilize the underlying Lotka-Volterra dynamics. Lastly, we show that the "naive" methods consisting in taking the original system (2.1.1), first averaging out (in time) the transfer operators, then performing a standard central manifold technique to obtain a limiting dynamics, may give the wrong prediction in terms of the qualitative behaviour of the whole system as ε goes to zero. In other words, the time-dependent manifold technique we develop here, followed by an averaging procedure in time, are the necessary steps in order to recover the correct qualitative prediction in terms of ecological aspects (e.g. extinction of survival of the species at hand).

Description of the model

We consider prey and predators evolving in a spatial domain that is the union of N ∈ N * subdomains, usually called "patches" in the ecological literature. We pick up a small dimensionless parameter ε expressing the ratio between the time-scales of migrations and of demographic evolution, and denote by

p ε i (t), i ∈ [[1, N ]],
the number of preys occupying the i th site at time t, and by

q ε i (t), i ∈ [[1, N ]
] the corresponding number of predators. Introducing the vectors

p ε (t) = (p ε 1 (t), • • • , p ε N (t)) T ∈ R N and q ε (t) = (q ε 1 (t), • • • , q ε N (t)) T ∈ R N ,
the initial-value problem can be written as

         dp ε (t) dt = 1 ε K p t ε p ε (t) + f (p ε (t), q ε (t)) , p ε (0) = p 0 , dq ε (t) dt = 1 ε K q t ε q ε (t) + g (p ε (t), q ε (t)) , q ε (0) = q 0 , (2.2.1)
where K p and K q are time-dependent transport matrices defined by

(K p (θ)) i,j = σ p i,j (θ) for i = j, (K p (θ)) i,i = - N k=1 σ p k,i (θ),
and by similar equations for K q . Here, the transfer rates σ p i,j (θ) and σ q i,j (θ) are the proportions of preys, respectively predators, moving from site j to site i at time θ. These functions are assumed to satisfy the following key assumption

∀i = j, σ p i,j (θ) is T -periodic with T = 2π, (2.2.2) ∀i = j, ∀θ, σ p i,j (θ) > 0,
(2.2.3) while σ p i,j is assumed smooth in θ. The similar assumptions are made for σ q i,j (θ). Note that, as a direct consequence of these assumptions, the vector 1 = (1, . . . , 1) T ∈ R N is a left-eigenvector of both K p and K q , namely 1 T K p (θ) = 1 T K q (θ) = 0. This translates the conservation of mass along spatial migrations: individuals do not die along their spatial migrations. Besides, the periodicity assumption (2.2.2) corresponds to diurnal variations of the migration coefficients, while the positivity assumption (2.2.3) ensures that populations actually migrate from any site i to any site j at any time, a fact that evenually ensures existence of a (unique, time-dependent) equilibrium repartition on the sites for each time θ, i.e. a unique right-eigenvector satisfying K p (θ) p eq (θ) = 0 and similarly with "q", for each θ, as proved below -see Section 2.3.1.

As for the functions f and g, they model non-linear interactions of Lotka-Volterra type between species, namely

f = (f 1 , . . . , f N ) T , g = (g 1 , . . . , g N ) T ,
where, for each site index 1 ≤ i ≤ N , we have

f i (p ε (t), q ε (t)) = a p,i p ε i (t) -b p,i p ε i (t) q ε i (t), g i (p ε (t), q ε (t)) = -a q,i q ε i (t) + b q,i p ε i (t) q ε i (t).
(2.2.4)

Here, coefficient a p,i is the birth-rate of preys on site i, while a q,i is the death-rate of predators on site i, and in the same way, b p,i is the death-rate of preys on site i caused by the predators, while b q,i is the birth-rate of predators due to the presence of preys. All those quantities are assumed non-negative. We also assume for simplicity that these quantities do not depend on time, a harmless assumption. We will sometimes write (2.2.4) in the shorthand form

f (p ε (t), q ε (t)) = a p p ε (t) -b p p ε (t) q ε (t), g (p ε (t), q ε (t)) = -a q q ε (t) + b q p ε (t) q ε (t),
where the notation a p p ε , b p p ε q ε , and so on, implicitely denote the componentwise vector products

a p p ε := (a p,i p ε i ) N i=1 ∈ R N , b p p ε q ε := (b p,i p ε i q ε i ) N i=1 ∈ R N
, and so on.

(2.2.5)

The rest of this chapter is devoted to the analysis of the limit ε → 0 in (2.2.1), under the above assumptions. Note in passing that the mere existence of solutions to (2.2.1) is obvious for each fixed value of ε, in the sense that for any ε > 0 there is a solution (p ε (t), q ε (t)) defined on some time interval [0, T ε ] with T ε > 0. The fact that there is a common T > 0 such that solutions to (2.2.1) exist on [0, T ] uniformly in ε is not clear at once. This fact comes below as a byproduct of our analysis -see Corollary 2.3.6.

Remark 2.2.1 For the sake of simplicity, interactions between species are modeled here with constant coefficients for the species interactions, namely a p,i does not depend on the slow time t nor on the fast periodic time t/ε, etc. However, our results would remain true mutatis mutandis with time-dependent coefficients, that is in the case where a p.i ≡ a p,i (t, t/ε), where a p.i (t, θ) is smooth on R × T, and similarly for the other birth/death rates.

In the similar spirit, note that we retain here and throughout this chapter the Lotka-Volterra interaction terms (2.2.4) as a paradigm which we fully analyze. Needless to say, any other basic nonlinear predator-prey model, based on fixed smooth nonlinear interaction-terms f (p, q) and g(p, q), whose underlying ecological coefficients are then assumed to both depend on the site i as well as on the slow time t and the fast time t/ε in a smooth and periodic fashion, could be analyzed along the similar techniques. We do not dwell on that aspect later.

Analysis and reduction of the system

In this section we transform the original equations (2.2.1) into a form that is amenable to the center manifold approach presented in Section 1.1, separating a fast variable for which relaxation to an almost-equilibrium occurs, and a slow variable which eventually carries the global demographic evolution.

Main properties of the linear part of the system

In this subsection, we present the spectral properties of the transport operators K p and K q . They translate the above mentionned relaxation towards an almost-equilibrium. Our main result in this paragraph is Proposition 2.3.2. It is obtained using the generalized entropy method [START_REF] Perthame | Transport equations in biology[END_REF].

Proposition 2.3.1 (Spectral decomposition associated with K p and K q ) For all θ ∈ T, the operator K p (θ) admits 0 as a simple eigenvalue, while all other eigenvalues of K p (θ) have a negative real part. The right-eigenvector p eq (θ) associated to 0 can be chosen as a smooth function w.r.t. θ, with positive entries, satisfying 1 T p eq = 1.

Moreover the following property holds true, namely

∃β > 0, ∀θ ∈ T, ∀λ ∈ Sp(K p (θ)) \ {0}, ℜ(λ) ≤ -2β.

Besides, denoting

E 0 = {z ∈ R N , 1 T z = 0} = {z ∈ R N , N i=1 z i = 0}, we have ∀θ ∈ T, E 0 ⊕ Span(p eq (θ)) = R N ,
and E 0 is stable upon the action of K p (θ).

The same properties hold for matrix K q , where the right-eigenvector associated with the simple eigenvalue 0 is denoted by q eq (θ), and E 0 is again a stable supplementary space.

Proposition 2.3.2 (Exponential decay of the resolvent of K p and K q ) Following the previous Proposition, define for each θ K p (θ) : E 0 → E 0 as the linear operator induced by K p (θ) on the stable subspace E 0 . Define R p (t, s) as the resolvent associated with the variable coefficients matrix K p (θ), namely as the solution to

d dt R p (t, s) = K p (t)R p (t, s), R p (s, s) = Id, (2.3.1)
Then, there exists µ p > 0 and C p ≥ 1 such that for any t ≥ s we have

R p (t, s) ≤ C p e -µp(t-s) . (2.3.2)
The similar statements hold true for matrix K q (θ) and its resolvent R q (t, s).

Remark 2.3.3 Note that Proposition 2.3.2 provides a reinforced version of Proposition 2.3.1, in that Proposition 2.3.1 asserts the fact that K p (θ) has a spectrum lying in {ℜ(z) ≤ -2β < 0} for any θ, while Proposition 2.3.2 asserts that the resolvent of K p (θ) has exponential decay. Note also that Proposition 2.3.1 by no means implies Proposition 2.3.2, despite the fact that there is a fixed supplementary space E 0 , independent of θ, which carries all negative eigenvalues of K p (θ). The oscillations of the equilibrium eigenspace Span(p eq (θ)) as θ varies may destroy the expected exponential decay of the resolvent R p (t, s) in general. Only the specific structure of the migration operator K p (θ) allows to recover Proposition 2.3.2, thanks to the generalized entropy method [START_REF] Perthame | Transport equations in biology[END_REF],

Proof:[of Proposition 2.3.1] The proof is classical. Fix θ ∈ T and denote K p ≡ K p (θ).

Relation 1 T K p = 0 ∈ R N (conservation of mass) proves that 1 is an eigenvector of K T p associated with the eigenvalue 0. On the other hand, taking λ > 0 large enough, matrix λId + K T p possesses only positive entries (due to the positivity of the σ p i,j 's). Besides, using again relation 1 T K p = 0 ∈ R N we deduce that for any i, the sum of the entries of K T p on line i vanishes, so that the (positive) entries of λId + K T p on line i sum up to λ. Hence the spectral radius ρ(λId + K T p ) is ≤ λ, and the previous observation provides ρ(λId + K T p ) = λ. The Perron-Frobenius Theorem then asserts that ρ(λId + K T p ) = λ is a simple eigenvalue of λId + K T p . As a consequence, removing the term λId, matrix K T p admits 0 as a simple eigenvalue, and all other eigenvalues of K T p have real part ≤ -2β < 0 for some β > 0. Removing the transposition provides that 0 is a simple eigenvalue of K p , and all other eigenvalues have real part ≤ -2β < 0 as well. Lastly, applying again the Perron-Frobenius Theorem, we observe that λId + K p only possesses positive entries, and we deduce that the eigenspace of K p associated with the eigenvalue 0 is generated by some vector p eq which only possesses positive entries, and which we may normalize so as to ensure 1 T p eq = 0.

This proves the first part of the above statement for each fixed value of θ.

The smooth dependence of p eq upon θ then stems form the smooth dependence of K p upon θ, together with the fact that the eigenvalue 0 has constant multiplicity 1 for any θ. The uniformity of the upper bound β above also comes from the continuity of K p upon θ ∈ T.

Lastly, the spectral decomposition

E 0 ⊕ Span(p eq (θ)) = R N is obvious, since relation 1 T K p = 0 ∈ R N implies that for any x ∈ R N we have K p x ∈ E 0 . Proof:[of Proposition 2.3.2]
We use the generalized entropy method, see [START_REF] Perthame | Transport equations in biology[END_REF]. The proof is in several steps.

First step. Construction of a Floquet eigenproblem.

We first apply the Floquet theorem to K p (θ), which is T -periodic in θ. The resolvent of K p (θ), denoted R(K p )(t, s) to avoid confusion with the resolvent R p (t, s) of K p (θ), is of the form P (ts) exp((ts)A p ) for some T -periodic and invertible matrix P (θ) satisfying P (0) = Id, and for some constant coefficients matrix A p . Besides, standard ODE considerations show that the resolvent of K p (θ) preserves componentwise positivity of vectors in R N (this is a standard fact concerning linear Boltzmann-like operators). Therefore the Perron-Frobenius theory applies and the spectral radius of R(K p )(T, 0), namely of the resolvent of K p (θ) over one period, is a simple eigenvalue of R(K p )(T, 0), for which there exists a unique associated eigenvector having positive entries which sum up to one. Let µ per be this spectral radius, and p per the associated eigenvector. We have by definition

µ per = ρ(R(K p )(T, 0)) = ρ(exp(T A p )).
Let us now prove that µ per = 1. From the equality 1 T K p (θ) = 0, valid for any θ, it comes 1 T exp(T A p ) = 1 T and 1 is an eigenvalue of exp(T A p ). Besides, taking λ > 0 large enough, the matrix K p (θ) + λ has positive entries, uniformly in θ, hence so does the resolvent R(K p + λId)(t, s) = R(K p )(t, s) exp(λ(ts)).

In particular exp(T A p ) = R(K p )(T, 0) has positive entries. Therefore by the Gershgorin Theorem we have µ per ≤ max i 1 T exp(T A p ) i and we deduce, using the equality 1 T exp(T A p ) = 1 T , that µ per = 1.

From this comes the existence of a (unique) T -periodic solution p(θ) to the direct problem

d dθ p(θ) = K p (θ) p(θ), with the normalization T 0 N i=1 p i (θ)dθ = 1.
Indeed p(θ) is (up to normalization) equal to R(K p )(θ, 0) p per , where p per is the above mentionned eigenvector associated with µ per = 1. The vector p(θ) has positive entries for any θ. Similarly, by a simple transposition, there exists a (unique) T -periodic solution φ(θ) to the dual problem

d dθ φ(θ) = φ(θ) K p (θ), with the normalization T 0 N i=1 p i (θ) φ i (θ) dθ = 1.
We actually have φ(θ) = 1 T for any θ, and φ(θ) has positive entries for any θ as well.

In particular and for later convenience note that there is a constant c 0 > 0 such that for any t ∈ [0, T ] and any index i we have

0 < c 0 ≤ p i (t) ≤ 1/c 0 and 0 < c 0 ≤ φ i (t) ≤ 1/c 0 .
(2.3.3)

Second step. The generalized entropy inequality.

The key point is to use the generalized entropy equality. It asserts, amongs others, that any solution x(t) to dx(t)/dt = K p (t) x(t), satisfies the equality

d dt N i=1 φ i (t) p i (t) x i (t) p i (t) 2 = - N i, j = 1, i = j φ i (t) (K p (t)) i,j p j (t) x j (t) p j (t) - x i (t) p i (t) 2 , (2.3.4) upon setting x i (t) = x i (t) -ρp i (t), and ρ = T 0 N i=1 x i (t) φ i (t) = T N i=1
x i (0) .

(2.3.5)

We refer to [START_REF] Perthame | Transport equations in biology[END_REF] for the proof, which turns out to be an explicit computation (the difficult point here is to introduce the correct convex functional of interest, here the weighted, time-dependent l 2 norm . 2 := i φ i (t)p i (t) (./p i (t)) 2 , that is based on the periodic solutions p(t) and φ(t)). In (2.3.5) we used φ(t) = 1 T . From (2.3.4) we deduce, setting c 1 = min i =j min t∈[0,T ] σ p i,j (t) > 0, (see assumptions (2.2.2)-(2.2.3)), the bound

d dt N i=1 φ i (t) p i (t) x i (t) p i (t) 2 ≤ -c 1 N i,j=1 φ i (t) p j (t) x j (t) p j (t) - x i (t) p i (t) 2 .
(2.3.6) Note that the vector x(t) in (2.3.6) satisfies the normalisation

T 0 N i=1
x i (t) dt = 0, thanks to the conditions satisfied by p(t) and φ(t).

Next, we assert that there is a c 2 > 0 such that any solution x(t) to dx(t)/dt = K p (t)x(t) satisfies, with the above notation, and for any t ∈ [0, T ] (i.e. over one period), the lower bound N i,j=1 φ i (t) p j (t)

x j (t)

p j (t) - x i (t) p i (t) 2 ≥ c 2 N i=1 φ i (t) p i (t) x i (t) p i (t) 2 .
(2.3.7)

We postpone the proof of (2.3.7) for the time being.

Putting together (2.3.6) and (2.3.7) provides, thanks to the Gronwall Lemma, that for any t ∈ [0, T ] and for any solution to dx(t)/dt = K p (t)x(t) we have We now pick up an arbitrary solution x(t) to dx(t)/dt = K p (t)x(t) with x(0) ∈ E 0 . Since by Proposition 2.3.1 the subspace E 0 is stable under the action of K p (θ) for any θ, it comes that x(t) actually satisfies dx(t)/dt = K p (t)x(t) with x(0) ∈ E 0 . Hence x(t) = R p (t, 0)x(0) with the notation of the Proposition.

N i=1 φ i (t) p i (t) x i (t) p i (t) 2 ≤ exp(-c 1 c 2 t) N i=1 φ i (0) p i (0) x i (0) p i (0) 2 . ( 2
On the other hand, inequality (2.3.8) implies, using the fact that x(t) ∈ E 0 , which in turn implies x(t) ≡ x(t), that for any t ≥ 0 we have

N i=1 φ i (t) p i (t) x i (t) p i (t) 2 ≤ exp(-c 1 c 2 t) N i=1 φ i (0) p i (0) x i (0) p i (0) 2 .
Using the uniform bound (2.3.3) then provides

N i=1 [x i (t)] 2 ≤ exp(-c 1 c 2 t) (c 0 ) 4 N i=1 [x i (0)] 2 .
Estimate (2.3.2) follows.

Fourth step. Proof of the Poincaré inequality (2.3.7).

We argue by contradiction. In the case where (2.3.7) is not valid, there is a sequence of solutions x (k) (t) to

dx (k) (t)/dt = K p (t)x (k) (t)
, and a sequence of times

t (k) ∈ [0, T ], with N i=1 φ i (t (k) ) p i (t (k) ) x (k) i (t (k) ) p i (t (k) ) 2 = 1
for each k, together with the convergence

N i,j=1 φ i (t (k) ) p j (t (k) ) x (k) j (t (k) ) p j (t (k) ) - x (k) i (t (k) ) p i (t (k) ) 2 → 0 as k → ∞.
Under these circumstances, and thanks to the bounds (2.3.3), it comes that the sequence of functions x (k) (t), defined over [0, T ], is bounded in C 0 ([0, T ], R N ), uniformly in k. On the other hand, an elementary computation shows that d x (k) (t)/dt = K p (t) x (k) (t), for any k. Hence x (k) (t) is uniformly bounded in C 1 ([0, T ], R N ), and actually in C n ([0, T ], R N ) for any integer n. Therefore, possibly taking subsequences, and using the Ascoli-Arzela Theorem, we may assume t (k) → t * for some t * , and that the sequence of x (k) (t) converges towards some function y(t) in, say, C 1 ([0, T ], R N ). Naturally, the function y(t) satisfies dy(t)/dt = K p (t) y(t). Besides, passing to the limit in the various relations satisfied by the sequence x (k) (t), we obtain

T 0 i y i (t) dt = 0, together with N i=1 φ i (t * ) p i (t * ) y i (t * ) p i (t * ) 2 = 1,

Reduction of the system

Starting from differential system (2.2.1), we are now in position to develop a time dependent center manifold approach to analyze the limit ε → 0. This is a natural procedure. Indeed, according to the previous paragraph, our transport matrices K p (θ) and K q (θ), which carry the stiff factor 1/ε in (2.2.1), possess a spectral decomposition such that the eigenspace associated with the zero eigenvalue has constant dimension, and there is a stable supplementary space (here E 0 , but this space could possibly depend on θ if necessary) such that the eigenvalues of the correspondingly induced matrices have all eigenvalues with real part ≤ -2β < 0. Even more, the resolvent of the associated differential equation is exponentially decaying. For that reason, it is natural to look for a θ-dependent change of variables (p, q) → (x, z) in (2.2.1), where variable x (the "slow" variable) is roughly associated with the zero eigenvalues of our transport matrices, while z (the "fast" variable) is roughly associated with the eigenvalues with negative real part, and such that the whole nonlinear system (2.2.1) may be recast in the form

dx dt = F x, z, t ε , dz dt = 1 ε B t ε z + G x, z, t ε , (2.3.9)
where the θ-dependent matrix B(θ) is periodic, has eigenvalues with negative real part, and possesses an exponentially decaying resolvent.

The main result of this paragraph is the following Proposition 2.3.4 (Reduction of the original system to a normalized form) Consider the system (2.2.1), with unknowns (p ε (t), q ε (t)) ∈ R 2N . Set the macroscopic variable (we drop the ε-dependence for convenience)

x(t) = (x p (t), x q (t)) T ∈ R 2 , where x p (t) = N i=1 p ε i (t) (=total number of preys),
and

x q (t) = N i=1 q ε i (t) (=total number of predators),
Set also the linear operator (we use the notation of Proposition 2.3.2)

B(θ) = K p (θ) 0 0 K q (θ) ,
acting on E 0 ×E 0 . According to Proposition 2.3.2, the resolvent of B(θ), denoted by R(t, s) for convenience, decays exponentially with time, in that there is a µ 0 > 0 and a C 0 ≥ 1 such that for any t ≥ s we have

R(t, s) ≤ C 0 e -µ0(t-s) .
With this notation, there exists a one-to-one linear mapping Φ t ε : R 2N → R 2 ×(E 0 ) 2 , that depends smoothly and periodically in the variable t/ε, and there exist smooth and explicitely computable functions F (x, z, θ) and G(x, z, θ), defined on R 2 × (E 0 ) 2 × T, such that, defining the change of variables

Φ t ε : p ε q ε ∈ R N × R N → x z ∈ R 2 × E 2 0 the system (2.2.1) transforms into        dx(t) dt = F x(t), z(t), t ε , dz(t) dt = 1 ε B t ε z(t) + G x(t), z(t), t ε .
(2.3.10)

Remark 2.3.5 Once system (2.2.1) is recast as (2.3.10), we shall use the results of Section 1.1.

Proof:[of Proposition 2.3.4] Our strategy of proof is the following. Given that matrix K p (θ) has a simple eigenvector p eq (θ) associated to 0 and other eigenvalues lower than a constant -β (Proposition 2.3.1), our first step consists in treating separately the equation projected on p eq (θ) and the equations projected on E 0 = {z ∈ R N , z T 1 = 0}, and similarly for the index "q". Our second step consists in removing the remaining stiff term through a time-dependent change of variables. This is the key non-trivial step.

First step. Projecting on the natural almost-equilibrium.

In view of the spectral decompositions Span(p eq (θ)) ⊕ E 0 = R N , together with Span(q eq (θ)) ⊕ E 0 = R N , we introduce the decomposition

x p (t) = 1 T p ε (t) ∈ R and y p (t) = p ε (t) -x p (t)p eq t ε ∈ E 0 ,
and similarly for q ε (t). Needless to say, x p (t) and x q (t) denote, respectively, the total number of preys resp. predators, over all sites i = 1, ..., N . The supplementary variables y p (t) and y q (t) simply denote the difference with the natural almost-equilibrium repartition p ε (t) = x p (t) p eq t ε , and similarly for the "q" index. We have the identities p ε (t) = x p (t)p eq t ε + y p (t) ∈ Span p eq t ε ⊕ E 0 , q ε (t) = x q (t)q eq t ε + y q (t) ∈ Span q eq t ε ⊕ E 0 ,

In the same spirit, we also introduce the decomposition f x 1 (x p , x q , y p , y q , θ) = 1 T f (x p p eq (θ) + y p , x q q eq (θ) + y q ) ∈ R, f y 1 (x p , x q , y p , y q , θ) = f (x p p eq (θ) + y p , x q q eq (θ) + y q ) f x 1 (x p , x q , y p , y q , θ) p eq (θ) ∈ E 0 , and similarly for g x 1 and g y 1 . With obvious shorthand notation we have

f (p ε (t), q ε (t)) = f x 1 p eq t ε + f y 1 ∈ Span p eq t ε ⊕ E 0 , g (p ε (t), q ε (t)) = g x 1 q eq t ε + q y 1 ∈ Span q eq t ε ⊕ E 0 .
Still in the same direction, we observe that the relation 1 T p eq (θ) = 1, valid for any θ, leads upon differentiation to ∀θ, ṗeq (θ) ∈ E 0 , and similarly for qeq (θ), where the dot denotes differentiation.

Projecting system (2.2.1) according to the spectral decompositions Span(p eq (θ))⊕E 0 = R N , and Span(q eq (θ))⊕

E 0 = R N , clearly leads to (see Proposition 2.3.2)                              ẋp = f x 1 x p , x q , y p , y q , t ε , ẏp = 1 ε K p t ε y p - 1 ε x p ṗeq t ε + f y 1
x p , x q , y p , y q , t ε , ẋq = g x 1 x p , x q , y p , y q , t ε ,

ẏq = 1 ε K q t ε y q - 1 ε x q qeq t ε + g y 1 x p , x q , y p , y q , t ε .
(2.3.11)

Second step. Getting rid of the additional stiff terms.

System (2.3.11) still is not of the desired form (2.3.10). We need to get rid of the stiff terms in 1 ε x p ṗeq and 1 ε x q qeq . This is the reason why we now introduce a change of variables y p → z p , of the form

z p (t) := y p (t) -h 0 p x p (t), t ε (2.3.12)
(and similarly for the index "q"), where the -still unknown -function h 0 p (x, θ) is required to be periodic in θ (and we assume for the time being that h 0 p (x, θ) ∈ E 0 for any x ∈ R and θ ∈ T). The aim is to obtain, for the variable z p , a differential equation of the form

żp (t) = 1 ε K p t ε z p (t) + f x p , x q , z p , z q , t ε ,
where f is a function without any pre-factor 1/ε. We readily stress the point that the change of variables (2.3.12) in fact is guided by a central manifold point of view, and the seeked function h 0 p may be interpreted as a first order expansion of a center manifold for the (fast) variable y p in (2.3.11).

In any circumstance, starting from the Ansatz (2.3.12), and differentiating w.r.t time t, leads to

żp (t) = ẏp (t) -∂ x h 0 p x p (t), t ε ẋp (t) - 1 ε ∂ θ h 0 p x p (t), t ε = 1 ε K p t ε z p (t) -∂ x h 0 p x p (t), t ε f x 1 x p , x q , y p , y q , t ε + f y 1 x p , x q , y p , y q , t ε - 1 ε ∂ θ h 0 p x p (t), t ε + 1 ε K p t ε h 0 p x p (t), t ε - 1 ε x p ṗeq ,
from which it becomes clear that h 0 p (x, θ) should be taken as a periodic solution of the following equation

∂ θ h 0 p (x, θ) = K p (θ) h 0 p (x, θ) -x ṗeq (θ).
In order to solve the previous equation, we use the resolvent R p (t, s) associated with the operator K p (θ) acting on E 0 . The Duhamel formula provides

h 0 p (x, θ) = R p (θ, 0) h 0 p (•, 0) - θ 0 R p (θ, ϕ)x ṗeq (ϕ)dϕ.
Since h 0 p is assumed T -periodic with period T , the following relation is necessary

(Id -R p (T, 0))h 0 p (x, 0) = - T 0 R p (T, ϕ)x ṗeq (ϕ)dϕ.
Since Id -R p (T, 0) is invertible (because of the spectral properties of R p (T, 0) obtained in Proposition 2.3.2), the only solution of the previous equation is given by

h 0 p (x, 0) = -(Id -R p (T, 0)) -1 T 0 R p (T, ϕ)x ṗeq (ϕ)dϕ ∈ E 0 .
This leads us to define h 0 , for any (x, θ) ∈ R × T as

h 0 p (x, θ) := -x R p (θ, 0)(Id -R p (T, 0)) -1 T θ-T R p (0, ϕ) ṗeq (ϕ)dϕ ,
or, in other words

h 0 p (x, θ) := -x S p ( ṗeq ) (θ)
, where we introduce the nonlocal operator S p

S p : v = v(θ) → S p (v)(θ) := R p (θ, 0)(Id -R p (T, 0)) -1 T θ-T R p (0, ϕ)v(ϕ)dϕ .
Note that S p maps C 0 (T) to C 0 (T), but also that S p maps C 0 (T; E 0 ) to C 0 (T; E 0 ). Since ṗeq (θ) ∈ E 0 for any θ, this provides that S p ( ṗeq )(θ) ∈ E 0 as well. In some sense we have

S p = ∂ θ -K p (θ) -1
, in that for each periodic function v, operator S p computes the unique periodic solution w to

∂ θ -K p (θ) w = v.
In any case, all this defines the quantity h 0 p (x, θ) ∈ E 0 and the operator S p and we obviously introduce the analogous quantities for the index "q". Eventually, the change of variables y p (t) → z p (t) and y q (t) → z q (t), from E 0 to E 0 , that we introduce now is z p (t) = y p (t) + x p (t) S p ( ṗeq ) t ε , z q (t) = y q (t) + x q (t) S q ( qeq ) t ε .

(2.3.13)

With this new set of variables, the system (2.3.11) becomes

                             ẋp = f x 2 x p , x q , z p , z q , t ε , żp = 1 ε K p t ε z p + f y 2 x p , x q , z p , z q , t ε , ẋq = g x 2 x p , x q , z p , z q , t ε , żq = 1 ε K q t ε z q + g y 2 x p , x q , z p , z q , t ε , (2.3.14)
and the following definitions are used f x 2 (x p , x q , z p , z q , θ) = f x 1 x p , x q , z p + h 0 p (x p , θ), z q + h 0 q (x q , θ), θ , f y 2 (x p , x q , z p , z q , θ) = f y 1 x p , x q , z p + h 0 p (x p , θ), z q + h 0 q (x q , θ), θ + f x 2 (x p , x q , z p , z q , θ) S p ( ṗeq ) (θ), together with the equivalent definitions for g x 2 and g y 2 . Lastly, denoting

x = (x p , x q ) T , z = (z p , z q ) T , and introducing the abuse of notation

f x 2 ≡ f x 2 (x, z, θ), f y 2 ≡ f y 2 (x, z, θ), g x 2 ≡ g x 2 (x, z, θ), g y 2 ≡ g y 2 (x, z, θ),
as well as the notation

F (x, z, θ) = f x 2 (x, z, θ) g x 2 (x, z, θ) ∈ R 2 , G(x, z, θ) = f y 2 (x, z, θ) g y 2 (x, z, θ) ∈ E 0 × E 0 ,
system (2.3.14) may be written in the shorter form

         ẋ = F x, z, t ε , x (t 0 ) = x 0 , ż = 1 ε B t ε z + G x, z, t ε , z (t 0 ) = z 0 .
Third step. Summarizing.

All in all, the t/ε-dependent change of variables

(p ε , q ε ) ∈ R 2N → (x, z) ∈ R 2 × (E 0 ) 2 that transforms (2.2.1)
into the reference system (2.3.10) reads, with the above notation x = (x p , x q ), z = (z p , z q ), with (2.3.15)

x p = 1 T p ε = N i=1 p ε i , z p = p ε -x p p eq t ε -S p ( ṗeq ) t ε ,
and similarly for "q". The inverse mapping reads

p ε = z p + x p p eq t ε -S p ( ṗeq ) t ε ,
and similarly for "q". We recall the definition, valid for any periodic function v(θ),

S p (v)(θ) := R p (θ, 0)(Id -R p (T, 0)) -1 T θ-T R p (0, ϕ)v(ϕ) dϕ , (2.3.16)
and similarly for "q", which in some sense means

S p = ∂ θ -K p (θ) -1
. These formulae entirely define the isomorphism Φ t ε : R 2N → R 2 × E 2 0 that is referred to in the statement of Proposition 2.3.4, as well as its inverse. It depends smoothly and periodically in the variable t/ε, and is entirely based on the linear operator B(θ) and on associated spectral objects.

With this notation the nonlinearities F and G in (2.3.10) simply are conjugated to the original nonlinearities f and g in (2.2.1) through the isomorphism Φ, through the relation

F G (x, z, t/ε) = Φ t ε -1 • f g • Φ t ε (p ε , q ε ).
Since the original nonlinearities are smooth, so are F and G. In a case where the original nonlinearities f and g would have limited C r smoothness, then so would have F and G as well.

Proposition 2.3.4 has the following immediate corollary, which we state here for convenience.

Corollary 2.3.6 (Existence of solutions to the original system on a uniform time interval) Fix R > 0. Then, there exists ε 0 = ε 0 (R) such that the following holds.

There exists a T (R) > 0, such that for any 0 < ε < ε 0 (R), and for any initial data (p ε (0), q ε (0)) in (2.2.1) satisfying (p ε (0), q ε (0)) ≤ R, the unique solution (p ε (t), q ε (t)) to (2.2.1) associated with the initial choice (p ε (0), q ε (0)) is well-defined on the time-interval [0, T (R)], independently of ε.

Proof:[of corollary 2.3.6]

Take R > 0. Proposition 2.3.4 asserts that we may equivalently argue on the reduced system (2.3.10). Therefore we may assume that the initial datum (x(0), z(0)) in (2.3.10) belongs to the ball B c1R of radius c 1 R in R 2 × E 2 0 , where c 1 > 0 is a constant that only depends on the function B(θ) through the isomorphism Φ t ε . Let assume for sake of brevity that c 1 = 1. Denote by M = M (R) a bound on the function (F, G)(x, z, θ) on B 2 C0 R × T, where C 0 is the constant in Proposition 2.3.4. As long as the solution (x(t), z(t)) to (2.3.10) belongs to the ball B 2 C0 R , we may write, using the Duhamel formula and the exponential smallness of the resolvent R(t, s),

z(t) ≤ C 0 e -µ0t/ε R + t 0 C 0 e -µ0(t-s)/ε M ds ≤ C 0 R + C 0 µ -1 0 ε M ≤ 2 C 0 R
as long as we choose ε ≤ µ 0 /M. On top of that, we may as well estimate

x(t) ≤ R + t 0 M ds ≤ R + M t ≤ 2 C 0 R
as long as we choose 0 ≤ t ≤ C 0 /M. The result follows.

Remark 2.3.7 A technically important consequence of the above corollary, besides the existence of a common time interval on which all solutions (p ε (t), q ε (t)) to (2.2.1) are well-defined, is the following. For any given R > 0, up to restricting our attention on a certain, fixed, possibly small, time interval [0, T (R)], we may always assume that the function (p ε (t), q ε (t)) to (2.2.1), or equivalently (x(t), z(t)) to (2.3.10), belongs to the ball of radius R for all times. Note that it may happen that T (R) = +∞ if the structure of the equations makes sure that the solution remains in the ball of radius R for all times.

In particular, given R > 0, we may always introduce a smooth truncation χ R on the space of variables (p ε (t), q ε (t)), or equivalently in the variables (x, z), a truncation which is one on the ball of radius R, and zero outside the ball of radius 2R, and consider the truncated nonlinearities (f, g) or (F, G) that coincide with the original nonlinearities on the ball of radius R, and vanish outside the ball of radius 2R. With that truncation in mind, we see that the functions (p ε (t), q ε (t)), or (x(t), z(t)) may always be seen as the restrictions, on the time interval [0, T (R)] for some 0 < T (R) ≤ +∞, of solutions to a system of the form (2.2.1) or (2.3.10) that are defined on the whole time interval [0, +∞[, and are associated with nonlinearities that are bounded as well as all their derivatives, globally.

A center manifold approach

In this section, we apply the theoretical results of Section 1.1 to the system (2.3.10):

         ẋ = F x, z, t ε , x (t 0 ) = x 0 , ż = 1 ε B t ε z + G x, z, t ε , z (t 0 ) = z 0 ,
where x = (x p , x q ) ∈ R 2 is the number of preys and predators, irrespective of the site. We have separated a slow and a fast variable, and the linear operator B has an exponentially small resolvent.

We apply Theorem 1.1.2 to (2.3.10). We refer to Remark 2.3.7: though the functions F and G stemming from our original equations (2.2.1) do not satisfy the boundedness assumptions listed in the present Theorem, we may nonetheless apply the Theorem in this case, on any restricted time interval [0, T (R)] where we can ensure that the solution belongs to a given ball of radius R. Hence, there exists ε 0 > 0, and a function εh(x, θ)

∈ C 1 (R n × T) ∩ C 0 (T; C r (R n )) (which is an O(ε))
, defined for all 0 < ε < ε 0 , with the following property. For all x 0 ∈ R n and θ 0 ∈ T, the solution (x(t), θ(t), z(t)) of (2.3.10) with initial conditions

x(0) = x 0 , θ(0) = θ 0 , z(0) = εh(x 0 , θ 0 ), satisfies the relation, for all t, z(t) = εh x(t), θ 0 + t ε .
(2.4.1)

Remark 2.4.1 In the previous section, we have first introduced an intuitive new variable y, and then defined z in (2.3.12) to get rid of additional stiff terms. Hence,

y(t) = z(t) + h 0 x(t), t ε .
Then, using (2.4.1) (with θ 0 = 0), we have :

y(t) = h 0 x(t), t ε + εh x(t), t ε .
This equality explains the choice of the name h 0 : this function may be interpreted as a zero order expansion of a center manifold.

Then, Theorem 1.1.4 shows the exponential convergence of (x(t), z(t)) towards the center manifold, up to a change of initial condition x 0 to x ε 0 . Indeed, for any fixed T ∞ > 0, the solution components of the reduced system

         ẋh = F x h , εh x h , θ 0 + t ε , θ 0 + t ε , x h (0) = x ε 0 z h (t) = εh x h (t), θ 0 + t ε satisfy the following error estimate ∀t ∈ [0, T ∞ ], z(t) -z h (t) + x(t) -x h (t) ≤ Ce -µt/ε .
Hence, we can reduce the stiff dynamics of the 2N unknowns (p ε (t), q ε (t)), to a non-stiff dynamics, with fast oscillating coefficients, of the 2 unknowns (x p (t), x q (t)), the total number of preys and predators.

Moreover, according to Theorem 1.1.6, the center manifold εh can be approximate to every order in ε. For all r ≥ 1, there exists

h [r] ε (x, θ) := εh 1 (x, θ) + ε 2 h 2 (x, θ) + • • • + ε r h r (x, θ), with h 1 (x, θ), . . . , h r (x, θ) periodic in θ, such that εh -h [r] ε ∞ ≤ C r ε r+1 .
Eventually, we apply Theorem (1.2.4) to average the fast dynamics.

Qualitative behaviour 2.5.1 Derivation of the first terms of the expansion

Applying the results of Section 1.1.4 to our initial Lotka-Volterra system with fast migrations (2.2.1), we come up with the following formulae.

The zeroth order reduced system associated with (2.2.1) is

         ẋp (t) = A 0 p t ε x p (t) -B 0 p t ε x p (t)x q (t), ẋq (t) = -A 0 q t ε x q (t) + B 0 q t ε x p (t)x q (t),
(2.5.1) with

y p (t) = -x p (t) S p ( ṗeq ) t ε , y q (t)(t) = -x q (t) S q ( qeq ) t ε .
Here we have defined the coefficients A 0 p (θ) = 1 T (a p p eq (θ)) , B 0 p (θ) = 1 T (b p p eq (θ) q eq (θ)) , A 0 q (θ) = 1 T (a q q eq (θ)) , B 0 q (θ) = 1 T (b q p eq (θ) q eq (θ)) , and we set p eq (θ) = p eq (θ) -S p ( ṗeq ) (θ), q eq (θ) = q eq (θ) -S q ( ṗeq ) (θ).

We also used as in (2.2.5) the convention that a p p eq etc. denote the componentwise vector products (a p,i p eq,i

) N i=1 ∈ R N etc.
The first order reduced system associated with (2.2.1) is

                           ẋp (t) = A 0 p + εA 1 p t ε x p (t) -B 0 p + εB 1 p t ε x p (t) x q (t) +ε C 1 p t ε x p (t) 2 x q (t) + ε D 1 p t ε x p (t) x q (t) 2 , ẋq (t) = -A 0 q + εA 1 q t ε x q (t) + B 0 q + εB 1 q t ε x p (t) x q (t) -ε C 1 q t ε x q (t) 2 x p (t) -ε D 1 q t ε x q (t) x p (t) 2 , (2.5.2) with          y p (t) = -x p (t) S p ( ṗeq ) t ε + ε x p (t) E 1 p t ε -x p (t) x q (t) F 1 p t ε , y q (t)(t) = -x q (t) S q ( qeq ) t ε + ε -x q (t) E 1 q t ε + x p (t) x q (t) F 1 q t ε .
(2.5.3)

Here we have introduced the following coefficients. On the one hand we defined

E 1 p (θ) = S p Π peq (θ) (a p p eq (θ)) (θ), F 1 p (θ) = S p Π peq (θ) (b p p eq (θ)
q eq (θ)) (θ), E 1 q (θ) = S q Π qeq (θ) (a q q eq (θ)) (θ), F 1 q (θ) = S q Π qeq (θ) (b q p eq (θ) q eq (θ)) (θ), where we have defined the projections on E 0 parallel to p eq (θ) resp. parallel to q eq (θ)

Π peq (θ) : v ∈ R N → v -1 T v p eq (θ), Π qeq (θ) : v ∈ R N → v -1 T v q eq (θ).
Note that the θ-dependent coefficients E 1 p , F 1 p , E 1 q and F 1 q all belong to E 0 for each θ. On the other hand we defined

A 1 p (θ) = 1 T a p E 1 p (θ) , B 1 p (θ) = 1 T b p -p eq (θ) E 1 q (θ) + q eq (θ) E 1 p (θ) + a p F 1 p (θ) , C 1 p (θ) = 1 T -b p p eq (θ) F 1 q (θ) , D 1 p (θ) = 1 T b p q eq (θ) F 1 p (θ) , as well as A 1 q (θ) = 1 T -a q E 1 q (θ) , B 1 q (θ) = 1 T b q -p eq (θ) E 1 q (θ) + q eq (θ) E 1 p (θ) -a q F 1 q (θ) , C 1 q (θ) = 1 T b q q eq (θ) F 1 p (θ) , D 1 q (θ) = 1 T -b q p eq (θ) F 1 q (θ)
. System (2.5.2) is a fast time dependent Lotka-Volterra system, wuth cubic corrective terms of size ε, in the variables (x p (t), x q (t)). As we shall see below, the cubic corrections carry the key qualitative feature of the asymptotic dynamics.

Application of the averaging results

In the case where we start from system (2.2.1), the formulae of Section 1.2.2 take the following form.

The zeroth order averaged reduced system associated with (2.2.1) is

ẋp (t) = A 0 p x p (t) -B 0 p x p (t)
x q (t), ẋq (t) = A 0 q x q (t) -B 0 q x p (t)x q (t),

(2.5.4) with

A 0 p = 1 T T 0 A 0 p (θ) dθ,
and similarly for the other coefficients. At zero order, the averaged system is thus simply the zero order reduced system (2.5.1), where the oscillatory coefficients A 0 p (θ) etc. are averaged out. Needless to say, equation (2.5.4) is a Lotka-Volterra system on the global predator-prey populations (x p , x q ), whose coefficients A 0 p etc. are obtained as appropriate averages, both in the sites i = 1, . . . , N and in the fast ocillating variable θ ∈ T, of the original Lotka-Volterra dynamics. It is implicit here that the actual repartition of preys over the sites is, at zero order, given by the relation y p = x p S p ( ṗeq ) t ε , and similarly for "q", see (2.5.3). The first order averaged reduced system associated with (2.2.1) is

             ẋp (t) = A 0 p + ε A p 1 x p (t) -B 0 p + ε B 1 p + ε B 1 p x p (t) x q (t) +ε C 1 p + C 1 p x p (t) 2 x q (t) + ε D 1 p x p (t) x q (t) 2 , ẋq (t) = -A 0 q + ε A q 1 x q (t) + B 0 q + ε B 1 q + ε B 1 q x p (t) x q (t) -ε C 1 q + C 1 q x q (t) 2 x p (t) -ε D 1 q x q (t) x p (t) 2 , (2.5.5)
where as before . means averaging out in θ, and we have defined the new coefficients

B 1 p = 1 2T T 0 θ 0 B 0 p (θ ′ ) A 0 q (θ) -A 0 q (θ ′ ) B 0 p (θ) dθ ′ dθ, B 1 q = - 1 2T T 0 θ 0 B 0 q (θ ′ ) A 0 p (θ) -A 0 p (θ ′ ) B 0 q (θ) dθ ′ dθ, C 1 p = 1 2T T 0 θ 0 B 0 p (θ ′ ) B 0 q (θ) -B 0 q (θ ′ ) B 0 p (θ) dθ ′ dθ, C 1 q = - 1 2T T 0 θ 0 B 0 q (θ ′ ) B 0 p (θ) -B 0 p (θ ′ ) B 0 q (θ) dθ ′ dθ = C 1 p
As in (2.5.2), the first order averaged and reduced system (2.5.5) is a Lotka-Volterra system, wuth cubic corrective terms of size ε, in the variables (x p (t), x q (t)), whose coefficients are obtained as appropriate averages of the original coefficients in both the index i and the variable θ. The cubic corrections carry the key qualitative feature of the asymptotic dynamics.

Qualitative analysis of the zero and first order averaged reduced systems

At zero order, the global dynamics of (2.2.1) is thus foreseen to be described by (2.5.5). Qualitatively, this means that the global predator-prey populations oscillate around an equilibrium point in a periodic fashion. The important point is now the effect of the first order corrective terms. The corrections intervening in terms of degree one and two in (2.5.5) do not modify the global dynamics (which still is of Lotka-Volterra type as far as these terms are concerned), as foreseen by the zero order asymptotic model (2.5.4). yet the cubic terms do rule out this qualitative property: though the implicit function theorem still foresees the existence of an equilibrium point, in general the cycles around the equilibrium point will be broken. And the question is: do we observe a stabilization or a destabilization of the cycles foreseen at zero order?

In other words, we need to qualify the stability of system (2.5.5) around its equilibrium. This we do in the spirit of Poggiale discussion in [START_REF] Poggiale | Lotka-Volterra's model and migrations: Breaking of the well-known center[END_REF]. Indeed, the system being posed in dimension two, a standard geometric criterion applies, for small values of ε > 0. Writing (2.5.4) as ẋ = F 0 (x) and system (2.5.5) as ẋ = F 0 + ε F 1 (x), since (2.5.4) has an equilibrium point F 0 (x 0 ) = 0 at some x 0 = (x 0 p , x 0 q ) T ∈ R * + 2 (whose exact value is inessential), it is clear that (2.5.5) has an equilibrium point at some x ε ∈ R * + 2 (whose exact value is inessential as well). Now, since x 0 is actually a center for the vector field F 0 , we have

Tr ∂ x F 0 (x 0 ) = 0.
This is beacause (2.5.4) is a Lotka-Volterra system. Therefore we know that

Tr ∂ x F 0 + ε F 1 (x ε ) = λ ε + O(ε 2 ),
for some given real coefficient λ. In two dimension the following criterion therefore applies:

if λ > 0, system (2.5.5) has an unstable focus at x ε , if λ < 0, system (2.5.5) has a stable focus at x ε .

We cannot conclude if λ = 0. On the other hand, an easy computation provides

λ = x 0 p x 0 q C 1 p + C 1 p -C 1 q + C 1 q =: x 0 p x 0 q σ.
This serves as a definition for

σ = C 1 p -C 1 q ,
and all in all we have the value σ = -1 T b p p eq S q Π qeq b q p eq q eq + b q q eq S p Π peq b p p eq q eq .

We conclude that if σ > 0, system (2.5.5) has an unstable focus at x ε , if σ < 0, system (2.5.5) has a stable focus at x ε .

And the above Theorems allow to conclude that this stability resp. instability criterion holds as well for the total populations (x p , x q ) in the original model (2.2.1).

It is important to note that the naive method to deal with the original model (2.2.1), consisting in first freezing the coefficients of the matrices K p (θ) and K q (θ), say to their mean value K p and K q , and next performing the natural central manifold analysis, leads to a system of the form (2.5.5) yet with the wrong value of σ. The original model (2.2.1) may very well be stable for small values of ε while the modified model with frozen coefficients may be unstable, and conversely. We give an example below.

One example with N = 2 sites

We apply the above method on one simple example with two sites. Our starting equation is

dp ε dt = 1 ε K p t ε p ε + f (p ε , q ε ), dq ε dt = 1 ε K q t ε q ε + g(p ε , q ε ),
an equation that is posed in two dimensions, i.e. p ε and q ε both belong to R 2 .

Stability issues.

Here we choose

K p = -1, K q = -1, p eq (θ) = 1 -a(θ) a(θ) , q eq = 1 -b b , which means K p (θ) = -a(θ) 1 -a(θ) a(θ) -(1 -a(θ)) , K q (θ) = -b 1 -b b -(1 -b)) ,
together with a(θ) = a 0 + a 1 cos(θ) + a -1 sin(θ),

and b p = b 1 p b 2 p , b q = b 1 q b 2 q .
An easy computation provides, in this particular case,

C 1 p = -b (1 -b) b 1 p (1 -a 0 ) -b 2 p a 0 b 1 q (1 -a 0 ) -b 2 q a 0 + (b 1 p + b 2 p ) (b 1 q + b 2 q ) 8 (a 2 1 + a 2 -1 ) , C 1 q = b 1 p (1 -b) -b 2 p b b 1 q (1 -b) -b 2 q b a 0 (1 -a 0 ) - (a 2 1 + a 2 -1 ) 4 .
Next we select the following simple values, leaving b 2 p as the only variable for the time being, namely

a 0 = b = 1/2, a 1 = 0.1, a -1 = 0.37, b 1 p = 0.2, b 1 q = 0.5, b 2 q = 1.5. We have σ = 0.0036 > 0 if b 2 p = 0.
4, meaning instability of the equilibrium point in that case, while σ = -0.0055 < 0 if b 2 p = 0.7, meaning stability. In the similar spirit, note that if one uses the wrong procedure consisting in first replacing K p and K q by their mean value, next performing the central manifold approach, one obtains a wrong prediction σ wrong whose value is

σ wrong = -b (1 -b) b 1 p (1 -a 0 ) -b 2 p a 0 b 1 q (1 -a 0 ) -b 2 q a 0 -b 1 p (1 -b) -b 2 p b b 1 q (1 -b) -b 2 q b a 0 (1 -a 0 ),
With the same choice of parameters as above, we have σ = -0.0055 < 0 if b 2 p = 0.7, meaning stability of the original system, while σ wrong = 0.0075 > 0 for the same value of b 2 p , meaning that freezing the oscillatory coefficients at once wrongly predicts instability of the system.

Here, we present an explicit computation of the sign of the coefficients σ and σ wrong in that case. According to the previous calculus:

σ = C 1 p -C 1 q = -b (1 -b) b 1 p -(b 1 p + b 2 p ) a 0 b 1 q -(b 1 q + b 2 q ) a 0 -a 0 (1 -a 0 ) b 1 p -(b 1 p + b 2 p ) b b 1 q -(b 1 q + b 2 q ) b + a 2 1 + a 2 -1 4 -2 b (1 -b) (b 1 p + b 2 p ) (b 1 q + b 2 q ) + b 1 p -(b 1 p + b 2 p ) b b 1 q -(b 1 q + b 2 q ) b .
We fix b 2 p = λb 1 p and b 2 q = µb 1 q . The equation becomes:

σ = -b (1 -b) [1 -(1 + λ) a 0 ] [1 -(1 + µ) a 0 ] -a 0 (1 -a 0 ) [1 -(1 + λ) b] [1 -(1 + µ) b] + a 2 1 + a 2 -1 4 [-2 b (1 -b) (1 + λ) (1 + µ) + [1 -(1 + λ) b] [1 -(1 + µ) b]] .
To simplify the notations, we introduce:

A = 1 -a 0 a 0 (donc a 0 = 1 1 + A ), B = 1 -b b , R 2 = a 2 1 + a 2 -1 < a 2 0 .
We obtain the following equivalence, where ∼ denotes the fact that σ is equal, up to positive factors, to the term on the right:

σ ∼ -B b 2 a 2 0 [A -λ] [A -µ] -A a 2 0 b 2 [B -λ] [B -µ] + b 2 R 2 4 [-2 B (1 + λ) (1 + µ) + [B -λ] [B -µ]] .
We rewrite it as:

σ ∼ -[λ -A] [µ -A] + A B [λ -B] [µ -B] + R 2 4Ba 2 0 [λ -B] [µ -B] -2 B (1 + λ) (1 + µ) .
To simplify the computation, we choose A = B and we introduce

α = R 2 /(4Ba 2 0 ) < 1/(4B) = 1/(4A).
Eventually, we obtain:

σ ∼ -(2 -α) [λ -A] [µ -A] -2αA (λ + 1) (µ + 1), σ wrong ∼ -2 [λ -A] [µ -A] .
Hence, the sign of σ wrong changes each time λ = A or µ = A are crossed. Concerning σ, we have σ < 0 when:

λ -A λ + 1 =x µ -A µ + 1 =y > 2αA α -2 if α < 2 < 2αA α -2 if α > 2 without condition if α = 2
We remark that:

-A ≤ x = 1 - 1 + A λ + 1 ≤ 1, -A ≤ y ≤ 1.
To find the conditions for a change of sign for σ, we plot the hyperboles: there is possibility when one branch of the hyperbole is in the square x ∈ [-A, 1], y ∈ [-A, 1]. In Fig. 2.1 (left), the case α < 2 is illustrated : there is a change of sign for σ when α < 2 3 . In Fig. 2.1 (right), the case α > 2 is illustrated : there is a change of sign for σ when α(1 -2A) > 2 or α(A -2) > 2A.

We illustrate this on an example :

A = B = 1, µ = 3, α = 1/10.
Here, we have σ > 0 if and only if λ < 15/23 when σ wrong > 0 if and only if λ < 1. Hence, for 15/23 < λ < 1, we have σ < 0 while σ wrong > 0. Numerical issues. We still work on this example, with the following parameters :

a p = 0.4 0.3 , b p = 0.2 0.1 , a q = 0.1 0.2 , b q = 0.5 0.3 , b = 0.3, a 0 = 0.4, a 1 = 0.3, a -1 = 0.2.
We compute the exact solution using a Runge-Kutta method of order 4 with initial condition p 0 = 0.1 0.2 and q 0 = 0.3 0.4 . We represent on Fig. 2.2 the exact solution and the corresponding phase portrait for ε = 0.1 and on Fig. 2.3 for ε = 0.01. Then, we compute explicitly all the coefficients A 0 p , A 1 p , etc, and we define h 0 and h 1 such that hε = h 0 + εh 1 is an approximation up to order ε 2 of the central manifold. The shape of the functions is given on Fig. 2.4 for x p = x q = 1. Figure 2.2: Direct approximation of the solution p 1 , p 2 , q 1 and q 2 for ε = 0, 1 and the corresponding phase portrait. Figure 2.3: Direct approximation of the solution p 1 , p 2 , q 1 and q 2 for ε = 0, 01 and the corresponding phase portrait. Evolution of h 1 q (x p ,x q ,θ) for x p =x q =1 θ Figure 2.4: The shape of h 0 and h 1 as functions of θ.

Introduction

In the neighborhood of an equilibrium point of a dynamical system, the center manifold is made of orbits which are neither attracted by the stable manifold nor repulsed by the unstable one. A preliminary step when studying equilibria is usually to linearize the system: the phase-space can then be decomposed as the direct sum of the stable and unstable eigenspaces of the linear operator. The former corresponds to eigenvectors associated with eigenvalues having negative real parts whereas the latter is formed by eigenvectors associated with eigenvalues having positive real parts. If the equilibrium under consideration is hyperbolic (i.e. if all eigenvalues have nonzero real parts), the behavior of the dynamical system is fully characterized. If not however, i.e. if some eigenvalues have null real parts, then the corresponding eigenvectors give rise to a center manifold. If these eigenvalues not only have null real parts but are zero, then the center manifold is called a slow manifold. This is the situation we consider in this paper and which appears, for instance, in a number of applications to population dynamics1 [AR94, CHL09]. More precisely, we are concerned with partitioned systems of differential equations of the form

ẋ = εf (x, z) ż = -Λz + εg(x, z) , (3.1.1)
with initial condition (x(0), z(0)) = (x 0 , z 0 ) ∈ R n × R m and where Λ ∈ R m×m is a diagonal matrix with strictly positive diagonal elements λ i , i = 1, . . . , m. In essence, existing theorems from the literature address the possibility that the λ i are different. In this first paper however, we analyze the technically less demanding situation where all λ i coincide. Besides and without additional loss of generality, we fix to 1 their common value.

A statement of the center manifold theorem

The celebrated center-manifold theorem -see for instance [START_REF] Carr | Applications of Centre Manifold Theory[END_REF]-assumes here the following wording:

Theorem 3.1.1 (Center manifold theorem and shadowing principle) Let B R ⊂ R n × R m be the ball of radius R centered at the origin. For all R > 0, there exists ε * > 0 and T > 0 such that the solution (x(t), z(t)) of (3.1.1) exists for all 0 < ε < ε * , all 0 ≤ t ≤ T /ε and all initial condition (x 0 , z 0 ) ∈ B R . Moreover, there exists a (ε-dependent) function h : R n → R m , defined for all 0 < ε < ε * , such that

M = {(x, z) ∈ R n × R m ; z = εh(x)}
is an invariant manifold of (3.1.1) in the following sense: if (x 0 , z 0 ) ∈ M ∩ B R , then (x(t), z(t)) ∈ M for all t ∈ [0, T /ε]. Denoting ϕ t the exact flow of the reduced equation

u = εf (u, εh(u)) (3.1.2)
one can assert that there exists µ > 0 such that for all

(x 0 , z 0 ) ∈ B R ∀t ∈ 0, T ε , z(t) = ε h (ϕ t (x 0 )) + O e -µt .
Furthermore, there exists a modified initial data

x ε 0 ∈ R n such that ∀t ∈ 0, T ε , x(t) = ϕ t (x ε 0 ) + O e -µt and z(t) = εh(ϕ t (x ε 0 )) + O e -µt . (3.1.3)
Finally, if T = +∞ for some R, then equations (3.1.3) are satisfied for all t ≥ 0.

The interest of the center manifold theorem is apparent: provided the function h and the value of x ε 0 are known, the dynamics of (3.1.1) is asymptotically described by the reduced system

ẋ∞ = εf x ∞ , εh(x ∞ ) z ∞ = εh(x ∞ )
with modified initial condition x ∞ (0) = x ε 0 . Besides, it is not hard in our setting to obtain a formal ε-expansion of h from the partial differential equation

ε∂ x h(x) f x, εh(x) = -h(x) + εg x, εh(x) .
(3.1.4)

Scope of the paper

Theorem 3.1.1 states the existence of a function h, and a perturbed initial condition x ε 0 , allowing for a reformulation of the dynamics of (3.1.1) with an asymptotically exponentially small discrepancy. Nevertheless, whereas h may be approached iteratively through (3.1.4), nothing is said on how to construct x ε 0 . Obtaining an expression of x ε 0 in terms of x 0 and z 0 is hence clearly part of our motivation, as it appears as a prerequisite for obtaining an exponentially-close x-approximation. An additional strong motivation for this work stems from the need for a better approximation of the transient phase, that is to say the "small" interval of time close to the initial time where the solution (x(t), z(t)) undergoes a rapid variation. From this point of view, Theorem 3.1.1 is indeed largely unsatisfactory given that e -µt ≈ 1 for small values of t. Our main contribution in this work is to show that it is possible to complement the reduced center-manifold equation in x ∞ with a second equation in ỹ, leading to a system

ẋ∞ = εf x ∞ , εh(x ∞ ) ẏ = εG(x ∞ , ỹ) ,
whose solution for the initial condition (x ∞ (0), ỹ(0)) = (x ε 0 , y 0 ) can then be used to compute the solution (x(t), z(t)) of (3.1.1) exactly for all time. All the transformations required to prove this result are obtained via B-series, introduced as such in [START_REF] Hairer | On the Butcher group and general multi-value methods[END_REF] and pioneered by J.C. Butcher [START_REF] Butcher | An algebraic theory of integration methods[END_REF][START_REF] Butcher | The numerical analysis of ordinary differential equations[END_REF]. B-series are series expansions in powers of ε, which allow for the effective explicit computation of, for instance, the exact solution of (3.1.1). They involve two types of terms: on the one hand, scalar coefficients2 which are universal (in the sense that they are independent of the specific functions f and g) and encode the intrinsic properties of the class of systems being studied, and on the other hand, so-called elementary differentials (composed of various derivatives of f and g and constructed in a very simple way) . This type of representation has proved to be very helpful to construct modified equations [START_REF] Chartier | Algebraic structure of B-series[END_REF] or to analyze highly-oscillatory differential equations [START_REF] Chartier | Higher order averaging, formal series and numerical integration i: B-series[END_REF][START_REF] Chartier | Higher order averaging, formal series and numerical integration ii: the quasi-periodic case[END_REF]. The B-series approach may also be applied to the derivation of estimates for the remainder (the terms O (e -µt ) of Theorem 3.1.1, see for instance [START_REF] Chartier | Higher order averaging, formal series and numerical integration iii: error bounds[END_REF][START_REF] Chartier | A formal series approach to averaging: exponentially small error estimates[END_REF]). This aspect will be the subject of a forthcoming paper and will not be to addressed any longer here. It is worth mentioning that word-series, though less general than B-series, constitute an appealing alternative to B-series, as they are much simpler to compose. In a recent series of papers [START_REF] Murua | Averaging and computing normal forms with word series algorithms[END_REF][START_REF] Murua | Computing normal forms and formal invariants of dynamical systems by means of word series[END_REF][START_REF] Murua | Word series for dynamical systems and their numerical integrators[END_REF], A. Murua and J.M. Sanz-Serna resort to word-series to compute normal forms of a large class of systems including (3.1.1). Their works share many similarities with the present one: in particular [START_REF] Murua | Averaging and computing normal forms with word series algorithms[END_REF] considers the same transport equation as in Section 3.2.4 of this paper, or Section 2.4 of [START_REF] Chartier | Higher order averaging, formal series and numerical integration ii: the quasi-periodic case[END_REF]. However, it differs in that the adjoined initial conditions are not prescribed at the same time. This choice manifests itself in the resulting normal equations (and will be discussed in greater details in Remark 3.2.26): the form obtained here in Theorem 3.2.25 and sketched above is close to the standard one of Theorem 3.1.1 and retains its main advantage, i.e. the reduction of dimensionality: starting from a problem posed in R n × R m , the normal form obtained here is a partially decoupled system.

The main ideas of the paper are exposed in Section 3.2. At first, we shall motivate in Subsection 3.2.1 the use of B-series, with a direct attempt at deriving the formal expansions of the solution of (3.1.1). We will then present, in Subsection 3.2.2, the trees and elementary differentials required in this context and introduce Taylor-indexed B-series in Subsection 3.2.3, together with some of their features which are essential to the subsequent analysis. Subsection 3.2.4 is devoted to the main properties of the B-series coefficients of the exact solution of (3.1.1). In particular, it is shown therein that they are amenable to a time-scale separation after which they obey a transport equation. The transformation that maps x 0 to x ε 0 is also explicitly defined in this part. Equations of the centermanifold dynamics are then derived in Subsection 3.2.5, where a theorem similar to Theorem 3.1.1 is given. We conclude this second section with the statement of the main result of the paper in Subsection 3.2.6.

In Section 3.4, we illustrate numerically the main outcomes of our analysis on two simple examples. For both systems, we derive a third-order approximation of all transformations considered and show that the results of our theorems are indeed valid up to errors of fourth-order with respect to (w.r.t.) ε.

Center manifold via B-series

In [START_REF] Chartier | Higher order averaging, formal series and numerical integration i: B-series[END_REF] and [START_REF] Chartier | Higher order averaging, formal series and numerical integration ii: the quasi-periodic case[END_REF], where highly-oscillatory systems are studied, the analysis leans on a Fourier expansion of the periodic vector field. Here, the corresponding ad-hoc tool is a Taylor expansion in the neighbourhood of z = 0. We hereafter explain how to proceed.

Through the change of variable z(t) = e -tΛ y(t), system (3.1.1) can be equivalently rewritten as ẋ = εf x, e -tΛ y := εF tΛ (x, y) ẏ = εe tΛ g(x, e -tΛ y) := εG tΛ (x, y) , (3.2.1) with initial condition (x(0), y(0)) = (x 0 , z 0 ). Assuming that both f and g are real-analytic w.r.t. the y-variable, we have the Taylor expansions

f (x, z) = f (x, 0) + k∈N m 1 k! D k z f (x, 0) z k , (3.2.2) g(x, z) = g(x, 0) + k∈N m 1 k! D k z g(x, 0) z k , (3.2.3)
where

k! = k 1 ! • • • k n ! and where D k z f (x, 0) z k i = ∂ |k| f i (x, 0) ∂z k1 1 • • • ∂z kn n z k1 1 . . . z kn n with |k| = k 1 + . . . + k n
and similarly for g. In particular, one has

D k z f (x, 0) (e -tΛ y) k = e -t(k•λ) D k z f (x, 0) y k with k • λ = n i=1 k i λ i ,
so that one can eventually write

F tΛ (x, y) = f (x, 0) + k∈N m e -t(k•λ) k! D k z f (x, 0) y k , G tΛ (x, y) = e tΛ g(x, 0) + k∈N m e -t(k•λ) k! e tΛ D k z g(x, 0) y k ,
expressions which can be further simplified by taking into account that all λ i are assumed to be equal to 1

F t (x, y) = f (x, 0) + k∈N m e -t|k| k! D k z f (x, 0) y k , G t (x, y) = e t g(x, 0) + k∈N m e -t(|k|-1) k! D k z g(x, 0) y k .
To sum up, the equations analyzed throughout this paper are of the following form

ẋ(t) = = εF t (x, y) := ε ∞ k=0 e -kt f k (x, y), x(0) = x 0 ẏ(t) = = εG t (x, y) := ε ∞ k=0 e -(k-1)t g k (x, y), z(0) = z 0 (3.2.4)
where we have denoted3

f k (x, z) = k∈N m , |k|=k 1 k! (D k z f )(x, 0) z k and g k (x, z) = k∈N m , |k|=k 1 k! (D k z g)(x, 0) z k .
Since we shall not study the convergence of the series manipulated in this paper4 , we will furthermore assume that the series in (3.2.4) are finite, i.e. that f and g are polynomials.

Expansion of the transient solution

In order to motivate the introduction of trees and B-series in next subsection, we first derive a few terms of the formal ε-expansion of the components (x(t), y(t)) of the solution and this is done by considering equations (3.2.4) in their integral form:

x(t) = x 0 + ε t 0 k≥0
e -ks f k (x(s), y(s))ds = x 0 + O(ε),

y(t) = y 0 + ε t 0 k≥0
e -(k-1)s g k (x(s), y(s))ds = y 0 + O(ε).

Introducing these expressions in the right-hand side of the equations, we now obtain

x(t) = x 0 + ε k≥0 t 0 e -ks ds f k (x 0 , y 0 ) + O(ε 2 ), y(t) = y 0 + ε k≥0 t 0 e -(k-1)s ds g k (x 0 , y 0 ) + O(ε 2 ).
Omitting the argument (x 0 , y 0 ) of the various functions, a third iteration then leads to It is clear that this procedure à la Picard can be iterated to obtain the ε 3 , ε 4 , ... terms of the expansions of x(t) and y(t). However, the growing complexity of the expressions arising in the process impedes a systematic construction. This is the reason why we shall use Taylor-indexed rooted trees and associated elementary differentials as a mean to derive the sought-after series with explicit inductions.

x(t) = x 0 + ε k≥0 t 0 e -ks ds f k + ε 2 k,r t 0 e -ks s 0 e -rσ dσds (∂ x f k )f r

Taylor-indexed bicoloured trees and elementary differentials

We consider bi-coloured rooted trees where black vertices refer to function f , white vertices to function g and where each vertex has been labelled with a index k ≥ 0 associated to the k-th terms f k and g k in the Taylor expansions of f and g. For the sake of simplicity, we use hereafter the word tree.

Definition 3.2.1 The set of Taylor-index bicoloured trees (or simply trees) T = T • ∪ T • is defined recursively as follows:

1. For any index k ∈ N, the tree with a single indexed vertex • k belongs to T • and the tree with a single indexed vertex • k belongs to T • .

2. For any index k ∈ N, any (p, q) ∈ N × N with q ≤ k, any (u 1 , . . . , u p ) ∈ T p • and (v 1 , . . . , v q ) ∈ T q • , the tree [u 1 , . . . , u p , v 1 , . . . , v q ] • k obtained by connecting the roots of u 1 , . . . , u p , v 1 , . . . , v q to a new root

• k , belongs to T • . Similarly, [u 1 , . . . , u p , v 1 , . . . , v q ] • k ∈ T • .
For homogeneity, the empty trees ∅ x and ∅ y will sometimes be used to denote F ∅x (x, y) = x and F ∅y (x, y) = y, and we shall write accordingly

T• = T • ∪ {∅ x }, T• = T • ∪ {∅ y } and T = T• ∪ T• .
Trees with a number of branches q in T • strictly greater than the index k carried by the root are not permitted. This comes from the fact that the q-th derivatives of f k and g k with respect to y vanish identically for q > k, given that both f k and g k are k-linear maps with respect to y. For instance, 1 1 ∈ T • but 0 1 / ∈ T • . This will become completely clear with Definition 3.2.3.

The order of a tree w ∈ T , denoted |w|, is its number of vertices. The symmetry factor measures how symmetric the tree looks like and is defined as follows: note that vertices with different labels are distinguished. Definition 3.2.2 The symmetry factor is defined recursively on T as follows:

1. For all k ∈ N, σ • k = σ • k = 1.

Let w be of the form either

[u µ1 1 , . . . , u µp p , v ν1 1 , . . . , v νq q ] • k or [u µ1 1 , . . . , u µp p , v ν1 1 , . . . , v νq q ]
• k , where trees u i and v j are assumed to be pairwise distinct, and where exponents µ i and ν j indicate that u i and v j are repeated µ i and ν j times. Then

σ w = p i=1 µ i !σ µi ui q j=1 ν j !σ νj vj .
Finally, to each tree we associate an elementary differential, i.e. a function from R n × R m to either R n or R m , depending on whether the root is black or white. The label of the root then determines which function is differentiated.

Definition 3.2.3 The elementary differentials associated to trees of T are defined recursively as follows:

1. For all k ∈ N,

F • k (x, y) = f k (x, y) and F • k (x, y) = g k (x, y). 2. If u = [u 1 , . . . , u p , v 1 , . . . , v q ] • k , then F u (x, y) = (∂ p x ∂ q y f k )(x, y) F u1 (x, y), .
. . , F up (x, y), F v1 (x, y), . . . , F vq (x, y) ,

and if v = [u 1 , . . . , u p , v 1 , . . . , v q ] • k , then F v (x, y) = (∂ p x ∂ q y g k )(
x, y) F u1 (x, y), . . . , F up (x, y), F v1 (x, y), . . . , F vq (x, y) .

According to previous definitions, the truncated expansions (3.2.5) and (3.2.6) can be rewritten as

x(t) = x 0 + u∈T•,|u|≤2 ε |u| α u (t) σ u F u (x 0 , y 0 ) + O(ε 3 ), y(t) = y 0 + v∈T•,|v|≤2 ε |v| α v (t) σ v F v (x 0 , y 0 ) + O(ε 3 ), with α • k (t) =

Taylor-indexed partitioned B-series

In this subsection, we now consider Taylor-indexed partitioned B-series, which will constitute the main tool employed in this paper. For brevity again, we shall simply call them B-series.

Definition 3.2.4 A Taylor-indexed partitionned B-series (or simply B-series) with coefficients a : T → C, is a formal expansion of the form

B(a, (x, y)) = a ∅x x + u∈T• ε |u| σ u a u F u (x, y), a ∅y y + v∈T• ε |v| σ v a v F v (x, y) .
We shall incidentally write

B • (a, (x, y)) = a ∅x x + u∈T• ε |u| σ u a u F u (x, y)
for the component x of B(a, (x, y)) and accordingly B • (a, (x, y)) for the component y. Although this is not reflected in the notations, a B-series depends on the functions f k and g k , and thus on f and g, through the elementary differentials F . According to this definition and to equations (3.2.4), we can write the function (x, y) → (εF t (x, y), εG t (x, y)) as a B-series

(εF t (x, y), εG t (x, y)) = B(β(t), (x, y)) (3.2.7)
with coefficients β depending on t and defined as follows:

β • k = e -kt , β • k = e -(k-1)t , β w = 0 for all w ∈ T \{• k , • k , k ∈ N}.
Two B-series with coefficients a and b such that a ∅x = a ∅y = 1 can be composed to form a new B-series with coefficients c, that is to say B(b, B(a, (x, y)) = B(c, (x, y)) with c = a * b.

More precisely, each c w is an explicitly known polynomial of the a w ′ and b w ′′ . The star product is non-commutative and the set of near-identity mappings a ∈ C T with a ∅x = a ∅y = 1 is a non-commutative group G, named Butcher group, with unit element 1 1, defined by 1 1 ∅x = 1 1 ∅y = 1 and 1 1 w = 0 for all w ∈ T . In particular, every element a ∈ G has an inverse a -1 ∈ G such that a * a -1 = 1 1 = a -1 * a. Note that more generally, the star product a * b is well defined for a ∈ G and b ∈ C T . For the sake of illustration, we give the terms of the star-product for some trees of order less than, or equal to, 3, where j, k, l are three positive integers:

c ∅x = b ∅x c ∅y = b ∅y c •j = b ∅x a •j + b •j c •j = b ∅y a •j + b •j c j k = b ∅x a j k + b •j a • k + b j k c j k = b ∅y a j k + b •j a • k + b j k c j k = b ∅x a j k + b •j a • k + b j k c j k = b ∅y a j k + b •j a • k + b j k
It is apparent that the color of vertices does not play a specific role, so for the trees of order 3, we content ourselves in this brief exposition with the following mono-coloured trees: To conclude this subsection, we now re-derive the expansion of the transient solution (x(t), y(t)) by using the star-product. Denoting (x(t), y(t)) = B(α(t), (x 0 , y 0 )), the differential equations (3.2.1) can be rewritten in terms of B-series as d dt B(α(t), (x 0 , y 0 )) = B(α(t) * β(t), (x 0 , y 0 )) B(α(0), (x 0 , y 0 )) = B(1 1, (x 0 , y 0 ))

c j k l = b ∅x a j k l + b •j a • k a • l + b j k a • l + b j l a • k + b j k l c j k l = b ∅x a j k l + b •j a j k + b j k a • l + b j l a • k + b j k l A
Of course, this may be translated as a Cauchy problem in terms of the coefficients α ∈ G R as

d dt α(t) = α(t) * β(t), α(0) = 1 1, (3.2.9)
where β is defined in (3.2.7). Note that since β vanishes for trees of orders greater than 1 and owing to the expression of the star-product, we hereby obtain for all k ∈ N

d dt α • k (t) = β • k (t), d dt α • k (t) = β • k (t),
and for all u = [u 1 , . . . , u p , v 1 , . .

. , v q ] • k ∈ T • and v = [u 1 , . . . , u p , v 1 , . . . , v q ] • k ∈ T • , d dt α u (t) = β • k (t) p i=1 α ui (t) q i=1 α vi (t), d dt α v (t) = β • k (t) p i=1 α ui (t) q i=1
α vi (t), (3.2.10) which, together with the initial conditions α w (0) = 0, give the formulae of Subsection 3.2.1. Since z 0 = y 0 , we can also write

x(t) = x 0 + u∈T• ε |u| σ u α u (t)F u (x 0 , z 0 ), (3.2.11) z(t) = e -t z 0 + v∈T• ε |v| σ v (e -t α v (t))F v (x 0 , z 0 ). (3.2.12)
The first coefficients α(t) for trees of orders less than or equal to 2 are given in Tables 3.1, 3.2.

u • 0 • k 0 0 0 r r ≥ 1 α u t 1-e -kt k t 2 2 -1+rt+e -rt r 2 u k r 1 0 k 1 k r r ≥ 1 (k, r) = (1, 0), r = 1 α u r-(r+k)e -kt +ke -(r+k)t rk(r+k) -1 + t + e -t
1-e -kt (1+kt) k 2 r-1-(k+r-1)e -kt +ke -(k+r-1)t (r-1)k(k+r-1)

≥ 1. v • 1 • k 1 0 1 r k r 0 1 r = 0 r = 0 α v t 1-e -k ′ t k ′ t 2 2 -1+rt+e -rt r 2 r-(k ′ +r)e -k ′ t +k ′ e -(k ′ +r)t rk ′ (k ′ +r) -1 + e t -t v 1 1 1 r 2 0 k 1 k r r = 1 k ≥ 2 r = 1, k ≥ 2, k + r = 2 α v t 2 2 -1+r ′ t+e -r ′ t (r ′ ) 2 t -1 + e -t 1-(1+k ′ t)e -k ′ t (k ′ ) 2 r ′ -(k ′ +r ′ )e -k ′ t +k ′ e -(k ′ +r ′ )t r ′ k ′ (k ′ +r ′ ) Table 3.2: Coefficients α for trees v of T • with |v| ≤ 2, k = 1, k ′ = k -1 and r ′ = r -1.

The transport equation

This subsection contains all the technical results used to state and prove the main results of the paper. In the spirit of Lemma 2.4 and Definition 2.5 of [START_REF] Chartier | Higher order averaging, formal series and numerical integration ii: the quasi-periodic case[END_REF], we define polynomial functions by separating the slow-time t and fast-time τ variables. This requires here a little bit more care than in [START_REF] Chartier | Higher order averaging, formal series and numerical integration ii: the quasi-periodic case[END_REF], since we wish to keep track of the fact that coefficients α(t) involve exponential terms of the form e -kt with k ≥ -1 only. Whereas in quasistroboscopic averaging as in [START_REF] Chartier | Higher order averaging, formal series and numerical integration ii: the quasi-periodic case[END_REF], the value at τ = 0 of coefficients γ(t, τ ) was playing a major role, this is the value at τ = +∞ which here becomes central to the analysis. where α is the solution of equation (3.2.9).

Proof: The proof proceeds by induction on the order of trees. For order 1, the assertion of the proposition can be straightforwardly checked. Now, consider u α vi (s)ds.

= [u 1 , . . . , u p , v 1 , . . . , v q ] • k ∈ T • and v = [u 1 , . . . , u p , v 1 , . . . , v q ] • k ∈ T • ,
Using the induction hypothesis and taking into account that the set of polynomial functions is an algebra which is stable by derivation and integration (w.r.t. both t and τ ), since q ≤ k the function

e -ks p i=1 α ui (s) q i=1 α vi (s) = e -(k-q) p i=1 α ui (s) q i=1 (e -s α vi (s))
is also of the form w(t, τ )| τ =t for some polynomial function w, and so is α u (t). A similar conclusion holds for α v (t) if q ≤ k -1. Now if q = k, we write w(t, τ ) = P (t) + e -τ w(t, τ ) where P (t) = lim τ →+∞ w(t, τ ) is polynomial in t and w another polynomial function. Denoting Q the unique polynomial such that Q ′ + Q = P , we have

α v (t) = t 0
(e s P (s) + w(s, s))ds = e t Q(t) + t 0 w(s, s)ds so that e -t α v (t) is again a polynomial function.

For instance, according to Table 3.1 and Table 3.2, we have

α k 0 (t) = 1 -e -kt k 2 - te -kt k and α 1 0 (t) = -t + e t -1,
so that by substituting t by τ in exponential terms, we obtain

γ k 0 (t, τ ) = 1 -e -kτ k 2 - te -kτ k and γ 1 0 (t) = -1 -t + e τ .
The values of γ for trees or orders less than or equal to 2 are given in tables 3.3 and 3.4.

u • 0 • k 0 0 0 r r ≥ 1 γ u t 1-e -kτ k t 2 2 -1+rt+e -rτ r 2 u k r 1 0 k 1 k r r ≥ 1 (k, r) = (1, 0), r = 1 γ u r-(r+k)e -kτ +ke -(r+k)τ rk(r+k) -1 + t + e -τ
1-e -kτ (1+kt) k 2 r-1-(k+r-1)e -kτ +ke -(k+r-1)τ (r-1)k(k+r-1)

Table 3.3: Coefficients γ for trees u ∈ T • with |u| ≤ 2 and k ≥ 1. v • 1 • k 1 0 1 r k r 0 1 r = 0 r = 0 γ v t 1-e -k ′ τ k ′ t 2 2 -1+rt+e -rτ r 2 r-(k ′ +r)e -k ′ τ +k ′ e -(k ′ +r)τ rk ′ (k ′ +r) -1 + e τ -t v 1 1 1 r 2 0 k 1 k r r = 1 k ≥ 2 r = 1, k ≥ 2, k + r = 2 γ v t 2 2 -1+r ′ t+e -r ′ τ (r ′ ) 2 t -1 + e -τ 1-(1+k ′ t)e -k ′ τ (k ′ ) 2 r ′ -(k ′ +r ′ )e -k ′ τ +k ′ e -(k ′ +r ′ )τ r ′ k ′ (k ′ +r ′ ) Table 3.4: Coefficients γ for trees v of T • with |v| ≤ 2, k = 1, k ′ = k -1 and r ′ = r -1. Proposition 3.2.8 If γ ∈ G R×R is defined as in (3.2.13), then γ(0, 0) = 1 1 and ∀(t, τ ) ∈ R × R, ∂ t γ(t, τ ) + ∂ τ γ(t, τ ) = γ(t, τ ) * β(τ ). ( 3
, 0) = α(0) = 1 1.
In contrast with the general situation where equation (3.2.14) may have infinitely many solutions with the mere initial condition γ(0, 0) = 1 1, the polynomial nature of γ ensures here that there is only one. Actually, uniqueness can be ensured by prescribing the value of γ at (0, 0), as in [START_REF] Chartier | Higher order averaging, formal series and numerical integration ii: the quasi-periodic case[END_REF], or at (0, τ 0 ), as in [START_REF] Murua | Averaging and computing normal forms with word series algorithms[END_REF]. Since the asymptotic dynamics of B • (γ(t, τ ), (x 0 , y 0 )) is here attained for τ = +∞, we have to address the question of uniqueness of the solution of (3.2.14) in the two situations of Lemma 3.2.9 and Lemma 3.2.10. Lemma 3.2.9 Given a polynomial function w : R × R → R, c 1 ∈ R and c 2 ∈ R, there exists a unique polynomial solution of For k = 1, . . . , j, there exists a unique polynomial solution of ϕ ′ kkϕ k = W k , while for k > j, ϕ ′ kkϕ k = 0 implies ϕ k = 0. As for k = 0, we get ϕ 0 (t) = t 0 W 0 (s)ds + C where C is then uniquely defined by ϕ(0, 0) = n k=1 ϕ k (0) + C = c 1 or by ϕ(0, +∞) = C = c 2 . Lemma 3.2.10 Given a polynomial function w : R × R → R and c ∈ R, there exists a unique solution ϕ(t, τ ) = e τ ψ(t, τ ), with ψ a polynomial function, of :

∂ t ϕ(t, τ ) + ∂ τ ϕ(t, τ ) = w(t, τ ), ϕ(0, 0) = c 1 or ϕ(0, +∞) = c 2 .
∂ t ϕ(t, τ ) + ∂ τ ϕ(t, τ ) = e τ w(t, τ ), ϕ(0, 0) = c Proof: Writing w(t, τ ) = j k=0 e -kτ W k (t) and ϕ(t, τ ) = e τ r k=0 e -kτ ϕ k (t), the differential equation becomes r k=0 e (1-k)τ (ϕ ′ k (t) + (1 -k)ϕ k (t)) = j k=0 e (1-k)τ W k (t).
For k = 0, 2, . . . , j, there exists a unique polynomial solution of ϕ

′ k + (1 -k)ϕ k = W k , while for k > j, ϕ ′ k -kϕ k = 0 implies ϕ k = 0. As for k = 1, we get ϕ 1 (t) = t 0 W 1 (s)ds + C, where C is then uniquely defined by ϕ(0, 0) = ϕ 0 (0) + C + n k=2 ϕ k (0) = c. Proposition 3.2

.11

1. There exists a unique polynomial map γ ∈ G R×R solution of (3.2.14) satisfying γ(0, 0) = 1 1.

There exists a unique polynomial map

δ ∈ G R×R solution of ∂ t δ(t, τ ) + ∂ τ δ(t, τ ) = δ(t, τ ) * β(τ ) (3.2.15) satisfying both conditions δ| T• (0, +∞) = 1 1| T• and δ| T• (0, 0) = 1 1| T• .
Remark 3.2.16 Trees with zero norm are distinctive for the corresponding elementary differentials do not depend on y. This stems from the very definition of f k (x, y) and g k (x, y) as k-linear maps w.r.t. the variable y, and from the fact that, in a tree with zero norm, all functions f k or g k appearing are differentiated exactly k times. It is also worth mentioning that in a tree of zero norm, the index carried by any vertex is exactly equal to the number of upcoming • branches. In this sense, the set of zero norm indexed partitioned trees is isomorphic to the set of partitioned trees with no label. For the sake of illustration, consider for instance the zero-norm tree u

= [• 0 , • 0 , • 0 , • 0 , • 0 ] •3 . Since ∂ 3 y f 3 (x, y) = 1 3! ∂ 3 y (∂ 3 y f )(x, 0) y 3 = (∂ 3 y f )(x, 0), we have F u (x, y) = (∂ 2 x ∂ 3 y f )(x, 0) F •0 (x, y), F •0 (x, y), F •0 (x, y), F •0 (x, y), F •0 (x, y) = (∂ 2 x ∂ 3 y f )(x, 0) f (x, 0), f (x, 0), g(x, 0), g(x, 0), g(x, 0) . Proposition 3.2.17 For all u = [u 1 , . . . , u p , v 1 , . . . , v q ] • k ∈ T • and all v = [u 1 , . . . , u p , v 1 , . . . , v q ] • k ∈ T • such that q = k, the following relations δ∞ u (t) = p i=1 δ ∞ ui (t) q j=1 δ ∞ vj (t), δ ∞ u (0) = 0 and δ∞ v (t) + δ ∞ v (t) = p i=1 δ ∞ ui (t) q j=1 δ ∞ vj (t) (3.2.24) are satisfied. Furthermore, if w ∈ T is such that w > 0, then δ ∞ w (t) = 0 for all t ∈ R.
Proof: By definition of δ, we have

(∂ t + ∂ τ )δ(t, τ ) = δ(t, τ ) * β(τ ) with δ ∞ | T• (0) = δ| T• (0, +∞) = 1 1| T• . Hence, for all u = [u 1 , . . . , u p , v 1 , . . . , v q ] • k ∈ T • , one has e -(k-q)τ p i=1 δ ∞ ui (t) q j=1 δ ∞ vj (t) + e -τ w(t, τ ) = δ∞ u (t) + e -τ w(t, τ ) (3.2.25)
where w and w are polynomial functions. Similarly, for all v = [u 1 , . . . , u p , v 1 , . .

. , v q ] • k ∈ T • , one has e -(k-1-q)τ p i=1 δ ∞ ui (t) q j=1 δ ∞ vj (t) + w(t, τ ) = e τ ( δ∞ v (t) + δ ∞ v (t)) + w(t, τ ). ( 3

.2.26)

For k = q, this clearly proves equations (3.2.24). Now, suppose that u > 0 and v > 0 respectively: then either k < q or k = q and at least one amongst the branches u i or v j has a strictly positive norm. In the first case, we obtain respectively δ∞ u (t) = 0 and δ∞ v (t) + δ ∞ v (t) = 0 (according to (3.2.25) and (3.2.26)), so that δ ∞ u (t) = 0 owing to the initial condition δ ∞ u (0) = 0, and δ ∞ v (t) = 0 owing to the fact that δ ∞ v (t) is polynomial in t. In the second case, relation (3.2.24) is satisfied with a right-hand side that vanishes owing to an induction argument, and we can conclude similarly. Remark 3.2.18 Note that the initial condition δ ∞ | T• ( 0) is not known a priori. It is actually determined by solving equation (3.2.24). This is accordance with the fact that δ ∞ are the B-series coefficients of the asymptotic dynamics of (x(t), z(t)), the solution of (3.2.1).

Relation (3.2.16) can be rewritten as

δ(t + t ′ , τ ) = γ-1 (0) * γ(t ′ ) * δ(t, τ ). Furthermore, it is clear that δ ∞ | T• (t) = δ| T• (t, +∞) = γ-1 (0) * γ| T• (t, ∞) = γ-1 (0) * γ| T• (t), and that δ ∞ | T• (t) = lim τ →∞ e -τ δ| T• (t, τ ).
Passing to the appropriate limits, we thus obtain

δ ∞ (t + t ′ ) = γ-1 (0) * γ(t ′ ) * δ ∞ (t).
Given now that δ ∞ w ≡ 0 if w > 0, the left term in the star product of the right-hand side is evaluated only for trees of T • , and since

(γ -1 (0) * γ(t ′ )) | T• = γ-1 (0) * γ | T• (t ′ ) = δ ∞ | T• (t ′ ),
we may write

δ ∞ (t + t ′ ) = δ ∞ (t ′ ) * δ ∞ (t), (3.2.27) or δ ∞ (t + t ′ ) = δ ∞ | T• (t ′ ) * δ ∞ (t), (3.2.28)
at the price of a slight abuse of notations (we have indeed here δ ∞ ∅y (t ′ ) = 0, in contrast with the definition of the * -product in Subsection 3.2.3). The perfectly rigorous way to write (3.2.27) would be

δ ∞ (t + t ′ ) = (1 1 z + δ ∞ (t ′ )) * δ ∞ (t)
where 1 1 z w = 0 for all w / ∈ T• ∪ T • and 1 1 z ∅y = 1, where the addition of the term 1 1 z is harmless and only technical, since it does not appear in the effective computations of the star product.

Remark 3.2.19 Note that the B-series B(δ ∞ (t), (x, y)) is of the form (x + O(ε), O(ε)), since for all t ≥ 0, δ ∞ ∅z (t) = 0. Besides, since δ ∞ w ≡ 0 for w > 0, that it does not depend on z, that is to say

B(δ ∞ (t), (x, y)) = B(δ ∞ (t), (x, 0)).
As a consequence, the asymptotic behaviour depends on z 0 = y 0 only through the modified intial condition x ε 0 .

In the sequel, we shall denote

(x ∞ (t), z ∞ (t)) = B(δ ∞ (t), (x ε 0 , 0)), (3.2.29)
where (x ε 0 , y 0 ) = Φ -1 0 (x 0 , y 0 ). It may be seen as the trajectory

(x(t, τ ), z(t, τ )) = B • (δ(t, τ ), (x ε 0 , z 0 )), B • (e -τ δ(t, τ ), (x ε 0 , z 0 ))
in the limit τ = +∞, i.e. as the shadow solution of (x(t), y(t)) on the center manifold. We are now in position to identify the center-manifold function εh and the vector field εF of the slow dynamics in the variable x.

Proposition 3.2.20 Consider the two functions h and F defined by their B-series expansions εh(x) := B • (δ ∞ (0), (x, 0)) and εF (x) := B • ( δ∞ (0), (x, 0)).

The following relations hold true εf (x, εh(x)) = B • δ∞ (0), (x, 0) and εg(x, εh(x))) = B • δ∞ (0) + δ ∞ (0), (x, 0) .

(3.2.30)

In particular, one has F (x) = f (x, εh(x)).

Proof: By definition of β, we have εf (x, z), εg(x, z) = B(β(0), (x, z)). Besides, we can write

(x, εh(x)) = B(δ ∞ (0), (x, 0)) = B(1 1 z + δ ∞ (0), (x, 0)),
so that, using the star-product, we get

(εf (x, εh(x)), εg(x, εh(x))) = B (1 1 z + δ ∞ (0)) * β(0), (x, 0) .
To sum up, we wish to prove that

B • (1 1 z + δ ∞ (0)) * β(0), (x, 0) = B • ( δ∞ (0), (x, 0)), B • (1 1 z + δ ∞ (0)) * β(0), (x, 0) = B • ( δ∞ (0) + δ ∞ (0), (x, 0)).
Now, given that F w (x, z)| z=0 = 0 whenever w > 0, all we have to prove is that for all w ∈ T with w = 0, one has

(1 1 z + δ ∞ (0)) * β(0) w = δ∞ w (0) if w ∈ T • , (1 1 z + δ ∞ (0)) * β(0) w = δ∞ w (0) + δ ∞ w (0) if w ∈ T •
Consider a tree w ∈ T either of the form [u 1 , . . . , u p , v 1 , . . . , v q ] •q or [u 1 , . . . , u p , v 1 , . . . , v q ] •q . According to Proposition 3.2.17 on the one hand, and by definition of the * -product on the other hand, we indeed have

δ∞ w (0) = p i=1 δ ∞ ui (0) q j=1 δ ∞ vj (0) = [(1 1 z + δ ∞ (0)) * β(0)] w if w ∈ T • , δ∞ w (0) + δ ∞ w (0) = p i=1 δ ∞ ui (0) q j=1 δ ∞ vj (0) = [(1 1 z + δ ∞ (0)) * β(0)] w if w ∈ T • ,
where we used the explicit definition of β (as in (3.2.10)).

Theorem 3.2.21 The solution (x ∞ (t), z ∞ (t)) defined in (3.2.29) satisfies the following center-manifold system

ẋ∞ (t) = εF (x ∞ (t)) = εf x ∞ (t), εh(x ∞ (t)) z ∞ (t) = εh(x ∞ (t)) , (3.2.31) with initial condition (x ∞ (0), z ∞ (0)) = (x ε 0 , εh(x ε 0 )).
Proof: For trees of T • , relation (3.2.28) becomes

δ ∞ | T• (t + t ′ ) = δ ∞ | T• (t) * δ ∞ | T• (t ′ ),
so that, by choosing t ′ = 0, one has

δ ∞ | T• (t) = δ ∞ | T• (t) * δ ∞ | T• (0), which
(t) = δ ∞ | T• (t) * δ∞ | T• (0). Example 3.2.22
The first terms of the expansion of εh(x) are computed in Table 3.8, where, for conciseness, we have omitted the argument (x, 0) of all functions and used obvious simplified notations for derivatives. Collecting

• 0 0 0 1 0 0 0 0 1 0 0 2 0 0 0 0 0 1 0 0 1 1 0 0 1 0 1 -1 1 2 -1 1 1 -2 1 -1 g g x f g z g g xx (f, f ) g xz (f, g) g zz (g, g) g x f x f g z g x f g z g z g g x f z g Table 3.8: Coefficients δ ∞ (0) for trees u of T • with |u| ≤ 3.
Notation 3.3.24 In Maple, a forest w 1 . . . w n will be denoted [w 1 , . . . , w n ].

To construct B -in Maple, the difficult part lies in the manipulation of indexes. Indeed, we have a direct access to the branches of the tree, but we need to associate each index to the right branch. For example, with the tree In fact, the definition of the index vector is well suited : we read the list of branches and the list of indexes (starting with the second index) and we associate to a branch the number of indexes corresponding to its order. We define B + as a function of a forest, a type of root and an index of root. The only point is to find the rank of the new tree in the global order. For that, we compare it with the list of ordered trees. Now, the coproduct ∆ can be defined. It is done recursively, as in (3.3.4). The map ∆ gives an element of H ⊗ H, with on each side a sum of forests. The point here is to take into account the sum of the elements forests⊗ forests. As ∆ is defined recursively, we have to multiply such sums of forest between them, a multiplication of forests being the concatenation of the forests. We create a function multiplication to do that. Hence, we have all the elements to compute the coproduct.

The composition of B-series α * β is computed using (3.3.3) and we create a function compose(α,β,arbre), with time coefficients if α and β are t or (t, τ ) dependent.

Compute the inverse of a B-series

The only difficulty that remains is to compute the inverse of a B-series. We start by rewriting the formula of the composition of two B-series (0.2.10) and the formula of the inverse of a B-series (3.3.2) for w ∈ T : The formula to compute the inverse looks like the composition formula for α * α, with the addition of (-1) #P .

(α * β) w = β u + p∈P(u)
Example 3.3.27 For the tree :

(α * β) = β + 2β α • + β • (α • ) 2 + α β ∅ α -1 = -α + 2α α • -α • (α • ) 2
Looking at the example, we realise that we have to take care of something else : in the inverse, the partitions , ∅ and ∅, are the same. So we have to check if the partition was not already used. We adapt the definition of the composition to add those two points and the function inverse is defined.

with initial conditions x 0 = 0.8 and z 0 = 0.05. Its exact solution6 for ε = 0.01 is drawn on the left of Figure 3.3. Functions f and g being polynoms of respective degrees 3 and 2 in z, all functions f k and g k vanish identically, except the following ones:

f 1 (x, z) = ∂ z f (x, 0)z = -x 3 z, f 3 (x, z) = 1 6 ∂ 3 z f (x, 0)z 3 = 1 3 x 3 z 3 , g 0 (x, z) = g(x, 0) = x, g 2 (x, z) = 1 2 ∂ 2 z g(x, 0)z 2 = - 1 2 xz 2 .
As a result, only ten trees in T • of orders less than or equal to 2 have to be considered. Their corresponding coefficients γ(0) can be read from Table 3.7 and are listed in Table 3.12 together with their associated elementary differentials. The second-order truncated expansion of x ε 0 , defined by (x ε 0 , y 0 ) = Φ -1 0 (x 0 , y 0 ), is thus of the form 1/2 1/18 -1 1/2 1/6 1/12 F u (x, y) -x 3 y 1 3 x 3 y 3 -x 5 y 4 -x 5 y 4 3x 5 y 2 1 3 x 5 y 6 x 4 1 2 x 4 y 2 x 4 y 2 -1 2 x 4 y 4 Table 3.12: Coefficients γ(0) and elementary differentials for trees of T • with order ≤ 2.

x 0 + εX 1 + ε 2 X 2 , with Getting the third-order term requires a few more calculations that we do not reproduce here: Much fewer terms are present in the expansion of εh(x), as coefficients δ ∞ (0) vanish for all trees with a non-zero norm. Truncating once again at order 3, we indeed obtain εh [3] (x) := εx + ε 3 -1 2 x 3 + x 4 = ε h(x) + O(ε 4 ). ∞ (t) (on the left), z(t) and ε h(x ∞ (t)) (on the right). Finally, the differences x(t)x ∞ , εh [3] (x

[3] ∞ ) ẏ[3] = ε G [3] x [3] ∞ , ỹ[3]
(3.4.3) with initial conditions x ∞ (0) = x ε 0 , ỹ(0) = z 0 , z[3] (t) = e -t ỹ[3] (t), and where εG [3] denotes the third-order truncation of εG. Note that, in order to write down εG [2] , only eight trees of orders less than or equal to 2 have to be considered and their corresponding coefficients δ(0) and elementary differentials computed: they are listed in Table 3.13. Getting the terms of order 3 is totally straightforward, though definitely more painful:

v ∈ T • • 0 • 2 0 1 0 3 2 1 2 3 2 0 2 2
δv (0) 0 0 -1 0 0 0 1 0 F v (x, y) x 1 2 xy 2 -x 3 y 1 3 x 3 y 3 1 2 x 3 y 3 -1 6 x 5 y 5 -x 2 y 1 2 x 2 y 3 Table 3.13: Coefficients δ(0) and elementary differentials for trees of T • with order ≤ 2.

εG [3] (x, y) = ε 2 (x 3 yx 2 y) + ε 3 2x 3 -1 2

x 3 y 2 -2x 4 = εG(x, y) + O ε 4 .

According to Theorem 3.2.25, it remains to approximate the change of variables Φ t up to order 3 in ε, a task that necessitates the values of coefficients δ(0, t), presented in Table 3.14 for trees of orders less than or equal to 2. We shall omit the other details of the calculations and define [3] (t) = e -t y [3] (t), we are finally in a position to represent the differences x(t)-x [3] (t) (left) and z(t)-z [3]

A slow manifold with oscillatory dynamics

The second system we consider has higher dimension (3) and is written previous example, all functions f k and g k vanish at the exception of

f 0 (x, z) = f (x, 0) = 0 -1 1 0 x, f 1 (x, z) = ∂ z f (x, 0)z = -z 0 -1 1 0 x g 0 (x, z) = g(x, 0) = x 2 1 x 2 2 .
Proceeding as for previous example, we obtain the elements of Table 3.15, where we have denoted J = 0 -1 1 0 , and upon computing third-order terms in a similar way we arrive at the following truncated expansion of x ε

0 u ∈ T • • 0 • 1 0 0 0 1 1 0 1 1 1 0 γu (0) 0 1 0 -1 1 1/2
-1 F u (x, y) Jx -y Jx -x y x y x -y 2 x -x 2 1 x 2 2 Jx Table 3.15: Coefficients γ(0) and elementary differentials for trees of T • with order ≤ 2.

x ε 0 = 1 -

ε 2 2 y 2 0 + ε 3 2
x 2 0,1 x 2 0,2 y 0 x 0 +εy 0 + ε 2 x 2 0,1 x 2 0,2 + ε 3 1 6 y 3 0 + 2(1 + 2y 0 )(x 0,1 x 3 0,2x 3 0,1 x 0,2 ) Jx 0 + O ε 4 .

The center manifold is simply here ε h [3] (x) := εx 2 1 x 2 2 -2ε 2 (x 3 1 x 2x 3 2 x 1 ) + ε 3 (-4x 1 x 2 + 2(x 3 1 x 2x 3 2 x 1 ) + 2(x 4 1 -4x 2 1 x 2 2 + x 4 2 )) = εh(x) + O(ε 4 ).

On the right of Figure 3.8 is represented the exact solution of the system ∞,1 (t) and x 2 (t)x

ẋ[3] ∞ = εf x [3] ∞ , εh [3] (x [3] ∞ ) z [3] ∞ = ε h [3] (x
[3] ∞,2 (t), while Figure 3.11 draws the difference z(t)εh [3] (x ∞,1 (t), and between x 2 (t) and x ∞ (t), z(t)) of (3.4.6) for ε = 0.01, x 0 = (0.5; 0.7) T and z 0 = 0.05 (right).

Computing the numerical solution of the truncated normal form of Theorem 3.2.25 leads to Figure 3.11, which represents the components of the following system 3.16 and further computations lead to ε G [3] (x, y) = -2ε 2 (x 1 x 3 2x 3 1 x 2 )(1 + y) + 2ε 3 (x 2 1 x 2 2 (2 + x 1 x 3 2x 3 1 x 2 )x 4 1 + 4x 2 1 x 2 2x 4 2 ).

ẋ[3] ∞ = ε f x [3] ∞ , εh [3] (x [3] ∞ ) ẏ[3] = ε G [3] (x
Finally, coefficients δ(0, t) and elementary differentials required for the expansion of Φ t are presented in Table 3.17 for trees of orders less than or equal to 2. The differences x 1 (t)x

  1

  = h(x h (t)) + O (e -µt ) . (0.1.6) Lemme 0.1.6 Soit (x(t), z(t)) solution de (0.1.2) avec |(x 0 , z 0 )| assez petit. Alors il existe c ≥ 0 et µ > 0 telles que : ∀t > 0, |z(t)h(x(t))| ≤ ce -µt |z 0h(x 0 )|.

0. 1 .

 1 UNE DYNAMIQUE DE VARI ÉT É CENTRALE xv 0.1.2 Une variété centrale périodique rapidement oscillante (Résumé du chapitre 1)

Figure 1 :

 1 Figure 1: x(t) ne décrit pas la dynamique asymptotique de x(t).

Figure 2 :

 2 Figure 2: Construction de x ε 0 avec x(t).

  INTRODUCTION où Ψε t est le flot de l'équation différentielle ẋ = F ε (x). (0.1.13) Le principe d'une méthode de moyennisation est donc celui représenté sur la Figure (4).

Figure 5 :

 5 Figure 5: La dynamique du problème d'interaction proie-prédateur avec migrations rapides.

  Théorème 0.2.10 Soient a, b : T ∪ {∅} → R deux applications, avec a ∅ = 1 et soient B(a, y) et B(b, y) les deux B-séries associées. Alors leur composition est encore une B-série B(b, B(a, y)) = B(a * b, y), avec a * b défini par : ∀u ∈ T , (a * b) u = b u + p∈Pu b p p∈u\p a p + b ∅ a u , (0.2.9) où P u est l'ensemble des sous-arbres 3 de u et pour p ∈ P u , u\p est l'ensemble des branches laissées de côté par ce découpage. Remarque 0.2.11 1 1 est l'élément neutre pour la composition. Exemple 0.2.12 Soit u = , il y a trois manières de le couper -Couper la branche de gauche, alors p = et u\p = {•}, -Couper la branche de droite, alors p = et u\p = {•}, -Couper les deux branches, alors p = • et u\p = {•, •}.

Figure 7 :

 7 Figure 7: Forme normale et équation exacte.

Figure 2 . 1 :

 21 Figure 2.1: The situation when α < 2 (left) and when α > 2 (right).

0 ee

 0 r-1)σ dσds (∂ y f k )g r + O(ε 3 ), (3.2.5) y(t) = y 0 + ε k≥0 t -(k-1)s ds g k + ε 2 -rσ dσds (∂ x g k )f r r-1)σ dσds (∂ y g k )g r + O(ε 3 ).(3.2.6)

  k-1)s-(r-1)σ dσds.

  property worth mentioning is the right-linearity of the * -product: if a ∈ G, and b and b ′ are in C T , then one has a * (b + b ′ ) = a * b + a * b ′ . (3.2.8) A further immediate property that one can easily infer from the formulae above and that we shall frequently use in the sequel is the fact that (a * b)| T• = a * b| T• and (a * b)| T• = a * b| T• .

Proof:

  Writing w(t, τ ) = j k=0 e -kτ W k (t) and ϕ(t, τ ) = r k=0 e -kτ ϕ k (t), the differential equation becomes r k=0 e -kτ (ϕ ′ k (t)kϕ k (t)) = j k=0 e -kτ W k (t).

  Example 3.3.25 With the example above, the tree is written as[[o, [o, o, 1, 1], [o, [o, o, 1, 1], 3, 2], 49, 4], [k, l, m, n]], so the branches are [o, o, 1, 1] and [o, [o, o, 1, 1], 3, 2]. The tree [o, o, 1, 1] is of order 1, so we associate the second index l. The tree [o, [o, o, 1, 1], 3, 2] is of order 2, so we have the indexes [m, n]. Hence, the result of B -is : [ [[o, o, 1, 1], [l]], [[o, [o, o, 1, 1], 3, 2], [m, n]] ].

  Notation 3.3.26 In Maple, the result of ∆ is denoted as :Forest 1 ⊗ Forest 2 + . . . + Forest n 1 ⊗ Forest n 2 = [ [[Forest 1], [Forest 2]], . . . , [[Forest n 1 ], [Forest n 2 ]] ].

Figure 3

 3 Figure 3.3: Solution (x(t), y(t)) of (3.4.1) (left) and solution (x

∞

  (t)) of (3.4.2) (right).

Figure 3 . 3 ∞

 33 Figure3.3 (right) then represents the solution of the following system, which is a third-order approximation of the center-manifold system (3.2.31)ẋ[3] ∞ = εf x [3] ∞ , εh [3] (x [3] ∞ ) ż[3] ∞ = εh [3] (x [3] ∞ )

∞

  (t) (left) and z(t)ε h(x

∞Figure 3 . 4 :

 34 Figure 3.4: Comparison between x(t) and x

∞

  (t) (left) and between z(t) and ε h(x

∞

  (t)) (right).

Figure 3 . 5 :

 35 Figure 3.5: Comparison between x(t) and x

∞

  (t) (left) and between z(t) and ε h(x

∞

  (t)) (right).Eventually, we simulate numerically the solution of the truncated normal form derived in Theorem 3.2.25. More precisely, Figure3.6 represents the solution of the following system 

  Figure 3.6: Solution (x

∞

  (t), z[3] (t)) of (3.4.3) for x 0 = 0.8, z 0 = 0.05 and ε = 0.01.

Figure 3 . 7 :

 37 Figure3.7: Comparison between x(t) and x[3] (t) (left) and between z(t) and z[3] (t) (right).

∞

  Figure 3.8: Exact solution for ε = 0.01, x 0 = (0.5; 0.7) T and z 0 = 0.05 (left) and solution (x

∞

  (t), z

∞

  (t)) of (3.4.5) (right).

Figure 3

 3 Figure3.9: Exact trajectory and its projection onto the (x 1 , x 2 )-plane.

  4.5) with initial condition x[3] ∞ (0) = x ε 0 . Finally, Figure 3.10 represents the quantities x 1 (t)x [3]

∞

  (t)), for several values of ε. They demonstrate a clear-cut numerical confirmation of Theorem 3.2.21.

Figure 3

 3 Figure 3.10: Error between x 1 (t) and x [3]

Figure 3 .

 3 Figure 3.11: Error between z(t) and ε h(x

∞

  (t)) (left) and solution (x

  x ∞ (0) = x ε 0 , ỹ(0) = z 0 and z[3] (t) = e -t ỹ[3]. The coefficients δ(0) associated with εG[2] are listed in Table

1

  (t) (left) and x 2 (t)x

2

  (t) (right) are drawn on Figure3.12, while the difference z(t)z[3] (t) is plot on Figure3.13, for several values of ε. Again, numerical experiments clearly support the assertion of Theorem 3.2.25.

Figure 3

 3 Figure 3.12: Comparison between x 1 (t) and x

  ) (left) and between x 2 (t) and x

Figure 3

 3 Figure 3.13: Comparison between z(t) and z [3] (t).

  e ikθ et β = 0 sinon. En écrivant y solution de (0.2.13) sous forme de B-série

	y(t) = B(α(t), y 0 ),	(0.2.15)
	on montre que α est solution de l'équation différentielle
	α(t) = α(t) * β, α(0) = 1 1.	(0.2.16)
	Une étude de (0.2.16) donne l'existence d'une fonction γ(t, θ) polynomiale en (t, e iθ ) telle que
	∀t ≥ 0, α(t) = γ(t, tω).

Table 3 .

 3

1: Coefficients α for trees u ∈ T • with |u| ≤ 2 and k

  Definition 3.2.6 We say that a function w : R × R → R is a polynomial function if there exists a real polynomialP ∈ R[X 1 , X 2 ] in 2 variables X 1 , X 2 , such that w(t, τ ) = P (t, e -τ ). Furthermore, η : R×R → G is a polynomial map, if (i) η w is a polynomial function for each w ∈ T • ;(ii) e -τ η w is a polynomial function for each w ∈ T • .

	Proposition 3.2.7 There exists a unique polynomial map γ : R × R → G such that
	α(t) = γ(t, t),	(3.2.13)

Lemma 3.2.5 Let w : (t, τ ) ∈ R × R → w(t, τ ) ∈ R be a continuous function which, for each fixed τ , is a polynomial in t and for each fixed t, is a polynomial in e -τ . If for all t ∈ R, w(t, t) = 0, then for all (t, τ ) ∈ R × R, w(t, τ ) = 0.

  .2.14) Proof: By virtue of the chain rule and equation (3.2.9) , equation (3.2.14) is satisfied for τ = t. The mapping γ being polynomial, Lemma (3.2.5) allows to assert that equation (3.2.14) is actually satisfied for all (t, τ ). Finally, γ(0

  proves the second part of (3.2.31). As for the first part of the statement, we consider equation (3.2.28) for trees of T • , differentiate it w.r.t. t ′ , then evaluate at t ′ = 0 and obtain δ∞ | T•

  Table 3.14: Coefficients δ(0, t) and elementary differentials for trees of T with order ≤ 2. • (x [3] ∞ (t), ỹ[3] (t)) = (x(t), y(t)) + O ε 4 , and z

	Denoting
						(x [3] (t), y [3] (t)) = Φ	[3]
		e -t	1 2	e -2t -e -t	1 3	e -3t -1 2	e -2t	1 4	e -4t -1 3	e -3t	e -t -1	1 2	+ 1 2	e -2t -e -t
	F w (x, y)	x 4		1 2	x 4 y 2		x 4 y 2				-1 2	x 4 y 4	-x 2 y	1 2	x 2 y 3
	v ∈ T •	• 0			• 2		0	1					0	3	2	1	2	3
	δ v (0, t)	-1 + e t		1 -e -t		0			-1 6	(1 -e -2t )	-1 2	(1 -e -2t )	-1 12	(1 -e -4t )
	F v (x, y)	x		1 2	xy 2		-x 3 y					1 3	x 3 y 3	1 2	x 3 y 3	-1 6	x 5 y 5

Φ [3] t = Φ t + O ε 4 . t

Dans cette introduction, les mots "fonction régulière" désigneront une fonction de classe C ∞ en toutes ses variables. Dans le reste de la thèse, on travaillera parfois avec des fonctions de classe C r pour r ∈ N * . Le contexte sera toujours précisé.

La théorie des word-series aurait pu être utilisée. Cependant, les B-series sont plus adaptées à l'étude de méthodes numériques générales. De plus, les particularités présentées par les arbres dits dans le chapitre

de norme nulle sont fondamentales pour notre étude et n'auraient pas d'équivalent avec les word-series.

Un sous-arbre de u est un arbre qui possède la même racine que u, obtenu en coupant une ou plusieurs branches selon les découpages admissibles, définis dans[START_REF] Connes | Hopf algebras, renormalization and noncommutative geometry[END_REF].

Il s'agit ici d'une définition un peu modifiée du terme "polynomiale", voir Définition (3.2.6).

La définition exacte distingue les arbres appartenant à T• de ceux de T•.

Throughout this text, the word "smooth" will refer to a C ∞ dependence in the relevant variables. In places however, we may sometimes deal with C r functions for some fixed integer r ≥ 1. This point will be discussed whenever it appears in the text.

Note that here we implicitely use the Floquet Theorem for the resolvent R.

A number of recent publications consider the more involved situation where the matrix Λ exhibits additionally a periodic dependency in time[START_REF] Poggiale | Fast oscillating migations in a predator-prey model[END_REF][START_REF] Castella | Analysis of a time-dependent problem of mixed migration and population dynamics[END_REF].

The coefficients we consider here are time-dependent.

Note that for all λ ∈ R, f k (x, λz) = λ k f k (x, z) and similarly for g k .

This will be done in a forthcoming paper.

This equation is documented in almost every book on the center manifold theorem, see for instance[START_REF] Carr | Applications of Centre Manifold Theory[END_REF].

It is actually obtained as the result of a very high-precision simulation.

In the same way that [CMSS10] has paved the way for [CCMSS11, CMMV15].

Remerciements

Error on q 2 Figure 2.5: Approximation of the solution p 1 , p 2 , q 1 and q 2 for ε = 0, 1 and the error with the exact solution. Error on q 2 Figure 2.6: Approximation of the solution p 1 , p 2 , q 1 and q 2 for ε = 0, 01 and the error with the exact solution.

Then, we use the center manifold to approximate the solution of the initial problem by the solution on the center manifold with initial condition (x ε 0 , εh(x ε 0 , 0)). The error is theoretically an O e -µ t ε + ε 2 . The approximate solution and the error with the exact solution are plotted on Fig. 2.5 for ε = 0.1 and 2.6 for ε = 0.01.

On Fig. 2.7, we zoom on small times to illustrate the phase of exponential decrease of the error. Eventually, we illustrate the evolution of the error as a function of ε. We compute the error between exact and approximate solution with an initial condition on the center manifold. Hence, the error is theoretically a O ε 2 . We illustrate this rate of convergence in Fig. 2.8. A formal series approach to the center manifold theorem Proof: Both assertions can be proved along the same lines and we thus concentrate on the second one. The proof proceeds by induction on the order of trees. For u = • k ∈ T • , k ≥ 0, equation (3.2.15) with the considered initial conditions gives

which, owing to Lemma 3.2.9, has a unique polynomial solution. For v = • k ∈ T • , k ≥ 0, we can conclude similarly by Lemma 3.2.10. Now, consider u = [u 1 , . . . , u p , v 1 , . . . , v q ] • k ∈ T • a tree of order n ≥ 2. As δ is a solution of (3.2.14), we have

where, by induction hypothesis, δ ui (t, τ ) and e -τ δ vj (t, τ ) are polynomial, so that w(t, τ ) is also a polynomial function.

The assumption on the initial conditions and Lemma 3.2.9 thus imply that δ u (t, τ ) is uniquely defined. For a tree v = [u 1 , . . . , u p , v 1 , . . . , v q ] • k ∈ T • of order n ≥ 2, we can conclude similarly using Lemma 3.2.10. As we shall see below, the map δ embeds the dynamics on the center manifold. Its values for trees or orders less than or equal to 2 are given in tables 3.5 and 3.6.

t + e -τ -e -kτ (1+kt) k 2

-(k+r-1)e -kτ +ke -(k+r-1)τ (r-1)k(k+r-1)

Table 3.5: Coefficients δ for trees u ∈ T • with |u| ≤ 2 and k ≥ 1.

Prior to the next proposition, which states the main result of this subsection, we introduce the following key change of variables of the center-manifold theory. Definition 3.2.12 Let γ : R → G be defined by

We denote by Φ 0 the map (x, y) → Φ 0 (x, y) := B(γ -1 (0), (x, y)).

Note that by definition, γ| T• (0) = 1 1| T• , so that y is left unchanged by Φ -1 0 . The values of γ| T• (0) for trees or order less or equal to 2 are given in Table 3.7.

Remark 3.2.13 Denoting x(t, τ ) = B • (γ(t, τ ), (x 0 , y 0 )), the solution x(t, +∞) in the limit τ → +∞ lies on the center manifold and x(0, +∞) is nothing but the modified initial condition x ε 0 . In other words, (x ε 0 , y 0 ) = Φ -1 0 (x 0 , y 0 ).

(3.2.21)

Proof: Premultiplying (3.2.14) by γ-1 (t ′ ) and using the right-linearity of the * -product, one sees that γ-1 (t ′ ) * γ(t + t ′ , τ ) also satisfies an equation of the form (3.2.15) with

By Proposition 3.2.11, the solution is unique and is thus independent of t ′ . In particular, one has

where Φ 0 has been introduced in Definition 3.2.12.

Dynamics on the center manifold

The map δ, being a polynomial map, can be decomposed as

The norm • is defined recursively on T as follows: for all k ∈ N,

Table 3.9: Coefficients δ∞ (0) for trees of T with order ≤ 3.

all terms up to order 3 in ε, we thus obtain

In next proposition, we prove that the center-manifold function εh obtained here satisfies the familiar partial differential equation (3.2.32) 5 .

Proposition 3.2.23

The function h satisfies the following partial differential equation

(3.2.32)

we have ( ẋ∞ (t), ż∞ (t)) = B( δ∞ (t), (x ε 0 , 0)). From Theorem 3.2.21, we also have

To sum up, the following three relations hold for all x

from which we can straightforwardly conclude.

Remark 3.2.24 Note that if one defines the B-series map

hold true. These two relations guarantee the two consistency requirements that: (i) the B-series B(δ ∞ (0), (x, z)), mapping (x, z) to (x, εh(x)) ∈ M, is a projection (onto the manifold, perpendicularly to the x-axis); (ii) the change of variables Φ 0 coincides with the identity map when acting on the manifold M.

The overall results of this subsection are sketched in Figure 3.1. The main novelty here is the map Φ -1 0 , whose expansion in powers of ε we derive, and which tells us how to transform the initial condition.

Figure 3.1: Center-manifold and exact equations.

Reduction to normal form

The center manifold theorem, whether in its standard enunciation 3.1.1 or in its B-series variation 3.2.21, decouples the asymptotic dynamics into a slow variable x -which obeys a nonstiff reduced model-and an enslaved variable z -which becomes a direct function of x. In the transient phase (for small values of t), this remains unsatisfactory since we have no way to recover the full exact solution from the reduced model. The information provided by the initial condition z 0 is indeed lost as soon as Π is applied (see Figure 3.1). If we wish to recover the exact solution of (3.2.1) (through a decoupled -although not reduced-model), we may further exploit Proposition 3.2.14. To this aim, we thus consider the map

which is, roughly speaking, an "unprojected" version of δ ∞ (t). On the one hand, composing (3.2.17) from the left by γ-1 (0), we obtain

so that, using Remark 3.2.19

.2.36)

Then the following relation holds true

(3.2.37)

Proof: As already noticed above, we have [START_REF] Murua | Averaging and computing normal forms with word series algorithms[END_REF] (see Example 1 pp. 28 and Theorem 9). They can be used as well to bring the original system to a normal form which however differs from the one obtained here: mutadis mutandis, system (3.2.36) transforms -via (3.2.19), (3.2.20) and (3.2.21)-into a coupled system in both variables and the gain of working in lower dimension is lost. The corresponding vector field is indeed the B-series with coefficients dγ(t,0) dt t=0

whose first terms can be written as follows

for component ẋ and

for component ẏ. These expansions make apparent the dependence w.r.t. both x and y of both components, in contrast with (3.2.36), which may be thought of as a triangularization.

In order to impart to Theorem 3.2.25 its practical value, we list coefficients δ(0, t), relative to Φ t , for all trees of orders less than or equal to two, in Table 3.10 and Table 3.11. We further illustrate Theorems 3.2.21 and 3.2.25

-(k+r-1)e -kt +ke -(k+r-1)t (r-1)k(k+r-1)

Table 3.10: Coefficients δ(0, t) for trees u ∈ T • with |u| ≤ 2 and k ≥ 1.

schematically in Figure 3.2.

NUMERICAL IMPLEMENTATION OF B-SERIES

57

Figure 3.2: Center-manifold and normal-form equations.

Numerical implementation of B-series

To implement numerically B-series, one key point is to compute the composition and the inverse among B-series. The definitions of the introduction (see Definitions 0.2.10 and 0.2.15) use partition of trees, which is not an easy notion to implement. In this section, the theory of Hopf algebras of trees is introduced in order to define the composition and the inverse of B-series in a computable way.

Hopf algebra of trees

Definition 3.3.1 A Hopf algebra H is a bialgebra , i.e an (unital associative) algebra and a (counital coassociative) coalgebra. Moreover, it is equipped with an antiautomorphism called the antipode S : H → H, such that the following diagram is commutative :

with ∆ : H → H ⊗ H the coproduct and µ : H ⊗ H → H the multiplication.

The composition law can be described using a Hopf algebra introduced by Connes and Kreimer. To explain the notions, the trees with one color and no index are considered. We are back in the initial case of the introduction, where we wanted to describe in terms of B-series the solution of the simple differential system : d t y = εf (y). Definition 3.3.2 A forest is a set of trees in T ∪ {∅}. The forest regrouping the elements w 1 , . . . , w n ∈ T ∪ {∅} (with potential repetitions) is denoted w 1 . . . w n .

Definition 3.3.3 H is the set of linear combinations of forests in T ∪ {∅}.

To attach H a commutative R-algebra structure, we define a scalar multiplication ., an addition + and a multiplication µ : H ⊗ H → H. Example 3.3.4 We have :

which is a linear combination of forests and so belongs to H.

Proposition 3.3.5 The unit element for the multiplication µ is the empty forest ∅ :

In order to make H a bialgebra, a coproduct compatible with the algebras laws is defined. This coproduct is denoted ∆ : H → H ⊗ H and it is coassociative :

The bialgebra H is a Hopf algebra if there exists in addition an antipode S : H → H satisfying :

The map S will be used to compute the inverse in the Butcher group.

The coproduct of Connes and Kreimer

The elements of the Butcher group a ∈ G are such that a ∅ = 1 and they can be extended to unital algebra maps by linearity and by the property :

∀ w 1 , . . . , w n ∈ T , a w1...wn = a w1 . . . a wn .

Hence, a : T ∪ {∅} → R is extended to a : H → R.

Example 3.3.6 Let a : T ∪ {∅} → R be an element of G. The map a is extended to H by :

In order to define the coproduct of Connes and Kreimer, the following two functions are introduced.

Definition 3.3.7 For w = [w 1 , . . . , w n ] • ∈ T , the functions B -: T → H and B + : H → T are defined as:

Using the extension of the Butcher-group elements to unital algebra maps, the definition of the coproduct is the following one. Definition 3.3.9 The coproduct is defined recursively by

(3.3.1)

Example 3.3.10 For trees of order lower or equal to 3, the coproduct gives :

Then, the antipode is defined to complete the construction of Connes and Kreimer Hopf-algebra.

Definition 3.3.11 The antipode S is given by : ∀w ∈ T , S w = p∈P(w) (-1) #P(w) ω∈P(w)

w.

Example 3.3.12 For trees of order lower or equal to 3, this yields to :

The inverse in the Butcher group is expressed using the antipode :

Composition of B-series

The main point of this section is the use of the following proposition. The coproduct of Connes and Kreimer is used in an alternative definition of the composition of B-series.

This characterization of the composition of two B-series, without partitions and subtrees, is the one that we implement on Maple.

The numerical implementation Define a global ordered set of trees

To work numerically with trees, they need to be ordered. In a first step, we work with unindexed bicoloured trees. Definition 3.3.15 The concatenation of two trees u and v, denoted by u.v, is the tree obtained by connecting the root of v to the root of u by a branch.

Example 3.3.16 The concatenation of and • is : .• = .

Definition 3.3.17 The global order ≺ is defined recursively as :

• The tree is of rank 1 and • of rank 2.

• We assume that we have created a set of global ordered trees with all the trees of order ≤ n. They all have a rank in the global order ≺ and the last tree has a rank r.

To create the first tree of order n + 1, we concatenate the tree of rank 1, namely , with the first tree if order n, namely w n : this gives .w n of rank r + 1 in the global order. Then, we concatenate with every other tree of order n, according to their rank in the global order and we attribuate the following rank to the new tree. Once it is done, we have created r n new trees. We define the r + r n + 1 th element of the forest by concatening • (of rank 2) with w n , etc.

Example 3.3.18 The first trees are the following :

If we create the trees as explained above, we will create the same tree several times. Indeed, we have :

which represent the same differential element. We add the condition that a tree is uniquely represented by the tree w = [w 1 , . . . , w n ] • or • with w 1 ≺ . . . ≺ w n . In the previous example, the tree is uniquely represented by :

We define a function GenerateTree(n) which creates all the trees of order lower or equal than n. To do so, we need to represent the trees in a way that Maple can understand. Notation 3.3.19 We write o for a root in • and x for a root in •. Each unindexed tree will be represented using the following model : w = [type of the root, list of branches, rank in the global order, order].

Example 3.3.20 Hence, we have the following definitions :

Compute the coproduct

To compute the composition of B-series, we need to calculate the coproduct ∆ defined as in (3.3.1), but with a small adaptation (due to the presence of bicoloured trees). We modify the Definition 3.3.7. This yields to the coproduct :

Compute all the coefficients

Step1:

We compute β(t):

Step 2: We define α(t) recursively as

Step 3: We define γ(t, τ ) the polynomial map such that α(t) = γ(t, t). To do so, we replace in the expression of α(t) all the e -kt by e -kτ .

Step 4: We define γ(t) as: ∀u ∈ T • , γu (t) = γ u (t, +∞), ∀v ∈ T • , γv (t) = γ u (t, 0).

Step 5: We define the coefficients γ-1 (t). Hence, we have access to Φ -1 0 (x, y) = B(γ(0), (x, y)), the change of variables which sends (x 0 , y 0 ) on (x ε 0 , y 0 ).

Step 6: We define δ(t, τ ) = γ-1 (0) * γ(t, τ ).

Step 7: We define δ ∞ (t) as:

Hence, we have access to (x ∞ (t), z ∞ (t)) = B(δ ∞ (t), (x ε 0 , 0)) and to the center manifold εh(x) = B(δ ∞ (0), (x ε 0 , 0)).

Step 7: We define δ(t) as:

Hence, we have access to (x ∞ (t), ỹ(t)).

Step 8: We define the change of variables Φ t (x, y) = B(δ(0, t), (x, y)). Hence, we can approximate the exact solution (x(t), y(t)) for any t.

Numerical illustration of the results

For the sake of illustration, we next compute the various expansions considered in Theorem 3.2.25 in two simple examples. All expansions are truncared at ε 3 , i.e. ε 4 -terms and smaller are neglected.

Two coupled scalar equations

The first system we consider is of the form

Table 3.16: Coefficients δ(0) and elementary differentials for trees of T • with order ≤ 2.

Table 3.17: Coefficients δ(0, t) and elementary differentials for trees of T with order ≤ 2.

Future work

In the process of reduction of (3.1.1) to its center-manifold form, it is of paramount importance to determine the correct new initial condition x ε 0 . Furthermore, if some information is to be extracted from this reduced system in the transient phase, it is necessary to add a complementary equation and then to apply an appropriate change of variables. All the transformations introduced to this aim in this paper have been developed in terms of B-series, the convergence of which is usually not ensured. It is thus our intention to analyze truncation errors and their impact in a future paper. Of much interest also, in our opinion, is the extension of these results to the situation where problem (3.1.1) involves different eigenvalues λ i . Finally, we think that a better understanding of the structure of the center-manifold equations will lead to new numerical schemes 7 and this is the reason why B-series, which can represent virtually all numerical methods, have been preferred here to word-series.

Abstract :

In this thesis, we study a differential system regulated by two phenomena: a center manifold dynamics and a periodic fast oscillating dynamics. We want to analyse the qualitative behaviour of the system, and to approximate the solution efficiently. In the first chapter, we prove the existence of a fast oscillating center manifold. Then, we prove that the asymptotic behaviour of its solution is given by the shadowed solution on the center manifold, and that it can be approximated up to every order. We use averaging results in order to handle the fast oscillating dynamics. Eventually, we derive a smooth approximated system, without fast oscillations, with the same asymptotic dynamics as the solution of the initial problem. In the second chapter, previous results are applied to a prey-predator system over N distinct sites. The model mixes long time prey-predator interaction (Lotka-Volterra) and fast migrations among sites, with periodic coefficients. A first step is to apply two changes of variables in order to bring this system back to the formalism of the first part. In a second step, we use the results of the first chapter and we derive explicit expansions of the first order approximated systems. At lowest order, it is still of Lotka-Volterra type, with average coefficients, and the terms of higher order perturb this equilibrium. Eventually, these results (both qualitative and quantitative) are illustrated on an example. In the last chapter, we adapt the B-series theory to the study of a simplified version of the system. Firstly, we obtain formal expansions for all the quantities related to the center manifold introduced in the first chapter : this gives informations about the asymptotic behaviour of the system. Secondly, we approximate the solution of the initial system for every time as the composition of a change of variables and the solution of a partially decoupled system. Eventually, we illustrate these results on two examples.