Direct and inverse Hopf bifurcation in a neutral delay differential equation
model of reed conical instrument
Tom Colinot(1) , Louis Guillot(1) , Jean Kergomard(1)(1) Aix-Marseille Univ, ECM, CNRS, LMA, UMR7031, Marseille, France, kergomard@lma.cnrs-mrs.fr

Abstract
In conical musical instrument, the self-sustained oscillations appear at a blowing pressure threshold. At thisthreshold, the mathematical model encounters a so-called Hopf bifurcation. The bifurcation can be either direct
or inverse which determines the stability of the arising oscillation. It is well known that for cylindrical instruments (clarinet-like instruments), the bifurcation is direct most of the time. In the case of conical instruments,
the nature of the bifurcation has been studied analytically before, for instance using geometrical approximations
based on cylinders. It is expected to ﬁnd an inverse bifurcation in some cases. The proposed study addresses thebifurcation of an idealized model of reed instrument with a lossless conical resonator and a compliant mouthpiece, using two methods : a direct time-integration algorithm and a Taylor series based continuation methodof periodic solutions, based on the harmonic balance method. It allows to draw the whole branch of periodic
solutions and to deduce the type of Hopf bifurcation from a full diagram. Super-critical Hopf bifurcations canbe found, with two fold bifurcations whose positions depend on the value of the geometrical parameters of theinstrument.
Keywords: Hopf bifurcation, conical instrument, neutral delayed differential equations, harmonic balancemethod, time-integration

1INTRODUCTION

The study of periodic solutions of physical models has shown to be an efﬁcient tool to understand the behaviorof the systems considered. A ﬁrst common application is to predict and avoid undesired instabilities. In musicalacoustics, the periodic solutions are sought for as they correspond to the regimes where an instrument produces
a note. A lot of methods are available in the literature to study these regimes. Among those, two essentiallydifferent approaches will be considered. The ﬁrst one is time-integration method. From the knowledge of thesystem before the starting time of integration, it reconstructs the behavior of the system for all greater times.
Without a priori knowledge on the solution type, it is then possible to reconstruct the whole signal or waveform.
For all the values of the parameters, a new time-integration is performed and thus the behavior of the system
is known on a range of the parameters. However, this can be more complex in the case of several co-existing
solutions and these methods are often unable to obtain unstable regimes. If the aim is the study of periodic
solutions, the transient regime before the steady state can be very long and difﬁcult to characterize. The second
type of approach is continuation methods. These methods allow to keep track of the waveform as one parameterof the system is varied. Only steady-state solutions are considered. Thus, from one starting periodic solution
a whole branch of periodic solution is determined. From the solution branch, it is possible to determine thecritical value of the parameter above which a sound arise (the Hopf bifurcation point). A study of the oscillationthresholds of the clarinet by continuation is proposed in [8].
In this short article, the two complementary approaches are applied to a simpliﬁed model of saxophone. Previous work treat the bifurcation of saxophone-like instruments using multi-cylindrical approximations of theresonator [5, 6] or approaches based on the input impedance [7]. The model we use here has a conical resonator modelled by a delayed system. It is taken from [3]. Although losses and reed dynamics are ignored,
this model has been shown to reproduce internal pressure waveform for the saxophone. A reformulation of themodel leads to a strongly nonlinear neutral delay differential equation, for which no analytical solution exists.

6446

However, the equation can be solved numerically using sound synthesis (time-integration) or harmonic balance.
From sound synthesis, the general features of the produced sound are determined in a two-dimensional spaceof parameters value. This study helps to give a hint of the behavior of the system but fails to state clearly
whether the Hopf bifurcation is super-critical of sub-critical. In practice, this information is interesting from themusician point of view. The super-critical case features a soft appearance of the sound in the instrument while
the sub-critical case is often linked with a sudden appearance of the sound and a hysteresis loop. The continuation study of the same model is performed with a Taylor series based method coupled with the harmonicbalance method [2]. It allows to have a continuous representation of the solution-branch. From the continuation
diagram, the nature (super- or sub-critical) of the Hopf bifurcation and the hysteresis loops can be determined.
The results obtained with the two methods are compared to validate the approaches.

2

Model and numerical methods

2.1 Model
The acoustical variables at the input of the resonator p and u are in dimensionless form, obtained from the thephysical values p̂ and û asûp̂p=, u = Zc
,
(1)
pM
pM
where pM is the static pressure necessary to close the reed channel completely, and Zc = ρc
S is the characteristicimpedance. The control parameters of the model are2pm, ζ = Zc wH,
(2)
γ=pM
ρ pM
where pm is the pressure in the mouth of the musician, w is the width of the reed channel and H is the heightof the reed channel at rest. The control parameter ζ represents dimensionless reed opening at rest.
The reed is considered as a spring without mass: its position is proportional to the pressure difference between
the input of the instrument and the mouth of the musician. The nonlinear characteristic giving the ﬂow u
depending on the pressure p and the control parameters is deduced from the Bernoulli law [4] and writesζ (1 − γ + p)sign(γ − p) |γ − p| if p ≥ γ − 1.
(3)
u = F(p) =0if p < γ − 1
The resonator of the saxophone is considered as a lossless conical waveguide with a lumped mouthpiece, whose
volume equals the volume of the missing apex of the cone [3]. The pressure p at the input obeys the equation tdp
6c −cdp(t) =(t − τ̂) −p(t  )dt  ,p (t) + p+ (t − τ̂) +(4)dt
dtx1
2x1 t−τ̂
where p+ = (p + u)/2 and p− = (p − u)/2 and the parameters are x1 the length of the missing apex of the cone,τ̂ = 2/c is the duration of a round trip in the resonator of length .
Deriving equation (4) leads to an equation in p only, that can be used for harmonic balance:√√(5)p̈ + 3(1 − F  (p)) ṗ + p = p¨τ − 3(1 + F  (pτ )) p˙τ + pτ ,where pτ is the delayed pressure and the time and delay are dimensionless:c√2 √t˜ =3t , τ =
3.
(6)
x1
x1
The dot notation signals derivation with respect to the dimensionless time. With the dimensionless time, there
are only 3 parameters : the control parameters ζ and γ, and the delay τ which is proportional to the ratio /x1 .
In this paper, ζ is ﬁxed at 1 and x1 at 0.4 m, which is a typical value for alto saxophones, as we concentrateon the effects of parameters γ and .

6447

2.2 Sound synthesis
Equation (4) may be discretized to serve for sound synthesis synthesis with time step ts :cts
(pn + pn−1 − pn−2M − pn−2M−1 )In = In−1 +4x1
6cts 11(pn − un ) + (pn−2M + un−2M ) + In
pn+1 = pn + pn−2M+1 − pn−2M −x1 22un+1 = F(pn+1 ),(7)(8)(9)

where M = /(cts ) and In is the discretized version of the integral of p.
In order to explore the oscillating regimes of the model, a cartography is calculated by time-integration, describing the peak-to-peak amplitude of the last 0.1 s of the pressure signal p in simulations of 1 second of soundeach, with different values of γ and . Figure 1 presents this cartography. The mouth pressure parameter γvaries from 0.25 to 0.5, around the threshold where oscillations appear, at about 1/3. The lengths  range from0.3 m, below x1 , to 1.1 m which is approximately the total length of an alto saxophone resonator. A changein behavior appears at the apparition of the oscillation (γ  0.35), depending on the value of the length : for < x1 the amplitude of the oscillations increases smoothly with γ, and for  > 2x1 the change in amplitude isbrutal. There is an intermediary region with a smooth increase of the amplitude followed by a brutal jump.
This could indicate a change in the nature of the Hopf bifurcation, from super-critical to sub-critical. To verify
this idea, we use harmonic balance, and continuation in function of γ for 3 particular values of  :  = x1 = 0.4,
 = 0.6 m and  = 1 m.

2
1.5
10.5Figure 1. Cartography of the ﬁnal peak-to-peak amplitude of the internal pressure signal p. Red dotted lines:
particular values of  where the bifurcation is studied by continuation (ﬁgures 3 and 2).

2.3 Continuation with harmonic balance method
A Taylor-series based continuation method, called asymptotic numerical method (ANM), is used. It is based on
the numerical continuation of a system of equations
R(X, λ ) = 0

where X, R(X, λ ) ∈ Rn and λ ∈ R.

(10)

The parameter λ is the continuation parameter and the vector of unknowns is X. The solution-set of equations (10) can be represented as a collection of (one-dimensional) curves in the total space (X, λ ). The ANM
has the speciﬁcity to require a speciﬁc treatment of the equations to obtain an equivalent system of equation
with at most quadratic nonlinearities [1]. This allows to compute an approximation of the solution branches(X(a), λ (a)) on the form of Taylor series
X(a) = X0 + X1 a + X2 a2 + . . .

6448

and λ (a) = λ0 + λ1 a + λ2 a2 + . . .

(11)

where a is a parametrization of the solution-branch which veriﬁes here a = X1 (X − X0 ) + λ1 (λ − λ0 ). In thepresent work, the equation of the model is discretized in the frequency domain using harmonic balance method(HBM), that is a truncated Fourier expansion of the unknown
Hp(t) = p0 + ∑ pc,k cos(kωt) + ps,k sin(kωt)

(12)

k=1

This ansatz is put in the quadratic recast of the model equation in order to obtain a system of the form (10)where X is a vector containing the all Fourier coefﬁcients and the unknown angular frequency ω. The parameter
λ is γ or  in the following study. The truncated Fourier expansion (at order H = 50 in the applications) ofp and the angular frequency ω are then known along the solution branch. The generalization of the coupling
between HBM and ANM in the case of systems with time-delay has been done in [2].

3Results

The results in this section are presented on the form of several ﬁgures. On these ﬁgures, the amplitude and theangular frequency of the acoustic pressure p computed with the continuation method are represented in blue
solid lines. The truncation order for the Fourier series is H = 50. The red crosses are the point-wise periodic
solutions given by the sound synthesis algorithm. The red circles are quasi-periodic steady state solutions. The
comparison is made for values of the parameters along the red dotted lines of the ﬁgure 1. The agreementof both methods is very good on most of the results. A surprising feature of the system at  ﬁxed is thateven if the Hopf bifurcation seems super-critical in all cases, the larger  is, the sooner the Hopf bifurcation isfollowed by a ﬁrst fold bifurcation (or limit point), before the second one where the periodic solution retrieves
its stability. For  sufﬁciently small (see right ﬁgure 3) the two fold bifurcations are getting closer and with lesseffect on the dynamics of the system. However, in this case the periodic solution seems to become unstable
quickly to the beneﬁt of a quasi-periodic solution obtained with sound synthesis. In terms of frequency, apitch-ﬂattening effect is observed for high γ values, but before the ﬁrst fold bifurcation, the pitch rises when γaugments.

4CONCLUSIONSFor this simple saxophone model, the Hopf bifurcation is never found to be sub-critical, but hysteretic behavior
can appear due to the presence of two fold bifurcations, which are increasingly far apart as the length of theresonator augments. Further work may include a study of the dependence of these phenomena on the volume of
the mouthpiece or the reed opening control parameter ζ , as well as a study on the stability of the equilibrium
and the oscillating solution in the region where they coexist.

ACKNOWLEDGEMENTS
The authors would like to thank Christophe Vergez and Philippe Guillemain for helpful comments and discussions. This work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of theA*MIDEX project (ANR-11-IDEX-0001-02), funded by the Investissements d’Avenir French Government program managed by the French National Research Agency (ANR). This study has been supported by the French
ANR LabCom LIAMFI (ANR-16-LCV2-007-01).

REFERENCES[1] Guillot, L, Cochelin, B, Vergez, C. (2019) A generic and efﬁcient Taylor series–based continuation method using a quadratic recast of smooth nonlinear systems. Int J Numer Methods Eng. ;https://doi.org/10.1002/nme.6049

6449

Amplitude of p
Angular frequency

Amplitude of p
Angular frequency

Angular frequency

Amplitude of p

Figure 2. Bifurcation diagram representing the ﬁnal peak-to-peak amplitude and the frequency of the internalpressure signal p, for  = 1. Blue line: continuation and harmonic balance. Red crosses: sound synthesis.

Figure 3. Bifurcation diagram representing the ﬁnal peak-to-peak amplitude and the frequency of the internalpressure signal p. Blue line: continuation and harmonic balance. Red crosses: sound synthesis. Left :  = 0.6.
Right :  = 0.4.
[2] Guillot, L, Vergez, C, Cochelin, B. (2019). Continuation of periodic solutions of various types of delaydifferential equations using asymptotic numerical method and harmonic balance method. Nonlinear Dyn.;6450

https://doi.org/10.1007/s11071-019-04958-y[3] Kergomard, J, Guillemain, P, Silva, F, Karkar, S. (2016) Idealized digital models for conical reed instruments, with focus on the internal pressure waveform. The Journal of the Acoustical Society of America,
139(2) :927–937.
[4] Wilson, T, Beavers, G. (1974) Operating modes of the clarinet. The Journal of the Acoustical Society ofAmerica, 56(2) :653–658.
[5] Ollivier, L, Dalmont, J-P, Kergomard, J. (2004) Idealized models of reed woodwinds. Part I: Analogy withthe bowed string, Acta acustica united with acustica, vol. 90, no. 6, pp. 1192- 1203.
[6] Dalmont, J-P, Gilbert, J, and Kergomard, J. (2000) Reed instruments, from small to large amplitude periodic
oscillations and the Helmholtz motion analogy, Acta Acustica united with Acustica, vol. 86, no. 4, pp. 671684.
[7] Ricaud, B, Guillemain, P, Kergomard, J, Silva, F, Vergez, C. (2009) "Behavior of reed woodwind instruments around the oscillation threshold," Acta Acustica united with Acustica, vol. 95, no. 4, pp. 733-743.
[8] Karkar, S, Vergez, C, Cochelin, B. (2012) Oscillation threshold of a clarinet model: A numerical continuation approach. The Journal of the Acoustical Society of America, 131(1), 698-707.

6451