$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, y), \qquad (2.3a) $$

$$ \dot{y} = \epsilon g(y), \qquad (2.3b) $$

où $y$ est maintenant une nouvelle variable d'état dont l'évolution dans le temps est régie par une équation différentielle que l'on suppose du premier ordre (c.f. Éq. (2.3b)). Cette équation fait ici intervenir la fonction $g$ qui définit la forme du profil temporel de $y$.

Comme évoqué en Section 1.2, la présence du paramètre $\epsilon$ dans le membre de droite de l'Éq. (2.3b) fait de l'Éq. (2.3) un système $(n, 1)$-rapide-lent dont les variables rapides sont les coordonnées du vecteur $\mathbf{x}$ qui correspondent aux variables d'états de l'Éq. (2.1) et la variable lente est $y$ le paramètre de bifurcation variant dans le temps.

**REMARQUE 2.1 (CAS STATIQUE ET DYNAMIQUE).** L'Éq. (2.1) représente une situation dite *statique* car le paramètre de bifurcation $\gamma$ est constant dans le temps et l'Éq. (2.3) (ou (2.2)) représente une situation dite *dynamique* car le paramètre de bifurcation $\gamma$ varie (lentement) dans le temps.

Ce vocabulaire provient de théorie de la bifurcation dynamique [15] dont les objets d'étude sont des systèmes de type (2.3) lorsque la variable lente franchit un point de bifurcation de l'Éq. (2.1). Ce dernier est donc parfois appelé point de bifurcation statique pour le différencier de la valeur de $y$ à laquelle émergent les oscillations de l'Éq. (2.3), appelée point de bifurcation dynamique.

De même, le diagramme de bifurcation de l'Éq. (2.1) peut être appelé *diagramme de bifurcation statique* et la représentation des variable rapides $x_i(t)$ ($i = 1, \dots, n$) de (2.3) en fonction de la variable lente $y(t)$ est parfois appelée *diagramme de bifurcation dynamique*.

Dans cette logique, les bassins d'attractions de (2.1) sont appelés *bassins d'attraction statiques* (SBA pour *static basin of attraction*).

Les travaux de recherche résumés dans ce chapitre se placent dans le cadre de la théorie de la bifurcation dynamique et plus généralement dans celui de l'analyse des systèmes dynamiques lents-rapides en étudiant le comportement de l'Éq. (2.3) quand la variable lente $y$ :

(i) croise la valeur du point de bifurcation de Hopf de l'Éq. (2.1). Dans ce cas, les oscillations apparaissent en général pour une valeur de pression dans la bouche (c.-à-d. le point de bifurcation dynamique) plus grande que la valeur statique. Ce phénomène, appelé *retard à la bifurcation* [15], est décrit en Section 2.4.

(ii) traverse une zone de multistabilité du diagramme de bifurcation de l'Éq. (2.1). Dans ce contexte, le travail précurseur de Schumacher et Woodhouse [178] sur les instrument à cordes frottées et plus récemment celui de Colinot et coll. [43] sur un modèle de saxophone montrent, à l'aide de simulations temporelles, que des profil temporels de paramètres de contrôle ayant la même valeur cible peuvent produire des sons différents ; ce qui est difficilement interprétable en termes de bassins d'attraction statique. C'est la raison pour laquelle les concepts de *bassin d'attraction dynamique* (DBA pour *dynamic basin of attraction*) et de séparatrice entre deux DBA sont introduits. Ceci est présenté en Section 2.5.

Pour aider à contextualiser et à comprendre les résultats de recherche présentés en Sections 2.4 et 2.5, un bref état de l'art sur le modèle physique élémentaire d'un instrument à anche est présenté en Section 2.2. Ensuite, l'étude de stabilité statique (locale et globale) du modèle dit « à un mode » (cf. Section 2.2.4.2) est effectuée en Section 2.3.
## 2.2 BREF ÉTAT DE L'ART SUR LE MODÈLE PHYSIQUE ÉLÉMENTAIRE D'UN INSTRUMENT À ANCHE

La Fig 2.1 présente le schéma simplifié d'un bec de type clarinette ou saxophone indiquant les variables et les orientations d'axes choisies utilisées.

<figure><img src="image_36.png" /><figcaption>FIGURE 2.1. Schéma de principe d'un bec à anche simple de type bec de clarinette ou de saxophone. Notations et orientations choisies.</figcaption></figure>

### 2.2.1 L'exciteur

#### 2.2.1.1 La dynamique de l'anche

Dans la littérature il est commun d'assimiler l'anche d'une clarinette (ou d'un saxophone) à un oscillateur harmonique amorti à un degré de liberté [10, 215] ou bien à un simple ressort sans masse ni amortissement [111].

Le premier cas, souvent appelé « modèle à une masse », considère donc que l'anche est un système masse-ressort dont les paramètres sont obtenus à partir de ceux de l'anche lorsque l'on se limite à son premier mode de flexion : la raideur statique (par unité de surface) de l'anche, $K_s$, plus ou moins élevée selon la « force » de l'anche; la première fréquence de résonance $f_r$ de l'anche et enfin la masse modale du premier mode $m_s = K_s/\omega_r^2$, elle aussi exprimée par unité de surface. La dynamique simplifiée de l'anche peut donc être décrite par l'équation suivante :

$$ \frac{d^2 z}{dt^2} + q_r \omega_r \frac{dz}{dt} + \omega_r^2 z = - \frac{\Delta P}{m_s}, \qquad (2.4) $$

où $z(t)$ est le déplacement vertical du bout de l'anche, $\omega_r = 2\pi f_r$ la première pulsation de résonance de l'anche et $q_r\omega_r$ le coefficient d'amortissement. Le terme $\Delta P$ représente la différence de pression s'exerçant sur les deux faces de l'anche et s'exprime comme suit :

$$ \Delta P = P_m - P. \qquad (2.5) $$

Enfin le signe « - » apparaissant dans le membre de droite de l'Éq. (2.4) signifie que l'anche de clarinette est une anche en-dedans [66], c.-à-d. qu'une différence $\Delta P$ positive entre les deux faces de l'anche tend à refermer le canal d'anche.

À part dans le registre suraigu, la fréquence de jeu est en général beaucoup plus basse que la fréquence de résonance $f_r$ de l'anche. Ceci motive l'utilisation du second modèle usuel d'anche qui ignore la dynamique de l'anche ce qui, à partir de (2.4), conduit à :

$$ K_s z = -\Delta P. \qquad (2.6) $$
### 2.2.1.2 Hydrodynamique de l'écoulement à l'entrée du bec

Le comportement de l'air entrant dans le bec de l'instrument est ce qui détermine la non-linéarité de l'instrument essentielle à la création d'auto-oscillations. La modélisation de cet élément consiste à trouver une relation entre le débit $U$ entrant dans l'instrument et la différence de pression $\Delta P$. Cette relation, non linéaire, est appelée caractéristique non linéaire de l'instrument. C'est Backus [10] qui, en 1963, fait les premières mesures de la caractéristique non linéaire. Suite à ses mesures, il propose une expression empirique de la caractéristique dans laquelle le débit est proportionnel à $(\Delta p)^{4/3}$ et à $(z + H)^{2/3}$, où $(z + H)$ représente l'évolution au cours du temps de la hauteur du canal d'anche. Cette expression est ensuite réfutée par Hirschberg et coll. [98, 99], qui proposent l'expression de la caractéristique aujourd'hui la plus couramment utilisée. Elle repose sur l'équation de Bernoulli pour l'écoulement dans le canal entre anche et bec, et la création d'un jet à sa sortie avec dissipation de son énergie cinétique par turbulence sans récupération de la pression. La vitesse $v_j$ du jet à la sortie du canal d'anche s'écrit dans ce cas :

$$v_j = \sqrt{\frac{2\Delta P}{\rho}}, \qquad (2.7)$$

où $\rho$ est la masse volumique de l'air.

Pour déterminer le débit $U$ entrant dans l'instrument à partir de l'expression (2.7) il nous faut connaître la section du jet $S_j$ à la sortie du canal d'anche. Si la séparation de l'écoulement est supposée localisée à la sortie du canal d'anche (en réalité la position du point de séparation dépend de la géométrie de la sortie du canal d'anche [98]), la section du jet est proportionnelle à la hauteur $(z + H)$ du canal. Le plus simple est de considérer que la section du jet est directement égale à la section transverse du canal d'anche, dans ce cas $S_j$ s'écrit :

$$S_j = (z + H)w, \qquad (2.8)$$

où $w$ est la largeur du canal d'anche. Cette proportionnalité entre la section du jet $S_j$ et l'ouverture du canal d'anche $(z + H)$ est mesurée par Dalmont et coll. [50]. Cependant, dans le cas d'un canal court ou pour de grandes ouvertures, un phénomène de *vena-cantracta* peut apparaître (voir [4] section 5.10 et annexes C.3 et C.5). Dans ce cas, la section $S_j$ du jet est inférieure à la section transverse du canal (cf. [73], Fig.I-5). Ceci a été observé sur la clarinette par mesure Doppler de vitesse et visualisation de l'écoulement [227]. Si nous considérons un bec suffisamment fermé et disposant d'un canal d'anche long, le débit à la sortie du canal peut finalement s'écrire :

$$U = (z + H)w \sqrt{\frac{2\Delta P}{\rho}}. \qquad (2.9)$$

L'Éq. (2.9) donne l'expression du débit lorsqu'il entre dans le bec. Pour la compléter il faut considérer le cas inverse d'un débit entrant dans la bouche de l'instrumentiste. Ce phénomène est rare numériquement et n'a encore jamais été observé expérimentalement. Il convient cependant de le prendre en compte de façon à avoir une expression la plus générale possible (dans le cadre des hypothèses énoncées précédemment). En admettant que les hypothèses précédentes restent valables dans le cas où la différence de pression est inversée, l'expression du débit $U$ en fonction de la différence de pression $\Delta P$ s'écrit finalement :

$$U = U_a \left(1 + \frac{z}{H}\right) \sqrt{\frac{|\Delta p|}{P_M}} \operatorname{sgn}(\Delta P), \qquad (2.10)$$

où sgn représente la fonction signe et la grandeur $U_a$ est proportionnelle à un débit et vaut :

$$U_a = wH \sqrt{\frac{2}{\rho} P_M}, \qquad (2.11)$$

$P_M$ est la valeur de la pression dans la bouche pour laquelle l'anche plaque ($z = -H$, voir la Fig. 2.1) et reste plaquée (la pression dans le bec devenant alors nulle). À partir de l'Éq. (2.6) on a donc
$P_M = K_s/H$; appelée pression de plaquage. Au delà de cette valeur le canal d'anche est fermé (ceci n'a en fait jamais été observé dans [50]) et le débit $U$ devient par conséquent égal à zéro, il vient donc l'expression finale du débit, fonction de la différence de pression $\Delta P$ et du déplacement du bout de l'anche $z$ :

$$U = U_A \left(1 + \frac{z}{H}\right) \sqrt{\frac{|\Delta P|}{P_M}} \operatorname{sgn}(\Delta P)\Theta(z + H), \qquad (2.12)$$

où $\Theta$ est la fonction de Heaviside. Un exemple de caractéristique non linéaire théorique calculée d'après l'Éq. (2.12) est montré en Fig. 2.2.

<figure><img src="image_37.png" /><figcaption>FIGURE 2.2. Caractéristique non linéaire calculée d'après (2.12). $H = 1\text{mm}$, $w = 1.3\text{cm}$, $\rho = 1.2\text{kg/m}^3$ et $P_M = 10\text{kPa}$.</figcaption></figure>

Nous devons maintenant connaître le débit total $U_{in}$ entrant dans le résonateur. Le débit total n'est pas directement le débit $U$ calculé précédemment, il faut lui ajouter le débit $U_r$ créé par l'anche qui agit comme un piston aspirant et repoussant de l'air à sa vitesse. En prenant en compte l'orientation choisie (cf. Fig. 2.1), le débit $U_r$ s'écrit alors :

$$U_r = -S_r \frac{dz}{dt}, \qquad (2.13)$$

où $S_r$ est la surface effective de l'anche. Finalement le débit total entrant dans le résonateur de l'instrument a pour expression :

$$U_{in} = U - S_r \frac{dz}{dt}. \qquad (2.14)$$

### 2.2.2 Le résonateur

Le corps d'un instrument à vent est un résonateur acoustique. Dans l'hypothèse d'une propagation linéaire et en ondes planes, son comportement en régime forcé peut être décrit dans le domaine fréquentiel par son impédance d'entrée $Z$ :

$$Z(\omega) = \frac{\hat{P}(\omega)}{\hat{U}_{in}(\omega)}, \qquad (2.15)$$

où $\hat{P}(\omega)$ et $\hat{U}_{in}(\omega)$ sont respectivement les transformées de Fourier de $P(t)$ et $U_{in}(t)$.

Le module $|Z|$ de l'impédance d'entrée mesurée sur un saxophone alto (doigté de Si grave) est représenté en Fig. 2.3.
<figure><img src="image_38.png" /><figcaption>**FIGURE 2.3.** Figure extraite de [42] (Fig. 2). Module de l'impédance d'entrée mesurée sur un saxophone alto (doigté de si grave). Le module de l'impédance est normalisé par l'impédance caractéristique à l'entrée de l'instrument.</figcaption></figure>

**DÉFINITION 2.1 (REGISTRE).** Dans un contexte musical, si la fréquence de la note jouée par l'instrument à anche pour un doigté donné est proche du premier pic de $|Z|$, on parle du premier registre de l'instrument. Si pour le même doigté la fréquence de la note jouée est proche d'un des pics suivants, on parle de registres d'ordre supérieurs : deuxième registre pour le deuxième pic, troisième registre pour le troisième pic, etc.

Dans le domaine temporel l'expression (2.15) s'écrit :

$$U_{in}(t) = [h * P](t), \qquad (2.16)$$

où $h(t)$ est la réponse impulsionnelle du résonateur, transformée de Fourier inverse de l'admittance d'entrée $Y(\omega)$ (inverse de l'impédance d'entrée) [37] (Chapitre 9). Il est souvent plus commode de travailler avec la variable $U$ qui, à partir des Éqs. (2.14) et (2.16), s'écrit :

$$U(t) = [h * P](t) + S_r \frac{dz}{dt}. \qquad (2.17)$$

### 2.2.3 Modèles adimensionnés d'un instrument à anche simple

#### 2.2.3.1 Modèle à trois équations.

Les Éqs. (2.4), (2.12) et (2.17) forment un système de trois équations à trois inconnues ($z$, $P$ et $U$) modélisant le comportement d'un instrument de musique à anche simple. Les paramètres de contrôle (ou de bifurcation) sont la pression dans la bouche $P_m$ et la force d'appui de la lèvre du musicien sur l'anche, reliée au débit $U_a$. Pour étudier ce modèle il est judicieux de le réécrire en variables adimensionnées [36, 111]. Les déplacements sont adimensionnés par l'ouverture au repos du bec $H$, les pressions par la pression de plaquage de l'anche $P_M$ et les débits par $P_M/Z_c$, où $Z_c = \rho c/S$ est l'impédance caractéristique des ondes planes dans ce résonateur. Les variables et paramètres adimensionnés s'écrivent donc :

$$ \begin{aligned} \Delta p &= \Delta P/P_M &; \quad p &= P/P_M \\ u &= Z_c U/P_M &; \quad \gamma &= P_m/P_M \\ x &= z/H + \gamma &; \quad \zeta &= Z_c U_A/P_M. \end{aligned} \qquad (2.18) $$
En utilisant les variables adimensionnées définies par (2.18), le système formé des Éqs. (2.4), (2.12) et (2.17) s'écrit maintenant :

$$ u = [h * p] + \frac{Z_c H S_r}{p_M} \frac{dz}{dt} \qquad (2.19a) $$

$$ \frac{1}{\omega_r^2} \frac{d^2 x}{dt^2} + \frac{q_r}{\omega_r} \frac{dx}{dt} + x = p \qquad (2.19b) $$

$$ u = \zeta(1 + x - \gamma) \sqrt{|\gamma - p|} \operatorname{sgn}(\gamma - p) \Theta(1 + x - \gamma). \qquad (2.19c) $$

Le modèle précédent est appelé communément « modèle à 3 équations ». Les variables de sortie sont $p$, $u$ et $z$ et les paramètres de contrôle sont $\gamma$ et $\zeta$.

### 2.2.3.2 Modèle à deux équations

Comme mentionné en Section 2.2.1, la fréquence du registre fondamental des notes jouées par un instrument à anche simple comme la clarinette ou le saxophone est toujours inférieure à la fréquence de résonance $f_r$ de l'anche. Ceci permet de s'affranchir de la dynamique de l'anche. Dans ce cas, en notant que $p_M = K_s H$, les Éqs. (2.6) et (2.12) s'écrivent respectivement :

$$ z = -H \frac{\Delta p}{P_M}, \qquad (2.20) $$

et

$$ U = U_A \left( 1 + \frac{\Delta P}{P_M} \right) \sqrt{\frac{|\Delta P|}{P_M}} \operatorname{sgn}(\Delta p) \Theta(P_M - \Delta P). \qquad (2.21) $$

L'Éq. (2.21) indique que la caractéristique non linéaire est à dérivée discontinue en $\Delta P = P_M$ ce qui d'un point de vue physique n'est pas tout à fait satisfaisant. Cependant, les résultats expérimentaux [48, 50] confirment la partie de la courbe correspondant à l'anche non battante ($\Delta P < P_M$).

Pour le cas d'un tuyau cylindrique, le débit $U_r$ créé par l'anche est supposé avoir pour effet une simple correction de longueur du résonateur qui est ajoutée à la longueur totale [36]. L'Éq. (2.19a) du résonateur s'écrit dans ce cas :

$$ U(t) = [h * P](t). \qquad (2.22) $$

Les Éqs. (2.21) et (2.22) constituent le « modèle à 2 équations ». En variables adimensionnées il s'écrit :

$$ u = [h * p] \qquad (2.23a) $$

$$ u = F(p) = \zeta(1 + p - \gamma) \sqrt{|\gamma - p|} \operatorname{sgn}(\gamma - p) \Theta(1 + p - \gamma) \qquad (2.23b) $$

La Fig. 2.4(a) représente la fonction $F$ pour les paramètres de contrôle $\gamma = 0.42$ et $\zeta = 0.6$.

### 2.2.3.3 Représentation en onde aller et onde retour

La solution générale de l'équation d'onde, qui régit la propagation des ondes acoustiques (pression et débit) dans le résonateur de l'instrument, s'exprime (en supposant une propagation à une dimension) sous la forme d'une onde se propageant vers la droite, l'onde aller, et d'une onde se propageant vers la gauche, l'onde retour. La pression acoustique est la somme de ces deux ondes progressives.

Le support temporel de la réponse impulsionnelle $h(t)$ est long, ce qui peut poser des problèmes si l'on veut résoudre numériquement le système. Pour remédier à ce problème il est judicieux de considérer une représentation en onde de pression aller $p^+$ et onde de pression retour $p^-$ [137, 177] plutôt que celle en pression $p$ et débit $u$. Comme nous travaillons avec des variables adimensionnées, les ondes aller et retour sont définies ainsi :

$$ p = p^+ + p^- ; \quad u = p^+ - p^-, \qquad (2.24) $$
<figure><img src="image_39.png" /><figcaption>FIGURE 2.4. Représentation de la caractéristique non linéaire en variables $[p, u]$ (a) et $[p^+, p^-]$ (b) avec comme paramètres de contrôle $\gamma = 0.42$ et $\zeta = 0.6$.</figcaption></figure>

soit,

$$p^+ = \frac{1}{2}(p+u); \quad p^- = \frac{1}{2}(p-u). \qquad (2.25)$$

Dans cette représentation l'équation du résonateur (2.23a) devient :

$$p^{-}(t) = [r * p^{+}](t), \qquad (2.26)$$

où $r(t)$ est la fonction de réflexion du résonateur qui n'est autre que la transformée de Fourier inverse du coefficient de réflexion $R(\omega)$, défini par :

$$\hat{P}^{-}(\omega) = R(\omega) \hat{P}^{+}(\omega). \qquad (2.27)$$

Quant à la relation non linéaire (2.23b), elle devient :

$$p^{+} = G(-p^{-}). \qquad (2.28)$$

La fonction $G$ est obtenue par rotation de 45° de la caractéristique non linéaire $F$ [133]. Il est aussi possible d'obtenir une expression analytique de la fonction $G$. L'obtention de cette expression analytique repose sur le changement de variable $X = \gamma - p^+ - p^-$ et $Y = \gamma - 2p^-$, et sur la résolution d'une équation du troisième degré (cf. Annexe A de [195]). En utilisant les variables $p^+$ et $p^-$ le modèle à deux équations devient :

$$\left\{ \begin{array}{l} p^{-} = [r * p^{+}] \\ p^{+} = G(-p^{-}) \end{array} \right. \qquad (2.29a)$$

$$\left\{ \begin{array}{l} p^{+} = G(-p^{-}) \\ p^{-} = [r * p^{+}] \end{array} \right. \qquad (2.29b)$$

La Fig. 2.4(b) représente la fonction $G$ pour les paramètres de contrôle $\gamma = 0.42$ et $\zeta = 0.6$.

### 2.2.4 Modèles simples utilisés

Les travaux présentés dans ce manuscrit s'appuient essentiellement sur deux modèles simplifiés d'instruments à anche simple décrits dans cette section.

#### 2.2.4.1 Modèle de « Raman »

Le modèle dit « de Raman » (pour la clarinette [49]) est l'un des modèles les plus simples d'instrument à anche. Il tient son nom de l'analogie avec le modèle de corde frottée proposé par le physicien indien C. V. Raman [166].

Dans un instrument de type clarinette, deux types de pertes sont à considérer<fnref n="3" /> : les pertes visco-

<footnote n="3">3. Les pertes volumiques dans l'air sont en effet négligeables par rapport à ces deux types de pertes.</footnote>
thermiques à l'intérieur du résonateur et les pertes par rayonnement à l'extrémité de ce dernier. Aux basses fréquences les pertes visco-thermiques sont prépondérantes et elles sont, dans le modèle, prises en compte dans la fonction de réflexion; dans [164] Polack et coll. proposent une expression analytique réaliste de la fonction de réflexion d'un tube cylindrique.

Le modèle de Raman fait l'hypothèse d'un résonateur cylindrique dont l'impédance est nulle à l'extrémité aval (réflexion totale avec changement de signe pour la pression, ce qui est à peu près valable aux basses fréquences) et considère que les pertes visco-thermiques à l'intérieur du résonateur sont indépendantes de la fréquence. Dans ce cas la fonction de réflexion $r(t)$ correspond à un retard $\tau$, dû à un aller-retour dans le résonateur, multiplié par un coefficient réel positif $\lambda$ qui traduit la diminution d'amplitude due aux pertes :

$$r(t) = -\lambda\delta(t - \tau), \qquad (2.30)$$

où $\delta$ est la distribution de Dirac et $\tau = 2l/c$ est le temps que met une onde de célérité $c$ à faire un aller-retour dans le cylindre de longueur $l$. Le paramètre de perte $\lambda$, dont une valeur typique est 0.9, a pour expression :

$$\lambda = e^{-2\alpha l}, \qquad (2.31)$$

où $\alpha$ est le coefficient d'amortissement [110] :

$$\alpha \approx 3 \cdot 10^{-5} \sqrt{f/R}. \qquad (\text{SI}) \qquad (2.32)$$

$R$ est le rayon du résonateur et $f$ est la fréquence considérée. Dans le cas du modèle de Raman le coefficient d'amortissement $\alpha$ est constant et communément calculé à la fréquence de jeu [48].

Dans ce cas l'Éq. (2.26) s'écrit simplement :

$$p^{-}(t) = -\lambda p^{+}(t - \tau), \qquad (2.33)$$

En remplaçant dans (2.28), on obtient :

$$p^{+}(t) = G(\lambda p^{+}(t - \tau)) = f(p^{+}(t - \tau)), \qquad (2.34)$$

qui, en discrétisant à tous les intervalles $\tau$ et en notant $p^{+}(n\tau) = p_{n}^{+}$, s'écrit :

$$p_{n}^{+} = f(p_{n-1}^{+}). \qquad (2.35)$$

L'expression de la fonction $f$ est obtenue dans [195]. Le modèle de clarinette de Raman se simplifie donc en l'équation aux différences (2.35) (à temps discret) qui décrit l'évolution au cours du temps de la variable $p^{+}$ et qui peut être résolue par un simple processus itératif. Connaissant $p^{+}$, il est ensuite possible d'obtenir $p^{-}$ par (2.33) puis $p$ et $u$ par (2.24).

Le modèle communément appelé « modèle sans pertes » peut être obtenu simplement en posant $\lambda = 1$ ce qui signifie que les pertes dans le résonateur sont elles aussi ignorées.

En supposant des paramètres de contrôle $\zeta$ et $\gamma$ constants, le modèle de Raman (ou sans pertes) permet d'expliquer la naissance des oscillations dans un instrument à anche simple. Il s'agit de la bifurcation à doublement de période (cf. [119] et la Définition 2.2 ci-dessous) de la solution d'équilibre de l'Éq. (2.35) (obtenue en résolvant $x = f(x)$ et correspondant au silence) vers une solution périodique à deux états (obtenus en résolvant $x = f \circ f(x)$ et correspondant à un son musical).

**DÉFINITION 2.2 (BIFURCATION À DOUBLEMENT DE PÉRIODE).** Soit l'équation au dif-
férence :

$$x_n = f(x_{n-1}), \qquad (2.36)$$

Si $x^*$ et $f'$ sont respectivement le point fixe (solution de $x^* = f(x^*)$) et la dérivée de la fonction d'itération $f$, alors la bifurcation associée à $f'(x^*) = -1$, est appelée bifurcation à doublement de période. En anglais, on parle de *period-doubling bifurcation* ou de *flip bifurcation*.

Une analyse de bifurcation statique (c.-à-d. à paramètre de bifurcation $\zeta$ et $\gamma$ constants) relativement détaillée du modèle de Raman est présentée en Section 1.3 du manuscrit de ma thèse de doctorat [t1].
### 2.2.4.2 Modèle à un mode

Nous considérons ici la décomposition modale de l'impédance d'entrée (2.15) du résonateur de l'instrument. Cette décomposition modale infinie s'écrit sous la forme suivante

$$Z(\omega) = Z_c \sum_n Z_n \frac{j\epsilon_n\omega_n\omega}{\omega_n^2 + j\epsilon_n\omega_n\omega - \omega^2}, \qquad (2.37)$$

où $Z_n$, $\omega_n$ and $\epsilon_n$ sont des paramètres modaux, respectivement le facteur modal, la pulsation de résonance et l'inverse du facteur de qualité du n-ième pic de l'impédance (correspondant au n-ième mode de la colonne d'air du résonateur).

En utilisant la transformée de Fourier de l'Éq. (2.22) et les variables adimensionnées (2.18), nous avons aussi $Z(\omega)/Z_c = \hat{u}(\omega)/\hat{p}(\omega)$, où $\hat{u}(\omega)$ et $\hat{p}(\omega)$ sont les transformées de Fourier de $u(t)$ et $p(t)$, respectivement. Par conséquent en effectuant la transformée de Fourier inverse de (H.22) nous obtenons les équations suivantes, écrites dans le domaine temporel :

$$\frac{d^2 p_n}{dt^2} + \epsilon_n \omega_n \frac{dp_n}{dt} + \omega_n^2 p_n = Z_n \epsilon_n \omega_n \frac{du}{dt}, \quad \forall n \in \mathbb{N}, \qquad (2.38)$$

Le modèle dit « à un mode » est obtenu en ne considérant que le premier mode du résonateur, c.-à-d. en prenant uniquement $n = 1$ dans l'Éq. (2.38). En notant $p_1 = p$, en utilisant l'Éq. (2.23b) et en passant au temps adimensionné $t' = \omega_1 t$, nous obtenons l'équation différentielle suivante :

$$\ddot{p} + \epsilon_1 (1 - Z_1 F'(p)) \dot{p} + p = 0 \qquad (2.39)$$

où $F'(p) = dF(p)/dp$ et où le temps adimensionné $t'$ est noté $t$ et $d()/dt = (\dot{)}$ pour alléger les notations.

## 2.3 RÉSULTATS PRÉLIMINAIRES : ANALYSE DE STABILITÉ STATIQUE (LOCALE ET GLOBALE) DU MODÈLE À UN MODE

Après avoir utilisé le modèle de Raman (à temps discret) pendant ma thèse, mes recherches actuelles se focalisent surtout sur des modèles à temps continu et notamment sur le modèle à un mode. Pour aider à contextualiser et à comprendre les résultats récents présentés en Sections 2.4 et 2.5, l'étude de stabilité statique (locale et globale) du modèle à un mode est effectuée dans cette section. Nous rappelons que le paramètre de bifurcation considéré ici, comme dans l'ensemble de ce document, est la pression dans la bouche adimensionnée $\gamma$. Sauf mention contraire, les paramètres du modèles de clarinette à un mode (2.39) utilisés dans les sections qui suivent sont donnés dans le Tab. H.1.

### 2.3.1 Analyse de stabilité locale : naissance des oscillations

Dans le cadre du modèle à un mode, l'apparition d'un son musical (c.-à-d. ici un régime périodique) dans un instrument à anche simple s'explique par la perte de stabilité de la solution d'équilibre triviale de l'Éq. (2.39) (c.-à-d. $p = 0$ et $\dot{p} = 0$) au profit d'une solution périodique.

Pour déterminer la stabilité de la position d'équilibre il est judicieux d'écrire l'Éq. (2.39) dans l'espace des phase ($p = x_1, \dot{p} = x_2$), ce qui donne :

$$\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \gamma), \qquad (2.40)$$

avec $\mathbf{x} = (x_1, x_2)^T$ et $\mathbf{f}(\mathbf{x}, \gamma) = (x_2, -\epsilon_1(1 - Z_1F'(x_1))x_2 + x_1)^T$ (avec $()^T$ le vecteur transposé du vecteur ($$)). À noter que l'Éq. (2.40) est de la même forme que l'Éq. (2.1) avec $\mathbf{x} \in \mathbb{R}^2$.

Le perte de stabilité de la solution d'équilibre s'effectue via une bifurcation de Hopf (voir [179] et la Définition 2.3 ci-dessous).
**DÉFINITION 2.3 (BIFURCATION DE HOPF).** Une bifurcation de Hopf est une bifurcation locale dans laquelle un point fixe d'un système dynamique perd sa stabilité tandis qu'une paire de valeurs propres complexes conjuguées de la linéarisation autour de ce point fixe franchissent l'axe imaginaire du plan complexe (c.-à-d. les parties réelles deviennent positives). Si nous considérons un système dynamique du plan de la forme suivante :

$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \lambda) \qquad (2.41) $$

avec $\mathbf{x} \in \mathbb{R}^2$, $\lambda \in \mathbb{R}$ le paramètre de bifurcation considéré et $\mathbf{f}$ une fonction vectorielle non linéaire suffisamment dérivable. Alors les conditions d'une bifurcation de Hopf sont :

1. $\mathbf{f}(\mathbf{x}_0, \lambda) = 0$, c.-à-d. $\mathbf{x}_0$ est une solution d'équilibre (ou point fixe) de (2.41);

2. $D_{\mathbf{x}}\mathbf{f}(\mathbf{x}_0, \lambda)$ (la matrice jacobienne de $\mathbf{f}$ évaluée en $\mathbf{x}_0$) à une paire de valeurs propres complexes conjuguées $\mu(\lambda) = \alpha(\lambda) \pm j\beta(\lambda)$ dont la partie réelle s'annule pour $\lambda = \lambda_0$, c.-à-d. $\alpha(\lambda_0) = 0$;

3. $\frac{d\alpha(\lambda_0)}{d\lambda} \neq 0$.

Si les conditions précédentes sont respectées, il y a alors naissance d'une solution périodique en $(\mathbf{x}_0, \lambda_0)$ dont la pulsation initiale est $\beta(\lambda_0)$. $\lambda_0$ est appelé point de bifurcation de Hopf statique

Lorsque l'amplitude de la pression $p$ dans le bec de l'instrument est proche de zéro, l'Éq. (2.23b) s'écrit s'implement :

$$ u = F(p) = \zeta(1 + p - \gamma)\sqrt{\gamma - p}. \qquad (2.42) $$

Pour expliquer l'apparition d'une solution périodique de l'Éq. (2.39) par déstabilisation de la solution d'équilibre triviale, il est courant d'utiliser le développement limité de (2.42) à l'ordre trois<fnref n="4" /> autour de $p = 0$. Dans ce cas la relation débit pression devient :

$$ u = F(p) \approx c_0 + c_1 p + c_2 p^2 + c_3 p^3, \qquad (2.43) $$

avec $c_0 = \zeta(1 - \gamma)\sqrt{\gamma}$, $c_1 = \zeta \frac{3\gamma - 1}{2\gamma^2}$, $c_2 = -\zeta \frac{3\gamma + 1}{8\gamma^2}$ et $c_3 = -\zeta \frac{\gamma + 1}{16\gamma^2}$. L'approximation (2.43) est pertinente si l'on considère des solutions périodiques d'amplitude modérée.

En utilisant l'Éq. (2.43), la matrice jacobienne $\mathbf{f}_{\mathbf{x}}(\mathbf{x}, \gamma)$ de (2.40) estimée au point fixe trivial $\mathbf{x}_0 = \mathbf{0}$ s'écrit :

$$ D_{\mathbf{x}}\mathbf{f}(\mathbf{0}, \gamma) = \begin{pmatrix} 0 & 1 \\ -1 & B(\gamma) \end{pmatrix} \qquad (2.44) $$

avec

$$ B(\gamma) = -\epsilon_1 \left( 1 - \frac{(3\gamma - 1)\zeta Z_1}{2\sqrt{\gamma}} \right). \qquad (2.45) $$

Le calcul des valeurs propres $\mu(\gamma)$ de (2.44) donne :

$$ \mu(\gamma) = \frac{1}{2} \left( B(\gamma) \pm j \sqrt{4 - B(\gamma)^2} \right), \qquad (2.46) $$

si $-2 < B(\gamma) < 2$, ce qui se vérifie dans le cas d'un amortissement faible. Selon la Définition 2.3, nous avons une bifurcation de Hopf quand $\gamma = \hat{\gamma}^{\text{st}}$, le point de bifurcation statique, solution de $B(\hat{\gamma}^{\text{st}}) = 0$. On montre facilement en résolvant $B(\hat{\gamma}^{\text{st}}) = 0$ que :

$$ \hat{\gamma}^{\text{st}} = \frac{1}{3} + \frac{2 + 2\sqrt{1 + 3\zeta^2 Z_1^2}}{9\zeta^2 Z_1^2}. \qquad (2.47) $$

<footnote n="4">4. L'ordre trois est nécessaire pour que le modèle possède une solution périodique en plus de la solution d'équilibre.</footnote>
À noter que si le pertes dans le résonateur sont ignorées (ce qui revient dans notre modèle à $Z_1 \to +\infty$) alors $\hat{\gamma}^{st} = 1/3$.

Toujours en utilisant l'approximation polynomiale (2.43) de la caractéristique non linéaire du débit, il est possible d'obtenir analytiquement (de façon approchée) le diagramme de bifurcation statique du modèle à un mode. Pour cela une méthode de moyennisation est employée. Il existe dans la littérature plusieurs manières d'aborder la moyennisation. La plus courante en mécanique est d'utiliser l'approximation de Krylov-Bogolyubov, telle qu'évoquée en Section 1.2.1.3. L'autre (voir Chapitre 7 de [172]) est d'utiliser la nature également lente-rapide de la représentation amplitude-phase du mouvement. Cette seconde méthode, appliquée au modèle à un mode (2.39), est détaillée en Section H.9 de ce manuscrit 5 et résumée très brièvement ici. D'abord, l'Éq. (2.39) est écrite comme suit :

$$ \ddot{p} + \epsilon_1 h(p, \dot{p}, \gamma) + p = 0. \quad (2.48) $$

Ensuite, nous utilisons une représentation complexe similaire à celle proposée en Section 1.2.1.3, c.-à-d. en posant $\xi e^{jt} = \dot{p} + j\bar{p}$ puis $\xi = x e^{j\varphi}$ où $x$ et $\varphi$ représentent respectivement l'amplitude et la phase de $p$. En combinant $\xi e^{jt} = \dot{p} + j\bar{p}$ et son complexe conjugué, $p$ et $\dot{p}$ s'écrivent en effet :

$$ p = \frac{\xi e^{jt} - \xi^* e^{-jt}}{2j} = x \sin(t + \varphi) \quad (2.49a) $$

$$ \dot{p} = \frac{\xi e^{jt} + \xi^* e^{-jt}}{2} = x \cos(t + \varphi) \quad (2.49b) $$

On montre alors, en supposant un faible amortissement (c.-à-d. $0 < \epsilon_1 \ll 1$), que $x$ et $\varphi$ sont découplés, que l'équation différentielle de l'amplitude $x$ peut être approchée par :

$$ \dot{x} = \epsilon_1 F(x, \gamma) = \epsilon_1 x \left( -\frac{1}{2} - \frac{\zeta Z_1 (32\gamma^2(1-3\gamma) + 3(\gamma+1)x^2)}{128\gamma^{5/2}} \right). \quad (2.50) $$

et que la phase $\varphi$ est constante. Nous montrons que $\partial_x F(0, \gamma) = \frac{2B(\gamma)}{\epsilon_1}$ (avec $B(\gamma)$ donné par l'Éq. (2.45)). La solution triviale de (2.50) subit donc une bifurcation *fourche* (voir par exemple [119]) correspondant à la bifurcation de Hopf de la solution triviale de (2.48) (cf. Éq. (2.46)), c.-à-d. pour $\gamma = \hat{\gamma}^{st}$. Les solutions d'équilibre non triviales de (2.50) correspondent aux solutions périodiques de (2.48). Les deux solutions non nulles de $F(x, \gamma) = 0$, notée $x_{-}^e$ et $x_{+}^e$ sont :

$$ x_{\pm}^{e} = \pm \frac{4\sqrt{\frac{2}{3}\gamma\sqrt{(3\gamma-1)\zeta Z_1 - 2\sqrt{\gamma}}}}{\sqrt{(\gamma+1)\zeta Z_1}} \quad (2.51) $$

qui sont réelles quand $\gamma > \hat{\gamma}^{st}$ et qui correspondent donc à deux solutions périodiques de (2.48) de même amplitude mais de phase opposée. Ces deux solutions sont associées à la même énergie représentée par $x^2$. De plus, d'un point de vue musical il s'agit d'une même note car ces deux solutions oscillent à la même pulsation, avec donc la même énergie.

La stabilité des solutions (2.51) est obtenue en regardant le signe de $\partial_x F(x_{\pm}^e, \gamma)$ : $x_{\pm}^e$ sont stables si $\partial_x F(x_{\pm}^e, \gamma) < 0$ et instables si $\partial_x F(x_{\pm}^e, \gamma) > 0$. On montre que :

$$ \partial_x F(x_{\pm}^e, \gamma) = - \frac{B(y)}{\epsilon_1} \quad (2.52) $$

qui est donc négatif quand $\gamma > \hat{\gamma}^{st}$.

Les solutions (2.51) sont donc stables quand elles existent (c.-à-d. quand elles sont réelles). Il s'agit donc d'un cas de bistabilité. Pour un tel système à une dimension les bassin d'attraction statiques (SBA encore pour *static basin of attraction*) des solutions stables sont triviaux. En effet, si $x_0 = x(0) < 0$ (resp. $x_0 > 0$), alors c'est $x_{-}^e$ (resp. $x_{+}^e$) qui sera atteinte. Le phénomène de bistabilité

<footnote n="5">5. Attention en Section H.9 le tempts $t$ est noté $t'$.</footnote>
<figure><img src="image_40.png" /><figcaption>**FIGURE 2.5.** Diagramme de bifurcation de (2.50). La solution d'équilibre non triviale $x = 0$ correspond également à la solution triviale de (2.48), et les solutions d'équilibre non triviales $x_+^e$ et $x_-^e$ (voir (2.51)) représentent l'amplitude des solutions périodiques de (2.48).</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**FIGURE 2.6.** Exemples d'intégration numérique des Éqs. (2.48) et (2.50) : series temporelle de $x(t)$ (vert) et $p(t)$ (en bleu). Mêmes paramètres qu'en Fig. 2.5 avec (a) $\gamma = 0.2 < \hat{\gamma}^{\text{st}}$ (c.-à-d. la solution triviale est stable) et (b) $\gamma = 0.6 > \hat{\gamma}^{\text{st}}$ (c.-à-d. la solution triviale est instable et la solution périodique est stable).</figcaption></figure>

sera abordé plus en détail sur un exemple plus interessant en Section 2.3.2. En effet, comme évoqué plus haut, il s'agit ici d'un cas de bistabilité entre deux solutions musicalement équivalentes. De plus, nous pourrions montrer que même les bassins d'attraction dynamiques, qui seront introduits en Section 2.5, sont également triviaux.

Le diagramme de bifurcation de (2.50) est représenté en Fig. 2.5. Nous y voyons la solution d'équilibre triviale $x = 0$, qui correspond également à la solution triviale de (2.48), et les solutions d'équilibre non triviales (2.51), qui représentent l'amplitude des solutions périodiques de (2.48). De plus, des exemples d'intégration numérique des Éqs. (2.48) et (2.50) sont montrés en Fig. (2.6) : pour $\gamma = 0.2 < \hat{\gamma}^{\text{st}}$ en Fig. 2.6(a) et $\gamma = 0.6 > \hat{\gamma}^{\text{st}}$ en Fig. 2.6(b). Nous pouvons remarquer un bon accord entre la simulation de (2.48) et celle de (2.50).

Le diagramme représenté en Fig. 2.5, obtenu avec l'approximation polynomiale (2.43) de la ca-
ractéristique non linéaire, est représentatif de celui obtenu avec (2.23b) pour de faibles amplitudes
d'oscillations, c.-à-d. quand $\gamma$ ne dépasse pas trop $\hat{\gamma}^{\text{st}}$. Quand $\gamma$ s'éloigne de $\hat{\gamma}^{\text{st}}$ il faut donc utiliser
(2.23b). Dans ce cas, les hypothèses de la moyennisation restent valables (c.-à-d. supposer que le terme
$\epsilon_1 h(p, \dot{p}, \gamma)$ dans (2.48) est petit devant $\ddot{p} + p$) mais il n'est plus possible d'obtenir une expression
analytique de $f(x, \gamma)$ dans l'Éq. (2.50). Par conséquent, il faut déterminer directement les solutions
périodiques de (2.39) (par équilibrage harmonique ou continuation par exemple) puis en calculer la
<figure><img src="image_42.png" /><figcaption>**FIGURE 2.7.** Diagramme de bifurcation de l'Éq. (2.39) avec $F(p)$ donné par (2.23b) obtenu par continuation numérique des solutions périodiques et application de la théorie de Floquet à l'aide du logiciel AUTO [57]. Les solutions stables sont représentées en bleu et les solutions instables en rouge. Les points H de la bifurcation de Hopf et le point SN de la bifurcation nœud-col (saddle-node) des solutions périodiques sont représentés par des points bleus et point rouge, respectivement. La zone grisée représente la zone de bistabilité. Les paramètres du modèle sont les mêmes que ceux utilisés en Fig. 2.5. *Figure produite et fournie par Soizic Terrien*<fnref n="6" />.</figcaption></figure>

stabilité par la théorie de Floquet (voir par exemple [145, 179] pour les fondements théoriques de ces méthodes). Ceci est présenté dans la section suivante.

### 2.3.2 Analyse de stabilité globale : compétition entre régimes

Le diagramme de bifurcation de l'Éq. (2.39) avec $F(p)$ donné par (2.23b) obtenu par continuation numérique des solutions périodiques et application de la théorie de Floquet à l'aide du logiciel AUTO [57] est représenté en Fig 2.7. Pour cela la fonction $F(p)$ est régularisée en utilisant les égalités suivantes :

$$ \Theta(1 - \gamma + p) = \frac{1 - \gamma + p + |1 - \gamma + p|}{2(1 - \gamma + p)} \quad (2.53a) $$

$$ \operatorname{sgn}(\gamma - p) = \frac{\gamma - p}{|\gamma - p|} \quad (2.53b) $$

et l'approximation de la fonction racine carrée :

$$ |x| \approx \sqrt{x^2 + \eta} \quad (2.54) $$

avec $\eta$ un paramètre de régularisation, pris ici tel que $\eta = 10^{-3}$.

Le diagramme de la Fig 2.7 montre un domaine de bistabilité avec présence d'une solution périodique et de la solution d'équilibre (voir zone grisée entre $\gamma = 1$ et $\gamma \approx 1.3$ en Fig 2.7). Comme précédemment, il s'agit d'un cas où les SBA des solutions stables sont triviaux. En effet, la séparatrice entre les deux SBA est la solution périodique instable. Ceci est illustré en Fig. 2.8. La Fig. 2.8(a) montre, dans l'espace des phases $(p, \dot{p})$ pour les mêmes paramètres qu'en Fig. 2.7 avec $\gamma = 1.1$, les solutions périodiques (cercle bleu pour la solution stable et cercle rouge pour la solution instable) et la position d'équilibre (stable et représentée par un point bleu). Les étoiles jaune, verte et rouges représentent les conditions initiales (CI) de la simulation numérique de l'Éq. (2.39) dont les séries temporelles sont reportées en Fig. 2.8(b), haut, milieu et bas, respectivement. La simulation dont la CI est à l'intérieur du cercle rouge (Fig. 2.8(b) en haut) converge bien vers la position d'équilibre. Pour les CI à l'extérieur du cercle rouge les simulations convergent vers la solution périodique, que la CI soit à l'intérieur du cercle bleu (Fig. 2.8(b) au milieu) ou à l'extérieur (Fig. 2.8(b) en bas).

<footnote n="6">6. Laboratoire d'Acoustique de l'Université du Mans (LAUM), UMR 6612, Institut d'Acoustique - Graduate School (IA-GS), CNRS, Le Mans Université</footnote>
<figure><img src="image_43.png" /><figcaption>**FIGURE 2.8.** Illustration de la bistabilité de l'Éq. (2.39). (a) Représentation, dans l'espace des phases $(p, \dot{p})$ et pour les mêmes paramètres qu'en Fig. 2.7 avec $\gamma = 1.1$, des solutions de l'Éq. (2.39) : solutions périodiques (cercle bleu bleu pour la solution stable et cercle rouge la solution instable), position d'équilibre stable (un point bleu). Les étoiles jaune, verte et rouges représentent les conditions initiales (CI) de simulation numérique de l'Éq. (2.39) dont les séries temporelles sont reportées en (b). La simulation dont la CI est à l'intérieur du cercle rouge (étoile jaune et (b) en haut) converge vers la position d'équilibre. Pour les CI à l'extérieur du cercle rouge les simulations convergent vers la solution périodique que la CI soit à l'intérieur du cercle bleu (étoile verte et (b) au milieu) ou à l'extérieur (étoile rouge et (b) en bas).</figcaption></figure>

Qu'advient-il de la notion de bassin d'attraction lorsque le paramètre de contrôle varie lentement dans le temps, c.-à-d. lorsque le modèle passe de la forme (2.1) à la forme (2.3) ? Nous tenterons de donner quelques éléments de réponse à cette question en Section 2.5 en introduisant la notion de *bassin d'attraction dynamique*. D'abord, le phénomène de retard à la bifurcation (qui s'avère être aussi une forme de réponse à la question posée ci-dessous) est abordé dans la Section (2.4) qui suit.

# 2.4 RETARD À LA BIFURCATION

## 2.4.1 Synthèse des travaux de recherche publiés

Une vision naïve qui consiste à penser que pendant un transitoire d'attaque l'apparition des oscillations se produit au point de bifurcation statique est souvent erronée. En effet, on observe en général l'apparition des oscillations pour une valeur plus grande que le point de bifurcation statique, appelée *point de bifurcation dynamique*. C'est ce que l'on appelle le *retard à la bifurcation* [15]. Un état de l'art relativement complet sur le sujet, qui n'est pas propre aux instrument de musique, est fourni dans l'introduction de [54].

Dans le cadre de ma recherche doctorale, le retard à la bifurcation a été mis en évidence expérimentalement sur une clarinette de laboratoire à l'aide de mesures *in vitro* : l'instrument joue, inséré
<figure><img src="image_44.png" /><figcaption>**FIGURE 2.9.** (a) Bouche artificielle du LMA dans laquelle est insérée une clarinette de laboratoire. (b) et (c) Résultats de mesure sur bouche artificielle. L'expérience consiste à augmenter linéairement, avec un coefficient directeur *k*, la pression dans la bouche jusqu'à apparition des oscillations puis à la diminuer jusqu'à observer le silence. La figure (c) représente l'amplitude de la pression dans le bec *P* en fonction de la pression dans la bouche *P*<sub>m</sub> pour différentes valeurs de *k* avec *k*<sub>1</sub> < *k*<sub>2</sub> < *k*<sub>3</sub> < *k*<sub>4</sub> < *k*<sub>5</sub> < *k*<sub>6</sub>.</figcaption></figure>

dans une bouche artificielle. Cette dernière est un dispositif expérimental permettant de mettre l'ins-
trument en auto-oscillation sans avoir recours au musicien (voir Fig. 2.9(a) <fnref n="7" />). L'expérience consistait
à augmenter linéairement, avec un coefficient directeur *k*, la pression dans la bouche jusqu'à ap-
parition des oscillations puis à la diminuer jusqu'à observer le silence. Un exemple est présenté en
Fig. 2.9(b), où l'on voit en noir la pression dans la bouche et en gris la pression dans le bec. Les
résultats, publiés dans [a4], sont résumés en Fig. 2.9(c) qui représente l'amplitude de la pression
dans le bec en fonction de la pression dans la bouche artificielle pour différentes valeurs de *k* (avec
*k*₁ < *k*₂ < *k*₃ < *k*₄ < *k*₅ < *k*₆). Pour une valeur de *k* donnée, on observe en effet que l'apparition
des oscillations en montée se produit pour une valeur de pression dans la bouche (qui correspond au
point de bifurcation dynamique) plus grande que la valeur où le son s'éteint en descente. Cette der-
nière étant proche du point de bifurcation statique. Une deuxième observation est que ce retard à la
bifurcation observé expérimentalement augmente avec la valeur de *k* : plus on fait varier rapidement
la pression dans la bouche, plus le retard est important.

<footnote n="7">7. La bouche artificielle utilisée est celle développée au LMA dans le cadre du projet ANR CONSONNES.</footnote>
Ce travail expérimental a été complété par une étude théorique basée sur des modèles d'instrument extrêmement simples. Tout d'abord, pendant ma thèse, le modèle de Raman (cf. Section (2.2.4.1)), qui suppose une propagation acoustique avec pertes non dispersives dans le résonateur de l'instrument, a été utilisé. Dans ce cas, les équations du modèle se réduisent à une simple équation aux différences ou suite récurrente (système à temps discret) non linéaire (cf. Éq. (2.35)) qui peut être résolue simplement par un processus itératif. Les phénomène de retard à la bifurcation dans le modèle de Raman est étudié en détail dans mon manuscrit de thèse [t1] dont les chapitres concernés sont les retranscriptions entre autres de trois articles [a1, a5, a6]. Le second modèle utilisé est le modèle « à un mode », décrit en Section (2.2.4.2) qui prend la forme d'une équation différentielle à temps continu proche de celle de l'oscillateur de Van der Pol (cf. l'Éq. (2.39)).

Nous présentons dans cette section quelques éléments d'analyse du retard à la bifurcation dans le modèle à un mode, l'analyse détaillée est présentée dans [a15], reproduit dans l'Annexe G de ce manuscrit. Il est à noter que dans [a15], par rapport aux notations introduites en Section 2.2.4.2, $Z_1$ est noté $F_1/(\epsilon_1\omega_1)$ et $\epsilon_1$ est noté $\alpha_1$.

Pour des raisons de commodité de calcul, on pose

$$ \gamma = \tilde{\gamma} + \hat{\gamma}^{\text{st}}. \qquad (2.55) $$

Par conséquent, le modèle à un mode, écrit dans l'espace des phases (cf. Éq. (2.40)), prend la forme suivante :

$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \tilde{\gamma}), \qquad (2.56) $$

où, pour alléger les notation, $\mathbf{f}(\mathbf{x}, \tilde{\gamma} + \hat{\gamma}^{\text{st}})$ est noté $\mathbf{f}(\mathbf{x}, \tilde{\gamma})$. À noter que par le biais de (2.55) ; le point de bifurcation de Hopf de (2.56) est $\hat{\gamma}^{\text{st}} = 0$.

Dans le cas dynamique avec un croissance linéaire de $\gamma$ et donc également de $\tilde{\gamma}$, le modèle à un mode, prend la forme de l'Éq. (2.3) avec $\mathbf{x} = (x_1 = p, x_2 = \dot{p}) \in \mathbb{R}^2$ et $g(y) = 1$, soit

$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, y), \qquad (2.57a) $$

$$ \dot{y} = \epsilon \qquad (2.57b) $$

où $\gamma'$ est notée $y$ et considérée comme une autre variable d'état.

Avant d'aborder les aspects théoriques, le phénomène de retard à la bifurcation est illustré à nouveau en Fig. (2.10), cette fois à l'aide de simulations numériques du modèle à un mode (2.57). Le diagramme de bifurcation de Éq. (2.56) (c.-à-d. le diagramme statique), identique à celui de la Fig. 2.5, est représenté en Fig. 2.10(a). Les séries temporelles $p(t)$ et $y(t)$ obtenues par intégration numérique de l'Éq. (2.57) sont tracées en Fig. 2.10(b). La représentation paramétrique, $p(t)$ en fonction de $y(t)$ de cette même intégration numérique superposés au diagramme de bifurcation statique est montrée en Fig. 2.10(c) et met en évidence le phénomène de retard à la bifurcation. On observe en effet que les oscillations de (2.57) apparaissent pour une valeur de $y$ bien supérieure au point de bifurcation statique de (2.56).

Il est connu [22] que pour décrire correctement le retard à la bifurcation il faut sortir du cadre déterministe. Pour cela un forçage stochastique est ajouté à l'Éq. (2.48). Avec en plus le changement de paramètre (2.55), cela donne :

$$ \ddot{p} + \epsilon_1 h(p, \dot{p}, \tilde{\gamma}) + p = \nu \xi(t) \qquad (2.58) $$

où $h(p, \dot{p}, \tilde{\gamma} + \hat{\gamma}^{\text{st}})$ est noté simplement $h(p, \dot{p}, \tilde{\gamma})$ et $\xi(t)$ est un bruit blanc unitaire de moyenne nulle, c.-à-d.

$$ \mathbb{E}[\xi(t)] = 0 \quad \text{et} \quad \mathbb{E}[\xi(t)\xi(t + \tau)] = \delta(\tau) \qquad (2.59) $$

avec $\nu$ le niveau de bruit. $\mathbb{E}[X]$ est l'espérance mathématique de la variable aléatoire $X$ et $\delta$ la distribution de Dirac.

La méthode de *moyennisation stochastique* [114, 191] permet comme en Section 1.2.6 d'obtenir, à partir de (2.58), l'équivalent stochastique de l'Éq. (2.50). Le système obtenu, appelé *stochastic slow*
<figure><img src="image_45.png" /><figcaption>**FIGURE 2.10.** (a) Diagramme de bifurcation de Éq. (2.40) (c.-à-d. le diagramme statique). (b) Séries temporelles $p(t)$ (bleu) et $y(t)$ (vert) obtenues par intégration numérique de l'Éq. (2.57). Le phénomène de retard à la bifurcation apparaît en Fig.2.10(c) où l'on observe en effet les oscillations de (2.57) apparaître pour une valeur de y bien supérieure au point de bifurcation statique de (2.40).</figcaption></figure>

dynamics dans [a15], est décrit par le système (1, 1)-lent-rapide stochastique suivant 8 :

$$dx_t = F(x_t, y_t)dt + \sigma dW_t \qquad (2.60a)$$

$$dy_t = \epsilon \qquad (2.60b)$$

avec $W_t$ le processus de Wiener, relié au bruit blanc par $\xi_t = dW_t/dt$ et $\sigma = \nu/\sqrt{2}$. À noter que le problème initial, décrit par (2.57), d'une bifurcation de Hopf dynamique devient, après moyennisation, celui d'une bifurcation fourche dynamique.

En divisant (2.60a) par (2.60b) et en utilisant les propriétés de décalage temporel et d'échelle du processus de Wiener (voir la Section 3.1 de [163] et/ou le Chapitre 2 de [153]) il est possible d'utiliser $y$ comme variable temporelle. Dans ce cas l'Éq. (2.60) devient :

$$dx_y = \frac{1}{\epsilon} F(x_y, y)dy + \frac{\sigma}{\sqrt{\epsilon}} dW_y. \qquad (2.61)$$

Les détails sur l'obtention des Éqs. (2.60) et (2.61) sont donnés en Section G.3.

Dans l'Éq. (2.61), $x_y$ représente un processus stochastique qui décrit en fonction de $y$ (devenue la variable temporelle du problème), l'évolution de la variable aléatoire « amplitude » $x$. Jusqu'à maintenant nous avons parlé du point de bifurcation dynamique comme de la valeur de $y$ où émergent les oscillations. Nous allons maintenant en donner la définition mathématique. De façon générale cette définition sera (i) une statistique obtenue sur un ensemble de réalisations du processus $x_y$ et (ii) aura un caractère arbitraire, contrairement au point de bifurcation statique. En effet, la définition de « l'émergence des oscillations » a elle-même un caractère arbitraire. Il existe deux définitions principales du point de bifurcation dynamique (ou, de façon équivalente, du retard à la bifurcation). La première définition (voir par exemple [11] ou [198]) consiste à raisonner sur $x_y^2$ (c.-à-d. l'énergie du signal) en définissant le point de bifurcation dynamique comme la valeur $y$ pour laquelle $E[x_y^2]$ atteint un certain seuil, par exemple la valeur initiale $x_{y0}^2$ avec $y_0$ la valeur initiale de $y$. Cette définition est illustrée en Fig. 2.11(a) qui représente $x_y^2$ (en bleu) pour six réalisations de (2.61) 9 et la moyenne $E[x_y^2]$ calculée sur ces six réalisations. Une échelle logarithmique est utilisée pour l'axe des ordonnées. Le point de bifurcation dynamique (repéré par une ligne verticale pointillée rouge)

<footnote n="8">8. Il est de coutume en calcul stochastique d'indiquer la dépendance temporelle des variables en indice, cette convention est utilisée dans [a15] et dans la fin de cette section.</footnote>

<footnote n="9">9. L'intégration numérique de l'Éq. (2.61) est réalisée avec la fonction ItoProcess du logiciel Mathematica [216]. En pratique c'est l'Éq.. (2.60) qui est simulée, on utilise ensuite une représentation dans le plan $(y, x)$ ou $(y, x^2)$ selon le besoin.</footnote>
<figure><img src="image_46.png" /><figcaption>**FIGURE 2.11.** Illustration des définitions du point de bifurcation dynamique. Six réalisations du processus $x_y$ décrit par (2.61). Une échelle logarithmique est utilisée pour l'axe des ordonnées. (a) Première définition : le point de bifurcation dynamique (repéré par une ligne verticale pointillée rouge) correspond à l'abscisse du point d'intersection entre la courbe de $\mathbb{E}[x_y^2]$ et la droite horizontale d'équation $x_y^2 = x_{y0}^2$. (b) Seconde définition : le point de bifurcation dynamique est la valeur moyenne $\mathbb{E}[Y]$ (repérée par une ligne verticale pointillée verte) des six valeurs de Y obtenues (points verts).</figcaption></figure>

correspond à l'abscisse du point d'intersection entre la courbe de $\mathbb{E}[x_y^2]$ et la droite horizontale d'équation $x_y^2 = x_{y0}^2$. Le retard à la bifurcation (ou bifurcation delay) est défini comme la différence entre le point de bifurcation dynamique et le point de bifurcation statique, ici $\hat{y}^{\text{st}} = 0$ (repéré par une ligne verticale pointillée grise). La deuxième définition (voir par exemple [106, 132] ou [190]) utilise la notion de temps de premier passage [168], c.-à-d. le temps $T$ (ici le temps est représenté par la variable $y$, donc le « temps » de premier passage est noté $Y$) pour que le processus $x_y^2$ atteigne pour la première fois un certain seuil (ici $x_{y0}^2$ comme pour la première définition). Le point de bifurcation dynamique est ici défini comme la valeur moyenne $\mathbb{E}[Y]$ de ces temps de premier passage obtenus sur un ensemble de réalisations de (2.61). Cette seconde définition est illustrée en Fig. 2.11(b), qui montre les six mêmes réalisations de (2.61) que la Fig. 2.11(a). Dans ce cas, le point de bifurcation dynamique est la valeur moyenne (repérée par une ligne verticale pointillée verte) des six valeurs de $Y$ obtenues (points verts) et la définition du retard à la bifurcation est la même que pour la première définition. En général, ces deux définitions aboutissent à deux résultats légèrement différents mais proches.

Pour ces deux définitions, c'est la valeur moyenne qui est considérée, c.-à-d. $\mathbb{E}[x_y^2]$ et $\mathbb{E}[Y]$ pour la première et la deuxième définition, respectivement. Ce choix est également arbitraire, d'autres statistiques pourraient être choisies, les valeurs les plus probables par exemple seraient aussi pertinentes physiquement.

**REMARQUE 2.2 (CHOIX DE LA DÉFINITION DU POINT DE BIFURCATION DYNAMIQUE).** Dans ce manuscrit, c'est la première définition mentionnée plus haut qui est adoptée. Le point de bifurcation dynamique est donc la valeur de $y > y_0$, notée $y^{\text{dyn}}$, telle que :

$$ \sqrt{\mathbb{E}[x_y^2]} = x_{y0}. \qquad (2.62) $$

Dans le cas déterministe l'Éq. (2.62) devient simplement $\sqrt{x_y^2} = x_{y0}$.

Le point de bifurcation dynamique correspond à des petites valeurs de l'amplitude. Par conséquent, il est pertinent d'étudier la forme linéairisée de l'Éq. (2.61) autour de $x_y = 0$. En notant $a(y) =$
$\partial_x F(0, y)$, cette dernière s'écrit :

$$dx_y = \frac{1}{\epsilon}a(y)x_y dy + \frac{\sigma}{\sqrt{\epsilon}}dW_y \quad (2.63)$$

dont la solution est (voir Annexe G.4.1) :

$$x_y = x^{\det}(y) + \frac{\sigma}{\sqrt{\epsilon}} e^{\frac{1}{\epsilon} A(y)} \int_{y_0}^{y} e^{-\frac{1}{\epsilon} A(y')} dW_{y'} \quad (2.64)$$

où $x^{\det}(y) = x_{y_0} e^{\frac{1}{\epsilon}(A(y)-A(y_0))}$ est la solution du problème déterministe (c.-à-d. (2.63) avec $\sigma = 0$) et $A(y)$ la primitive de $a(y)$. À partir de (2.64) l'espérance de $x_y^2$ est obtenue comme suit :

$$\mathbb{E}[x_y^2] = \mathcal{D}(y) + \mathcal{S}(y) \quad (2.65)$$

où

$$\mathcal{D}(y) = (x^{\det}(y))^2 = x_{y_0}^2 e^{\frac{2}{\epsilon}(A(y)-A(y_0))} \quad (2.66a)$$

$$\mathcal{S}(y) = \frac{\sigma^2}{\epsilon} e^{\frac{2}{\epsilon} A(y)} \int_{y_0}^{y} e^{-\frac{2}{\epsilon} A(y)} dy. \quad (2.66b)$$

D'après (G.39) et selon [19], pour une équation différentielle de la forme de (2.61), on peut distinguer trois comportements ou régimes :

* **Regime I :** $S(y) \ll D(y)$. Dans ce cas, le niveau de bruit est si faible qu'il peut être ignoré et le problème est identique au cas déterministe (G.32) qui subit le plus grand retard à la bifurcation possible.

* **Regime II :** $S(y) \gg D(y)$ avec un niveau de bruit pas trop élevé. Le bruit ne peut plus être ignoré mais il reste suffisamment faible pour que le modèle linéarisé (G.30) reste valide. Le retard à la bifurcation existe toujours mais il est réduit par rapport au régime I.

* **Régime III :** $S(y) \gg D(y)$ avec un niveau de bruit élevé. Dans cette situation, le comportement du système est dominé par le bruit et la trajectoire de $x_y$ quitte le voisinage de zéro avant que le point de bifurcation statique $y = 0$ ne soit atteint. Dans un tel cas, l'approximation linéaire n'est plus valable et, comme le soulignent Berglung et Gentz [19], la notion de retard à la bifurcation n'a plus de sens.

Seuls les régimes I et II sont étudiés dans ce présent travail.

À partir de l'expression analytique de $\mathbb{E}[x_y^2]$ donnée en (2.65), il est possible d'obtenir :

(i) les domaines d'existence des régimes (cf. Section G.4.4, notamment la Fig. G.7). Ceci montre que la frontière entre les régimes II et III dépend surtout du niveau de bruit $\sigma$. En revanche, la frontière entre régimes I et II dépend du triplet $(\epsilon, \sigma, y_0)$.

(ii) en traitant les régimes I et II séparément, une expression théorique du point de bifurcation dynamique dans chacun de ces régimes est obtenue. Les détails sont donnés en Sections G.4.2 et G.4.3 pour les cas déterministe et stochastique, respectivement. Le point de bifurcation dynamique est noté $x^{\det}(y)$ dans le cas déterministe et $x^{\text{stoch}}(y)$ dans le cas stochastique (il y a en fait deux expressions $x^{\text{stoch},a}(y)$ et $x^{\text{stoch},b}(y)$ correspondant à deux niveaux d'approximation différents). Les résultats obtenus permettent de retrouver les caractéristiques bien connues du retard à la bifurcation. Dans le cas déterministe (Régime I) le retard à la bifurcation dépend de $y_0$ mais est indépendant de $\epsilon$. Ce comportement déterministe est illustré en Fig. 2.12 où l'Éq. (2.61) est simulée avec $\sigma = 0$ pour différentes valeurs de $y_0$ en fixant $\epsilon$ (cf. Fig 2.12(a)) et pour différentes valeurs de $\epsilon$ en fixant $y_0$ (cf. Fig 2.12(a)). Dans le cas stochastique (Régime II), le retard à la bifurcation dépend du niveau de bruit $\sigma$ et de la pente $\epsilon$ et ne dépend pas de
<figure><img src="image_47.png" /><figcaption>**FIGURE 2.12.** Caractéristiques du retard à la bifurcation dans le cas déterministe. Intégration numérique de l'Éq. (2.61) avec $\sigma = 0$ (a) pour différentes valeurs de $y_0$ en fixant $\epsilon$ et (b) pour différentes valeurs de $\epsilon$ en fixant $y_0$.</figcaption></figure>

la valeur initiale $y_0$. Ceci est montré en Fig. 2.13 qui montre $\sqrt{\mathbb{E}[x_y^2]}$ (calculée à chaque fois sur vingt réalisations de l'Éq. (2.61)) en fonction de $y$ pour différentes valeurs de $\sigma$ avec $\epsilon$ et $y_0$ fixés (cf. Fig. 2.13(a)); pour différentes valeurs de $\epsilon$ avec $\sigma$ et $y_0$ fixés (cf. Fig. 2.13(b)); et pour différentes valeurs de $y_0$ avec $\sigma$ et $\epsilon$ fixés (cf. Fig. 2.13(c)). Il interessant de remarquer que le comportement observé pour le Régime II en Fig 2.13(b) est qualitativement le même sur celui observé expérimentalement dans [a4] (cf. Fig. 2.9(c)).

(iii) la densité de probabilité $\rho(x, y)$ :

$$ \rho(x, y) = \frac{1}{\sqrt{2\pi S(y)}} e^{-\frac{(x-x_{\text{det}(y)})^2}{2S(y)}}. \quad (2.67) $$

telle que, pour une valeur donnée de $y$, la probabilité de trouver $x_y$ dans l'intervalle $[x, x + dx]$ est $\rho(x, y)dx$. La comparaison entre l'expression théorique (2.67) de $\rho(x, y)$ et la densité obtenue numériquement est présentée en Fig. G.11.

Ces travaux ont également fait l'objet de communications dans des congrès internationaux [ci1-ci5] et un congrès national [cn3].

## 2.4.2 Note sur la seconde définition du point de bifurcation dynamique

La probabilité $Q(y)$ que $x_y^2 \le x_{y0}^2$ est égale à la probabilité que $-x_{y0} \le x_y \le x_{y0}$, notée $\text{Pr}(-x_{y0} \le x_y \le x_{y0})$. Or, à partir de (2.67) on a par définition :

$$ \begin{aligned} \text{Pr}(-x_{y0} \le x_y \le x_{y0}) &= Q(y) \\ &= \int_{-x_{y0}}^{x_{y0}} \rho(x, y)dx \\ &= \frac{1}{2} \left[ \text{erf}\left(\frac{x_{y0} - x^{\text{det}}(y)}{2\sqrt{S(y)}}\right) + \text{erf}\left(\frac{x_{y0} + x^{\text{det}}(y)}{2\sqrt{S(y)}}\right) \right] \end{aligned} \quad (2.68) $$

où erf est la fonction erreur.

Si comme précédemment le temps de premier passage $Y$ est la valeur de $y$ quand $x_y^2$ atteint $x_{y0}^2$ pour la première fois, alors la probabilité que $y \le Y$, notée $\text{Pr}(y \le Y)$, vaut également $1 - P(y)$. La
<figure><img src="image_48.png" /><figcaption>**FIGURE 2.13.** Caractéristiques du retard à la bifurcation dans le cas stochastique. $\sqrt{\mathbb{E}[x_y^2]}$ (calculée à chaque fois sur vingt réalisations de l'Éq. (2.61)) en fonction de y pour (a) différentes valeurs de $\sigma$ avec $\epsilon$ et $y_0$ fixés; (b) pour différentes valeurs de $\epsilon$ avec $\sigma$ et $y_0$ fixés; et (c) pour différentes valeurs de $y_0$ avec $\sigma$ et $\epsilon$ fixés.</figcaption></figure>

probabilité que le temps de premier passage Y soit le temps d'observation y est donc $Pr(Y \le y) = 1 - Pr(y \le Y) = 1 - Q(y)$ noté $P(y)$. Par conséquent, la probabilité que Y se situe entre y et $y + dy$ est $\tilde{R}(y)dy = P(y + dy) - P(y)$. Si l'on utilise un développement limité au premier ordre, on obtient :

$$ \tilde{R}(y) = \frac{dP(y)}{dy} \qquad (2.69) $$

qui est la densité de probabilité du temps de premier passage.

Comme précédemment, les Régimes I et II sont étudiés séparément. Pour le Régime I, comme $D(y) = x^{\text{det}}(y)^2 \gg S(y)$, on a $\tilde{R}(y) \approx \delta(x_{y_0} - x^{\text{det}}(y))$. En introduisant $\ell(y) = x_{y_0} - x^{\text{det}}(y)$ dont les racines sont $y_0$ and $\hat{y}_{\text{det}}^{\text{dyn}}$ (cf. Eq. (G.41)) et en utilisant la propriété de composition de la distribution de Dirac<fnref n="10" />, il vient :

$$ \tilde{R}(y) \approx \tilde{R}_I(y) = \frac{\delta(y - y_0)}{|\ell'(y_0)|} + \frac{\delta(y - \hat{y}_{\text{det}}^{\text{dyn}})}{|\ell'(\hat{y}_{\text{det}}^{\text{dyn}})|} \qquad (2.70) $$

Dans le cas du Régime II, on a $D(y) \ll S(y)$ et donc l'Éq. (2.68) devient

$$ \Pr(-x_{y_0} \le x_y \le x_{y_0}) = 1 - P(y) = \erf\left(\frac{x_{y_0}}{2\sqrt{S(y)}}\right) \qquad (2.71) $$

$$ \delta(g(x)) = \sum_i \frac{\delta(x - x_i)}{|g'(x_i)|} $$

où les $x_i$ sont les racines (supposées simples) de la fonction g.

<footnote n="10">10. La propriété de composition de la distribution de Dirac correspond à l'identité suivante :</footnote>
<figure><img src="image_49.png" /><figcaption>**FIGURE 2.14.** Comparaison, dans un cas stochastique (c.-à-d. en Régime II) entre l'expression analytique (2.72) de la densité de probabilité du temps de premier passage (courbe bleue) et son estimation numérique obtenue à partir de 500 réalisations des processus (2.58) (Full order, histogramme orange) (2.63) (Averaged, histogramme bleu). Les valeurs du point de bifurcation dynamique obtenues avec la première définition (cf. Remarque (2.2)) et la seconde définition (cf. Éq. (2.73)) sont également représentées par des lignes verticales rouge et verte, respectivement. Les valeurs du niveau de bruit $\sigma$ sont (a) $\sigma = 10^{-5}$ et (b) $\sigma = 10^{-7}$.</figcaption></figure>

et

$$ \tilde{R}(y) \approx \tilde{R}_{II}(y) = \frac{dP(y)}{dy} = x_{y0} \frac{S'(y)}{\sqrt{2\pi} S(y)^{3/2}} e^{-\frac{x_{y0}^2}{2S(y)}}. \quad (2.72) $$

Selon la seconde définition évoquée plus haut, le point de bifurcation dynamique vaut :

$$ \mathbb{E}[Y] = \int_{-\infty}^{+\infty} y \tilde{R}(y) dy. \quad (2.73) $$

Dans le cas déterministe (2.73) redonne bien la valeur $\hat{y}_{det}^{dyn}$ et dans le cas stochastique il n'y a pas d'expression analytique simple de (2.73) qui est donc calculée numériquement.

Ces derniers résultats sont illustrés en Fig. 2.14 qui compare, dans un cas stochastique (c.-à-d. en Régime II), l'expression analytique (2.72) de la densité de probabilité du temps de premier passage à son estimation numérique obtenue à partir de 500 réalisations des processus (2.58) (appelé Full order) (2.63) (appelé Averaged). Les valeurs du point de bifurcation dynamique obtenues avec la première définition (cf. Remarque (2.2)) et la seconde définition (cf. Éq. (2.73)) sont également représentées par des lignes verticales rouge et verte, respectivement. Nous pouvons ainsi observer que ces deux définitions ne fournissent pas des valeurs égales mais proches. Le niveau de bruit $\sigma$ est plus important en Fig. 2.14(a) qu'en Fig. 2.14(b) où le retard est donc plus important comme cela a déjà été observé en Fig. 2.13(a).

## 2.5 BASSINS D'ATTRACTION DYNAMIQUES

Cette section reprend les résultats principaux de [p1] (reproduit en Annexe H) qui introduit la notion de *bassin d'attraction dynamique*.

### 2.5.1 Introduction

Les instruments de musique à anche sont connus pour produire une multitude de régimes sonores. Dans les modèles physiques, cette complexité se manifeste par la coexistence de régimes d'équilibre, périodiques ou quasi-périodiques et même chaotiques. Nous avons vu par exemple en Fig. 2.7 que le diagramme de bifurcation statique du modèle à un mode présentait un domaine de bistabilité entre
la solution d'équilibre et une solution périodique. Les diagrammes de bifurcation de modèles plus raffinés montrent en général aussi la présence de zones de multistabilité où coexistent des solutions périodiques stables (voir par la Fig. 3 de Colinot et coll. [43]). En pratique, du point de vue du jeu instrumental, la multistabilité peut être reliée à la facilité avec laquelle une note donnée peut être jouée. Par exemple, dans le cas d'un instrument à anche comme le saxophone, les musiciens savent que pour un doigté donné correspondant à une note grave, les registres d'ordre supérieurs (deuxième et troisième, cf. Définition 2.1) sont parfois plus faciles à produire que le premier registre. C'est également sur ce principe de multistabilité que fonctionnent les cuivres, comme les trompettes et les trombones.

En plus de repérer l'apparition du son, pendant un transitoire d'attaque il est important de savoir quel régime sera effectivement atteint lorsque différents régimes stables coexistent pour les mêmes paramètres. Pour investiguer cela, le modèle à un mode est à nouveau utilisé avec ici la « vraie » caractéristique non linéaire de l'instrument, c.-à-d. $u = F(p)$ donnée par l'Éq. (2.23b). Le diagramme de bifurcation statique de ce modèle est donc représenté en Fig. 2.7.

En général, la multistabilité peut être reliée au phénomène de transition critique (critical transition), c'est-à-dire un changement soudain et de grande amplitude de la sortie d'un système dynamique causé par un faible changement de l'entrée. L'article de revue [63] examine les résultats récents sur le lien entre multistabilité et transition critique et présente quelques applications en physique, neurosciences, sciences du climat, biologie et écologie. Également appelées points de basculement (tipping points) en sciences du climat, les transitions critiques peuvent se produire de différentes manières. Ashwin et coll [7] proposent une classification en trois catégories basée sur la description mathématique des mécanismes conduisant à leur apparition. Les premières sont les transitions induites par le bruit<fnref n="11" /> (TIBr). En effet, les systèmes multistables sont connus pour être potentiellement très sensibles au bruit qui peut les faire basculer d'un état stable à un autre. Les TIBr n'impliquent pas forcément de variations dans le temps d'un paramètre de bifurcation du système étudié et ne sont pas spécifiquement étudiées dans ce manuscrit. Il existe ensuite les transitions induites par une bifurcation<fnref n="12" /> (TIBi) qui ne sont autres que les bifurcations dynamiques évoquées en Section 2.4. Enfin, il y a les transitions induites pas la variation dans le temps d'un paramètre de bifurcation du système<fnref n="13" /> (TIV). Les TIBi et les TIV considèrent spécifiquement des systèmes dynamiques de la forme de (2.3).

De façon générale, les solutions d'un système de la forme de l'Éq. (2.3) (à la fois transitoires et permanentes) ne sont plus solutions de l'Éq. (2.1) puisque $y$ n'est plus constant dans le temps. Cependant, les solutions permanentes (solutions d'équilibre, périodiques ou quasi-périodiques) stables de l'Éq. (2.1) sont associées à des variétés<fnref n="14" /> invariantes localement attractives de l'Éq. (2.3) qui évoluent dans l'espace des phases de dimension ($n + 1$) de ce dernier et attirent les orbites (ou trajectoires) qui se trouvent dans leur voisinage. Dans le cas des TIBi, on étudie le comportement de (2.3) quand la variable lente $y$ croise un point de bifurcation $\hat{y}^{st}$ de l'Éq. (H.1). En général, la transition vers une autre variété attractive ne se produit pas à $\hat{y}^{st}$ mais à une valeur plus grande, c'est le phénomène de retard à la bifurcation décrit en Section 2.4. Les TIV, identifiées par Wieczorek et coll. [212], sont des transitions entre variétés attractives de l'Éq. (2.3) qui se produisent sans la présence de bruit ou de la traversée de point de bifurcation.

Les travaux cités dans le paragraphe précédent (voir l'introduction de [p1] en Annexe H.1 pour davantage de références) se focalisent sur les conditions provoquant la sortie d'une orbite de (2.3) du voisinage d'une variété qu'elle longeait. Ces travaux portent donc sur le début de la transition critique. Cependant, la question de la fin de la transition, lorsque plusieurs autres variétés attractives peuvent être potentiellement atteintes par l'orbite, est peu abordée dans la littérature. Dans ce contexte, la question que nous nous posions en fin de Sect. 2.3.2 peut être reformulée ainsi : que se passe-t-il pour une orbite de l'Éq. (2.3) lorsqu'elle quitte le voisinage d'une certaine variété attractive (ou une variété

<footnote n="11">11. Traduction de *noise-induced transitions*.</footnote>

<footnote n="12">12. Traduction de *bifurcation-induced transitions*.</footnote>

<footnote n="13">13. Traduction de *rate-induced transitions*.</footnote>

<footnote n="14">14. Intuitivement, une variété peut être vue comme la généralisation de la classification qui établit qu'une courbe est une variété de dimension 1 et une surface est une variété de dimension 2 (voir par exemple [69] et [40]).</footnote>
qui a été localement attractive) et que plusieurs variétés attractive sont présentes ?

Nous fournissons des éléments de réponse à ces questions en introduisant la notion de bassin d'attraction dynamique que nous illustrons d'abord à l'aide du modèle à un mode (cas d'une TIBi) en Section 2.5.2.1, qui synthétise les résultats de [p1]. Ensuite afin de souligner le caractère général de la méthodologie, son application à un cas d'une TIV est présenté en Section 2.5.3.

## 2.5.2 Bassin d'attraction dynamique du modèle à un mode

Pour investiguer cette thématique complexe sur un modèle académique d'instrument sujet à la multistabilité, le modèle à un mode (2.39) est à nouveau utilisé avec ici la « vraie » caractéristique non linéaire de l'instrument, c.-à-d. $u = F(p)$ (cf. Éq. (2.23b)). Le diagramme de bifurcation statique de ce modèle est donc représenté en Fig. 2.7. Dans les faits, nous utilisons une approximation polynomiale de la fonction $F(p)$ afin de pouvoir utiliser la méthode de moyennisation évoquée en Section 2.3.1. Par conséquent, le modèle obtenu a une forme similaire à (2.1) avec $x \in \mathbb{R}$, c.-à-d. :

$$ \dot{x} = \epsilon_1 f(x, \gamma) \qquad (2.74) $$

où $x$ est à nouveau l'amplitude des oscillations de la pression acoustique $p$ dans le bec de l'instrument. La dynamique de la pression de souffle (maintenant notée $y$) est ajoutée comme suit : $\dot{y} = \hat{e}g(y)$. Cela mène à :

$$ \dot{x} = \epsilon_1 f(x, y) \qquad (2.75a) $$

$$ \dot{y} = \hat{e}g(y). \qquad (2.75b) $$

Le changement de temps $t \to \epsilon_1 t$ conduit finalement à l'équation suivante :

$$ \dot{x} = f(x, y) \qquad (2.76a) $$

$$ \dot{y} = \epsilon g(y). \qquad (2.76b) $$

avec

$$ \epsilon = \frac{\hat{e}}{\epsilon_1} \qquad (2.77) $$

qui est le petit paramètre pertinent à utiliser ici.

L'Éq. (2.76) peut aussi être écrit par rapport au temps lent $\tau = \epsilon t$ comme suit

$$ \epsilon x' = f(x, y) \qquad (2.78a) $$

$$ y' = g(y), \qquad (2.78b) $$

où, comme au Chapitre 1, $(.)' = \frac{d(.)}{d\tau}$.

Les étapes permettant d'obtenir (2.76) sont détaillées en Annexe H.8.2 (attention encore aux quelques inhomogénéités de notations entre le corps du texte et cette annexe).

Le diagramme de bifurcation statique du problème (c.-à-d. celui de l'Éq. (2.74)) est représenté en Fig. 2.15. On observe une allure proche de celle du diagramme de la Fig. 2.7, avec un domaine de bistabilité entre l'équilibre trivial et un équilibre non trivial.

### 2.5.2.1 Cas d'une d'une pression de souffle qui croît dans le domaine de bistabilité

Nous considérons ici une augmentation linéaire de la pression de souffle, c.-à-d. $g(y) = 1$, et $y$ croît dans le domaine de bistabilité. La problématique est illustrée en Fig. 2.16 à l'aide de simulations numériques de (2.78) pour deux conditions initiales légèrement différentes (notées $N_1$ et $N_2$). Les séries temporelles obtenues sont représentées en Figs. 2.16(a) et 2.16(b) avec en bleu $y(t)$ et en vert $x(t)$. Dans le premier cas, on observe $x(t)$ partir de la condition initiale puis atteindre rapidement le voisinage de l'axe des abscisses et y rester jusqu'à la fin de la simulation. Dans le second cas, le début de la simulation est identique, mais une première transition rapide de $x(t)$ du voisinage de zéro vers
<figure><img src="image_50.png" /><figcaption>**FIGURE 2.15.** Diagramme de bifurcation de l'Éq. (2.74). Une zone de bistabilité entre l'équilibre trivial et un équilibre non trivial est présente entre $\gamma = 1.11$ and $\gamma = 1.34$ (zone bleutée).</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>**FIGURE 2.16.** Simulations numériques de l'Éq. (2.78) pour deux conditions initiales légèrement différentes (notées $N_1$ et $N_2$). Les séries temporelles obtenues sont représentées en (a) et (b) avec $y(t)$ en bleu et $x(t)$ en vert. (c) Superposition de ces mêmes trajectoires (en bleu pour celle initialisée en $N_1$ et en vert celle initialisée en $N_2$), représentées dans le plan $(y, x)$, et du diagramme de bifurcation de l'Éq. (2.74).</figcaption></figure>

une valeur non nulle se produit ; $x(t)$ évolue ensuite lentement puis une seconde transition rapide se produit où $x(t)$ redescend proche de zéro et y reste jusqu'à la fin de la simulation. La superposition de ces mêmes trajectoires représentées dans le plan $(y, x)$ avec le diagramme de bifurcation statique (voir Fig. 2.16(c)) permet de mieux comprendre les comportements observés. En effet, au départ des deux simulations, les orbites rejoignent l'unique branche stable du diagramme correspondant à l'équilibre trivial. En effet, cette branche du diagramme de bifurcation de (2.74) est également une variété invariante de (2.78) (de dimension 1) : attractive quand $y < \hat{\gamma}^{\text{st}}$, répulsive quand $y > \hat{\gamma}^{\text{st}}$ et attractive à nouveau pour $\gamma > 1.12$. Un phénomène de retard à la bifurcation est ensuite observé : les deux orbites longent pendant un certain temps la variété devenue répulsive après avoir dépassé le point de bifurcation statique. L'orbite débutant en $N_1$ reste sur la variété répulsive jusqu'a ce qu'elle redevienne attractive vers $y = 1.11$. Quant à l'orbite débutant en $N_2$, elle saute de la variété répulsive pour rejoindre l'autre variété attractive de (2.78) également de dimension 1 et « associée » à la branche de l'équilibre stable non trivial du diagramme de bifurcation de (2.74) ; une fois le domaine de bistabilité franchi l'orbite redescend sur la droite $x = 0$, la seule variété invariante disponible.

Les variétés invariantes attractives de (2.78) associées aux branches stables du diagramme de bifurcation de (2.74) mentionnées plus haut peuvent être définies rigoureusement dans le cadre de
l'approche géométrique des perturbations singulières [62, 107] déjà évoquée en Section 1.2 et dont les éléments importants sont donnés en Annexe H.4.

Ainsi, l'Éq. (2.74) (cas statique) peut être vue comme le sous-système rapide associé à l'Éq. (2.78), soit,

$$ \dot{x} = f(x, y), \tag{2.79a} $$

$$ \dot{y} = 0. \tag{2.79b} $$

Le sous-système lent s'écrit :

$$ 0 = f(x, y) \tag{2.80a} $$

$$ y' = 1. \tag{2.80b} $$

La résolution de l'Éq. (2.80a) permet de définir la variété critique $M_0$ de l'Éq. (2.78) telle que

$$ M_0 = \{ (x, y) \in \mathbb{R}^2 \mid x = x^*(y) \} \tag{2.81} $$

avec $x^*(y)$ les branches de solution de

$$ f(x^*(y), y) = 0 \tag{2.82} $$

dans chaque intervalle de $y$ où $\partial_x f(x, y)$ ne s'annule pas et donc où l'on peut écrire que $f(x, y) = 0$ est équivalent à $x = x^*(y)$ par le théorème des fonctions implicites (cf. Théorème 2.1).

**THÉORÈME 2.1 (DES FONCTIONS IMPLICITES, VERSION $\mathbb{R}^2$).** Soit $f$ une fonction de classe $C^k$ définie sur un ouvert $\Omega \subset \mathbb{R}^2$, à valeurs dans $\mathbb{R}$. Soit $(a, b) \in \Omega$ tel que $f(a, b) = 0$. On suppose que la dérivée partielle de $f$ par rapport à la seconde variable est non nulle en $(a, b)$ :

$$ \partial_y f(a, b) \neq 0. $$

Alors il existe un voisinage ouvert $U$ de $(a, b)$ dans $\mathbb{R}^2$, un intervalle ouvert $I$ de $\mathbb{R}$ contenant $a$ et une fonction $\phi : I \to \mathbb{R}$ de classe $C^k$ telle que, pour tout $(x, y) \in U$ on ait :

$$ f(x, y) = 0 \Longleftrightarrow y = \phi(x). $$

$$ \text{De plus, pour tout } x \in I, \text{ on a } \phi'(x) = -\frac{\partial_x f(x, \phi(x))}{\partial_y f(x, \phi(x))}. $$

Chaque point de $M_0$ étant une solution d'équilibre du sous-système rapide (2.79), la variété critique n'est autre que ce nous appelions jusqu'à maintenant le diagramme de bifurcation statique du problème. Notons que la variété critique est indépendante de la manière dont $y$ varie dans le temps et donc indépendante de la fonction $g(y)$ du cas général (2.76).

La Définition H.2, réécrite ci-dessous en français, précise la définition d'une variété critique normalement hyperbolique.

**DÉFINITION 2.4.** Soit $a^*(y)$ la linéarisation de (2.79a) (de dimension 1) sur $M_0$, c.-à-d. en $x = x^*(y)$, définie comme suit :

$$ a^*(y) = \partial_x f(x^*(y), y) \tag{2.83} $$

Une valeur $x^*(y)$ de la variable rapide $x$ est une solution d'équilibre hyperbolique de (2.79) si $a^*(y) \neq 0$; stable si $a^*(y) < 0$ et instable si $a^*(y) > 0$. La variété critique est attractive (resp. répulsive) si $a^*(y) < 0$ (resp. $a^*(y) > 0$) pour $y \in I$ avec $I$ un intervalle de $\mathbb{R}$. Un sous-ensemble $M_0^{nh}$ de $M_0$ est normalement hyperbolique si pour chaque point $P = (x^*(y), y) \in M_0^{nh}$ on a $a^*(y) \neq 0$, c.-à-d. $M_0^{nh}$ est soit attractive, soit répulsive $^a$.

<footnote>a. Normalement signifie ici que chaque point de $M_0$ doit être hyperbolique uniquement dans la direction normale (c'est-à-dire non tangente) à elle-même (voir par exemple la définition 2.3.4 et le texte en dessous de [40]).</footnote>
Le Théorème de Fenichel (cf. Theorem H.1) stipule que si la variété critique $M_0$ est attractive (resp. répulsive), alors il existe une variété invariante $M_\epsilon$ qui se trouve à une distance de l'ordre de $O(\epsilon)$ de $M_0$. Cette variété invariante $M_\epsilon$ attire (resp. repousse) les orbites de (2.76) (ou (2.78)) qui sont dans son voisinage, c.-à-d. que $M_\epsilon$ est également attractive (resp. répulsive).

Sur la base des éléments de l'approche géométrique des perturbations singulières fournis en An-
nexe H.4 et pour certains rappelés ici, la définition de ce que nous appelons un *bassin d'attraction dynamique* va pouvoir être formalisée. Avant cela, la définition du domaine de bistabilité est donnée ci-dessous.

**DÉFINITION 2.5 (DOMAINE DE BISTABILITÉ DE $M_0$).** Un intervalle

$$D =]y_1, y_u[ \tag{2.84}$$

de $\mathbb{R}$ est appelé *domaine de bistabilité* de la variété critique, si pour $y \in D$, l'Éq. (2.82) possède deux branches attractives $x_1^*(y)$ et $x_2^*(y)$, c.-à-d. $a_1^*(y) = \partial_x f(x_1^*(y), y) < 0$ et $a_2^*(y) = \partial_x f(x_2^*(y), y) < 0$. Nécessairement, l'Éq. (2.82) possède aussi une branche répulsive $x_3^*(y)$ pour laquelle $a_3^*(y) = \partial_x f(x_3^*(y), y) > 0$. La variété critique n'est pas normalement hyperbolique en $y = y_1$ et $y = y_u$.

Comme énoncé plus haut, la variété critique $M_0$ représente le diagramme de bifurcation statique du sous-système rapide. Par conséquent, le domaine de bistabilité de $M_0$ est identique au domaine de bistabilité du problème statique montré en Fig. 2.15.

Il est maintenant possible de donner la définition d'un bassin d'attraction dynamique. Notez que les définitions qui suivent sont valables pour les systèmes (1, 1)-rapides-lents de la forme de l'équation (2.78) et dont la variété critique $M_0$ possède un domaine de bistabilité.

**DÉFINITION 2.6 (BASSIN D'ATTRACTION DYNAMIQUE.).** Dans le domaine de bistabilité $D$ (cf. Définition 2.5) la variété critique $M_0$ de (H.4) possède deux branches attractives $M_{0,a_1}$ et $M_{0,a_2}$ telles que

$$M_{0,a_i} = \{x \in \mathbb{R}, y \in D \mid x = x_i^*(y)\}, \quad i = 1, 2 \tag{2.85}$$

et séparée par une branche répulsive $M_{0,r}$

$$M_{0,r} = \{x \in \mathbb{R}, y \in D \mid x = x_3^*(y)\} \tag{2.86}$$

Par le théorème de Fenichel Fenichel H.1, l'Éq.(2.78) possède deux variétés invariantes localement attractives $M_{\epsilon,a_i}$ dans un voisinage $O(\epsilon)$ de $M_{0,a_i}$ ($i = 1, 2$) et une variété invariante localement répulsive $M_{\epsilon,r}$ dans un voisinage $O(\epsilon)$ de $M_{0,r}$.

Le *bassin d'attraction dynamique* (DBA pour *dynamic basin of attraction*) d'une branche attractive donnée $M_{\epsilon,a_i}$ ($i = 1, 2$) est le sous-ensemble de l'espace des phases pour lequel les orbites provenant de conditions initiales dans le DBA finissent par longer $M_{\epsilon,a_i}$ ($i = 1, 2$) lorsque le domaine de bistabilité $D$ est traversé.

Nous écrivons dans la Définition 2.6 « ... finissent par longer ... » car des phénomènes transitoires peuvent exister avant que les branches $M_{\epsilon,a_i}$ ($i = 1, 2$) ne soient atteintes.

La séparatrice entre deux DBA, ainsi que sa nature, sont enfin définies.

**DÉFINITION 2.7 (SÉPARATRICE ENTRE DEUX DBA).** La séparatrice entre deux DBA (SDBA pour *separatrix between DBA*) est la frontière dans l'espace des phases qui sépare deux DBA.
<figure><img src="image_52.png" /><figcaption>**FIGURE 2.17.** Simulations numériques de l'Éq. (2.78) pour deux jeux de conditions initiales légèrement différents (notées $N_1$ et $N_2$) et de l'Éq. (2.87) initialisé en $N_3$ défini par (2.88). (a) Superposition des séries temporelles de l'Éq. (2.78) (en bleue pour celle initialisée en $N_1$, en vert pour celle initialisée $N_2$) et de l'Éq. (2.87) (c.-à-d. la SDBA, en rouge), représentées dans le plan $(y, x)$, et de la variété critique $M_0$ (2.81) (qui n'est autre que la diagramme de bifurcation de l'Éq. (2.74)). (b) Mêmes séries temporelles qu'en (a) avec une échelle logarithmique pour l'axe des ordonnées. $\epsilon = 0.0183$.</figcaption></figure>

**DÉFINITION 2.8 (NATURE DE LA SDBA).** La SDBA est la variété invariante répulsive $M_{\epsilon,r}$ considérée au-delà de la limite inférieure $y_1$ du domaine de bistabilité D jusqu'à un point de l'espace des phases ayant la plus petite valeur possible de la variable lente y.

De façon pratique, $M_{\epsilon,r}$ n'est autre qu'une orbite particulière de l'Éq. (2.76) (ou (2.78)). À la manière de Berglund et Landon [26], cette solution particulière est ici approchée numériquement au moyen d'une procédure de renversement temporel, c.-à-d. par une intégration temporelle numérique de

$$ \epsilon x' = -f(x, y) \qquad (2.87a) $$

$$ y' = -1 \qquad (2.87b) $$

en choisissant une condition initiale $(x_0, y_0)$ appropriée, c.-à-d. sur la borne supérieure du domaine de bistabilité D et sur $M_{0,r}$, c.-à-d. sur le point $N_3$ tel que :

$$ N_3 = (x_3^*(y_u), y_u). \qquad (2.88) $$

Nous utilisons le fait que, selon le Théorème de Fenichel H.1, $M_{0,r}$ se trouve dans un voisinage de l'ordre de $O(\epsilon)$ de $M_{\epsilon,r}$ et que cette dernière est attractive pour les orbite de (2.87). L'intégration numérique est ensuite arrêtée à la valeur requise de la variable lente y. Cette méthode de renversement temporel a l'avantage de pouvoir être facilement généralisée à un espace des phases de dimension supérieure, même si dans ce cas des précautions sont à prendre (voir Section 2.5.2.3).

La méthodologie décrite ci-dessus est illustrée en Fig. 2.17. La Fig. 2.17(a) montre le même résultat de simulation numérique de l'Éq. (2.78) que celui montré en Fig. 2.16(c). Nous y voyons en plus l'orbite du système en temps inversé (2.87) initialisé en $N_3$ selon (2.88) et des flèches ont été rajoutées sur les courbes pour montrer le sens du temps. Cette orbite représente la SDBA entre les DBA des variétés invariantes $M_{\epsilon,a_1}$ et $M_{\epsilon,a_2}$, notés DBA$_1$ et DBA$_2$, respectivement. Au vu de ces résultats, la Fig. 2.16(c) peut maintenant s'interpréter en observant que, même si $N_1$ et $N_2$ sont très proches dans le plan $(y, x)$, on a $N_1 \in$ DBA$_1$ et $N_2 \in$ DBA$_2$. C'est la raison pour laquelle l'orbite initialisée en $N_2$ finit par longer $M_{\epsilon,a_2}$ lorsque le domaine de bistabilité D est traversé alors que celle initialisée en $N_1$ longe $M_{\epsilon,a_1}$ (qui correspond à $x = 0$) pendant toute la simulation. Les lignes de courant de (2.78)
sont également représentées (lignes grises fléchées). Ces dernières illustrent effectivement que, pour $y < y_u$, les orbites commençant avec une condition initiale dans DBA₁ (resp. DBA₂), c'est-à-dire en-dessous (resp. au-dessus) de la SDBA, finissent par longer $M_{\epsilon,a_1}$ (resp. $M_{\epsilon,a_2}$). Quand $y$ dépasse le domaine de bistabilité, toutes ces orbites rejoignent $M_{\epsilon,a_1}$, la seule variété attractive restante.

Pour mieux comprendre ce qui se passe au voisinage de $x = 0$, une échelle logarithmique pour les ordonnées est utilisée en Fig. 2.17(b). Nous pouvons voir que de très petites amplitudes sont rencontrées. Il est difficile d’envisager cette situation dans un problème réel où un bruit résiduel empêcherait d’atteindre des amplitudes aussi faibles. Ceci sera évoqué dans les perspectives en Section 4.2.2.

### 2.5.2.2 Cas d'une d'une pression de souffle qui converge dans le domaine de bistabilité

En situation de jeux lors de l'attaque d'une note, le musicien ne va pas faire croître indéfiniment la pression dans sa bouche mais plutôt l'augmenter jusqu'à une valeur cible (ceci n'est qu'une stratégie de contrôle parmi d'autres). Si cette valeur cible est supérieure à la borne supérieure du domaine de bistabilité, l'étude présentée dans les sections précédentes reste pertinente. Dans cette section nous étudions le cas où la valeur cible se trouve dans le domaine de bistabilité. Pour cela, un exemple impliquant un profil exponentiel de la variables lente $y$ convergeant vers une valeur $y_{\text{targ}}$ située dans le domaine de bistabilité $D$, c.-à-d. $g(y) = y_{\text{targ}} - y$ (avec $y_{\text{targ}} \in D$) est considéré. À l'échelle de temps lente $\tau$ l'équation du modèle à un mode devient donc :

$$ \epsilon x' = f(x, y) \qquad (2.89a) $$

$$ y' = y_{\text{targ}} - y. \qquad (2.89b) $$

Cette situation diffère du cas étudié jusqu'à présent, dans lequel la pression de souffle ne sature pas dans le domaine de bistabilité. En effet, la recherche des points fixes de l'Éq. (2.89), c.-à-d. la résolution de $f(x, y) = 0$ et $y_{\text{targ}} - y = 0$ conduit à $f(x, y_{\text{targ}}) = 0$ et donc $x = x_i^*(y_{\text{targ}})$ ($i = 1, 2, 3$). Par conséquent, les points $P_i^* = (x_i^*(y_{\text{targ}}), y_{\text{targ}})$ ($i = 1, 2, 3$) sont les points fixes de l'Éq. (2.89). La matrice jacobienne $\mathbf{J}(P_i^*)$ de (2.89) évaluée en chaque point fixe $P_i^*$ s'écrit :

$$ \mathbf{J}(P_i^*) = \begin{pmatrix} \frac{a_i^*(y_{\text{targ}})}{\epsilon} & \frac{\partial_y f(x_i^*(y_{\text{targ}}), y_{\text{targ}})}{\epsilon} \\ \frac{0}{\epsilon} & -1 \end{pmatrix} \qquad (2.90) $$

dont les valeurs propres sont $\mu_1 = \frac{a_1^*(y_{\text{targ}})}{\epsilon}$ et $\mu_2 = -1$. Par conséquent, comme $a_1^*(y_{\text{targ}}) < 0$, $a_2^*(y_{\text{targ}}) < 0$ et $a_3^*(y_{\text{targ}}) > 0$ (voir Définition 2.5), $P_1^*$ et $P_2^*$ sont des noeuds (nodes) - c.-à-d. $\mu_1\mu_2 > 0$ - et $P_3^*$ est un point-selle (saddle) - c.-à-d. $\mu_1\mu_2 < 0$ - . Ces définitions classiques des nœuds et des points-selles sont tirées de [179].

Le problème devient donc un problème classique de recherche des SBA des points fixes stables d'un système de dimension 2. Dans ce cas, les DBA des branches attractives $M_{\epsilon,a_i}$ ($i = 1, 2$) et les BAS des points fixes $P_1^*$ et $P_2^*$ sont identiques. Cependant, nous choisissons de conserver le vocabulaire DBA et SDBA car le problème (2.89) reste un problème dynamique tel que défini en Section 2.1. La séparatrice entre DBA est la variété stable du point fixe instable $P_3^*$ [38] qui est aussi une variété invariante localement répulsive $M_{\epsilon,r}$ de (2.89). Cette dernière est à nouveau calculée à l'aide d'une simulation en temps inversé. L'équation

$$ \epsilon x' = -f(x, y) \qquad (2.91a) $$

$$ y' = -(y_{\text{targ}} - y) \qquad (2.91b) $$

est donc intégrée numériquement avec une condition initiale choisie ici comme une petite perturbation de $P_3^*$.

Un exemple est donné en Fig. 2.18 dont le mode de représentation est similaire à celui de la Fig. 2.17(a). La valeur cible $y_{\text{targ}}$ de la pression dans la bouche est mise en évidence par une ligne verticale orange d'équation $y = y_{\text{targ}}$. Les points fixes de l'Éq. (2.89) apparaissent donc comme les
<figure><img src="image_53.png" /><figcaption>**FIGURE 2.18.** Figure similaire à la Fig. 2.17(a). En rouge, une orbite obtenue par l'intégration numérique du système a temps inversé (2.87) avec $g(y) = y_{\text{targ}} - y$ et une condition initiale (CI) choisie comme une petite perturbation du point fixe non trivial $P_3^*$. En bleu et vert, deux orbites obtenues par l'intégration de l'Éq. (2.89) : avec une CI d'abord dans le DBA$_1$ (en bleu) et dans le DBA$_2$ (en vert). Les lignes de courant de l'Éq. (2.89) (lignes fléchées grises) sont également représentées.</figcaption></figure>

points d'intersection entre cette dernière et la variété critique $M_0$ qui pour rappel est indépendante de la fonction $g(y)$ et donc identique à celle obtenue précédemment en Section 2.5.2.1. Un zoom au niveau du point fixe instable $P_3^*$ est effectué, on peut observer la forme typique des lignes de courant d'un système dynamique au niveau d'un de ses points-selles. Ici, c'est donc la variété stable de $P_3^*$ qui représente la SDBA.

### 2.5.2.3 Quid de l'extension de la méthode à des modèles de dimension supérieure ?

Les concepts de DBA et de SDBA ont été introduits pour des systèmes (1, 1)-rapides-lents. Cependant, les modèles plus réalistes d'instruments à anche sont de dimension supérieure. Travailler dans le plan a permis d'illustrer facilement ces concepts et donc de les comprendre plus aisément. Les concepts de DBA et de SDBA restent cependant valables pour les systèmes $(m, n)$-rapides-lents. Cependant, il y a deux choses auxquelles il faut faire attention lorsque le sous-système rapide est de dimension supérieure à 1. La première est que la variété critique peut être *saddle-type*, c.-à-d. ni attractive ni répulsive (dans ce cas, la matrice jacobienne du sous-système rapide évaluée sur la variété critique a des valeurs propres dont certaines ont des parties réelles positives et d'autres néga-tives, cf. Définition 1.2). Si la SDBA doit être calculée à partir d'une variété invariante *saddle-type* (au lieu d'une variété répulsive, comme c'était le cas dans les sections précédentes), la procédure de renversement temporel ne peut plus être utilisée. En effet, une variété *saddle-type* en temps direct reste *saddle-type* en temps inversé. La deuxième chose à considérer est que le sous-système rapide peut avoir d'autres types de solutions permanentes que les solutions d'équilibre, comme des solutions périodiques ou même quasi-périodiques. Ainsi, le diagramme de bifurcation complet du sous-système rapide (incluant les solutions périodiques et si possible les solutions quasi-périodiques) n'est plus la variété critique du système $(m, n)$-rapide-lent étudié qui n'inclut que les solutions d'équilibre. Les variétés invariantes d'un système $(m, n)$-rapide-lent associées à des solutions quasi-périodiques du sous-système rapide correspondant semblent très compliquées à étudier. En revanche, les solutions périodiques peuvent elles être traitées de la même manière que les solutions d'équilibre. Tout d'abord, parce que Berglund [22] propose un analogue du théorème de Fenichel quant à l'existence de variétés invariantes qui longent des familles d'orbites périodiques$^{15}$ et aussi parce que les solutions périodiques

<footnote n="15">15. La démonstration repose sur une procédure de moyennisation.</footnote>
diques sont faciles à calculer numériquement. Comme dans le cas des solutions d'équilibre, pour que la méthode du retournement temporel puisse être utilisée, la variété invariante considérée pour le calcul de la SDBA doit être répulsive. Pour les solutions périodiques, cela signifie que $m - 1$<fnref n="16" /> des multiplicateurs de Floquet associés à la solution périodique considérée ont un module strictement supérieur à un (voir par exemple [179], Chap. 7, concernant les détails sur la stabilité locale des solutions périodiques).

En passant en variables amplitude-phase $x$ et $\varphi$ (cf. Section 2.3.1) et en prenant en compte la variation dans le temps de la pression de souffle $y$, le modèle à un mode non moyenné (2.39) devient d'un système (2, 1)-rapides-lent de la forme

$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, y), \qquad (2.92a) $$

$$ \dot{y} = \epsilon g(y), \qquad (2.92b) $$

avec $\mathbf{x} \in \mathbb{R}^2$. Les sous-systèmes rapide et lent associés à (2.92) sont :

$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, y), \qquad (2.93a) $$

$$ \dot{y} = 0, \qquad (2.93b) $$

et

$$ \mathbf{0} = \mathbf{f}(\mathbf{x}, y), \qquad (2.94a) $$

$$ y' = g(y), \qquad (2.94b) $$

respectivement.

Dans cette situation où $\mathbf{x} \in \mathbb{R}^2$, le sous-système rapide (qui pour rappel représente le cas statique) possède deux solutions périodiques (une stable et une instable) en plus de la solution d'équilibre triviale. La variété critique de (2.92), solution de l'Éq. (2.94a), ne représente donc plus toutes les solutions permanentes de (2.93). Le domaine de bistabilité du diagramme de bifurcation de (2.93) incluant les solutions périodiques (qui n'est donc plus la variété critique) est en revanche encore D avec toujours une branche attractive correspondant à l'équilibre trivial (2.93). Les équilibres non triviaux, stable et instable, du système moyenné (2.76) correspondent à des solutions périodiques (2.93) de même stabilité. On peut donc en déduire qu'il existe des variétés invariantes (de dimension 2) associées aux solutions périodiques stable et instable de (2.93), l'une attractive et l'autre répulsive, respectivement. En effet, pour un système de dimension 2, une solution périodique stable (resp. instable) est associée à deux multiplicateurs de Floquet, l'un est égal à 1¹⁶ et l'autre de module inférieur (resp. supérieur) à un. Par conséquent, la variété invariante associée à la solution périodique instable de l'Éq. (2.93) est effectivement répulsive et la procédure de renversement temporal peut être utilisée.

L'Annexe H.6.2 illustre les éléments évoqués dans le paragraphe précédent. Dans ce cas, comme énoncé dans le paragraphe précédent, la méthode de renversement temporelle reste pertinente.

### 2.5.2.4 Calcul de séparatrice par apprentissage supervisé

Nous abordons aussi cette thématique relativement récente sur la notion de bassin d'attraction dynamique appliquée à des modèles d'instruments de musique à l'aide d'une méthode d'apprentissage supervisé appelée *machine à vecteurs de support* (ou *support vector machine*) [46, 92]. Toujours en utilisant le modèle à un mode (2.39) avec cette fois un profil de type rampe pour la pression dans la bouche, c.-à-d. une croissance linéaire stoppée à un temps $T$ après lequel la pression de souffle reste constante, égale $\gamma_t$ :

$$ \begin{aligned} \gamma(t) &= \gamma_0 + \epsilon \quad \text{if} \quad t < t_c \\ \gamma(t) &= \gamma_t \quad \text{if} \quad t \ge t_c. \end{aligned} \qquad (2.95) $$

<footnote n="16">16. Une solution périodique d'un sous-système rapide de dimension $m$ est associée à $m$ multiplicateurs de Floquet, mais l'un d'eux est toujours égal à 1 (voir le lemme 7.3 de [179]).</footnote>
<figure><img src="image_54.png" /><figcaption>**FIGURE 2.19.** Simulations numérique de (2.39) avec (2.95) pour deux profils de $\gamma(t)$ variant dans le temps (selon (2.95)) et un constant $\gamma(t) = \gamma_t$. La valeur initiale $\gamma_0$, la valeur cible $\gamma_t = 1.05$ sont les mêmes pour les deux simulations se faisant à $\gamma$ variant dans le temps. Seule la valeur de $\hat{\epsilon}$ change : $\hat{\epsilon} = 0.006$ (courbes jaunes) et 0.01 (courbes vertes). La simulation effectuée à $\gamma = \gamma_t$ est représentée en bleu. (a) Profils temporels de $\gamma$ et (b) de la pression acoustique simulée $p(t)$. Figure produite et fournie par Soizic Terrien$^6$.</figcaption></figure>

Noter ici que même si la pression de souffle varie dans le temps la notation $\gamma(t)$ est utilisée (et pas $y(t)$ comme précédemment). La valeur cible est choisie pour être dans la zone de bistabilité (zone grisée en Fig. 2.7).

L'effet d'un tel profil temporel de $\gamma$ sur les régimes observés est illustré en Fig. 2.19 qui montre montre la simulation numérique de (2.39) avec (2.95) pour deux profils de $\gamma(t)$ variant dans le temps (selon (2.95)) et un constant $\gamma(t) = \gamma_t$. Les trois profils de $\gamma(t)$ en fonction du temps sont représentés en Fig. 2.19(a) et les séries temporelles simulées de la pression acoustique $p(t)$ sont représentées en Fig. 2.19(b). Les conditions initiales $p(0)$ et $w(0)$ (avec $w(t) = \dot{p}(t)$) sont identiques pour chaque simulation. La valeur initiale $\gamma(0) = \gamma_0$, la valeur cible $\gamma_t = 1.05$ sont les mêmes pour les deux simulations se faisant à $\gamma$ variant dans le temps. La seule différence est la valeur de la pente $\hat{\epsilon}$ : $\hat{\epsilon} = 0.006$ et 0.01, ce qui correspond respectivement à des valeurs de $\epsilon$ de 0.22 et 0.37 (cf. Éq. (2.77)). Lorsque $\gamma$ et égal à $\gamma_t$, les conditions initiales choisies se trouve dans le bassin d'attraction de la solution d'équilibre, le régime permanent observé est donc non oscillant. Pour le $\hat{\epsilon}$ le plus grand (en vert en Fig. 2.19), on obtient la solution d'équilibre, comme dans le cas où $\gamma$ est constant (en bleu Fig. 2.19). Inversement, pour le plus petit $\hat{\epsilon}$ (en jaune en Fig. 2.19), la solution périodique est finalement observée. Ceci démontre que, pour des conditions initiales identiques sur $p, w$ et $\gamma$ et une même valeur cible $\gamma_t$, la valeur de $\hat{\epsilon}$ peut modifier la nature du régime permanent.

Les DBA sont déterminés directement à partir de simulations temporelles du système (en temps direct). Le principe est simplement de juger *a posteriori* si un jeu de conditions initiales ($p_0, \dot{p}_0, \gamma_0$) est situé dans le DBA de la solution d'équilibre ou dans celui de la solution périodique. Au lieu de calculer directement les DBA en « force brute » en explorant tout l'espace des conditions initiales, une classification par *machine à vecteurs de support* (ou *support vector machine*) [46, 92] est utilisée pour déterminer directement et à moindre coup de calcul (la frontière dans l'espace des phases du modèle) entre les simulations qui finissent pas osciller et celles qui s'éteignent, c.-à-d. ce que nous avons défini précédemment comme la SDBA.

La Fig. 2.20 montre un example de SDBA (surface bleue) obtenue par *machine à vecteurs de support*. Pour des conditions initiales choisies à l'intérieur de la surface bleue, le système converge vers la solution d'équilibre. Au contraire, pour des conditions initiales choisies à l'extérieur de la surface bleue, le système converge vers le régime périodique stable. La SDBA est comparée à l'extension suivant $\gamma(0)$ de la séparatrice entre SBA calculé pour $\gamma = \gamma_t$ (cette séparatrice est similaire à celle montrée en Fig. 2.8, c.-à-d. la solution périodique instable du cas statique à $\gamma = \gamma_t$). Si la dynamique de $\gamma$ n'avait aucun effet, les surfaces bleue et jaune seraient confondues.

Ces résultats mettent en évidence que la dynamique de la pression de souffle $\gamma$ affecte considérablement la nature du régime permanent observé pour une valeur cible donnée. En pratique, pour certaines valeurs de $\gamma_0$ un régime peut devenir très difficile à jouer. En effet, comme le montre la Fig. 6, il devient très difficile en pratique d'accéder au régime périodique quand $\gamma_0$ s'approche de 0, car presque toute condition initiale sur $p$ et $w$ conduit à la solution d'équilibre. En revanche, pour
<figure><img src="image_55.png" /><figcaption>**FIGURE 2.20.** SDBA obtenue par machine à vecteurs de support (surface bleue) comparé à l'extension suivant $\gamma(0)$ de la séparatrice entre SBA calculé pour $\gamma = \gamma_t$ (surface jeune). $\hat{\epsilon} = 0.006$ et $\gamma_t = 1.05$. Figure produite et fournie par Soizic Terrien$^6$.</figcaption></figure>

$\gamma_0 \in [0.1, 1]$ il est plus facile d'accéder au régime périodique dans le cas dynamique que dans le cas statique.

Cette approche, type « boite noire », permet d'aborder la question des bassins d'attraction dynamiques avec moins d'apriori sur le modèle que l'approche analytique détaillée dans ce manuscript. Un effort de comparaison entre les résultats obtenus avec ces deux approches complémentaires devra être réalisé à l'avenir.

Ce travail a fait l'objet d'un article [p2] (en fin d'écriture) et d'une communication dans un congrès international [ci12].

### 2.5.3 Bassin d'attraction dynamique dans de la cas d'un « rate-induced transitions »

La méthodologie présentée dans les sections précédentes est appliquée au cas d'une TIV ou *rate-induced transition* dans cette section. Pour cela, un exemple tiré de Ashwin et coll. [6] est traité. Nous n'allons pas analyser les résultats théoriques de cet article, mais uniquement d'appliquer la méthodologie présentée en Section 2.5.2 pour montrer qu'elle permet d'interpréter des résultats de simulations d'un modèle académique de la littérature subissant une TIV.

Le modèle étudié est un système (1, 1)-rapide-lent de la forme de (2.76). Ashwin et coll. utilisent $r$ (pour *rate*) à la place de $\epsilon$ utilisé dans ce manuscrit, inversent les noms des fonctions $f$ et $g$ et notent $\lambda$ la variable lente à la place du $y$ utilisé ici. Par conséquent les équations du modèle prennent la forme suivante :

$$ \dot{x} = g(x, \lambda) \qquad (2.96a) $$

$$ \dot{\lambda} = rf(\lambda). \qquad (2.96b) $$

avec

$$ g(x, \lambda) = - \left( (x + a + b\lambda)^2 + c \tanh(\lambda - d) \right) \left( x - K / \cosh(e\lambda) \right) \qquad (2.97) $$

où $a, b, c, d, e$ et $K$ sont des constantes, et

$$ f(\lambda) = -(\lambda - 1)(\lambda + 1). \qquad (2.98) $$
<figure><img src="image_56.png" /><figcaption>**FIGURE 2.21.** Fig. 4 de Ashwin et coll. [6] modifiée pour en faciliter la lecture. Variété critique attractive (lignes noires) et répulsive (ligne pointillée noire) et *pullback attractors*<fnref n="18" /> (ligne rouge avec une flèche noire pour montrer le sens du temps). (a) $r = r_a$, (b) $r = r_b$ et (c) $r = r_c$ avec $r_a < r_b < r_c$.</figcaption></figure>

La solution de (2.96b) est

$$ \lambda(t) = \frac{\lambda_0 \cosh(rt) + \sinh(rt)}{\lambda_0 \sinh(rt) + \cosh(rt)}. \qquad (2.99) $$

Le profil temporel de $\lambda$ est donc un profil exponentiel.

Nous notons $x = x_i^*(\lambda)$ ($i = 1, 2, 3$) les solutions de :

$$ g(x^*(\lambda), \lambda) = 0. \qquad (2.100) $$

Dans ce cas, la variété critique $M_0$ de (2.96) possède deux branches attractives $M_{0,a_1}$ et $M_{0,a_2}$ telles que

$$ M_{0,a_1} = \{x \in \mathbb{R}, \lambda \in \mathbb{R} \mid x = x_1^*(\lambda)\}, \qquad (2.101a) $$

$$ M_{0,a_2} = \{x \in \mathbb{R}, \lambda \in D =]\lambda_1, +\infty[ \mid x = x_2^*(\lambda)\}, \qquad (2.101b) $$

et une branche répulsive $M_0$ définie comme suit :

$$ M_{0,r} = \{x \in \mathbb{R}, \lambda \in D \mid x = x_3^*(\lambda)\}. \qquad (2.102) $$

L'intervalle $D$ représente donc le domaine de bistabilité de $M_0$.

Pour ce système académique, les branches $x = x_i^*(\lambda)$ ($i = 1, 2, 3$) ont les expressions analytiques suivantes :

$$ x_1^*(\lambda) = K \operatorname{sech}(e\lambda), \qquad (2.103a) $$

$$ x_2^*(\lambda) = -a - b\lambda + \sqrt{c \tanh(d - \lambda)}, \qquad (2.103b) $$

$$ x_3^*(\lambda) = a - b\lambda - \sqrt{c \tanh(d - \lambda)}. \qquad (2.103c) $$

Ces branches sont représentées en Fig. 2.21, tirée de Ashwin et coll. [6] puis modifiée pour en faciliter la lecture. Les paramètres utilisés sont :

$$ a = -0.25, \quad b = 1.2, \quad c = -0.4, \quad d = -0.3, \quad e = 3, \quad K = 2. \qquad (2.104) $$

La Fig. 2.21 montre aussi (en rouge avec une flèche noire pour montrer le sens du temps) des orbites particulières de l'Éq. (2.96) initialisées sur $M_{0,a_1}$, c.-à-d. sur un point $N_1$ tel que :

$$ N_1 = (\lambda_0, x_1^*(\lambda_0)) \qquad (2.105) $$

avec $\lambda_0 \to -1$<fnref n="17" />. Ces orbites particulières sont appelées *pullback attractors*<fnref n="18" /> par les auteurs. De gauche à droite, la Fig. 2.21 montre des résultats obtenus pour trois valeurs croissantes du paramètre

<footnote n="17">17. La valeur $\lambda_0 = -1$ rend $\lambda$ constant et donc égal -1 (cf. Éq. (2.99)).</footnote>

<footnote n="18">18. La notion de *pullback attractor* (voir par exemple [117]) n'est pas étudiée dans ce manuscrit. Ashwin et coll., dans le contexte de leur étude, définissent un *pullback attractor* par rapport à une branche stable donnée, ici $x_1^*(\lambda)$, comme une orbite $(\lambda(t), x(t))$ de l'Éq. (2.96) telle que $x(t) \to x_1^*(\lambda)$ quand $t \to -\infty$. Le *pullback attractor* est donc pour $M_{0,a_1}$ ce que notre SDBA est pour $M_{0,r}$.</footnote>
<figure><img src="image_57.png" /><figcaption>**FIGURE 2.22.** Reproduction (qualitative) de la Fig. 2.21 avec ajout de la SDBA et des lignes de courant de l'Éq. (2.96). Orbites de l'Éq. (2.96) (en rouge) et SDBA (en bleu), (c.-à-d. orbites de (2.96), dont on a inversé le temps, initialisées sur $M_{0,r}$ avec $\lambda_0 \to 1$). (a) $r = r_a$, (b) $r = r_b$ et (c) $r = r_c$ avec $r_a < r_b < r_c$.</figcaption></figure>

Pour la première valeur de $r$ (voir Fig. 2.21(a)) l'orbite commence donc sur $M_{0,a_1}$, s'en éloigne, puis la rejoint. Ensuite, pour la seconde valeur (voir Fig. 2.21(b)), l'orbite commence également sur $M_{0,a_1}$, s'en éloigne de nouveau, puis rejoint $M_{0,a_2}$. Enfin, quand on augmente encore $r$ (voir Fig. 2.21(c)), l'orbite commence toujours sur $M_{0,a_1}$, s'en éloigne davantage encore qu'en Fig. 2.21(a), puis finit par la rejoindre à nouveau.

Les résultats montrés en Fig. 2.21 peuvent s'interpréter en termes de bassin d'attraction dynamique comme montré en Fig. 2.22 qui reproduit des situations qualitativement similaire à celle montrée en Fig. 2.21¹⁹. Nous y voyons des orbites de l'Éq. (2.96) (en rouge) équivalentes à celles montrées en Fig. 2.21. En initialisant le système (2.96), dont on a inversé le temps, sur $M_{0,r}$, c.-à-d. sur un point $N_2$ tel que :

$$N_2 = (\lambda_0, x_3^*(\lambda_0)) \qquad (2.106)$$

avec $\lambda_0 \to 1$<fnref n="20" /> on obtient la SDBA (représentée ici en bleu) entre les DBA des variétés $M_{0,a_1}$ et $M_{0,a_2}$. Nous observons que, selon la valeur du paramètre $r$, une même condition initiale $N_1$ ne se situe pas dans le même DBA. Ainsi, en Fig. 2.22(a), pour $r = r_a$, $N_1$ se trouve d'abord dans le DBA de $M_{0,a_1}$ puis dans celui de $M_{0,a_2}$ pour $r = r_b$ en Fig. 2.22(b) et à nouveau dans le DBA de $M_{0,a_1}$ pour $r = r_c$ en Fig. 2.22(c).

Nous illustrons avec cet exemple que la SDBA peut varier fortement avec la vitesse de variation du paramètre, avec des implications importantes sur le comportement du système (à conditions initiales fixées).

<footnote n="19">19. En effet, Ashwin et coll. 2.21 ne fournissent pas les valeurs numériques des paramètres utilisés pour leurs simulations numériques.</footnote>

<footnote n="20">20. La valeur $\lambda_0 = 1$ rend $\lambda$ constant et égal 1 (cf. Éq. (2.99), en remplaçant $t$ par $-t$).</footnote>
## TRAVAIL DE RECHERCHE PÉRIPHÉRIQUE
### SUR LE CHAOS POLYNOMIAL

Ce travail de recherche s'est réalisé au cours de la thèse de Duc Thinh Kieu intitulée « Méthodes de réduction de modèles pour l'analyse du comportement dynamique des structures à paramètres incertains » [e6] dirigée par Sébastien Berger (co-directeur de thèse), Jean-Mathieu Mencik (co-directeur de thèse), Marie-Laure Gobert (co-encadrante de thèse) et moi-même (co-encadrant de thèse) et soutenue et soutenue le 10 février 2020.

Le manuscrit de thèse est disponible au lien suivant : https://www.theses.fr/2020ISAB0002

La conception mécanique dispose d'outils de calcul puissants capables de prédire le comportement dynamique de systèmes mécaniques. Cependant, afin d'assurer les niveaux vibratoires imposés par les cahiers des charges tout au long du cycle de vie de ces systèmes, il est impératif de prendre en compte les incertitudes dans leur cycle de conception. L'objectif est de concevoir des systèmes mécaniques optimaux et robustes, c.-à-d. des structures qui gardent leurs performances tout le long de leur cycle de vie malgré le caractère incertain de leurs paramètres. La première approche pour décrire le comportement dynamique de structures à paramètres incertains consiste à modéliser ces structures par éléments finis (EF) et à réaliser des simulations de ces modèles en effectuant un ensemble de tirages des valeurs des paramètres suivant leur loi de probabilité. Cette approche est connue sous le nom de méthode de Monte Carlo (MC). Elle s'avère néanmoins prohibitive en termes de temps de calcul dès lors que le nombre de ddl et de variables aléatoires est élevé. L'objectif de ce travail de thèse a été de développer des approches numériques s'appuyant d'une part sur des méthodes de prise en compte d'incertitudes basées sur des développements en chaos polynomiaux non intrusifs, et d'autre part, sur des modèles numériques « réduits » basées sous une sous-structuration de type Craig-Bampton, c'est-à-dire impliquant un nombre réduit de ddl. Dans une première partie, plusieurs méthodes de prise en compte d'incertitudes basées sur le chaos polynomial creux non intrusif sont développées afin de tenir compte d'un nombre significatif de paramètres incertains. Ces méthodes sont appliquées à l'étude du comportement dynamique d'un modèle d'embrayage et plus particulièrement à l'étude de la stabilité à partir de la méthode indirecte de Lyapunov. Une partie de ces travaux a fait l'objet d'une publication [a16] et d'une communication dans un congrès national [cn4]. Dans la deuxième partie (à laquelle ma participation a été moindre), la méthode basée sur le chaos polynomial creux dit anisotropique est associée à la méthode de sous-structuration dynamique de Craig-Bampton afin d'étudier les fonctions de réponse en fréquence d'une structure avec paramètres incertains. Cette approche est appliquée à une structure académique composée de trois plaques reliées par des interfaces de raideurs incertaines. Dans un premier temps, seule la méthode de réduction de Craig-Bampton est utilisée. L'impact de la troncature de la base modale est alors évalué. Dans un deuxième temps, le chaos polynomial creux anisotropique est utilisé seul puis associé à différents modèles réduits de type Craig-Bampton. Il est montré que les coûts en temps de calcul obtenus sont très largement inférieurs à ceux requis par la technique de Monte Carlo, pour des niveaux de précision et de confiance similaires.

En plus de ce travail de thèse, cette thématique a aussi été abordée dans le cadre de partenariats industriels, notamment par le biais de trois stages de Master 2 : [e1, e8] financés par l'entreprise Air Liquide dans le cadre des projets [Caline 1] et [Caline 2] (cf. Section 5.4.2) et [e3] financé par le CEA dans le cadre du projet [Sure-Methonum] (cf. Section 5.4.1).
## CONCLUSION GÉNÉRALE ET PERSPECTIVES

### Sommaire

<table><tr><td>4.1 Bilan des activités de recherche</td><td>97</td></tr><tr><td>4.2 Perspectives à court et moyen terme</td><td>97</td></tr><tr><td>4.2.1 Perspectives concernant le contrôle passif non linéaire de vibrations</td><td>97</td></tr><tr><td>4.2.2 Perspectives concernant les phénomènes transitoires dans les instruments de musique à anche</td><td>103</td></tr></table>

### 4.1 BILAN DES ACTIVITÉS DE RECHERCHE

C'est manuscrit a présenté mes travaux de recherche, obtenus depuis ma thèse au LAUM (2010-2013), en passant par le LMA (2013-2015) et aujourd'hui en 2023 au LaMé. J'ai eu la chance, dans le cadre de deux application différentes, de pouvoir développer une recherche avec un fil conducteur bien identifiable : celui de la dynamique non linéaire de systèmes rapides-lents. Je pense que dans ce manuscrit a pu transparaître mon goût pour la recherche fondamentale, avec un ancrage académique fort. Les perspectives énoncées dans la suite s'inscrivent encore dans cette logique même si le lecteur remarquera également mon souhait de donner plus de place à l'expérimentation dans mon travail, tout en l'orientant davantage vers le monde socio-économique.

### 4.2 PERSPECTIVES À COURT ET MOYEN TERME

#### 4.2.1 Perspectives concernant le contrôle passif non linéaire de vibrations

##### 4.2.1.1 Contrôle et récupération d'énergie vibratoire nuisible dans les éoliennes
*Coll. : S. Berger, M. Bavencoffe et G. Poulin-Vittrant*

Cette recherche future s'effectuera dans le cadre du projet CoREVE.

#### PROJET ACADÉMIQUE RÉGIONAL CoREVE (2024-2026)

Il s'agit d'un appel à « projets de recherche d'initiative académique » (APR IA) de la région Centre-Val de Loire remporté en juillet 2023 et qui débutera début 2024.

**Titre du projet** : Dispositifs passifs pour le contrôle et la récupération de l'énergie vibratoire nuisible dans les systèmes mécaniques tournants : application aux éoliennes

**Acronyme** : CoREVE (Contrôle et Récupération d'Énergie Vibratoire dans les Éoliennes)

**Descriptif du projet.** Dans le contexte de la transition écologique, les sources d'énergie renouvelables comme les éoliennes sont devenues incontournables. Cependant, ces dernières peuvent subir des vibrations structurelles nuisibles. En effet, l'industrie éolienne utilise des rotors de plus en plus grands pour capter davantage d'énergie et réduire les coûts. L'agrandissement de la structure rend les pales beaucoup plus flexibles et donc plus susceptibles de subir des vibrations de grande amplitude
causées notamment par des chargements aérodynamiques importants. Cela nuit aux performances de l'éolienne [2] et cause des dommages par fatigue [192]. Certaines vibrations, dites auto-entretenues, peuvent même endommager la structure [100, 157]. Un enjeu économique et industriel consiste donc à concevoir des dispositifs permettant d'éliminer ou du moins réduire ces vibrations.

Dans le passé des absorbeurs dynamiques linéaires passifs appelés TMD (Tuned Mass Damper) ont été utilisés pour tenter de réduire ces vibrations [44, 120, 141]. Les TMD sont efficaces lorsque les vibrations à réduire sont stationnaires et ne présentent qu'une composante fréquentielle principale. Ils deviennent cependant inefficaces dans la plupart des situations réelles où l'éolienne subit des variations temporelles de ses paramètres modaux et de son chargement aérodynamique (vibrations instationnaires et multi-fréquentielles).

Ces limitations conduisent les chercheurs et les industriels à développer des alternatives. Dans cette optique, les dispositifs actifs ont été largement étudiés. Fitzgerald et coll. [65] ont proposé un TMD actif. Une variante de ce dernier a été proposée par Fitzgerald et Basu [187]. Le contrôle des vibrations de rotors à vitesse variable a également été proposé dans [186]. Des systèmes de contrôle semi-actifs ont également été investigués par Arrigan et coll. [5]. Des TMD multiples [55] ou ayant plusieurs degrés de liberté [105] ont aussi été étudiés numériquement.

Les absorbeurs non linéaires bistable de type BNES sont choisis dans ce projet comme solution alternative aux TMD et aux NES monostables pour leur capacité à atténuer les vibrations de faibles amplitudes de façon passive et sur de larges gammes de fréquence. De plus, les réseaux d'absorbeurs sont également des dispositifs permettant de remédier à certaines limitations des TMD. Dans le cas de la réponse forcée harmonique d'une structure élastique, l'existence d'un réseau périodique de TMD peut notamment donner lieu à un phénomène de bandes de fréquences interdites [218] ou permettre une atténuation multimodale [167]. Il existe relativement peu de travaux dans la littérature sur les réseaux d'absorbeurs non linéaires. Un réseau de NLTVA peut permettre l'atténuation multimodales de structures primaires non linéaires [128]. Les réseaux de NES sont encore moins étudiés. L'existence du projet européen (Horizon 2020) MATASINK (https://doi.org/10.3030/896942) montre cependant qu'il s'agit d'une problématique académique et industrielle d'actualité. Il est intéressant de noter que dans [128, 167] l'utilisation de patchs piézoélectriques peut jouer le même rôle qu'un absorbeur purement mécanique.

Ainsi, l'objectif principal du projet CoREVE est d'étudier la capacité d'absorbeurs non linéaires de type BNES et/ou d'un réseau d'absorbeurs non linéaires à atténuer significativement les vibrations de systèmes mécaniques tournants avec une application aux éoliennes. La réduction des vibrations nuisibles dans les éoliennes permettra in fine d'en améliorer les performances. Il s'agira ensuite d'associer ce dispositif à des transducteurs piézoélectriques afin de récupérer l'énergie vibratoire (pour par exemple alimenter en énergie des capteurs positionnés sur les pales et utiles pour le contrôle de la santé de ces structures) et d'améliorer si possible l'efficacité du dispositif d'atténuation vibratoire.

> Par le biais d'une collaboration avec Maxime Bavencoffe (Maître de conférences, INSA Centre Val de Loire) et Guylaine Poulin-Vittrant (Directrice de recherche, CNRS), pour les parties du projet CoREVE traitant de la récupération d'énergie vibratoire, je m'appuierai sur l'expertise du GREMAN (UMR 7347) dans le domaine des matériaux piézoélectriques et de leurs applications

**Résultats préliminaires.** Dans ce contexte, trois stages Master 2 on permis d'obtenir des résultats préliminaires. Dans les deux premiers, nous avons étudions les vibrations causées entre autre par le couplage de modes entre le rotor et le mât d'une éolienne. Durant le stage de Riccardo Gega [e2], un modèle analytique linéaire d'éolienne reproduisant ce phénomène de couplage de modes et couplé à un NES a été proposé. L'étude a permis de montrer que le NES permet d'atténuer les vibrations causées par le couplage de modes. En effet, le modèle sans NES étant linéaire, l'amplitude des oscillations ne sature pas vers un cycle limite mais diverge vers l'infini et on observe que l'ajout du NES permet l'existence de régimes d'oscillations d'amplitude finie, dits atténués. Cependant, bien que finies, les
amplitudes observées restent trop importantes. Pendant le stage de Joris Nguyen [e7], un autre modèle analytique d'éolienne incluant dans pales non linéaires a été proposé permettant d'obtenir, lorsqu'une instabilité est observée, non plus des régimes divergents mais des solutions périodiques ou quasi-périodiques. Ceci permet d'évaluer l'amplitude des régimes oscillants que le NES devra atténuer. Ce travail de modélisation focaliser sur le phénomène de couplage de mode a permis de nous familiariser avec les modèles d'éoliennes. Il s'avère que bien qu'extrêmement dangereux, ce phénomène d'instabilité vibratoire se produit en général à des vitesse de rotation du rotor plus élevées que les vitesses nominales d'une éolienne. Par conséquent, dans le projet CoREVE un plus large éventail de problèmes vibratoires dans les éoliennes sera considéré. Nous nous orienterons vers des modèles plus réalistes (obtenu par la méthodes des éléments finis par exemple). Concernant le comportement du BNES, nous étudierons d'abords des modèles de petites tailles dont l'analyse nous permettra l'acquérir de l'intuition sur leur comportement très riche et complexe. C'est dans cette optique qu'une première étude a été réalisée pendant le stage d'Helmi Chaabene [e4]. Cette étude numérique a traité d'un système académique comportant un système primaire linéaire (masse-ressort-amortisseur) forcé harmoniquement et connecté à un BNES. À l'aide d'intégration temporelle numérique de ce système, six régimes de fonctionnement ont été recensés. Ces régimes sont de même nature que ceux observé dans [a13] lorsque que le système primaire est un oscillateur de Van der Pol (cf. Annexe C). Ce travail a également permis de localiser les régions de forte efficacité dans l'espace des paramètres du BNES.

**Principales étapes du projet.** Les principales étapes du projet CoREVE seront :

* **Phase (1).** Compléter l'étude numérique/analytique préliminaire [e4]. Dans cette étape il s'agira d'associer des études paramétriques obtenues par l'intégration temporelle numérique de modèle académiques de petite taille à des études analytiques basées sur la méthodologie générale proposée dans [a13]. Ces modèles académiques auront pour but de représenter un phénomène vibratoire donné rencontré dans les éoliennes : forçage (harmonique, large bande et/ou aléatoire et/ou pouvant être d'amplitude relativement faible) et instabilités dynamiques. Le cas d'un d'un système primaire linéaire forcé harmoniquement et connecté à un BNES sera traité en premier, la démarche sera ensuite étendue au cas :

(a) des autres phénomènes vibratoires énoncés plus haut ;

(b) d'un modèle académique mécano-électrique incluant un patch piézoélectrique pour la récupération d'énergie vibratoire. Il s'agira de trouver un compromis entre atténuation de vibration et récupération d'énergie ;

(c) au cas d'un réseau de BNES seuls puis associés à des patchs piézoélectriques.

* **Phase (2).** À partir de l'expertise acquise dans la phase précédente cette Phase (2) consistera à établir un modèle réaliste mécanique d'éolienne pouvant reproduire plusieurs des phénomènes vibratoires étudiés dans le Phase (1). Il s'agira ensuite de développer un modèle mécano-électrique d'une éolienne couplée à un dispositif d'atténuation non linéaire passif de vibrations (à base d'un ou plusieurs BNES) incluant un système récupération d'énergie vibratoire à base de patchs piézoélectriques. Cette phase de modélisation se déroulera en trois étapes :

(a) Étape I (éolienne seule) : (i) Établissement, à partir de l'étude bibliographique, d'un modèle d'éolienne incluant le couplage mât-pale et les principales sources de chargement aérodynamique. (ii) Simulation numériques temporelle du modèle pour la caractérisation de son comportement vibratoire. Il s'agira ici d'être en mesure de reproduire des résultats choisis de la littérature.

(b) Étape II (éolienne + dispositif BNES) : (i) Couplage du modèle précédemment établi à un dispositif constitué d'un ou plusieurs BNES. (ii) Simulation numériques temporelle du modèle couplé. Il s'agira aussi de comparer les résultats obtenus à ceux de la littérature afin de mettre en évidence l'apport des méthodes proposées dans CoREVE.
<figure><img src="image_58.png" /><figcaption>**FIGURE 4.1.** Photo du banc d'essai de l'INSA Centre Val de Loire. Le système primaire est une poutre sollicitée harmoniquement en flexion à l'aide d'un pot vibrant. Le BNES y est constitué d'une poutre en flambage sur laquelle est attachée une masse supposée ponctuelle en son milieu.</figcaption></figure>

(c) Étape III (éolienne + dispositif BNES + dispositif de récupération d'énergie) : (i) L'exper-
tise du GREMAN sur les matériaux piézoélectriques et leur modélisation va être sollicitée
afin de proposer un modèle mécano-électrique réaliste et d'en caractériser son rendement
électrique. (ii) Dans cette dernière étape de l'établissement du modèle complet, un bilan
sera effectué afin d'envisager une application des résultats obtenus dans le cadre d'une
collaboration future avec un partenaire industriel.

* **Phase (3).** Effectuer une campagne de mesure afin de corréler des résultats expérimentaux avec les résultats théoriques/numériques obtenus dans les phases précédentes. La validation expérimentale ne se fera dans un premier temps que sur le modèle poutre-BNES en forçage harmonique. Pour cela, le dispositif expérimental utilisé s'appuiera sur celui développé en 2023 lors du projet DPT2MA (cf. Section 5.4 et Fig. 4.1). Il s'agit d'un banc d'essai, inspiré de la littérature [103], incluant une poutre sollicitée harmoniquement en flexion à l'aide d'un pot vibrant. Le BNES y est constitué d'une poutre en flambage sur laquelle est attachée une masse supposée ponctuelle en son milieu. Les mesures sont effectuée à l'aide de vibromètres laser (non présents sur la Fig. 4.1). Il faudra dans un premier temps retrouver expérimentalement le comportement typique d'un BNES pour ensuite y associer un patch piézoélectrique.

4.2.1.2 Analyse lente-rapide approfondie et collaboration avec la communauté des ma-
thématiques appliquées

Coll. : N. Berglund

Analyse des solutions de la dynamique rapide d'un oscillateur auto-entretenu couplé à un BNES. L'étude effectuée dans [a13] (cf. Section 1.2.4 et Annexe C montre que les solutions de la dynamique rapide du flot de modulation ne sont pas toujours des solutions d'équilibre. La fig. 1.20 montre en effet que l'on peut observer des régimes régimes périodiques, quasi-périodiques et chaotique. La connaissance de la stabilité de la variété critique du flot de modulation (dont pour rappel chaque point correspond aux solutions d'équilibre) n'est plus suffisante pour effectuer l'analyse de stabilité globale complète des systèmes étudiés. Par exemple, en Fig. 1.19 nous observons que même si la structure de la variété critique nous permet de comprendre l'évolution de l'orbite du flot de
<figure><img src="image_59.png" /><figcaption>**FIGURE 4.2.** Similaire à la Fig. 1.19 sauf que l'intégration numérique du flot de modulation (1.68) est effectuée pour une valeur de σ dépassant la limite d'atténuation $\sigma_{ml}$.</figcaption></figure>

modulation (en bleu), elle ne fournit pas de critère permettant de prédire l'apparition de chaque phase rapide où l'orbite saute d'une branche de variété attractive à une autre. La Fig. 1.19 montre que c'est la phase lente 3 (FE3) qui détermine si la situation va rester inoffensive où si elle va devenir dangereuse. En effet lorsque l'on augmente le paramètre de bifurcation $\sigma$ (cf. Fig. 4.2), on observe que l'orbite du flot de modulation saute sur la branche attractive droite de $M_{01}$ (correspondant à la ligne rouge la plus à droite et notée $M_{01}^{a,2}$ en Annexe C) et arrive au dessus du point fixe instable (point noir), l'orbite longe $M_{01}$ vers le haut ce qui produit une situation dangereuse. Plus rarement, il arrive aussi que le saut rapide vers $M_{01}^{a,2}$ ne soit pas précédé d'une phase lente chaotique mais d'une phase lente périodique ou quasi-périodique

Dans ce contexte je projette les perspectives de recherche suivantes :

**À court terme.** Une étude systématique des solutions du sous-système rapide et des bassins d'attraction des solutions stables. En traçant ces solutions (avec leur stabilité et les bassins d'attraction des solutions stable) dans l'espace des phases pour un jeu de paramètres du système primaire et du BNES donné on obtiendrait l'équivalent de ce qu'est la variété critique pour les solutions d'équilibre. Même s'il parait peu probable d'obtenir ce type d'informations analytiquement (à part éventuellement s'il s'agit de solutions périodiques), même obtenues numériquement elles permettront d'avancer encore davantage sur la compréhension du système. Notamment en mettant en évidence les conditions d'apparition des phases rapides et leurs liens avec la dynamique lente. Il s'agirait par exemple dans le cas de la Fig. 1.19 de connaître la valeur de la variable lente $a$ (qui est un paramètre pour le sous-système rapide) à partir de laquelle les oscillations chaotiques de la dynamique rapide emmènent l'orbite du flot de modulation dans le bassin d'attraction de la branche attractive $M_{01}^{a,2}$. Puis d'observer la position relative de cette valeur par rapport au point fixe instable situé sur $M_{01}^{a,2}$.

**À moyen (et long) terme.** L'étude du sous-système rapide mentionnée ci-dessus permet d'obtenir beaucoup d'informations sur le comportement du système. Cependant, le comportement asymptotique quand $\epsilon \to 0$ ainsi décrit n'est pas toujours totalement représentatif du comportement du flot de modulation. Nous avons par exemple vu en Section 1.2.5 (cf. Fig. 1.21) qu'au passage des points-col de la variété critique une loi d'échelle était nécessaire pour décrire correctement la trajectoire du flot de modulation. Dans ce cas la loi d'échelle permet en effet d'approcher localement la variété invariante $M_{\epsilon}$ du flot de modulation.

Cette étude asymptotique ne prévoit pas non plus le phénomène de retard à la bifurcation
<figure><img src="image_60.png" /><figcaption>**FIGURE 4.3.** Zoom de Fig. 1.19 effectué au niveau du point de bifurcation de Hopf ($a^{\text{FP}}$, $b^{\text{FP}}$, $c^{\text{FP}}$) du sous-système rapide (1.69) (cf. (C.51)).</figcaption></figure>

évoquée en Section 2.4. La Fig. 4.3 montre que ce phénomène se produit également pendant le régime CR1 au niveau du point de bifurcation de Hopf ($a^{\text{FP}}$, $b^{\text{FP}}$, $c^{\text{FP}}$) du sous-système rapide (1.69) (cf. (C.51)).

Dans ces situations il est donc nécessaire de connaître les variétés invariantes associées aux solutions du sous-système rapide. Cela lorsqu'il s'agit de solutions d'équilibre (comme en Section 1.2.5) mais également de solutions périodiques, quasi-périodiques et chaotiques. Dans le cas de solutions périodiques les résultats généraux proposés par Berglund et coll. [21, 22], avec qui j'entame une collaboration (cf. paragraphe suivant), pourraient être mis au œuvre.

**Effet d'un forçage stochastique faible sur le phénomène de pompage énergétique.** Le travail de recherche décrit dans ce paragraphe s'effectuera en partie dans le cadre de la thèse d'Israa Zogheib que je co-encadrerai avec Nils Berglund (directeur de thèse) de l'Institut Denis Poisson (UMR CNRS 7013, Orléans) et spécialiste des systèmes dynamiques rapides-lents stochastiques [20, 22]. Le fil conducteur de ces travaux de thèse sera d'appliquer des résultats théoriques connus sur les équations différentielles ordinaires lentes-rapides perturbées par un bruit à des systèmes mécaniques concrets. Ces applications permettront également de mettre en évidence des extensions nécessaires de la théorie afin de décrire correctement les phénomènes physiques mis en jeu, conduisant à une meilleure compréhension de ces derniers et à des résultats mathématiques originaux.

L'objectif initial sera d'obtenir des résultats théoriques permettant d'expliquer les résultats numériques obtenus dans [a3] (cf. Section E). Dans ce cadre, l'idée de départ (cf. Fig. 1.24) sera d'estimer la probabilité *Pr* qu'après chaque passage au niveau du point-col gauche de la variété critique une réalisation du processus stochastique (1.94) (en rouge en Fig. 1.24) arrive sur la partie attractive droite de $M_0$ au dessus du plus grand point fixe du flot de modulation (point vert en Fig. 1.24).

L'étude se focalisera d'abord sur la forme normale d'une bifurcation col-nœud dynamique (cf. Éq. (1.77)) avec un forçage stochastique faible sur la variable lente. En effet, en appliquant la méthode de réduction à la variété centrale décrite en Section 1.2.5.2, nous pouvons montrer qu'au voisinage du point-col de la variété critique le flot de modulation stochastique (1.94) peut s'écrire sous la forme
suivante :

$$ \hat{\epsilon}q' = h(q, v) = -q^2 - v \qquad (4.1a) $$

$$ v' = 1 + \sqrt{\hat{\epsilon}\hat{\sigma}\xi(\hat{\tau})}. \qquad (4.1b) $$

La différence de signe entre les membres de droite de (1.77a) et (4.1a) s'explique des hypothèses de départ différentes des méthodes de moyennisation utilisées, c.-à-d. (1.5) et (E.6)-(E.7), respectivement. Le bruit n'agissant que sur la variable lente (1.94) et la réduction à la variété centrale ne portant que sur les variables rapides, l'obtention de (4.1) ne pose pas de difficultés supplémentaires par rapport au cas déterministe.

Dans le plan de phase $(q, v)$, si le processus (4.1) commence à $(q_0, v_0)$ près d'un point de la partie attractive $M_{0,a}$ de la variété critique, lorsque la trajectoire atteint une valeur donnée $Q$ de la variable rapide $q$, nous détectons la valeur obtenue pour la variable lente $v$, notée $V$. Nous définissons ensuite la densité de probabilité (PDF pour probability density function) $R_Q(v)$ comme suit : pour une valeur donnée $Q$ de $q$, la probabilité que $V$ se situe entre $v$ et $v + dv$ est $R_Q(v)dv$. La fonction de répartition (CDF pour cumulative distribution function) correspondante est notée $P_Q(v)$. Ceci est illustré en Fig. 4.4 pour $\hat{\epsilon} = 0.012$, $\hat{\sigma} = 0.75$ et $Q = -5$. On notera que $R_Q(v)$ est ici représenté schématiquement sur la Fig. 4.4 par une fonction gaussienne (ce qui n'est pas forcément le cas en réalité).

<figure><img src="image_61.png" /><figcaption>**FIGURE 4.4.** Position du problème. Dans le plan de phase $(q, v)$ : la variété critique $M_{0,c.}$-à-d. $q = \pm\sqrt{-v}$ (ligne continue verte) ; la solution déterministe similaire (D.52) au signe près (ligne pointillée grise) ; 11 réalisations du processus stochastique (4.1) (lignes continues avec différentes teintes de rouge). La PDF $R_{-5}(v)$ est également représentée de manière schématique (ligne bleue continue).</figcaption></figure>

En revenant au flot de modulations initiale (1.94), nous pensons que la connaissance de $R_Q(v)$ et $P_Q(v)$ permettra de caractériser (du moins en partie) la probabilité $Pr$ évoquée plus haut.

### 4.2.2 Perspectives concernant les phénomènes transitoires dans les instruments de musique à anche

#### 4.2.2.1 Bassins d'attraction dynamiques de modèles d'instruments complexes

*Coll. : S. Terrien et C. Vergez*

Par rapport au modèle simpliste à un mode utilisé dans les sections précédentes, la complexification par étape du modèle passerait par :
(i) **Le raffinement du modèle physique de l'instrument (c.-à-d. de l'Éq. (2.3a))**. Une première étape serait d'utiliser un modèle à deux modes tel que proposé par Doc et coll. [56]. Le diagramme de bifurcation de ce dernier pouvant posséder des domaines de bistabilité entre régimes périodiques, ce modèle serait un outil pertinent, c.-à-d. simple mais suffisant, pour étudier les séparatrices entre les DBA de deux solutions périodiques. Dans ce cas, une analyse de stabilité locale et globale statique sera nécessaire au préalable afin de déterminer le diagramme de bifurcation statique complet et les bassins d'attraction statiques du modèle.

Nous observons en Annexe H.6.2 (voir aussi la Section 2.5.2.3) la limite de la méthode de renversement temporel pour le calcul de le SDBA associée à des variétés invariantes *saddle-type*. Dans ce cas, des méthodes de continuation d'orbites, déjà initiées dans [cn6] à l'aide du logiciel AUTO [57], pourraient être utilisées. L'intérêt des méthodes de continuation, par rapport à l'intégration temporelle, est qu'elles peuvent suivre des solutions instables et sont donc pertinentes pour le calcul de variétés *saddle-type* le cas échéant. Dans ce contexte, des travaux issus de la communautés des mathématiques appliquées [61, 89], proposent des algorithmes de calcul de telles variétés avec des applications dans le domaine des neurosciences.

Par ailleurs, nous avons observé en Fig. 2.17(b) que la séparatrice entre deux DBA pouvait impliquer de très petites amplitudes. Il est difficile d'envisager cette situation dans un problème réel (ou même dans un calcul numérique réalisé avec une précision finie) où un bruit résiduel empêcherait la séparatrice d'atteindre des amplitudes aussi faibles (comme dans le cas de la retard à la bifurcation décrit en Section 2.4). D'un point vue numérique la problématique du bruit peut être abordée naïvement par une approche de type Monte-Carlo [64] qui consiste obtenir des données statistiques (moyenne, écart-type, ...) et/ou à reconstruire des fonctions de densité de probabilité de grandeurs d'intérêt à partir d'un ensemble de réalisations d'un processus stochastique. Avec le même modèle que celui étudié en Section 2.5.2 mais perturbé par un bruit blanc (tel que présenté en Section 2.4), le plus simple serait par exemple une approche en temps direct<fnref n="1" /> de type « force brute » qui consisterait à calculer pour chaque point de l'espace des conditions initiales ($p_0, \dot{p}_0, \gamma_0$) la probabilité d'observer des oscillations quand le domaine de bistabilité est traversé. L'approche de Monte-Carlo est cependant très couteuse en temps de calcul car il faut un nombre important de réalisations pour obtenir des résultats précis en termes de grandeurs statistique et encore plus lorsque l'on souhaite reconstruire des fonctions de densité de probabilité. De plus c'est une approche de type « boite noire » et donc qui ne permet pas l'explication des mécanisme sous-jacents aux comportements observés. Par conséquent, de façon complémentaire, un approche analytique doit aussi être envisagée (voir Section 4.2.2.2).

(ii) **La prise en compte de la dynamique temporelle d'autres paramètres de contrôle.**
Dans un premier temps, il s'agirait de prendre en compte la dynamique temporelle du paramètre $\zeta$ en remplaçant $g(y)$ dans l'Éq. (2.3a) par $g(y)$ avec $y \in \mathbb{R}^2$. Les profils temporels utilisés seraient déterminés par l'intuition ou préférablement à partir de mesures *in vivo* réalisées à l'aide de becs instrumentés.

### 4.2.2.2 Bassins d'attraction dynamiques et analyse lente-rapide approfondie : études déterministe et stochastique

*Coll. : N. Berglund*

Le travail de recherche décrit dans cette section s'amorcera également dans le cadre de la thèse d'Israa Zogheib (c.f. § *Effet d'un forçage stochastique faible sur le phénomène de pompage énergétique* en Section 4.2.1.2). L'objectif initial de la thèse sera dans ce contexte d'écrire le modèle d'instrument

<footnote n="1">1. En effet, dans le cas stochastique le renversement temporel n'est pas trivial comme dans le cas déterministe (c.-à-d. uniquement remplacer $t$ pat $-t$) et même s'il existe des travaux sur les équations différentielles rétrogrades [158], c'est un domaine de recherche à part entière en mathématique.</footnote>
dont la pression de souffle varie dans le temps, c.-à-d. l'Éq. (2.76)), réécrite ici comme suit :

$$ \frac{dx}{dy} = \frac{1}{\epsilon} \frac{f(x, y)}{g(y)} = \frac{1}{\epsilon} F(x, y) \qquad (4.2) $$

sous la forme d'un système dynamique ayant la forme générale suivante (similaire à l'Éq. (G.32)) :

$$ \frac{du}{ds} = \frac{1}{\epsilon} \hat{F}(u, s) \qquad (4.3) $$

où s pourrait être vu comme une abscisse curviligne paramétrant la SDBA et $u(s)$ la projection orthogonale en $s$ du point $(y, x)$ d'une orbite quelconque de (4.2) (voir Fig. 4.5). Sans perte de généralité, l'orientation est choisie telle que $u(s) > 0$ si $(y, x) \in \text{DBA}_2$. La difficulté réside dans le fait que nous ne possédons pas de solution analytique pour la SBDA (et nous n'en aurons probablement pas) et dans la définition même de $u$. Dans ce contexte la méthodologie proposée par Berglund et Gentz [21, 25] sur la paramétrisation d'orbites périodiques (qui peut être vue comme une généralisation de la méthode de moyennisation) semble pouvoir être applicable.

<figure><img src="image_62.png" /><figcaption>**FIGURE 4.5.** Schéma du changement de variables $(y, x) \to (s, u)$. $s$ est l'abscisse curviligne paramétrant la SDBA (en rouge) et $u(s)$ est la projection orthogonale en $s$ du point $(y, x)$ d'une orbite quelconque de (4.2) (en noir). Sans perte de généralité, l'orientation est choisie telle que $u(s) > 0$ si $(y, x) \in \text{DBA}_2$.</figcaption></figure>

Le problème mis sous la forme (4.3) devient trivial d'un point de vue déterministe. En effet, si $u(0) < 0$ (resp. $u(0) > 0$) alors $u(0) \in \text{DBA}_1$ (resp. $u(0) \in \text{DBA}_2$) et la branche $M_{\epsilon,a_1}$ (resp. $M_{\epsilon,a_2}$) sera nécessairement atteinte. D'un point de vue stochastique par contre, même si la condition initiale est dans l'un des DBA, le bruit peut très bien faire basculer la trajectoire dans l'autre et faire que le comportement attendu d'un point vue déterministe n'est pas observé. En s'inspirant des résultats décrits en Section 2.4, le problème stochastique peut s'écrire sous la forme suivante :

$$ dx_y = \frac{1}{\epsilon} F(x_y, y) dy + \frac{\sigma}{\sqrt{\epsilon}} dW_y. \qquad (4.4) $$

en fonction des variables $y$ et $x$, et

$$ du_s = \frac{1}{\epsilon} \hat{F}(u_s, s) ds + \frac{\sigma}{\sqrt{\epsilon}} dW_s. \qquad (4.5) $$

en fonction des variables $s$ et $u$.

Dans un premier temps, dans une approche locale, nous souhaitons étudier le comportement du système dans le voisinage de la SDBA. Le problème devient donc similaire à celui traité en Section 2.4 dans le cadre de l'étude du retard à la bifurcation. À la différence qu'ici nous ne nous intéressons pas
à déterminer un point de bifurcation dynamique mais à savoir si une trajectoire commençant à une
certaine position initiale ($s_0, u_0$) (avec par exemple $u_0 < 0$) atteindra la branche attractive $M_{\epsilon,a_2}$
(comme attendu dans le cas déterministe) ou l'autre branche attractive $M_{\epsilon,a_1}$. L'intérêt de la forme
(4.5) est que dans le voisinage de la SDBA nous pouvons justifier la linéarisation $\hat{F}(u_s, s)$ autour de
$u_s = 0$ (voisinage de la SDBA). Ceci n'est pas possible avec (4.4) car le voisinage de la SDBA ne
correspond pas forcément à $x$ ou $y$ petit. En notant $\hat{a}(s) = \partial_u \hat{F}(0, s)$, la linéarisation de (4.5) s'écrit :

$$
du_s = \frac{1}{\epsilon} \hat{a}(s) u_s ds + \frac{\sigma}{\sqrt{\epsilon}} dW_s. \qquad (4.6)
$$

Les résultats existants sur un système de la forme de (4.6) dans le cas d’une bifurcation fourche dynamique en présence de bruit ([22], Section 3.4) pourront dans un premier temps être adaptés à la structure plus complexe du modèle de clarinette. La question est de savoir si cela pourra être généralisé à des des systèmes de dimension supérieure, par exemple si l’on souhaite étudier un cas de bistabilité entre deux solutions périodiques.

4.2.2.3 Sourdine pour instruments à anche basée sur le pompage énergétique

Coll. : C. Vergez et P.-O. Mattei

Les applications déjà existantes du pompage énergétique en vibroacoustique laissent en effet penser
que les NES pourraient être utilisés pour concevoir des sourdines efficaces pour instruments à anche et
à cordes frottées. En effet, les solutions existantes ne sont à l'heure actuelle pas satisfaisantes alors que
les produits efficaces pour les autres instruments (piano, cuivres, ...) suscitent un grand intérêt parmi
les musiciens amateurs et professionnels. On peut illustrer ceci par le succès de la gamme « *Silent Instruments* » de YAMAHA. D'autre part, la nécessité qu'ont certains instrumentistes citadins de
s'exercer quotidiennement les rendent prêts à acquérir des dispositifs très couteux, comme les cabines
acoustiques par exemple.

Cette thématique fait le lien entre mes deux activités de recherche principales. Elle envisage une
application des NES en acoustique musicale qui me tient particulièrement à coeur. Initiée il y a un
moment lors du Stage de Master 2 de Barthélémy Jeandel [e5] qui avait montré qu’un NES purement
acoustique (tel que proposé par Gourdon et coll. [85]) était difficilement envisageable, cette recherche
a été continuée de façon épisodique en considérant une membrane pour la conception du NES. Des
verrous scientifiques ont été mis en évidence, notamment liés aux dimensions des membranes qui
devraient être utilisées mais également à la nature des modèles d’instrument à anche.

En effet, comme nous l'avons évoqué plus haut, dans les modèles d'instruments à anche l'appa-
rition d'un son musical correspond à la perte de stabilité de la solution d'équilibre au profit d'une
solution oscillante (en général périodique). Dans les modèles plus complexes que le modèle à un mode,
comme le modèle à deux modes mentionné en Section 4.2.2.1, cette perte de stabilité s'effectue par
les bifurcations de Hopf successives des modes du modèle. C'est donc un cas d'instabilités multiples
(c.-à-d. que plus d'une paire de valeurs propres complexes conjuguées de la matrice jacobienne du
modèle évalué en la solution d'équilibre ont des parties réelles positives). Les travaux présentés en
Section 1.2.3 qui synthétisent les résultats de [a9] peuvent être vus comme une preuve de concept
quant à l'utilisation d'un NES pour atténuer des régimes oscillants résultant d'instabilité double (c.-
à-d. deux modes instables). Qu'en est-il lorsqu'il y a plus que deux modes instables ? Comment ces
résultats théoriques peuvent-ils aider pour la conception d'un NES ? Cette conception mènera-t-elle
à un NES réalisable pratiquement ? Les regimes du système couplé instrument-NES sont-ils perti-
nents musicalement ? Des questions auxquelles il faudra répondre si l'on veut espérer une application
industrielle.

Afin d'aborder cette thématique également d'un point de vue expérimental, une collaboration
avec le LMA semble opportune étant donnée l'expertise de ce dernier dans la conception de NES à
base membrane de petite taille.
4.2.2.4 Recherche de profil optimaux et prise en compte de l'interaction instrumentiste-instrument

La notion d'optimisation me paraît importante dans le domaine de la pratique musicale. En effet, on peut supposer que l'instrumentiste, quel que soit son niveau, essaye de jouer au mieux selon ses capacités et son instrument. Il s'agit ici d'optimisation objective, puisque sont recherchés des profils temporels des paramètres de contrôle (ou la forme de la fonction g(y) dans l'Éq. (2.3)) permettant d'accomplir une tâche donnée de façon optimale (une attaque de note la plus rapide possible par exemple) et ce, sans aucune contrainte liée au niveau d'apprentissage de l'instrumentiste ou à la qualité de l'instrument.

à plus long terme, il s'agirait d'être en mesure de modéliser des situations de contrôle avec contre-réaction, c.-à-d. quand le modèle d'instrument (2.3) prend la forme suivante :

$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, y), \qquad (4.7a) $$

$$ \dot{y} = \epsilon g(\mathbf{x}, y), \qquad (4.7b) $$

et devient un modèle instrument-instrumentiste.

La contre-réaction du musicien serait définie comme les corrections apportées par ce dernier « en temps réel » aux paramètres de contrôle pour arriver au résultat souhaité de manière optimale pour lui. En effet, comme énoncé précédemment, on part du principe que l'instrumentiste essaye de jouer au mieux selon ses capacités à un instant donné de son apprentissage mais aussi de son instrument. Dans cette partie on supposera que les caractéristiques de l'instrument sont connues et on se concentrera sur le contrôle du musicien.

La question de la modélisation de la contre-réaction du musicien sur l'instrument, c.-à-d. la forme de la fonction g(x, y) sera abordée dans le cadre général de la théorie du contrôle et de l'automatique (En particulier leurs applications à la description des mouvements humains [29, 104, 139]) et/ou en s'inspirant de mesure in vivo. Une première approche consistera, comme précédemment, à résoudre un problème d'optimisation sans contrainte. Cette démarche ne prend cependant en compte ni les limites physiologiques de l'être humain - et donc peut mener à des résultats irréalisables même pour un musicien expert - ni le niveau d'apprentissage de l'instrumentiste et par conséquent ne permet pas de faire la différence entre un instrumentiste débutant et un musicien expert. Une seconde étape serait encore de résoudre un problème d'optimisation mais en y ajoutant des contraintes liées à la physiologie et au niveau du musicien. Tout le problème réside alors dans la définition de ces contraintes.

Cette recherche nécessiterait un rapprochement avec la communauté des sciences du mouvement.
# SELECTION D'ARTICLES
<figure><img src="image_63.png" /></figure>
*Mechanical Systems and Signal Processing 131 (2019) : 728-50 [a8]*

# STEADY-STATE REGIMES PREDICTION OF A MULTI-DEGREE-OF-FREEDOM UNSTABLE DYNAMICAL SYSTEM COUPLED TO A SET OF NONLINEAR ENERGY SINKS

Baptiste BERGEOT and Sergio BELLIZZI

## ABSTRACT

A general method to predict the steady-state regimes of a multi-degree-of-freedom unstable vibrating system (the primary system) coupled to several nonlinear energy sinks (NESs) is proposed. The method has three main steps. The first step consists in the diagonalization of the primary underline linear system using the so-called biorthogonal transformation. Within the assumption of a primary system with only one unstable mode the dynamics of the diagonalized system is reduced ignoring the stable modes and keeping only the unstable mode. The complexification method is applied in the second step with the aim of obtaining the slow-flow of the reduced system. Then, the third step is an asymptotic analysis of the slow-flow based geometric singular perturbation theory. The analysis shows that the critical manifold of the system can be reduced to a one dimensional parametric curve evolving in a multidimensional space. The shape and the stability properties of the critical manifold and the stability properties of the fixed points of the slow-flow provide an analytical tool to predict the nature of the possible steady-state regimes of the system. Finally, two examples are considered to evaluate the effectiveness and advancement of the proposed method. The method is first applied to the prediction of the mitigation limit of a breaking system subject to friction-induced vibrations coupled to two NESs, and next an airfoil model undergoing an aeroelastic instability coupled to a NESs setup (from one to four) is discussed. Theoretical results are compared, for validation purposes, to direct numerical integration of the system. The comparisons show good agreement.

**Keywords:** Multi-degree-of-freedom unstable system, set of nonlinear energy sinks, passive mitigation, relaxation oscillations, mitigation limit, asymptotic analysis

## A.1 INTRODUCTION

In the framework of passive vibration control, it is now recognized that the nonlinear absorbers based on the concept of Targeted Energy Transfer (TET) also known as Nonliner Energy Sink (NES) are good candidate to consider especially when low frequency and high level are concerned. A basic NES generally consists of a light mass, an essentially nonlinear spring and a viscous linear damper. Due to the essentially nonlinear stiffness (with no linear part), a NES can tune to any frequency content displayed by the primary vibrating structure even if its mass is small. The TET concept results from nonlinear mode bifurcations and it can be described as a 1:1 resonance capture [76, 203]. A first review of the state-of-the-art can be found in [205] and more recently in [129] as a part of the state-of-the-art about technologies for nonlinear dissipative devices including for example (without being exhaustive) polynomial, non-polynomial, piecewise linear and non-smooth stiffness NES, and also vibro-impact NES, track NES, bistable NES and magnetic NES. The generation of new technologies is always in progress. Recent studies include for example a novel prototype of NES enhanced by an inerter for vibration reduction [223, 225], an acoustic NES [85] and a virbo-acoustic NES with a controlled acoustic load [33] for noise attenuation and a grounded NES for rotor system vibration suppression [221].

NESs are also used to control dynamic instabilities. The possible suppression of the limit cycles of a Van der Pol oscillator utilizing a NES was first demonstrated numerically in [122] and next considered theoretically in [75]. A series of papers was also dedicated to aeroelastic mitigation [78, 123, 131, 208], helicopter groud resonance instability mitigation [a10, a12] and mode coupling instability mitigation in a friction system [a14]. As a recent works, aeroelastic suppression of an airfoil with control surface using a nonlinear energy sink is
considered in [93]. Panel flutter suppression with NES is discussed in [156]. Note that in these two recent works only complete Limit Cycle Oscillations (LCO) suppression or make the LCOs less intense are considered.

A single NES (coupled to the primary system) was the NES configuration the most often studied. Multi-Degree-Of-Freedom (multi-DOF) NES in series [83, 87, 193, 200] or in parallel [a7, 30, 126, 208, 224] were also considered. The use of several NESs allow to increase the range of excitations (level and frequency) within which the NES can absorb energy efficiently and/or to extract vibration energy simultaneously from multiple linear modes of the primary system.

To capture the complexity of the dynamic of the coupled system (primary system + NES) due to the nonlin-ear behaviour of the NES, various analysis have been developed including direct time integration method [221, 224], harmonic balance method [222, 225] or complex-averaging method [223]. The complex-averaging me-thod [135] gives also access to the slow-flow model of the nonlinear system. This slow-flow is governed by different time scales and analytical methods such as multiple scale analysis [144], or Geometric Singular Per-turbation Theory (GSPT) [107] can be used to extract the complex dynamics. In [188], for example harmonic forcing was considered where response regimes are characterized in terms of periodic and relaxation oscilla-tions of the slow-flow also named Strongly Modulated Responses (SMR). An alternative method known as mixed multiple scale/harmonic balance method has been proposed by Luongo and Zulli [130] leading to an equivalent result without the complexification-averaging being necessary. This approach was used by the same authors [131] in case of NES control of aeroelastic instability.

To go further, a discussion on the relationship between the dimensionality of the critical manifold, the nature and the distribution of the fixed points of the slow dynamics and the observed response regimes is explored in [74] where it is apparent that the analysis is difficult for systems with large dimension. Hence, it seems interesting to investigate the possibility of implementing a procedure able to deal with a multi-DOF system coupled to a large number of NESs.

In this work, auto-oscillating systems are considered and without loss of generality we focus on multi-degree-of-freedom unstable vibrating systems undergoing cubic nonlinearities and having only one mode which can become unstable through Hopf bifurcation. Our objective is to predict the steady-state regimes (i.e com-plete suppression of instability, mitigation through periodic response or mitigation through SMR) when the (primary) system is coupled to a family of single NESs. The proposed procedure starts with the diagonaliza-tion of the primary underlying linear system written in the state-space form using the so-called biorthogonal transformation. Then, the dynamics of the diagonalized system is reduced ignoring the stable modes and keeping only the unstable mode of the primary system. The complexification-averaging method is next ap-plied to the resulting coupled system, leading to slow-flow of the system. Based on the GSPT, it is shown that the critical manifold of the system can be reduced to a one dimensional parametric curve evolving in a multidimensional space. This form of critical manifold is similar to that obtained considering a network of parallel single NESs [a7]. Knowing the stability properties of the critical manifold and the fixed point of the slow-flow (position and stability), the response regimes can be predicted. Finally as in [a7] we introduce and predict the mitigation limit of the set of NESs, defined as the value of a chosen bifurcation parameter which separates harmless situations (corresponding to responses in which the set of NESs acts) from harmful situations (corresponding to responses in which the set of NESs does not act).

The system under study is introduced in the next section starting from the physical model, applying successively a rescaling step, a diagonalization step and a reduction step. In Sect. A.3, the slow-flow dynamics is first obtained applying the complexification-averaging method and then analyzed using the geometric singular perturbation theory. Sects. A.3.1 to B.3.4 describe in detail the mathematical developments which allow to obtain the critical manifold and the fixed points (position and stability) of the slow-flow. From the knowledge of that, the Sect. A.3.5 provides the procedure to predict the steady-state response regimes and consequently the mitigation limit. Finally, in Sect. D.5, two examples are considered. In Sect. A.4.1, the method is first applied to the prediction of the mitigation limit of a breaking system subject to friction-induced vibrations coupled to two NESs. A airfoil model undergoing an aeroelastic instability coupled to a set of NESs is next considered in Sect. A.4.2.

## A.2 THE SYSTEM UNDER STUDY

### A.2.1 The initial model

One considers in this work the following primary nonlinear system

$$ \tilde{M}\ddot{\tilde{x}} + \tilde{C}\dot{\tilde{x}} + \tilde{K}\tilde{x} + \tilde{g}^{\text{NL}}(\tilde{x}) = 0, \quad (A.1) $$
where $\tilde{\mathbf{x}} = (\tilde{x}_1, \dots, \tilde{x}_N)^T$ with $(\cdot)^T$ denotes the transpose operator, the dot represents the differentiation with respect to the time $t$ (already assumed dimensionless) and

$$ \tilde{\mathbf{M}} = \begin{bmatrix} \tilde{m}_{1,1} & \dots & \tilde{m}_{1,N} \\ \vdots & \ddots & \vdots \\ \tilde{m}_{N,1} & \dots & \tilde{m}_{N,N} \end{bmatrix}, \quad \tilde{\mathbf{C}} = \begin{bmatrix} \tilde{c}_{1,1} & \dots & \tilde{c}_{1,N} \\ \vdots & \ddots & \vdots \\ \tilde{c}_{N,1} & \dots & \tilde{c}_{N,N} \end{bmatrix} \quad \text{and} \quad \tilde{\mathbf{K}} = \begin{bmatrix} \tilde{k}_{1,1} & \dots & \tilde{k}_{1,N} \\ \vdots & \ddots & \vdots \\ \tilde{k}_{N,1} & \dots & \tilde{k}_{N,N} \end{bmatrix} \quad (A.2a) $$

are constant matrices. The nonlinear vector-valued function $\tilde{\mathbf{g}}^{\text{NL}}$ is called nonlinearity of the primary system. We assume that each component of $\tilde{\mathbf{g}}^{\text{NL}}$ is a linear combination of monomial terms of order 3.

We also assume that Eq. (D.1) can undergo a single dynamic instability of the trivial solution and that the nonlinear function $\tilde{\mathbf{g}}^{\text{NL}}$ allows the existence of Limit Cycle Oscillations (LCOs) on which the system can saturate.

System (D.1) can undergo a single dynamic instability of the trivial solution. In order to protect the primary system from this instability, $M$ purely cubic ungrounded NESs with masses $\tilde{m}_{h_m}$, damping coefficients $\tilde{c}_{h_m}$ and cubic stiffness $\tilde{\lambda}_{h_m}$ ($m = 1, \dots, M$) are used. Taking into account the NESs displacements $\tilde{h}_m(t)$ ($m = 1, \dots, M$), the equations of the whole coupled system result on the following form

$$ \tilde{\mathbf{M}}\ddot{\mathbf{x}} + \tilde{\mathbf{C}}\dot{\mathbf{x}} + \tilde{\mathbf{K}}\tilde{\mathbf{x}} + \tilde{\mathbf{g}}^{\text{NL}}(\tilde{\mathbf{x}}) + \tilde{\mathbf{B}}(\operatorname{diag}(\tilde{c}_{h_m})(\mathbf{T}\dot{\mathbf{x}} - \dot{\tilde{\mathbf{h}}}) + \operatorname{diag}(\tilde{\lambda}_{h_m})\mathbf{f}^{\text{NL}}(\mathbf{T}\tilde{\mathbf{x}} - \tilde{\mathbf{h}})) = \mathbf{0} \quad (A.3a) $$

$$ \operatorname{diag}(\tilde{m}_{h_m})\ddot{\mathbf{h}} - (\operatorname{diag}(\tilde{c}_{h_m})(\mathbf{T}\dot{\mathbf{x}} - \dot{\tilde{\mathbf{h}}}) + \operatorname{diag}(\tilde{\lambda}_{h_m})\mathbf{f}^{\text{NL}}(\mathbf{T}\tilde{\mathbf{x}} - \tilde{\mathbf{h}})) = \mathbf{0}, \quad (A.3b) $$

where $\tilde{\mathbf{h}} = (\tilde{h}_1, \dots, \tilde{h}_M)^T$, and the constant matrices

$$ \tilde{\mathbf{B}} = \begin{bmatrix} \tilde{b}_{1,1} & \dots & \tilde{b}_{1,M} \\ \vdots & \ddots & \vdots \\ \tilde{b}_{N,1} & \dots & \tilde{b}_{N,M} \end{bmatrix} \quad \text{and} \quad \mathbf{T} = \begin{bmatrix} t_{1,1} & \dots & t_{1,N} \\ \vdots & \ddots & \vdots \\ t_{M,1} & \dots & t_{M,N} \end{bmatrix} \quad (A.4) $$

are the influence coefficient matrices (depending on the position of the NESs). The nonlinear vector function $\mathbf{f}^{\text{NL}}$, due to the NESs attachments, is given by

$$ \mathbf{f}^{\text{NL}}(\mathbf{y}) = (y_1^3, \dots, y_M^3)^T. \quad (A.5) $$

Finally $\operatorname{diag}(d_m)$ denotes the diagonal matrix with the main diagonal $(d_1, d_2, \dots, d_M)$.

System (D.2) is the system under study (a $N$ multi-DOFs the primary system coupled with $M$ one-DOF NESs).

We now assume that the masses, $\tilde{m}_{h_m}$, and the damping coefficients, $\tilde{c}_{h_m}$, of the $M$ NESs are small and we introduce a small, dimensionless parameter $\epsilon$ ($0 < \epsilon \ll 1$) and the associated rescaled coefficients, $m_{h_m}$ and $c_{h_m}$, as

$$ \tilde{m}_{h_m} = \epsilon m_{h_m} \quad \text{and} \quad \tilde{c}_{h_m} = \epsilon c_{h_m}, \quad \text{for } m = 1, \dots, M. \quad (A.6) $$

Then, rescaling the variables $\tilde{\mathbf{x}}$ and $\tilde{\mathbf{h}}$ through $\epsilon$ as

$$ \mathbf{x} = \frac{\tilde{\mathbf{x}}}{\sqrt{\epsilon}} \quad \text{and} \quad \mathbf{h} = \frac{\tilde{\mathbf{h}}}{\sqrt{\epsilon}} \quad (A.7) $$

and inserting the rescaled parameters $m_{h_m}$ and $c_{h_m}$ and the rescaled variables $\mathbf{x}$ and $\mathbf{h}$ into (D.2) lead to

$$ \tilde{\mathbf{M}}\ddot{\mathbf{x}} + \tilde{\mathbf{C}}\dot{\mathbf{x}} + \tilde{\mathbf{K}}\mathbf{x} + \epsilon\tilde{\mathbf{g}}^{\text{NL}}(\mathbf{x}) + \epsilon\tilde{\mathbf{B}}\operatorname{diag}(m_n)\ddot{\mathbf{h}} = \mathbf{0} \quad (A.8a) $$

$$ \epsilon \operatorname{diag}(m_{h_m})\ddot{\mathbf{h}} - (\epsilon \operatorname{diag}(c_{h_m})(\mathbf{T}\dot{\mathbf{x}} - \dot{\mathbf{h}}) + \epsilon \operatorname{diag}(\tilde{\lambda}_{h_m})\mathbf{f}^{\text{NL}}(\mathbf{T}\mathbf{x} - \mathbf{h})) = \mathbf{0}. \quad (A.8b) $$

Finally multiplying (A.8a) by $\tilde{\mathbf{M}}^{-1}$ and (A.8b) by $\operatorname{diag}(1/m_{h_m})$ we obtain

$$ \ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} + \epsilon\mathbf{g}^{\text{NL}}(\mathbf{x}) + \epsilon\mathbf{B}\ddot{\mathbf{h}} = \mathbf{0} \quad (A.9a) $$

$$ \epsilon \ddot{\mathbf{h}} - (\epsilon \operatorname{diag}(\gamma_{h_m})(\mathbf{T}\dot{\mathbf{x}} - \dot{\mathbf{h}}) + \epsilon \operatorname{diag}(\alpha_{h_m})\mathbf{f}^{\text{NL}}(\mathbf{T}\mathbf{x} - \mathbf{h})) = \mathbf{0}, \quad (A.9b) $$

with

$$ \mathbf{C} = \tilde{\mathbf{M}}^{-1}\tilde{\mathbf{C}}, \quad \mathbf{K} = \tilde{\mathbf{M}}^{-1}\tilde{\mathbf{K}}, \quad \mathbf{B} = \tilde{\mathbf{M}}^{-1}\tilde{\mathbf{B}}\operatorname{diag}(m_n) $$

and $\gamma_{h_m} = c_{h_m}/m_{h_m}$, $\alpha_{h_m} = \tilde{\lambda}_{h_m}/m_{h_m}$ and $\mathbf{g}^{\text{NL}}(\mathbf{x}) = \tilde{\mathbf{M}}^{-1}\tilde{\mathbf{g}}^{\text{NL}}(\mathbf{x})$.
### A.2.2 Model reduction

In this section, the system (D.3) is reduced taking into account that the primary system (D.1) can undergo only a single instability of the trivial solution.

First, to simplify asymptotic analysis in next section, it is convenient to introduce new coordinates as

$$ \mathbf{v} = \mathbf{x} + \epsilon \mathbf{B} \mathbf{h} \qquad (A.10a) $$

$$ \mathbf{w} = \mathbf{h} - \mathbf{T} \mathbf{x}, \qquad (A.10b) $$

giving reciprocally,

$$ \mathbf{x} = (\mathbf{I}_N + \epsilon \mathbf{B} \mathbf{T})^{-1}(\mathbf{v} - \epsilon \mathbf{B} \mathbf{w}) \qquad (A.11a) $$

$$ \mathbf{h} = (\mathbf{I}_M + \epsilon \mathbf{T} \mathbf{B})^{-1}(\mathbf{w} + \epsilon \mathbf{T} \mathbf{v}), \qquad (A.11b) $$

where $\mathbf{I}_n$ is the identity matrix of size $n$.

Using Eqs. (D.4) with the inverse matrices $(\mathbf{I}_N + \epsilon \mathbf{B} \mathbf{T})^{-1}$ and $(\mathbf{I}_M + \epsilon \mathbf{T} \mathbf{B})^{-1}$, approximated respectively by $(\mathbf{I}_N - \epsilon \mathbf{B} \mathbf{T})$ and $(\mathbf{I}_M - \epsilon \mathbf{T} \mathbf{B})$, Eq. (D.3) is transformed into the following form

$$ \ddot{\mathbf{v}} + \mathbf{C} \dot{\mathbf{v}} + \mathbf{K} \mathbf{v} - \epsilon \left[ \mathbf{C} \mathbf{B} (\mathbf{T} \dot{\mathbf{v}} + \dot{\mathbf{w}}) + \mathbf{K} \mathbf{B} (\mathbf{T} \mathbf{v} + \mathbf{w}) \right] + \epsilon \mathbf{g}^{\mathrm{NL}}(\mathbf{v}) = \mathbf{0} \qquad (A.12a) $$

$$ \epsilon \left[ \mathbf{T} \mathbf{B} \operatorname{diag}(\gamma_{h_m}) \dot{\mathbf{w}} + \mathbf{T} \mathbf{C} \mathbf{B} (\mathbf{T} \dot{\mathbf{v}} + \dot{\mathbf{w}}) + \mathbf{T} \mathbf{K} \mathbf{B} (\mathbf{T} \mathbf{v} + \mathbf{w}) + \mathbf{T} \mathbf{B} \operatorname{diag}(\alpha_{h_m}) \mathbf{f}^{\mathrm{NL}}(\mathbf{w}) \right] = \mathbf{0}, \qquad (A.12b) $$

where only the first order term in $\epsilon$ has been retained in the nonlinearity of the primary underline linear system.

To capture the essential features of the single instability, model reduction is only performed on Eq. (D.5a), keeping the whole second equation. To achieve this, Eq. (D.5a), is written in state-space form as follows

$$ \dot{\mathbf{y}} = \mathbf{A} \mathbf{y} + \epsilon [\mathbf{D}_1 \mathbf{y} + \mathbf{D}_2 \mathbf{w} + \mathbf{D}_3 \dot{\mathbf{w}}] + \epsilon \mathbf{D}_4 \mathbf{g}^{\mathrm{NL}}(\mathbf{y}) \qquad (A.13) $$

where $\mathbf{y} = (v_1, \dots, v_N, \dot{v}_1, \dots, \dot{v}_N)^T$, $\mathbf{g}^{\mathrm{NL}}$ trivially redefined as $\mathbf{g}^{\mathrm{NL}}(\mathbf{y}) = \mathbf{g}^{\mathrm{NL}}(\mathbf{v})$,

$$ \mathbf{A} = \left[ \begin{array}{c|c} \mathbf{0} & \mathbf{I}_N \\ -\mathbf{K} & -\mathbf{C} \end{array} \right], \qquad (A.14) $$

$$ \mathbf{D}_1 = \left[ \begin{array}{c|c} \mathbf{0} & \mathbf{0} \\ \mathbf{K} \mathbf{T} \mathbf{B} & \mathbf{C} \mathbf{T} \mathbf{B} \end{array} \right], \quad \mathbf{D}_2 = \left[ \begin{array}{c} \mathbf{0} \\ \mathbf{K} \mathbf{B} \end{array} \right], \quad \mathbf{D}_3 = \left[ \begin{array}{c} \mathbf{0} \\ \mathbf{C} \mathbf{B} \end{array} \right] \quad \text{and} \quad \mathbf{D}_4 = \left[ \begin{array}{c} \mathbf{0} \\ -\mathbf{I}_N \end{array} \right]. \qquad (A.15) $$

The matrix $\mathbf{A}$ characterizes the linear dynamics of the primary system.

Because the matrix $\mathbf{A}$ is not a symmetric matrix, its diagonalization must necessarily be carried out using biorthogonality property of the right eigenvectors $\mathbf{r}_i$ ($i = 1, \dots, 2N$) and the left eigenvectors $\mathbf{l}_j$ ($j = 1, \dots, 2N$) of $\mathbf{A}$ hereafter defined. The method is briefly recalled in this section. The general procedure is provided for example in [140].

Let us consider the right and left eigenvector matrices

$$ \mathbf{R} = [\mathbf{r}_1 \ \mathbf{r}_1^* \ \dots \ \mathbf{r}_N \ \mathbf{r}_N^*] \quad \text{and} \quad \mathbf{L} = [\mathbf{l}_1 \ \mathbf{l}_1^* \ \dots \ \mathbf{l}_N \ \mathbf{l}_N^*] \qquad (A.16) $$

and the associated eigenvalue diagonal matrix $\mathbf{\Lambda} = \operatorname{diag}(\lambda_1, \lambda_1^*, \dots, \lambda_N, \lambda_N^*)$ corresponding to the following right and left eigenvalue problems

$$ \mathbf{A} \mathbf{R} = \mathbf{R} \mathbf{\Lambda} \quad \text{and} \quad \mathbf{A}^T \mathbf{L} = \mathbf{L} \mathbf{\Lambda}. \qquad (A.17) $$

The right and left eigenvectors satisfy the biorthogonality properties stating that $\mathbf{L}^T \mathbf{R}$ and $\mathbf{L}^T \mathbf{A} \mathbf{R}$ are diagonal matrices. It is convenient to normalize the two sets of eigenvectors in order to obtain

$$ \mathbf{L}^T \mathbf{R} = \mathbf{I}_N, \qquad (A.18) $$

that results in

$$ \mathbf{L}^T \mathbf{A} \mathbf{R} = \mathbf{\Lambda}. \qquad (A.19) $$

The biorthogonal transformation consists in introducing the biorthogonal coordinates which are constituted of $N$ pairs of complex conjugates, $q_n$ and $q_n^*$ ($n = 1, \dots, N$), and defined by the following relations

$$ \mathbf{y} = \mathbf{R} \mathbf{q} \Leftrightarrow \mathbf{q} = \mathbf{L}^T \mathbf{y}, \qquad (A.20) $$
where

$$ \mathbf{q} = (q_1, q_1^*, \dots, q_N, q_N^*)^T. \qquad (A.21) $$

Introducing Eq. (D.8) in Eq. (D.6), the equations of motion take the form of the following system

$$ \dot{\mathbf{q}} = \mathbf{\Lambda}\mathbf{q} + \epsilon\mathbf{L}^T[\mathbf{D}_1\mathbf{R}\mathbf{q} + \mathbf{D}_2\mathbf{w} + \mathbf{D}_3\dot{\mathbf{w}}] + \epsilon\mathbf{L}^T\mathbf{D}_4\mathbf{g}^{\text{NL}}(\mathbf{R}\mathbf{q}), \qquad (A.22) $$

and therefore Eqs. (D.5) are equivalent to

$$ \begin{aligned} & \dot{\mathbf{q}} - \mathbf{\Lambda}\mathbf{q} - \epsilon\mathbf{L}^T[\mathbf{D}_1\mathbf{R}\mathbf{q} + \mathbf{D}_2\mathbf{w} + \mathbf{D}_3\dot{\mathbf{w}}] - \epsilon\mathbf{L}^T\mathbf{D}_4\mathbf{g}^{\text{NL}}(\mathbf{R}\mathbf{q}) = 0 \\ & \ddot{\mathbf{w}} + \text{diag}(\gamma_{h_m})\dot{\mathbf{w}} - \text{TC}(\mathbf{R}_{\text{dl}}\mathbf{q}_u + \mathbf{R}_{\text{dr}}\mathbf{q}_d) - \text{TK}(\mathbf{R}_{\text{ul}}\mathbf{q}_u + \mathbf{R}_{\text{ur}}\mathbf{q}_d) + \\ & \text{diag}(\alpha_{h_m})\mathbf{f}^{\text{NL}}(\mathbf{w}) + \epsilon[\text{TB diag}(\gamma_{h_m})\dot{\mathbf{w}} + \text{TCB}(\mathbf{T}(\mathbf{R}_{\text{dl}}\mathbf{q}_u + \mathbf{R}_{\text{dr}}\mathbf{q}_d) + \dot{\mathbf{w}}) + \\ & \qquad \text{TKB}(\mathbf{T}(\mathbf{R}_{\text{ul}}\mathbf{q}_u + \mathbf{R}_{\text{ur}}\mathbf{q}_d) + \mathbf{w}) + \text{TB diag}(\alpha_{h_m})\mathbf{f}^{\text{NL}}(\mathbf{w})] = \mathbf{0}, \end{aligned} \qquad (A.23b) $$

where the matrix **R** (respectively the vector **q**) has been split into a $N \times N$-block matrix (respectively $N \times 1$-block vector) as

$$ \mathbf{R} = \begin{bmatrix} \mathbf{R}_{ul} & \mathbf{R}_{ur} \\ \mathbf{R}_{dl} & \mathbf{R}_{dr} \end{bmatrix} \quad \text{and} \quad \mathbf{q} = (\mathbf{q}_{u}^{T}, \mathbf{q}_{d}^{T})^{T} \qquad (A.24) $$

respectively.

From now we assume that one (and only one) mode (without loss of generality the first mode) of the primary system become unstable through Hopf bifurcation, when the chosen bifurcation parameter, denoted $\sigma$, crosses a particular parameter value $\sigma_0$ called *bifurcation point*, i.e

* if $\sigma < \sigma_0$, Re{$\{\lambda_n\} < 0 \forall i \in [1, N]$} and the trivial solution of the primary system is stable,

* if $\sigma > \sigma_0$, Re{$\{\lambda_1\} > 0$, Re{$\{\lambda_n\} < 0 \forall i \in [2, N]$}} and the trivial solution of the primary system is unstable.

Because we assume that only the first pair of eigenvalues ($\lambda_1, \lambda_1^*$) may have positive real part, after some exponentially decaying transients the components $q_n$ and $q_n^*$ ($n = 2, \dots, N$) become small. Therefore, all terms related to $q_n$ and $q_n^*$ ($n = 2, \dots, N$) are omitted from further consideration and Eq. (D.9) is reduced to

$$ \dot{q}_1 = \lambda_1 q_1 + \epsilon l_1^T (\mathbf{D}_1(\mathbf{r}_1 q_1 + \mathbf{r}_1^* q_1^*) + \mathbf{D}_2 \mathbf{w} + \mathbf{D}_3 \dot{\mathbf{w}}) + \epsilon l_1^T \mathbf{D}_4 \mathbf{g}^{\text{NL}} (\mathbf{r}_1 q_1 + \mathbf{r}_1^* q_1^*), \qquad (A.25) $$

in which we stated that $q_n = q_n^* = 0$ ($n = 2, \dots, N$).

Finally, grouping Eqs. (D.12) and (D.5b), we obtain the reduced system

$$ \begin{aligned} & \dot{q}_1 - (\epsilon\rho + j\omega)q_1 - \epsilon l_1^T (\mathbf{D}_1(\mathbf{r}_1 q_1 + \mathbf{r}_1^* q_1^*) + \mathbf{D}_2 \mathbf{w} + \mathbf{D}_3 \dot{\mathbf{w}}) - \epsilon l_1^T \mathbf{D}_4 \mathbf{g}^{\text{NL}} (\mathbf{r}_1 q_1 + \mathbf{r}_1^* q_1^*) = 0 \\ & \ddot{\mathbf{w}} + \text{diag}(\gamma_{h_m})\dot{\mathbf{w}} - \text{TC}(\mathbf{r}_1^d q_1 + \mathbf{r}_1^d q_1^*) - \text{TK}(\mathbf{r}_1^u q_1 + \mathbf{r}_1^u q_1^*) + \text{diag}(\alpha_{h_m})\mathbf{f}^{\text{NL}}(\mathbf{w}) + \\ & \epsilon(\text{TB diag}(\gamma_{h_m})\dot{\mathbf{w}} + \text{TCB}(\mathbf{T}(\mathbf{r}_1^d q_1 + \mathbf{r}_1^d q_1^*) + \dot{\mathbf{w}}) + \text{TKB}(\mathbf{T}(\mathbf{r}_1^u q_1 + \mathbf{r}_1^u q_1^*) + \mathbf{w}) + \\ & \qquad \text{TB diag}(\alpha_{h_m})\mathbf{f}^{\text{NL}}(\mathbf{w})) = \mathbf{0} \end{aligned} \qquad (A.26b) $$

where the vector $\mathbf{r}_1$ has been split as $\mathbf{r}_1 = (\mathbf{r}_1^u, \mathbf{r}_1^d)^T$.

Note that $\lambda_1$ has also been replaced by

$$ \lambda_1 = \epsilon\rho + j\omega, \qquad (A.27) $$

with $j^2 = -1$ assuming a weak instability i.e. that Re{$\{\lambda_1\}$} is in the order of $O(\epsilon)$ (Re{$\{\lambda_1\} = \epsilon\rho$} with $\rho$ is in the order of $O(1)$).

## A.3 ASYMPTOTIC ANALYSIS

We focus on the motion of the system in the vicinity of a 1 : 1 resonance associated with a frequency close to the natural frequency ($\omega$ see Eq. (A.27)) of the primary system.
### A.3.1 The slow-flow

The slow-flow describes the evolution of the amplitude and phase of the system which evolve slowly compare to the oscillations at the frequency $\omega$. To obtain it, we use the Complexification-Averaging method introduced by [135] and discussed in detail in [205].

The first step is the complexification<fnref n="1" /> which consists in introducing the following complex vector

$$ \zeta = \dot{\mathbf{w}} + j\omega\mathbf{w} \qquad (A.28) $$

or equivalently

$$ \mathbf{w} = \frac{\zeta - \zeta^*}{2j\omega}, \qquad \dot{\mathbf{w}} = \frac{\zeta + \zeta^*}{2} \qquad \text{and} \qquad \ddot{\mathbf{w}} = \dot{\zeta} - \frac{j\omega}{2}(\zeta + \zeta^*) \qquad (A.29) $$

and expressing the complex variables $q_1$ and $\zeta$ as

$$ q_1 = \phi e^{j\omega t} \quad \text{and} \quad \zeta = \xi e^{j\omega t} \qquad (A.30) $$

where the scalar $\phi$ and the vector $\xi$ are the complex (assumed slow modulated) amplitude of the fast component $e^{j\omega t}$.

The second step is the averaging which consists in, first, substituting Eqs. (D.16) and (D.17) into Eqs. (A.26), and next, averaging the resulting equations over one period $T = 2\pi/\omega$. We obtain the slow-flow described by the following differential equations

$$ \dot{\phi} = \epsilon f(\phi, \xi), \qquad (A.31a) $$

$$ \dot{\xi} = g(\phi, \xi, \epsilon) \qquad (A.31b) $$

where the functions $f$ and $g$ are defined by

$$ f(\phi, \xi) = (\rho + l_1^T D_1 r_1 + l_1^T D_4 c^{NL}(r_1) |\phi|^2) \phi + \left( \frac{1}{2\omega j} l_1^T D_2 + \frac{1}{2} l_1^T D_3 \right) \xi \qquad (A.32a) $$

$$ = (A_0 + A_1 |\phi|^2) \phi + A_1^T \xi, \qquad (A.32b) $$

where the vector $c^{NL}(r_1)$ is obtained from the nonlinear function $g^{NL}(r_1)$ replacing each monomial $r_{1i}r_{1j}r_{1k}$ by $r_{1i}^*r_{1j}r_{1k} + r_{1i}r_{1j}^*r_{1k} + r_{1i}r_{1j}r_{1k}^*$ (where $r_{1i}$ is the $i$-th component of $r_1$) and

$$ \mathbf{g}(\phi, \xi, \epsilon) = (\mathbf{T}\mathbf{K}\mathbf{r}_1^u + \mathbf{T}\mathbf{C}\mathbf{r}_1^d)\phi - \frac{1}{2}\mathbf{diag}\left(j\omega + \gamma_{h_m} - j\frac{3\alpha_{h_m}}{4\omega^3}|\xi_m|^2\right)\xi \\ - \epsilon\left[\left(\mathbf{T}\mathbf{K}\mathbf{B}\mathbf{T}\mathbf{r}_1^u + \frac{1}{2}\mathbf{T}\mathbf{C}\mathbf{B}\mathbf{T}\mathbf{r}_1^d\right)\phi \right. \\ \left. + \frac{1}{2}\left(\mathbf{T}\mathbf{C}\mathbf{B} + \frac{1}{j\omega}\mathbf{T}\mathbf{K}\mathbf{B} + \mathbf{T}\mathbf{B}\mathbf{diag}\left(\gamma_{h_m} + j\frac{3\alpha_{h_m}}{4\omega^3}|\xi_m|^2\right)\right)\right] \xi, \qquad (A.33a) $$

$$ = \mathbf{g}_0(\phi, \xi) + \epsilon \mathbf{g}_1(\phi, \xi). \qquad (A.33b) $$

The scalars $A_0$ and $A_1$, the vector $A_1$ and the vector functions $g_0$ and $g_1$ have been introduced to simplify the analysis in next Section.

Finally, it is convenient to rewrite Eqs. (A.31) considering real variables. To achieve this, polar coordinates are introduced as

$$ \phi = se^{j\delta}, \qquad (A.34a) $$

$$ \xi = \mathrm{diag}(e^{j\theta_m}) \mathbf{r}, \qquad (A.34b) $$

where $\mathbf{r} = (r_1, \dots, r_M)^T$ and $s$ and $\delta$ (respectively $r_m$ and $\theta_m$) characterized the modulus and the argument of $\phi$ (respectively $\xi_m$). Substituting Eqs. (A.34) into Eqs. (A.31) and separating in real and imaginary parts, Eqs. (A.31) reduce to

$$ \dot{s} = \epsilon F(s, \mathbf{r}, \vartheta) \qquad (A.35a) $$

<footnote n="1">1. This step is not necessary for the variable $q_1$ because it is already a complex variable.</footnote>
$$ \dot{\mathbf{r}} = \mathcal{G}(s, \mathbf{r}, \vartheta, \epsilon) \qquad (A.35b) $$

$$ \dot{\vartheta} = \mathcal{H}(s, \mathbf{r}, \vartheta, \epsilon), \qquad (A.35c) $$

where the argument differences $\vartheta_m = \theta_m - \delta$ have been considered and written in vector form as $\vartheta = (\vartheta_1, \dots, \vartheta_M)^T$. The functions $\mathcal{F}$, $\mathcal{G}$ and $\mathcal{H}$ are deduced from functions $f$ and $g$ in Eqs. (A.31) as

$$ \mathcal{F} = \mathrm{Re} \{ f e^{-j\delta} \}, \qquad (A.36a) $$

$$ \mathcal{G} = \mathrm{Re} \{ \mathrm{diag} (e^{-j\theta_m}) \mathbf{g} \}, \qquad (A.36b) $$

$$ \mathcal{H} = \mathrm{Im} \left\{ \mathrm{diag} \left( \frac{e^{-j\theta_m}}{r_m} \right) \mathbf{g} \right\} - \epsilon \mathrm{Im} \left\{ \mathbf{U} \frac{f e^{-j\delta}}{s} \right\} \qquad (A.36c) $$

where $\mathbf{U} = (1, \dots, 1)^T$ is a vector of size $M$.

### A.3.2 The Critical Manifold

The slow-flow described by Eqs. (A.35) has a slow-fast nature. Indeed, Eq. (A.35a) contains only $\mathcal{O}(\epsilon)$ terms and Eqs. (A.35b) and (A.35c) contain both $\mathcal{O}(1)$ and $\mathcal{O}(\epsilon)$ terms. Consequently $s$ is the slow variable whereas the vectors $\mathbf{r}$ and $\vartheta$ contain the fast variables $r_m$ and $\vartheta_m$ ($m = 1, \dots, M$).

Equations (A.35) can be reformulated by switching from the fast time scale $t$ to the slow time scale $\tau = \epsilon t$ as

$$ s' = \mathcal{F}(s, \mathbf{r}, \vartheta) \qquad (A.37a) $$

$$ \epsilon \mathbf{r}' = \mathcal{G}(s, \mathbf{r}, \vartheta, \epsilon) \qquad (A.37b) $$

$$ \epsilon \vartheta' = \mathcal{H}(s, \mathbf{r}, \vartheta, \epsilon), \qquad (A.37c) $$

where $\epsilon' = \frac{d}{d\tau}$. Solutions of the slow-fast system (A.35) (or (A.37)) can exhibit slow and fast epochs characterized by the speed at which the solution advances.

Stating $\epsilon = 0$, the following subsystems are derived from (A.35) and (A.37) respectively

$$ \dot{s} = 0 \qquad (A.38a) $$

$$ \dot{\mathbf{r}} = \mathcal{G}(s, \mathbf{r}, \vartheta, 0), \qquad (A.38b) $$

$$ \dot{\vartheta} = \mathcal{H}(s, \mathbf{r}, \vartheta, 0), \qquad (A.38c) $$

which is the *fast subsystem*, and

$$ s' = \mathcal{F}(s, \mathbf{r}, \vartheta) \qquad (A.39a) $$

$$ 0 = \mathcal{G}(s, \mathbf{r}, \vartheta, 0), \qquad (A.39b) $$

$$ 0 = \mathcal{H}(s, \mathbf{r}, \vartheta, 0), \qquad (A.39c) $$

which is the *slow subsystem*.

The algebraic equations (A.39b) and (A.39c) of the slow subsystem define the so-called *Critical Manifold S* [107] as

$$ S := \left\{ (s, \mathbf{r}, \vartheta) \in \mathbb{R}^{2M+1} \mid \mathcal{G}(s, \mathbf{r}, \vartheta, 0) = 0 \text{ and } \mathcal{H}(s, \mathbf{r}, \vartheta, 0) = 0 \right\} \qquad (A.40) $$

corresponding to the subspace where the slow motions take place and also to the fixed points of the fast subsystem (A.38).

Looking at Eqs. (A.33b) and (A.36), one can deduce that

$$ \mathcal{G}(s, \mathbf{r}, \vartheta, 0) = \mathrm{Re} \{ \mathrm{diag}(e^{-j\theta_m}) \mathbf{g}_0 \} \quad \text{and} \quad \mathcal{H}(s, \mathbf{r}, \vartheta, 0) = \mathrm{Im} \left\{ \mathrm{diag} \left( \frac{e^{-j\theta_m}}{r_m} \right) \mathbf{g}_0 \right\}, \qquad (A.41) $$

and after some algebraic manipulation, Eqs (A.39b) and (A.39c) can be written in the reduced complex form as

$$ \mathrm{diag}(e^{-j\vartheta_m}) (\mathbf{T}\mathbf{K}\mathbf{r}_1^u + \mathbf{T}\mathbf{C}\mathbf{r}_1^d) s - \frac{1}{2} \mathrm{diag}\left(j\omega + \gamma_{h_m} - j\frac{3\alpha_{h_m}}{4\omega^3}r_m^2\right) \mathbf{r} = 0. \qquad (A.42) $$

Introducing the complex functions $F_m$ ($m = 1, \dots, M$) of a real variable as

$$ F_m(x) = \frac{\omega}{B_m} \left( 1 - j \frac{\gamma_{h_m}}{\omega} - \frac{3\alpha_{h_m}}{4\omega^4} x^2 \right) = R_m(x) + j I_m(x), \qquad (A.43) $$
where

$$B_m = -2j (\mathbf{T}\mathbf{K}\mathbf{r}_1^u + \mathbf{T}\mathbf{C}\mathbf{r}_1^d)_m \quad (A.44)$$

with $(.)_m$ stands for the $m$-th coordinate of the vector $(.)$, the Critical Manifold $S$ is characterized by

$$s^2 = H_m(r_m), \quad (m = 1, \dots, M), \qquad (A.45a)$$

$$\vartheta_m = -\arg(F_m(r_m)), \quad (m = 1, \dots, M) \qquad (A.45b)$$

where the real functions $H_m$ ($m = 1, \dots, M$) of a real variable are defined as

$$H_m(x) = x^2 [R_m(x)^2 + I_m(x)^2]. \qquad (A.46)$$

It is interesting to note that the general form (E.48) of the Critical Manifold $S$ is the same as found in [a7] in the case of a one-DOF unstable primary system coupled to several parallel NESs. Indeed, in the latter case, the equivalent functions $F_m$ can be obtained from (A.43) stating $B_m = 1$. From this form, as suggested in [a7], specific properties of the Critical Manifold can be deduced which can be advantageously used to characterize the response regimes of the system under study.

### A.3.3 Critical Manifold properties and stability

Due to Eqs. (E.48), the critical manifold $S$ appears as a one dimensional parametric curve evolving in $\mathbb{R}^{+M+1}$ as a solution of the $M$ nonlinear equations (E.48a) with respect to the $M+1$ unknown $(r_1, \dots, r_M, s)$. A classical continuation method can be used to obtain $S$. Equivalently $S$ can be obtained solving the following differential-algebraic equations

$$\sum_{m=1}^{M} \left( \frac{dr_m}{du}(u) \right)^2 + \left( \frac{ds}{du}(u) \right)^2 = 1, \quad \text{for } u_1 < u < u_2 \qquad (A.47a)$$

$$s(u)^2 - H_m(r_m(u)) = 0, \quad (m = 1, \dots, M) \qquad (A.47b)$$

$$\frac{ds}{du}(0) = 1, \quad \frac{dr_m}{du}(0) = 0, \quad (m = 1, \dots, M) \qquad (A.47c)$$

$$s(0) = s^0, \quad r_m(0) = r_m^0, \quad (m = 1, \dots, M) \qquad (A.47d)$$

where $u$ is arclength variable in the range $u_1(< 0)$ to $u_2(> 0)$ and $(r_1^0, \dots, r_M^0, s^0) \in \mathbb{R}^{+M+1}$ is a root of the algebraic equations (A.47b), chosen to be close to trial solution 0. An example of typical critical manifold for a system containing two NESs are plotted Fig. A.1. $S$ starts at the origin by a straight line followed by a finite number of successive simple curved lines curving in opposite directions (with respect to $s$) and ended by an unbounded straight line. Each simple curve is characterized by a curving point defined by a horizontal tangency in at least one $(r_m, s)$-plane (i.e $ds/dr_m = 0$).

The critical manifold $S$ can also be viewed as the intersection of the $M$ manifolds

$$M_m = \{(\mathbf{r}, s) \in \mathbb{R}^{+M+1} \mid s^2 = H_m(r_m)\}, \quad (m = 1, \dots, M) \qquad (A.48)$$

in the $(r_1, \dots, r_M, s)$-space.

A typical projection of $M_m$ onto the subspace $\mathbb{R}^2$ corresponding to the variables $(r_m, s)$ is depicted in Fig. A.2. Its geometrical form is characterized by the properties of the function $H_m$. From Eq. (A.46), it can be shown that, for each $m$, if the following relation holds

$$\gamma_{h_m} < \frac{\omega}{\sqrt{3}}, \qquad (A.49)$$

then $H_m$ admits a local maxima $r_m^{\max}$ and a local minima $r_m^{\min}$ given as the positive roots of its derivative $dH_m/dx$. An easy calculus shows that

$$r_m^{\max} = \frac{2}{3}\omega^2 \sqrt{\frac{2 - \sqrt{1 - 3\left(\frac{\gamma_{h_m}}{\omega}\right)^2}}{\alpha_{h_m}}} \quad \text{and} \quad r_m^{\min} = \frac{2}{3}\omega^2 \sqrt{\frac{2 + \sqrt{1 - 3\left(\frac{\gamma_{h_m}}{\omega}\right)^2}}{\alpha_{h_m}}} \qquad (A.50)$$

with $r_m^{\max} < r_m^{\min}$. Each extrema, defined by $(r_m^p, s_m^p)$ with $s_m^p = \sqrt{H_m(r_m^p)}$ where the upper-script $p$ stands for min or max (see Fig. A.2 where condition (E.50) is satisfied) is associated to a curving point of $S$.
<figure><img src="image_64.png" /><figcaption>FIGURE A.1. (a) Example of typical critical manifold S of a system containing two parallel NESs ($N = 2$) in the ($r_1, r_2, s$)-space (i.e. $\mathbb{R}^{+3}$). S (in red) is the intersection between the blue surface $s = \sqrt{H_1(r_1)}$ and the green surface $s = \sqrt{H_2(r_2)}$. The parameters of the primary system are those used in Sect. A.4.2 resulting in $\omega = 0.82$. The NESs parameters are: $\gamma_{h_1} = 0.3$, $\gamma_{h_2} = 0.2$, $\alpha_{h_1} = \alpha_{h_2} = 6$, $B_1 = 0.028 - j0.923$ and $B_2 = 0.030 - j0.819$ (corresponding to the second and the third NESs in the example of Sect. D.5). (b) Projection of the CM in the ($r_1, r_2$)-plane.</figcaption></figure>

As already mentioned, each point of S is also a fixed point for the fast subsystem (A.38) and the stability can be determined by examining the sign of the eigenvalue real parts of the Jacobian matrix of the differential system (A.38b)(A.38c). Hence the stability range of S is defined by the points ($r_1, \dots, r_N, s$) on S that satisfy

$$ \forall m \in [1, M], \quad \frac{dH_m}{dx}(r_m) > 0. \qquad (A.51) $$

Stability condition can be also written equivalently in terms of arclength variable u at the point

$$ (r_1(u), \dots, r_N(u), s(u)) $$

of S as

$$ \forall m \in [1, M], \quad \frac{ds}{du}(u) \left( \frac{dr_m}{du}(u) \right)^{-1} > 0. \qquad (A.52) $$

We assume from now that Eq. (E.50) is satisfied for the M NESs.

For each $m$ ($1 \le m \le M$), the extrema points ($r_m^{\text{max}}, s_m^{\text{max}}$) and ($r_m^{\text{min}}, s_m^{\text{min}}$) characterize two bounds on S where S ceases to be hyperbolic<fnref n="2" /> connecting stable or attractive (continuous line) and unstable or repulsive (dashed line) parts of S (see Fig. A.2). It follows that ($r_m^{\text{max}}, s_m^{\text{max}}$) (respectively ($r_m^{\text{min}}, s_m^{\text{min}}$)) defines a point $\mathbf{f}^{\text{ma},m}$ (respectively $\mathbf{f}^{\text{mi},m}$) on S called *fold point* by

$$ \mathbf{f}^{\text{ma},m} = (r_1^{\text{ma},m}, \dots, r_{m-1}^{\text{ma},m}, r_m^{\text{max}}, r_{m+1}^{\text{ma},m}, \dots, r_M^{\text{ma},m}, s_m^{\text{max}}) \qquad (A.53a) $$

$$ (\text{respectively } \mathbf{f}^{\text{mi},m} = (r_1^{\text{mi},m}, \dots, r_{m-1}^{\text{mi},m}, r_m^{\text{min}}, r_{m+1}^{\text{min},m}, \dots, r_M^{\text{min},m}, s_m^{\text{min}})) \qquad (A.53b) $$

with $r_n^{\text{ma},m}$ (respectively $r_n^{\text{mi},m}$) satisfied for $n = 1, \dots, M, n \ne m$, $H(r_n^{\text{ma},m}) = H(r_m^{\text{max}})$ (respectively $H(r_n^{\text{mi},m}) = H(r_m^{\text{min}})$). The number of fold points is equal to $2M$ (2 by NES).

<footnote n="2">2. S is hyperbolic if all eigenvalues of the Jacobian matrix of the differential systems (A.38b)(A.38c) have nonzero real part.</footnote>
<figure><img src="image_65.png" /><figcaption>**FIGURE A.2.** Part of the *Critical Manifold* (CM) corresponding to the *m*-th NES, Eq (A.46). The parameters of the primary system are those used in Sect. A.4.2 resulting in $\omega = 0.82$. The NESs parameters are: $\gamma_{h_m} = 0.2$, $\alpha_{h_m} = 6$, $B_m = 0.030 - j0.819$ (corresponding to the third NES in the example of Sect. D.5).</figcaption></figure>

Moreover in each $(r_m, s)$-plane, the dynamics of any point not on $S$ is entirely defined by the direction of the fast variable $r_m$ as given by Eq. (A.38a) (i.e. $s = \text{Cte}$) and which is indicated by the horizontal arrows in Fig. A.2. Hence at a fold point $f^{\text{ma},m}$ (respectively $f^{\text{mi},m}$) on $S$, a jump occurs (fast dynamic) to a point on $S$ named arrival point and denoted $j^{\text{ma},m}$ (respectively $j^{\text{mi},m}$). During the jump, $s = s_m^{\text{max}}$ and for $n = 1, \dots, M, n \neq m$ the trajectory lies in the manifold $M_n$ leading to

$$j^{\text{ma},m} = (r_1^{\text{ma},m}, \dots, r_{m-1}^{\text{ma},m}, r_m^{\text{up}}, r_{m+1}^{\text{ma},m}, \dots, r_M^{\text{ma},m}, s_m^{\text{max}}) \quad (\text{A.54a})$$

$$(\text{respectively } j^{\text{mi},m} = (r_1^{\text{mi},m}, \dots, r_{m-1}^{\text{mi},m}, r_m^{\text{do}}, r_{m+1}^{\text{mi},m}, \dots, r_M^{\text{mi},m}, s_m^{\text{min}})) \quad (\text{A.54b})$$

with the ordinate $r_m^{\text{up}}$ (respectively $r_m^{\text{do}}$) is obtained solving $H_m(r_m^{\text{up}}) = H_m(r_m^{\text{max}})$ (respectively $H_m(r_m^{\text{do}}) = H_m(r_m^{\text{min}})$) that leads to

$$r_m^{\text{up}} = \frac{2\sqrt{2}}{3}\omega^2 \sqrt{\frac{1 + \sqrt{1 + 3\left(\frac{\gamma_{h_m}}{\omega}\right)^2}}{\alpha_{h_m}}} \quad (\text{respectively } r_m^{\text{do}} = \frac{2\sqrt{2}}{3}\omega^2 \sqrt{\frac{1 - \sqrt{1 - 3\left(\frac{\gamma_{h_m}}{\omega}\right)^2}}{\alpha_{h_m}}}). \quad (\text{A.55})$$

It follows that the stability domain of $S$ is only characterized in the $(r_1, r_2, \dots, r_M)$-space $(\mathbb{R}^{+M})$ by the subspace $D$ given in terms of Cartesian product as

$$D = \prod_{m=1}^{M} D_m, \quad (\text{A.56})$$

where

$$D_m = [0 \ r_m^{\text{max}}] \cup [r_m^{\text{min}} + \infty) \quad (\text{A.57})$$

denotes the stability range associated to the *m*-th NES (see Fig. A.2).

Eq. (A.56) may be expanded as the union of $2^M$ disjunct Cartesian products $I_k$ (with $k = 1, \dots, 2^M$) as

$$D = \bigcup_{k=1}^{2^M} I_k, \quad (\text{A.58})$$

where $I_k$ appears as the product of close or left-close intervals resulting in an isolated stable or attractive part of the critical manifold $S$ as

$$I_k = \prod_{m \in C_k} [0, r_m^{\text{max}}] \prod_{m \in C_k^c} [r_m^{\text{min}}, +\infty), \quad (\text{A.59})$$
with $C_k$ is a part of the finite set $\{1, 2, \dots, M\}$ and $C_k^c$ its complementary. The first and the last Cartesian products are chosen as $I_1 = \prod_{m=1}^M [0, r_m^{\max}]$ and $I_{2M} = \prod_{m=1}^M [r_m^{\min}, +\infty)]$ respectively. With this choice, the first branch of $S$ is defined on $I_1$ whereas the last branch of $S$ is defined on $I_{2M}$ (see Fig. A.1(b)).

Finally a more compact representation of the CM can be used obtained projecting $S$ in the $(r_1, r_2, \dots, r_M)$-space resulting in a one dimensional parametric curve in $\mathbb{R}^{+M}$. The intersection of this curve with $D$ defines the stable zone of the CM which can be decomposed as successive branches, one branch by subset $I_k$ named hereafter $S_k$. The first branch $S_1$ starts at origin of $\mathbb{R}^{+M}$ denotes $\mathbf{b}_1^1 = \mathbf{0}$ and ends at the fold point $\mathbf{b}_1^2 = \mathbf{f}^{\text{ma}, \bar{m}}$ where $\bar{m} = \text{argmin}_{1 \le m \le M} r_m^{\max}$ (i.e $\bar{m}$ is the subscript value at which $r_m^{\max}$ is minimal). The last branch $S_{2M}$ starts at the fold point $\mathbf{b}_{2M}^1 = \mathbf{f}^{\text{mi}, \bar{m}}$ where $\bar{m} = \text{argmax}_{1 \le m \le M} r_m^{\min}$ (i.e $\bar{m}$ is the subscript value at which $r_m^{\min}$ is maximal) and it is unbounded (i.e $\mathbf{b}_{2M}^2 = \infty$). The other branches $S_k$ are ended by two fold points defined by $\mathbf{b}_k^1 = \mathbf{f}^{\text{mi}, \bar{m}}$ where $\bar{m} = \text{argmax}_{m \in C_k^c} r_m^{\min}$ and $\mathbf{b}_k^2 = \mathbf{f}^{\text{ma}, \bar{m}}$ where $\bar{m} = \text{argmin}_{m \in C_k} r_m^{\max}$. Note that for each of the points $\mathbf{b}_k^1$ and $\mathbf{b}_k^2$ (except for $\mathbf{b}_k^1$ and $\mathbf{b}_{2M}^2$) we can also associated the arrival points (following the jump) named $\mathbf{a}_k^1$ and $\mathbf{a}_k^2$ defined by $\mathbf{a}_k^1 = \mathbf{j}^{\text{ma}, \bar{m}}$ and $\mathbf{a}_k^2 = \mathbf{j}^{\text{mi}, \bar{m}}$ with the corresponding value of $\bar{m}$.

Such a multi-S-shape of the critical manifold $S$ provides a possibility for relaxation oscillations [86] of the slow-flow characterized by fast transitions (jumps) of the dynamics during each cycle (the possible fast transitions between the two stable branches of $S$ are denoted by arrows in Fig. A.2). Because several NESs are considered, complex relaxation oscillations scenarios may be contemplated. Such relaxation oscillations of the slow-flow explain the existence of *Strongly Modulated Responses* [75, 78, 188] (SMRs) of the non-averaged system (A.26).

We proceed to a more detailed analysis of the possible steady-state regimes in next section computing fixed points of the slow-flow and their stability.

### A.3.4 The fixed points of the slow-flow

Following the Geometric Singular Perturbation Theory (GSPT) [54, 62, 107], the nontrivial fixed points of the slow-flow (A.31) can be approximated computing the fixed points of Eq. (A.39).

Substituting first Eq. (E.48a) with $n=1$ (any $n \in [1, M]$ can be chosen) into Eq. (A.39a), we obtain

$$ \frac{1}{2\sqrt{H_1(r_1)}} \frac{dH_1}{dx}(r_1) r_1' = \mathcal{F}\left(\sqrt{H_1(r_1)}, r_1, \dots, r_M, \vartheta_1, \dots, \vartheta_M\right) \quad (A.60) $$

which can be rewritten, substituting Eq. (A.45b) into the resulting equation, as

$$ \frac{dH_1}{dx}(r_1) r_1' = f_{r_1}(r_1, \dots, r_M), \quad (A.61) $$

where

$$ f_{r_1}(r_1, \dots, r_M) = 2\sqrt{H_1(r_1)} \mathcal{F}\left(\sqrt{H_1(r_1)}, r_1, \dots, r_M, -\arg(F_1(r_1)), \dots, -\arg(F_M(r_M))\right). \quad (A.62) $$

The function $f_{r_1}$ will be explicitly given here after. Note that $f_{r_1}$ is a function of a real variable (here $r_1$) because all variables $r_n$ (with $n \in [2, M]$) are linked to $r_1$ ($H_n(r_n) = H_1(r_1)$, see Eq. (E.48a)).

From Eq. (B.40), it is possible to detect fixed points as

$$ f_{r_1}(r_1, \dots, r_M) = 0, \quad H_m(r_m) = H_1(r_1) \quad (m=2, \dots, M) \quad (A.63a) $$

$$ \frac{dH_1}{dx}(r_1) \neq 0 \quad (A.63b) $$

and folded singularities as

$$ f_{r_1}(r_1, \dots, r_M) = 0, \quad H_m(r_m) = H_1(r_1) \quad (m=2, \dots, M) \quad (A.64a) $$

$$ \frac{dH_1}{dx}(r_1) = 0. \quad (A.64b) $$

Folded singularities correspond to situations for which fixed points and fold points coincide, they are hints of possible *canard explosions* [16].

The remaining of the section is dedicated to the analysis of the fixed points which, for $0 < \epsilon \ll 1$, are assumed to be the fixed points of (A.35).

Combining Eq. (A.62) with Eqs. (E.48), Eqs. (A.43) and Eq. (A.46), and using the notations introduced Eq. (A.32b), the function $f_{r_1}$ reduces to

$$ f_{r_1}(r_1, \dots, r_M) = 2(A_0^R + A_1^R H_1(r_1)^2) H_1(r_1) + 2 \sum_m^M r_m^2 (A_{1m}^R R_m(r_m) + A_{1m}^I I_m(r_m)), \quad (A.65) $$