Continuation of periodic solutions of various types of delay differential equations using asymptotic numerical method and harmonic balance method

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.


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-delay systems.
The HBM, a frequency-based method to search for periodic solutions of dynamical systems, has its roots in Fourier theory [13]. The works [20,22] developed the HBM as we know it today. Its application to delay differential equations (DDE) has been done in [21], using an alternating frequency-time scheme. The originality of our work is that all the steps of resolution are treated directly in the frequency domain without any incursion in the time domain and thus without possible aliasing errors. From this point of view, the order of truncation of the Fourier series is the only source of error.
The ANM, a Taylor series-based continuation method, was first introduced in [4,5] for the computation of equilibrium of nonlinear structures [30]. This is an alternative to more common predictor-corrector methods [23], and it can be seen as a high-order predictor that does not need any correction step. A few developments based on this method are given below.
For dynamical systems without delay, the efficiency of the ANM coupled with the HBM has already been shown in several works (see [8,18] and [27] for some examples). The work [16] treats the case of quasi-periodic solutions using the same framework. To achieve better efficiency, the choice to recast the origi-nal system into a quadratic version has been done. The procedure to obtain such a recast has been explained in great detail in [15] and is briefly recalled in this paper. The way to recast non-polynomial nonlinearities was explained for the first time in [18], and a way to recast fractional order derivatives is already known from [28]. 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. For smooth systems, this often is an advantage as the convergence properties of the frequency representation (the Fourier series) are very good. In addition, the robustness of the ANM will be shown to overcome some inherent difficulties of DDE.
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.

Harmonic balance method
As described in [8], the periodic solutions t → x(t) of a dynamical systems are represented by their trun-cated Fourier decomposition which is found through the harmonic balance method: where ω 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 us seek periodic solutions of The system is rewritten quadratically as in [8] : The next step is to compute products and derivatives directly on the complex form of Fourier series (1). The products can be written as a convolution : The derivatives can be written : Now it is possible to write (3) using Fourier series with formula (4) and (5). For example, the first equation gives Balancing the harmonics yields This is a set of algebraic equations in the unknowns x −H , x −H +1 , . . . , x 0 , x 1 , . . . , x H and possibly ω, μ and F.

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 [29] 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 [15]. It is based on the numerical continuation of algebraic systems of the form where V ∈ R n+1 and R(V ) ∈ 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 [18]. System (8) can then be writteñ whereṼ ∈ R N +1 , N ≥ n and C, L and Q are respectively a constant, a linear and a quadratic operator. For 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 nonquadratic terms are due to derivatives (5) that lead to cubic terms x k ω 2 or x k μω depending on whether ω or/and μ is/are unknown/unknowns as can be seen on algebraic system (7) obtained from HBM. The introduction of the auxiliary variables Ω = ω 2 and Λ = μω solves this problem.
The coupling between ANM and HBM in quadratic framework (9) can be automatized and is detailed in [8]. Again, the key idea of the ANM is to compute a highorder Taylor series development of the solution branch. Let consider thatṼ 0 is a regular solution of system (9). LetṼ 1 be a tangent vector atṼ 0 . The classical arclength parameter a =Ṽ t 1 (Ṽ 0 −Ṽ ) is introduced, and the solution branch aroundṼ =Ṽ 0 is written as a power series with respect to ã Introducing series (10) in Eq. (9) and equating the terms of same power yields : where J 0 is the Jacobian matrix of the system at the pointṼ 0 . The order 0 equation states thatṼ 0 is a solution of Eq. (9), and the order 1 equation states thatṼ 1 is a tangent vector at the solution pointṼ 0 . Notice that for orders p ≥ 2, solving linear systems (11) require only one matrix J 0 to be inverted, which is an obvious advantage from a computational cost point of view. 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.
where f : R × R N 4 × R → R N is an analytic function that is periodic with respect to its first variable, x : R → R N is an unknown periodic function of unknown pulsation ω, τ is a real number called the delay and λ is a control parameter of interest called the continuation parameter. The continuation parameter λ may be equal to the delay τ as shown in Sect. 4.1.
A systematic way to recast Eq. (12) quadratically is first to introduce three auxiliary variables : The system has now 4N + 1 unknowns for 4N equations. The equations on delayed variables appear isolated, and a classical quadratic recast [8,15,18] can be applied to the function f to write system (13) quadratically : where v : R → 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,15,18]; thus, the Fourier coefficients of x τ are now derived from the Fourier coefficients of x with their quadratic recast.

Harmonic balance method applied to delayed variables
The order of truncation of the Fourier series is noted H . The truncated Fourier series of t → x(t) is: The Fourier series of t → x(t − τ ) is: Therefore, if a τ k and b τ k are the Fourier coefficients of t → x(t − τ ), their expression is given by:

Asymptotic numerical method applied to delayed variables
For k ≥ 1, the auxiliary variables c τ k = cos(kωτ ) and s τ k = sin(kωτ ) are defined to get a quadratic recast for (17): where in the last two equations, the differentiated forms are used as explained in [15] to compute the Taylor series expansions (10) of s τ k and c τ k . The last equations are defined as "initial conditions" plus a "differential system". For the computation of Taylor series (10), the "initial conditions" are used to compute the zerothorder coefficientṼ 0 and the differential system to compute higher-order coefficientsṼ p , p ≥ 1. This type of computation of Taylor series development is explained in detail in the book [14]. The main ideas are recalled below.
Let x = x 0 + ax 1 + a 2 x 2 + · · · 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 + · · · ) and c (= c 0 + ac 1 + a 2 c 2 + · · · ), the following procedure is applied: s 0 = sin(x 0 ) and c 0 = cos(x 0 ). -For p ≥ 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 is written: For p = 1, the identification of the constant term of system (19) gives s 1 = c 0 x 1 and c 1 = −s 0 x 1 . For p ≥ 1, the identification of the terms of order p − 1 in system (19) gives the following recursive formula:

Applications
In this section, four examples are given which have been explored with our method. The first one, a delay 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 [25]. 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 delay Duffing oscillator [17], exhibits a lot of resonances (both superharmonic and subharmonic). Its bifurcation diagram is represented with a focus on two symmetry-breaking bifurcations. The third example is a toy model of clarinet [26] 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 [9,19] using transcendental functions in its writing. A solution is compared to a finite-difference scheme from the literature [19] and a bifurcation diagram is drawn, showing an expected behaviour for this type of system according to [11].

A delay Van der Pol oscillator
In this part, we consider the Van der Pol oscillator (dimensionless) with delay feedback that has been studied in [1]: The continuation is performed with respect to τ . Introducing the auxiliary variables v(t) =u(t), In order to minimize the number of auxiliary variables, the variable F nl is the sum of the delayed variable ku(t − τ ) (with k constant) and the product v × r . Hence the delayed variable is not isolated like in (13) and (14), but as it appears linearly, this allows to use the approach explained in the previous section. The following quadratic recast is used: To validate the method, Fig. 1 shows the difference between time integration and harmonic balance for the periodic solution of Eq. (21) obtained with ε = 2, τ = 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 timeintegration reference is probably farther from the true  Figure 2 shows a part of the bifurcation diagram of system (21) with the delay τ 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]: 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 system (23) 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 Fig. 2, some parts of the diagram are of null amplitude : these are just some parts where the solution found is the equilibrium u = 0.

A forced delay Duffing oscillator
The example presented in this section is taken from [17]. It is an example of a forced delay differential equation. The dimensionless equation of motion (Eq. (3) in the reference) is : The continuation is performed with respect to λ. The following quadratic recast is used : where λ is the continuation parameter and the set of parameters used is : ζ = ε, μ = ε, u = ε, v = −ε, τ = π and f = 5 with ε = 0.05. Again, the time-delayed terms are not isolated but as they appear linearly in the equations, the recast of Sect. 3 is applicable. In [17], 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 [17] 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 Fig. 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 ω = λ 3 .
The zoom on the top of Fig. 3 shows the many superharmonic resonances and the symmetry-breaking (SB) bifurcations. The superharmonic resonances appear for λ 1 2 p+1 , p ∈ N and in these areas the solutions become more complex as the harmonics 2 p + 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,ẋ) are no longer symmetric.

A clarinet toy model
The following example is a simplified model of clarinet-like instrument, taken from [26] [Eqs. (4)-(5)-(6)-(7)-(10)]. This system has been chosen to highlight that our approach copes even with a degenerate form of general system (12). 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 nonlinearity and an equation with time-delayed variable. As in the reference, the reed dynamics is neglected and the dimensionless equations can be written : Fig. 3 Bifurcation diagram of the amplitude of the solution with respect to the forcing pulsation λ. On the primary resonance solution branch, there are several superharmonic resonances appearing for λ = 1 3 , 1 5 , 1 7 , 1 9 , . . .. 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 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 − . is the length of the instrument and c is the speed of sound. The delay 2 c in Eq. (26) represents the time for a pressure wave to go back and forth the cylindrical bore of the clarinet. λ is the reflection coefficient at the end of the instrument which can be seen as a damping term. γ is the static blowing pressure in the mouth of the player and ζ is the reed parameter. All the variables but the time are dimensionless in this model. In the following, the time is normalized by 2 c so that the value of the delay τ is 1. For simplicity, the nonlinear characteristics (the first equation of system (26)) is written under the assumption that −1 < p − γ < 0 which means according to [26] that the reed channel is open and air is allowed to enter the clarinet. System (26) becomes Its quadratic recast is 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 λ = 0.95 and ζ = 0.8. The continuation is performed with respect to γ , the dimensionless blowing pressure. The bifurcation diagram is plotted in Fig. 4. The values for which the system encounters a Hopf bifurcation γ H = 0.35 or a period doubling bifurcation γ PD = 0.43 are the same as the one given in [26] (table 1) up to the third significant digit. On Fig. 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 [26] were used as starting points for the continuation process, using their truncated decomposition in Fourier series.

A neutral saxophone model
In this section, a neutral system that is a dimensionless model of saxophone taken from [19] is studied. It has also been the subject of the conference talk [9]. The derivative of Eq. (40) in [19] is taken to get rid of the integrals. This will be taken into account when considering the Fourier developments of the unknowns. The following equation is obtained : where p is the acoustic pressure in the mouthpiece of the instrument; γ and ζ have the same meaning as in Sect. 4.3. Everytime there is a τ index, it means that the variable is delayed, i.e. p τ (t) = p(t − τ ). In order to have a smooth function, the nonlinear characteristics F( p) is an exponential regularized law [9] : To recast Eq. (29) quadratically, the idea is to write it on the form I It is then necessary to write F ( p) quadratically. This is done using the same idea as for the recast in Eq. (18) ; setting p = p − γ , the exponential is written using A = e 3( p−γ ) (3 p) and its dif-ferentiated form d A = 3Ad( p). These remarks lead to introduce the following auxiliary variables : which allows to write Eq. (29) in a quadratic format: The definition of A (third line, left-hand side in (32)) 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 (32)) is used to compute higher-order terms. This has been explained in [18]. The dimensionless delay τ can be expressed in terms of the length of the conical resonator 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 [19].  Below, the value of τ is set to 2 √ 3 x 1 with = 17 15 and x 1 = 1 5 . ζ is set to 1.2 to draw the bifurcation diagram with respect to γ in Fig. 5. As said above, Eq. (29) is the derivative of another 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 Fourier decomposition (15)] of p.
The waveform of each unknown over one period is a by-product of the method and is available for each value of γ , i.e. for each point of bifurcation diagram 5. The waveform of p is compared to the one obtained with the numerical integration algorithm presented in [19] for γ = 0.5 and is in very good agreement as can be seen in Fig. 6. On the same figure, the pressure obtained by integration with the built-in solver ddensd [24] of MAT-LAB is presented. The three methods give very close results. However, as expected in [19], their scheme is unstable and the Hopf bifurcation is detected around γ = 0.285, that is, below the theoretical value of γ = 1 3 which is found with our method.
The scheme used in [19] 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 [19]. For the signals represented in Fig. 6, the acoustic pressure obtained by Harmonic Balance has a mean value around the numerical precision 10 −16 while the time integration solver from [19] gives a mean around 10 −4 and ddensd gives 10 −3 .

Conclusion
An extension of the HBM, coupled with the ANM continuation, to time-delay systems is developed. Various types of systems exhibiting constant time delays is explored. Especially, delay 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 timeintegration 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.
The main perspective of this work would be to implement a robust frequency-based method for the stability analysis of delay differential equations. For equilibrium as well as for periodic solution, this would allow full bifurcation analysis of delay systems.