Baptiste Bergeot
email: baptiste.bergeot@centrale-marseille.fr

André Almeida

Christophe Vergez

Bruno Gazengel

Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure

Keywords: Musical acoustics, Clarinet-like instruments, Iterated maps, Dynamic Bifurcation, Bifurcation delay, Transient processes

teaching and research institutions in France or abroad, or from public or private research centers.

Prediction of the dynamic oscillation threshold in a

clarinet model with a linearly increasing blowing pressure

Introduction

One of the interests of mathematical models of musical instruments is to be able to predict certain characteristics of the produced sound given the gesture performed by the musician. In the case of a clarinet for instance, the amplitude, frequency or spectral content (the sound parameters) can be to a certain extent, determined as a function of the blowing pressure and lip force applied to the reed (the control parameters). A basic model, such as the one introduced by Wilson and Beavers [START_REF] Wilson | Operating modes of the clarinet[END_REF], allows to compute the amplitude of the oscillating resonator pressure from the knowledge of these two control parameters, giving results that follow the major tendencies observed in experiments. Several degrees of refinement can be added to this model, usually aiming at realistic sound and mechanical behavior. Well known simplifications of this model allow to study analytically the behavior of the clarinet. Simplified models, of course, are unable to describe or predict with refinement the exact harmonic content of the sound, or the influences of such important details as the reed geometry and composition or the vocal tract of the player. However, they can provide an understanding of the factors essential for the production of sound.

The highest degree of simplification of the model (introduced in Section 2) considers a straight, lossless (or losses independent of frequency) resonator and the reed as an ideal spring [START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF][START_REF] Kergomard | Elementary considerations on reedinstrument oscillations[END_REF][START_REF] Chaigne | Instruments à anche[END_REF]. With these assumptions, the system can be simply described by an iterated map [START_REF] Mcintyre | On the oscillations of musical instruments[END_REF]. Iterated maps often describe a succession of different regimes with variable periodicity. By analyzing the asymptotic values of these regimes it is possible to estimate: thresholds of oscillation, extinction, beating regimes, etc. [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF], amplitudes and stability of the steady state regime [START_REF] Ollivier | Idealized models of reed woodwinds. part 2 : On the stability of two-step oscillations[END_REF] and phenomena of period doubling [START_REF] Kergomard | Period doubling on cylindrical reed instruments[END_REF][START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF].

These characteristics arise from the so-called static bifurcation theory assuming that control parameters are constant. For example, these studies allow to find a static oscillation threshold γ st [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF] such that a constant regime is stable if the blowing pressure is below γ st and a periodic regime is stable if it is above γ st . More precisely, the oscillation emerges through a flip bifurcation [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF]. This behavior is static, obtained by choosing a constant blowing pressure, letting the system reach its final state, and repeating the procedure for other constant blowing pressures. Therefore, most studies using iterated map approach are restricted to a steady state analysis of the oscillation, even if transients are studied. They focus on the asymptotic amplitude regardless of the history of the system.

During a note attack transient the musician varies the pressure in her/his mouth before reaching a quasi-constant value. During this transient the blowing pressure cannot be regarded as constant. In a mathematical point of view increasing the control parameter (here the blowing pressure) makes the system non-autonomous and results from static bifurcation theory are not sufficient to describe its evolution. Indeed, it is known that, when the control parameter varies, the bifurcation point -i.e. the value of the blowing pressure where the system begins to oscillate -can be considerably delayed [START_REF] Kapral | Bifurcation structure of the nonautonomous quadratic map[END_REF][START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Fruchard | Sur le retard à la bifurcation[END_REF]. Indeed, the bifurcation point is shifted from γ st to a larger value γ dt called dynamic oscillation threshold. This phenomenon called bifurcation delay is not predicted by the static theory. Therefore, when the control parameter varies, results from the so-called dynamic bifurcation theory are required to properly analyze the system.

The purpose of this paper is to use results from dynamic bifurcation theory to describe analytically a simplified clarinet model taking into account a blowing pressure that varies linearly with time. In particular we propose a theoretical estimation of the dynamic oscillation threshold. Section 2 introduces the simplified mathematical model of a clarinet and the iterated map method used to estimate the existence of the oscillations inside the bore of the clarinet. Some results related to the steady state are presented in this section. Section 3 is devoted to the study of the dynamic system that takes into account a linearly increasing blowing pressure. The phenomenon of bifurcation delay is demonstrated using numerical simulations. A theoretical estimation of the dynamic oscillation threshold is also presented and compared with numerical simulations. In Section 4 the limits of this approach are discussed. It is shown, when the model is simulated, that the precision (the number of decimal digits used by the computer) has a dramatic influence on the bifurcation delay. The influence of the speed at which the blowing pressure is swept is also discussed.

y(t) -H 0 U (t) U r (t) U in (t) P (t)

Mouthpiece

State of the art 2.1 Elementary model

The model of the clarinet system used in this article follows an extreme simplification of the instrument, which can be found in other theoretical works [START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF][START_REF] Chaigne | Instruments à anche[END_REF].

This basic model separates the instrument into two functional elements. One of these is the bore, or resonator, a linear element where the pressure waves propagate without losses. The other is the reed-mouthpiece system, which is considered as a valve controlled by the pressure difference between the mouth and the mouthpiece. It is often called the generator and is the only nonlinear part of the instrument. A table of notation is provided in Appendix A.

The reed-mouthpiece system

The reed-mouthpiece system is depicted in Fig. 1. The reed is assumed to behave as an ideal spring characterized by its static stiffness per unit area K s . So, its response y to the pressure difference ∆P = P m -P is linear and is given by:

y = - ∆P K s . (1)
From (1) we can define the static closing pressure P M which corresponds to the lowest pressure that completely closes the reed channel (y = -H):

P M = K s H. (2)
The reed model also considers that the flow created by the motion of the reed U r is equal to zero, so that the only flow entering the instrument is created by the pressure imbalance between the mouth and the bore:

U in = U. (3)
The non-linearity of the reed-mouthpiece system is introduced by the Bernoulli equation which relates the flow U to the acoustic pressure P [START_REF] Hirschberg | A quasi-stationary model of air flow in the reed channel of single-reed woodwind instruments[END_REF][START_REF] Hirschberg | Aero-acoustics of wind instruments[END_REF]. This relation is the nonlinear characteristics of the exciter, given by:

U =                U A 1 - ∆P P M |∆P | P M sgn(∆P ) if ∆P < P M ; (4a) 0 if ∆P > P M . (4b)
The flow U A is calculated using the Bernoulli theorem:

U A = S 2P M ρ , ( 5
)
where S is the opening cross section of the reed channel at rest and ρ the density of the air.

Introducing the dimensionless variables and control parameters [START_REF] Chaigne | Instruments à anche[END_REF]:

∆p = ∆P /P M p = P/P M u = Z c U /P M γ = P m /P M ζ = Z c U A /P M . (6)
Z c = ρc/S res is the characteristic impedance of the cylindrical resonator of cross-section S res (c is the sound velocity). Equation (4) becomes:

F (p) =          ζ (1 -γ + p) |γ -p| sgn(γ -p) if γ -p < 1 ; (7a) 0 if γ -p > 1. (7b)
The parameters γ and ζ are the control parameters of the system. An example of the function F is shown in Fig. 2(a).

The resonator

Assuming that only plane waves exist in the resonator and propagate linearly, the resonator can be characterized by its reflection function r(t). The general expression relating p(t) to u(t) through r(t) is:

p(t) -u(t) = [r * (p + u)](t). ( 8
)
The resonator is modeled as a straight cylinder. Reflections at the open end of the resonator are considered perfect (no radiation losses) and viscous and thermal losses are ignored. In this case the reflection function becomes a simple delay with sign inversion:

r(t) = -δ(t -τ ), ( 9
)
where δ is the Dirac generalized function and τ = 2l/c is the round trip time of the sound wave with velocity c along the resonator of length l.

With the reflection function ( 9), equation ( 8) becomes:

p(t) -u(t) = -[p(t -τ ) + u(t -τ )] . (10)
Assuming that the blowing pressure γ skips instantaneously from 0 to a finite value and remains constant, p and u remain constant during the first half-period and hence during each forthcoming half-period. Therefore, p and u are square waves.

Using a discrete time formulation (the discretization is done at regular intervals τ ) and noting p(nτ ) = p n and u(nτ ) = u n , we obtain the following difference equation:

p n -u n = -(p n-1 + u n-1 ) . ( 11
)

Iterated map: outgoing and incoming wave representation

In linear acoustics any planar wave can be expanded into an outgoing wave p + and an incoming wave p -. Using the dimensionless variables defined in equation ( 6), the acoustic pressure p and flow u are given by:

p = p + + p - ; u = p + -p -, (12)
Replacing in equation [START_REF] Fruchard | Sur l'équation aux différences affine du premier ordre unidimensionnelle[END_REF],

p + = 1 2 (p + u) ; p -= 1 2 (p -u). (13)
By combining equations ( 7) and ( 12) a nonlinear relation G between p + and p -can be obtained:

p + = G -p -. (14)
An explicit expression of the function G was determined, for ζ < 1, by Taillard et al. [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF]. Fig. 2(b) shows an example of the function G. Using equations [START_REF] Fruchard | Bifurcation delay and difference equations[END_REF], the relation [START_REF] Fruchard | Sur l'équation aux différences affine du premier ordre unidimensionnelle[END_REF] becomes:

p - n = -p + n-1 . (15)
Finally, equations ( 14) and ( 15) define the iterated map [START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF][START_REF] Mcintyre | On the oscillations of musical instruments[END_REF]:

p + n = G -p - n = G p + n-1 . (16)
In the following, the variable p + will be used preferentially. The variable p can easily be calculated using equations ( 15) and (12).

Results from static bifurcation theory

The difference equation ( 16) can be analyzed using the static bifurcation theory, which assumes that the control parameters are constant. This will be hereafter referred to as the static case. The parameter γ will be specifically introduced as a subscript in the definition of the nonlinear characteristics [START_REF] Kapral | Bifurcation structure of the nonautonomous quadratic map[END_REF], stressing that this will be the parameter of interest in the current study (ζ will always consider to be constant):

p + n = G γ p + n-1 . (17)
Some of the predictions of the static bifurcation theory that are important to this work are recalled in the following sections while applying them to the map of equation ( 17) [START_REF] Kergomard | Elementary considerations on reedinstrument oscillations[END_REF][START_REF] Chaigne | Instruments à anche[END_REF][START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF].

Expression of the static regime and static oscillation threshold

For all values of the control parameter γ below a particular value of the parameter γ called static oscillation threshold and noted γ st the series p + n converges to a single value (the static regime), also referred to as the fixed point of G γ . It can be found by solving the following equation:

p + * = G γ p + * . ( 18
)
After solving the equation we obtain :

p + * (γ) = ζ 2 (1 -γ) √ γ. ( 19
)
When the static regime is reached p

+ n = p + n-1 = -p - n .
Therefore, for the variable p = p + + p -, the static regime is equal to zero.

The static regime exists for all values of the parameter γ but it is stable when γ < γ st and unstable when γ > γ st . The condition of stability of the static regime [START_REF] Kergomard | Elementary considerations on reedinstrument oscillations[END_REF] allowing to obtain the value of the static oscillation threshold is:

G γ p + * < 1, (20)
where G γ is the first derivative of the function G γ . The value of the static oscillation threshold is finally:

γ st = 1 3 . ( 21
)
Beyond the oscillation threshold, other bifurcations occur, the 2-valued oscillating regime becoming unstable and giving rise to a 4-valued oscillating state. This cascade is the classical scenario of successive period doublings, leading eventually to chaos [START_REF] Feigenbaum | The universal metric properties of nonlinear transformations[END_REF][START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF]. The values of the parameter γ for which appear the different 2n-valued oscillating regimes depend on the value of the parameter ζ: the smaller is ζ, the earlier the 2n-valued oscillating regimes appear. When γ = 1/2, whatever the value of ζ, a 2-valued oscillating regime reappears, the beating-reed regime . This is a particularity of model of the clarinet, it is due to the fact that when γp > 1 (equation (7b)) the reed presses against the mouthpiece lay. It can be shown [START_REF] Chaigne | Instruments à anche[END_REF] that in this permanent regime p = ±γ (c.f. Fig. 3).

Static bifurcation diagrams

Common representations of the static bifurcation diagram for clarinets usually show the steady state of the pressure inside the mouthpiece p or that of its amplitude (corresponding in the lossless model to the absolute value of p) with respect to the control parameter γ [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF]. In this paper, calculations are based on p + , so that most bifurcation diagrams will represent the steady state of the outgoing wave [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF].

Fig. 3 shows an example of these three representations of the static bifurcation diagram for ζ = 0.5. Fig. 3 represents only the two first branches of the diagrams. The first branch corresponds to the fixed points of the function G γ and the second branch represents the fixed points of the function (G γ • G γ ). On Fig. 3(b) and Fig. 3(c) the dashed line represents the curve of the static regime. For the variable p, the static regime is equal to zero and for the variable p + it is a function of the parameter γ, noted p + * (γ). Oscillating regimes with higher periodicities which may appear between γ = 1/3 and γ = 1/2 are not represented.

3 Time-varying blowing pressure 3.1 Problem statement

Definitions

Before presenting the problem, some definitions are introduced in order to avoid ambiguity in the vocabulary used hereafter. In the remainder of this paper, all simulations and calculations will be performed considering that the parameter ζ is a constant and equal to 0.5. The definitions presented below, used commonly in works dealing with bifurcation theory, can present some conflicts with that of musical acoustics. The terms that will be used in the remaining discussions are clarified in the following paragraphs:

Static case The control parameter γ is constant and the system is described by:

p + n = G γ p + n-1 . (22)
The steady state of the series p + n depends on the value of the control parameter γ. If γ is smaller than γ st , the series tends to a static regime. To avoid confusion, the static regime will now be called non-oscillating static regime. If γ is larger than γ st the steady state of p + n is an oscillating regime. This regime is called oscillating static regime. This behavior is still static, obtained by choosing a value of γ, letting the system reach its steady state, and repeating the procedure for each value of γ. Note that, even if the system tends to a steady state, the initial condition p + 0 often induces a transient regime.

Dynamic case As pointed in the introduction, in a musical context, the blowing pressure cannot always be considered constant. The dynamic case take this into account considering that the control parameter γ is variable and now written as γ n . When γ is a linear function of time, the system is described by the following difference equations:

p + n = G p + n-1 , γ n (23a) γ n = γ n-1 + . ( 23b
)
Since γ is changed only at each multiple of τ , the solution of equation ( 23) is still a square signal, i.e. two-state oscillating regime.

A slowly varying parameter implies that is arbitrarily small ( 1). The hypothesis of an arbitrarily small could be questioned in the context of the playing of a musical instrument. However, this hypothesis is required in order to use the framework of dynamic bifurcation theory (see forthcoming sections). An example of a numerical simulation performed on the system ( 23) is shown in Fig. 4 for ζ = 0.5, = 10 -3 and an initial condition γ 0 = 0.2. The initial value of the outgoing wave is p + 0 = G(-p - 0 = 0, γ 0 ). Indeed, for n = 0 the incoming wave p -is clearly zero, otherwise sound would have traveled back and forth with an infinite velocity.

The series p + n first shows a short oscillating transient, which will be called transient oscillating dynamic regime. This oscillation decays into a non-oscillating dynamic regime. Beyond a certain threshold, a new oscillation grows, giving rise to the final oscillating dynamic regime.

This paper will focus on the transition (i.e. the bifurcation) from the non-oscillating dynamic regime to the final oscillating dynamic regime. The value of the parameter γ for which the bifurcation occurs is called dynamic oscillation threshold, noted γ dt .

Bifurcation delay

Bifurcation delay occurs in nonlinear-systems with time varying control parameters. Fruchard and Schäfke [START_REF] Fruchard | Sur le retard à la bifurcation[END_REF] published an overview of the problem of bifurcation delay.

In fig. 5, the system (23) was simulated numerically, showing the time evolution of the series p + n and of the control parameter γ n (cf. Fig. 5(a)). To better understand the consequence of a time-varying parameter, the orbit of the series p + n is plotted as a function of the parameter γ n -in this case the evolution of the system can be interpreted as a dynamic bifurcation diagram. This is compared to the static bifurcation diagram in Fig. 5(b). We can observe that the static and the dynamic bifurcation diagrams coincide far from the static oscillation threshold γ st . However, in the dynamic case, we can see that the orbit continues to follow closely the branch of the fixed point of function G throughout a remarkable extent of its unstable range, i.e. after γ st : the bifurcation point is shifted from the static oscillation threshold γ st to the dynamic oscillation threshold γ dt . The term bifurcation delay is used to state the fact that the static oscillation threshold γ st is smaller than the dynamic oscillation threshold γ dt .

Non-standard analysis has been used in the past to study the phenomenon of bifurcation delay [START_REF] Fruchard | Sur l'équation aux différences affine du premier ordre unidimensionnelle[END_REF][START_REF] Fruchard | Bifurcation delay and difference equations[END_REF], explaing that one of the causes of the bifurcation delay is the exponential proximity between the orbit of the series p + n and the curve of the the fixed point of G. Other studies of bifurcation delay using standard mathematical tools -mathematics [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Baesens | Gevrey series and dynamic bifurcations for analytic slow-fast mappings[END_REF] or physics publications [START_REF] Kapral | Bifurcation structure of the nonautonomous quadratic map[END_REF][START_REF] Tredicce | Critical slowing down at a bifurcation[END_REF] -explain bifurcation delay as an accumulation of stability during the range of γ for which the fixed point of G is stable (i.e. 0 < γ n < γ st ). The dynamic oscillation threshold therefore appears as the value of the parameter γ at which the stability previously accumulated is compensated.

In musical acoustics literature some papers present results showing the phenomenon of bifurcation delay without never making a connection to the concept of dynamic bifurcation. For example this phenomenon is observed in simulations of clarinetlike systems using a slightly more sophisticated clarinet model (Raman's model) [START_REF] Atig | Saturation mechanism in clarinet-like instruments, the effect of the localised nonlinear losses[END_REF]. Raman's model takes losses into account although they are assumed to be independent of frequency (see [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF] for further explanation). Bifurcation delay can also explain the difficulty in estimating the static oscillation threshold by using a slowly variable blowing pressure [START_REF] Ferrand | Seuils d'oscillation de la clarinette : validité de la représentation excitateurrésonateur[END_REF]. In a preliminary work [START_REF] Bergeot | Measurement of attack transients in a clarinet driven by a ramp-like varying pressure[END_REF], bifurcation delays were experimentally observed in a clarinet-like instrument.

Analytical study of the dynamic case

This section presents an analytical description of a clarinetlike system in a dynamic case. The notion of invariant curve (φ(γ, )), invariant under the mapping [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF], will be needed for this study. The study of the stability of the invariant curve allows to define an analytical estimation of the dynamic oscillation threshold. A generic method to calculate the invariant curve is given by Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF]1, based on a perturbation method [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF].

Invariant curve

The invariant curve φ(γ, ) is invariant under the mapping ( 23), satisfying the following equation:

φ(γ, ) = G (φ(γ -, ), γ) . (24)
This curve plays a similar role for the dynamic system as fixed points for the static system, attracting or repelling the orbits. It is independent of the initial condition.

First of all, the invariant curve is expanded into a power series of , here truncated to the first order:

φ(γ, ) ≈ φ 0 (γ) + φ 1 (γ). (25)
Fig. 5 shows that, during the dynamic phase, the orbit of the series p + n closely follows the curve of the fixed points of G. This allows to linearize function G around the curve of the fixed points p + * (γ):

G(x, γ) ≈ G p + * (γ), γ + x -p + * (γ) ∂ x G p + * (γ), γ , (26) using the notation ∂ x G (x, y) = ∂G(x, y) ∂x , (27)
and knowing that G (p + * (γ), γ) = p + * (γ) (cf. equation ( 18)). Finally, using a Taylor expansion of φ(γ -, ) equation ( 24) is successively solved for the functions φ 0 (γ) and φ 1 (γ), yielding:

φ(γ, ) ≈ p + * (γ) + dp + * (γ) dγ ∂ x G (p + * (γ), γ) ∂ x G (p + * (γ), γ) -1 . (28)
Using the explicit expressions of p + * and dp + * /dγ we have:

φ(γ, ) ≈ ζ 2 (1 -γ) √ γ - ζ (3 γ -1) 4 √ γ ∂ x G (p + * (γ), γ) ∂ x G (p + * (γ), γ) -1 . (29)
More details about the calculation of the invariant curve are given in Appendix Appendix B.

To simplify the notation, in the rest of the document the invariant curve will be noted φ(γ). Its dependency on parameter is not explicitly stated.

Stability of the invariant curve and theoretical estimation of the dynamic oscillation threshold

A theoretical estimation of the dynamic oscillation threshold is done by identifying the value of γ for which the invariant curve looses its stability. The invariant curve is said to be unstable when the orbit of of the series p + n escapes from the neighborhood of the invariant curve φ(γ, ).

To investigate the stability of the invariant curve φ(γ, ), the function G in equation ( 23a) is expanded in a first-order Taylor series around the invariant curve [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF]:

p + n = G(p + n-1 , γ n ) ≈ G (φ(γ n -), γ n ) + p + n-1 -φ(γ n -) ∂ x G (φ(γ n -), γ n ) . (30)
A new variable is defined that describes the distance between the actual orbit and the invariant curve:

w n = p + n -φ(γ n ), (31)
and using equation [START_REF] Tredicce | Critical slowing down at a bifurcation[END_REF], equation (30) becomes:

w n = w n-1 ∂ x G (φ(γ n -), γ n ) . ( 32
)
The solution of equation ( 32) is formally:

w n = w 0 n i=1 ∂ x G (φ(γ i -), γ i ) , (33)
for n ≥ 1 and where w 0 is the initial value of w n . The absolute value of w n can be written as follow:

|w n | = |w 0 | exp n i=1 ln |∂ x G (φ(γ i -), γ i )| . (34)
Finally, using Euler's approximation the sum is replaced by an integral:

|w n | ≈ |w 0 | exp γn+ γ0+ ln |∂ x G (φ(γ -), γ )| dγ . ( 35
)
Equation ( 35) shows that the variable p + starts to diverge from the invariant curve φ(γ, ) when the argument of the exponential function changes from negative to positive. Therefore, the analytical estimation of the dynamic oscillation threshold γ th dt is defined by: where γ 0 is the initial value of γ. This result can be deduced from [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF] (equation (2.18)), it may also be obtained in the framework of non-standard analysis [START_REF] Fruchard | Canards et râteaux[END_REF].

γ th dt + γ0+ ln |∂ x G (φ(γ -), γ )| dγ = 0, (36)
The theoretical estimation γ th dt of the dynamic oscillation threshold depends on the initial condition γ 0 and the increase rate , it is therefore written γ th dt (γ 0 , ). A numerical solution γ th dt (γ 0 , ) of the implicit equation ( 36) is plotted in Fig. 6 as a function of the initial condition γ 0 and for = 10 -2 , 10 -3 and 10 -4 . γ th dt can be much larger than static oscillation threshold γ st = 1/3 for small initial conditions γ 0 . When the initial condition value γ 0 increases, γ th dt approaches the static threshold. Fig. 6(d) shows that the bifurcation delay seems to be independent of the increase rate if this value is sufficiently small (typically ≤ 10 -3 ). Equation (36) states that when γ = γ th dt we have |w n | ≈ |w 0 |, providing a good estimation of the dynamic oscillation threshold γ dt if |w 0 | is sufficiently small, i.e. if p + 0 is sufficiently close to φ(γ 0 ). γ 0 = 0 can be problematic since φ(0, ) = -∞, but a single iteration is sufficient to bring the orbit to a neighborhood of the invariant curve. Therefore, we make the assumption that

γ th dt (0, ) ≈ γ th dt ( , ). (37)
A non-exhaustive study done by running a few simulations shows that for = 10 -4 the error in γ dt due to this approximation is under 10 -8 , rising to 10 -7 when = 10 -3 and 2 × 10 -5 when = 10 -2 .

Benchmark of theoretical estimators for the dynamic threshold

Multiple criteria can be associated to the beginning of the oscillating regime. For instance, the oscillations can start before the series departs from the vicinity of the invariant curve as described in equation (36). Moreover, because of the approximation used between equations (35) and (36), the value of γ = γ th dt may not be an accurate estimation of the value at which the orbit departs from this vicinity.

For comparison, a dynamic oscillation threshold (noted γ num dt ) is calculated by simulating system (23) and compared with γ th dt . When the orbit of the series p + n is periodic, the sign of the second order difference of p + n changes sign at each iteration (i.e. the curve of p + n changes from upward to downward concave). In discrete time formulation the second order difference is given by:

δ 2 p + i = p + i -p + i-1 -p + i-1 -p + i-2 . (38)
Therefore, γ num dt , the oscillation threshold measured in numerical simulations, is reached when

δ 2 p + i-1 δ 2 p + i < 0, (39)
is satisfied. Then, in Fig. 6, γ num dt is compared with γ th dt (gray " * " depend on the increase rate, as previously shown in Fig. 6(d). Then, still in the case of a precision of 5000 decimals, γ num dt decreases for increasing , and γ th dt also decreases but to a lesser extent. For intermediate precisions (15, 100 and 500 decimals) the curve of γ num dt first increases before stabilizing close to the curve of γ th dt . For a given value of the precision, the larger the , the smaller is the accumulation of round-off errors created by the computer to reach a certain value of γ. This explains why the bifurcation delay first increases if the precision is not sufficiently high to simulate an analytic case. Beyond a certain value, all curves coincide with the one corresponding to the highest precision. That means that the system has reached the pair of parameters [precison ; ] needed to simulate an analytic case.

Conclusion

When considering mathematical models of musical instruments, oscillation threshold obtained through a static bifurcation analysis may be possibly very different from the threshold detected on a numerical simulation of this model.

For the first time for musical instruments, the differences between these two thresholds have been interpreted as the appearance of the phenomenon of bifurcation delay in connection with the concept of dynamic bifurcation.

Theoretical estimations of the dynamic bifurcation provided in this paper have to be compared with care to numerical simulations since the numerical precision used in computations plays a key role: for numerical precisions close to standard machine precision, the bifurcation towards the oscillating regime can occur at significantly lower mouth pressure values (while different most of the time from the threshold obtained through static bifurcation theory). Moreover, in that case, the threshold at which the oscillations start becomes more dependent on the increase rate of the mouth pressure.

The dependency on precision can be linked to the influence of noise generated by turbulence as the musician blows into the instrument. This would explain why the delays observed in artificially blown instruments are shorter than the predicted theoretical ones [START_REF] Bergeot | Measurement of attack transients in a clarinet driven by a ramp-like varying pressure[END_REF]. This will be the subject of further work on this subject, as well as the validity of these results for smoother curves of variation of the mouth pressure.

Moreover, in the light of results presented here for a basic model of wind instruments, varying the blowing pressure (even slowly) does not appear as the best way to experimentally determine Hopf bifurcations (static). In a musical context, since the blowing pressure varies through time, the dynamic threshold is likely to give more relevant informations than the static threshold, even if, in a real situation the influence of noise must be considered.

As a final remark, the simplistic model used in this work only describes one point per half-period of the sound played by the instrument. It is thus not suitable to describe different regimes (whose frequencies are harmonics of the fundamental one) that can be obtained by the instrument. However a simple extension of this model calculating the orbits of different instants within the half-period may be able to provide some insight on this subject.

Finally the expression of the invariant curve is:

Figure 1 :

 1 Figure 1: Schematic diagram of a single-reed mouthpiece. Presentation of variables, control parameters and choice of axis orientation. U is the flow created by the pressure imbalance Pm -P between the mouth and the bore, Ur is the flow created by the motion of the reed, Uin is the flow entering the instrument, y represents the position of the tip of the reed and H is the opening of the reed channel at rest.

Figure 2 :

 2 Figure 2: Nonlinear characteristics in u = F (p) representation (a) and p + = G(-p -) representation (b) for γ = 0.42 and ζ = 0.6.

Figure 4 :

 4 Figure 4: Numerical simulation performed on the system (23). (a) complete orbit of the series and (b) zoom near the non-oscillation dynamic regime. ζ = 0.5, = 10 -3 , γ0 = 0.2 and p + 0 = G(0, γ0).

Figure 5 :

 5 Figure 5: (a) Time evolution of the series p + n and of the control parameter γn.(b) Comparison between the series p + n and the static bifurcation diagram as a function of γn. ζ = 0.5, = 10 -4 , γ0 = 0 and p + 0 = G(0, γ0).

Figure 6 :

 6 Figure 6: Plot of γ dt a function of initial condition γ0, for different values of the slope . (a), (b) and (c): solid lines are the theoretical prediction γ th dt calculated from equation (36). Gray " * " markers represent the value γ num dt for which the system begins to oscillate. (d): combination of the previous theoretical predictions. " " represent the highest γ0 for which the system has enough time to reach a non-oscillating dynamic regime.

Figure 3 :

 3 Figure 3: Graphical representation of the static bifurcation diagrams for ζ = 0.5. Diagrams based on variables (a) |p|, (b) p and (c) p + . The dashed line represents the curve of the static regime, corresponding to the fixed point p + * (γ) of the function Gγ when the diagram is base on variable p + .

Figure 7 :

 7 Figure 7: Representation of the series p + n as a function of γn for ζ = 0.5, = 10 -3 , γ0 = 0.3 and p + 0 = G(0, γ0).

Figure 10 :

 10 Figure 10: Graphical representation of γ num dt as a function of for ζ = 0.5, γ0 = 0 and using five different precisions. A logarithmic scale is used in abscissa.

  φ(γ, ) ≈ p + * (γ) + ∂p + * (γ) ∂γ ∂ x G (p + * (γ), γ) ∂ x G (p + * (γ), γ) -1 .(52)

1In[START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF], the invariant curve is called adiabatic invariant manifold.

Acknowledgements

We wish to thank Mr. Jean Kergomard for his valuable comments on the manuscript. This work was done within the framework of the project SDNS-AIMV "Systèmes Dynamiques Non-Stationnaires -Application aux Instruments à Vent" financed by Agence Nationale de la Recherche (ANR).

markers). In some cases the series p + n never reaches the nonoscillating dynamic regime. An example of such situations is shown in Fig. 7. The values of γ num dt corresponding to the last initial values γ 0 for which the system has enough time to reach the non-oscillating dynamic regime are circled.

Fig. 6 shows that for = 10 -4 the theoretical result γ th dt provides a good estimation of the observed dynamic oscillation threshold. For = 10 -3 , the estimation is also good if the the initial condition is sufficiently small but as γ 0 gets closer to the static threshold γ st the system begins to oscillate before γ = γ th dt . Finally, for = 10 -2 , γ num dt is always smaller than γ th dt .

Limit of the model: influence of the precision

The phenomenon of bifurcation delay is very sensitive to noise: either numerical noise (round-off errors of the computer) or experimental noise (due to turbulence for instance). Indeed, when the static threshold is exceeded the system is very unstable. As a result, to observe bifurcation delay with numerical simulations and compare to theoretical results, it is necessary to perform calculations using a very high precision, as was done previously in this paper. For lower precisions the bifurcation delay can be considerably reduced (see [START_REF] Fruchard | Bifurcation delay and difference equations[END_REF] for an example in the logistic map).

Fig. 8 shows an example of numerical simulation performed on system [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF]. Fig. 8(a) and Fig. 8(b) differ only in the numerical precision (i.e. the number of decimal digits) used to calculate the orbit. The choice of the precision is possible using mpmath, the arbitrary precision library of Python. Fig. 8(a) was obtained using a precision of 5000 decimal digits, in this case γ th dt gives a good estimation of the bifurcation point, as it has already been shown in Fig. 6. On the other hand, using a precision of 15 decimal digits (Fig. 8(b)), the bifurcation delay is considerably reduced and the theoretical estimation of the dynamic oscillation threshold is not valid.

To highlight the influence of the precision γ num dt is calculated for different precisions. Results are plotted in Fig. 9 and compared to the analytical values γ st and γ th dt . The first thing to observe in Fig. 9 is the very high sensitivity of γ num dt to precision, yet all the values of γ num dt appear between γ st and γ th dt . For the lowest precision (7 decimal digits) the bifurcation delay disappears and γ num dt = γ st . If the precision is very high (typically ≥ 5000 decimals) γ num dt = γ th dt . Therefore, γ th dt can be interpreted as the limit of the bifurcation delay when precision tends to infinity. In cases with intermediate precisions (prec. = 15, 100 and 500) the bifurcation delay increases with the precision.

The sensitivity to the precision depends on the value of the increase rate : Fig. 10 plots γ num dt with respect to for different values of the numerical precision. Results are also compared with γ st and γ th dt . As above, for the lowest precision (7 decimals) the bifurcation delay disappears when is sufficiently small. Indeed, γ num dt is constant and equal to γ st . Then bifurcation delay occurs and increases with . The case of the highest precision (5000 decimals) is identical to an analytical case which would correspond to infinite precision. When is sufficiently small, the curves of γ num dt and γ th dt overlap. In this case when is small γ num dt is almost constant suggesting that the bifurcation delay does not

Appendix A Table of notation

Appendix B Invariant curve

The invariant curve φ(γ, ) is invariant under the mapping [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF], it therefore satisfies the following equation:

First of all, the invariant curve is expanded into a power series of and only he first-order is retained:

Secondly, the function G is linearized around the curve p + * (γ) of the fixed points:

where

Then, we make a Taylor expansion of φ(γ -, ): To obtain the approximate analytical expression of the invariant cure φ, equation (47) is successively solved for the functions φ 0 (γ) and φ 1 (γ).

As expected, to order 0 we find:

To order 1, we have to solve:


B Bergeot

A Almeida

C Vergez

B Gazengel

Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure: influence of noise

Keywords: Musical acoustics, Clarinet-like instruments, Iterated maps, Dynamic Bifurcation, Bifurcation delay, Transient processes, Noise, Finite precision

This paper presents an analysis of the effects of noise and precision on a simplified model of the clarinet driven by a variable control parameter.

When the control parameter is varied the clarinet model undergoes a dynamic bifurcation. A consequence of this is the phenomenon of bifurcation delay: the bifurcation point is shifted from the static oscillation threshold to an higher value called dynamic oscillation threshold.

In a previous work [8], the dynamic oscillation threshold is obtained analytically. In the present article, the sensitivity of the dynamic threshold on precision is analyzed as a stochastic variable introduced in the model. A new theoretical expression is given for the dynamic thresholds in presence of the stochastic variable, providing a fair prediction of the thresholds found in finite-precision simulations. These dynamic thresholds are found to depend on the increase rate and are independent on the initial value of the parameter, both in simulations and in theory.

Introduction

In classical (or static) bifurcation theory, all the parameters are constant, including the bifurcation parameter. The dynamic bifurcation theory focuses on systems where the bifurcation parameter is varying slowly over time. For a given system, the location of the bifurcation can be significantly different in the latter case.

A simple illustration is the flip bifurcation undergone by many one-dimensional discrete time nonlinear systems (among which the well known logistic map [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF] or a clarinet model [START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF][START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF]). When the bifurcation parameter is constant, the static bifurcation diagram summarizes the behavior of the system around the bifurcation : below the critical value of the parameter, the fixed point is stable (thus attractive), and above the critical value of the parameter the fixed point is unstable (thus repulsive) whereas a 2-valued cycle, born at the bifurcation, is stable (thus attractive). When the bifurcation parameter is varied over time, a bifurcation delay may appear : when the static bifurcation point is passed, the orbit remains in the neighborhood of the branch of the fixed points. After a certain time, the dynamic bifurcation dynamic oscillation threshold γ th dt theoretical estimation of the dynamic oscillation threshold of the clarinet model without noise or in "deterministic" situation γth dt theoretical estimation of the dynamic oscillation threshold of the noisy clarinet model in "sweepdominant" situation Γ th dt general theoretical estimation of the dynamic oscillation threshold of the noisy clarinet model, both for a "sweep-dominant regime" or a "deterministic regime" γ num dt dynamic oscillation threshold calculated on numerical simulations point is reached: the system escapes from the branch of the fixed points and moves abruptly to the 2-valued cycle. This behavior may be depicted in a dynamic bifurcation diagram. Fruchard and Schäfke [START_REF] Fruchard | Sur le retard à la bifurcation[END_REF] published an overview of the problem of bifurcation delay.

A previous article by the authors [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] analyzed the behavior of a simplified model of a clarinet when one of its control parameters (the blowing pressure) increases slowly linearly with time. Oscillations corresponding to the production of sound start at a much higher threshold than the one obtained in a static parameter case (i.e. higher than static bifurcation point of the system). The dynamic threshold (i.e. the dynamic bifurcation point) was described by an analytical expression, predicting that it does not depend on the increase rate of the blowing pressure (within the limits of the theory, i.e. slow enough increases), but that it is very sensitive to the starting value of the linear increase. This is a known behavior of such kind of system which timevarying parameter, shown by Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF] and, in the framework of nonstandard analysis, by Fruchard [START_REF] Fruchard | Canards et râteaux[END_REF].

These results are reproduced by simulations of the model, but only when very high precisions are used in the simulations. Running the simulation with the normal double-precision of a CPU results in much lower thresholds, although higher than the static ones.

The problem of the precision had already been mentioned in the seminal article [START_REF] Benoît | Chasse au canard ("duck hunting[END_REF] Chasse au canard (Duck Hunting in english).The canard phenomenon have similarities with the bifurcation delay. However, it can also appear in static situations: if the control parameter is higher than the static bifurcation point a stable limit cycle appears but, in particular cases, a delay can be observed in the limit cycle itself (see [START_REF] Benoît | Chasse au canard ("duck hunting[END_REF][START_REF] Fruchard | Sur le retard à la bifurcation[END_REF] for canards of forced Van der Pol equation). The shape of the resulting canard cycle in the phase space resembles that of a duck . This phenomenon can only exist in a very narrow interval of the parameter. Consequently, numerical simulations have to be performed with high precision and it was impossible in the beginning of the 80ies.

For the dynamic bifurcation, in contrast with the theory and simulations using high precision, when numerical simulations are running with finite precision, the dynamic threshold depends on the parameter increase rate, but doesn't depend on the starting value of the parameter. These properties have been observed on numerical simulations of the logistic map by Kapral and Mandel [START_REF] Kapral | Bifurcation structure of the nonautonomous quadratic map[END_REF] and in [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] in the case of a clarinet model.

To explain this discrepancy, round-off errors of the computer must be taken into account. In general this is done by introducing an ad hoc additive white noise in the model. For continuous-time systems we can cite Benoît [START_REF] Benoît | Linear dynamic bifurcation with noise[END_REF] and more recently Berglund and Gentz [START_REF] Berglund | The effect of additive noise on dynamical hysteresis[END_REF][START_REF] Berglund | Pathwise description of dynamic pitchfork bifurcations with additive noise[END_REF]. For discrete-time systems Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Baesens | Noise effect on dynamic bifurcations: The case of a period-doubling cascade[END_REF] propose a general method which is followed in the present paper.

Therefore, the aim of the present article is to formulate analytically an estimation of the dynamic bifurcation threshold in simulations performed with finite precision. The effect of finite numerical precision in simulations is modeled as an ad hoc additive white noise with uniform distribution. This hypothesis is tested in section 3. In section 4, a mathematical relation is derived for the behavior of the model affected with noise. The resulting theoretical expression of the dynamic oscillation threshold is compared to numerical simulations and its range of validity is discussed. The clarinet model and major results from [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] are first briefly recalled in section 2.

A table of notation is provided in table 1.

2 Dynamic oscillation threshold of the clarinet model without noise and problem statement

Clarinet model

The instrument is divided into two functional elements: the exciter and the resonator. The exciter of the clarinet is the reedmouthpiece system described by a nonlinear characteristics relating, by the Bernoulli equation, the instantaneous values of the flow u(t) across the reed entrance to the pressure difference ∆p(t) = p m (t)p(t) between the mouth of the musician and the clarinet mouthpiece [START_REF] Hirschberg | A quasi-stationary model of air flow in the reed channel of single-reed woodwind instruments[END_REF][START_REF] Hirschberg | Aero-acoustics of wind instruments[END_REF]. The reed is simplified into an ideal spring without damping or inertia. The resonator is approximated by a straight cylinder, described by its reflection function r(t). Considering that the resonator is a perfect cylinder in which the dispersion is ignored and the losses are assumed to be frequency independent [START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF][START_REF] Kergomard | Elementary considerations on reedinstrument oscillations[END_REF]. The reflection function r(t) becomes a simple delay with sign inversion (multiplied by an loss parameter λ) and is written:

r(t) = -λδ(t -τ ), (1)
where τ = 2L/c is the travel time for waves to propagate to the end of the resonator of length L at speed c and to return to the input.

The loss parameter λ takes into account the visco-thermal losses along the resonator, which at low frequencies are dominant over the radiation losses. It can be approximated by the expression:

λ = e -2αL , (2)
where α is the damping factor [START_REF] Keefe | Acoustical wave propagation in cylindrical ducts: transmission line parameter approximations for isothermal and non-isothermal boundary conditions[END_REF]:

α ≈ 3 • 10 -5 f /R. (3)
R is the resonator radius and f is chosen to be the fundamental playing frequency. A realistic value of the loss parameter is λ = 0.95.

The solutions p(t) and u(t) of the model depend on two control parameters: γ representing the blowing pressure and ζ related to the opening of the embouchure at rest. In this work, the control parameter ζ is always constant and equal to 0.5. Using the variables p + = 1 2 (p + u) and p -= 1 2 (pu) (outgoing and incoming pressure waves respectively) instead of the variables p and u, the nonlinear characteristic of the exciter is written:

p + = f -p -, γ . (4)
Outgoing and incoming pressure waves are also related through the reflection function r(t):

p -(t) = r * p + (t) = -λp + (t -τ ). (5)
Finally, by combining equations ( 4) and ( 5) and using a discrete time formulation (the discretization is done at regular intervals τ ) and noting p + (nτ ) = p + n and p -(nτ ) = p - n , we obtain the following iterated map [START_REF] Maganza | Bifurcations, period doublings and chaos in clarinet-like systems[END_REF][START_REF] Mcintyre | On the oscillations of musical instruments[END_REF][START_REF] Kergomard | Elementary considerations on reedinstrument oscillations[END_REF]:

p + n = f λp + n-1 , γ = G p + n-1 , γ , (6)
with, by definition: G(x) ≡ f (λx). The function G can be written explicitly for ζ < 1 (see Taillard et al. [START_REF] Taillard | Iterated maps for clarinet-like systems[END_REF]). When the control parameter γ is constant, for low values of γ the solution of eq. ( 6) stabilizes at an equilibrium point which corresponds to the fixed point p + * (γ) of the iterated function G. For a critical value γ st , namely the static bifurcation point (also called the static oscillation threshold) a flip bifurcation [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF] occurs, i.e.

G p

+ * (γ st ) = -1, (7)
leading to a 2-valued periodic regime that corresponds to sound production. For the lossless model (i.e. λ = 1) the static oscillation threshold is equal to γ st = 1/3. If λ < 1, γ st is larger than 1/3, an expression of γ st is given by Kergomard et al. [START_REF] Kergomard | Calculation of the spectrum of self-sustained oscillators using a variable troncation method[END_REF].

Dynamic bifurcation

For a linearly increasing control parameter γ, eq. ( 6) is replaced by eq. ( 8a) and (8b) :

p + n = G p + n-1 , γ n (8a) γ n = γ n-1 + . (8b)
The theory derived in section 4 requires that the parameter γ increases slowly, hence is considered arbitrarily small ( 1). Because of the time variation of the control parameter γ, the system (8) undergoes a bifurcation delay: the bifurcation point corresponding to the birth of the oscillations is shifted from the static oscillation threshold γ st [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF] to the dynamic oscillation threshold γ dt [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF]. The previous article by the authors [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] provides an analytical study of the dynamic flip bifurcation of the clarinet model (i.e. system (8)) in the case where λ = 1. The method is based on applications of dynamic bifurcation theory proposed by Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF]. The main focus of this work is on the properties of the dynamic oscillation threshold, recalled hereafter.

The trajectory of the system in the phase space (here constituted of a single variable p + ) through time is called the orbit. The dynamic oscillation threshold is defined as the value of γ for which the orbit escapes from the neighborhood of the invari-ant curve φ(γ, ). This definition is different from the one used in [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] where the dynamic threshold was defined as the value of γ for which the orbit starts to oscillate.

The invariant curve is the nonoscillating solution of the system [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF]. It plays the role of an attractor for variable parameters similarly to the role of the fixed point in a static case. The invariant curve is written as a function of the parameter, invariant under the mapping (8) and thus satisfying the following equation:

φ(γ, ) = G (φ(γ -, ), γ) . (9)
The procedure to obtain the theoretical estimation γ th dt of the dynamic oscillation threshold is as follows: a theoretical expression of the invariant curve is found for a particular (small) value of the increase rate (i.e.

1). The system ( 8) is then expanded into a first-order Taylor series around the invariant curve and the resulting linear system is solved analytically. Finally, γ th dt is derived from the analytic expression of the orbit. The dynamic oscillation threshold γ th dt is defined by [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF]:

γ th dt + γ0+ ln |∂ x G (φ(γ -), γ )| dγ = 0, (10)
where γ 0 is the initial value of γ (i.e. the starting value of the linear ramp). Two important remarks can be made on this expression (Fig. 6 of [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF]):

• γ th dt does not depend on the slope of the ramp , provided that is small enough,

• γ th dt depends on the initial value γ 0 of the ramp.

These properties are also observed in numerical simulations with very high precision.

Problem statement

The above theoretical predictions converge towards the observed simulation results for very high numerical precision (typically when thousands of digits are considered in the simulation). 1Figure 1(a) is a plot similar to figure 10 of [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF]. The only difference is that the bifurcation point γ num dt estimated on the simulation results is here defined by the point where the orbit leaves the neighborhood of the invariant curve. The motivation for this choice will appear clearly in section 3 where random variables are considered. [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] and noisy case [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF] with a noise of level σ = 10 -30 . The dynamic threshold of oscillation obtained over an average of 20 runs is plotted against the precision used in the simulations, showing that beyond a precision of about σ, the system affected with noise is insensitive to the precision. Moreover, figure 2 reveals that for a low numerical precision (though even significantly higher than typical precisions used in numerical simulations), the dependence of the bifurcation point on the initial value γ 0 is lost over a wide range of γ 0 .

The minimum precision for which round-off errors do not affect the behavior of the system depends on the precision itself and on the relative magnitude of the slope and the initial condition γ 0 . Indeed, figure 1 shows that, beyond a certain value of all curves join the one with highest precision. Curves for even higher precisions would overlap, allowing to conclude that they are representative of an infinitely precise case. As shown in figure 2, for given values of and of the numerical precision beyond a certain value of γ 0 , the theoretical result γ th dt allows to obtain a good prediction of the bifurcation delay.

As a conclusion, the theoretical results obtained in [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] are not able to predict the behavior of numerical simulations carried out at usual numerical precision. The aim of this paper is to show how the numerical precision can be included in a theoretical model that correctly describes numerical simulations. Firstly, it is shown that the model computed with a finite precision behaves similarly to the model with an ad-hoc additive white noise. This is done in the next section. Then, using theoretical results given by Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF], a modified expression describing the behavior of the model affected by noise (section 4) is proposed.

In figure 1(b), system (8) is simulated with λ = 0.95, a typical value to take into account losses in the cylindrical clarinet considered in this paper. The effect of the losses is to increase the dynamic threshold, as for the static one. However, the behavior of the lossless model and that of system with losses are qualitatively the same. Therefore, for sake of simplicity and without loss of generality, following analytical calculation and numerical simulations are performed using λ = 1.

Finite precision versus additive white noise

Differences between theoretical predictions and numerical simulations highlighted in the previous section are due to round-off errors that accumulate for finite precisions. The aim of this section is to test whether round-off errors of the computer can be modeled as an additive independent and identically distributed random variable (referred to as an additive white noise). This result is used in section 4 to derive theoretical predictions of the dynamic bifurcation point γ dt .

Results

Two numerical results are compared. The first is the simulation of the system (8) using a numerical precision pr 1 (hereafter referred as a finite precision case). The second one (hereafter referred as a noisy case) is the simulation of the following stochastic system of difference equations:

p + n = G p + n-1 , γ n + ξ n (11a) γ n = γ n-1 + , ( 11b
)
where ξ n is a uniformly distributed stochastic variable with an expected value equal to zero (i.e. E [ξ n ] = 0) and a level σ defined by:

E [ξ m ξ n ] = σ 2 δ mn , (12)
where δ mn is the Kronecker delta. The definition of the expected value E is provided in [START_REF] Ross | Introduction to Probability Models, 9ème édition, chap. 2 "Random variables[END_REF]. For comparison with the finite precision case the noise level σ is equal to 10 -pr1 . The bifurcation point γ num dt estimated on the simulations is defined as the value of γ for which the orbit leaves the neighborhood of the invariant curve. Since the mean value of the white noise ξ n is zero, the relevant quantity to study is the mean square deviation of the orbit from the invariant curve. Therefore, γ num dt is reached when:

E [w 2 n ] = , (13)
where w n = p + nφ(γ n , ) describes the distance between the actual orbit and the invariant curve. Among other possible criteria, the condition ( 13) is chosen because it is also used in the analytical calculation made in section 4.

To simplify the notation, in the rest of the document the invariant curve will be noted φ(γ). Its dependency on parameter is no longer explicitly stated.

In figures 3 and 4, γ num dt is estimated in the finite precision case and in the noisy case. In both cases an average is made on w 2 n obtained in 20 different simulations. Then, γ num dt is calculated on the mean signal using equation [START_REF] Fruchard | Sur le retard à la bifurcation[END_REF].

In figure 3, γ num dt is plotted with respect to the numerical precision for which both systems [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] and [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF] are simulated. The noise level σ of the noisy case modeled by the system (11) is equal to 10 -30 . For numerical precision below log 10 (σ) = 30, the noise level is smaller than round-off errors of the computer. In these situations, the effect of the additive noise in system [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF] is hidden by the effect of the round-off errors of the computer. The consequence is that the thresholds computed in finite precision case and in noisy case are equals.

For numerical precisions higher than 30, γ num dt computed on system [START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF] is constant because the influence of the round-off errors is now hidden by the additive noise which have a fixed level. Figure 3 shows that the transition between the regime for which the round-off error effect prevails over the additive noise affect and the regime for which the opposite occurs is abrupt. Therefore, the region where mixed effects of both round-off errors and additive noise play a role is very narrow. However, to avoid any influence of the numerical precision, the system (11) is simulated using a precision pr 2 = 2pr 1 .

Figure 4 confirms that the kind of noise introduced in the stochastic system can correctly describe the influence of a finite precision. Indeed, with the exception of the smallest precision (pr 1 = 7), the curves are nearly superimposed. Hence, in the next section, the stochastic system ( 11) is studied theoretically in order to predict results of numerical simulations of system (8) with finite precision.

Relevance of numerical results

To investigate the relevance of the numerical results, several indicators are calculated. First, the standard deviation STD w 2 n of the signal w 2

n is calculated at the dynamic threshold γ num dt and compared to E w 2 n , also calculated at γ num dt . Secondly, the dynamic threshold is calculated on each run. We obtain 20 values, noted γ num,i dt (i ∈ [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Kergomard | Calculation of the spectrum of self-sustained oscillators using a variable troncation method[END_REF]). The mean value E γ num,i dt is compared to the value γ num dt , estimated on the mean signal E [w 2 n ] (see section 3.1 where this numerical estimation method of γ num dt is used because it also used in analytical calculations in section 4.). The standard deviation STD γ num,i dt is calculated too.

Results are presented in table 2. The mean value E w 2 n and the standard deviation STD w 2 n at the dynamic threshold have the same order of magnitude. This suggests a bad repeatability of the numerical simulations. However, at the dynamic threshold, w n diverges sharply and a large deviation of it does not necessarily imply a large deviation of the dynamic threshold. The standard deviation STD γ num,i dt , in table 2, shows precisely a good repeatability of γ num,i dt .

Analytical study of the noisy dynamic case 4.1 General solution of the stochastic clarinet model

This section introduces a formal solution of the stochastic model that is valid when the orbit is close to the invariant curve. Function G in equation ( 11a) is expanded into a first-order Taylor series around the invariant curve. Using the variable w n = p + nφ(γ n ), the system (11) becomes:

w n = w n-1 ∂ x G (φ(γ n -), γ n ) + ξ n (14a) γ n = γ n-1 + . ( 14b
)
The solution of equation ( 14a) is [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF]:

10 -3
w n = w 0 n k=1 ∂ x G (φ(γ k -), γ k ) + n k=1 ξ k n m=k+1 ∂ x G (φ(γ m -), γ m ) , ( 15
)
where w 0 is the initial value of w n .

Because the additive white noise ξ n has a zero-value mean, as in section 3, the relevant indicator is the mean square deviation of the orbit from the invariant curve: E [w 2 n ]. Equation ( 15) squared becomes:

w 2 n = w 0 n k=1 ∂ x G (φ(γ k -), γ k ) 2 + n k=1 ξ k n m=k+1 ∂ x G (φ(γ m -), γ m ) 2 + 2w 0 n k=1   n j=1 ∂ x G (φ(γ j -), γ j )   ξ k × n m=k+1 ∂ x G (φ(γ m -), γ m ) . ( 16
)
Averaging has no effect on the first term of the right-hand side of equation ( 16) because the stochastic variable ξ n is not present. Using eq. ( 12), the average of the second term is simplified to:

σ 2 n k=1 n m=k+1 ∂ x G (φ(γ m -), γ m ) 2 . ( 17
)
Because E [ξ n ] = 0, the average of the third term of the righthand side of equation ( 16) is also equal to zero. Using the fact that a product can be expressed as an exponential of a sum of logarithms, the final expression of E w 2 n is given by:

E w 2 n = w 2 0 exp n k=1 ln |∂ x G (φ(γ k -), γ k )| 2 An + σ 2 n k=1 exp n m=k+1 ln |∂ x G (φ(γ m -), γ m )| 2 Bn . ( 18
)
The two terms of the right-hand side of equation ( 18) are denoted A n and B n .

Finally, using Euler's approximation, sums are replaced by integrals and the expressions of A n and B n become:

A n ≈ w 2 0 exp γn+ γ0+ 2 ln |∂ x G (φ(γ -), γ )| dγ , (19)
B n ≈ σ 2 γn+ γ0+ exp γn+ γ+ 2 ln |∂ x G (φ(γ -), γ )| dγ dγ. (20)
A n corresponds to the precise case studied in [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] which leads to the theoretical estimation γ th dt of the dynamic oscillation threshold for the system without noise (cf. equation [START_REF] Berglund | Pathwise description of dynamic pitchfork bifurcations with additive noise[END_REF]). B n is the contribution due to the noise.

The transform from discrete sums to integral can be questioned. Indeed, to transform the term A n in equation [START_REF] Keefe | Acoustical wave propagation in cylindrical ducts: transmission line parameter approximations for isothermal and non-isothermal boundary conditions[END_REF] to its integral form [START_REF] Kergomard | Elementary considerations on reedinstrument oscillations[END_REF], we assume that:

n k=1 ln |∂ x G (φ(γ k -), γ k )| ≈ γn+ γ0+ ln |∂ x G (φ(γ -), γ )| dγ , (21)
but if ∂ x G (φ(γ i -), γ i ) crosses over zero, one term in the sum comes close to ln(0) = -∞ and the equality ( 21) is respected only for small enough values of . Otherwise, in equation [START_REF] Keefe | Acoustical wave propagation in cylindrical ducts: transmission line parameter approximations for isothermal and non-isothermal boundary conditions[END_REF], we have: A n → 0 and B n → σ 2 , and consequently E [w 2 n ] → σ. In this case, for small noise levels, the difference between the orbit and the invariant curve comes close to zero, and as a result, the orbit needs more time to escape from the invariant curve neighborhood, i.e. the bifurcation delay is lengthened. This phenomenon is mentioned by Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Baesens | Gevrey series and dynamic bifurcations for analytic slow-fast mappings[END_REF] and Fruchard [START_REF] Fruchard | Sur le retard à la bifurcation[END_REF]. It can be observed for example in figure 1 where peaks (i.e. larger bifurcation delay) are there on curves in the above right part of the figure. In the rest of the paper, we use integral form, because it allows analytical integrations of noise contribution B n which will be considered in the remaining of this section.

A first glance on equations ( 19) and ( 20) allows to explain observation made in Section 2.3. Indeed, comparing the expressions of A n and B n , it possible to distinguish [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Baesens | Noise effect on dynamic bifurcations: The case of a period-doubling cascade[END_REF] two operating regimes, which, for a given value of w 0 , depends on , σ and γ 0 :

• A n B n (deterministic regime)
In this case the noise does not affect the bifurcation delay and the dynamic oscillation threshold can be determined by eq. ( 10).

• A n B n (sweep-dominant regime)
In this case, the bifurcation delay is affected by the noise. This regime is studied in the following section.

In Section 2.3, figures 1 and 2 represent two different cases distinguished by the parameter values: in certain areas of the figures, the dynamic bifurcation threshold does not depend on but depends on γ 0 , while in other areas the dynamic bifurcation threshold depends on but is not dependent on γ 0 . This observation may be interpreted as the existence of the two regimes detailed above: a sweep-dominant regime and a deterministic regime. The transition between the two regimes occurs abruptly as observed in figures 1 and 2.

Theoretical expression of the dynamic oscillation threshold of the stochastic model

The next step is to find an approximate expression of the standard deviation E [w 2 n ] for the sweep-dominant regime. In this regime, the term A n is negligible with respect to the contribution B n due to the noise, i.e. E [w 2 n ] ≈ √ B n . The purpose is to obtain a statistical prediction of the dynamic oscillation threshold for the stochastic system, hereafter referred as γth dt . It is assumed that 1, which implies that the invariant curve φ(γ) and the curve p + * (γ) of the fixed points in eq. ( 6) are close [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF], and allows the approximation:

∂ x G (φ(γ -), γ) ≈ ∂ x G p + * (γ), γ . (22)
Moreover, because of the noise, the bifurcation delay is expected to occur earlier, so that the dynamic oscillation threshold γ dt is assumed to be close2 to the static oscillation threshold γ st . The term ∂ x G (p + * (γ), γ) is then expanded in a first-order Taylor series around the static oscillation threshold γ st :

∂ x G p + * (γ), γ ≈ ∂ x G p + * (γ st ), γ st -1 : flip bifurcation + (γ -γ st ) ∂ xy G p + * (γ st ), γ st noted -K , (23)
finally we have:

∂ x G p + * (γ), γ ≈ -1 -K (γ -γ st ) , (24)
which is used in equation [START_REF] Kergomard | Calculation of the spectrum of self-sustained oscillators using a variable troncation method[END_REF]. Figure 5 shows the comparison between ∂ x G (p + * (γ),γ) and its tangent function -1 -2This hypothesis could be questioned because according to figures 1 and 2, even in the presence of noise, the bifurcation delay can be large. However, this hypothesis is required to carry out following calculations. Using expression (24) the integral

I 1 = γn+ γ+ 2 ln |∂ x G (φ(γ -), γ )| dγ , (25)
contained in the expression (20) of B n becomes:

I 1 = 2K γn+ γ+ (γ -γ st ) dγ = K (γ -γ st ) 2 γn+ γ+ .
(26) The small correction in the domain of integration can be neglected since 1. Therefore, we obtain:

I 1 = K (γ -γ st ) 2 γn γ = K (γ n -γ st ) 2 -(γ -γ st ) 2 . ( 27
)
By combining equations ( 20) and ( 27), B n is now written as:

B n ≈ σ 2 γn+ γ0+ exp K (γ n -γ st ) 2 -(γ -γ st ) 2 dγ = σ 2 exp K (γ n -γ st ) 2 × γn+ γ0+ exp - K (γ -γ st ) 2 dγ I2 . (28)
The function which appears in the integral I 2 is a Gaussian function with standard deviation

ν = 2K . ( 29
)
Integral I 2 is then [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]:

I 2 = 1 2 π K erf K (γ -γ st ) γn γ0 , (30)
where erf(x) is the error function. The initial condition γ 0 is supposed to be much lower than the static threshold γ st , so that equation (30) can be written:

I 2 = 1 2 π K erf K (γ n -γ st ) + 1 . (31)
The dependence on the initial condition γ 0 is now lost.

Since 1, for γ n > γ st the error function quickly becomes equal to 1 and the integral I 2 is simplified to I 2 = π K . Finally the expression of B n is:

B n ≈ σ 2 √ π K exp K (γ n -γ st ) 2 . ( 32
)
From equation (32) it is possible to obtain the expression of

E [w 2 n ] ≈ √ B n : E [w 2 n ] ≈ σ -1/4 π K 1/4 exp K 2 (γ n -γ st ) 2 . ( 33
)
The dynamic oscillation threshold γth dt is defined [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Baesens | Noise effect on dynamic bifurcations: The case of a period-doubling cascade[END_REF] as the value of γ n for which the standard deviation E [w 2 n ] leaves the neighborhood of the invariant curve. More precisely, the bifurcation occurs when E [w 2 n ] becomes equal to the increase rate , as defined in eq. ( 13). Finally, using equation (33), we have:

γth dt = γ st + - 2 K ln π K 1/4 σ 5/4 , (34)
which is the theoretical estimation of the dynamic oscillation threshold of the stochastic systems (11) (or of the system (8) computed using a finite precision) when it evolves in a sweepdominant regime. The bifurcation delay is a by-product of eq. (34) since it is simply γth dtγ st . The method presented in this section is based on a first-order Taylor expansion of the system (11) around the invariant curve φ(γ n ), leading to the linear system [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]. Using an asymptotic expansion of the error function it is possible to investigate the behavior of √ B n before γ n enters the neighborhood of the static oscillation threshold γ st . This study allows to define the domain of validity of this linear approximation, as done by Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Baesens | Noise effect on dynamic bifurcations: The case of a period-doubling cascade[END_REF]. This is σ √ (more details on obtaining the domain of validity are given in Appendix A). Otherwise, if σ √ , the orbit of the series p + n leaves the neighborhood of invariant curve φ(γ) before the static oscillation threshold is reached. In this case, the linear approximation is no longer valid. This situation is called by Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Baesens | Noise effect on dynamic bifurcations: The case of a period-doubling cascade[END_REF] noise-dominant regime and it is not investigated in the present paper. However, figure 6 shows the domain of existence of the different regimes in a plan [ ;log 10 (σ)]. The frontier between deterministic and sweep- However, when slowly increasing one of the control parameters, the distances between the state of the system and the invariant curve can become smaller than the round-off errors of the calculation, with dramatic effects on the time required to trigger an oscillation. In these cases, the inclusion of a stochastic variable in the theory allows to correctly estimate the threshold observed in simulations, which lies between the static and dynamic thresholds found for precise cases.

As a final remark, the present theoretical study is probably not restricted to describe numerical simulations. Indeed, the noise level σ measured in an artificially blown instrument is typically of the order of magnitude of 10 -3 . The domain of validity of the results: σ √ suggests that the comparison with experiment using blowing pressure with increase rates > 10 -6 (typically for usual clarinets that corresponds to ≈ 5Pa/s), seems to be possible although the noise level usually increases with the pressure applied to the instrument.

It is known that the static oscillation of the clarinet is difficult to measure by increasing, even slowly, the mouth pressure. The phenomenon of dynamic bifurcation is a possible reason. If that were proven experimentally, we could imagine to inverse the equation (34) to deduce the static threshold from the measurement of the noise level, the increase rate of the blowing pressure and the dynamic threshold.

B Asymptotic expansion of error function

The asymptotic expansion of the error function erf(x) for large negative x (x → -∞) is [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF]

Figure 1 :

 1 Figure 1: Graphical representation of γ num dt for different precisions (prec. = 7, 15, 30, 100, 500 and 5000) with respect to the slope and for γ0 = 0. Results are also compared to analytical static and dynamic thresholds: γst and γ th dt . (a) lossless model: λ = 1 and (b) typical losses in a cylindrical clarinet, λ = 0.95.

  ( =1 •10 -4 ) γ num dt ( =3 •10 -4 ) (a)Numerical precision is fixed (prec. = 50). γ num dt is computed for = 10 -4 and 3 • 10 -4 . 0.00 0.05 0.10 0.15 0.20 0.25 0.30 (prec. =15) γ num dt (prec. =100) (b) The increase rate of γ is fixed ( = 3 • 10 -4 ). γ num dt is computed for numerical precisions equal to 15 and 100.

Figure 2 :

 2 Figure 2: Plot of γ dt as a function of the initial condition γ0. Solid lines are the theoretical prediction γ th dt calculated from equation (10). Dashed line represent the values γ num dt .

  Fig-ure1 1(a) shows that for the usual double-precision of CPUs (around 15 decimals), theoretical predictions of the dynamic bifurcation point γ th dt are far from thresholds estimated on the numerical simulation results γ num dt . In particular, the numerical bifurcation point γ num dt depends on the slope , in contrast with the theoretical predictions γ th dt .

Figure 3 :

 3 Figure 3: Comparison of the dynamic threshold γ num dt obtained in numerical simulations of a clarinet model in finite precision case[START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF] and noisy case[START_REF] Dalmont | An analytical prediction of the oscillation and extinction thresholds of a clarinet[END_REF] with a noise of level σ = 10 -30 . The dynamic threshold of oscillation obtained over an average of 20 runs is plotted against the precision used in the simulations, showing that beyond a precision of about σ, the system affected with noise is insensitive to the precision.

Figure 4 :

 4 Figure 4: Comparison between γ num dt computed for finite precision cases and for noisy cases. For both cases and for each value of we compute the average of the signals wn = p + n -φ(γn) obtained over 20 runs. Then, γ num dt is calculated on the resulting signal. The numerical precisions used to simulate the finite precision cases are 7, 15, 30 and 100 decimal digits. γ0 = 0.

Figure 6 :

 6 Figure 6: Domain of existence of the deterministic, sweep-dominant and noise-dominant regimes in a plane ( , -log 10 (σ)). For finite precision cases, -log 10 (σ) corresponds to the value of the precision.

8 Figure 5 :

 85 Figure 5: Graphical representation of ∂xG p + * (γ),γ and its tangent function -1 -K (γ -γst) around the static oscillation threshold for ζ = 0.2, 0.5 and 0.8.

Figure 7 :

 7 Figure 7: Graphical representation of γ num dt for different precisions (prec. = 7, 15, 30, 100, 500 and 5000) with respect to the slope and for γ0 = 0. Results are also compared to analytical static and dynamic thresholds: γst, γ th dt and γth dt . (a) γ num dt and only finite precision cases are represented. (b) Both finite precision cases and noisy cases are represented for prec. = 7 and 15.

Figure 9 :

 9 Figure 9: Relative errors: RE [γst], RE γ th dt and RE Γ th dt for numerical precisions equal to 15 (a) and 100 (b).

Table 1 :

 1 Table of notation. All quantities are dimensionless.

		Table of Notation
	G	iterative function
	γ	musician mouth pressure (control parameter)
	ζ	control parameter related to the opening of the em-
		bouchure at rest
	p +	outgoing wave
	p -	incoming wave
	p + *	non-oscillating static regime of p + (fixed points of
		the function G)
	φ	invariant curve
	w	difference between p + and φ
		increase rate of the parameter γ
	σ	level of the white noise
	γst	static oscillation threshold
	γ dt

Table 2 :

 2 Mean value E w 2 n and standard deviation STD w 2 n of the signal w 2 n , calculated at the dynamic threshold γ num[START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Kergomard | Calculation of the spectrum of self-sustained oscillators using a variable troncation method[END_REF]) calculated on each run. All results are calculated for σ = 10 -7 and 10 -15 and for = 10 -4 , 10 -3 and 10 -3 . 10 -8 1.18 • 10 -6 2.83 • 10 -4 1.20 • 10 -8 1.46 • 10 -6 5.46 • 10 -4 10 -8 1.56 • 10 -6 3.02 • 10 -4 2.07 • 10 -8 1.86 • 10 -6 1.25• 10 -4

	dt	. Mean values E [w 2 n ]

Acknowledgements

This work was done within the framework of the project SDNS-AIMV "Systèmes Dynamiques Non-Stationnaires -Application aux Instruments à Vent" financed by Agence Nationale de la Recherche (ANR).

dominant regime corresponding to A n = B n is determined numerically using the equality γ th dt = γth dt . The condition σ √ is respected in this work since σ = 10 -pr with 7 ≤ pr ≤ 5000 and 8.10 -5 ≤ ≤ 10 -2 .

Discussion

In figure 7, γth dt defined by equation ( 34) is plotted against the increase rate . It is compared with γ num dt for different values of the precision and for γ 0 = 0. In figure 7(a), γ num dt is represented for finite precision cases. The differences between finite precision cases and stochastic cases observed for prec. = 7 and 15 are shown in figure 7(b). The theoretical result γth dt provides a good estimation of the dynamic oscillation threshold as long as the system remains in the sweep-dominant regime (with a better estimation when the bifurcation delay is small3). Otherwise, γ th dt is a better approximation of γ num dt , as expected in the deterministic regime.

Figure 8 shows the comparison between γth dt and γ num dt (only for finite precision cases) plotted against the initial condition γ 0 . In figure 8(a), variables are plotted for several values of and for a fixed numerical precision. The opposite is done in figure 8(b). As in figure 7, γth dt provides a good estimation of the dynamic oscillation threshold in the sweep-dominant regime, as well as γ th dt in the deterministic regime. Finally, to predict theoretically the dynamic bifurcation threshold Γ th dt of the stochastic system (11) (as well as of the system (8) when it is computed with a finite precision) the following procedure is proposed:

• compute the theoretical estimation γth dt of the stochastic system through eq. (34)

3This is an expected result because of the initial assumption of a small bifurcation delay in the presence of noise, leading to first-order Taylor expansions γst in previous calculation (see equation ( 23)).

• compute the theoretical estimation γ th dt of the system without noise through eq. ( 10)

• if γth dt < γ th dt the system remains in the "sweep-dominant regime" and the dynamic threshold Γ th dt is equal to γth dt , otherwise the "deterministic regime" is attained and the dynamic threshold Γ th dt is equal to γ th dt . Figure 9 compares the relative error RE of the three theoretical predictions of the oscillation threshold (γ st , γ th dt and Γ th dt ) with respect to γ num dt , as a percent value:

where X takes successively the values of γ st , γ th dt and Γ th dt . For standard double-precision (fig. 9(a), prec.=15), the sweepdominant regime is prevalent throughout most of the range of increase-rates studied in this article. Higher precisions (for instance prec.=100) imply the appearence of the deterministicregime for lower increase-rates. In this case, Γ th dt provides a better estimation of the oscillation threshold of the clarinet with a linearly increasing blowing pressure. Indeed, in situations represented in figure 9, RE Γ th dt never exceeds 15% while RE [γ st ] and RE γ th dt can reach 60% and 145% respectively. At slightly lower values of than the limit between the two regimes, γ th dt still provides a better estimation of γ num dt than Γ th dt , a situation that occurs for all values of the precision, according to figure 7.

Conclusion

In many situations, the finite precision used in numerical simulations of the clarinet system does not produce major errors in the final results that are sought. Such is the case, for instance, when estimating the amplitudes for a given regime.

A Limit of the linear calculation

The method presented in Section 4 is based on a first-order Taylor expansion of the system (11) around the invariant curve φ(γ n ) leading to define the linear system [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]. Following Baesens [START_REF] Baesens | Slow sweep through a period-doubling cascade: Delayed bifurcations and renormalisation[END_REF][START_REF] Baesens | Noise effect on dynamic bifurcations: The case of a period-doubling cascade[END_REF], we give here the upper bound of the domain of validity of this linear approximation.

Using equations ( 20) and (31), the expression B n is given by:

We investigate the behavior of E w 2 n before γ n enters in the neighborhood of the static oscillation threshold γ st . More precisely, we compute an approximate expression of E w 2 n when γ n < γ stν, where ν is defined by equation (29). To do this, the error function in equation ( 36) is expanded in a firstorder asymptotic series [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF] (the asymptotic expansion of the error function erf(x) for large negative x is recalled in Appendix B):

which is simplified in:

Using the explicit form of γ n , solution of equation (8b):

and (38), we have:

where n st is the iteration for which γ st is reached. Equation (40) means that when γ n < γ stν, the standard deviation E [w