Convergence Analysis of the LDG Method for Singularly Perturbed Reaction-Diffusion Problems

We analyse the local discontinuous Galerkin (LDG) method for two-dimensional singularly perturbed reaction–diffusion problems. A class of layer-adapted meshes, including Shishkin- and Bakhvalov-type meshes, is discussed within a general framework. Local projections and their approximation properties on anisotropic meshes are used to derive error estimates for energy and “balanced” norms. Here, the energy norm is naturally derived from the bilinear form of LDG formulation and the “balanced” norm is artificially introduced to capture the boundary layer contribution. We establish a uniform convergence of order k for the LDG method using the balanced norm with the local weighted L2 projection as well as an optimal convergence of order k+1 for the energy norm using the local Gauss–Radau projections. The numerical method, the layer structure as well as the used adaptive meshes are all discussed in a symmetry way. Numerical experiments are presented.


Introduction
Over the past few decades, singularly perturbed problems have attracted considerable attention in the scientific community. Such problems arise in many applications, including the modelling of viscous fluid flows, semiconductor devices, and more [1]. For reaction-diffusion problems, difficulties arise owing to the presence of boundary layers in the solution. Unless the meshes are sufficiently refined, traditional finite-difference or finite-element methods on uniform or quasi-uniform meshes yield oscillatory and inaccurate numerical solutions. Consequently, three common approaches have been proposed in the literature. The first is to use traditional numerical methods on strongly refined, layer-adapted meshes, such as the Shishkin-type (S-type) or Bakhvalov-type (B-type) mesh [2][3][4][5]. Various parameter-uniform convergence results have been established in this way; notably, the order of convergence and error constant are independent of the singular perturbation parameters. The second approach is to use a stabilised numerical method, such as the streamline diffusion finite-element method, interior-penalty discontinuous Galerkin method, or local discontinuous Galerkin (LDG) method [6][7][8]; well-behaved local error estimates have been investigated using uniform or quasi-uniform meshes. The third approach [2,[9][10][11] is to combine the aforementioned stabilised numerical methods with layer-adapted meshes. From a practical perspective, the third approach is preferable because it is more stable and less sensitive to the choice of transition point on the layer-adapted mesh.
The LDG method is a form of finite-element method; it was first proposed as a generalisation of the discontinuous Galerkin (DG) method for a convection-diffusion • For singularly perturbed reaction-diffusion problems, the boundary layer structure is considerably more complicated because of the parabolic layers along all boundary edges [21]. As a result, the regularity of the solution is complex. This adds many difficulties to the theoretical analysis, such as in the construction of layer-adapted meshes and the estimates of various approximation errors. • Better than the convergence order obtained in [9] for the convection-diffusion problem, we can establish an optimal convergence of order k + 1 for the LDG method in the energy norm through a more elaborate analysis for the two-dimensional Gauss-Radau projections on anisotropic meshes. • In the reaction-diffusion region, the balanced-norm is more suitable to reflect the contribution of the boundary layer component [22]. To date, balanced-norm error estimates are only available for the Galerkin finite-element method (FEM) [23,24], mixed FEM [22], and hp-FEM [25], but not for the LDG method. For the first time, we establish the uniform convergence of the LDG method for the balanced norm.
Note that the numerical method, layer structure as well as the used adaptive meshes are all discussed in a symmetry way, which simplify our analysis to some extent.
The remainder of this paper is organised as follows. In Section 2, we describe the LDG method. In Section 3, we introduce a class of layer-adapted meshes and state some elementary lemmas for them. In Section 4, we establish uniform convergence for the balanced and energy norms. In Section 5, we present numerical experiments.
Let Ω N = {K ij } j=1,2,...,N y i=1,2,...,N x be a rectangle partition of Ω with element K ij = I i × J j , where I i = (x i−1 , x i ) and J j = (y j−1 , y j ). We set h x,i = x i − x i−1 , h y,j = y j − y j−1 , and h = min K ij ∈Ω N min{h x,i , h y,j }. We let V N = {v ∈ L 2 (Ω) : v| K ∈ Q k (K), K ∈ Ω N }, (2) Symmetry 2021, 13, 2291 3 of 26 be the discontinuous finite-element space, where Q k (K) represents the space of polynomials on K with a maximum degree of k in each variable. V N is contained in a broken Sobolev space, expressed as whose function is allowed to have discontinuities across element interfaces. For v ∈ H 1 (Ω N ) and y ∈ J j , j = 1, 2, . . . , express the traces evaluated from the four directions. We denote as the jumps on the vertical and horizontal edges, respectively. Rewrite (1) into an equivalent first-order system: Next, we replace the smooth functions v, s, r by test functions v, s, r, respectively, in the finite element space V N and the exact solution w = (u, p, q) by the approximate solution w = (u, p, q). As this function is discontinuous in each of its components, we must also replace the fluxes u i,y , u x,j , p i,y , q x,j by the numerical fluxes u i,y , u x,j , p i,y , q x,j .
Then, the LDG scheme is defined as follows. Find w = (u, p, q) ∈ V 3 N such that in each element K ij , the variational forms for y ∈ J j and j = 1, 2, . . . , N y . Here, λ i,y ≥ 0 (i = 0, 1, · · · , N x ) are stabilisation parameters to be determined later. Analogously, for x ∈ I i and i = 1, 2, . . . , N x , we can define q x,j and u x,j for j = 0, 1, . . . , N y .
Write w, v = ∑ K ij ∈Ω N w, v K ij . Then, we write the above LDG method into a compact form: where with

Layer-Adapted Meshes
To introduce the layer-adapted meshes, we extract some precise information from the exact solution of (1) and its derivatives in [21] (Theorem 2.2) and in [26] (Theorem 4.2).

Proposition 1.
Assume that the solution u of (1) can be decomposed by where S is a smooth part, W k is a boundary layer part and Z k is a corner layer part. More precisely, for 0 ≤ i, j ≤ k + 2, there exists a constant C independent of ε such that and so on for the remaining terms. Here, The layer-adapted mesh is constructed as follows: For notational simplification, we assume that N x = N y = N. Let N ≥ 4 be a multiple of four. We introduce the mesh points and consider a tensor-product mesh with mesh points (x i , y j ). Because both meshes have the same structure, we only describe the mesh in the x-direction.
Note that for these meshes and under the assumption
Note that for these meshes and under the assumption   In the following, we state two preliminary lemmas that will be frequently employed in the subsequent analysis. Because h x,i = h y,i , i = 1, 2, . . . , N, we simply use h i to denote one of them.

Lemma 1. Suppose that
Then, there exists a constant C > 0 independent of ε and N such that Because of the symmetry of the mesh and layer function, (14) is also valid for the cases with y i , 1 − y i−1 and 1 − x i−1 .

Convergence Analysis
In this section, we perform uniform convergence analysis for the LDG method on layer-adapted meshes. Two related norms are considered. The first is the energy norm, which is naturally derived from the formulation of the LDG method; that is, w 2 E ≡ B(w; w). Therefore, we obtain by using integration by parts and some trivial manipulations. Here, z 2 = ∑ K∈Ω N z 2 K and z 2 K = z, z K . However, this norm is inadequate for reaction-diffusion problems because the layer contributions are not "seen" In fact, letting , which vanishes as ε → 0. Thus, the following "balanced" norm is introduced: This norm is balanced in the sense that both the regular and layer parts are bounded away from zero uniformly in ε. We drop the penalty parameter λ because they will be taken as O( √ ε) in the subsequent analysis and computation of balanced-norm error. In the following subsections, we perform convergence analysis for these two norms. Different projections are introduced, and the related approximation properties are investigated.

Convergence of Balanced Norm
First, we analyse the LDG method for the balanced norm (20). Let ω ∈ C 1 (Ω N ) and ω ≥ ω 0 > 0 be a general weight function. We define the piecewise local weight L 2 projection Π ω as follows: For each element K ∈ Ω N and for any z ∈ L 2 (Ω N ), Π ω z ∈ V N satisfies In the special case of ω = 1, this operator reduces to the classical local L 2 projection, which is denoted by Π.
where C > 0 is independent of ε and N. A similar procedure applies for u and q in other spatial directions.
Proof. The proof is provided in the Appendix A.1.
N be the numerical solution of the LDG scheme (5) on layer-adapted meshes (13) when σ ≥ k + 1.5. Here, k ≥ 0 is the degree of the piecewise polynomials used in the discontinuous finite element space. Then, there exists a constant C > 0, independent of ε and N, such that From Proposition 1 and the consistency of numerical flux, we obtain the error equation: and we use (23c) to obtain Consequently, we have Analogously, we have Here, (23f) and (23g) are used. Finally, we use the Cauchy-Schwarz inequality, (23d), and the assumption that λ i,y = λ x,j = √ ε for i, j = 0, 1, · · · , N, to obtain From (27)-(29), we have which implies from the definitions (19) and (20) and the choice of λ i,y = λ x,j = √ ε (i, j = 0, 1, · · · , N). Using (23) and a trivial inequality, we derive (24).
Remark that for the meshes listed in Table 1, one has balanced-norm error estimate for BS-,B-type mesh.

Improvement of Convergence in Energy Norm
By the analysis in Theorem 1, one can only obtain the energy-error bound of order O( 4 √ εN −k (max |ψ |) k+1/2 ). In this subsection, we perform an elaborate analysis and establish a more sharp convergence result in the energy norm. The following local Gauss-Radau projections are required.
For each element K ij ∈ Ω N and for any Lemma 5. [9] (Inequality (4.4)) There exists a constant C > 0, independent of the element size and z, such that

Lemma 6.
Let σ ≥ k + 1.5. Then, there exists a constant C > 0 independent of ε and N such that where is given by (17). Similarly, we can obtain the same conclusions in another spatial direction.
Proof. The proof is provided in the Appendix A.2.

Numerical Experiments
In this section, we present some numerical experiments. All calculations are conducted in MATLAB R2015B. The system of linear equations resulting from the discrete problems are solved by the lower-upper (LU)-decomposition algorithm. All integrals are evaluated using the 5-point Gauss-Legendre quadrature rule.
The LDG method (5) is applied to the layer-adapted meshes presented in Table 1, where σ = k + 1, k = 0, 1, 2, 3. We let e N be the error in either e E or e B for an N-element. As assumed in Theorems 1-2, we take the flux parameter λ 0 = λ N = √ ε and λ i = 0 (i = 1, 2, . . . , N − 1) to compute the energy-error e E , and choose λ i = √ ε (i = 0, 1, . . . , N) to compute the balanced-error e B . The corresponding convergence rates are computed by the following formulae: with p = 2 ln N/ ln(2N). Here, we use r S to reflect the convergence rate with respect to the power of N −1 ln N, see the error estimates in (32) and (36).

Example 1. Consider a linear reaction-diffusion problem
where f (x, y) is suitably taken such that the exact solution is u(x, y) = g(x)g(y) with We set ε = 10 −8 , small enough to bring out the singularly perturbed nature of (44). In Table 2, we list the balanced-norm errors and their convergence rates. The convergence rate is plotted in Figure 2. We observe convergence of order k + 1/2, which is a half-order superior to the estimate from (24). In Table 3, we present the energy norm errors and their convergence rates. The convergence rate is plotted in Figure 3, which agrees with our estimate (36).
We show the relevance of these errors to the small parameter ε. We let N = 256, k = 1, and vary the values of ε. We see that the errors in the balanced norm are almost unchanged, whereas the errors in the energy norm change slightly. For a visual understanding, we plot the energy errors via ε on log-log coordinates. In Figure 4, we observe the subtle influence of the ε 0.25 -factor on the energy errors; the results agree with our predictions.

Example 2. Consider the nonlinear reaction-diffusion problem
where f (x, y) is suitably taken such that the exact solution is u(x, y) = h(x)h(y) and We let ε = 10 −8 . In Tables 4 and 5, we list the error and convergence rates for the balanced and energy norms, respectively. We still observe convergence of orders k + 1/2 and k + 1 for the balanced-norm and energy norm errors. See the convergence rates for the balanced-norm and energy-norm errors in Figures 5 and 6.
Moreover, we test the dependence of these two types of errors on ε. We observe that the errors in the balanced norm are almost constant, whereas the errors in the energy norm change slightly. In Figure 7, we confirm the influence of the factor ε 0.25 on the energy norm errors obtained for the three layer-adapted meshes. Note that for this example, the ε 0.25 -factor is clearly observed on both the S-type and B-type meshes. This may be due to the fact that the regular part of the exact solution belongs to V N , as described in [10] (Theorem 3.5).

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A
In this part, we provide the technical proofs for Lemmas 4 and 6.
(2) Prove (23b) and (23c). It can be seen that For K ij ∈ Ω XL , we obtain from the L ∞ -stability (22a), the L ∞ -approximation (22b), and Lemma 1 that where we have used (10) and h i / √ ε ≥ CN −1 ≥ h j . For K ij ∈ Ω XM ∪ Ω XR , we obtain from the L ∞ -stability (22a) and σ ≥ k + 1 that Consequently, we obtain from (14) that using σ ≥ k + 1. In a similar fashion, we obtain From the solution decomposition and similar arguments for the other terms, we arrive at (23b) and (23c).
(3) Prove (23d). We start from the following inequality: For the first term, we notice that and proceed as before. For the second term, we use (23c). Thus, (23d) follows. The remaining inequalities of (23) can be proved analogously; we omit the details here.

Appendix A.2 Proof of Lemma 6
The conclusions are more precise than Lemma 4.1 of [9]. The proof proceeds similarly to Lemma 4. We mention several differences and use the same notations to prevent confusion.
Using the L 2 -stability, the uniformity of the mesh in Ω XM , and Proposition 1, we have To bound B 3 , we decompose it into two parts:
(2) Note that (η u ) − N,y = u N,y − π − y (u N,y ), where π − y is a one-dimensional Gauss-Radau projection regarding y and satisfies analogous stability and approximation conditions to that in Lemma 5. From the solution decomposition, we express u N,y = S N,y + E N,y + F N,y , where S N,y , E N,y and F N,y are functions of one variable y and satisfy |S N,y | ≤ Cε −j/2 e −β(1−y)/ √ ε . Following the similar line to that used to prove (34a), we obtain Using the triangle inequality leads to (35b). The proofs of (35c) and (35d) are similar and therefore omitted.