Continuous Galerkin methods on quasi-geometric meshes for delay differential equations of pantograph type

We analyze the optimal global and local convergence properties of 
 continuous Galerkin (CG) solutions on quasi-geometric meshes for delay differential equations with 
proportional delay. It is shown that with this type of meshes the attainable order of nodal superconvergence 
 of CG solutions is higher than of the one for uniform meshes. The theoretical 
results are illustrated by a broad range of numerical examples.

The optimal order of convergence of DG methods on uniform meshes for pantograph-type DDEs (1) have been studied in [7]. There it is shown the attainable order of the DG solution U in the space S (−1) m (J h ) of (discontinuous) piecewise polynomial of degree m ≥ 0 is O(h m+1 ), and that the optimal order of superconvergence at the mesh points cannot exceed O(h m+2 ) when m ≥ 1. This is in sharp contrast to DG solutions for ODEs [10] and for DDEs of constant delay [15], where the optimal orders of superconvergence at the mesh points are O(h 2m+2 ) and O(h 2m+1 ), respectively.
For DDEs with proportional delay, Xu, Huang and Chen [19] located all the superconvergence points of CG solutions based on uniform meshes. Their analysis is based on the supercloseness between the CG solution and the interpolant of the true solution.
The fact that for 0 < q < 1 the image of a subinterval (t n−1 , t n ) under the mapping θ(t) = qt will in general lie in the union (t i−1 , t i )∪[t i , t i+1 ) (i < n−2) of two adjacent previous subintervals makes the theoretical analysis and the computation of numerical solutions rather complex. In order to avoid this problem, different kinds of strategies have been proposed. A very prominent one consists in replacing the uniform meshes by a so-called "quasi-geometric meshes" (see for example [1,2,3,6]).
Thus, there arises the question as to what the orders of global convergence and local superconvergence of CG solutions will be under such quasi-geometric meshes. It is the aim of this paper to show that when employing certain quasi-geometric meshes a higher order of convergence accuracy can be attained for CG solutions at the nodal points.
The outline of the paper is as follows. In Section 2, we introduce the CG method for (1) on quasi-geometric meshes and discuss the existence and uniqueness of the CG solution. The main results on the optimal order of global convergence and local superconvergence of the CG solution are stated in Section 3. Section 4 shows that form of the discretized CG method and the collocation method are identical for delay differential equations with proportional delay. In Section 5, we provide numerical experiments to illustrate our theory. Possible extensions and future research work are described in the final Section 6.
2. The continuous Galerkin method for delay differential equations with proportional delay. In this section, we introduce the CG method for DDE (1) under quasi-geometric mesh. Then we give the computational form, and the existence and uniqueness of the CG equation. We assume that the given functions a, b and f in (1) are continuous on J.
2.1. The CG method. Suppose that on a given (small) initial subinterval J 0 = [0, t 0 ] of [0, T ], t 0 = q k T for a suitable value of k, the approximation φ(t) of the exact solution u is known. For instance, φ(t) can be obtained by the CG method or by the truncation of the Taylor expansion of the exact solution u(t). Subsequently, we solve the following equation:

2.2.
The computational form of the CG equation. After selecting the basis functions, the CG equation (5) can be rewritten as a system of linear algebraic equations for the vector U n = (u n,1 , · · · , u n,m+1 ) T .
The structure of these N systems changes for each value of n as we pass the following phase 1 to phase 2. To make this more precise, we first define g 0 = (L 1 (1), · · · , L m+1 (1)) , and introduce the matrices reflecting contributions corresponding to the nondelay terms in CG equations. The contributions of the delay term will introduced below with two distinct phases.
• Phase 1: 1 ≤ n ≤ l. In this phase, for t ∈ J n , the images qt < qt l = t l−l = t 0 lies in [qt 0 , t 0 ]. We let U n is given by the solution of the linear algebraic system: • Phase 2: l < n ≤ N . In this phase, for any t ∈ J n , the images qt are in J n−l . We define U n is given by the solution of the linear algebraic system: 2.3. Existence and uniqueness of the CG equation. We will briefly discuss the existence and uniqueness of the CG solution (defined by the solutions of the linear algebraic systems (8), (9)). Proof. Noting that the left-hand sides of the linear algebraic systems (8) and (9) are allÃ + h nÃn , and the given functions a(t) and b(t) are in C(J), we just need to show thatÃ is nonsingular.
Suppose that the basis functions on [−1, 1] are given by (see subsection 3.1 for details) From the properties of Legendre polynomials, we havẽ which meansÃ is nonsingular. Therefore, for any q ∈ (0, 1), there exists a positivẽ h, so that for all h ∈ (0,h) and 1 ≤ n ≤ N , (5) defines a unique CG solution 3. Global convergence and local superconvergence of the CG solution. In this section, we first list some preliminaries, then give the global convergence and local superconvergence of the CG solution, respectively.
3.1. Preliminaries. The follow preliminaries are given in this subsection which will play a key role in the proofs of the following Theorems 3.2 and 3.4.
We then introduce the Legendre polynomials on the interval [−1, 1]. They are defined by P 0 (s), P 1 (s) = s, · · · , P n (s) = 1 2 n n! d n ds n (s 2 − 1) n , n ≥ 2. The Legendre polynomials satisfy the following orthogonality relation Define M 0 (s) = 1, M 1 (s) = s and M k+1 (s) = s −1 P k (t)dt (k = 1, 2, · · · ), we get another family of polynomials, which we call the M-type series: The M-type series has the following properties: The zeros of M n+1 (s) are called the Gauss-Lobatto points of degree n + 1.

3.2.
Global error estimate of the CG solution.
Theorem 3.2. Assume the following.
is the exact solution of the initial value problem for the pantograph DDE (2). 3. U ∈ S Then for u ∈ W m+1,∞ ([t 0 , T ]), the following optimal global convergence estimate holds: Proof. For t ∈ J n , we take When j = 1, d n,1 = 1 2 (u(t n ) − u(t n−1 )). Integration u ′ (t n−1/2 + h n s) from −1 to s, we get an M-type series as following It is easy to see that For (11), we use multi-integration by parts and have (see [9] for details) The following idea is motivated by [9]. We construct a new m-degree polynomial interpolation u I of the exact solution u in J n , We let u I (t) = u(t), qt 0 ≤ t ≤ t 0 , and d * n,j is to be determined below. It is obvious that Then, the remainder term of J n is: From the formula (7) we see that

QIUMEI HUANG, XIUXIU XU AND HERMANN BRUNNER
In (15) and (16), we let the test functions η(t) and use the transformation t = t n−1/2 + h n s to get where By the orthogonality properties of P i (s) (i = 0, 1, · · · , m) and M j+1 (s) = 1 Thus, we obtain the following estimates: In order to determine the coefficients d * n,j in (17) and (18), we require that d * n,j satisfy the following equations: We then show that To prove (25), we proceed as follows: Step 1. We prove that We combine (12) and (22) to obtain Here, it is easy to see from (19), (20) and (21) Step 2. We show that Combining (23) and (26), we have This leads to When h n is sufficiently small, the coefficient matrices of the systems (23) are diagonally dominant. Thus, (27) holds.

QIUMEI HUANG, XIUXIU XU AND HERMANN BRUNNER
For n = l + 1, it is easy to show that the right-hand side term of (24) satisfies Similarly to the proof of (27) we find that Assume that (28) is valid for n = l + 2, · · · , N − 1. Considering the case n = N , we have This leads to Therefore, we have So (25) holds. The above arguments show that B n (R, η) can be expressed as By the inverse estimate |β 0 | ≤ C 1 −1 |η(s)|ds (see [9] for details), (29) and (30), there exists a constant C > 0 such that Now, we let From (7), we have When t 0 ≤ t ≤ t l , i.e., 1 ≤ i ≤ l, we have Ji |η(t)|dt u m+1,∞ .
Summing from 1 to n for J i and taking η(t) = θ ′ (t), we obtain Where we use the ab ≤ εa 2 + b 2 4ε andā := max t∈[t0,tn] a(t). Choose ε < 1 2 and omitting the first term of the right-hand side, there yields Here and in the following analysis, C denotes a generic constant that may assume different values at different places.When t ≥ t l , i.e., l < i ≤ N
By the inverse property θ L∞(In) ≤ Ch From the interpolation error estimates and (13) and (32), we see that (10) is true.
Observe that the maximum stepsize of the partition J h for [t 0 , T ] is attained in the last interval I k , where the distance is If we want the method described above to converge on the original interval [0, T ], we must also consider how to choose t 0 . It is suggestive that t 0 is chosen as to depend on l, such that that is, is the minimum integer for which (36) holds. Hence, the following convergence theorem holds.  is true.

Local superconvergence of the CG solution.
In this subsection, we get the superconvergence at the mesh points and other character points under the quasigeometric mesh.
The following theorem gives the result of the local superconvergence.

Remark 1.
We can also constant the auxiliary problem of the original problem to prove these results. In this framework, the superconvergence of the nodal points will be given and then the supercloseness between U and some interpolation Π h u of u will be obtained, and the superconvergence points are located by combining the two above results. We leave this to another paper.
The optimal superconvergence analysis can also be extended to DDEs with nonlinear vanishing delays θ(t).
3.4. Delay differential equations with nonlinear vanishing delay. We assume that the delay item θ is no longer qt, but more general delay function which is subject to the following conditions 1. θ(0) = 0 and θ(t) < t for t > 0, 2. min t∈J θ ′ (t) =: q 0 > 0. The delay differential equation (1) becomes It is easy to verify that the same results are true for the CG solution in S (0) m (J h ), where the quasi-geometric mesh is replaced by a quasi-graded mesh, where the ξ µ are now given by and the nodes t 0 = ξ 0 < t 1 < · · · < t l = ξ 1 < · · · < t 2l = ξ 2 < · · · < t (k−1)l = ξ k−1 are defined recursively by t n−l = θ(t n ), n = l, · · · , kl, with the last l + 1 nodes being assigned arbitrarily (see [2] for details). The computational form of CG solutions will require an obvious modification. We leave the details to the reader.
Example 1. We first use the CG method to solve the following DDE with proportional delay: Its exact solution is u(t) = sin(t) for any 0 < q < 1.
In Tables 1, 2 and 3 we present the errors and convergence orders of the piecewise quadratic CG solution for (37), with a = −1, b = 0.5, q = 0.1, 0.5, 0.9.  The above tables reveal that In Tables 4, 5 and 6 we present the errors and convergence orders of the piecewise cubic CG solution for (37), with q = 0.1, 0.5, 0.9. The numerical results reveal that the induced errors behave like For the nodal superconvergence, we also make a comparison of the numerical results between the quasi-geometric mesh and the uniform mesh. As it is not very easy to choose the same step-size for these two kinds of meshes, we select the relatively similar step-size. Numerical results are listed in the following Tables 7-8.
on the quasi-geometric mesh and on the uniform mesh.
In Tables 9 and 10 we show the error behavior of the piecewise quadratic and cubic CG solutions for (38), with a = −2, b = 1. The results in Tables 9 and 10 with θ(t) = t 2 . The exact solution is u(t) = cos(t). In initial subinterval J 0 = [0, t 0 ], the approximation of the exact solution u(t) is provided by Taylor expansion are φ(t) = 1 − t 2 2! + t 4 4! with m = 2 and φ(t) = 1 − t 2 2! + t 4 4! − t 6 6! with m = 3 . We choose quasi-graded meshes J h with last l + 1 nodes (l = 2, 4, 8, 16, 32, 64) being chosen equally spaced,Ñ = N + 1 = kl + 1 (k = κ + 1). Figure 1 exhibits the behavior of the error of the piecewise quadratic CG solution for (39). We find that superconvergence occurs at the Lobatto points and the mesh points. In Tables 12 and 13, we illustrate the error behavior of the piecewise    This illustrates the correctness of the theory.
6. Concluding remarks. The following three problems remain to be addressed in future research work.
• Analysis of the attainable order of local superconvergence of the DG method for pantograph-type DDEs under geometric meshes. • Postprocessing of the CG solutions for pantograph-type DDEs under geometric meshes. • Postprocessing of the DG solutions for pantograph-type DDEs under geometric meshes.