# Continuation of periodic solutions of various types of Delay Differential Equations using Asymptotic Numerical Method and Harmonic Balance Method

Louis Guillot, Christophe Vergez, Bruno Cochelin

► To cite this version:

Louis Guillot, Christophe Vergez, Bruno Cochelin. Continuation of periodic solutions of various types of Delay Differential Equations using Asymptotic Numerical Method and Harmonic Balance Method. Nonlinear Dynamics, 2019, 10.1007/s11071-019-04958-y. http://hal-01929349

HAL Id: hal-01929349

https://hal.science/hal-01929349

Submitted on 21 Nov 2018

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Continuation of periodic solutions of various types of Delay Differential Equations using Asymptotic Numerical Method and Harmonic Balance Method

Louis Guillot · Christophe Vergez · Bruno Cochelin

Received: date / Accepted: date

**Abstract** This article presents an extension of the Asymptotic Numerical Method combined with the Harmonic Balance Method to the continuation of periodic orbits of Delay Differential Equations. The equations can be forced or autonomous and possibly of neutral type. The approach developed in this paper requires the system of equations to be written in a quadratic formalism which is detailed. The method is applied to various systems, from Van der Pol and Duffing oscillators to toy models of clarinet and saxophone. The Harmonic Balance Method is ascertained from a comparison to standards time-integration solvers. Bifurcation diagrams are drawn which are sometimes intricate, showing the robustness of this method.

**Keywords** Nonlinear dynamics · Delay equations · Harmonic balance method · Numerical continuation · Periodic solutions · Quadratic recast · Asymptotic Numerical Method.

## 1 Introduction

The aim of this article is to extend the framework of the numerical continuation of periodic solutions of differential systems using the Harmonic Balance Method (HBM) coupled with the Asymptotic Numerical Method (ANM) to time-delayed systems.

For dynamical systems without delay, the efficiency of this method has already been shown in several works, see [8], [17] and [22] for some examples. The work [15] treats the case of quasi-periodic solutions using the same framework. To achieve better efficiency, the choice to recast the original system into a quadratic version has been done. The procedure to obtain such a recast has been explained in great details in [14] and is briefly recalled in this paper. The way to recast non-polynomial non-linearities was explained for the first time in [17] and a way to recast fractional order derivatives is already known from [23]. As opposed to already existing methods of continuation of delay differential equations (DDE), see [12] and its adaptation for neutral systems [2], the method developed in this paper is based on the frequency representation of the periodic solutions.

This article contains three sections. In the first section the ANM and the HBM are recalled. It is admitted that the method can be applied to virtually any explicit or implicit differential-algebraic equations. It is extended in a second section to the case of explicit or implicit delay algebraic-differential equations which can be forced or autonomous. Neutral systems can also be treated. The effect of time-delay on Fourier series is derived and the use of auxiliary variables to write the delayed variable in the frequency domain is detailed. This makes the quadratic recast of a delayed variable clear. In the last section, the method is applied to several examples, showing its versatility and its robustness particularly when applied to neutral systems. The convergence of the Fourier series with the order of truncation is shown and the solutions are compared to standard time-integration solvers. The continuation with respect to the delay is also explored in the first example.
## 2 Continuation of periodic solutions

### 2.1 Harmonic Balance Method

As described in [8], the periodic solutions $t \mapsto x(t)$ of a dynamical systems are represented by their truncated Fourier decomposition which is found through the Harmonic Balance Method :

$$x(t) = \sum_{k=-H}^{H} x_k \exp(ik\omega t) = a_0 + \sum_{k=1}^{H} a_k \cos(k\omega t) + b_k \sin(k\omega t) \quad (1)$$

where $\omega$ is the pulsation of the oscillations and $H$ is the number of harmonics i.e the order of truncation of the Fourier series. Given that $x_{-k}$ is the complex conjugate of $x_k$, the complex form and the real form of the Fourier series are equivalent, hence both can be used depending on which is more convenient.

For example, let's seek periodic solutions of

$$\ddot{x} + 2\mu\dot{x} + x + x^3 = F \cos(\omega t) \quad (2)$$

The system is rewritten quadratically as in [8] :

$$\ddot{x} + 2\mu\dot{x} + x + xy = F \cos(\omega t) \\ y = x^2 \quad (3)$$

The next step is to compute products and derivatives directly on the complex form of the Fourier series (1). The products can be written as a convolution :

$$\left(\sum_{k=-H}^{H} x_k \exp(ik\omega t)\right) \left(\sum_{k=-H}^{H} y_k \exp(ik\omega t)\right) = \sum_{k=-H}^{H} \exp(ik\omega t) \sum_{\substack{h_1+h_2=k \\ |h_1| \le H, |h_2| \le H}} x_{h_1} y_{h_2} \quad (4)$$

The derivatives can be written :

$$\begin{aligned} \dot{x}(t) &= \sum_{k=-H}^{H} ik\omega x_k \exp(ik\omega t) \\ \ddot{x}(t) &= \sum_{k=-H}^{H} -k^2\omega^2 x_k \exp(ik\omega t) \end{aligned} \quad (5)$$

Now it is possible to write (3) in the frequency domain by replacing with the formula (4) and (5) and by balancing the harmonics. This leads to a set of algebraic equations.

### 2.2 Continuation through Asymptotic Numerical Method

The Asymptotic Numerical Method (ANM) that has first been described in [4,5,6], relies on a high-order Taylor series representation of the solution-branch. This technique has already proven its efficiency for a lot of applications in engineering [10], mechanics [24] or acoustics for example. While some implementations relying on automatic differentiation do exist, see [3] for example, the choice is made here to work with a quadratic framework. A generic implementation of this latest approach which minimizes problem-dependent implementation has been developed [14]. It is based on the numerical continuation of algebraic systems of the form

$$R(V) = 0 \quad (6)$$

where $V \in \mathbb{R}^{n+1}$ and $R(V) \in \mathbb{R}^n$ is an analytic function of its argument. This system is always written in a quadratic format as a prerequisite of the method. Again, this formalism is not a constraint that we suffer but a choice that allows to treat a very wide range of problems as shown in [17]. The system (6) can then be written

$$\tilde{R}(\tilde{V}) = C + L(\tilde{V}) + Q(\tilde{V}, \tilde{V}) = 0 \quad (7)$$

where $\tilde{V} \in \mathbb{R}^{N+1}$, $N \ge n$ and $C, L$ and $Q$ are respectively a constant, a linear and a quadratic operator. For the small example (3), the introduction of an auxiliary variable $y = x^2$ in the time domain allows to have almost directly a quadratic algebraic system. The only non-quadratic terms
are due to the second order derivative (5) that leads to cubic terms $x_k\omega^2$. The introduction of the auxiliary variable $\Omega = \omega^2$ solves this problem.

The coupling between ANM and HBM in the quadratic framework (7) can be automatized and is detailed in [8]. Again, the key idea of the ANM is to compute a high-order Taylor series development of the solution-branch. Let consider that $\tilde{V}_0$ is a regular solution of the system (7). Let $\tilde{V}_1$ be a tangent vector at $\tilde{V}_0$. The classical arc-length parameter $a = \tilde{V}_1^t(\tilde{V}_0 - \tilde{V})$ is introduced and the solution-branch around $\tilde{V} = \tilde{V}_0$ is written as a power series with respect to $a$

$$ \tilde{V}(a) = \tilde{V}_0 + a\tilde{V}_1 + a^2\tilde{V}_2 + a^3\tilde{V}_3 + \dots \qquad (8) $$

More on the implementation of continuation in the case of algebraic systems coming from the harmonic balance method is available in [8]. The next section generalizes this approach for systems with time-delayed variables in their formulation.

### 3 Continuation of periodic solution of systems with time delay

In this section a system of delay differential equations, forced or autonomous and possibly neutral, of the following general form is considered:

$$ f(t, x(t), \dot{x}(t), x(t-\tau), \dot{x}(t-\tau), \lambda) = 0 \qquad (9) $$

Where $f: \mathbb{R} \times (\mathbb{R}^N)^4 \times \mathbb{R} \mapsto \mathbb{R}^N$ is an analytic function that is periodic with respect to its first variable, $x: \mathbb{R} \mapsto \mathbb{R}^N$ is an unknown periodic function of unknown pulsation $\omega$, $\tau$ is a real number called the delay and $\lambda$ is a control parameter of interest called the continuation parameter. The continuation parameter $\lambda$ may be equal to the delay $\tau$ as shown in section 4.1.

A systematic way to recast the equations (9) quadratically is first to introduce three auxiliary variables :

$$ \left\{ \begin{array}{l} f(t, x(t), y(t), x_{\tau}(t), y_{\tau}(t), \lambda) = 0 \\ y(t) - \dot{x}(t) = 0 \\ x_{\tau}(t) - x(t - \tau) = 0 \\ y_{\tau}(t) - y(t - \tau) = 0 \end{array} \right. \qquad (10) $$

The system has now $4N+1$ unknowns for $4N$ equations. The equations on delayed variables appear isolated and a classical quadratic recast [8, 17, 14] can be applied to the function $f$ to write the system (10) quadratically :

$$ \left\{ \begin{array}{l} g(t, x(t), y(t), x_{\tau}(t), y_{\tau}(t), \lambda, v(t)) = 0 \\ y(t) - \dot{x}(t) = 0 \\ x_{\tau}(t) - x(t - \tau) = 0 \\ y_{\tau}(t) - y(t - \tau) = 0 \end{array} \right. \qquad (11) $$

where $v: \mathbb{R} \mapsto \mathbb{R}^{N_{aux}}$ is a function of the time that can be expressed as quadratic expressions of all the variables. The system has now $4N+1+N_{aux}$ unknowns for $4N+N_{aux}$ quadratic equations. The two last equations are the only one that are not treated in [8, 17, 14], thus, the Fourier coefficients of $x_{\tau}$ are now derived from the Fourier coefficients of $x$ with their quadratic recast.

#### 3.1 Harmonic Balance Method applied to delayed variables

The order of truncation of the Fourier series is noted $H$. The truncated Fourier series of $t \mapsto x(t)$ is:

$$ x(t) = a_0 + \sum_{k=1}^{H} a_k \cos(k\omega t) + b_k \sin(k\omega t) \qquad (12) $$
The Fourier series of $t \mapsto x(t - \tau)$ is:

$$
\begin{align}
x(t - \tau) &= a_0 + \sum_{k=1}^{H} a_k \cos(k\omega(t - \tau)) \nonumber \\
&\quad + \sum_{k=1}^{H} b_k \sin(k\omega(t - \tau)) \nonumber \\
&= a_0 + \sum_{k=1}^{H} a_k (\cos(k\omega t) \cos(k\omega\tau) + \sin(k\omega t) \sin(k\omega\tau)) \nonumber \\
&\quad + \sum_{k=1}^{H} b_k (\sin(k\omega t) \cos(k\omega\tau) - \cos(k\omega t) \sin(k\omega\tau)) \tag{13} \\
&= a_0 + \sum_{k=1}^{H} [a_k \cos(k\omega\tau) - b_k \sin(k\omega\tau)] \cos(k\omega t) \nonumber \\
&\quad + \sum_{k=1}^{H} [b_k \cos(k\omega\tau) + a_k \sin(k\omega\tau)] \sin(k\omega t) \nonumber
\end{align}
$$

Therefore, if $a_k^\tau$ and $b_k^\tau$ are the Fourier coefficients of $t \mapsto x(t - \tau)$, their expression is given by:

$$
\begin{align}
a_0^\tau &= a_0 \nonumber \\
a_k^\tau &= a_k \cos(k\omega\tau) - b_k \sin(k\omega\tau), \quad k \ge 1 \tag{14} \\
b_k^\tau &= b_k \cos(k\omega\tau) + a_k \sin(k\omega\tau), \quad k \ge 1 \nonumber
\end{align}
$$

3.2 Asymptotic Numerical Method applied to delayed variables

For $k \ge 1$, the auxiliary variables $c_k^\tau = \cos(k\omega\tau)$ and $s_k^\tau = \sin(k\omega\tau)$ are defined to get a quadratic
recast for (14):

$$
\begin{align*}
\tilde{\tau} &= \omega\tau \\
a_0^\tau &= a_0 \\
a_k^\tau &= a_k c_k^\tau - b_k s_k^\tau \\
b_k^\tau &= b_k c_k^\tau + a_k s_k^\tau \\
s_k^\tau &= \sin(k\tilde{\tau}) \quad \text{and} \quad \mathrm{d}(s_k^\tau) = kc_k^\tau \times \mathrm{d}(\tilde{\tau}) \\
c_k^\tau &= \cos(k\tilde{\tau}) \quad \text{and} \quad \mathrm{d}(c_k^\tau) = -ks_k^\tau \times \mathrm{d}(\tilde{\tau})
\end{align*}
\tag{15}
$$

where in the last two equations, the differentiated forms are used as explained in [14] to compute
the Taylor series expansions (8) of $s_k^\tau$ and $c_k^\tau$. The last equations are defined as "initial conditions"
plus a "differential system". For the computation of the Taylor series (8), the "initial conditions"
are used to compute the zero-th order coefficient $\tilde{V}_0$ and the differential system to compute higher
order coefficients $\tilde{V}_p$, $p \ge 1$. This type of computation of Taylor series development is explained in
details in the book [13]. The main ideas are recalled below.

Let $x = x_0 + ax_1 + a^2 x_2 + \dots$ be a given Taylor series and let $s = \sin(x)$ and $c = \cos(x)$ be its
sine and its cosine. To compute the coefficients of the Taylor series of $s$ ($= s_0 + as_1 + a^2 s_2 + \dots$)
and $c$ ($= c_0 + ac_1 + a^2 c_2 + \dots$), the following procedure is applied:

$$
-s_0 = \sin(x_0) \text{ and } c_0 = \cos(x_0).
$$

- For $p \ge 1$, identify the term of order $p-1$ in the differentiated equations $ds = cdx$ and $dc = -sdx$.

For the second step, the development in Taylor series of the differentiated equations are written:

$$
\begin{align}
s_1 + 2as_2 + 3a^2 s_3 + \dots &= (c_0 + ac_1 + a^2 c_2 + \dots) (x_1 + 2ax_2 + 3a^2 x_3 + \dots) \tag{16} \\
c_1 + 2ac_2 + 3a^2 c_3 + \dots &= -(s_0 + as_1 + a^2 s_2 + \dots) (x_1 + 2ax_2 + 3a^2 x_3 + \dots) \nonumber
\end{align}
$$

For $p=1$, the identification of the constant term of the system (16) gives $s_1 = c_0x_1$ and $c_1 = -s_0x_1$.
For $p \ge 1$, the identification of the terms of order $p-1$ in the system (16) gives the following
recursive formula:

$$
\begin{align}
s_p &= \frac{1}{p} [c_0 \times px_p + c_1 \times (p-1)x_{p-1} + \cdots + c_{p-2} \times 2x_2 + c_{p-1} \times x_1] \notag \\
c_p &= -\frac{1}{p} [s_0 \times px_p + s_1 \times (p-1)x_{p-1} + \cdots + s_{p-2} \times 2x_2 + s_{p-1} \times x_1] \tag{17}
\end{align}
$$
## 4 Applications

In this section, four examples are given which have been explored with our method. The first one, a delayed Van der Pol oscillator [1], is autonomous. It is used to show the convergence of the Fourier series to a reference solution of the system, computed with a standard solver [20]. Moreover, a bifurcation diagram with respect to the delay is drawn and the location of the many Hopf bifurcations is ascertained with the results of a linear stability analysis. The second example, a forced delayed Duffing oscillator [16], exhibits a lot of resonances (both super-harmonic and sub-harmonic). Its bifurcation diagram is represented with a focus on two symmetry-breaking bifurcations. The third example is a toy model of clarinet [21] that is purely algebraic. Its periodic-solutions are squared shaped and their Fourier series converge rather slowly. This example shows the robustness of our approach. The fourth and last example is a neutral model of saxophone [18,9] using transcendental functions in its writing. A solution is compared to a finite-difference scheme from the literature [18] and a bifurcation diagram is drawn, showing an expected behaviour for this type of system according to [11].

### 4.1 A delayed Van der Pol oscillator

In this part we consider the Van der Pol oscillator (dimensionless) with delayed feedback that has been studied in [1]:

$$\ddot{u}(t) - \varepsilon(1 - u(t)^2)\dot{u}(t) + u(t) = \varepsilon ku(t - \tau) \quad (18)$$

The continuation is performed with respect to $\tau$. Introducing the auxiliary variables $v(t) = \dot{u}(t)$, $r(t) = 1 - u(t)^2$ and $F_{nl}(t) = v(t)r(t) - ku(t - \tau)$. In order to minimize the number of auxiliary variables, the variable $F_{nl}$ is the sum of the delayed variable $ku(t - \tau)$ (with $k$ constant) and the product $v \times r$. Hence the delayed variable is not isolated like in (10) and (11), but as it appears linearly, this allows to use the approach explained in the previous section. The following quadratic recast is used:

$$\begin{aligned}
\dot{v}(t) + u(t) - \varepsilon F_{nl}(t) &= 0 \\
v(t) - \dot{u}(t) &= 0 \\
r(t) - 1 + u(t)^2 &= 0 \\
F_{nl}(t) - v(t)r(t) - ku(t - \tau) &= 0
\end{aligned} \quad (19)$$

To validate the method, the figure 1 shows the difference between time integration and harmonic balance for the periodic solution of equation (18) obtained with $\varepsilon = 2$, $\tau = 0.2$ and $k = 0.2$. The number of harmonics is varied to show the convergence of the harmonic balance towards the time integration solution. The reference solution is a time integration performed with dde23 built-in solver of MATLAB. Its absolute and relative tolerance thresholds have been set to $10^{-8}$.

With a low number of harmonics, i.e for $H < 10$, the solution is not accurate and even the period of the oscillation is not computed correctly. Indeed, the plot $H = 5$ stops before $t = 10$ s, conversely to the other plots. From $H = 25$ to $H = 50$, the solution gets much closer to the reference. Solutions for $H = 100$ and $H = 200$ look superimposed. In fact, the tolerance at which the reference has been computed, i.e $10^{-8}$, is reached for $H = 50$. When computing the solution with a higher number of harmonics, the time-integration reference is probably farther from the true solution than the solution obtained with the harmonic balance method. In addition, the time-integration solver used takes approximately ten times the time needed to compute the solution with $H = 50$. Moreover, once a solution has been computed the computation of the bifurcation diagram through ANM is very fast: for the same computation time as one time-integration over twenty periods (to get rid of the transient), all the steady state solutions between $\varepsilon = 2$ and $\varepsilon = 5$ are available with our method.

The figure 2 shows a part of the bifurcation diagram of the system (18) with the delay $\tau$ as a continuation parameter. There are infinitely many Hopf bifurcations in this system. The analytic position of the Hopf bifurcations is determined by the stability analysis of the equilibrium $u = 0$ and is given by the solutions of the system [1]:

$$\left\{ \begin{array}{l} 1 - \omega^2 = k\varepsilon \cos(\tau\omega) \\ \varepsilon\omega = k\varepsilon \sin(\tau\omega) \end{array} \right. \quad (20)$$
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Differences over one period between several solutions of equation (18) obtained with the harmonic balance method with increasing *H* and a reference solution obtained with time integration. Here *ε* = 2, *τ* = 0.2 and *k* = 0.2.</figcaption></figure>

where *ω* is the pulsation of the arising periodic solution at the corresponding Hopf point. The continuation of the periodic solution branches has been initialized with an analytic approximation with *H* = 1 and the pulsation *ω* solution of the system (20) for the corresponding value of *τ*. In this particular case, the robustness of the ANM allows to follow all the solution branches from just one starting point. In the general case, it would probably necessitate to initialize with an analytic approximation at every Hopf bifurcations. On the figure 2 some parts of the diagram are of null amplitude : these are just some parts where the solution found is the equilibrium *u* = 0.

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Bifurcation diagram of equation (18) with respect to the delay $\tau$, $\varepsilon = 0.1$ and $k = 2$. There are as many Hopf bifurcations as branches starting from a null amplitude. The red dots are the analytical position of Hopf bifurcations given by the solutions of system (20).</figcaption></figure>
## 4.2 A forced delayed Duffing oscillator

The example presented in this section is taken from [16]. It is an example of a forced delayed differential equation. The dimensionless equation of motion (Eq.(3) in the reference) is:

$$ \ddot{x}(t) + 2\zeta\dot{x}(t) + x(t) + \mu x^3(t) = 2ux(t - \tau) + 2v\dot{x}(t - \tau) + f \cos(\lambda t) \quad (21) $$

The continuation is performed with respect to $\lambda$. The following quadratic recast is used:

$$ \begin{aligned} \dot{x}(t) - y(t) &= 0 \\ \dot{y}(t) + 2\zeta y(t) + x(t) + \mu r(t)x(t) - 2ux(t - \tau) - 2vy(t - \tau) &= f \cos(\lambda t) \\ r(t) - u(t)^2 &= 0 \end{aligned} \quad (22) $$

where $\lambda$ is the continuation parameter and the set of parameters used is: $\zeta = \varepsilon, \mu = \varepsilon, u = \varepsilon, v = -\varepsilon, \tau = \pi$ and $f = 5$ with $\varepsilon = 0.05$. Again, the time-delayed terms are not isolated but as they appear linearly in the equations, the recast of section 3 is applicable. In [16], the primary and the 1/3 subharmonic resonances were studied with the method of multiple scales. In this present article a whole bifurcation diagram is drawn, showing both primary and 1/3 subharmonic resonances. To find the starting point of the continuation, the approximations computed in [16] are used. The forcing amplitude has been chosen high to have a rich dynamics as can be seen in the bifurcation diagram which is represented in figure 3. To draw this diagram, the order of truncation of the Fourier series has been set to $H = 25$ which is sufficient away from superharmonic resonances. The 1/3 subharmonic resonance is disconnected from the main branch. It corresponds to a solution-branch where the pulsation of the oscillation verifies $\omega = \frac{\lambda}{3}$.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Bifurcation diagram of the amplitude of the solution with respect to the forcing pulsation $\lambda$. On the primary resonance solution-branch, there are several super-harmonic resonances appearing for $\lambda = \frac{1}{3}, \frac{1}{5}, \frac{1}{7}, \frac{1}{9}, \dots$. There are also two symmetry-breaking bifurcation, noted SB on the zoom of the diagram, detected with the method explained in [7] for which even harmonics appear.</figcaption></figure>

The zoom on the top of figure 3 shows the many superharmonic resonances and the symmetry-breaking (SB) bifurcations. The superharmonic resonances appear for $\lambda \gtrsim \frac{1}{2p+1}, p \in \mathbb{N}$ and in these areas the solutions become more complex as the harmonics $2p + 1$ are close to the natural pulsation of the system. In addition, two symmetry-breaking bifurcations are detected and the emerging solution is then continued with the method explained in [7]. On the bifurcated branches, the solutions have strong even harmonics appearing and their phase diagrams $(x, \dot{x})$ are no longer symmetric.
## 4.3 A clarinet toy model

The following example is a simplified model of clarinet-like instrument, taken from [21] [Eq.(4)-(5)-(6)-(7)-(10)]. This system has been chosen to highlight that our approach copes even with a degenerate form of the general system (9). In this example, there are no time derivatives in the equations of the motion. The unknowns are functions of the time and the dynamics comes from the coupling between a non-linearity and an equation with time delayed variable. As in the reference, the reed dynamics is neglected and the dimensionless equations can be written:

$$ 
\begin{aligned}
u(t) &= \begin{cases} 0 & \text{if } p(t) \le \gamma - 1 \\ \zeta(p(t) - \gamma + 1)\sqrt{|\gamma - p(t)|}\operatorname{sign}(\gamma - p(t)) & \text{if } p(t) > \gamma - 1 \end{cases} \\
p^{-}(t) &= -\lambda p^{+}(t - \frac{2\ell}{c})
\end{aligned} 
\quad (23) $$

where $p$ and $u$ are respectively the acoustic pressure and the volume flow in the mouthpiece of the instrument. $p^+$ and $p^-$ are respectively the forward and the backward progressive waves, linked to $p$ and $u$ through $p = p^- + p^+$ and $u = p^+ - p^-$. $\ell$ is the length of the instrument and $c$ is the speed of sound. The delay $\frac{2\ell}{c}$ in equation (23) represents the time for a pressure wave to go back and forth the cylindrical bore of the clarinet. $\lambda$ is the reflection coefficient at the end of the instrument which can be seen as a damping term. $\gamma$ is the static blowing pressure in the mouth of the player and $\zeta$ is the reed parameter. All the variables but the time are dimensionless in this model. In the following, the time is normalized by $\frac{2\ell}{c}$ so that the value of the delay $\tau$ is 1. For simplicity, the nonlinear characteristics (the first equation of the system (23)) is written under the assumption that $-1 < p - \gamma < 0$ which means according to [21] that the reed channel is open and air is allowed to enter the clarinet. The system (23) becomes

$$ 
\begin{aligned}
u(t) &= \zeta(p(t) - \gamma + 1)\sqrt{\gamma - p(t)} \\
p^{-}(t) &= -\lambda p^{+}(t - 1)
\end{aligned} 
\quad (24) $$

Its quadratic recast is

$$ 
\begin{aligned}
p^{-}(t) &= -\lambda p^{+}(t - 1) \\
v(t)^2 &= \gamma - p(t) \\
u(t) &= \zeta(p(t) - \gamma + 1)v(t)
\end{aligned} 
\quad (25) $$

The variable $v(t)$ is defined implicitly by the second equation. To ensure that the results are correct, we have to check that $v(t)$ is non-negative. The parameters chosen are $\lambda = 0.95$ and $\zeta = 0.8$. The continuation is performed with respect to $\gamma$, the dimensionless blowing pressure. The bifurcation diagram is plotted in figure 4.

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Bifurcation diagram of the toy model of clarinet (24) with respect to the dimensionless blowing pressure $\gamma$, with $H = 200$. The period doubling bifurcation is visible as the connection of the branches. The red dotted lines correspond to section of branches where the variable $v(t)$ is negative, i.e where the simplified model is not valid anymore to represent the clarinet functioning.</figcaption></figure>
The values for which the system encounters a Hopf bifurcation $\gamma_H = 0.35$ or a period doubling bifurcation $\gamma_{PD} = 0.43$ are the same as the one given in [21] (table 1) up to the third significant digit. On the figure 4, the Hopf bifurcation occurs at the lowest point of the solution-branch. As the mean value of $p$ is not zero, the $L^2$ norm of $p$ is not 0 at the bifurcation point. The exact solution of this system in the time-domain is a square and the decomposition in Fourier series is obviously not very well suited for this discontinuous type of signal. However, the numerical method is so robust that it is able to continue this solution on the whole domain which is remarkable. This also explains why the diagram represents the $L^2$ norm of $p$ and not its amplitude : Gibbs phenomenon leads to an overshoot at each of the two discontinuous point of the signal. This gives wrong amplitude for the acoustic pressure. The analytical solutions given in [21] were used as starting points for the continuation process, using their truncated decomposition in Fourier series.

## 4.4 A neutral saxophone model

In this section, a neutral system that is a dimensionless model of saxophone taken from [18] is studied. It has also been the subject of the conference talk [9]. The derivative of Eq.(40) in [18] is taken to get rid of the integrals. This will necessitate to be taken into account when considering the Fourier developments of the unknowns. The following equation is obtained :

$$ \ddot{p} + \sqrt{3}(1 - F'(p))\dot{p} + p = \ddot{p}_{\tau} + \sqrt{3}(-1 - F'(p_{\tau}))\dot{p}_{\tau} + p_{\tau} \quad (26) $$

where $p$ is the acoustic pressure in the mouthpiece of the instrument; $\gamma$ and $\zeta$ have the same meaning as in section 4.3. Everytime there is a $\tau$ index, it means that the variable is delayed, i.e $p_{\tau}(t) = p(t - \tau)$. In order to have a smooth function, the nonlinear characteristics $F(p)$ is an exponential regularized law [9] :

$$ F(p) = \zeta \frac{\exp(1)}{\sqrt{3}} \exp(3(p - \gamma))(\gamma - p) \quad (27) $$

To recast the equation (26) quadratically, the idea is to write it on the form $I + \sqrt{3}\dot{p} = I_{\tau} - \sqrt{3}\dot{p}_{\tau}$ with $I = \ddot{p} - \sqrt{3}F'(p)\dot{p} + p$. It is then necessary to write $F'(p)$ quadratically. This is done using the same idea as for the recast in equation (15) ; setting $\Delta p = p - \gamma$, the exponential is written using $A = \exp(3\Delta p)$ and its differentiated form $dA = 3Ad(\Delta p)$. These remarks lead to introduce the following auxiliary variables :

$$ \begin{aligned} \Delta p &= p - \gamma \\ A &= \exp(3\Delta p) \\ F'(p) &= \zeta \frac{\exp(1)}{\sqrt{3}} A(-1 - 3\Delta p) \\ q &= \dot{p} \\ I &= \ddot{p} + p - \sqrt{3}F'(p)q \end{aligned} \quad (28) $$

which allows to write the equation (26) in a quadratic format :

$$ \begin{aligned} I - I_{\tau} + \sqrt{3}(q + q_{\tau}) &= 0 \\ \Delta p - (p - \gamma) &= 0 \\ A - \exp(3\Delta p) &= 0 \quad \text{and} \quad dA - 3A \, d(\Delta p) = 0 \\ F'(p) - \zeta \frac{\exp(1)}{\sqrt{3}} A(-1 - 3\Delta p) &= 0 \\ q - \dot{p} &= 0 \\ I - (\ddot{p} + p - \sqrt{3}F'(p)q) &= 0 \end{aligned} \quad (29) $$

The definition of $A$ (third line, left hand side in (29)) is used only to compute the mean value (or the constant term) of its Fourier series and its differentiated form (third line, right hand side in (29)) is used to compute higher order terms. This has been explained in [17]. The dimensionless delay $\tau$ can be expressed in term of the length of the conical resonator $\ell$ of the saxophone, truncated at its input to plug the mouthpiece, and the length $x_1$ of the missing part of the cone (from the apex to the input of the truncated cone). This is not our purpose here to understand the refinements of this model, the interested reader can find the details in [18].

Below, the value of $\tau$ is set to $\frac{2\sqrt{3\ell}}{x_1}$ with $\ell = \frac{17}{15}$ and $x_1 = \frac{1}{5}$. $\zeta$ is set to 1.2 to draw the bifurcation diagram with respect to $\gamma$ in figure 5. As said above, the equation (26) is the derivative of another
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Bifurcation diagram of a model of saxophone with respect to the dimensionless blowing pressure $\gamma$ and for $\zeta = 1.2$. The sub-critical Hopf bifurcations can be seen clearly. It is not surprising for the saxophone to show such sub-critical bifurcations [11]. This diagram has been computed using $H = 50$ harmonics.</figcaption></figure>

equation. This implies that after the Harmonic Balance Method (HBM) has been performed, the constant coefficient equation (i.e its mean value) states $0 = 0$. As this leads to a singular system, this equation is replaced by $p^0 = 0$ where $p^0$ is the constant term [the $a_0$ of the Fourier decomposition (12)] of $p$.

The waveform of each unknown over one period is a by-product of the method and is available for each value of $\gamma$, i.e for each point of the bifurcation diagram 5. The waveform of $p$ is compared to the one obtained with the numerical integration algorithm presented in [18] for $\gamma = 0.5$ and is in very good agreement as can be seen in figure 6. On the same figure, the pressure obtained by integration with the built-in solver `ddensd` [19] of MATLAB is presented. The three methods give very close results. However, as expected in [18], their scheme is unstable and the Hopf bifurcation

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Comparison of the periodic solutions of equations (26) over one period for $\zeta = 1.2$ and $\gamma = 0.5$. In blue, the solution obtained with Harmonic Balance Method with $H = 50$ and extracted from the bifurcation diagram shown in figure 5. In red dashed, the solution obtained with the time-integration algorithm of [18] with 5000 samples per period. In yellow dotted, the solution obtained with `ddensd` Matlab's solver with absolute and relative thresholds set to $10^{-3}$. A zoom on the zone of beating reed is shown to enhance the differences between the waveforms.</figcaption></figure>
is detected around $\gamma = 0.285$, that is below the theoretical value of $\gamma = \frac{1}{3}$ which is found with our method.

The scheme used in [18] is very fast but requires a high sampling frequency (more than 1000 samples per period) to converge. The *ddensd* solver takes a longer time especially when the tolerance chosen is decreased to the minimum value possible of $10^{-5}$. In our approach, the transient is not available but the shape of the steady state is obtained on all the solution-branch instantly through the Fourier series construction, as well as the unstable periodic regimes. In addition, the Hopf bifurcation and the limit point can be deduced from the bifurcation diagram.

The model of saxophone presented should have a zero mean over one period by physical consideration [18]. For the signals represented in figure 6, the acoustic pressure obtained by Harmonic Balance has a mean value around the numerical precision $10^{-16}$ while the time integration solver from [18] gives a mean around $10^{-4}$ and *ddensd* gives $10^{-3}$.

## 5 Conclusion

An extension of the HBM, coupled with the ANM continuation, to time-delayed systems is developed. Various types of systems exhibiting constant time-delays is explored. Especially, delayed differential equations of neutral type do not require a specific setting in the framework presented. The available tools of the ANM continuation allow to compute bifurcated branches without the need of a specific implementation.

The coupling between HBM and ANM shows its efficiency and its robustness compared to time-integration methods. Especially in some cases where there are not many numerical methods available (neutral systems) and even when the Fourier series of the periodic solution have bad convergence properties. It is ascertained by a comparison to time-domain solvers. The algorithms used for this article are freely available online at http://manlab.cnrs-mrs.fr.

**Acknowledgements** The authors want to thank Tom Colinot for the help with the saxophone model and its associated time-integration method. This work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX* project (ANR-11-IDEX-0001-02), funded by the Investissements d'Avenir French Government program managed by the French National Research Agency (ANR).

## References

1. Atay, F.M.: Van der pol's oscillator under delayed feedback. Journal of Sound and Vibration **218**(2), 333-339 (1998)

2. Barton, D.A., Krauskopf, B., Wilson, R.E.: Collocation schemes for periodic solutions of neutral delay differential equations. Journal of Difference Equations and Applications **12**(11), 1087-1101 (2006)

3. Charpentier, I., Cochelin, B.: Towards a full higher order ad-based continuation and bifurcation framework. Optimization Methods & Software (2018)

4. Cochelin, B.: A path-following technique via an asymptotic-numerical method. Computers & structures **53**(5), 1181-1192 (1994)

5. Cochelin, B., Damil, N., Potier-Ferry, M.: Asymptotic-numerical methods and pade approximants for non-linear elastic structures. International journal for numerical methods in engineering **37**(7), 1187-1213 (1994)

6. Cochelin, B., Damil, N., Potier-Ferry, M.: The asymptotic-numerical method: an efficient perturbation technique for nonlinear structural mechanics. Revue européenne des éléments finis **3**(2), 281-297 (1994)

7. Cochelin, B., Medale, M.: Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations. Journal of Computational Physics **236**, 594-607 (2013)

8. Cochelin, B., Vergez, C.: A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions. Journal of sound and vibration **324**(1), 243-262 (2009)

9. Colinot, T., Kergomard, J.: Formulation analytique de la pression interne d'un modèle idéalisé d'instrument conique à anche. In: CFA 2018/VISHNO (2018)

10. Cottanceau, E., Thomas, O., Véron, P., Alochet, M., Deligny, R.: A finite element/quaternion/asymptotic numerical method for the 3d simulation of flexible cables. Finite Elements in Analysis and Design **139**, 14-34 (2018)

11. Dalmont, J.P., Gilbert, J., Kergomard, J.: Reed instruments, from small to large amplitude periodic oscillations and the helmholtz motion analogy. Acta Acustica united with Acustica **86**(4), 671-684 (2000)

12. Engelborghs, K., Luzyanina, T., Roosevelt, D.: Numerical bifurcation analysis of delay differential equations using dde-biftool. ACM Transactions on Mathematical Software (TOMS) **28**(1), 1-21 (2002)

13. Griewank, A., Walther, A.: Evaluating derivatives: principles and techniques of algorithmic differentiation, vol. 105. Siam (2008)

14. Guillot, L., Cochelin, B., Vergez, C.: A generic and efficient Taylor series based continuation method using a quadratic recast of smooth nonlinear systems (2018). URL https://hal.archives-ouvertes.fr/hal-01827832. Submitted ; https://hal.archives-ouvertes.fr/hal-01827832
15. Guillot, L., Vigué, P., Vergez, C., Cochelin, B.: Continuation of quasi-periodic solutions with two-frequency harmonic balance method. Journal of Sound and Vibration **394**, 434–450 (2017)

16. Hu, H., Dowell, E.H., Virgin, L.N.: Resonances of a harmonically forced duffing oscillator with time delay state feedback. *Nonlinear Dynamics* **15**(4), 311–327 (1998)

17. Karkar, S., Cochelin, B., Vergez, C.: A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities. *Journal of Sound and Vibration* **332**(4), 968–977 (2013)

18. Kergomard, J., Guillemain, P., Silva, F., Karkar, S.: 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 (2016)

19. Shampine, L.F.: Dissipative approximations to neutral ddes. Applied Mathematics and Computation **203**(2), 641–648 (2008)

20. Shampine, L.F., Thompson, S.: Solving ddes in matlab. Applied Numerical Mathematics **37**(4), 441–458 (2001)

21. Taillard, P.A., Kergomard, J., Laloë, F.: Iterated maps for clarinet-like systems. *Nonlinear dynamics* **62**(1-2), 253–271 (2010)

22. Vigué, P., Vergez, C., Karkar, S., Cochelin, B.: Regularized friction and continuation: Comparison with coulomb's law. *Journal of Sound and Vibration* **389**, 350–363 (2017)

23. Vigué, P., Vergez, C., Lombard, B., Cochelin, B.: Continuation of periodic solutions for systems with fractional derivatives. *Nonlinear Dynamics* pp. 1–15 (2018)

24. Zahrouni, H., Cochelin, B., Potier-Ferry, M.: Computing finite rotations of shells by an asymptotic-numerical method. *Computer methods in applied mechanics and engineering* **175**(1-2), 71–85 (1999)