A LEGENDRE-GAUSS COLLOCATION METHOD FOR NONLINEAR DELAY DIFFERENTIAL EQUATIONS

. In this paper, we introduce an eﬃcient Legendre-Gauss collocation method for solving nonlinear delay diﬀerential equations with variable delay. We analyze the convergence of the single-step and multi-domain versions of the proposed method, and show that the scheme enjoys high order accuracy and can be implemented in a stable and eﬃcient manner. We also make numerical comparison with other methods.

1. Introduction.Delay differential equations (DDEs) naturally arise in diverse science and engineering applications [25], and in recent years, a large body of literatures has been devoted to numerical solutions of DDEs (see, e.g., [3]- [6] and [9,10,14,26]).Among the existing methods, numerical schemes based on Taylor's expansions or quadrature formulas have been frequently used (cf.[11,12,25,26,28] and the references therein), e.g., the implicit Runge-Kutta methods, which can be systematically designed and often provide accurate approximations.Over the years, spectral method has become increasingly popular and been widely used in spatial discretizations of PDEs owing to its high order of accuracy (cf.[7,8,13,16,17,18]).The solution of a DDE globally depends on its history due to the delay variable, so a global spectral method could be a good candidate for numerical DDEs.Some work has been done along this line, and we particularly point out that Ito, Tran and Manitius [27] proposed and analyzed a Legendre-tau method for linear DDEs with one constant delay, where the solution was approximated by a truncated Legendre series, and the unknown expansion coefficients was solved from a Lanczostau formulation.Moreover, Guo and Wang [21], and Guo and Yan [23] developed the Legendre-Gauss collocation methods for ODEs based on Legendre polynomial expansions.Meanwhile, Guo and Wang [20], and Guo et al. [22] discussed the Laguerre-Gauss-type collocation methods for ODEs.However, it is more interesting but challenging to develop and analyze such type of high-order methods for nonlinear DDEs with variable delay of the form where W (t) = U (t − θ(t)), f, θ and V are given functions and the delay variable: θ(t) ≥ 0. We start with a single-step scheme motivated by [21].Basically, we approximate the solution by a finite Legendre series, and collocate the numerical scheme at Legendre-Gauss points to determine the coefficients.Due to the delay variable, it's inevitable to solve nonlinear systems likewise for the implicit Runge-Kutta methods, but the iterative solvers can be implemented efficiently thanks to the fast transform [2].The scheme has an infinite order of accuracy both in time and the delay variable, and is particularly attractive for DDEs with highly oscillatory solutions and/or stiff behavior.Therefore, it enjoys some remarkable advantages over the Runge-Kutta-type methods.
For a more effective implementation, we suggest a multi-domain scheme to advance in time subinterval by subinterval due to the following twofold considerations.Firstly, the resultant system for the expansion coefficients can be solved more stably and efficiently for a small or modest number of unknowns.In particular, for large T, it is desirable to partition the solution interval (0, T ) and solve the subsystems successively.We also want to make a point that we merely need to store the Legendre coefficients of the solution in relevant "delay" subintervals to recover the solution in the subinterval that needs to resolve.Hence, the scheme can be implemented efficiently and economically.On the other hand, the multi-domain scheme provides us a flexibility to handle DDEs involving non-smooth initial data and/or solutions.In such cases, the partition of the interval can be adapted to the evolution process as the adaptive Runge-Kutta methods, and we are able to place more points in the subintervals that are needed.
The rest of the paper is organized as follows.In the next section, we present and analyze the single-step Legendre-Gauss collocation method, and provide some numerical results to justify our theoretical analysis.The multi-domain version is described, and the convergence result is also derived and numerically illustrated in Section 3. The final section is for some concluding discussions.
2. Single-step Legendre-Gauss collocation method.In this section, we describe and analyze a single-step numerical integration process for the DDE (1) using the Legendre-Gauss interpolation, which serves as a base for the multi-domain scheme to be presented in the forthcoming section.
2.1.Preliminaries.Let L l (•) be the Legendre polynomial of degree l defined on [−1, 1].The shifted Legendre polynomials L T,l (t), t ∈ [0, T ], are defined by (cf.[21]) According to the properties of the standard Legendre polynomials, we have Moreover, there holds the recursive relation The set of L T,l (t) forms a complete L 2 (0, T )−orthogonal system, namely, where δ l,m is the Kronecker symbol.Thus for any v ∈ L 2 (0, T ), we write We now turn to the Legendre-Gauss interpolation.Let {t N j , ω N j } N j=0 be the Legendre-Gauss quadrature nodes (in (−1, 1) and arranged in ascending order) and weights.Accordingly, we define the shifted Legendre-Gauss interpolation nodes and weights as Let P N (0, T ) be the set of all real polynomials of degree at most N.We recall that for any φ ∈ P 2N +1 (0, T ), ω N T,j φ(t N T,j ).
(6) Denote by (u, v) T and v T the inner product and the norm of the space L 2 (0, T ), respectively.Define the following discrete inner product and norm associated with the Legendre-Gauss quadrature: T,N .
2.2.The single-step scheme.We next present the single-step scheme for the delay differential equation (1).For this purpose, we denote the grid set by Λ N := t N T,k : 0 ≤ k ≤ N ⊂ (0, T ).The single-step Legendre-Gauss collocation approximation to (1) is to find u N (t) ∈ P N +1 (0, T ), such that where the delay term Here, we recall that θ(t) ≥ 0 and V (•) is a known function.Denote by , and set w N (t) = u N (t − θ(t)).The collocation scheme ( 14)-( 15) can be reformulated as: Find u N (t) ∈ P N +1 (0, T ) such that supplemented with the initial condition u N (0) = V (0).We notice that it is an implicit scheme.
An important problem is how to resolve (16).Indeed, we may follow Lambert [28] to design an algorithm (mainly for ODEs) to resolve the discrete system (16) based on the Lagrange interpolation.However, it is known that Lagrange interpolation is not stable for large N. Hence, we propose a stable approach by expanding u N (t) directly in terms of the shifted Legendre polynomials and determining the unknown coefficients of the collocation scheme (16).This approach is stable for large N, and much easier to implement particularly for DDEs, since one only needs to store the coefficients of the numerical solution at each step.One may refer to Remark 2.1 of [21] for some more features of this method for ODEs.We next describe the numerical implementation of the collocation scheme (16).For this purpose, we expand the collocation solution as Furthermore, let [l] be the integer part of l.According to [21], we have On the other hand, L T,l (0) = (−1) l .Hence, ( 16) is equivalent to where and and A N T be the matrix with the entries a N T,k,l , 0 ≤ k ≤ N, 1 ≤ l ≤ N + 1.Then we can rewrite (18) into the following matrix form: In actual computations, we solve u N T,l N +1 l=0 from (19), and recover the collocation solution u N (t), 0 < t ≤ T from (17).
We tabulate in Table 1 the condition numbers of A N T with T = 1 and various N , which indicates that the condition numbers grow like N 2 /4 as with the normal collocation scheme using Lagrangian basis.

Error analysis.
In this subsection, we analyze the convergence of the scheme (16).In particular, we shall prove the spectral accuracy of numerical solution u N (t).Let I T,N be the Legendre-Gauss interpolation operator as defined before.Denote Lemma 2.1.Let U and u N be respectively solutions of (1) and (16). where Proof.Let Then we have from (1) that Subtracting ( 22) from ( 16) yields Clearly, t −1 (E N (t) − E N (0)) ∈ P N (0, T ).Thereby, by (7) and integration by parts, we deduce that On the other hand, we have from ( 23) that where Since t −1 (E N (t) − E N (0)) ∈ P N (0, T ) and G N T,2 (t) ∈ P N (0, T ), we use (7) to obtain that for any ǫ > 0, T .Inserting the above estimate and ( 24) into (25) gives that We next estimate G N T,2 T .Clearly, by (12) with dU dt and r − 1 instead of u and r, we have that for integer r ≥ 2, Moreover, by (13), Therefore On the other hand, by (7), Substituting the above and ( 29) into (26), we obtain (20).
We now consider several typical f and analyze the numerical errors.Hereafter, β denotes any positive number less than 1  4 .Case I. Consider ( 16) with the linear delay: Assume that f (x, y, t) satisfies the following Lipschitz conditions in x and y.That is, there exists real numbers γ 1 ≥ 0 and γ 2 ≥ 0 such that and Theorem 2.2.If the conditions (30)-(32) hold, U ∈ H r (0, T ), with integer 2 ≤ r ≤ N + 1, and for certain δ > 0, In particular max where c β is a positive constant depending only on β.
Remark 2.1.The condition (33) is necessary for the proof, but it should not be essential.In fact, some numerical examples do not meet this condition, but the numerical scheme still converges.This remark also applies to multiple-domain cases.
Case II.Assume that the delay function satisfies: Moveover, f (x, y, t) satisfies the Lipschitz condition (31).
In particular, if 1 N ≤ T ≤ 1, then we may take r = N + 1 in (60) and (61), to reach that 2.4.Numerical results.In this subsection, we present some numerical results to illustrate the efficiency of our single-step algorithm.
Linear variable delay (Case I): As pointed out in [15], the exact solution is u(t) = e t .Obviously, the conditions (31) and (32) hold with γ 1 = 1 2 and γ 2 = 1 2 e T 2 .Moreover, the inequality (33) is satisfied for T = 0.3.But for T = 1, the inequality (33) is no longer valid.In Fig. 1, we plot the numerical errors at t = T with T = 0.3, 1 and various values of N.They indicate that the numerical errors decay exponentially as N increases.In particular, we can observe from Fig. 1 that even if the condition (33) is not satisfied (for instance, T=1), our algorithm is still valid.
This is a nonlinear DDE with the exact solution u(t) = sin(t), see [15].In Fig. 2, we plot the numerical errors at t = 1 using the single-step scheme (16) and Newton-Raphson iteration method with various values of N.They indicate that the numerical errors decay exponentially as N increases.We also note that the condition (32) is not satisfied for problem (64), but our algorithm is still valid.Linear constant delay (Case II): The exact solution is u(t) = sin(t).Obviously, the condition (31) hold with γ 1 = 0.Moreover, the inequality (51) is satisfied.In Fig. 3, we plot the numerical errors at t = π 2 with various values of N.They indicate that the numerical errors decay exponentially as N increases.3. Multiple-domain Legendre-Gauss collocation method.In the last section, we investigated the single-step Legendre-Gauss collocation method.The numerical errors decay very rapidly as N and r increase.While the foregoing singlestep collocation method provide accurate results, in actual computation, it is not convenient to resolve the discrete system (19) with very large mode N. As commented in the introduction, it is advisable to partition the interval (0, T ) into a finite number of subintervals and solve the equations subsequently on each subinterval.Firstly, replacing T and N by τ 1 and N 1 in ( 16) and all other formulas in Subsection 2.2, we can derive an alternative algorithm, with which we obtain the numerical solution u N1 1 (t) ∈ P N1+1 (0, τ 1 ).Then we evaluate the numerical solutions u Nm m (t) ∈ P Nm+1 (0, τ m ), 2 ≤ m ≤ M, step by step.Finally, the global numerical solution of ( 1) is given by We now present the numerical scheme for u Nm m (t).Denote by t Nm τm,k and ω Nm τm,k , 0 ≤ k ≤ N m the nodes and the corresponding Christoffel numbers of the shifted Legendre-Gauss interpolation on the interval (0, τ m ).Let where where ).We see from (67) and (68) that the local numerical solution u Nm m (t) is actually an approximation to the local exact solution U m (t), with the approximate initial data u Nm m (0 Error analysis.We next analyze the numerical errors.Denote Let U m and u Nm m be respectively solutions of (68) and (67). where Then, following the same line as in the derivation of ( 23), we obtain from ( 67) and (68 Like ( 24), we have Next, we use (70) to derive that , j ≥ 0. Hence, by (7) with τ m and N m instead of T and N , we deduce from the above inequality that According to (29) with the parameters N m and τ m , we have A combination of ( 71)-(73) leads to (69).
We now consider several typical f, and analyze the numerical errors.Hereafter, let τ = max 1≤j≤M τ j and N = min 1≤j≤M N j .Moreover, we always assume that for any 1 ≤ i ≤ j ≤ M, τi τj is bounded above.Case I. Consider (67) with the linear delay: Assume that f (x, y, t) satisfies the Lipschitz conditions (31) and (32).

Numerical results.
In this subsection, we present some numerical results to illustrate the efficiency of our multiple-domain algorithm.
Linear variable delay (Case I).
We consider the example (63) with T = 1 and uniform step-size grid.As pointed out before, the conditions (31) and (32) hold with γ 1 = 1 2 and γ 2 = 1 2 e 1 2 .Furthermore, the inequality (75) is satisfied for τ m ≡ 0.25, but for τ m ≡ 0.5, the inequality (75) is no longer valid.Moreover, (1 − λ)T m ≤ T m−1 with λ = 1 2 .In Fig. 4, we plot the numerical errors at t = 1 with uniform τ m and N m .They indicate that the numerical errors decay exponentially as N m increases and τ m decreases.In particular, we can observe from Fig. 4 that even if the condition (75) is not satisfied, our multiple-domain algorithm is still valid (see the case τ m ≡ 0.5).

Nonlinear variable delay (Case I).
We consider the example (64).In Fig. 5, we plot the numerical errors at t = 1 using the multiple-domain scheme (67) and Newton-Raphson iteration method with uniform τ m and N m .They indicate that the numerical errors decay exponentially as N m increases and τ m decreases.Again, we observe from Fig. 5 that our algorithm is still valid even if the condition (32) is not satisfied.Linear variable delay (Case II).
), for t > 0,     As pointed out in [24], the exact solution is The first derivative of the solution is not continuous at t = −0.5 and t = 0, therefore the second derivative also has a jump at the points where the time lag is equal to −0.5 and 0, for instance, t = ξ 1 = 1 and t = ξ 2 = √ 2. Denote by n 0 = 2k the total number of steps.Obviously, the conditions (31) and (32) hold with γ 1 = 0 and γ 2 = 1.Furthermore, by choosing suitably T m , we can ensure that the condition (91) holds.In Fig. 6, we plot the numerical errors at As pointed out, in actual computation, even if the conditions for Theorems 3.1 and 3.2 are not satisfied, the numerical solutions of our method match the exact solutions very well.We next present some other numerical examples to illustrate the efficiency of our method, which can be found in [5].
Nonlinear constant delay (Case II): This is a well-known equation from biology, solution of which has the breaking points at the integers 0, 1, • • • .
In Table 1, we compare the errors of our method LGC3(5) (N = 3, 5 and the step-size is constant) with the method UCIRK4(6) (the uniformly corrected implicit Runge-Kutta method of order 4(6) presented in [5]) and the method IRKSR4(6) (the implicit Runge-Kutta superconvergence rate 4(6) presented in [4]), at the point t = 20 with respect to the reference solution u(20) = 4.671437497500.Since the three methods have the same mesh and the same number of interpolation nodes, we can observe that our method provides more accurate numerical results.In particular, our method is still valid even for large time step-size.In Fig. 7, we plot the numerical errors at t = 20 with τ m = 1, 1 ≤ m ≤ 20 and uniform N m .They indicate that the numerical errors decay exponentially as N m increases.(95) The solution has the breaking points: ξ 0 = 0, ξ 1 = 2.1461932206205825852, ξ 2 = 4.9254498245082464926, etc..
In Table 2, we compare the numerical errors of various methods at the point t = ξ 2 with respect to the reference solution: u(ξ 2 ) = 76.3734726693768056269.
Since the three methods have the same mesh and the same number of interpolation nodes, we can observe that our method provides again much better numerical results.
In particular, our method is still valid even for large time step-size.In Fig. 8, we plot the numerical errors at t = ξ 2 with τ 1 = ξ 1 , τ 2 = ξ 2 − ξ 1 and uniform N m .They indicate that the numerical errors decay exponentially as N m increases.
4. Concluding discussions.In this paper, we proposed the single-step and multiple-domain Legendre-Gauss collocation integration processes for nonlinear DDEs.These approaches have several attractive features.• The single-step Legendre-Gauss collocation method is easy to be implemented for nonlinear DDEs, and possesses the spectral accuracy.It enjoys computational efficiency and accuracy over some popular existing methods.• By using the multiple-domain Legendre-Gauss collocation method, we can use moderate mode N to evaluate the numerical solutions more stably and effectively.Moreover, it can be adapted to solutions with jump-discontinuities.For any fixed mode N, the numerical errors decay as τ → 0. For any fixed step size in time, the numerical error decays very rapidly as N increases.The reason of this phenomena might be that there exists the factor N −r in the error estimations of our new method, whereas there is no such factor in the error estimations of the corresponding implicit Runge-Kutta method.In fact, in the derivation of algorithm of our method, we benefit from the orthogonality of Legendre polynomials, while in the derivation of algorithm of the implicit Runge-Kutta method, one used the Lagrange interpolation on the Legendre-Gauss interpolation nodes, which is not stable for large N .In particular, our methods are much easier to be implemented than the implicit Runge-Kutta method for DDEs, since we only need to save the coefficients of numerical solutions in each step.• The numerical errors of our methods are characterized by the semi-norms of exact solutions in certain Sobolev spaces.These sharp norms are in particular necessary for problems with degenerated initial data.
The numerical results demonstrated the spectral accuracy of proposed algorithms and coincided with the theoretical analysis very well.

Remark 2 . 2 . 3 2
We see from (34), (35), (53) and (54) that the errors |U (T ) − u N (T )| and U −u N T decay rapidly as N and r increase.The convergence rate is O(N −r ).Thus, the smoother the exact solution, the smaller the numerical errors.In other words, the scheme (16) possesses the spectral accuracy.Remark 2.3.If d r U dt r ∈ L ∞ (0, T ), N ≥ r − 1 and T ≤ 1, then we use Theorems 2.1 and 2.2 to deduce that

Figure 1 .Figure 2 .
Figure 1.The numerical errors of Legendre-Gauss collocation method at t = T .

3. 1 .
The multiple-domain scheme.We now describe the multiple-domain scheme.Let M and N m , 1 ≤ m ≤ M be any positive integers.We first decompose the interval (0, T ] into M subintervals (T m−1 , T m ], 1 ≤ m ≤ M,such that the set of T m includes all breaking points, where T 0 = 0 and T M = T. Denote by τ m = T m −T m−1 , 1 ≤ m ≤ M. We shall use u Nm m (t) ∈ P Nm+1 (0, τ m ) to approximate the solution U in the subinterval (T m−1 , T m ].

Figure 4 .
Figure 4.The numerical errors of Legendre-Gauss collocation method at t = 1.

Figure 5 .
Figure 5.The numerical errors of Legendre-Gauss collocation method at t = 1.

Figure 7 .
Figure 7.The numerical errors of Legendre-Gauss collocation method at t = 20.
and uniform N m .They indicate that the numerical errors decay exponentially as N m increases and τ m decreases.

Table 1 .
The condition numbers with T = 1.