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 Diﬀerential 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 Duﬃng os-cillators 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-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 quasiperiodic 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 nonpolynomial 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 differentialalgebraic 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 truncated 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's 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 the Fourier series (1). The products can be written as a convolution : The derivatives can be written :ẋ 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.

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 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 [17]. The system (6) 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 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 ω 2 . The introduction of the auxiliary variable Ω = ω 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Ṽ 0 is a regular solution of the system (7). LetṼ 1 be a tangent vector atṼ 0 . The classical arc-length parameter a =Ṽ t 1 (Ṽ 0 −Ṽ ) is introduced and the solution-branch aroundṼ =Ṽ 0 is written as a power series with respect to ã 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.

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: 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 section 4.1.
A systematic way to recast the equations (9) 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,17,14] can be applied to the function f to write the system (10) quadratically : where v : R → R Naux 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 τ 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 (14) 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 and c τ k . 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Ṽ 0 and the differential system to compute higher order coefficientsṼ p , p ≥ 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 + . . . 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: For the second step, the development in Taylor series of the differentiated equations are written: For p = 1, the identification of the constant term of the system (16) 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 the system (16) gives the following recursive formula:

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 periodicsolutions 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].

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]:ü 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 (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:v 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 ε = 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 time-integration reference is probably farther from the true solution than the solution obtained with the harmonic balance method. In addition, the timeintegration 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 ε = 2 and ε = 5 are available with our method.
The figure 2 shows a part of the bifurcation diagram of the system (18) 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 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.

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 : 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 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 solutionbranch where the pulsation of the oscillation verifies ω = λ 3 . 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 super-harmonic 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.
The zoom on the top of figure 3 shows the many superharmonic resonances and the symmetrybreaking (SB) bifurcations. The superharmonic resonances appear for λ 1 2p+1 , p ∈ 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,ẋ) are no longer symmetric.

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 : 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 equation (23) 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 the system (23)) is written under the assumption that −1 < p − γ < 0 which means according to [21] that the reed channel is open and air is allowed to enter the clarinet. The system (23) 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 figure 4. The values for which the system encounters a Hopf bifurcation γ H = 0.35 or a period doubling bifurcation γ P D = 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.

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 : where p is the acoustic pressure in the mouthpiece of the instrument; γ and ζ have the same meaning as in section 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 the equation (26) 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 equation (15) ; setting ∆p = p − γ, the exponential is written using A = exp(3∆p) and its differentiated form dA = 3Ad(∆p). These remarks lead to introduce the following auxiliary variables : which allows to write the equation (26) in a quadratic format : 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 τ can be expressed in term 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 [18]. Below, the value of τ is set to 2 √ 3ℓ x1 with ℓ = 17 15 and x 1 = 1 5 . ζ is set to 1.2 to draw the bifurcation diagram with respect to γ in figure 5. As said above, the equation (26) is the derivative of another 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.
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 γ, 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 γ = 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 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. is detected around γ = 0.285, that is below the theoretical value of γ = 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 .

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.