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

HAL
open science

# Nuancing stability for self-sustained musical instruments: transient durations and basins of attraction

Martin Pégeot, Christophe Vergez, Jean-Baptiste Doc, Vincent Fréour, Tom Colinot

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

► To cite this version:

Martin Pégeot, Christophe Vergez, Jean-Baptiste Doc, Vincent Fréour, Tom Colinot. Nuancing stability for self-sustained musical instruments: transient durations and basins of attraction. Forum Acusticum, Sep 2023, Turin, Italy. 10.61782/fa.2023.0180. hal-04800154

HAL Id: hal-04800154

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

Submitted on 23 Nov 2024

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# NUANCING STABILITY FOR SELF-SUSTAINED MUSICAL INSTRUMENTS: TRANSIENT DURATIONS AND BASINS OF ATTRACTION

Martin Pégeot¹*

Christophe Vergez¹

Jean-Baptiste Doc²

Vincent Fréour³

Tom Colinot⁴

¹ Aix Marseille Univ, CNRS, Centrale Marseille, LMA UMR7031, Marseille, France

² Laboratoire de Mécanique des Structures et Système Couplés, Conservatoire National des Arts et Métiers, Paris, France

³ YAMAHA Corporation, Research and Development Division, Hamamatsu, Japan

⁴ Buffet Crampon, Mantes-La-Ville, France

## ABSTRACT

We are interested in the behavior of self-sustained musical instruments, such as wind instruments and bowed string instruments. Those dynamical systems can be modeled by autonomous equations. The quasi-static analysis of those equations is thoroughly represented in the literature. In this framework, stationary solutions are calculated, along with their local stability. However, those systems may be multistable: multiple locally stable solutions coexist for identical parameter values. Transient analysis notions are then needed to predict the actual regime obtained. This work aims at proposing analysis tools and graphical representations to complement the quasi-static analysis. On the one hand, we nuance the notion of stability using basins of attraction (determined with Support Vector Machines). On the other hand, we propose to enrich bifurcation diagrams with information about transient durations. Those methods are applied to a fifth order Van der Pol oscillator, which is an prototype of a self-sustained musical instrument model.

**Keywords:** Self-sustained musical instruments, Stability analysis, Basins of attraction, Transient durations

## 1. INTRODUCTION

Self-sustained musical instruments sometimes have multistable behavior, that is to say that multiple locally stable solutions can coexist for a given set of parameters [1, 2]. As a result, it may be both hard to obtain the desired sound for the musician and to predict the stationary regime of the system for the scientist.

In this paper, we propose to quantify the stability of a solution – related to its ease of production – by determining the size of its basin of attraction. This notion is known in the non-linear dynamics community as the basin stability [3]. To estimate the basin stability of a solution, we propose a two-step method. First, we determine the geometry of its basin of attraction using Explicit Design Space Decomposition (EDSD). It is a method based on Support Vector Machines (SVM) developed by [4]. Then, we use a Monte-Carlo method relying on the trained SVM as surrogate model, to compute the size of the basin of attraction. In the end, this indicator can be used to classify the regimes of a musical instrument by ease of production.

The playability of a regime also depends on the duration of its transient regime [5]. Indeed, a stable solution which takes too long to establish will hardly be usable musically. An instrument with long transient durations would be considered poorly made by experienced musicians. Thus, we present preliminary ideas for studying the transient regime.

<footnote>*Corresponding author: pegeot@lma.cnrs-mrs.fr.
Copyright: ©2023 Martin Pégeot et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 Unported License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.</footnote>
Those methods are illustrated on a fifth order Van der Pol oscillator [6], which presents similarities with models of self-sustained musical instruments. Its main advantage is that its phase space is of dimension 2, which facilitates the exploration of new analysis methods and concepts.

This paper is organized as follows. The system studied is presented in Sec.2. The notion of basin stability and its computation are then described in Sec.3. Basin stability and transient analysis are illustrated for this system in Sec.4 and conclusions are given in Sec.5.

## 2. STUDY SYSTEM

### 2.1 Choice of the system

In order to better highlight the interest of our analysis methods on self-sustained musical instruments, we study a fifth order Van der Pol oscillator [6]. It is a normal form of any system of ODE's which starts oscillating through the loss of stability of the equilibrium solution (Hopf bifurcation) [7]. It means that even a model of self-sustained musical instrument including many degrees of freedom (for instance many resonance modes) can be represented around the Hopf bifurcation by this Van der Pol oscillator. Since we are interested in hysteresis and bi-stability phenomenon at the oscillation threshold, this normal form has to be order-five so that it can show for adequate parameter values an inverse Hopf bifurcation followed by a saddle-node bifurcation [6].

The authors precise that this system has been chosen for its simplicity in this exploratory phase. In the long term, the goal is to apply them on complete self-sustained musical instrument models. Such models are presented in [1,2] (reed instrument models with a modal description of the resonator).

### 2.2 Equation of motion

The oscillator considered is driven by the following differential equation

$$ \ddot{x} = [\mu - \sigma(x^2 + \dot{x}^2) - \nu(x^2 + \dot{x}^2)^2] \dot{x} - x, \quad (1) $$

where $x$ and $\dot{x}$ denote the position and the velocity of the oscillator. $\sigma$ and $\nu$ are two parameters which are fixed $\sigma = -1.5$ and $\nu = 0.1$ in order to reproduce the behavior exhibited by musical instruments (bistability induced by an inverse Hopf bifurcation and a fold). $\mu$ is the control

parameter in this study. In the context of wind musical instruments, it can be interpreted as the blowing pressure that causes sound emergence trough linear instability [1].

Analytical solutions of Eqn. (1) are given in [6]. They write

$$ x = X \cos(t + \varphi), \quad (2) $$

with

$$ X = 0 \quad \text{or} \quad X = \sqrt{\frac{-\sigma \pm \sqrt{\sigma^2 + 4\mu\nu}}{2\nu}}. \quad (3) $$

The equilibrium is always defined. The periodic solution of largest amplitude is defined above a saddle-node bifurcation and the periodic solution of smallest amplitude is defined between the saddle-node bifurcation and an inverse Hopf bifurcation. The saddle-node bifurcation occurs at $\mu_{SN} = -\frac{\sigma^2}{4\nu}$ and the Hopf bifurcation at $\mu_H = 0$.

### 2.3 Bifurcation diagram

To study the behavior of such a system, it is possible to plot its bifurcation diagram. A bifurcation diagram depicts the different solutions in steady state, their amplitude and their local stability, with respect to the control parameter. For systems where no analytical solution is known, the authors use the software Matlab [8] to compute equilibrium and periodic solutions. This software combines the harmonic balance method approximation with Taylor expansion over the curvilinear abscissa. With this software, a great number of harmonics can be used to represent periodic solutions, thus leading to good approximations. More details can be found in [8].

For the system described by Eqn. (1), Manlab solutions are identical to the analytical solutions (thus only the first harmonic is not null). Its bifurcation diagram is given in Fig.1. The amplitude of the oscillations are represented by the $L_2$ norm of the state vector $||X||_{L_2}$, which also represents the mechanical energy of the system. The local stability of the solutions is given by the lifestyle. It is computed with linear stability analysis for the equilibrium [9], and with Floquet theory for periodic solutions [10]. The two bifurcations can be observed: an inverse Hopf bifurcation between stable equilibrium, unstable equilibrium and unstable periodic solution, and a saddle-node bifurcation between stable and unstable periodic solutions. This behavior is typical of what can be observed on the saxophone, for instance, especially for
low frequency notes [1].

Between those two bifurcation points, the system is multistable: the equilibrium and the largest periodic solution are both locally stable (but none is globally stable). This region could be interesting for a musician because the periodic solution has a lower amplitude than for greater values of $\mu$. In order to play pianissimo, the musician might thus target this region.

To produce sounds in this area, two strategies are possible. The first one is to set $\mu > 0$ until a sound appears and then to reduce $\mu < 0$ to its desired value. The second one is to set $\mu < 0$ to its desired value and then to perturb the system by imposing specific initial conditions. For a saxophone, initial position and velocity of the reed $(x, \dot{x})$ could be imposed with a flick of the tongue. In both cases, once the pianissimo regime is reached, a large enough perturbation could set the system back to the equilibrium.

This musical example illustrates the importance of quantifying the stability of a solution. Such a quantification would indicate how easy it is to produce and maintain the desired regime in a multistability situation. In the following section, a basin stability approach will be introduced to quantify the stability of a solution.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Bifurcation diagram of the system studied. In abscissa is the control parameter $\mu$ and in ordinate is the 2-norm of the state vector $||X||_{L_2}$. Solid lines represent locally stable solutions and dotted lines represent unstable solutions.</figcaption></figure>

# 3. BASIN STABILITY

## 3.1 Principle

The basin stability method was first introduced under those terms by [3] and has widely been used since then in the non-linear dynamics community [11–14]. It has been developed to overcome the shortcomings of the linear stability analysis regarding multistable situations i.e. its incapacity to estimate the resilience of solutions to a non-small perturbation. The idea is to characterize the stability of a solution by the size of its basin of attraction. The basin stability of a solution indicates the probability to asymptotically return to this solution given a random not necessarily small perturbation.

Considering $\mathcal{Q}$ a subset of the state space with a finite volume, $A$ the solution of interest ($A \subset \mathcal{Q}$), $\mathcal{B}$ its basin of attraction ($\mathcal{B} \subset \mathcal{Q}$) and $\rho$ a probability density function, the basin stability $S_B$ is defined as follows

$$S_B = \int_{\mathcal{Q}} \mathbf{I}_B(\mathbf{X}) \rho(\mathbf{X}) d\mathbf{X}, \quad (4)$$

where

$$\mathbf{I}_B(\mathbf{X}) = \begin{cases} 1, & \text{if } \mathbf{X} \in \mathcal{B} \\ 0, & \text{otherwise.} \end{cases} \quad (5)$$

If there is knowledge on the expected perturbation, the probability density function $\rho$ can follow a specific distribution. A uniform distribution should be used otherwise. In any case $\int_{\mathcal{Q}} \rho(\mathbf{X})d\mathbf{X} = 1$, so $S_B \in [0, 1]$.

To estimate $S_B$ numerically, which will be noted $\hat{S}_B$, one can use a Monte-Carlo technique to sample randomly the subset $\mathcal{Q}$. The basin stability estimation $\hat{S}_B$ is then given by

$$\hat{S}_B = M/N, \quad (6)$$

where $M$ is the number of samples belonging to $\mathcal{B}$ and $N$ is the total number of samples.

As pointed out in [3, 13, 14], this experiment corresponds to $N$ independent Bernoulli trials with probability of success $S_B$, leading to the absolute standard error for $\hat{S}_B$ due to sampling

$$\text{err}(\hat{S}_B) = \sqrt{\frac{S_B(1 - S_B)}{N}}. \quad (7)$$

Error due to misclassification of the samples should be added to this. The classification method used in this study is described in the remaining of this section.
## 3.2 Classification of the initial conditions

In this paragraph, we present the criteria chosen to classify the samples, that is to say to determine to which basin of attraction belongs a sample. For the system of this study, the basins of attraction are analytically known and their boundary is the unstable periodic solution [6]. However, this method is meant to be applied on more complex systems for which basins of attraction might be unknown. Usually, a criteria based on a time integration is used [15]. In that respect, we have chosen that type of criteria to test our method.

For the scaled Van der Pol oscillator of Eqn. (1), the time evolution of the $L_2$ norm of the system is monotonic. While the system get attracted by one stable solution or another, its $L_2$ norm tends toward a characteristic value of the attractor. This behavior is visible on Fig.4. All the red trajectories tend toward the stable periodic solution which norm is approximately 3.5, and all the blue trajectories tend toward the stable equilibrium which norm is equal to 0.

To classify a set of initial conditions, we compute the evolution of the system with finite differences (we use the Matlab function ode45). The integration time is arbitrarily chosen to reach steady state. By the end of the numerical integration, if the 2-norm is smaller than an arbitrarily small value ($\epsilon = 0.1$ in this study), then the sample is classified as belonging to the basin of attraction of the equilibrium. Otherwise, it is classified as belonging to the basin of attraction of the stable periodic solution.

Depending on the dynamics of the system and on its transient durations, the integration time needed to reach the steady state might be long. Classifying a great number of samples would thus be costly. In order to increase the number of samples without increasing the calculation time, we propose to build a low-cost classification function with a method based on Support Vector Machines (SVM). This classification function captures the geometry of the basins of attraction.

## 3.3 Computation of the basin geometry

The computation of the basin geometry is a classification problem and Support Vector Machines (SVM) are very adapted for that kind of problematic. Thus different methods using SVM to determine basin boundaries can be found in the literature [16, 17]. Here, we propose to

use the Explicit Design Space Decomposition (EDSD) method [4].

This is an iterative method. At each iteration, the boundary is first estimated out of few samples, using SVM. New samples are then added on the estimated boundary while being as far as possible from other samples. This two steps are then repeated until the estimated boundary has converged (see [4] for details on the convergence criteria). In the end, one has a function $f: Q \to R$ which is positive if the sample to evaluate belongs to the basin of attraction of interest, and negative otherwise. This evaluation has a very low computational cost and a high number of samples can thus be used to estimate the size of the basin. For this size estimation, we use a Latin hypercube sampling (LHS). To train the SVM, the more costly evaluation function presented in the previous sub-section is used.

In general, there is no error measure available for the boundary estimated with EDSD [4]. Later on, the error of the estimated basin stability will only take into account the error coming from sub-sampling. Its expression is given in Eqn. (7).

# 4. ILLUSTRATION

## 4.1 Basin stability analysis

The basin stability of the equilibrium point of the system has been undertaken using the methods presented in Sec.3 (namely, the EDSD to compute the basin geometry and a LHS to compute its relative size). The results are given in Fig.2, for different values of $\mu$ inside the bistability area. The basin stability of the periodic solution is the complementary of the basin stability of the equilibrium because they are the only two stable solutions. The region of interest is $x \in [-5, 5]$ and $\dot{x} \in [-5, 5]$ because the basin of attraction of the equilibrium point is included in it for all considered values of $\mu$. Standard error bars due to down-sampling are represented above result bars (reminder: error coming from EDSD is not taken into account). This error remains low because $N = 10^5$ samples have been used to calculate the size of the basin of attraction. That is the point of having a two-step method (EDSD first and LHS with surrogate model then).

For this system, the size of the basin of attraction of the equilibrium point can be analytically calculated
because its boundary coincides with the unstable periodic solution. This unstable solution is a circle in the phase space. Its radius is given in Eqn. (3) (smallest non zero value of X). The actual basin stability $S_B$ is superimposed over the estimated one $\hat{S}_B$.

One can observe that the estimation method is correct in this situation. This means that the error due to EDSD remains small here. Moreover, one notices that the equilibrium point is less stable when $\mu$ increases. Eventually, $S_B$ reaches 0 at the Hopf bifurcation (for $\mu_H = 0$), because the equilibrium point becomes unstable. Similarly, $S_B$ would jump to 1 at the saddle-node bifurcation because the equilibrium point becomes the only stable solution. Its basin of attraction suddenly transforms from a circle included into the region of interest to the whole region. This explains the discontinuity of $S_B$ at the saddle-node bifurcation (which is not represented to keep a readable figure size in the bistability region).

To give a little bit of context, this $S_B$ curve shape means that musicians would need to disturb greatly their instrument (with a powerful tongue strike for instance) to produce a sound when they blow softly ($\mu \approx -5$ for instance). On the contrary, only a small disturbance would be needed if their mouth pressure is high ($\mu \approx 0$). One speculative interpretation is to conclude that it is harder to attack a note pianissimo than mezo-piano on this "instrument". This conclusion must be treated with care, as we do not know which gesture – and therefore which initial conditions – is easy for a musician to execute.

With this metric, the stable solutions of a musical instrument can be classified from the easiest to produce to the hardest. Moreover, its value can be interpreted as the probability to produce one regime or another, without specific control from the musician in terms of initial conditions.

## 4.2 Transient analysis

This section presents preliminary ideas to study the transient regime.

Now that the likelihood of appearance of stationary solutions has been estimated, one can wonder how the system will behave during its transient regime and how long it will last. To study the global features of the transient behavior of the system, we propose to observe

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Basin stability estimation of the equilibrium (bars) from Eqn. (6) versus actual basin stability (dotted line) from Eqn. (4). For the estimation, basin geometry is computed with EDSD and basin size is computed with a latin hypercube ($10^5$ samples).</figcaption></figure>

its state evolution starting from initial states close to the boundary of the two basins of attraction.

In Fig.3, we plot in the phase space the three stationary solutions along with multiple trajectories initialized from one side and the other of the unstable limit cycle (boundary of the basins of attraction). It appears that some regions of the phase space are more likely visited than others.

In Fig.4, the time evolution of the 2-norm of the state vector $||X||_{L_2}$ is given for all those trajectories. It highlights the influence of the initial phase conditions over the transient duration. Indeed, for an identical mechanical energy given to the system ($||X||_{L_2}$ is the mechanical energy of the system), different initial phase conditions lead to different transient durations.

## 5. CONCLUSION

In this paper, we quantified the stability of the solutions of a multistable system, using the notion of basin stability. This metric allows us to rank the solutions in order of stability – i.e. in order of ease of production. Moreover, its value has a meaning: it corresponds to the probability of the solution to happen, given ordinary initial conditions. From the authors perspective, this notion seems particularly promising for describing musical instrument
<figure><img src="image_6.png" /></figure>

**Figure 3.** Phase portrait of the system studied. The thick black lines and the black dot in the center correspond to the solutions. Stability is given by the linestyle (straight if stable and dotted if unstable). The thin lines are trajectories initialized from one side (red) and the other (blue) of the unstable periodic solution.

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

**Figure 4.** Time evolution of the 2-norm of the state vector $||X||_{L^2}$ for trajectories depicted in Fig.3. The unstable periodic solution is represented by the dotted line and the two stable solutions correspond to the upper and lower boundaries of the graph.

behaviors.

To compute basin stability, we combined SVM methods – for basin geometry calculation – with Monte-Carlo techniques – for basin size calculation. Furthermore,

knowledge on basin boundary geometry helped us to have an overview of the transient regime by initializing time integrations along those boundaries. Main tendencies such as most likely transient durations and regions of the phase space have thus been highlighted.

## 6. REFERENCES

[1] 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.

[2] K. Takahashi, H. Kodama, A. Nakajima, and T. Tachibana, "Numerical Study on Multi-Stable Oscillations of Woodwind Single-Reed Instruments," *Acta Acustica united with Acustica*, vol. 95, pp. 1123–1139, Nov. 2009.

[3] P. J. Menck, J. Heitzig, N. Marwan, and J. Kurths, "How basin stability complements the linear-stability paradigm," *Nature Physics*, vol. 9, pp. 89–92, Feb. 2013. Number: 2 Publisher: Nature Publishing Group.

[4] 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, Oct. 2010.

[5] S. M. Logie, S. Bilbao, J. P. Chick, and D. M. Campbell, "The Influence of Transients on the Perceived Playability of Brass Instruments," 2010.

[6] D. Dessi, F. Mastroddi, and L. Morino, "A fifth-order multiple-scale solution for Hopf bifurcations," *Computers & Structures*, vol. 82, pp. 2723–2731, Dec. 2004.

[7] A. Algaba, E. Freire, E. Gamero, and C. García, "New aspects of the orbital normal form of the Hopf singularity: The Rayleigh and the van der Pol forms," *International Journal of Non-Linear Mechanics*, vol. 105, pp. 20–26, Oct. 2018.

[8] L. Guillot, B. Cochelin, and C. Vergez, "A Taylor series-based continuation method for solutions of dynamical systems," *Nonlinear Dynamics*, vol. 98, pp. 2827–2845, Dec. 2019.

[9] L. Velut, C. Vergez, J. Gilbert, and M. Djahanbani, "How Well Can Linear Stability Analysis Predict the Behaviour of an Outward-Striking Valve Brass Instrument Model?," *Acta Acustica united with Acustica*, vol. 103, pp. 132–148, Jan. 2017.
[10] L. Guillot, A. Lazarus, O. Thomas, C. Vergez, and B. Cochelin, "A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems," *Journal of Computational Physics*, vol. 416, p. 109477, Sept. 2020.

[11] P. Brzeski, M. Lazarek, T. Kapitaniak, J. Kurths, and P. Perlikowski, "Basin stability approach for quantifying responses of multistable systems with parameters mismatch," *Meccanica*, vol. 51, pp. 2713-2726, Nov. 2016.

[12] P. Brzeski, J. Wojewoda, T. Kapitaniak, J. Kurths, and P. Perlikowski, "Sample-based approach can outperform the classical dynamical analysis - experimental confirmation of the basin stability method," *Scientific Reports*, vol. 7, p. 6121, July 2017. Number: 1 Publisher: Nature Publishing Group.

[13] P. Schultz, P. J. Menck, J. Heitzig, and J. Kurths, "Potentials and limits to basin stability estimation," *New Journal of Physics*, vol. 19, p. 023005, Feb. 2017. Publisher: IOP Publishing.

[14] Y. Che, C. Cheng, Z. Liu, and Z. J. Zhang, "Fast basin stability estimation for dynamic systems under large perturbations with sequential support vector machine," *Physica D: Nonlinear Phenomena*, vol. 405, p. 132381, Apr. 2020.

[15] M. Stender and N. Hoffmann, "bSTAB: an open-source software for computing the basin stability of multi-stable dynamical systems," *Nonlinear Dynamics*, vol. 107, pp. 1451-1468, Jan. 2022.

[16] A. R. Armiyoon and C. Q. Wu, "A novel method to identify boundaries of basins of attraction in a dynamical system using Lyapunov exponents and Monte Carlo techniques," *Nonlinear Dynamics*, vol. 79, pp. 275-293, Jan. 2015.

[17] D.-W. Lee and J.-W. Lee, "Estimating Basin of Attraction for Multi-Basin Processes Using Support Vector Machine," *Management Science and Financial Engineering*, vol. 18, pp. 49-53, May 2012.