© The Author(s), Published by EDP Sciences, 2024

1 Introduction

Musical instruments are complex dynamical systems. Some of them are self-sustained oscillators, meaning that a continuous energy supply can make them oscillate. Wind instruments and bowed string instruments belong to this category. These instruments admit a wide variety of regimes, which sometimes coexist for certain values of control parameters. This is known as multistability. Multistability is a common phenomenon in self-sustained oscillators. It has been observed experimentally, theoretically and numerically on a wide variety of instruments, e.g., single and double reed instruments [1–3], flutes [4], brass instruments [5], vocal folds [6] and bowed string instruments [7]. There are multiple works focusing on mapping the operating regimes of these instruments [8]. We find in [9, 10] among the first curves representing the evolution of the amplitude of a solution as a function of a control parameter. The stability of these solutions is also studied few years later [11]. A seminal work of Dalmont et al. [12] represents bifurcation diagrams of self-sustained musical instruments. These diagrams show the evolution of known solution properties as a function of one or more control parameters. The stability of these solutions is given, as well as some key features, such as amplitude and fundamental frequency for periodic solutions. Thanks to these diagrams, one can predict the behavior of an instrument whose control parameters vary in a quasi-static manner – i.e., slowly in relation to the system’s rate of evolution. In the case of multistability, bifurcation diagrams predict hysteretic behaviors, such as those of reed and brass instruments [13] or flute-like instruments [14]. Some unexpected regimes are also predicted with this method, such as the ghost note in brass instruments [15] or the wolf tone of the cello [16]. Nevertheless, it is much more difficult to predict the regime of a multistable instrument when the control parameters vary rapidly. In particular, the blowing pressure dynamics are shown to have an influence over the regime selection in flute-like instruments [17]. Similar results are shown in [18] for the saxophone.

Similarly, once the control parameters are fixed, the time required to reach the steady state depends on the previous evolution of the control parameters. This phenomenon was observed by [19, 20] on the clarinet. Generally speaking, the duration of attack transients [7, 21] or transitions between notes [22] are important topics, since the quality of a transient perceived by a musician or a listener depends greatly on its characteristics and duration [23]. The duration of the transients that gives rise to a sound is thus studied on multiple self-sustained instruments, from the trombone [24] to the cristal Baschet [25], using linear stability analysis. Like the prediction of the playing regime, the duration of the transient is estimated for control parameters varying in a quasi-static manner. Linear stability analysis cannot predict dynamic phenomena such as those observed by [19, 20]. Finally, despite their high efficiency in the quasi-static regime, current methods for analyzing self-oscillating instruments are limited when considering rapidly varying control parameters. By taking these rapid variations into account, the dimensionality of the problem becomes infinite since the parameters can follow any time series. Consequently, it motivates researchers to limit their studies to specific transient scenarios. For example, several authors consider pressure ramps of variable slope and study the influence of the pressure rise rate on the selected regime [26, 27] or the transient duration [20]. In this article, we study rapid variations in control parameters through stochastic initial conditions, in cases of multistability between an equilibrium and a periodic regime. We are interested in situations where the musician moves quickly from one quasi-static configuration to another, as is the case during an attack [21, 28], a break or a change of note [22, 29, 30]. Quasi-statistical analysis is not valid during this transition and the new quasi-static configuration is studied under unknown initial conditions. With this in mind, we study the statistical influence of initial conditions on the selected regime and on the duration of the transient. The probability of obtaining each regime is estimated with an already existing method called basin stability [31], and an original approach is proposed to determine the probability of obtaining each transient duration. Finally, the distributions of initial conditions considered for these statistical methods are discussed. In particular, initial conditions generated by specific transient control scenarios are studied.

The remainder of the document is organized as follows. The system chosen to introduce the methods is given in Section 2. Section 3.1 focuses on the basin stability metric and its use for musical instruments. The transient duration is investigated in Section 3.2 and the influence of the distribution of initial conditions is explored in Section 4. To highlight the interest of the methods on more sophisticated and high dimensional systems, an application to a physical model of trumpet is presented in Section 5. The different methods and results are discussed in Section 6. Perspectives and conclusions are drawn in Section 7.

2 Minimal multistable system

Most physical models of self-sustained musical instruments have no exact analytical solutions for transient or steady-state regimes. Therefore, numerical methods are often used to study them, such as finite differences [32] or harmonic balance method [33]. A number of authors have proposed less detailed models with exact analytical solutions and only few state variables. These simplified models help to investigate the mechanisms behind self-oscillations [34]. In order to compare our statistical approach with exact analytical results, we have chosen to illustrate our methods on such a system. In this article, we focus on a Van der Pol oscillator with fifth-order nonlinearity, as described in [35]. This oscillator is characterized by a region of multistability, where both an equilibrium and a periodic solution are stable. This feature can also be observed on several musical instruments, such as saxophones [27] or brass instruments [5]. Its phase space is of dimension two and can therefore be displayed simply, enabling more detailed dynamic analysis. A technical publication [36] presents the musical interest of this system and provides a demonstrator solving it in real time (https://zenodo.org/records/8413627). Note that even though the global behavior of musical instruments is varied and rich, the system studied here can be obtained by change of variables and polynomial expansions from any system presenting this type of multistable behavior. It is only used to illustrate the interest of a statistical approach to predict the behavior of a multistable system under transient excitations.

The motion of the considered oscillator is given by the following equation:

$$\ddot{x}+\left[-\mu +\sigma (x^2+{\dot{x}}^2)+\nu (x^2+{\dot{x}}^2{)}^2\right]\dot{x}+x=0.$$(1)

We impose σ = −1.5 and ν = 0.1 in order to obtain a multistable saxophone-like or brass-like behavior. The parameter μ is the control parameter of this simple “musical instrument” and its effect is analogous to the blowing pressure for wind instruments [36]. A harmonic balance method with only one harmonic gives the amplitude X of solutions of the form:

$$x=X\cos (t+\phi ).$$(2)

There are up to three solutions, depending on the value of μ:

$$\{\begin{array}{cc}{X}_{\mathrm{Eq}}=0,& \\ X_{2\pi +}=\sqrt{\frac{-\sigma +\sqrt{{\sigma }^2+4\mu \nu }}{2\nu }}& \mathrm{if} \mu \ge {\mu }_{\mathrm{SN}},\\ X_{2\pi -}=\sqrt{\frac{-\sigma -\sqrt{{\sigma }^2+4\mu \nu }}{2\nu }}& \mathrm{if} {\mu }_{\mathrm{SN}}\ge \mu \ge {\mu }_H,\end{array}$$(3)

with μ_H = 0 and ${\mu }_{\mathrm{SN}}=-\frac{{\sigma }^2}{4\nu }$. These solutions consist of an equilibrium (X_Eq) and two 2π-periodic solutions (X_2π+ and X_2π−), which are represented on a bifurcation diagram in Figure 1. This oscillator shows an inverse Hopf bifurcation at μ = μ_H and a saddle-node bifurcation at μ = μ_SN. Between these two bifurcations, the system is multistable: both the equilibrium X_Eq and the largest periodic solution X_2π+ are stable. Saxophones [27] and brass instruments [5] are likely to present this bifurcation sequence and therefore this multistability behavior.

Figure 1 Bifurcation diagram of a fifth order Van der Pol oscillator given in equation (1) with respect to the control parameter μ, with σ = −1.5 and ν = 0.1. Stable solutions are represented with continuous lines, unstable solutions in broken lines and bifurcations are indicated with star markers. The region of multistability is shaded in grey.

A traditional bifurcation diagram only depicts the steady state solutions. The influence of initial conditions on the transient and the steady-state behavior is concealed. This influence can be seen on a phase plane for a given value of the control parameters. In Figure 2, we show the phase plane of the fifth order Van der Pol oscillator for μ = −3, which is a situation of multistability. For the remainder of the paper, we use this value each time we need to choose a specific μ for illustration. The equation of motion (1) is numerically integrated with the function ode45 of Matlab [37], for three initial states very close to one another. In that case, the slight change in initial conditions leads the system either to the equilibrium (in blue) or to the stable limit cycle (in red). The two oscillating trajectories have different transient durations: one is almost instantaneous while the other lasts approximately one period.

Figure 2 Illustration of the sensitivity to the initial conditions, for the Van der Pol oscillator, regarding the steady state regime and the transient duration. On the left side, the phase plane of the system for μ = −3, with trajectories coming from three close initial states (indicated with star markers). On the right side, the evolution of the position x with respect to time (dimensionless).

In short, Figures 1 and 2 illustrate that the steady state and the transient duration of this minimal self-sustained musical instrument are sensitive to the initial conditions and to the control parameter value around a bifurcation. This sensitivity is merely a mark of non linearity and is shared by most non linear dynamical systems. In the Sections 3.1 and 3.2, we present two statistical approaches to predict the selected regime and the transient duration for such a system.

3 Methods

3.1 Basin stability

Some self-sustained musical instruments are multistable. In such a case, once their control parameters are settled, their steady state depends exclusively on their initial state. The subset of initial conditions leading to a specific attractor is called its basin of attraction [38, 39]. In Figure 3, we depict the basins of attraction of the fifth order Van der Pol oscillator, for μ = −3. The basin of the equilibrium is shaded in blue and the basin of the periodic solution is shaded in red. In that case, the two basins are separated by the unstable limit cycle of amplitude X_2π−, plotted in broken line and given in equations (1) and (3). In the general case, however, the geometry of a basin of attraction can be more complex and its boundaries difficult to find. We discuss this problem in Section 6. Moreover, in the scope of this article, we consider that the initial conditions resulting from a quick change in the control parameters are unknown. They are randomly taken within a subset of the initial state space. Since the initial state is unknown, the precise geometry of the basins of attraction does not suffice to predict its playing regime.

Figure 3 Basins of attraction of the fifth order Van der Pol oscillator (μ = −3). Initial conditions taken in the red shaded area give rise to the periodic solution whereas initial conditions taken in the blue shaded area give rise to the equilibrium solution.

To take into account the uncertainty of the initial state, we propose to use the notion of basin stability, a probabilistic metric introduced by [31]. The principle is to evaluate the size of the basins of attraction with respect to the total size of the considered space of initial states. The result gives the probability to reach each attractor, considering a random initial state. Formally, considering a subset of interest of the space of initial conditions $\mathcal{Q}$, an attractor A, its basin of attraction B and a probability density function ρ, the basin stability ${\mathcal{S}}_{\mathcal{B}}\left(A\right)$ is defined as follows:

$${\mathcal{S}}_{\mathcal{B}}\left(A\right)={\int }_{\mathcal{Q}}\mathrm{}{\mathbf{I}}_B\left(\mathbf{X}\right)\rho \left(\mathbf{X}\right)\mathrm{d}\mathbf{X},$$(4)

where

$${\mathbf{I}}_B\left(\mathbf{X}\right)=\{\begin{array}{ll}1,& \mathrm{if} \mathbf{X}\in B,\\ 0,& \mathrm{otherwise}.\end{array}$$(5)

The probability density function ρ can follow a specific distribution if the random choice of initial conditions is not uniform. This aspect is discussed in Sections 4 and 6. Due to the nature of ρ as a probability density function over $\mathcal{Q}$, the basin stability metric is always comprised between 0 and 1. ${\mathcal{S}}_{\mathcal{B}}\left(A\right)=0$ indicates that the basin of attraction of the solution is statistically never reached. ${\mathcal{S}}_{\mathcal{B}}\left(A\right)=1$ means that the basin of attraction of A occupies all of the region of interest $\mathcal{Q}$, and that A is reached by every trajectory. All the illustrations of Section 3 are realized with a uniform distribution. The volumic integral in the expression of ${\mathcal{S}}_{\mathcal{B}}$, in equation (4), is rarely calculable exactly and we rather compute the corresponding discrete sum. In other words, we take N random samples in $\mathcal{Q}$ and we evaluate the proportion of samples belonging to the studied basin of attraction. The basin stability estimation, which is noted ${\widehat{\mathcal{S}}}_{\mathcal{B}}$, then writes:

$${\widehat{\mathcal{S}}}_{\mathcal{B}}\left(A\right)=M/N,$$(6)

where M is the number of samples belonging to B, the basin of attraction of A. As pointed out in [31, 40, 41], this computation corresponds to N independent trials with probability of success ${\mathcal{S}}_{\mathcal{B}}$. The resulting standard error due to sub sampling writes:

$$\mathrm{err}\left({\widehat{\mathcal{S}}}_{\mathcal{B}}\right)=\sqrt{\frac{{\mathcal{S}}_{\mathcal{B}}(1-{\mathcal{S}}_{\mathcal{B}})}{N}}.$$(7)

In practice, ${\mathcal{S}}_{\mathcal{B}}$ is unknown so it is replaced by its estimation ${\widehat{\mathcal{S}}}_{\mathcal{B}}$ in equation (7). This error only accounts for sub-sampling. It is interesting to notice that the error does not depend on the dimension of the state space. As a result, the basin stability is a metric suitable for low dimension problems as well as for high dimension ones. In general, for more complex systems, the basin boundaries are unknown and the classification method would consist in time integrating the system and observing toward which attractor it would converge. Such strategy is illustrated in Figure 4. However, for this minimal multistable system, the basin boundary is analytically known: the boundary is the unstable periodic solution whose ${\mathcal{L}}_2$-norm X_2π− is given in equation (3). As a result, it is much faster to classify the samples with this analytical expression rather than with time integrations.

Figure 4 Illustration of the classification method in the general case. The plotted trajectories are obtained with time integrations starting from the stars as initial conditions. Red trajectories are inside the basin of attraction of the stable periodic solution X2π+ whereas the blue ones belong to the basin of the equilibrium XEq.

Figure 5 gives the basin stability of the Van der Pol oscillator. To compute it, we arbitrarily choose a subset of the phase space $\mathcal{Q}$ which includes all the attractors and which is independent of μ. For μ ∈ [−7, 1] as in Figure 1, the subset $\mathcal{Q}=\left(\right[-4,\mathrm{ }4]\mathrm{ }\times \mathrm{ }[-4,\mathrm{ }4\left]\right)$ meets these conditions. To uniformly sample $\mathcal{Q}$, we use a Latin Hypercube Sampling of 100 samples. The basin stability of the equilibrium is plotted in blue and the one corresponding to the periodic solution in red. The error bars correspond to the standard error given in equation (7). In addition, within the multistability region, i.e., for μ_SN ≥ μ ≥ μ_H, the exact basin stabilities write:

$$\{\begin{array}{l}{S}_{\mathcal{B}}\left(X_{\mathrm{Eq}}\right)=\frac{\pi X_{2\pi -}^2}{\mathrm{area}\left(\mathcal{Q}\right)},\\ S_{\mathcal{B}}\left(X_{2\pi +}\right)=1-S_{\mathcal{B}}\left(X_{\mathrm{Eq}}\right).\end{array}$$(8)

Figure 5 Basin stability of the two stable solutions computed from 100 samples. In red: basin stability of the periodic solution. In blue: basin stability of the equilibrium. Error bars correspond to the standard error of the size estimation of the Monte-Carlo method due to sub-sampling equation (7). In black: basin stability computed with the analytic expression of the basin boundary (circle of known diameter X2π− given in (3)). The region of multistability is shaded in grey.

These exact solutions are represented with black lines. They mainly remain inside the error bars estimated with ${\widehat{\mathcal{S}}}_{\mathcal{B}}$ (i.e., ${\mathcal{S}}_{\mathcal{B}}$ replaced by ${\widehat{\mathcal{S}}}_{\mathcal{B}}$ in Eq. (7)). This validates the error estimation of the sample-based method given in equation (7). Outside the multistability area, the only stable solution has a basin stability equal to 1 whereas the other one has a basin stability equal to 0. Regarding this metric, non existing solutions and unstable solutions are identical and have a zero basin stability. Inside the multistability area, both solutions have a non zero basin stability. The stability of the periodic solution increases with μ whereas the stability of the equilibrium decreases. One can observe that ${\mathcal{S}}_{\mathcal{B}}$ is continuous at the Hopf bifurcation but discontinuous at the saddle-node bifurcation. Indeed, in the first case, the basin of attraction of the equilibrium shrinks progressively until it becomes a dot in the phase space for μ = μ_H. In the second case, the periodic solution suddenly disappears at μ_SN, even though its basin of attraction had previously occupied a large part of $\mathcal{Q}$. The value of ${\mathcal{S}}_{\mathcal{B}}$ corresponds to the probability to obtain the related solution, given a random initial state inside $\mathcal{Q}$. If a musician were to impose uniformly random initial conditions in $\mathcal{Q}=\left(\right[-\mathrm{4,4}]\times [-\mathrm{4,4}\left]\right)$, Figure 5 indicates that he/she would have at least a 60% chance of reaching the periodic regime in the steady state. However, the assumption of a musician peeking a random initial condition uniformly seems oversimplified. This topic is discussed in Sections 4 and 6.

Overall, the basin stability indicates which values of the parameters (here, the only parameter is μ) fosters one solution or another. It gives the probability to produce a specific regime, considering a subset $\mathcal{Q}$ of the phase space and a probability density function ρ of $\mathcal{Q}$. The result depends greatly on the choice of $\mathcal{Q}$ and ρ. The standard error of this metric depends on the number of samples and on the relative size of the studied basin of attraction. It does not depend on the phase space dimension however and thus is very well adapted to high but finite dimensional systems.

3.2 Transient duration

In music, transient phenomena play a crucial role in sound and instrument discrimination [42–44]. Musicians may also have expectations concerning the characteristics of these transients, and in particular concerning their duration. For example, Guettler and Askenfelt [23] highlight the importance of the type and duration of violin attacks on their quality as perceived by professional string players. Several authors [21, 22, 25, 45] assume that musicians and instrument makers look for short transitions, whether for the onset or for sequences of linked notes. In this paper, we thus consider the reaction time of an instrument – i.e., its tendency to produce short transients – as a playability criterion. In self-sustained instruments, the transient depends on the parameters of the system but also on the initial conditions, as can be seen on Figure 6. In this figure, the evolution of the fifth order Van der Pol oscillator has been computed for a large number of random initial states, all other things being equal. In Figure 6a, the trajectories are represented in the phase plane whereas in Figure 6b, the same trajectories are represented with respect to time. The transient part of a trajectory is defined between its initial state and the moment it reaches the neighborhood of an attractor. The transient duration is denoted τ hereafter. The neighborhoods of the two attractors are shaded in grey and are defined at a distance ε in ${\mathcal{L}}_2$-norm from the attractor. Although the size of these neighborhoods has an impact on the transient duration results, for now it is arbitrarily set to ε = 0.1. Figure 6b highlights the diversity of transient durations that can be obtained only by picking different initial conditions. Some tendencies seem to appear. For instance, for this distribution of initial conditions, trajectories leading to the equilibrium are usually longer than those leading to the periodic solution.

Figure 6 (a) Phase trajectories and (b) time evolution of the associated norm of the system’s state for 300 random initial conditions. Blue trajectories lead to the equilibrium and red trajectories lead to the periodic solution. The neighborhoods of the two stable solutions are shaded in grey. They are defined at a distance ε = 0.1 in ${\mathcal{L}}_2$-norm from the stable solutions.

These tendencies stand out when the distributions of the transient durations are represented graphically. Information on a distribution can be displayed as in Figure 7. The upper panel gives the cumulative distribution function, i.e., the proportion of trajectories with a transient duration shorter than the value τ_i represented in abscissa. The two lower panels depict the corresponding violin plots [46, 47], where the box plots give the extreme values (tips of the whiskers), the median and the two quartiles, and the shaded curves represent the probability density function of the distributions. To compute the latter, we used the kernel density estimator of Matlab (ksdensity). This statistic information is calculated for each attractor separately. As in Section 3.1, we use a Latin Hypercube Sampling to uniformly explore the space of initial states. With this sampling, we do not control the exact number of samples in each basin of attraction, but it can be estimated with the basin stability ${\mathcal{S}}_{\mathcal{B}}$ previously computed. As a result, for lowest values of ${\mathcal{S}}_{\mathcal{B}}$, the statistical transient analysis relies on only few samples and may not be very representative. It should be remembered that such misrepresentations occur only on highly improbable regimes. If a specific case calls for in-depth study of these regimes, the region of interest $\mathcal{Q}$ can be adapted, leaving the method otherwise unchanged. The cumulative distribution function lends itself well to interpretation in musical context, where one could define a longest acceptable transient duration. The probability of the transient being shorter than this upper limit can be read directly on the cumulative distribution function, under the hypothesis of random initial conditions following a uniform law.

Figure 7 (Top) Cumulative probability function of the transient duration distribution. In ordinate, the proportion of trajectories with a transient duration below the value given in abscissa. (Bottom) Corresponding violin plots. The box plots indicate the extreme values, the median, the upper and lower quartiles whereas the curves around represent the probability density function of the transient duration distribution. A latin hypercube of 1000 samples is used and the transient durations are defined with ε = 0.1.

We use the violin plot representation in Figure 8 to depict the evolution of the transient duration distribution with respect to the control parameter μ. Plotting the probability density function over the boxplot is particularly interesting when there are multiple maxima in the distribution since multiple maxima are not visible on a box plot. For instance, some examples of bi-modal distributions can be observed around the saddle-node bifurcation (for μ ≈ μ_SN).

Figure 8 Violin plots of the transient duration with respect to the control parameter μ. In blue: equilibrium. In red: periodic solution. Number of samples: N = 1000. Panel (b) is a magnified version of the transients leading to the periodic solution. The region of multistability is shaded in grey.

Figure 8 shows that most transients leading to the periodic solution are very short compared with its period (T = 2π). Indeed, the violin plots are centered around low values of τ (a magnified view for μ = −3 is visible in Fig. 8b). However, the maximum transient durations are markedly larger than the median duration, up to 800 times longer for μ = −0.26. This indicates that, for most values of μ, a few initial conditions have extremely long transients. In fact, these particular initial states are very close to the unstable limit cycle. In that case, the trajectory follows the unstable limit cycle for some time while slowly diverging from it, which lengthens the transient duration. These rare and very long durations have particularly large values in the vicinity of the Hopf bifurcation. For the equilibrium, transient durations are longer and distributed more uniformly for all the values of μ. Outside the vicinity of the two bifurcations, the transient durations leading to both solutions statistically decrease as μ increases. This is coherent with the fact that −μ is the linear part of the damping. As μ increases, the system becomes less and less damped until μ eventually becomes positive and energy is added to the system. Near the saddle-node bifurcation, this trend does not apply. The equilibrium becomes the only stable solution but its transient duration globally increases around μ_SN. This is a usual phenomenon around a saddle-node bifurcation. It is interpreted by Strogatz in [38] as the remaining ghost of a nearby attractor that does not exist anymore. Around the Hopf bifurcation, the periodic solution becomes the only stable solution but its maximum transient duration significantly increases. However, apart from these very few long transients, Figure 8b shows that the overall trend does not change significantly for the periodic solution around the Hopf bifurcation. Considering initial conditions chosen from a uniform distribution, these long transients are highly improbable.

4 Influence of the distribution of initial conditions

In order to provide a more meaningful analysis in a context of musical performance, we consider another distribution for the initial conditions, informed by temporal evolutions of the control parameter. More precisely, we now consider initial conditions resulting from a transient variation of the control parameter μ, just before it reaches a constant target value. Our goal is therefore to consider what we call hereafter a “transient-informed” distribution and to analyze how it affects the basin stability results. To obtain a transient-informed distribution, we compute the time evolution of the system under several transient variations of μ. These transient-informed distributions could also be used to compute the transient duration statistics, but we only focus on the basin stability here, since the results are more remarkable.

4.1 Control scenarios

Systems such as the fifth order Van der Pol oscillator presenting an inverse Hopf bifurcation followed by a saddle-node bifurcation produce a sound of non-zero amplitude when the control parameter exceeds the Hopf bifurcation (this is the case for the trumpet and the saxophone for certain sets of parameters). Nevertheless, it is possible to reduce the sound amplitude afterward by decreasing the control parameter under the Hopf bifurcation value, into the multistability region. To produce a low-volume sound at the onset, one can attempt to realize this gesture quickly by reducing the control parameter, and move the system into the multistable region before the limit cycle is attained. However, this strategy may fail if the system remains too close to the equilibrium during the transient. It would become trapped in the equilibrium’s basin of attraction when the parameter stops varying, which would result in no stable sound being produced. We consider piecewise linear control scenarios, as depicted in Figure 9. These scenarios are defined by four parameters: their initial value μ₀, their maximal value μ_max that corresponds to what can be called the overshoot value, their final value μ_end at which the basin stability is computed and their duration τ_att that corresponds to the “attack duration”. For simplicity, we choose the same slope (in absolute value) for the increasing and the decreasing phases (before and after the overshoot).

Figure 9 Control scenarios used to compute the transient-informed initial conditions. The multistability region is shaded in gray. μ0 is the initial value of μ, μmax is the overshoot, μend is the final value of μ, at which the basin stability is computed, μH and μSN are the value of μ at the Hopf and at the saddle-node bifurcations and τatt is the attack duration.

The order of magnitude of μ₀ is based on the bifurcation diagrams of the trumpet models presented in [5]. The objective is to give as much importance to the multistability as in these trumpet models. Therefore μ₀ is set to have a similar ratio between the size of the multistability region [μ_SN, μ_H] and the size of the stability region of the equilibrium [μ₀, μ_H] as in the models of [5]. μ₀ is then set to μ₀ = −25. At the end of the attack, μ = μ_end and the basin stability is computed for this value of μ. The overshoot μ_max and the duration of the attack τ_att are the two remaining parameters of the control scenarios. In the following, we study the influence of a given overshoot μ_max on the regime selection, while the attack duration τ_att remains free to take different values. We arbitrarily choose μ_max ∈ [1, 10] and τ_att ∈ [1, 20] (dimensionless values). In [21], Ernoult and Fabre measured on three recorder players pressure rise times going approximately from 5 to 150 ms, for all fingerings between C5 (1046.5 Hz) and C6 (2093 Hz). Therefore, by choosing τ_att ∈ [1, 20], we keep a ratio between the longest and the shortest attacks with the same order of magnitude as in [21].

4.2 Distribution generation

To generate a transient-informed distribution of initial conditions for a specific μ_end, we apply 10 control scenarios with different durations (τ_att ∈ [1, 20]) and with a unique value of μ_max, to N₀ = 10 pre-initial conditions chosen close to the equilibrium (Latin Hypercube Sampling taken inside x ∈ [−0.1, 0.1] and $\dot{x}\in [-\mathrm{0.1,0.1}]$). As a result, we obtain a transient-informed distribution of N = 100 initial conditions depicted in Figure 10. Firstly, it can be observed that this distribution is far from being uniform; most transient-informed initial conditions are very close to either the equilibrium or the periodic solution. It is expected for such a system to leave an unstable equilibrium through its fast eigendirection (i.e., associated with the eigenvalue of highest modulus) and to converge towards a stable equilibrium through its slow eigendirection (i.e., associated with the eigenvalue of smallest modulus). Such behavior is explained by Strogatz in example 5.2.3, p.133 of [38]. Precisely, around the equilibrium in Figure 10, the states of the system at μ = μ_end are gathered along a specific direction whereas the states of the system at μ = μ_max are regrouping along another direction. In addition, we plot on the same figure the fast eigendirection of the equilibrium at μ = μ_max and its slow eigendirection at μ = μ_end. At the overshoot, the system clearly leaves the equilibrium through its fast eigendirection. At the end of the attack, some of the trajectories converge toward the equilibrium. If the simulation was performed over a longer duration, they would eventually gather along the slow eigendirection of the equilibrium.

Figure 10 Example of a transient-informed distribution generated with 10 control scenarios of different durations. (•) Pre-initial Latin Hypercube of size N0 = 10; () states of the system when μ = μmax; () transient-informed distribution at the end of the attack. The black dotted line is the fast eigendirection of the equilibrium at μ = μmax and the red one is its slow eigendirection at μ = μend.

4.3 Basin stability with transient-informed distribution

The basin stability is then computed using these transient-informed distributions. The results are shown in Figure 11 for three different overshootvalues. For the sake of readability, only the basin stability of the periodic solution is represented (we recall that the basin stability of the equilibrium is the complementary to one).

Figure 11 Basin stability of the periodic solution computed with different distributions. The red curves are computed with transient-informed distributions of initial conditions, each curve corresponding to a specific overshoot μmax. The black curve is computed with a uniform distribution (it is the same as in Fig. 5 without the error bars).

The result computed with the uniform distribution is represented in black. Firstly, this Figure 11 highlights that the basin stability depends greatly on the considered distribution, even though some features are preserved: ${\mathcal{S}}_{\mathcal{B}}\left(\mu \right)$ is discontinuous at μ = μ_SN, continuous at μ = μ_H and monotonously increasing for the cases where μ_max = 3 and 4. In the case where μ_max = 2, the overshoot scenario does not suffice to leave the basin of the equilibrium and no oscillations are observed in the multistable region. Secondly, for the situation considered, this Figure 11 shows that the large-overshoot attacks are more likely to make the system end up on the periodic solution. Moreover, the basin stability increases faster when μ_end gets closer to μ_H. In terms of playability, a slowly varying ${\mathcal{S}}_{\mathcal{B}}$ can be interpreted as a region where the difficulty to produce a sound does not depend much on the final value of the control parameter. For the situation under study, this means that for a given value of μ_max, all the attacks with μ_end in this region of slow varying ${\mathcal{S}}_{\mathcal{B}}$ have almost the same probability to converge towards the periodic solution. As a result, if we ignore the difficulty for a musician to maintain μ > μ_SN, there is no additional difficulty to reducing the sound amplitude to its minimal value.

5 Application to a trumpet model

In order to evaluate the interest of the methods introduced in previous sections to a more advanced and high dimensional model of musical instrument, the analysis of basin stability and transient durations is applied in this section to a physical model of B_♭ trumpet.

5.1 Presentation of the model

The model used is identical to the one considered in [5]. It is detailed in Appendix A. The lips are represented by a one-degree-of-freedom oscillator. The air column inside the body of the instrument, also known as the resonator, is modeled by a modal truncation of its input impedance. These two elements are coupled by a nonlinear function that models the air jet through the lip channel. In this article, we consider an 11-mode truncation of the resonator and fix the resonance frequency of the lips around the second impedance peak (f_L = 200 Hz and f₂ = 232.7 Hz). This peak corresponds to the lowest note of the instrument, in “normal” playing conditions (the first peak being associated with the “pedal note”). The variables of this model are the position and velocity of the lip as well as the real and imaginary parts of the 11 modal pressures. This results in a system of dimensionality 24. The parameter values used in this article are given in Appendix A.

5.2 Bifurcation diagram

Unlike the fifth-order Van der Pol oscillator, this trumpet model has no analytical solutions. To determine its solutions, numerical continuation can be used. Here, we use the Manlab software [33]. The bifurcation diagram of the model is shown in Figure 12. The amplitude of the mouthpiece pressure p and the frequency of the sound produced f_play are plotted with respect to the blowing pressure p₀. In this zone, the system exhibits an inverse Hopf bifurcation followed by a saddle-node bifurcation, like the Van der Pol oscillator presented in previous sections. The periodic regime is associated with the second mode of the resonator, represented on the bottom panel of Figure 12 by a horizontal dotted line.

Figure 12 Bifurcation diagram of the 11-mode trumpet model, calculated with Manlab. Top pane: amplitude of the oscillations ${\Vert p\Vert }_{{\mathcal{L}}_2}$ versus blowing pressure p0. Bottom pane: fundamental frequency of the oscillations fplay versus blowing pressure p0. Stable solutions are shown as solid lines, unstable solutions as dotted lines. Bifurcations are represented by stars.

5.3 Basin stability and transient durations computed with bSTAB

To calculate the system’s basin stability and transient durations, we use the bSTAB toolbox [48]. For our application, we have made a number of modifications to this toolbox, which are detailed in Appendix B. These modifications rely on convergence criteria in order to identify when the system has reached steady state. The numerical integration lasts until these criteria are met, rather than for an arbitrary time which would be identical for all trajectories. These modifications allow us both to reduce the duration of the time integrations required to compute basin stability, and to calculate the duration of transients (which is not provided in the original toolbox). These criteria must be defined carefully, in order to avoid samples missclassification, leading to errors in the estimation of basin stability and transient durations.

5.3.1 Initial conditions from a uniform distribution

As for the Van der Pol oscillator, we first consider a uniform random distribution of initial conditions ${\mathcal{Q}}_u$ selected inside a hyper-rectangle. This hyper-rectangle is chosen large enough to contain all the solutions, whatever the value of p₀. Its bounds are set as the extrema of the stable periodic solution at p₀ = 0.8 kPa, rounded up to the nearest integer value.

The basin stability calculated with ${\mathcal{Q}}_u$ for N = 300 samples is given in Figure 13a, the transient duration distributions are given in Figure 13b and a zoomed-in view is given in Figure 13c. The basin stability shares some features with the Van der Pol oscillator. Indeed, it is discontinuous at the saddle-node bifurcation and remains continuous at the Hopf bifurcation. The observations and comments made about this feature on the simple Van der Pol oscillator translate directly to this more complicated system. In addition, ${\mathcal{S}}_{\mathcal{B}}$ is monotonous. However, unlike the Van der Pol oscillator, for dynamical systems of dimensionality three or more, the basins’ boundaries are different than the unstable periodic solutions. These objects even have different dimensionalities (the basins’ boundaries are hypersurfaces and the limit cycles are hyperlines). Consequently, the basin stability of the equilibrium cannot be inferred from a size estimate of the unstable periodic solution and the sampling approach is then required.

Figure 13 Basin stability (a) and transient duration distributions (b) of the trumpet model calculated with a uniform ${\mathcal{Q}}_u$ distribution of N = 300 initial conditions. The dotted lines in (b) indicate the limit of the scale used for the zoomed view in (c). Error bars are given by the equation (7).

Concerning transient durations, the Hopf bifurcation induces a few long transients in its vicinity, leading both to the equilibrium and to the periodic solution. This is probably related to the increasing timescale of the slowest eigendirection whose stability is reversed at the Hopf bifurcation. Eventually, the timescale of this eigendirection goes to the infinity at the Hopf bifurcation. The ghost effect around the saddle-node bifurcation, where the periodic solutions do not exist yet, is also visible but less than for the Van der Pol oscillator. This may be due to the fact that in a phase space of dimensionality two, for a value of the control parameter slightly below the value of the saddle-node, the trajectories initiated far away from the equilibrium necessarily have to enter the region of influence of this ghost limit cycle. On the contrary, in a phase space of higher dimensionality, these trajectories can pass away from this region. Moreover, the median and the interquartile range of both solutions increase around the Hopf bifurcation (even though it is hardly noticeable for the periodic solution with that scale).

Finally, this trumpet model sometimes exhibits transients of more than several seconds, which is extremely long compared to usual musical timescales. Beyond a certain duration, these asymptotic regimes can be considered unplayable or at least difficult to play. However, the regime that is heard during one of these extremely long transients can be treated as a playable regime by the musician, even if it is unstable. In this respect, one could see a percussion instrument as a dynamic system with only one asymptotic solution, the equilibrium, but whose transients are long enough to be used to make music. We could then consider a “transient” category for samples with a transient duration exceeding a certain “musical duration”, rather than classifying them according to their asymptotic solution.

5.3.2 Initial conditions from a distribution based on archetypal time variations of p₀

In the following, we study the influence of attack transients on basin stability. We apply the same method as in Section 4 to generate distributions ${\mathcal{Q}}_{\mathrm{tb}}$ based on transient evolutions of p₀. Again, we consider pressure rises with overshoot (see Fig. 9), in the manner of a strong tonguing attack or a sforzando, as can be found in [49]. The configuration and the parameters chosen here differ from [49]. Consequently, we consider slightly different control transients than those measured in that article, so that p_0,max ∈ [1, 10] kPa and τ_att ∈ [10, 500] ms. Ten attack transients are applied to each sample of a uniform distribution of N₀ = 100 pre-initial conditions. This leads to a final distribution of N = 1000 initial conditions. Each variable of the pre-initial distribution is selected in ranges 0.3 times smaller than in Section 5.3.1. These ranges are chosen arbitrarily and have a significant influence on the results.

In Figure 14a, we present the influence of the attack duration τ_att on the behavior of the system. Each curve has a fixed value of τ_att and 10 different values of p_0,max ∈ [1, 10] kPa. It is striking to see that the basin stability does not evolve monotonously with τ_att. Indeed, one would expect that the longer the system stays in the region of monostability of the periodic solution, the higher the chances it has to end up in its basin of attraction. However, when τ_att increases from 10 to 50 ms, ${\mathcal{S}}_{\mathcal{B}}$ decreases. ${\mathcal{S}}_{\mathcal{B}}$ does not evolve much for τ_att ∈ [50, 100] ms, and it then increases, as expected, for τ_att > 100 ms. These observations can be analyzed by considering recent results concerning dynamical bifurcations [50]: how close the system has come to the equilibrium when the Hopf bifurcation is crossed is crucial to predict the dynamics above the bifurcation. This closeness is the result of the time spent by the system below the bifurcation but also of the number of significant digits used in the simulation (numerical noise is therefore inevitable). This makes direct interpretation not straightforward. In Figure 14b, we study instead the influence of the overshoot amplitude p_0,max for various attack durations τ_att ∈ [10, 500] ms. In that case, ${\mathcal{S}}_{\mathcal{B}}$ increases with the overshoot value, until it reaches its maximal value for p_0,max ∈ [5, 7] kPa. Hence, it seems that there are optimal values of transient control parameters to maximize ${\mathcal{S}}_{\mathcal{B}}$. Further study would be needed to understand the dynamical reasons leading to this non trivial behavior.

Figure 14 Basin stability of the trumpet model calculated with distributions ${\mathcal{Q}}_{\mathrm{tb}}$ of N = 1000 initial conditions based on attack transients (10 scenarios applied to each of the N0 = 100 pre-initial conditions). (a) Each curve corresponds to a specific attack duration τatt. (b) Each curve corresponds to a specific overshoot value p0,max. Only the basin stability of the periodic regime is shown.

6 Discussion

6.1 Computational cost

The methods presented in this paper allow to compute the basin stability and the transient duration distributions of a system. These metrics describe the global, statistical behavior of a system and could be very useful to study musical instruments, notably their playability. However, in order to take into account all the potential behaviors of a system, their computation can be costly. As a result, the computation cost limits the possible applications of these methods, for parametric studies for example.

Initially, [31] proposed to use Monte-Carlo techniques to numerically estimate the basin stability. Some limitations of these techniques have been drawn [40] when applied on strange basins of attraction (e.g., fractal and riddled/intermingled). However, such basins remain to be exhibited on self-sustained musical instruments. Some authors also proposed enrichment to the metric, to take into account uncertainties and variations of the parameters [51, 52]. Additionally, improvements on classical Monte Carlo approaches, notably using machine learning tools, show promise to accelerate basin stability computation. This could be interesting for more complex musical instrument models. For example, [41] proposed a method based on support vector machines. The idea is to find the basins’ boundaries with a limited number of samples. Work still needs to be done to determine which method should be used depending on the situation and on the classification cost of a single sample. Indeed, if the classification of each sample is fast, it may be slower to find basins’ boundaries rather than to apply a simple Monte-Carlo size estimation. For instance, it is the case for the Van der Pol oscillator presented in this article.

If the transient duration statistics and the basin stability are both studied, the two metrics can be computed simultaneously, using the same time integrations (as in Sect. 5). The cost needed to compute the basin stability with Monte-Carlo methods only depends on the size of the basin as indicated by the absolute standard error given by equation (7). For systems with higher dimensions, that cost could increase if the classification cost does. Moreover, the relative standard error, which writes ${\mathrm{err}}_{\mathrm{rel}}\left({\widehat{\mathcal{S}}}_{\mathcal{B}}\right)=\sqrt{\frac{1/{\mathcal{S}}_{\mathcal{B}}-1}{N}}$, increases as ${\mathcal{S}}_{\mathcal{B}}$ decreases. Higher dimension systems may have more stable solutions which lead to smaller basin stability values. More samples would thus be needed to keep the relative standard error low. This standard error formula stands for independent trials with only two outcomes – the sample is either inside or outside the studied basin of attraction. The calculation of the transient duration distributions does not enter into this category since the outcome is a continuous random variable: a transient duration. With a sparse sampling, some transient behaviors that appear for specific initial conditions might be missed. Consequently, it could be interesting to estimate the convergence of the probability density functions in order to chose an adequate number of samples.

The statistical methods presented here can be applied applied at low computational cost on the fifth order Van der Pol oscillator, because the two attractors and the boundary of the basins of attraction are analytically known and characterized by their ${\mathcal{L}}_2$-norm. On more complex systems, the attractors and basins’ boundaries usually do not have analytic expressions. The main difficulty encountered in carrying out the statistical analysis of the trumpet model presented in Section 5 was to define adequate classification and stopping criteria. Indeed, non-linear systems sometimes evolve with very varied dynamics, making it difficult to distinguish between a steady state and a transient regime with slow dynamics. This kind of behavior leads to significant differences in transient duration between different initial conditions, especially at a bifurcation, where one eigendirection becomes infinitely slow. These discrepancies in transient durations are particularly visible in Figure 13b, around the Hopf bifurcation. In this regard, the Manlab solutions were very helpful to define efficient classification and convergence criteria (see Appendix B). For systems with a wide range of transient durations, we would strongly recommend to set a convergence criteria rather than to use a fixed integration time, as this is the case in the version of bSTAB currently online. This would improve performances and limit missclassification.

6.2 Initial conditions and musical gestures

In Section 3, the presented statistical approach is based on a uniform distribution of initial conditions. This is a general distribution and it allows different instrument configurations to be compared. However, the resulting probabilities probably fail to translate the musician’s playing experience. As a first attempt to tackle this issue, we study in Section 4 initial conditions resulting from an archetypal control gesture. This approach highlights the importance of the considered distribution of initial conditions and it can be interpreted in musical terms more directly than random initial conditions. We would like to draw the reader’s attention to the fact that the results obtained with this transient-based approach not only depend on the transient control parameters, but also on the choice of the pre-initial conditions. Overall, applicative studies on playability would greatly benefit from a preliminary description of the control scenarios and, where relevant, initial conditions distributions. In that perspective, it could be of great interest to measure transient blowing scenarios provided by humans (musicians and non-musicians) in an experimental study. Moreover, as illustrated by Figure 10, dynamical systems are more likely to cross specific regions of their phase space during transients. Choosing distributions of initial conditions based on the slow and fast eigendirections (depending on the solutions that exist during the transient) could also be an interesting idea.

In this paper, we focus on initial conditions induced by the fast evolution of one control parameter. However, musicians may also impact the initial conditions by acting directly on a state variable. One example among many others is the use of the tongue in reed instruments to impose initial reed positions and velocities. Measurements of tongue-induced reed positions in a clarinet are presented in [53] for example. To improve the statistical approaches, these initial conditions should also be taken into account.

6.3 Stability, transient duration and playability

Other studies were interested in predicting the transient duration of a self-sustained musical instrument [24, 25]. A common approach is to linearize the system around an unstable equilibrium and to consider the exponential growth of a perturbation. This exponential growth depends on the positive real part of the eigenvalues of the matrix of derivatives (i.e., the jacobian matrix). This method is better adapted for studying initial conditions close to an equilibrium solution. The Floquet theory might be used in the same way for initial conditions next to a periodic solution, although it has yet to be applied in that manner to a musical instrument model (for more details on the Floquet theory, refer to [54, 55]). In this article, we consider initial conditions not necessarily close to any solution. Hence, these linear analyses might not be valid for every sample. On more complex systems with several competing multistable solutions, its local nature prevents it from giving a complete analysis of the transients. However, the systematic and fast nature of linear analysis makes it a good candidate to bolster certain sections of a complete statistical study.

Finally, the statistical study of the trumpet model presented in Section 5 highlights several interests of the methods. The long transients observed in Figure 13 raise the question of the playability of a regime as a function of its transient duration. In a standard musical context, we argue that a regime that takes several seconds to establish cannot really be considered playable. Therefore, certain stable regimes could actually be unplayable due to transient behaviors, which shows the relevance of nuancing the notion of stability using transients when studying playability. The transient-based analysis of Section 5.3.2 also indicates that the transient control parameters of a musical gesture may have nontrivial optimal values. Future works could analyze these types of transient control scenarios into more details, relying on the statistical approaches as playability guidelines.

7 Conclusion and perspectives

In this article, we present a sample-based approach that can be used to enrich a bifurcation diagram. In this study we use it to investigate two playability issues: the prediction of the steady state in multistable situations and the prediction of the transient duration. The method used for the first issue is called the basin stability [31], whereas the method proposed for describing the responsiveness of a dynamical system is new, to the authors’ knowledge. These two methods rely on distributions of initial conditions that represent the transient action of a musician. The results strongly depend on this choice of distribution and we proposed a method based on time integration to generate transient-informed distributions.

The methods are then applied to a physical model of trumpet with the aim to evaluate their interest in a more practical case. The transient-based basin stability of this model present non trivial tendencies regarding some transient control parameters. Moreover, very long transients are highlighted, which raises the question of the playability of such asymptotic regime. Overall, these statistical methods show interest for the analysis of musical instruments, but they can also be extended to a large variety of systems.

Future works will focus on two aspects. First, these methods will be applied on musical instrument models in configurations showing multiple oscillating regimes and rich transients. Then, measurements on musicians will be conducted to evaluate the initial conditions that they induce, depending on the desired musical effect.

Conflicts of interest

The authors declare that they no conflicts of interest in relation to this article.

Data availability statement

No new data were created or analyzed in this study.

Appendix A

Trumpet model

We consider the brass instrument model described in [5]. The convention used to represent the lip position is given in Figure A1. The dimensioned and unregularized equations are as follows:

$$\{\begin{array}{l}\ddot{x}+\frac{{\omega }_L}{Q_L}\dot{x}+{\omega }_L^2(x-x_0)=\frac{p_0-p}{{\mu }_L},\\ {\dot{p}}_n-s_n{p}_n=Z_c{C}_{n}u \forall n\in [1,N_m],\\ p=2\sum _{n=1}^{N_m}\mathrm{}\mathfrak{R}\left(p_n\right),\\ u=W\sqrt{\frac{2|p_0-p|}{\rho }}\cdot \mathrm{sign}(p_0-p)\cdot \max (x,0),\end{array}$$(A1)

with

$$\{\begin{array}{l}{\omega }_L=\sqrt{\frac{k}{m}},\\ Q_L=\frac{\sqrt{\mathrm{km}}}{c},\\ {\mu }_L=\frac{m}{S}.\end{array}$$(A2)

S is the surface area of the lip to which the pressures p₀ and p are applied, μ_L is the surface mass of the lip, Q_L its quality factor and ω_L its angular eigenfrequency. The resonator is represented by its input impedance, which is treated as a sum of N_m modes with poles s_n and residues C_n. The lip allows air to pass through a rectangular surface of height x and width W. The system given in (A1) is scaled as follows:

$$\{\begin{array}{l}\ddot{\stackrel{̃}{x}}+\frac{{\omega }_L}{Q_L}\dot{\stackrel{̃}{x}}+{\omega }_L^2(\tilde{x}-1)={\omega }_L^2(\gamma -\tilde{p}),\\ {\dot{\stackrel{̃}{p}}}_n-s_n{\tilde{p}}_n=C_n\tilde{u} \forall n\in [1,N_m],\\ \tilde{p}=2\sum _{n=1}^{N_m}\mathrm{}\mathfrak{R}\left({\tilde{p}}_n\right),\\ \tilde{u}=\zeta \sqrt{|\gamma -\tilde{p}|}\cdot \mathrm{sign}(\gamma -\tilde{p})\cdot \max (\tilde{x},0).\end{array}$$(A3)

Figure A1 Schematic diagram of the one-degree-of-freedom lip model. Assume x = 0 when the lip is in the closed position and note x0 the position of the lip at rest (the situation shown here).

Variables affected by this scaling are marked with a ~ symbol and are defined as follows:

$$\{\begin{array}{l}\tilde{x}=\frac{x}{x_0},\\ p_M={\mu }_L{\omega }_L^2{x}_0,\\ \tilde{p}=\frac{p}{p_M},\\ \gamma =\frac{p_0}{p_M},\\ {\tilde{p}}_n=\frac{p_n}{p_M},\\ u_M=\frac{p_M}{Z_c},\\ \tilde{u}=\frac{u}{u_M},\\ \zeta =Z_{c}W{x}_0\sqrt{\frac{2}{\rho p_M}}.\end{array}$$(A4)

Irregular functions appearing in the flow term are regularized:

$$\{\begin{array}{l}|•{|}_r=\sqrt{{•}^2+ϵ},\\ \mathrm{sign}(•{)}_r=\frac{•}{|•{|}_r},\\ \max (•,0{)}_r=\frac{•+|•{|}_r}{2}.\end{array}$$(A5)

The regularization term ϵ is arbitrarily set to ϵ = 10⁻⁶.

Time is scaled by the first modal angular frequency ${\omega }_1=\mathfrak{I}\left(s_1\right)$ and dimensionless time derivatives are denoted ′. The time dimensionless poles, residues and lip angular frequency respectively write ${\widehat{C}}_n,{\widehat{s}}_n \mathrm{and} {\widehat{\omega }}_L$. Finally, we separate the real R_n and imaginary I_n parts of the modal pressures p_n, and the system we solve is as follows:

$$\{\begin{array}{c}\tilde{x}″+\frac{{\widehat{\omega }}_L}{Q_L}\tilde{x}\text{'}+{\widehat{\omega }}_L^2(\tilde{x}-1)={\widehat{\omega }}_L^2(\gamma -\tilde{p}),\\ \begin{array}{c}{\tilde{R}}_n^{\text{'}}=\mathfrak{R}\left({\widehat{s}}_n\right){\tilde{R}}_n-\mathfrak{I}\left({\widehat{s}}_n\right){\tilde{I}}_n+\mathfrak{R}\left({\widehat{C}}_n\right)\tilde{u}\\ \begin{array}{cc}& \hspace{1em}\forall n\in [1,N_m],\end{array}\end{array}\\ \begin{array}{c}\begin{array}{c}{\tilde{I}}_n^{\text{'}}=\mathfrak{R}\left({\widehat{s}}_n\right){\tilde{I}}_n+\mathfrak{I}\left({\widehat{s}}_n\right){\tilde{R}}_n+\mathfrak{I}\left({\widehat{C}}_n\right)\tilde{u}\\ \begin{array}{cc}& \forall n\in [1,N_m],\end{array}\end{array}\\ \begin{array}{c}\tilde{p}=2\sum _{n=1}^{N_m}\mathrm{}\mathfrak{R}\left({\tilde{p}}_n\right),\\ \tilde{u}=\zeta \sqrt{|\gamma -\tilde{p}{|}_r}\cdot \mathrm{sign}(\gamma -\tilde{p}{)}_r\cdot \max (\tilde{x},0{)}_r.\end{array}\end{array}\end{array}$$(A6)

The parameter values chosen for this article are as follows:

$$\{\begin{array}{l}{\omega }_L=2\pi \times 200 \mathrm{rad}·{\mathrm{s}}^{-1},\\ {\mu }_L=2 \mathrm{kg}·{\mathrm{m}}^{-2},\\ Q_L=3,\\ x_0=1·{10}^{-4} \mathrm{m},\\ W=8·10^{-3} \mathrm{m},\\ \mathrm{Zc}=1.83 \mathrm{kg}·{\mathrm{s}}^{-1}·{\mathrm{m}}^{-4},\\ ϵ=10^{-6}.\\ \end{array}$$(A7)

Appendix B

Modifications to bSTAB

The bSTAB toolbox [48] has been designed to calculate the basin stability of any dynamic system as automatically as possible. It generates initial conditions, performs a time integration for each of them, then classifies the obtained regime by comparing it to reference signals. In our case, we compare the peak-to-peak amplitude of R₂, the real part of the second modal pressure, to that of the MANLAB solutions. Since classification must be performed on the steady part of the regime, the time integration continues until the following convergence criteria are met.

During the last duration t^⋆:

1. The amplitude of R₂ is within ε of that of a reference solution;

2. The difference between the amplitude of R₂ and that of the reference decreases.

The transient duration is defined as the duration after which the envelope of R₂ stays within ε of the amplitude of the reference solution. In this article, we choose ε = 0.1.

The computation of the transient duration is illustrated in Figure B1a. The corresponding trajectory is also represented in Figure B1b, with its projection along x, R₂ and R₄.

Figure B1 (a) Computation of the transient duration (represented by the dotted line) and (b) projection of the corresponding trajectory along x, R2 and R4. The red dot is the initial state, the blue cross is the equilibrium, the red line is the stable periodic solution and the black broken line is the unstable periodic solution.

References

All Tables

All Figures

	Figure 1 Bifurcation diagram of a fifth order Van der Pol oscillator given in equation (1) with respect to the control parameter μ, with σ = −1.5 and ν = 0.1. Stable solutions are represented with continuous lines, unstable solutions in broken lines and bifurcations are indicated with star markers. The region of multistability is shaded in grey.
In the text

	Figure 2 Illustration of the sensitivity to the initial conditions, for the Van der Pol oscillator, regarding the steady state regime and the transient duration. On the left side, the phase plane of the system for μ = −3, with trajectories coming from three close initial states (indicated with star markers). On the right side, the evolution of the position x with respect to time (dimensionless).
In the text

	Figure 3 Basins of attraction of the fifth order Van der Pol oscillator (μ = −3). Initial conditions taken in the red shaded area give rise to the periodic solution whereas initial conditions taken in the blue shaded area give rise to the equilibrium solution.
In the text

	Figure 4 Illustration of the classification method in the general case. The plotted trajectories are obtained with time integrations starting from the stars as initial conditions. Red trajectories are inside the basin of attraction of the stable periodic solution X2π+ whereas the blue ones belong to the basin of the equilibrium XEq.
In the text

Figure 5 Basin stability of the two stable solutions computed from 100 samples. In red: basin stability of the periodic solution. In blue: basin stability of the equilibrium. Error bars correspond to the standard error of the size estimation of the Monte-Carlo method due to sub-sampling equation (7). In black: basin stability computed with the analytic expression of the basin boundary (circle of known diameter X2π− given in (3)). The region of multistability is shaded in grey.

In the text

	Figure 6 (a) Phase trajectories and (b) time evolution of the associated norm of the system’s state for 300 random initial conditions. Blue trajectories lead to the equilibrium and red trajectories lead to the periodic solution. The neighborhoods of the two stable solutions are shaded in grey. They are defined at a distance ε = 0.1 in ${\mathcal{L}}_2$-norm from the stable solutions.
In the text

Figure 7 (Top) Cumulative probability function of the transient duration distribution. In ordinate, the proportion of trajectories with a transient duration below the value given in abscissa. (Bottom) Corresponding violin plots. The box plots indicate the extreme values, the median, the upper and lower quartiles whereas the curves around represent the probability density function of the transient duration distribution. A latin hypercube of 1000 samples is used and the transient durations are defined with ε = 0.1.

In the text

	Figure 8 Violin plots of the transient duration with respect to the control parameter μ. In blue: equilibrium. In red: periodic solution. Number of samples: N = 1000. Panel (b) is a magnified version of the transients leading to the periodic solution. The region of multistability is shaded in grey.
In the text

	Figure 9 Control scenarios used to compute the transient-informed initial conditions. The multistability region is shaded in gray. μ0 is the initial value of μ, μmax is the overshoot, μend is the final value of μ, at which the basin stability is computed, μH and μSN are the value of μ at the Hopf and at the saddle-node bifurcations and τatt is the attack duration.
In the text

	Figure 10 Example of a transient-informed distribution generated with 10 control scenarios of different durations. (•) Pre-initial Latin Hypercube of size N0 = 10; () states of the system when μ = μmax; () transient-informed distribution at the end of the attack. The black dotted line is the fast eigendirection of the equilibrium at μ = μmax and the red one is its slow eigendirection at μ = μend.
In the text

	Figure 11 Basin stability of the periodic solution computed with different distributions. The red curves are computed with transient-informed distributions of initial conditions, each curve corresponding to a specific overshoot μmax. The black curve is computed with a uniform distribution (it is the same as in Fig. 5 without the error bars).
In the text

	Figure 12 Bifurcation diagram of the 11-mode trumpet model, calculated with Manlab. Top pane: amplitude of the oscillations ${\Vert p\Vert }_{{\mathcal{L}}_2}$ versus blowing pressure p0. Bottom pane: fundamental frequency of the oscillations fplay versus blowing pressure p0. Stable solutions are shown as solid lines, unstable solutions as dotted lines. Bifurcations are represented by stars.
In the text

	Figure 13 Basin stability (a) and transient duration distributions (b) of the trumpet model calculated with a uniform ${\mathcal{Q}}_u$ distribution of N = 300 initial conditions. The dotted lines in (b) indicate the limit of the scale used for the zoomed view in (c). Error bars are given by the equation (7).
In the text

Figure 14 Basin stability of the trumpet model calculated with distributions ${\mathcal{Q}}_{\mathrm{tb}}$ of N = 1000 initial conditions based on attack transients (10 scenarios applied to each of the N0 = 100 pre-initial conditions). (a) Each curve corresponds to a specific attack duration τatt. (b) Each curve corresponds to a specific overshoot value p0,max. Only the basin stability of the periodic regime is shown.

In the text

	Figure A1 Schematic diagram of the one-degree-of-freedom lip model. Assume x = 0 when the lip is in the closed position and note x0 the position of the lip at rest (the situation shown here).
In the text

	Figure B1 (a) Computation of the transient duration (represented by the dotted line) and (b) projection of the corresponding trajectory along x, R2 and R4. The red dot is the initial state, the blue cross is the equilibrium, the red line is the stable periodic solution and the black broken line is the unstable periodic solution.
In the text