USING AUTOMATIC DIFFERENTIATION TO COMPUTE PERIODIC ORBITS OF DELAY DIFFERENTIAL EQUATIONS

. In this paper we focus on the computation of periodic solutions of Delay Diﬀerential Equations (DDEs) with constant delays. The method is based on deﬁning a Poincar´e section in a suitable functional space and looking for a ﬁxed point of the ﬂow in this section. This is done by applying a Newton method on a suitable discretisation of the section. To avoid computing and storing large matrices we use a GMRES method to solve the linear system because in this case GMRES converges very fast due to the compactness of the ﬂow of the DDE. The derivatives of the Poincar´e map are obtained in a simple way, by applying Automatic Diﬀerentiation to the numerical integration. The stability of the periodic orbit is also obtained in a very eﬃcient way by means of Arnoldi methods. The examples considered include temporal and spatial Poincar´e sections.


(Communicated by Yang Kuang)
Abstract. In this paper we focus on the computation of periodic solutions of Delay Differential Equations (DDEs) with constant delays. The method is based on defining a Poincaré section in a suitable functional space and looking for a fixed point of the flow in this section. This is done by applying a Newton method on a suitable discretisation of the section. To avoid computing and storing large matrices we use a GMRES method to solve the linear system because in this case GMRES converges very fast due to the compactness of the flow of the DDE. The derivatives of the Poincaré map are obtained in a simple way, by applying Automatic Differentiation to the numerical integration. The stability of the periodic orbit is also obtained in a very efficient way by means of Arnoldi methods. The examples considered include temporal and spatial Poincaré sections.
1. Introduction. Delay Differential Equations (DDEs) theory have been studied in the mathematical literature by many authors, such as [8,10] and many others. However, many mathematical models with constant time delays or with state delays require complicated mathematical analysis to be treated only partially. As a result, such models require numerical methods to analyse the dynamical features; steadystate points, periodic orbits, bifurcation behaviours, etc. The development of these methods is, in general, difficult due to the properties of a delay dynamical system and it may be much more complicated than the non-delay case. Consequently, research on numerical techniques for DDEs has been mainly focused not only on time integration but also in the computation of invariant objects ( [18]).
The existence of periodic solutions is an important topic of interest in many applications of DDEs. There are many results and methods to study their existence, stability and dependence on parameters of certain classes of DDEs, [20]. However, many of these results and the corresponding methods are essentially theoretical and they cannot be easily applied to a general nonlinear system with several delays.
In this paper we focus on the computation of periodic orbits and their stability. The main idea is to compute the periodic orbit as a fixed point of a suitable Poincaré map by means of the Newton method. The main novelty is that the linear system that appears at each Newton iteration is solved by a GMRES method, which do not require the computation of the full Jacobian matrix, but only its action on given vectors (directional derivatives). These derivatives are computed by means of Automatic Differentiation (AD), a technique that is very suitable to propagate derivatives along a computer code, since it is very simple to implement [6,7,15].
We recall that the periodicity problem for DDEs is an infinite-dimensional problem because such a problem is defined in an infinite-dimensional space. Indeed, while only an initial point is required to characterise periodic solutions of an Ordinary Differential Equation (ODE), computing a delay periodic orbit requires an initial function to be found. Such computation will be done using a shooting approach [14] discretising an initial function on the delay interval and using a discrete version of the Poincaré operator whose eigenvalues will be approximations of the Floquet multipliers of the periodic solution of a likely large system. They will be clustered to zero and usually only a few of them have large modulus [8].
1.1. Delay differential equations. A Delay Differential Equation (DDE) is a functional equation of the formẋ where g is a suitable function defined on a domain contained in R × C to R n , where C denotes the vector space of continuous functions from the closed interval [−1, 0] to R n which, endowed with the sup norm, is a Banach space. As usual, the function x t ∈ C is defined as 2) It is well known ( [8]) that, if g is Lipschitz, then for any initial condition x t0 ≡ u ∈ C there exists t f > t 0 and a unique solution of (1.1), x u , defined on [t 0 −1, t f ) verifying this initial condition.
It is also well known the lack of smoothness of the solution, no matter how smooth the initial data is. For instance, (1.2) can have discontinuities in the derivatives at t 0 + k being k a positive integer. However, it becomes smoother in such points when k becomes larger. More precisely, if f is of class C p and x u is of class C q (q < p) at t 0 + k, it will be of class C q+1 at t 0 + k + 1. Finally, notice that if x u is periodic, then such lack of smoothness will not appear and the periodic orbit is as smooth as f .
2. Numerical methods. There are several well-known methods for the numerical integration of DDEs (see, for instance, [3]), and some are available as public domain codes. In this section we explain how to modify an existing numerical integrator (explicit or implicit) of DDEs such that it also propagates efficiently the derivatives of the flow. As an example, we have applied it to two algorithms.
The first numerical integrator has been developed by ourselves, as a modification of the standard Runge-Kutta-Fehlberg (7)8, to deal with DDEs with a single constant delay that it can be assumed to be 1. We have used Barycentric Rational Interpolation (based on the use of Chebyshev points) for the interpolation of the solution on the interval of delay. This interpolation scheme provides convergence on the whole delay interval ( [2]) which makes it very suitable in this situation. See Appendix A.2 for more details.
The second numerical integrator is the well known retard code explained in [9] and provided in E. Hairer's website. It is based on a Runge-Kutta pair of orders 4 and 5 with step size control and dense output, in such a way that, after every successful integration step the computed solution is stored in a convenient way such that it is available in the next steps. Let us remark that this code can handle DDEs with several delays.
We have modified both codes to use Automatic Differentiation to propagate the derivatives of the flow, as discussed in the next section.
2.1. Automatic differentiation. Automatic Differentiation (AD) is a computational tool to obtain derivatives of the output of an algorithm w.r.t. initial data and/or parameters [6,7,15]. In this paper we are only interested in the computation of first order derivatives so we focus on this case (an extended discussion of the applications of AD with higher derivatives can be found in [11,5]).
To summarize AD technique for first derivatives, assume that we have a computer program that, given a data x ∈ R n , produces a result y ∈ R m , Here, Ψ denotes the implemented algorithm. Of course, we are assuming that Ψ is of class C 1 . Then, given a set of data x 0 = (x 0 1 , . . . , x 0 n ) we replace x 0 by x 0 + s = (x 0 1 + s 1 , . . . , x 0 n + s n ), where s is a vector of symbols. The coefficient of the symbol s j is the partial derivative of the actual value w.r.t x 0 j so it is 1 at the beginning. Then, we modify the arithmetic operations in Ψ such that every operation is also performed on the first derivatives. The result of the evaluation of Ψ on x 0 is so we have also produced the Jacobian of the function. This technique is easily implemented in C++ since this language allows for the replacement of the basic arithmetic types and operations by user defined ones. For this reason we have translated the original FORTRAN version of retard into C++ in order to work with our AD library.
As it has been mentioned in the Introduction, we plan to use a GMRES iterative solver. A characteristic of this method is that it requires from the user a function that, given a vector v, returns the vector Av (here A is the matrix of the linear system to be solved). In our case, this implies that we will need to compute the directional derivative of the flow (in the direction v). To do this, we will use a single symbol, say s 1 , and replace the initial condition x 0 by x 0 + s 1 v. Using the same notation as before, and assuming that Ψ denotes the flow at a time, say, t f , we have where u is the directional derivative of the flow in the direction v (in other words, u = DΨ(x 0 )v). As it has been mentioned before, we could also produce the complete Jacobian by using as many symbols as the dimension of x 0 .

2.2.
Poincaré sections. Given a dynamical system, a Poincaré section is defined as a codimension 1 smooth manifold transversal to the flow. In many cases, Poincaré sections are not global, in the sense that not all the trajectories of the system must cross the section. They are a standard tool to study the dynamics in regions of the space where the orbits have some kind of recurrence. In our case (the flow of a DDE), the Poincaré section is a functional space of "initial conditions" for the DDE.
In this paper we use Poincaré sections to decrease the complexity of the computations: instead of applying the Newton method to a discretisation of the full periodic orbit, we apply it to the piece of trajectory contained in the section. This reduces significantly the amount of computations, specially when the period is large. To simplify the presentation, let us discuss the two kind of sections we have used.
2.2.1. Spatial section. We will focus first on the most common kind of section: the one defined by a codimension one hyperplane. They are enough for most situations and, at the same time, very easy to work with. In this paper, a C p -section is a hyperplane such that l is a bounded linear operator and α ∈ R. The degree of smoothness p depends on the example at hand. As we are looking for periodic orbits, we can use the same smoothness as the DDE. Moreover, we need the flow of the DDE to be transversal to the section, at least locally. A usual transversality condition is l(ẋ) > 0. For a more detailed discussion on Poincaré sections for DDEs see, for instance, [22]. To compute the Poincaré map P , we start from some initial data in the section, we propagate the corresponding orbit till it crosses the section in the right direction and we compute the intersection point with the section. Under generic conditions, it is possible to see that the map P is continuous and compact in C p ( [22]).
As it is usual in numerical methods, the elements of the Poincaré section S have to be discretised. Here we have used Barycentric Rational Interpolation, based on Chebyshev points. Its main advantage is that it converges on the whole interval of interpolation. See Appendix A.1 for more details. So, if these functions are discretised as a table of values (t i , u i ) ∈ [−1, 0] × R n , i = 0, . . . , N , the section is given by an expression σ(u 0 , . . . , u N ) = 0, where being a and a i fixed vectors in R n defining the desired section. As usual, the integration starts from an initial condition in the section and, during the integration, the condition σ(u 0 , . . . , u N ) = 0 is checked. When σ changes sign in the right direction, the equation σ(u 0 (t), . . . , u N (t)) = 0 is solved (by means of a Newton or a secant method) to find the return time to the section, t , and the image of the initial condition by the Poincaré map.
We note that the return time t is not a multiple of the delay. This means that we have to interpolate the trajectory on the time interval [t − 1, t ] which contains a point on which some derivatives of the orbit are not smooth (see Section 1.1). This implies that the interpolating function does not need to be as accurate as we would expect. However, note that our goal is to compute periodic orbits, which are known to be as smooth as the DDE. In all examples considered, the DDE is C ∞ (or even analytic) so, although we may have some relevant interpolation error at the beginning of the iterations to compute the periodic orbits, this error will disappear when the iterations converge to the fixed point corresponding to a periodic orbit.
Let us comment on the computation of the differential of the Poincaré map. As it has been mentioned before, we plan to use GMRES to solve the linear systems that appear at each step of a suitable Newton method, so let us discuss first the computation of a directional derivative. In Section 2.1 we have discussed the propagation along the flow of a directional derivative by means of AD. Once the integration is finished because is has arrived to the Poincaré section at a time t = t (see above), the directional derivative is not contained in the tangent space of the section (in this case, as the section is an hyperplane, the tangent space is the section itself). As it is usual in this situation, the directional derivative has to be projected onto the section following the direction of the vector field of the orbit at t = t (the intersection point with the section). If the complete Jacobian is needed, it can be computed in the same way, applying the projection to each column.

2.2.2.
Temporal section. This is the standard section used for a system that depends periodically on time, with a known period ρ. Then, generally speaking, periodic orbits have a period t * which is a multiple of ρ. Hence, it is natural to use the Poincaré section given by the time (t = 0 mod ρ) and the Poincaré map is the time-ρ flow. The differential of this Poincaré map is computed by Automatic Differentiation, as explained in Section 2.1. The computations are done in the same way as in the previous section, but note that here they are simpler: there is no need to solve an equation to arrive to the section, and the directional derivative of the flow is already the directional derivative of the map (without the need of any projection).
3. Periodic orbits. We assume that we have a DDE with a periodic orbit, with a period larger than the delay. We also assume that we have a suitable Poincaré section in the sense that the periodic orbit is transversal to it. We look for the periodic orbit as a fixed point of this Poincaré section. This means that we have to solve a nonlinear equation of the form G(u) = Φ(u) − u = 0, where Φ denotes the Poincaré map. This equation is solved by Newton iteration.
3.1. The Newton method. There are two relevant ingredients in a Newton's method. The first one is to choose a suitable initial condition and the second one is to provide the differential of the function that one wants to obtain the zero. The former is usually obtained using a priory information on the model. The computation of the differential strongly depends on each situation. Clearly, the differential computation has been reduced to the computation of D u Φ. As in the ODE approach, this computation could be done by using variational equations.
Here we have chosen to use Automatic Differentiation.
As we have mentioned before, the initial condition u of a DDE is discretised as a table of values, let us say (s i , u i ), u i = u(s i ), for i = 0, . . . , N . After some integration steps, we obtain a new point, namely v = x u t * +s * . This numerical integration can be seen as a sequence of elementary operations that, given an initial table of values (s i , u i ), produces another table of values (s i , v i ). Then Automatic Differentiation (AD) can be used to compute the derivatives of the result w.r.t. the initial data, by simply overloading the arithmetic used during the numerical integration. To produce the Jacobian matrix, we add a different symbol to each component of the vectors u i , and we propagate it (to first order) along all the operations of the numerical integration. With this, we obtain the final results v i plus linear polynomials in the initial symbols. The coefficients of these polynomials give the Jacobian of the map. Note that we can use the same ideas to compute a directional derivative: given the initial direction, we add to the initial data a single symbol times the direction. Then, we propagate this symbol to obtain the final result plus the symbol times a vector. This vector is the directional derivative. Note that the computation of a directional derivative requires much less operations than to compute a full Jacobian.
3.1.1. Using GMRES. In some cases, the dimension of this Jacobian matrix can be large (several state variables with several delays, for instance). An option to deal with very large linear systems is to use an iterative solver like GMRES [17]. These solvers do not require to obtain the full matrix of the system, they only require to compute the product of this matrix by a vector. Note that this product is the directional derivative of the Poincaré map w.r.t. this vector. Hence, instead of computing the matrix, we perform a numerical integration to compute the required directional derivatives. If the spectrum of the matrix is clustered, methods such as GMRES have a fast convergence [1,4]. We note that the differential of the flow of the DDEs considered here is given by a compact operator, which guarantees that the spectrum of this differential is clustered.

3.2.
Stability. Once a delay periodic orbit, with initial condition u * and period t * , has been found by the single-shooting method explained before, its stability is determined by the spectrum of the linearised monodromy operator. According to [8], if the DDE is smooth, autonomous and t * ≥ 1, then this linearised operator is a compact operator whose spectrum has zero as its only cluster point.
Since only a few eigenvalues are important, the computation of the respective eigenvalues and eigenfunctions should be done for some of the available methods for solving large-scale eigenvalue problems in an iterative way, in a similar way as the GMRES method. In this direction, the ARnoldi PACKage, also called ARPACK [13], is capable of solving large-scale Hermitian, non-Hermitian, standard or generalised eigenvalue problems. It is designed to compute a few, say k, eigenvalues with user-specified features using M · O(k) + O(k 2 ) storage and no auxiliary extra storage required. This software is quite standard in numerical computation and it is also well-known that it is based upon an algorithmic variant of the Arnoldi process called the Implicitly Restarted Arnoldi Method (IRAM) [21,13].
3.3. Computer implementation. As we have already explained, Automatic Differentiation requires the use of a new arithmetic. Numerically, it can be expressed as a new data type, that is, we consider a pair (x, x ) instead of just the value x where here x denotes the derivative that is being propagated. Hence in our approach if the initial condition u of a DDE is discretised as a table of values, let us say u = (u 0 , . . . , u N ) with s i the abscissae and u i = u(s i ) for i = 0, . . . , N . To each u i we associate a collection of symbols u . Thus, the initial data is (u, u ) with u = (u 0 , . . . , u N ). After some integration steps, we obtain a new point, namely v = x u t * +s * which is discretised in v at the same abscissae and the propagated values v w.r.t. the initial data (u, u ). This numerical integration can be seen as a sequence of elementary operations that, given an initial condition (u, u ), produces another value (v, v ). The size of our execution is larger because of this data type. Indeed, u i is an n-dimensional vector and u i has either dimension n in the directional case or n(N + 1) in the Jacobian case. It is clear that this technique does not depend on the integrator of the DDE as long as it accepts the new data type and its arithmetic. Notice that, in fact, we are computing the derivative of the algorithm that computes (at the same the execution) the Poincaré map and one must be aware that if a spatial section is used, v must be projected to the section.

Numerical examples.
As first test, we use the model studied in [14]. The equation is defined bẏ Clearly, x ≡ 0 is an equilibrium solution. We are going to illustrate the computation of two different cases. The first one consists in a temporal section and the second one in a spatial section.

4.1.
A first example with temporal Poincaré section. Let us consider a periodic perturbation of period 2π, such as, The Implicit Function Theorem in Banach spaces tells us that there is a periodic orbit of period 2π close to ε = 0 and x ≡ 0. Thus, a temporal section of time 2π leads to a delayed stroboscopic mapping. The solution is continued using a pseudoarc-parameter method with respect to the three parameters in (4.2), the parameter α and ε and the delay τ . In this last case, the (4.2) is modified by a change of t, that is, Because of the change of time, the temporal section is also suitably modified with the value of 2π/τ . The initial values of these parameters has been α = 1.5, ε = 10 −4 and τ = 1 with a discretisation size N = 32.

4.2.
An example with a spatial Poincaré section. Studying the stability of the equilibrium point of the system (4.1) with respect to the parameter, one shows a family of periodic solutions bifurcates at α = π 2 because of a Hopf bifurcation. These periodic orbits have an unknown period and a spatial section will allow us to compute one of them. If the initial condition is discretised by a table of values (s i , u i ), i = 0, . . . , N . The section is given by an affine hyperplane, let us say, The values a and a i must be fixed at each step of the continuation in such a way that the section must be transverse to the flow. We have selected the easiest one  In the left hand side, the orbit is displayed with the initial condition in −1 ≤ t ≤ 0 and final lag-segment once the second has been crossed two times. The phase space of the periodic orbit is shown in the right hand side.
σ(x 0 , . . . , x N ) = x 0 and α close to the bifurcation point, e.g. α = 1.57. Figure 4.2 shows the periodic orbit around the solution x ≡ 0 and α = 1.57. Such a computation has been done with a Newton tolerance of 10 −12 , a Section tolerance 10 −14 and an integrator tolerance of 10 −15 . Once an initial periodic orbit has been computed for α = 1.57 a pseudo-arc-length method ( [19]) can be done to perform a continuation with respect to the parameter α, see

4.3.
An example with several delays. The next example illustrates that our method is independent of the model whenever one has a suitable integrator. In particular, it can be applied for several delays, such as, the equatioṅ  x'(t) x(t) It has been studied extensively [16,12] and existence of periodic orbits has been proved for some values of the parameters. We take concrete values: λ 1 = λ 2 = 2.5, λ 3 = 0.25, τ 1 = 1.65, τ 2 = 0.35 and τ 3 = 1 such that the existence of a periodic orbit is guaranteed. The initial condition u at time t 0 = 0 is the integration of the orbit through cos(t) at t 0 = 0 until a final time of 4. This initial condition u is discretised with N = 64, i.e. 65 points, and (s i , u i ) with s i Chebyshev abscissae of second type and u i = u(s i ). We have applied the spatial section x N = u N to detect a periodic orbit after 4 crosses using the Hairer's integrator [9] with a tolerance 10 −14 . In this case, we have executed the two approaches of AD, that is, a full Jacobian implementation, that uses an LU solver, and a Jacobian-free implementation that uses a GMRES solver. The maximum dimension of the Krylov subspace has been 5  for each of the six different linear systems to be solved during the Newton's method with a tolerance of 10 −12 . The use of the iterative solver allows us to consider 10 −3 b 2 as stopping criterion of the GMRES of the linear system Ax = b. The CPU time to compute the periodic orbit is 4 seconds in the case of the full Jacobian while the Jacobian-free only needs 0.73 seconds. The Figure 4.6 shows the orbit which is stable after computing the three first eigenvalues.

5.
Conclusions and future work. This paper illustrates, in essence, that the same techniques applied in ordinary differential equations to compute a periodic orbit and to perform a continuation with respect to parameters can be done in a delay differential equation context. Although these techniques are almost the same, the main contribution is how the differential required by the Newton's method is computed. This is done by changing (enlarging) the standard double precision type to a new data type that also contains derivatives (which are also double precision types). The computer arithmetic is then enlarged to operate on this new data type so that it propagates not only values but also derivatives along the algorithms (this is done quite easily in languages like C++). The advantage of this point of view is that these techniques can be applied very easily to any numerical integrator. We note that, to compute a periodic orbit, we only need to discretise a segment of the orbit (of the size of the delay), which implies that the number of unknowns of the linear system is smaller than the number of unknowns obtained when discretising the full orbit. Moreover, when the problem becomes more complex and the Newton method leads to high dimensional systems, we have shown that iterative methods such as GMRES work very effectively (see Section 4.3). Finally, we have seen that the stability of the orbits can also be studied efficiently by means of iterative methods such as those based on Arnoldi's schemes.
To illustrate these methods, we have considered several DDEs models with families of periodic orbits and Hopf bifurcations. We have continued these families and computed their stability. In the near future we expect to extend the methods being used for quasi-periodic solutions in ODEs to DDEs by means of similar techniques.
Appendix A. Another delay time integrator. Let us assume a Delay Differential Equation (DDE) with constant delay. It can be seen as a chain of Ordinary Differential Equations (ODEs). Therefore a reasonable first approach is to consider a standard integrator and whenever an unknown value was required, it would be interpolated. However, such an interpolation should be as far accurate as the integrator is. We propose a Runge-Kutta-Fehlberg as ODE integrator and barycentric rational interpolation.
A.1. Barycentric rational approximation. Let ϕ : [a, b] → C be a function, φ be the affinity defined by Therefore it is enough to store the first N 2 Chebyshev points and the N function values.
On the other hand, if x k are the Chebyshev points of second type defined by N.B.: II 0 > · · · > II N and II k = −II N −k for any 0 ≤ k < N 2 . (A.2) becomes R[ϕ; II](y) So it is enough to store the first N 2 + 1 Chebyshev points and the N + 1 function values. In both cases, R[ϕ; I] and R[ϕ; II] does not have any pole in [−1, 1]. Furthermore, under analytic hypothesis one proves that the error with respect to the exact value is O(ρ −N ) for some ρ > 0, [2]. The error behaviour of each of them are really similar although R[ϕ; I] has an error more similar to that computed by the classical Chebyshev approximation. Hence, one may think that it does not matter whether I or II are used but each of them have their advantages. For instance, • I does not contain the values −1 and 1.
• II contains the values −1 and 1.
• II aN ak = II N k for any 0 ≤ k ≤ N and positive integer a. • If N is even, then II N/2 = 0. Notice also that if ϕ gives values on an n-dimensional space, one can apply the same scheme in each of the n coordinates so the complexity will be Θ(nN ).
A.2. Delay-Runge-Kutta-Fehlberg. The scheme stores the values at the Chebyshev points in the lag interval I k = [t 0 + k, t 0 + (k + 1)], it computes and stores the values at the Chebyshev points of the next lag interval I k+1 = [t 0 +(k+1), t 0 +(k+2)] and the evaluation of the DDE is done by the Chebyshev approximation or the Barycentric Rational Approximation on I k . Let us assume that the integration step is performed by a Runge-Kutta-Fehlberg method, namely rkf. The Algorithm 1 sketches the method. Note that the main idea is that the step h, which is modified during the integration, must be bounded with a maximum step whose value is just the difference between Chebyshev points, the last point of the interval or the final time.