# Basins of attraction in a dynamical system with a time-varying control parameter: the case of attack transients in a simple model of reed musical instrument

Soizic Terrien, Baptiste Bergeot, Christophe Vergez, Samy Missoum

► To cite this version:

Soizic Terrien, Baptiste Bergeot, Christophe Vergez, Samy Missoum. Basins of attraction in a dynamical system with a time-varying control parameter: the case of attack transients in a simple model of reed musical instrument. Journal of Sound and Vibration, 2025, 618, pp.119241. 10.1016/j.jsv.2025.119241. hal-05264904

HAL Id: hal-05264904

https://hal.science/hal-05264904v1

Submitted on 24 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Basins of attraction in a dynamical system with a time-varying control parameter: the case of attack transients in a simple model of reed musical instrument

Soizic Terrien<sup>a</sup>, Baptiste Bergeot<sup>b</sup>, Christophe Vergez<sup>c</sup>, Samy Missoum<sup>d</sup>

<sup>a</sup>Laboratoire d'Acoustique de l'Université du Mans (LAUM), UMR 6613,
Institut d'Acoustique - Graduate School (IA-GS), CNRS, Le Mans Université,
Avenue Olivier Messiaen, 72085 Le Mans, France

<sup>b</sup>INSA CVL, Univ. Orléans, Univ. Tours, LaMé EA 7494,
3 Rue de la Chocolaterie, CS 23410, 41034 Blois, France

<sup>c</sup>Aix Marseille Univ, CNRS, Centrale Med, LMA UMR7031,
4 impasse Nikola Tesla, 13013 Marseille, France

<sup>d</sup>Department of Aerospace and Mechanical Engineering, University of Arizona,
Tucson, AZ 85721, United States of America

Abstract

Multistability, a common phenomenon in dynamical systems in general and in musical instruments in particular, corresponds to the coexistence of several stable regimes for a given set of parameters. When a system is multistable, which regime is observed depends on the initial conditions on the state variables, and therefore relates to the basins of attraction. Here we consider a simple model of reed musical instrument. A bifurcation analysis shows that the system is bistable on a range of the blowing pressure, which is one of the main control parameters for the musician. In the bistability region, the boundary between basins of attraction is computed for the classical case of constant control parameters in the one hand, and for saturating ramps of the blowing pressure on the other hand. These profiles of the blowing pressure are chosen as prototypical representations of the so-called attack transients of the blowing pressure, during which the musician starts to blow in the instrument to create a musical note from silence. Attack transients are known to be of crucial importance in a musical context, yet mostly overlooked in physical models of the musician-instrument system. The separatrix between basins of attraction is computed numerically by combining time-domain simulations with a machine learning technique relying on support vector machine. Our results demonstrate that attack transients of the blowing pressure significantly affect the practical *observability* of the different stable regimes, suggesting that they might be used by musicians to navigate between basins. From a more general point of view, these results relate to rate-induced tipping phenomena and might therefore be of interest beyond the particular case of musical instruments.

# 1. Introduction

Self-sustained musical instruments, such as bowed string instruments, flutes, brass instruments or reed instruments are complex dynamical systems, known to produce a wide diversity of dynamical regimes. This includes equilibrium solutions where no sound is produced, different periodic regimes corresponding to different musical notes [1, 2, 3, 4], but also quasiperiodic and chaotic regimes [5, 6, 7, 8]. The existence and stability properties of these regimes depend on both design
parameters, which are fixed at the making stage by the instrument maker, and control parameters which are adjusted continuously by the musician while playing.

Similarly to many dynamical systems in applications ranging from biology to population dynamics, neurosciences, optics or climate models [9], musical instruments have been shown to display multistability, that is to say coexistence of several stable regimes (attractors) for identical parameters. In the context of multistability, which regime is observed depends on the initial conditions on the system variables. From a dynamical systems point of view, this relates to the basins of attraction of the different stable regimes, defined as the set of initial conditions for which the system converges towards a particular regime after any transient dynamics has died out; see for example [10]. In practice, basins of attraction relate to the likelihood to observe a given regime [11, 12], and therefore in a musical context to the *playability* of sound regimes: although stable, a regime with small basin can be very difficult to observe as this requires to set very specific initial conditions.

In the last decades, basins of attraction have been investigated in many different dynamical systems. This demonstrated that even simple systems can display a particularly complex structure of basins of attraction [13, 14, 15]. In parallel, different numerical methods have been developed that allow to access the boundaries (separatrices) between basins. This includes automated methods based on a large number of simulations from different initial conditions [11] and advanced continuation methods to compute the separatrix between basins in the phase space, [16, 17]. The latter method has been applied mainly to two-dimensional separatrices (i.e. surfaces); see for example [18, 19, 20]. In the context of nonlinear dynamics of musical instruments, although several studies addressed the question of multistability [2, 21, 22, 12], basins of attraction are mentioned in very few studies [23, 1, 12] despite their crucial importance.

Here, basins of attraction are investigated in a simple model of single-reed instrument. A bifurcation analysis shows that this model is bistable on a range of values of the blowing pressure, which is one of the main control parameters for the musician. Classically, basins of attraction are investigated for fixed, constant values of the system parameters. Nevertheless, in the context of wind musical instruments, investigating the influence of a time-varying parameter on practical playability of the different regimes is particularly relevant. First, it is unlikely that musicians can control initial conditions on the different state variables independently. Secondly, attack transients of the control parameters, which correspond here to the transient dynamics of the blowing pressure to its target value, are known to be of crucial importance in a musical context [24]. Finally, recent experimental and numerical investigations have demonstrated the major role played by the dynamics of the blowing pressure on the oscillation threshold and regime change thresholds [25, 26, 22, 27], a manifestation of the general phenomenon of dynamic bifurcations [28]. Here, we investigate numerically in a simple model of a reed musical instrument the influence of attack transients of the mouth pressure on the observed regimes. We aim at comparing the basins of attraction estimated for a constant value of the blowing pressure chosen in a bistability region of the bifurcation diagram with those obtained when an attack transient of the blowing pressure is considered. More precisely, we use ramp-like profiles of the blowing pressure in which the blowing pressure linearly increases over a transient period until it reaches the value used when $\gamma$ is constant, where it saturates. In practice, we investigate if attack transients can make a particular regime easier or more difficult to observe than predicted by the basins investigated when the parameters are constant.

From a more general point of view, the question of the influence of a time-varying control parameter on the observability of coexisting stable regimes is of interest beyond the specific case
of musical instruments. In particular, so-called critical transitions and rate-induced tipping have attracted considerable attention in the recent years, mainly in the context of climate modeling [29, 30, 31]. This phenomenon corresponds to an abrupt regime shift caused by the rate of change of a parameter rather than by its absolute value and can occur, for example, when a parameter changes too quickly for the system to adapt, thus causing the system to leave the basin of attraction of the initial state, even if it remains stable. Importantly, such transitions cannot necessarily be fully explained by the classical bifurcation analysis and basins of attraction, which are static in nature. Taking into account possible tipping phenomena is crucial to understand and predict regime transitions in real-life systems, where parameters evolve in time. Nevertheless, to the best of our knowledge, this has been barely investigated in the context of engineering systems in general and acoustical systems in particular.

This article is organised as follows. The model equations are detailed in section 2.1. In section 2.2, a bifurcation analysis in the blowing pressure is performed. This shows that the system is bistable on a range of the parameter of interest. Because the considered system is low-dimensional, the boundary between the two basins of attraction for constant parameters is obtained directly from the bifurcation analysis, as discussed in section 2.3. These basins for constant parameters are compared with results of time-domain simulations performed for time varying profiles of the blowing pressure parameter. In particular, we investigate the basins of attraction of a model that embeds a transient of the blowing pressure. As detailed in section 3, the surface separating the basins in the case where a transient of the control parameter is taken into account is estimated in time-domain simulations using a machine learning technique based on an explicit parameter space decomposition and a support vector machine (SVM) classifier. This approach has been used already in musical acoustics to determine the parameters leading to quasiperiodic sound regimes in a minimal model of single reed instruments [6]. This method allows us to focus on the exploration of the phase space close to the boundary between basins and, as such, to avoid a time-consuming exploration of the space of possible initial conditions. The effect of linear saturating ramps of the blowing pressure on the boundary between basins of attraction is explored in section 4. This is done for different values of the slope of the ramp (corresponding to different attack transients duration) and different target values of the blowing pressure.

## 2. Problem statement

### 2.1. Background on the model of reed musical instrument

Sound production in self-sustained musical instruments results from the nonlinear coupling between an exciter and a resonator [32]. In reed musical instruments, such as clarinets or saxophones, the sound production mechanism can be described as follows. As represented schematically in Figure 1, the pressure difference between the musician's mouth (i.e. the blowing pressure) and the instrument mouthpiece results in an air flow that makes the reed oscillate. This modulates the air flow entering the instrument, which constitutes a source for the air column contained in the instrument. The acoustic field thus created sustains the reed motion close to the resonance frequencies of the instrument resonator.

The state-of-the-art model of reed instruments includes three main elements (see for example chapter 9 in [33]): a one degree-of-freedom oscillator accounting for the reed motion, a nonlinear function that links the difference of pressure between the instrument mouthpiece and the musician's mouth to the flow entering the instrument, and the acoustical input impedance of the instrument's resonator which describes the response of the resonator in the frequency domain.
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic representation of the exciter of a single-reed musical instrument, with $\tilde{h}(\tilde{t})$ the opening of the reed channel, $h_0$ its value at rest (i.e. for zero blowing pressure), $p_m$ the pressure in the musician's mouth, $\tilde{p}(\tilde{t})$ the acoustical pressure in the instrument mouthpiece and $\tilde{u}(\tilde{t})$ the modulated flow entering the instrument.</figcaption></figure>

In the literature, the reed is classically modelled by a mass - spring - damper oscillator [33]. Because its resonance frequency is large compared to the typical playing frequencies, inertia and damping can be neglected in first approximation [33]. The reed is therefore represented through its sole stiffness, and the reed displacement $\tilde{h}$ depends linearly on the pressure difference between the instrument mouthpiece and the musician's mouth:

$$K_r(\tilde{h}(\tilde{t}) - h_0) = \tilde{p}(\tilde{t}) - p_m, \qquad (1)$$

with $h_0$ the height of the reed channel at rest, $K_r$ the equivalent surfacic stiffness of the reed, $p_m$ the pressure in the musician's mouth and $\tilde{p}$ the pressure in the mouthpiece (see Figure 1). Note that the dimensional variables, including time, are denoted with a tilde in anticipation of the upcoming non-dimensionalization of the model. It is assumed here that any contribution of the jet slowdown to the pressure in the mouthpiece can be neglected in first approximation due to turbulent mixing at the reed channel exit.

Following Bernoulli law, the flow $\tilde{u}$ entering the mouthpiece is written as follows:

$$\tilde{u}(\tilde{t}) = \operatorname{sgn}(p_m - \tilde{p}(\tilde{t})) W \tilde{h}^+(\tilde{t}) \sqrt{\frac{2 | p_m - \tilde{p}(\tilde{t}) |}{\rho}}, \qquad (2)$$

with $W$ the width of the reed channel and $\rho$ the air density. As shown in Figure 1, the amplitude of the reed motion is limited by the instrument mouthpiece: when $h \le 0$, the reed channel is closed and the flow $\tilde{u}(\tilde{t})$ is zero. Following the so-called ghost reed model [34], this is accounted for in equation (2) by the positive part $\tilde{h}^+(\tilde{t})$ of $\tilde{h}(\tilde{t})$, defined by: $\tilde{h}^+(\tilde{t}) = \max(\tilde{h}(\tilde{t}), 0)$. The response of the resonator is assumed to be linear, and is modelled by the input impedance $Z_{in}$ of the resonator. It is written in the frequency domain as a (truncated) sum of $n$ resonance modes [35]:

$$\tilde{Z}_{in}(\tilde{\omega}) = \frac{\tilde{P}(\tilde{\omega})}{\tilde{U}(\tilde{\omega})} = Z_c \sum_{i=1}^{n} \frac{Z_i \frac{\omega_i}{Q_i} j \tilde{\omega}}{\omega_i^2 - \tilde{\omega}^2 + j \tilde{\omega} \frac{\omega_i}{Q_i}}, \qquad (3)$$

where $\tilde{P}(\tilde{\omega})$ and $\tilde{U}(\tilde{\omega})$ are the Fourier transforms of the pressure and flow variables $\tilde{p}(\tilde{t})$ and $\tilde{u}(\tilde{t})$,
respectively. Parameters $Z_i$, $\omega_i$ and $Q_i$ are the modal amplitude, the resonance angular frequency and the quality factor of the $i^{th}$ mode, respectively, and $Z_c = \frac{\rho c}{S}$ is the characteristic impedance of plane waves in the resonator. Here, $S$ in the cross-section of the resonator and $c$ the sound velocity.

We aim at investigating the simplest reed instrument model displaying multistable dynamics. The motivation is twofold. First, an in-depth investigation of basins of attraction becomes practically complicated in high-order dynamical systems, with a large number of initial conditions to test. Secondly, our purpose here is to illustrate the phenomenon of rate-induced tipping which, to the best of our knowledge, has been poorly explored in engineering and never in musical acoustics, as well as its possible implications for musical instruments playing. For sake of clarity, we aim at exploring the simplest configuration where only two different stable regimes (and thus two basins of attractions) coexist. Here, the dimension of the model relates directly to the number of acoustic resonance modes. We consider here a single acoustic mode (i.e. $n = 1$ in equation (3)). Although this is a drastic simplification, this single mode model has been shown to reproduce the main features of the dynamics of a clarinet, apart from the spectral content of periodic regimes which is significantly affected by higher-order resonance modes [33]. For the sake of numerical conditioning, system (1)-(3) is written in a dimensionless form by defining the following variables (see for example [6]):

$$ 
\begin{aligned}
p(\tilde{t}) &= \frac{\tilde{p}(\tilde{t})}{P_M}; \\
u(\tilde{t}) &= \frac{Z_c}{P_M} \tilde{u}(\tilde{t}); \\
h(\tilde{t}) &= \frac{\tilde{h}(\tilde{t})}{h_0}; \\
t &= \omega_1 \tilde{t};
\end{aligned} 
\quad (4) $$

with $P_M = K_r h_0$ the minimal value of the mouth pressure required to close the reed channel (see equation (1)). The dimensionless mouth pressure $\gamma$ and reed channel height $\zeta$ are the two main control parameters for the musician and are defined as follows:

$$ 
\begin{aligned}
\gamma &= \frac{p_m}{P_M}; \\
\zeta &= Z_c W h_0 \sqrt{\frac{2}{\rho P_M}}.
\end{aligned} 
\quad (5) $$

The flow equation (2) is written as follows with respect to the dimensionless variables:

$$ u(t) = \zeta \operatorname{sgn}(\gamma - p(t)) [1 - \gamma + p(t)]^+ \sqrt{|\gamma - p(t)|} \quad (6) $$

For sake of numerical conditioning, this equation is regularised to avoid non-differentiability [12]. First, the function $h^+(t)$ and the sign function are rewritten as follows:

$$ 
\begin{aligned}
h^+(t) &= \frac{h(t) + |h(t)|}{2}, \\
\operatorname{sgn}(x) &= \frac{x}{|x|},
\end{aligned} 
\quad (7) $$
and the absolute value function is regularised as follows:

$$|x| \approx \sqrt{x^2 + \eta}, \qquad (8)$$

with $\eta$ a small regularisation parameter. Following [34], we consider $\eta = 10^{-3}$ throughout this article. Overall, the regularised flow is written as follows:

$$u(t) \approx \frac{1}{2} \zeta(\gamma - p) \left[ (\gamma - p)^2 + \eta \right]^{-\frac{1}{4}} \left[ 1 - \gamma + p + \sqrt{(1 - \gamma + p)^2 + \eta} \right]. \qquad (9)$$

Using equation (3) with $n = 1$, and using the dimensionless time $t = \omega_1 \tilde{t}$, one obtains, for the dimensionless variables:

$$P(\omega) = \frac{Z_1 \frac{1}{Q_1} j \omega}{1 - \omega^2 + j \omega \frac{1}{Q_1}} U(\omega), \qquad (10)$$

where $\omega$ is the dimensionless angular frequency. Finally, through an inverse Fourier transform of (10), one obtains a system of two first-order ordinary differential equations for the dimensionless acoustical pressure $p(t)$ in the instrument mouthpiece and its first order time derivative $w(t)$:

$$\begin{aligned} \dot{p}(t) &= w(t); \\ \dot{w}(t) &= \frac{Z_1}{Q_1} \dot{u}(t) - \frac{1}{Q_1} w(t) - p(t). \end{aligned} \qquad (11)$$

where $u(t)$ is given by (9) and now indicates the derivative with respect to the dimensionless time $t$. Unless stated otherwise, the parameter values considered throughout the article are $\zeta = 0.1$, $Z_1 = 50$, $\omega_1 = 1440$, $Q_1 = 36.6$, as considered in [36]. Throughout this article, the bifurcation parameter is chosen as the dimensionless mouth pressure $\gamma$.

## 2.2. One-parameter bifurcation diagram

A bifurcation analysis of model (11) is performed using AUTO-07p [37], a software for numerical continuation. Figure 2 shows a bifurcation diagram of (11) representing the amplitude of the dimensionless acoustic pressure $p(t)$ along the different coexisting branches of solutions, with respect to the dimensionless blowing pressure $\gamma$. Blue and red curves represent stable and unstable solutions, respectively. This shows that an equilibrium (non-oscillating) solution exists on the whole represented range of $\gamma$. It is stable for the lower values of $\gamma$ up to $\gamma = 0.42$ where it destabilises in a supercritical Hopf bifurcation $H$, before stabilising back through a second (subcritical) Hopf bifurcation at $\gamma = 1$. The branch of stable periodic solutions that emerges at the supercritical Hopf bifurcation and the branch of unstable periodic solutions emerging from the subcritical Hopf bifurcation collide and disappear in a saddle-node bifurcation of periodic orbits, as indicated by a SN point in Figure 2.

Between the subcritical Hopf bifurcation and the saddle-node bifurcation of periodic orbits, i.e. for $\gamma \in [1; 1.32]$ (shaded area in Figure 2), the system shows bistability between a periodic solution (corresponding to a musical note) and an equilibrium solution (where no sound is produced). Which stable regime is observed in practice depends on the considered initial conditions $p(0) = p_0$ and $w(0) = w_0$ of the state variables $p$ and $w$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Bifurcation diagram of model (11), showing the amplitude of the acoustical pressure $p(t)$ with respect to the dimensionless mouth pressure $\gamma$. Stable and unstable solutions are represented by thick blue parts and thin red parts of the curves, respectively. Points H of Hopf bifurcation and point SN of saddle-node bifurcation of periodic solutions are represented by blue dots and red dots, respectively. The shaded area indicates the bistability region.</figcaption></figure>

### 2.3. Basins of attraction for constant parameters

We focus on this bistable region of the bifurcation diagram, and we investigate the basins of attraction. These are defined as the set of initial conditions, in the $(p, w)$-phase plane, leading to each solution in the long-term (i.e. after any transient dynamics has died out). We first consider here the basins of attraction of model (11) where no time evolution of $\gamma$ is taken into account.

Because the model is here two-dimensional, the boundary (separatrix) between the basins of the equilibrium and periodic solution is the unstable periodic orbit that emerges from the subcritical Hopf bifurcation point (Figure 2); see for example [38]. Figure 3(a) shows the phase portrait of (11) for $\gamma = 1.05$, as obtained as a direct product of the numerical bifurcation analysis performed with AUTO-07p. The closed blue curve and the blue dot represent the stable periodic solution and the stable equilibrium, respectively. The red closed curve is the unstable periodic solution which separates the basins of attraction of the two stable solutions. In practice, in time-domain simulation, the equilibrium solution is observed for initial conditions $(p(0); w(0))$ chosen inside the closed red curve [38, 39]. All other initial conditions lead to the stable periodic solution. This is illustrated in Figure 3(b), which shows simulated time series of the dimensionless acoustical pressure $p$ for three different sets of initial conditions $(p_0, w_0)$, which are indicated in panel (a) by stars with corresponding colors. Simulations are performed with the MATLAB solver ode45 [40] with a maximal (dimensionless) time step $T_e = \omega_1 \times 10^{-5}$. The time series in green is (top subpanel in Figure 3(b)) is simulated from initial conditions chosen inside the unstable periodic solution in panel (a). As expected, the system converges towards the equilibrium solution. On the other hand, for initial conditions chosen either between the unstable and stable periodic orbits or outside the stable periodic orbit in panel (a), the trajectory approaches the stable periodic solution, as shown by the yellow and purple curves (middle and bottom subpanels in Figure 3 (b)). Figure 3(c) shows the evolution of the two dimensional phase portrait with respect to the bifurcation parameter $\gamma$, thus illustrating the evolution of the basins separatrix within the bistability region. In particular, the size of the equilibrium basin relates directly to the amplitude of the unstable periodic solution.

From a practical point of view, properties of the basins of attraction, including their relative sizes and geometry, relate to the *observability* of the different regimes: if the initial conditions are chosen randomly, a stable regime with a small basin is less likely to be observed than a regime
with a large basin [11]. In the context of musical instruments, basins of attraction might relate to
the playability, that is to say to the fact that a regime is easy to play or not from the musician
point of view [39]. In Figure 3(a) for example, the equilibrium basin is small compared to the
periodic solution ones: in practice, observing the equilibrium regime requires a careful selection of
the initial conditions and would be less easily accessible for the musician. It is important to note
that in higher dimensional dynamical systems – for example a reed instrument model that would
include a larger number *n* of resonance modes in equation (3) – the basins boundary can be much
more complex [41, 15]. In the general case in particular, determining the basins separatrix from a
bifurcation analysis is not necessarily straightforward.

<figure><img src="image_4.png" /><figcaption>Figure 3: (a): Phase portrait of system (11) for $\gamma = 1.05$, showing a stable and an unstable periodic solutions (blue and red closed curves, respectively) and a stable equilibrium (blue dot). Colored stars indicate initial conditions considered in time-domain simulations. (b) Time series of the acoustical pressure $p$ obtained in time-domain simulations of (11), for $\gamma = 1.05$ and for the different sets of initial conditions indicated in panel (a). (c): Evolution of the 2D phase portrait of (11) in the plane $(p, w)$, with respect to the continuation parameter $\gamma$.</figcaption></figure>

2.4. *Effect of the dynamics of* $\gamma$

<figure><img src="image_5.png" /><figcaption>Figure 4: Time-domain simulations of model (11), with a linear saturating ramp with slope $s$ of the mouth pressure $\gamma$. The initial value of $\gamma$ is $\gamma_0 = 0.065$, and its target value is $\gamma_t = 1.05$, initial conditions on the state variables are $(p_0, w_0) = (0, 0.06)$. Two different values of the slope $s = 0.006$ (yellow curves) and $s = 0.01$ (green curves) are considered. The case of a constant $\gamma = 1.05$ is represented in blue for sake of comparison. (a) Temporal profiles of $\gamma$ and (b) simulated acoustical pressure $p$ with respect to (dimensionless) time $t$.</figcaption></figure>

The dynamics of wind musical instruments in general and reed instruments in particular is most
often investigated considering constant values of the blowing pressure, the main control parameter
for the musician.

Recent experimental, numerical and theoretical studies [25, 26, 22, 12, 27] have shown that a
time-varying blowing pressure can significantly affect the observed dynamics. This includes the
value of the parameter at which bifurcations are observed, but also which regime is observed, among the different stable solutions, after the value of $\gamma$ has reached its target value.

Here, we aim at investigating the effect, on the observed regime, of transient dynamics of $\gamma$ from an initial value to a target (fixed) value $\gamma_t$ in the bistability range of the bifurcation diagram (see Figure 2). From an application point of view, this relates to attack transients, which are known to be of crucial importance both for musicians and instrument makers. From a more general point of view, we aim here at paving the way towards a better understanding of the relation between basins of attraction and musical practice: it is likely that the musician does not control directly and independently the initial conditions on the different state variables (here, the modal pressures and their derivatives). Rather, we investigate whether or not the musician might use the transient dynamics of the blowing pressure as an additional control parameter to navigate between basins and obtain the desired sound regime.

The effect of transient dynamics of $\gamma$ on the observed regime is illustrated in Figure 4. This shows time-domain simulations of (11) for a fixed constant value of $\gamma$ (i.e. when no attack transient is considered), and for two linearly increasing and saturating ramps of $\gamma$, as defined by the following equation:

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

Here, $\gamma_0$ is the starting value of $\gamma$ at $t = 0$, $s$ is the slope of $\gamma$ on its transient part, $\gamma_t$ is the target (saturating) value of $\gamma$ and $t_c = \frac{\gamma_t - \gamma_0}{s}$ the rising time from $\gamma_0$ to $\gamma_t$. In Figure 4, $\gamma_0 = 0.065$, $\gamma_t = 1.05$, and the initial conditions $(p_0, w_0) = (0, 0.06)$ are in the equilibrium basin of the model with no time variation of $\gamma$ (see Figure 3). Two different values of $s$ are considered, corresponding to different attack transient durations. Figure 4(a) shows the time profiles of $\gamma$ and Figure 4(b) show the corresponding simulated time series of the acoustical pressure $p$.

As expected, for constant parameters (blue curves), the equilibrium solution is observed after the transient on the acoustical variable $p(t)$ has died out. For a time-varying $\gamma$, the system also converges towards the equilibrium regime for the shorter attack transient (green curves). For the longer attack transient (yellow curves), it settles on the stable periodic solution. This demonstrates that, for identical initial conditions and bounds for $\gamma$, changing the rate of increase $s$ can push the system to a different basin of attraction than the one predicted for a constant $\gamma$, resulting in the observation of a different regime. These observations can be interpreted as a rate-induced tipping phenomenon, where the observed dynamics depends on a non trivial manner on the rate of the time-varying parameter, and not only on its absolute value.

In fact, which regime is observed in the long term depends on which basin of attraction (as determined for constant parameters with no parameter transient) the system is when $\gamma$ stops evolving (i.e. when the target value $\gamma_t$ of $\gamma$ is reached). This is illustrated in Figure 5, which represents the phase portrait obtained for a fixed constant value $\gamma = 1.05$, superimposed with results of one thousand time-domain simulations. For all these simulations, the target value of $\gamma$ is $\gamma_t = 1.05$ but different starting values $\gamma_0$ and different transition times between $\gamma_0$ and $\gamma_t$ are considered. For each simulation, the values of the state variables when the saturating value $\gamma_t$ is reached is represented by a circle, whose color indicates the observed regime (yellow for the periodic solution and green for the equilibrium solution). This shows that different time profiles of $\gamma$ can bring the system to different locations of the phase portrait and, as such, to one or the other of the basins of attraction determined for a fixed, constant value $\gamma_t$ of the control parameter.
<figure><img src="image_6.png" /><figcaption>Figure 5: Phase portrait of (11) for $\gamma = 1.05$, showing a stable periodic solution (blue closed curve), an unstable periodic solution (red closed curve) and a stable equilibrium (blue dot). Results of time-domain simulations of system (11,12) are superimposed, for $\gamma_t = 1.05$ and different values of $s, \gamma_0, p_0$ and $w_0$: shown is, for each simulation, the value of $(p,w)$ at $t = t_c$ when $\gamma$ reaches its target value. Yellow circles indicate simulations converging to the stable periodic solution, and green dots indicate simulations converging to the stable equilibrium.</figcaption></figure>

### 3. Numerical Method: explicit design space decomposition and support vector machine

A systematic investigation of the geometry of the basins of attraction in the 3D space $(p_0, w_0, \gamma_0)$, for different values of $\gamma_t$ and $s$ would rely on a very large number of simulations and would therefore be practically difficult. An alternative approach is to determine, in the $(p_0, w_0, \gamma_0)$-space, the boundary between the basins of attraction. This is performed here using the Matlab toolbox CODES, based on an Explicit Design Space Decomposition (EDSD) technique, and on a support vector machine classifier (SVM), see for example [42, 43] and [44, 45] as entry points to the extensive literature. SVM is a machine learning technique which has been widely used for classification tasks. The underlying idea is to generate an explicit (linear or nonlinear) decision function that corresponds to the boundary between two different *classes*, for example two different regimes of a model, in a space of relevant parameters or, as here, in the space of initial conditions.

Here, the two classes correspond to the two stable regimes and the decision function is an estimate of the basins separatrix. The classification criterion relies on the value of the maximal amplitude $p_{max}$ of the simulated acoustical pressure $p(t)$ after any transient dynamics has died out: $p_{max} < 0.1$ for the equilibrium regime and $p_{max} > 0.1$ for the periodic regime. This criterion is system-dependent and has been chosen here based on the results of the bifurcation analysis detailed above.

An initial SVM decision function, corresponding to an approximation of the boundary between both classes in the design space, is constructed using a set of uniformly distributed *training samples*. For each sample (corresponding to a set of values in the space of initial conditions), the model is simulated to determine to which regime (which *class*) the sample belongs. The first decision function is subsequently refined using an adaptive sampling scheme: additional simulations are performed for well-chosen samples, selected in sparse regions of the considered 3D space, and in regions where the probability of misclassification by the current SVM decision function is the highest (i.e. on the estimated boundary). This is achieved by selecting a new sample (i.e. a new set
of initial conditions) on the previously estimated boundary while maximizing the minimum distance
from existing samples. A single additional simulation is performed for this set of initial conditions.
The sample is then classified according to the pressure-based criterion previously defined. The
decision function is then updated using this new training sample and the process repeats for the
next iteration. Following [42], the convergence of this iterative calculation process of the basins
boundary is checked based on a set of 100 × 5<sup>d</sup> uniformly distributed convergence points, with d = 3
the dimension of the design space. Importantly, no simulation is performed from these convergence
points and they are not taken into account to determine the boundary between basins of attraction.
Rather, at each step *k* of the iterative process, one determines to which class these convergence
points belong to based solely on the SVM decision function. The fraction of convergence points
that changed class compared to iteration *k* − 1 is calculated. The convergence criterion is based on
both the fraction ε<sub>1</sub> of convergence points that change class between two successive iterations, and
on its rate of change ε<sub>2</sub>. Convergence is considered to be reached when ε<sub>1</sub> and ε<sub>2</sub> are simultaneously
below small positive user-defined numbers, typically chosen as 10<sup>−8</sup>.

This method allows for a systematic investigation of the influence of a time-varying $\gamma$ on the
observed regime, for different slopes of $\gamma$ and target values $\gamma_t$. Throughout this article, the basins
of attractions observed when a time-variation of $\gamma$ is considered are compared to the basins of
attraction of model (11) which considers a constant value of $\gamma$. In practice, for each value of $\gamma_t$ and
$s$, the boundary between the basins of attraction is estimated in the ($p_0, w_0, \gamma_0$)-space through time-
domain simulations with a Matlab Runge-Kutta solver of 4th order [40], for two thousand sets of
(randomly chosen) initial conditions that form the training samples. The maximal (dimensionless)
time step $T_s = \omega_1 \times 10^{-5}$ used for simulations has been chosen so that the results are not affected
by a further decrease of the time step. The duration of each simulation depends on the considered
slope $s$ (and consequently on the duration of the transient of $\gamma$) such that the steady-state regime
is reached. One thousand additional simulations are performed for samples chosen in an adaptive
manner to refine the basins boundary. We verified that, after this refinement iterative process, the
convergence criterion is satisfied.

The result of the iterative refinement process is illustrated in Figure 6 for a target value $\gamma_t = 1.05$
and a slope $s = 0.006$. Figure 6(a) depicts the initial SVM decision function trained with two
thousand uniformly distributed samples in the ($\gamma_0, p_0, w_0$) space along with their corresponding
classes. Changes in the predicted boundary after 100 and 1000 additional adaptive samples are
also depicted in Figure 6(b) and (c). Finally, panel (c) shows the refined boundary between basins
after one thousand adaptive samples have been taken into account. Here, the chosen numbers of
training and adaptive samples bring the two convergence criteria discussed above to $\epsilon_1 = 10^{-8}$ and
$\epsilon_2 < 10^{-9}$, respectively. For further illustration of the method, the interested reader can refer to
[6, 46] where this approach is illustrated on related models of wind musical instruments.

4. Influence of a transient of $\gamma$ on the basins boundaries

4.1. Linear saturating profiles of $\gamma$

The simplest case of linear saturating ramps of $\gamma$, as illustrated above in Figures 5 and 6, is
explored further. In Figure 6(b), for initial conditions chosen inside the blue surface, the system
converges towards the equilibrium regime. On the contrary, for initial conditions chosen outside
the blue surface, the system converges towards the stable periodic regime. For sake of comparison,
the yellow surface shows the separatrix between the basins of attraction obtained for a constant
value $\gamma(t) = \gamma_t = 1.05 \ \forall t$, extended here over the whole range of $\gamma_0$. Said otherwise, the yellow
<figure><img src="image_7.png" /><figcaption>Figure 6: Estimated boundary between the basins of attraction (blue surfaces) of system (11 - 12) in the presence of an initial transient of $\gamma$, represented in the $(\gamma_0, p_0, w_0)$-space, for a saturating linear ramp of $\gamma$ with target value $\gamma_t = 1.05$ and slope $s = 0.006$. (a) Initial estimation with an SVM based on two thousand training samples, represented by dots whose color encodes the class (i.e. the basin of attraction) they belong to (red for the equilibrium and green for the periodic solution). (b) Refined SVM after one hundred additional adaptive samples represented by circles. (c) Final boundary with one thousand additional adaptive samples. The yellow surface shows the boundary between the basins of attraction obtained without transient of $\gamma$, for a constant value $\gamma = 1.05$, extended here over the whole range of $\gamma_0$ for sake of comparison.</figcaption></figure>

surface represents what the boundary between the basins computed in the presence of an initial transient of $\gamma$ would be if this transient of $\gamma$ (before the target value is reached) had no effect. This demonstrates that the transient dynamics of $\gamma$ considerably affects the shape of the boundary between basins and, as such, the observability of both regimes. In particular, for very small initial value of $\gamma_0$, it becomes impossible in practice to reach the periodic solution: no matter what the initial values of the state variables $p_0$ and $w_0$ are, the system always converges towards the equilibrium solution. On the other hand, for values of $\gamma_0$ between 0.1 and 0.9, it is considerably easier to reach the periodic solution compared to the case with a constant $\gamma$: here, the region of the $(p_0, w_0)$-plane of initial conditions on the state variables leading to the equilibrium is drastically reduced compared to the case with no transient of the parameter (constant $\gamma$).

### 4.1.1. Influence of the slope $s$

The influence of the slope $s$ is investigated for the same target value of the dimensionless mouth pressure ($\gamma_t = 1.05$) than considered in Figure 6. Figure 7 shows the boundary between the basins in the case of a time-varying mouth pressure for $s = 0.01$ (panel (a)), $s = 0.05$ (panel (b)), $s = 0.1$ (panel (c)) and $s = 0.5$ (panel (d)).

The extreme values $s = 0.006$ and $s = 0.5$ considered in Figures 6 and 7(d) were selected so that the transition times from an initial $\gamma$ of zero to the target value correspond to realistic durations of attack transients observed in musicians' performances [24, 47]. Figures 7 show that, in the case with initial transients of $\gamma$, the basins of attraction strongly depend on the slope $s$ of $\gamma$. In particular, as $s$ increases, the basins for the model with initial transient of $\gamma$ get closer and closer to the basins obtained without the transient, for a constant value of $\gamma$. For $s = 0.5$, in Figure 7(d), the basins obtained with and without initial transient of $\gamma$ are mostly similar, showing that the dynamics of $\gamma$ has very little effect on the observed regime in this case. Overall, the fact that the dynamics of $\gamma$ has or not an influence on the observed regime depends on the timescales
<figure><img src="image_8.png" /><figcaption>Figure 7: Boundary between the basins of attraction (blue surfaces) of (11,12), represented in the $(\gamma_0, p_0, w_0)$-space, for saturating linear ramps of $\gamma$ with target value $\gamma_t = 1.05$. The slope is $s = 0.01$ (panel (a)), $s = 0.05$ (panel (b)), $s = 0.1$ (panel (c)) and $s = 0.5$ (panel (d)).</figcaption></figure>

separation between the state variables and $\gamma$. More precisely, when $\gamma$ evolves quickly compared to the internal timescales of model (11), the state variables $p$ and $w$ of the original system have very little time to evolve during the attack transient of $\gamma$: when $\gamma$ reaches $\gamma_t$, the values of the state variables $p$ and $w$ are very close to their initial values $p_0$ and $w_0$. It is thus very unlikely for the system to be brought, during the transient of $\gamma$, in a different basin of attraction. Consequently, when an initial transient of $\gamma$ is considered and for a very large slope $s$, the basins of attraction are very close to the basins observed when no dynamics of $\gamma$ is taken into account. On the other hand, if the timescale on which $\gamma$ evolves is of the same order of magnitude (or slower) than the internal timescales of system (11), the state variables have time to evolve during the transient of $\gamma$. They can therefore be far from their initial values when $\gamma$ saturates. In this case, the system can be brought in a different basin of attraction during the transient of $\gamma$. This explains why the basins of attraction computed for a constant $\gamma$ equal to the target value $\gamma_t$ – which is used in the time-varying $\gamma$ case – are not able to predict which regime is observed in this latter case.

From this perspective, the shape of the basins of attraction in Figure 6 (for a small value of $s$) can be interpreted as follows. When $\gamma_0 \in [0.4; 1]$, that is to say in the range of $\gamma$ where the
equilibrium solution is unstable (as shown by the bifurcation diagram in Figure 2), the system evolves at first towards the stable periodic regime, which is the only stable solution, regardless the initial conditions $p_0$ and $w_0$ on the state variables. Consequently, when $\gamma$ reaches $\gamma_t$ (which is chosen here slightly above the value of $\gamma$ at which the system becomes bistable), most trajectories have reached the periodic regime, and stay on this regime. This explains why, in the case of a time-varying $\gamma$, the basin of the equilibrium solution is drastically reduced compared to the case of a constant $\gamma$.

## 4.2. Influence of the target value $\gamma_t$ of $\gamma$

In the previous section, the influence of the dynamics of $\gamma$ was investigated for a target value $\gamma_t = 1.05$ of $\gamma$. As shown in the bifurcation diagram in Figure 2, this is slightly above the value of $\gamma$ at which the system becomes bistable. This implies that, in the case of a constant $\gamma$, the basin of attraction of the equilibrium solution is significantly smaller than the one of the periodic solution. Moreover, in the case of a time-varying $\gamma$, as $\gamma$ increases from a low value to this target value, the periodic solution is the only stable solution during the major part of this transient (as the bifurcation diagram shows).

Figures 8 and 9 show the boundaries between the basins of attraction for the model with an initial transient of $\gamma$ (11-12) for saturating linear ramps of $\gamma$ with four different values of the slope $s$, for target values $\gamma_t = 1.1$ (Figure 8) and $\gamma_t = 1.2$ (Figure 9). These values are chosen within the range of $\gamma$ for which the system is bistable.

As shown above in Figure 3(c), in the case where $\gamma$ is constant, the basin of attraction of the equilibrium solution becomes larger as $\gamma$ is increased within the bistability region of the bifurcation diagram. Indeed, this basin is bounded by the unstable periodic solution in the $(p, w)$-plane of state variables (shown in red in Figure 3(c), whose amplitude increases monotonically between the Hopf bifurcation at $\gamma = 1$ (where the system becomes bistable) and the saddle-node bifurcation SN at $\gamma = 1.4$ after which the system is not bistable anymore.

The comparison between Figures 6, 7, 8 and 9 shows that the shape of the separatix between the basins of attraction, as well as the way it depends on the value of the slope $s$ of $\gamma$, are qualitatively similar for all the different target values $\gamma_t$ considered here. In particular, it is worth noting that, for very small values of $\gamma_0$, the basin of the equilibrium solution is always larger when there is an initial transient in $\gamma$ than when $\gamma$ is constant regardless the value of $\gamma_t$. On the other hand, the equilibrium basin is always smaller when there is an initial transient in $\gamma$ than when $\gamma$ is constant on a range of $\gamma_0 \in [0.1; \gamma_t]$. The only exceptions are for the larger values of the slope $s$, for which the basin in the case with initial transient of $\gamma$ tends to be identical to the one obtained for constant $\gamma$ regardless the value of $\gamma_0$, as explained in the previous section. Nevertheless, it can be noted that the size of the basin of the equilibrium obtained with a transient on $\gamma$ is related to the size of the corresponding basin computed for constant $\gamma$. This can be seen for example, by comparing the basins obtained for a given value of the slope $s$ but different target mouth pressures $\gamma_t$. In the case with constant blowing pressure, the equilibrium basin becomes larger when $\gamma$ increases within the bistability region of the bifurcation diagram (see Figure 3(c)). Similarly, in the case with transients of $\gamma$, when the target value of the blowing pressure $\gamma_t$ increases (while staying in the bistability region of the bifurcation diagram), the basin of the equilibrium also enlarges. When comparing Figures 7(a), 8(b) and 9(b), this manifests itself by the fact that the bottleneck of the separatrix is less and less pronounced as $\gamma_t$ increases.
<figure><img src="image_9.png" /><figcaption>Figure 8: Boundary between the basins of attraction (green surfaces) of (11,12), represented in the $(\gamma_0, p_0, w_0)$-space, for saturating linear ramps of $\gamma$ with target value $\gamma_t = 1.1$. The slope is $s = 0.006$ (top left), $s = 0.01$ (top right), $s = 0.05$ (bottom left) and $s = 0.1$ (bottom right).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Boundary between the basins of attraction (dark blue surfaces) of (11,12), represented in the $(\gamma_0, p_0, w_0)$-space, for saturating linear ramps of $\gamma$ with target value $\gamma_t = 1.2$. The slope is $s = 0.006$ (top left), $s = 0.01$ (top right), $s = 0.05$ (bottom left) and $s = 0.1$ (bottom right).</figcaption></figure>
5. Conclusion

We investigated the influence of the time-variation of a control parameter on the dynamics of a simple model of reed musical instrument, written as a system of nonlinear ordinary differential equations. A single acoustic mode of the resonator has been taken into account, in order to keep the dimension of the dynamical system as low as possible. Despite this drastic simplification, a bifurcation analysis has shown that, when control parameters are constant in time, the system displays multistability in a range of the blowing pressure, which is one of the two main control parameters for the musician. We focused specifically on the question of the access to a particular regime in the context of multistability, that is to say to the question of the basins of attraction. In practice, it is unlikely that the musician can control directly the initial conditions on the different state variables independently. Because attack transients are known to be a crucial part of the musician expertise, we focus on their influence on the observed regime. Saturating ramps of the blowing pressure have been considered as prototypical profiles of attack transients, and the boundaries between the basins of attraction have been calculated through time domain simulations when such initial transients are accounted for. These are represented in the space of initial conditions on both the state variables and time-varying blowing pressure, for different target values and slopes of the blowing pressure.

Our results show that basins obtained in the cases without and with initial transients of the control parameter can differ significantly. In practice, this means that a given regime can become more difficult or easier to play when a well-chosen attack transient of $\gamma$ is performed compared to the case with no transient (or with an infinitely fast transient). In the phase space, this corresponds to a change of the size of the set of initial conditions leading to the considered regime. In the case with initial transients of parameters we showed that the basins strongly depends on the rate of change of the time-varying parameter. This influence is mediated by the interplay between the internal timescales of the system and the timescale on which the control parameter evolves. This strongly suggests that studying this model through a slow-fast dynamical systems point of view would be particularly valuable. However, this is beyond the scope of this paper, and is discussed elsewhere [27]. Additionally, investigating the influence of more realistic profiles of the control parameter might also be of interest for future work.

From a practical point of view, our results suggest that attack transients, and in particular the rate of change of control parameters, are an important aspect of the instrument control for the musicians, that allows them to access in a controlled manner to the desired sound regime among the different possible ones. This might be of practical importance for example for virtual prototyping of musical instrument. In particular, we argue here that basins of attraction might be strongly related to the *playability* or *ease of playing* of the different notes. These notions have been widely investigated in the literature and often related, for example, to oscillation thresholds [36, 48, 49]. Because the basin of attraction of a stable regime describes the set of initial conditions for which a regime is observed in practice, we expect this notion to be closely related to the probability of observing this regime in practice [11] and consequently to the fact that this regime is or not easy to play for the musician. In this context and as our results show, overlooking initial transients and investigating the stability of solutions and their basins while considering the control parameters as constant might not be sufficient to characterise the playability of regimes.

Our results provide a preliminary insight on the effect of time-varying parameters on the observability of regimes in a multistable model of musical instrument. Nevertheless, extending these investigations to a higher-dimensional model would be particularly valuable to unveil the impli-
cations of these observations in an actual instrument. Indeed, the model considered here takes
into account a single resonance mode of the instrument. Although this model has been shown
to reproduce the qualitative dynamics of a clarinet-like instrument as far as the first register is
concerned, recent studies, through extensive bifurcation analysis, have shown that higher-order
resonance modes induce a high degree of multistability, including between several periodic regimes
[50, 12, 36, 51]. The method considered here can easily handle higher-dimensional systems and
the same approach could therefore be considered. Nevertheless, visualizing the basins of attraction
when the dimension of the space of initial conditions exceeds three remains challenging [51]. More-
over, the attack of a note is also expected to induce variations in ζ, an effect that will be examined
in future work.

Finally, we expect that these results might be of interest beyond the specific case of musical
instruments considered in this article. In particular, the question of the influence of a time-varying
parameter relates to rate-induced dynamics which has attracted a lot of attention recently in a
variety of dynamical systems including biological and climate models [31, 30, 29]. The fact that
rate-induced tipping has such a strong impact on such a simple model illustrates the very general
character nature of this phenomenon which is nevertheless most often overlooked. Compared to
climate systems where tipping has been primarily investigated, acoustical systems and in particular
musical instruments are much easier to investigate and characterize experimentally. As such, these
numerical observations pave the way to future experimental works that would aim at characterising
and demonstrating rate-induced tipping in an actual physical system.
References

[1] G. Mueller and W. Lauterborn, "The bowed string as a nonlinear dynamical system," *Acta Acustica united with Acustica*, vol. 82, no. 4, pp. 657-664, 1996.

[2] R. Auvray, B. Fabre, and P.-Y. Lagrée, "Regime change and oscillation thresholds in recorder-like instruments," *The Journal of the Acoustical Society of America*, vol. 131, no. 2, pp. 1574-1585, 2012.

[3] T. Colinot, P. Guillemain, C. Vergez, J.-B. Doc, and P. Sanchez, "Multiple two-step oscillation regimes produced by the alto saxophone," *The Journal of the Acoustical Society of America*, vol. 147, no. 4, pp. 2406-2413, 2020.

[4] R. Mattéoli, J. Gilbert, S. Terrien, J.-P. Dalmont, C. Vergez, S. Maugeais, and E. Brasseur, "Diversity of ghost notes in tubas, euphoniums and saxhorns," *Acta Acustica*, vol. 6, p. 32, 2022.

[5] J. W. Coltman, "Jet offset, harmonic content, and warble in the flute," *The Journal of the Acoustical Society of America*, vol. 120, no. 4, pp. 2312-2319, 2006.

[6] J.-B. Doc, C. Vergez, and S. Missoum, "A minimal model of a single-reed instrument producing quasi-periodic sounds," *Acta Acustica united with Acustica*, vol. 100, no. 3, pp. 543-554, 2014.

[7] S. Terrien, C. Vergez, P. de La Cuadra, and B. Fabre, "Experimental analysis of non-periodic sound regimes in flute-like musical instruments," *The Journal of the Acoustical Society of America*, vol. 149, no. 3, pp. 2100-2108, 2021.

[8] E. Gourc, C. Vergez, P.-O. Mattei, and S. Missoum, "Nonlinear dynamics of the wolf tone production," *Journal of Sound and Vibration*, vol. 516, p. 116463, 2022.

[9] U. Feudel, A. N. Pisarchik, and K. Showalter, "Multistability and tipping: From mathematics and physics to climate and brain—minireview and preface to the focus issue," *Chaos: An Interdisciplinary Journal of Nonlinear Science*, vol. 28, no. 3, 2018.

[10] C. Grebogi, E. Ott, and J. A. Yorke, "Metamorphoses of basin boundaries in nonlinear dynamical systems," *Physical Review Letters*, vol. 56, no. 10, p. 1011, 1986.

[11] G. Datseris and A. Wagemakers, "Effortless estimation of basins of attraction," *Chaos: An Interdisciplinary Journal of Nonlinear Science*, vol. 32, no. 2, 2022.

[12] T. Colinot, C. Vergez, P. Guillemain, and J.-B. Doc, "Multistability of saxophone oscillation regimes and its influence on sound production," *Acta Acustica*, vol. 5, p. 33, 2021.

[13] S. W. McDonald, C. Grebogi, E. Ott, and J. A. Yorke, "Fractal basin boundaries," *Physica D: Nonlinear Phenomena*, vol. 17, no. 2, pp. 125-153, 1985.

[14] X. Xu, M. Wiercigroch, and M. Cartmell, "Rotating orbits of a parametrically-excited pendulum," *Chaos, Solitons & Fractals*, vol. 23, no. 5, pp. 1537-1548, 2005.

[15] S. Terrien, B. Krauskopf, and N. G. Broderick, "Bifurcation analysis of the yamada model for a pulsing semi-conductor laser with saturable absorber and delayed optical feedback," *SIAM Journal on Applied Dynamical Systems*, vol. 16, no. 2, pp. 771-801, 2017.

[16] J. P. Chávez, J. Schreiter, S. Siegmund, and C. Mayr, "A continuation approach for computing parameter-dependent separatrices in sram cells," *Applied Mathematical Modelling*, vol. 64, pp. 106-120, 2018.

[17] B. Krauskopf, H. M. Osinga, E. J. Doedel, M. E. Henderson, J. Guckenheimer, A. Vladimirsky, M. Dellnitz, and O. Junge, "A survey of methods for computing (un) stable manifolds of vector fields," *International Journal of Bifurcation and Chaos*, vol. 15, no. 03, pp. 763-791, 2005.

[18] J. Guckenheimer and C. Kuehn, "Computing slow manifolds of saddle type," *SIAM Journal on Applied Dynamical Systems*, vol. 8, no. 3, pp. 854-879, 2009.

[19] S. Farjami, V. Kirk, and H. M. Osinga, "Computing the stable manifold of a saddle slow manifold," *SIAM Journal on Applied Dynamical Systems*, vol. 17, no. 1, pp. 350-379, 2018.

[20] C. R. Hasan, B. Krauskopf, and H. M. Osinga, "Saddle slow manifolds and canard orbits in $r^4$ and application to the full hodgkin-huxley model," *The Journal of Mathematical Neuroscience*, vol. 8, pp. 1-48, 2018.

[21] S. Terrien, C. Vergez, and B. Fabre, "Flute-like musical instruments: A toy model investigated through numerical continuation," *Journal of sound and vibration*, vol. 332, no. 15, pp. 3833-3848, 2013.

[22] S. Terrien, R. Blandin, C. Vergez, and B. Fabre, "Regime change thresholds in recorder-like instruments: Influence of the mouth pressure dynamics," *Acta Acustica united with Acustica*, vol. 101, no. 2, pp. 300-316, 2015.

[23] R. Schumacher and J. Woodhouse, "The transient behaviour of models of bowed-string motion," *Chaos: An Interdisciplinary Journal of Nonlinear Science*, vol. 5, no. 3, pp. 509-523, 1995.

[24] A. Ernoult and B. Fabre, "Temporal characterization of experimental recorder attack transients," *The Journal of the Acoustical Society of America*, vol. 141, no. 1, pp. 383-394, 2017.
[25] B. Bergeot, A. Almeida, C. Vergez, and B. Gazengel, "Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure," *Nonlinear Dynamics*, vol. 73, pp. 521-534, 2013.

[26] B. Bergeot, A. Almeida, B. Gazengel, C. Vergez, and D. Ferrand, "Response of an artificially blown clarinet to different blowing pressure profiles," *The Journal of the Acoustical Society of America*, vol. 135, no. 1, pp. 479-490, 2014.

[27] B. Bergeot, S. Terrien, and C. Vergez, "Predicting transient dynamics in a model of reed musical instrument with slowly time-varying control parameter," *Chaos: An Interdisciplinary Journal of Nonlinear Science*, vol. 34, no. 7, 2024.

[28] J. K. Hale, "Introduction to dynamic bifurcation," in *Bifurcation Theory and Applications: Lectures given at the 2nd 1983 Session of the Centro Internationale Matematico Estivo (CIME) held at Montecatini, Italy, June 24-July 2, 1983*, pp. 106-151, Springer, 2006.

[29] P. D. Ritchie, H. Alkhayuon, P. M. Cox, and S. Wieczorek, "Rate-induced tipping in natural and human systems," *Earth System Dynamics*, vol. 14, no. 3, pp. 669-683, 2023.

[30] P. Ashwin, C. Perryman, and S. Wieczorek, "Parameter shifts for nonautonomous systems in low dimension: bifurcation-and rate-induced tipping," *Nonlinearity*, vol. 30, no. 6, p. 2185, 2017.

[31] A. Vanselow, S. Wieczorek, and U. Feudel, "When very slow is too fast-collapse of a predator-prey system," *Journal of theoretical biology*, vol. 479, pp. 64-72, 2019.

[32] M. E. McIntyre, R. T. Schumacher, and J. Woodhouse, "On the oscillations of musical instruments," *The Journal of the Acoustical Society of America*, vol. 74, no. 5, pp. 1325-1345, 1983.

[33] A. Chaigne and J. Kergomard, *Acoustics of musical instrument*, chapter 9. Springer, 2016.

[34] T. Colinot, L. Guillot, C. Vergez, P. Guillemain, J.-B. Doc, and B. Cochelin, "Influence of the ghost reed simplification on the bifurcation diagram of a saxophone model," *Acta Acustica united with Acustica*, vol. 105, no. 6, pp. 1291-1294, 2019.

[35] P.-A. Taillard, F. Silva, P. Guillemain, and J. Kergomard, "Modal analysis of the input impedance of wind instruments. application to the sound synthesis of a clarinet," *Applied acoustics*, vol. 141, pp. 271-280, 2018.

[36] J. Gilbert, S. Maugeais, and C. Vergez, "Minimal blowing pressure allowing periodic oscillations in a simplified reed musical instrument model: Bouasse-benade prescription assessed through numerical continuation," *Acta Acustica*, vol. 4, no. 6, p. 27, 2020.

[37] E. J. Doedel and B. Oldeman, "Auto-07p: continuation and bifurcation software," *Montreal, QC: Concordia University Canada*, 1998.

[38] S. H. Strogatz, *Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering*, chapter 6. CRC press, 2018.

[39] M. Pégeot, T. Colinot, J.-B. Doc, V. Fréour, and C. Vergez, "Playability of self-sustained musical instrument models: statistical approaches," *Acta Acustica*, vol. 8, p. 74, 2024.

[40] L. F. Shampine and M. W. Reichelt, "The matlab ode suite," *SIAM journal on scientific computing*, vol. 18, no. 1, pp. 1-22, 1997.

[41] E. Gwinn and R. Westervelt, "Fractal basin boundaries and intermittency in the driven damped pendulum," *Physical Review A*, vol. 33, no. 6, p. 4143, 1986.

[42] A. Basudhar and S. Missoum, "A sampling-based approach for probabilistic design with random fields," *Computer Methods in Applied Mechanics and Engineering*, vol. 198, no. 47-48, pp. 3647-3655, 2009.

[43] A. Basudhar and S. Missoum, "An improved adaptive sampling scheme for the construction of explicit boundaries," *Structural and Multidisciplinary Optimization*, vol. 42, pp. 517-529, 2010.

[44] N. Cristianini and B. Scholkopf, "Support vector machines and kernel methods: the new generation of learning machines," *AI Magazine*, vol. 23, no. 3, pp. 31-31, 2002.

[45] S. R. Gunn et al., "Support vector machines for classification and regression," *ISIS technical report*, vol. 14, no. 1, pp. 5-16, 1998.

[46] S. Missoum, C. Vergez, and J.-B. Doc, "Explicit mapping of acoustic regimes for wind instruments," *Journal of Sound and Vibration*, vol. 333, no. 20, pp. 5018-5029, 2014.

[47] W. Li, A. Almeida, J. Smith, and J. Wolfe, "How clarinettists articulate: The effect of blowing pressure and tonguing on initial and final transients," *The Journal of the Acoustical Society of America*, vol. 139, no. 2, pp. 825-838, 2016.

[48] R. Mattéoli, J. Gilbert, C. Vergez, J.-P. Dalmont, S. Maugeais, S. Terrien, and F. Ablitzer, "Minimal blowing pressure allowing periodic oscillations in a model of bass brass instruments," *Acta Acustica*, vol. 5, p. 57, 2021.

[49] E. Ukshini and J. J. Dirckx, "Objective characterization of saxophone mouthpiece playability," *Applied Acoustics*, vol. 206, p. 109309, 2023.

[50] J.-P. Dalmont, R. Mattéoli, S. Maugeais, S. Terrien, and C. Vergez, "The bichromator as the quintessence of a
saxophone and the question of harmonicity,” in *Forum Acusticum*, 2023.
[51] T. Passa, S. Terrien, S. Maugeais, and B. Gazengel, “Computing the basins of attraction of periodic solutions in a 4d model of a musical instrument,” in *European Nonlinear Dynamics Conference*, 2024.