On the Local Discontinuous Galerkin Method for Linear Elasticity

Abstract

Following Castillo et al. (2000) and Cockburn (2003), a general framework of constructing discontinuous Galerkin (DG) methods is developed for solving the linear elasticity problem. The numerical traces are determined in view of a discrete stability identity, leading to a class of stable DG methods. A particular method, called the LDG method for linear elasticity, is studied in depth, which can be viewed as an extension of the LDG method discussed by Castillo et al. (2000) and Cockburn (2003). The error bounds in 𝐿²-norm, 𝐻¹-norm, and a certain broken energy norm are obtained. Some numerical results are provided to confirm the convergence theory established.

1. Introduction

This paper is focused on systematically studying discontinuous Galerkin (DG) methods for the linear elasticity problem. Since the DG method was first introduced in 1970s, these methods have been applied for solving a variety of mathematical-physical problems including linear and nonlinear hyperbolic problems, Navier-Stokes equations, convection-dominated diffusion problems, and so on. The DG method may be viewed as high-order extensions of the classical finite volume method. Since no inter-element continuity is imposed, such methods can be defined on very general meshes including nonconforming meshes. Moreover, polynomials of arbitrary degree can be used on each element, making these methods suitable for ℎ𝑝-adaptivity. We refer to [1] for an excellent historical survey. As a generalization of the DG method proposed in [2] for the solution of the compressible Navier-Stokes equations, some local discontinuous Galerkin (LDG) methods were introduced and analyzed in [3, 4] for time-dependent convection-diffusion systems and second-order elliptic problems, respectively. After that, a unified theory was developed for conducting error analysis for DG methods of elliptic problems (cf. [5]). And a new framework was proposed in [6] for designing and analyzing DG methods, where stabilization mechanism frequently used in DG methods was also investigated. In [7], an important framework was proposed for constructing stable DG methods.

On the other hand, linear elasticity discusses how solid objects deform and become internally stressed on prescribed loading conditions, and is encountered extensively in structural analysis and engineering design. By use of usual displacement-based finite-element methods for the elasticity equation, we are able to numerically determine the displacement field directly. However, in many engineering applications, the stress field is a quantity of more interest. For this, we may apply mixed finite-element methods to solve the linear elasticity system, from which the aforementioned two physical quantities can be obtained simultaneously. Over the past four decades, there have been many efforts along this line. But due to the symmetry constraint on the stress tensor, it is extremely difficult to construct stable stress-displacement finite elements. In two spatial dimensions, the first stable finite element with polynomial shape functions is presented in [8]. For the lowest-order element, the discrete stress space is composed of certain piecewise cubic polynomials, while the discrete displacement space consists of piecewise linear polynomials. In three dimensions, a piecewise quadratic stress space is constructed with 162 degrees of freedom on each tetrahedron (cf. [9]). Another approach in this direction is to use composite elements (macroelements), in which the discrete displacement space consists of piecewise polynomials with respect to one triangulation of the domain, while the discrete stress space consists of piecewise polynomials with respect to a more refined triangulation (cf. [10–13]). It is mentioned that for solving the previous problem, several mixed elements with weakly imposed symmetry have also been developed (cf. [14–16]).

Regarding the complexity of mixed elements given above, the discontinuous Galerkin method is naturally a suitable alternative for numerically solving linear elasticity problems. To the best of our knowledge, a local DG method and an interior penalty DG method for linear elasticity are presented in [17] and [18, 19], respectively. A mixed DG method is given in [20], which one may find is covered by our general formulation below. Moreover, a mixed formulation is also extended to the case of nonsymmetric stress tensors (cf. [21]).

In this paper, we are going to look for new DG methods for the linear elasticity problem. Following [4, 7], we build up a framework to construct our DG methods. Then a discrete stability identity is established, from which we derive feasible choices of numerical traces and get a class of stable DG methods for linear elasticity. With a parameter taken to be zero, namely, 𝐶₂₂=0, we get an LDG method. Following [5] and using some technical arguments, we get optimal-order error estimates for the LDG method proposed in a certain broken energy norm, 𝐻¹-norm, and 𝐿²-norm, respectively. It should be emphasized that in our formulation the symmetry of the stress tensor is preserved automatically.

The rest of this paper is organized as follows. The basic framework of DG methods and the determination of numerical traces are presented in Section 2. An error analysis for the LDG method is given in Section 3. And in Section 4, some numerical results are included to confirm our theoretical convergence orders.

2. The DG Method for Linear Elasticity

2.1. Basic Framework for the DG Method

Assume that Ω⊂ℝ^𝑑(𝑑=2,3) is a bounded polygon or polyhedron. Let 𝝈=(𝜎_𝑖𝑗)_𝑑×𝑑 be the stress, let 𝐮=(𝑢₁,…,𝑢_𝑑)^𝑡 be the displacement, and let 𝐟=(𝑓₁,…,𝑓_𝑑)^𝑡 be the body force. Denote by 𝜺(𝐮)∶=(𝜀_𝑖𝑗(𝐮))_𝑑×𝑑 the linearized strain tensor with 𝜀_𝑖𝑗(𝐮)∶=(𝜕𝑢_𝑖/𝜕𝑥_𝑗+𝜕𝑢_𝑗/𝜕𝑥_𝑖)/2, tr the trace operator, and ∇⋅ the divergence operator. Then, the mathematical model of linear elasticity reads (cf. [22])𝒜𝝈−𝜺(𝐮)=𝟎inΩ,−∇⋅𝝈=𝐟inΩ,𝐮=𝟎on𝜕Ω, (2.1) where 𝒜 is the fourth-order compliance tensor defined by1𝒜𝝈=𝜆2𝜇(𝝈−𝑑𝜆+2𝜇(tr𝝈)𝜹), (2.2) where 𝜹∶=(𝛿_𝑖𝑗)_𝑑×𝑑 is the Kronecker tensor, and the positive constants 𝜆 and 𝜇 are the Lamé coefficients of the elastic material under consideration.

For defining our DG methods for problem (2.1), we introduce some notations first of all. For any Banach space 𝐵, the subspace of symmetric matrix-valued function is denoted by (𝐵)^𝑠_𝑑×𝑑. Given a bounded domain 𝐺⊂ℝ^𝑑 and a nonnegative integer 𝑚, let 𝐻^𝑚(𝐺) be the usual Sobolev space consisting of all functions in 𝐿²(𝐺) whose weak derivatives up to degree 𝑚 are also 𝐿²(𝐺)-integrable (cf. [23]). The corresponding norm is denoted by ‖⋅‖_𝑚,𝐺. Let 𝐻^𝑚₀(𝐺) be the closure of 𝐶^∞₀(𝐺) with respect to the norm ‖⋅‖_𝑚,𝐺. Moreover, we simply write ‖⋅‖_𝑚 for ‖⋅‖_𝑚,Ω.

Let {𝒯_ℎ}_ℎ>0 be a regular family of triangulations of Ω (cf. [23, 24]); ℎ∶=max_{𝐾∈𝒯ℎ}ℎ_𝐾 and ℎ_𝐾∶=diam(𝐾). Let ℰ_ℎ be the union of all faces (edges) of the triangulation 𝒯_ℎ and ℰ^𝑖_ℎ the union of all interior faces (edges) of the triangulation 𝒯_ℎ. For any 𝑒∈ℰ_ℎ, ℎ_𝑒 is set to be the diameter of 𝑒. Based on the triangulation 𝒯_ℎ, let 𝚺∶={𝝉∈(𝐿²(Ω))^𝑠_𝑑×𝑑;𝜏_𝑖𝑗∣_𝐾∈𝐻¹(𝐾),∀𝐾∈𝒯_ℎ,𝑖,𝑗=1,…,𝑑},𝐕∶={𝐯∈(𝐿²(Ω))^𝑑;𝑣_𝑖∣_𝐾∈𝐻¹(𝐾),∀𝐾∈𝒯_ℎ,𝑖=1,…,𝑑}. (2.3) The corresponding finite element spaces are given by 𝚺_ℎ∶={𝝉∈(𝐿²(Ω))^𝑠_𝑑×𝑑;𝜏_𝑖𝑗∣_𝐾∈𝒮₁(𝐾),∀𝐾∈𝒯_ℎ𝐕,𝑖,𝑗=1,…,𝑑},_ℎ∶={𝐯∈(𝐿²(Ω))^𝑑;𝑣_𝑖∣_𝐾∈𝒮₂(𝐾),∀𝐾∈𝒯_ℎ,𝑖=1,…,𝑑}, (2.4) where, for each 𝐾∈𝒯_ℎ, 𝒮₁(𝐾) and 𝒮₂(𝐾) are two finite-dimensional spaces of polynomials in 𝐾 containing 𝑃_𝑙(𝐾) and 𝑃_𝑘(𝐾), respectively, with 𝑘,𝑙≥0. Here, for a nonnegative integer 𝑚, 𝑃_𝑚(𝐾) stands for the set of all polynomials in 𝐾 with the total degree no more than 𝑚.

To guarantee uniqueness of the solution to the LDG method to be proposed, we always assume that𝜺(𝐕_ℎ)⊂𝚺_ℎ,𝒜𝚺_ℎ⊂𝚺_ℎ (2.5) For a function 𝑣∈𝐿²(Ω) with 𝑣∣_𝐾∈𝐻^𝑚(𝐾) for all 𝐾∈𝒯_ℎ, let ‖𝑣‖_𝑚,ℎ be the usual broken 𝐻^𝑚-type norm of 𝑣 defined by‖𝑣‖_𝑚,ℎ=(_{𝐾∈𝒯ℎ}‖𝑣‖²_𝑚,𝐾)^1/2. (2.6) If 𝐯 is a vector-valued or matrix-valued function, the corresponding term ‖𝐯‖_𝑚,ℎ is defined in a similar manner. For a vector or a matrix 𝐯, denote by |𝐯| the quantity (𝐯⋅𝐯)^1/2 or (𝐯∶𝐯)^1/2. Here the symbol ∶ stands for the double dot product operation of matrices. Throughout this paper, we use the notation “≲⋯” to indicate “≤𝐶⋯”, where 𝐶 is a generic positive constant independent of ℎ and other parameters, which may take different values at different appearances.

Let 𝐾⁺ and 𝐾⁻ be two adjacent elements of 𝒯_ℎ. Let 𝐱 be an arbitrary point of the set 𝑒^′=𝜕𝐾⁺∩𝜕𝐾⁻, and let 𝐧⁺ and 𝐧⁻ be the corresponding outward unit normals at that point. For a vector-valued function 𝐯 smooth inside each element 𝐾^±, let us denote by 𝐯^± the trace of 𝐯 on 𝑒^′ from the interior of 𝐾^±. Then we define averages and jumps at 𝐱∈𝑒^′ as follows: 1{𝑝}=2(𝑝⁺+𝑝⁻),[𝑝]=𝑝⁺𝐧⁺+𝑝⁻𝐧⁻,1{𝐯}=2(𝐯⁺+𝐯⁻),[𝐯]=𝐯⁺⋅𝐧⁺+𝐯⁻⋅𝐧⁻,1{𝝉}=2(𝝉⁺+𝝉⁻),[𝝉]=𝝉⁺𝐧⁺+𝝉⁻𝐧⁻. (2.7) If 𝐱 is on an edge/face 𝑒 lying on the boundary 𝜕Ω, the above terms are defined by {𝑝}=𝑝,[𝑝]=𝑝𝐧,{𝐯}=𝐯,[𝐯]=𝐯⋅𝐧,{𝝉}=𝝉,[𝝉]=𝝉𝐧, (2.8) where 𝐧 is the unit outward normal vector on 𝜕Ω. We define a matrix-valued jump ⟦⋅⟧ of a vector 𝐯 as follows: 1⟦𝐯⟧=2(𝐯⁺⊗𝐧⁺+𝐧⁺⊗𝐯⁺+𝐯⁻⊗𝐧⁻+𝐧⁻⊗𝐯⁻),if𝐱∈𝑒∈ℰ^𝑖_ℎ,1⟦𝐯⟧=2(𝐯⊗𝐧+𝐧⊗𝐯),if𝐱∈𝑒∈ℰ_ℎ∩𝜕Ω, (2.9) where 𝐯⊗𝐧 denotes the matrix whose (𝑖,𝑗)th component is 𝑣_𝑖𝑛_𝑗 for two vectors 𝐯 and 𝐧.

Now, we are ready to introduce a framework to derive DG methods for problem (2.1). Following [4, 7], we first derive a variational formulation for problem (2.1). Taking a double dot product with a matrix-valued function 𝝉 on both sides of the first equation of (2.1) and integrating by parts over 𝐾, we have_𝐾𝒜𝝈∶𝝉𝑑𝑥=−_𝐾𝐮⋅(∇⋅𝝉)𝑑𝑥+_𝜕𝐾𝐮⋅(𝝉𝐧)𝑑𝑠. (2.10) Multiplying the second equation of (2.1) by a vector-valued function 𝐯 and integrating by parts over 𝐾 yields_𝐾𝝈∶𝜺(𝐯)𝑑𝑥=_𝜕𝐾𝐯⋅(𝝈𝐧)𝑑𝑠+_𝐾𝐟⋅𝐯𝑑𝑥. (2.11)

Motivated by the above two identities, we may define our DG method as follows. Find an approximate solution (𝝈_ℎ,𝐮_ℎ)∈𝚺_ℎ×𝐕_ℎ such that_𝐾𝒜𝝈_ℎ∶𝝉𝑑𝑥=−_𝐾𝐮_ℎ⋅(∇⋅𝝉)𝑑𝑥+_𝜕𝐾𝐮_ℎ⋅(𝝉𝐧)𝑑𝑠,_𝐾𝝈_ℎ∶𝜺(𝐯)𝑑𝑥=_𝜕𝐾𝐯⋅(𝝈_ℎ𝐧)𝑑𝑠+_𝐾𝐟⋅𝐯𝑑𝑥, (2.12)(2.13) for all (𝝉,𝐯)∈𝚺_ℎ×𝐕_ℎ and all 𝐾∈𝒯_ℎ. Note that any function with the hat superscript is only defined over all faces of the triangulation 𝒯_ℎ, which is called a numerical trace. Since 𝝈_ℎ is symmetric, it is natural to choose 𝝈_ℎ as a symmetric matrix-valued function. Moreover, we only consider the case where the numerical traces are single valued over all faces (conservation).

2.2. Numerical Traces and the LDG Method

We begin by producing a stability identity for the continuous problem (2.1), a crucial step in constructing feasible numerical traces to get a stable DG method from (2.12)-(2.13). For this, taking a double dot product with 𝝈 on both sides of the first equation of (2.1) and then integrating over Ω, we have_Ω𝒜𝝈∶𝝈𝑑𝑥−_Ω𝜺(𝐮)∶𝝈𝑑𝑥=0. (2.14) Multiplying the second equation of (2.1) by 𝐮 and then integrating over Ω again, we find from the homogeneous boundary condition of 𝐮 that_Ω𝜺(𝐮)∶𝝈𝑑𝑥=_Ω𝐟⋅𝐮𝑑𝑥. (2.15) Now adding these two equations gives_Ω𝒜𝝈∶𝝈𝑑𝑥=_Ω𝐟⋅𝐮𝑑𝑥. (2.16) This is the stability identity corresponding to the continuous problem (2.1).

Next, we mimic the above derivation to get a discrete analogue of the stability identity (2.16) for the DG method (2.12)-(2.13). Taking 𝝉=𝝈_ℎ in (2.12) and summing over all 𝐾∈𝒯_ℎ, we have_Ω𝒜𝝈_ℎ∶𝝈_ℎ𝑑𝑥−_{𝐾∈𝒯ℎ}(−_𝐾𝐮_ℎ⋅(∇⋅𝝈_ℎ)𝑑𝑥+_𝜕𝐾𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠)=0. (2.17) Similarly, summing up (2.13) over all 𝐾∈𝒯_ℎ with 𝐯=𝐮_ℎ, we come to_{𝐾∈𝒯ℎ}(_𝐾𝝈_ℎ∶𝜺(𝐮_ℎ)𝑑𝑥−_𝜕𝐾𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠)=_Ω𝐟⋅𝐮_ℎ𝑑𝑥. (2.18) Adding the last two equations gives_Ω𝒜𝝈_ℎ∶𝝈_ℎ𝑑𝑥+Θ_ℎ=_Ω𝐟⋅𝐮_ℎ𝑑𝑥, (2.19) where Θ_ℎ∶=_{𝐾∈𝒯ℎ𝐾}(𝐮_ℎ⋅(∇⋅𝝈_ℎ)+𝝈_ℎ∶𝜺(𝐮_ℎ))𝑑𝑥−_{𝐾∈𝒯ℎ𝜕𝐾}(𝐮_ℎ⋅(𝝈_ℎ𝐧)−𝐮_ℎ⋅(𝝈_ℎ𝐧))𝑑𝑠. (2.20)

Lemma 2.1. Assume 𝝈_ℎ and 𝝈_ℎ are both symmetric, and the numerical traces 𝝈_ℎ and 𝐮_ℎ are single valued over all 𝑒∈ℰ_ℎ. Then _{𝐾∈𝒯ℎ𝜕𝐾⧵𝜕Ω}𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠=_{𝑒∈ℰ𝑖ℎ𝑒}[𝝈_ℎ]⋅{𝐮_ℎ}+{𝝈_ℎ}∶⟦𝐮_ℎ⟧𝑑𝑠,_{𝐾∈𝒯ℎ𝜕𝐾⧵𝜕Ω}𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠=_{𝑒∈ℰ𝑖ℎ𝑒}[𝝈_ℎ]⋅𝐮_ℎ𝑑𝑠,_{𝐾∈𝒯ℎ𝜕𝐾⧵𝜕Ω}𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠=_{𝑒∈ℰ𝑖ℎ𝑒}𝝈_ℎ∶⟦𝐮_ℎ⟧𝑑𝑠. (2.21)(2.22)(2.23)

Proof. We only prove the first identity, and the other two ones can be obtained in the similar manners. The left-hand side of (2.21) can be rewritten as _{𝐾∈𝒯ℎ𝜕𝐾⧵𝜕Ω}𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠=_{𝑒∈ℰ𝑖ℎ𝑒}(𝐮⁺_ℎ⋅(𝝈⁺_ℎ𝐧⁺)+𝐮⁻_ℎ⋅(𝝈⁻_ℎ𝐧⁻))𝑑𝑠. (2.24) On the other hand, _{𝑒∈ℰ𝑖ℎ𝑒}[𝝈_ℎ]⋅{𝐮_ℎ}𝑑𝑠=_{𝑒∈ℰ𝑖ℎ𝑒}12(𝐮⁺_ℎ⋅(𝝈⁺_ℎ𝐧⁺)+𝐮⁺_ℎ⋅(𝝈⁻_ℎ𝐧⁻)+𝐮⁻_ℎ⋅(𝝈⁺_ℎ𝐧⁺)+𝐮⁻_ℎ⋅(𝝈⁻_ℎ𝐧⁻))𝑑𝑠,_{𝑒∈ℰ𝑖ℎ𝑒}{𝝈_ℎ}∶⟦𝐮_ℎ⟧𝑑𝑠=_{𝑒∈ℰ𝑖ℎ𝑒}12(𝐧⁺⋅(𝝈⁺_ℎ𝐮⁺_ℎ)+𝐧⁻⋅(𝝈⁺_ℎ𝐮⁻_ℎ)+𝐧⁺⋅(𝝈⁻_ℎ𝐮⁺_ℎ)+𝐧⁻⋅(𝝈⁻_ℎ𝐮⁻_ℎ))𝑑𝑠. (2.25) Hence, the identity (2.21) is a direct consequence of (2.24)-(2.25).

By Lemma 2.1 and using integration by parts twice, we see that Θ_ℎ=_{𝐾∈𝒯ℎ}(_𝐾∇⋅(𝝈_ℎ𝐮_ℎ)𝑑𝑥−_𝜕𝐾𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠−_𝜕𝐾𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠)=_{𝐾∈𝒯ℎ}(_𝜕𝐾𝐧⋅(𝝈_ℎ𝐮_ℎ)𝑑𝑠−_𝜕𝐾𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠−_𝜕𝐾𝐮_ℎ⋅(𝝈_ℎ𝐧)𝑑𝑠)=_{𝑒∈ℰ𝑖ℎ𝑒}([𝝈_ℎ]⋅({𝐮_ℎ}−𝐮_ℎ)+({𝝈_ℎ}−𝝈_ℎ)∶⟦𝐮_ℎ+⟧)𝑑𝑠_𝜕Ω(𝐮_ℎ⋅((𝝈_ℎ−𝝈_ℎ)𝐧)−𝐮_ℎ⋅(𝝈_ℎ𝐧))𝑑𝑠. (2.26) Thus, if 𝑒∈ℰ^𝑖_ℎ, we take𝝈_ℎ={𝝈_ℎ}−𝐶₁₁⟦𝐮_ℎ𝐮⟧,_ℎ={𝐮_ℎ}−𝐶₂₂[𝝈_ℎ], (2.27) and if 𝑒∈ℰ_ℎ∩𝜕Ω, we take𝝈_ℎ=𝝈_ℎ−𝐶₁₁⟦𝐮_ℎ𝐮⟧,_ℎ=0, (2.28) where 𝐶₁₁ and 𝐶₂₂ are two nonnegative continuous functions on 𝑒. When 𝐶₂₂=0, the corresponding method is called the LDG method for linear elasticity, viewed as a generalization of the LDG method for second-order elliptic problems in [4, 7].

For the above choice of numerical traces, we have by some direct manipulation thatΘ_ℎ=_{𝑒∈ℰ𝑖ℎ𝑒}(𝐶₂₂[𝝈_ℎ]²+𝐶₁₁⟦𝐮_ℎ⟧²)𝑑𝑠+_𝜕Ω𝐶₁₁⟦𝐮_ℎ⟧²𝑑𝑠≥0. (2.29) Then, the combination of (2.19) and (2.29) allows a discrete stability estimate:‖𝝈_ℎ‖²₀≲_Ω𝐟⋅𝐮_ℎ𝑑𝑥, (2.30) which is essential in constructing a reliable DG method (cf. [7]).

Let us further show the unique solvability of problem (2.12)-(2.13) with the numerical traces given by (2.27) and (2.28), whenever 𝐶₁₁>0 and the finite element spaces 𝚺_ℎ and 𝐕_ℎ satisfy condition (2.5). In fact, it suffices to verify that this DG method only has zero solution when 𝐟=𝟎. By setting 𝐟=𝟎, the stability identity (2.19) gives_Ω𝒜𝝈_ℎ∶𝝈_ℎ𝑑𝑥+_{𝑒∈ℰ𝑖ℎ𝑒}(𝐶₂₂[𝝈_ℎ]²+𝐶₁₁⟦𝐮_ℎ⟧²)𝑑𝑠𝑑𝑠+_𝜕Ω𝐶₁₁⟦𝐮_ℎ⟧²𝑑𝑠=0, (2.31) which implies that 𝝈_ℎ=𝟎 in Ω, ⟦𝐮_ℎ⟧=𝟎 on ℰ^𝑖_ℎ, and 𝐮_ℎ=𝟎 on 𝜕Ω, provided that 𝐶₁₁>0. Therefore, from (2.12) the definition of 𝐮_ℎ (cf. (2.27) and (2.28)), and Lemma 2.1, we know that_𝐾𝝉∶𝜺(𝐮_ℎ)𝑑𝑥=0,∀𝝉∈𝚺_ℎ, (2.32) which, due to (2.5), implies that 𝜺(𝐮_ℎ)=𝟎. By Korn's inequality (2.34) on the discontinuous finite element space, given below, it is easy to see that ‖𝐮_ℎ‖_1,ℎ=0. Then with 𝐮_ℎ=𝟎 on 𝜕Ω and ⟦𝐮_ℎ⟧=𝟎, we conclude that 𝐮_ℎ=𝟎 in Ω, as required.

Remark 2.2. The standard Korn's inequality states that (cf. [25]) ‖𝐯‖²₁≲‖𝜺(𝐯)‖²₀+_𝜕Ω|𝐯|²𝑑𝑠,∀𝐯∈𝐻¹(Ω)^𝑑. (2.33) The following Korn's inequality on the discontinuous finite-element space is given in [26]: ‖𝐯‖²_1,ℎ≲‖𝜺(𝐯)‖²_0,ℎ+_{𝑒∈ℰℎ}ℎ_𝑒⁻¹‖⟦𝐯⟧‖²_0,𝑒,∀𝐯∈𝐕_ℎ. (2.34)

As in [4, 7], the DG method (2.12)-(2.13) with the numerical traces (2.27)-(2.28) can also be written in a mixed formulation (cf. [27]). After some direct manipulation, the approximate solution (𝝈_ℎ,𝐮_ℎ) can be characterized as the unique solution of the following variational problem. Find (𝝈_ℎ,𝐮_ℎ)∈𝚺_ℎ×𝐕_ℎ such that 𝑎(𝝈_ℎ,𝝉)+𝑏(𝐮_ℎ,𝝉)=0,−𝑏(𝐯,𝝈_ℎ)+𝑐(𝐮_ℎ,𝐯)=𝐹(𝐯), (2.35)(2.36) for all (𝝉,𝐯)∈𝚺_ℎ×𝐕_ℎ, where 𝑎(𝝈,𝝉)∶=_Ω𝒜𝝈∶𝝉𝑑𝑥+_ℰ𝑖ℎ𝐶₂₂[𝝈]⋅[𝝉]𝑑𝑠,𝑏(𝐯,𝝉)∶=_{𝐾∈𝒯ℎ𝐾}𝐯⋅(∇⋅𝝉)𝑑𝑥−_ℰ𝑖ℎ{𝐯}⋅[𝝉]𝑑𝑠,𝑐(𝐮,𝐯)∶=_ℰℎ𝐶₁₁⟦𝐮⟧∶⟦𝐯⟧𝑑𝑠,𝐹(𝐯)∶=_Ω𝐟⋅𝐯𝑑𝑥. (2.37)

Remark 2.3. If we take 𝐶₁₁=0 and 𝐶₂₂=𝛼ℎ^𝑠_𝑒, the mixed formulation (2.35)-(2.36) will be reduced to the one adopted in [20].

3. Error Analysis for the LDG Method

In this section, we provide an error analysis for the DG method (2.35)-(2.36) in the case 𝐶₂₂=0, which is named as the LDG method in the last section. For this purpose, we are going to derive a primal formulation for the approximate method. First we introduce a global lifting operator 𝐫∶(𝐿²(ℰ_ℎ))^𝑠_𝑑×𝑑→𝚺_ℎ defined by_Ω𝐫(𝝓)∶𝝉𝑑𝑥=−_ℰℎ𝝓∶{𝝉}𝑑𝑠,∀𝝉∈𝚺_ℎ,𝝓∈(𝐿²(ℰ_ℎ))^𝑠_𝑑×𝑑. (3.1) Moreover, for each 𝑒∈ℰ_ℎ, we introduce a local lifting operator 𝐫_𝑒∶(𝐿²(𝑒))^𝑠_𝑑×𝑑→𝚺_ℎ defined by_Ω𝐫_𝑒(𝜙)∶𝝉𝑑𝑥=−_𝑒𝝓∶{𝝉}𝑑𝑠,∀𝝉∈𝚺_ℎ,𝝓∈(𝐿²(𝑒))^𝑠_𝑑×𝑑. (3.2) Then it is easy to check that 𝐫_𝑒(𝝓) is only nonzero in the triangles with 𝑒 as one edge, and there holds the identity𝐫(𝝓)=_{𝑒∈ℰℎ}𝐫_𝑒(𝝓∣_𝑒),∀𝝓∈(𝐿²(ℰ_ℎ))^𝑠_𝑑×𝑑. (3.3) Now we express 𝝈_ℎ in terms of 𝐮_ℎ. From the first equation in Lemma 2.1 and (2.35), it follows that _Ω𝒜𝝈_ℎ∶𝝉𝑑𝑥=−_{𝐾∈𝒯ℎ𝐾}𝐮_ℎ⋅(∇⋅𝝉)𝑑𝑥+_ℰ𝑖ℎ{𝐮_ℎ}⋅[𝝉]𝑑𝑠=_{𝐾∈𝒯ℎ𝐾}𝜺(𝐮_ℎ)∶𝝉𝑑𝑥−_ℰℎ⟦𝐮_ℎ⟧∶{𝝉}𝑑𝑠=_{𝐾∈𝒯ℎ𝐾}𝜺(𝐮_ℎ)∶𝝉𝑑𝑥+_Ω𝐫(⟦𝐮_ℎ⟧)∶𝝉𝑑𝑥, (3.4) for all 𝝉∈𝚺_ℎ. Then by (2.5), we get𝒜𝝈_ℎ=𝜺(𝐮_ℎ)+𝐫(⟦𝐮_ℎ⟧), (3.5) that is, 𝝈_ℎ=2𝜇(𝜺(𝐮_ℎ)+𝐫(⟦𝐮_ℎ⟧))+𝜆tr(𝜺(𝐮_ℎ)+𝐫(⟦𝐮_ℎ⟧))𝜹. (3.6) Substituting 𝝈_ℎ from the last equation into (2.36), we get a primal formulation for the LDG method as follows. Find 𝐮_ℎ∈𝐕_ℎ such that𝑎_ℎ(𝐮_ℎ,𝐯)=_Ω𝐟⋅𝐯𝑑𝑥,∀𝐯∈𝐕_ℎ, (3.7) where 𝑎_ℎ(𝐰,𝐯)∶=_{𝐾∈𝒯ℎ𝐾}2𝜇(𝜺(𝐰)+𝐫(⟦𝐰⟧))∶(𝜺(𝐯)+𝐫(⟦𝐯⟧))𝑑𝑥+_{𝐾∈𝒯ℎ𝐾}𝜆tr(𝜺(𝐰)+𝐫(⟦𝐰⟧))tr(𝜺(𝐯)+𝐫(⟦𝐯⟧))𝑑𝑥+_ℰℎ𝐶₁₁⟦𝐰⟧∶⟦𝐯⟧𝑑𝑠. (3.8)

Next, we consider the consistency of the method (3.7). Assume that (𝝈,𝐮)∈(𝐻^𝑠+1(Ω))^𝑠_𝑑×𝑑×𝐻^𝑠+2(Ω)^𝑑, where 𝑠≥0 is some nonnegative integer, and 𝐮 is the exact solution of (2.1). It is easy to see that ⟦𝐮⟧=0 on ℰ_ℎ. Therefore, for any 𝐯∈𝐕_ℎ, we have 𝑎_ℎ(𝐮,𝐯)=_{𝐾∈𝒯ℎ𝐾}(2𝜇𝜺(𝐮)∶(𝜺(𝐯)+𝐫(⟦𝐯⟧))+𝜆tr(𝜺(𝐮))tr(𝜺(𝐯)+𝐫(⟦𝐯⟧)))𝑑𝑥=_{𝐾∈𝒯ℎ𝐾}𝝈∶(𝜺(𝐯)+𝐫(⟦𝐯⟧))𝑑𝑥. (3.9) From (3.9), the definition of lifting operator 𝐫 (see (3.1)), and the fact that {𝝈}=𝝈, we further find that𝑎_ℎ(𝐮,𝐯)=𝑎(𝐮,𝐯)+_ℰℎ⟦𝐯⟧∶{𝝈−𝐏_ℎ𝝈}𝑑𝑠+_{𝐾∈𝒯ℎ𝐾}(𝝈−𝐏_ℎ𝝈)∶𝐫(⟦𝐯⟧)𝑑𝑥, (3.10) where𝑎(𝐮,𝐯)∶=_{𝐾∈𝒯ℎ𝐾}𝝈∶𝜺(𝐯)𝑑𝑥−_ℰℎ⟦𝐯⟧∶𝝈𝑑𝑠. (3.11) Noting that [𝝈]=0 in ℰ^𝑖_ℎ and using integration by parts and the same technique for deriving (2.21), we deduce from (2.1) that𝑎_ℎ(𝐮,𝐯)=−_{𝐾∈𝒯ℎ𝐾}∇⋅𝝈𝑑𝑥=_Ω𝐟⋅𝐯𝑑𝑥. (3.12) Hence, we know that the LDG method (3.7) is not consistent with respect to the bilinear form 𝑎_ℎ(⋅,⋅), but it admits the following identity:𝑎_ℎ(𝐮−𝐮_ℎ,𝐯)=_ℰℎ⟦𝐯⟧∶{𝝈−𝐏_ℎ𝝈}𝑑𝑠+_{𝐾∈𝒯ℎ𝐾}(𝝈−𝐏_ℎ𝝈)∶𝐫(⟦𝐯⟧)𝑑𝑥,∀𝐯∈𝐕_ℎ. (3.13)

It is worth noting that the LDG method does not contain any parameter which needs to be quantified a priori. In what follows, we choose 𝐶₁₁=𝜂_𝑒ℎ_𝑒⁻¹ on each 𝑒∈ℰ_ℎ with {𝜂_𝑒}_{𝑒∈ℰℎ} having a uniform positive bound from above and below.

Next, we present a useful estimate for the local lifting operator 𝐫_𝑒. For this, let 𝐕(ℎ)∶=𝐕_ℎ+𝐻¹₀(Ω)^𝑑 and define the mesh-dependent energy norm for 𝐯∈𝐕(ℎ) by|‖𝐯‖|²=‖𝜺(𝐯)‖²_0,ℎ+_{𝑒∈ℰℎ}ℎ_𝑒⁻¹‖⟦𝐯⟧‖²_0,𝑒. (3.14)

Lemma 3.1. For any 𝐯∈𝐕(ℎ) and 𝑒∈ℰ_ℎ, one has ‖𝐫_𝑒(⟦𝐯⟧)‖²_0,ℎ≲ℎ_𝑒⁻¹‖⟦𝐯⟧‖²_0,𝑒. (3.15)

Proof. Since 𝐯∈𝐻¹₀(Ω)^𝑑 implies that ⟦𝐯⟧=0, it suffices to verify the result for 𝐯∈𝐕_ℎ. Taking 𝝓=⟦𝐯⟧ and 𝝉=𝐫_𝑒(⟦𝐯⟧)in (3.2), we get ‖𝐫_𝑒(⟦𝐯⟧)‖²_0,ℎ=−_𝑒⟦𝐯⟧∶{𝐫_𝑒(⟦𝐯⟧)}𝑑𝑠≤‖⟦𝐯⟧‖_0,𝑒‖{𝐫_𝑒(⟦𝐯⟧)}‖_0,𝑒. (3.16) On the other hand, by the scaling arguments, the trace theorem, and the inverse inequality, ‖{𝐫_𝑒(⟦𝐯⟧)}‖²_0,𝑒≲ℎ_𝑒⁻¹‖𝐫_𝑒(⟦𝐯⟧)‖²_0,ℎ+ℎ_𝑒|𝐫_𝑒(⟦𝐯⟧)|²_1,ℎ≲ℎ_𝑒⁻¹‖𝐫_𝑒(⟦𝐯⟧)‖²_0,ℎ. (3.17) Therefore, (3.15) is a direct consequence of (3.16) and (3.17).

Lemma 3.2 (boundedness). For any (𝐰,𝐯)∈𝐕(ℎ)×𝐕(ℎ), it holds that 𝑎_ℎ(𝐰,𝐯)≲|‖𝐰‖||‖𝐯‖|. (3.18)

Proof. According to the Cauchy-Schwarz inequality and the basic properties of 𝐫 and 𝐫_𝑒, we have 𝑎_ℎ≲(𝐰,𝐯)_{𝐾∈𝒯ℎ}(_𝐾|𝜺(𝐰)+𝐫(⟦𝐰⟧)|²𝑑𝑥⋅_𝐾|𝜺(𝐯)+𝐫(⟦𝐯⟧)|²𝑑𝑥)^1/2+_ℰℎ𝜂_𝑒ℎ_𝑒⁻¹⟦𝐰⟧∶⟦𝐯⟧𝑑𝑠≲(‖𝜺(𝐰)‖²_0,ℎ+‖𝐫(⟦𝐰⟧)‖²_0,ℎ)^1/2⋅(‖𝜺(𝐯)‖²_0,ℎ+‖𝐫(⟦𝐯⟧)‖²_0,ℎ)^1/2+_ℰℎ𝜂_𝑒ℎ_𝑒⁻¹⟦𝐰⟧∶⟦𝐯⟧𝑑𝑠≲(‖𝜺(𝐰)‖²_0,ℎ+_{𝑒∈ℰℎ}‖𝐫_𝑒(⟦𝐰⟧)‖²_0,ℎ)^1/2(‖𝜺(𝐯)‖²_0,ℎ+_{𝑒∈ℰℎ}‖𝐫_𝑒(⟦𝐯⟧)‖²_0,ℎ)^1/2+_ℰℎ𝜂_𝑒ℎ_𝑒⁻¹⟦𝐰⟧∶⟦𝐯⟧𝑑𝑠, (3.19) which, together with (3.15), yields 𝑎_ℎ(𝐰,𝐯)≲|‖𝐰‖||‖𝐯‖|+_ℰℎ𝜂_𝑒ℎ_𝑒⁻¹⟦𝐰⟧∶⟦𝐯⟧𝑑𝑠≲|‖𝐰‖||‖𝐯‖|, (3.20) as required.

Lemma 3.3 (stability). For any 𝐯∈𝐕_ℎ, it holds that 𝑎_ℎ(𝐯,𝐯)≳|‖𝐯‖|². (3.21)

Proof. Using the Cauchy-Schwarz inequality and (3.15), we get 𝑎_ℎ(𝐯,𝐯)≥2𝜇‖𝜺(𝐯)+𝐫(⟦𝐯⟧)‖²_0,ℎ+_ℰℎ𝜂_𝑒ℎ_𝑒⁻¹|⟦𝐯⟧|²𝑑𝑠≥2𝜇‖𝜺(𝐯)‖²_0,ℎ+2𝜇‖𝐫(⟦𝐯⟧)‖²_0,ℎ+4𝜇_{𝐾∈𝒯ℎ𝐾}𝜺(𝐯)∶𝐫(⟦𝐯⟧)𝑑𝑥+𝜂_0ℰℎℎ_𝑒⁻¹|⟦𝐯⟧|²𝑑𝑠≥2𝜇(1−𝜖)‖𝜺(𝐯)‖²_0,ℎ1+2𝜇(1−𝜖)‖𝐫(⟦𝐯⟧)‖²_0,ℎ+𝜂_0ℰℎℎ_𝑒⁻¹|⟦𝐯⟧|²𝑑𝑠≥2𝜇(1−𝜖)‖𝜺(𝐯)‖²_0,ℎ+(2𝜇𝐶₁+𝜂₀−2𝜇𝐶₁𝜖)_ℰℎℎ_𝑒⁻¹|⟦𝐯⟧|²𝑑𝑠. (3.22) Here 𝜂₀=min_{𝑒∈ℰℎ}𝜂_𝑒, 0<𝜖<1, and 𝐶₁ is the constant in (3.15). Therefore (3.21) is true if we choose 𝜖 such that 2𝜇𝐶₁+𝜂₀−2𝜇𝐶₁/𝜖>0, that is, 2𝜇𝐶₁/(2𝜇𝐶₁+𝜂₀)<𝜖<1.

Now, we are in a position to give error analysis for the LDG method (3.7). The main idea of our derivation is based on the framework on error analysis of DG methods for second-order elliptic problems (cf. [5]). Let 𝐐_ℎ be an 𝐿²-projection operator from 𝐕 onto the finite-element space 𝐕_ℎ. Let 𝐏_ℎ denote the usual 𝐿²-orthogonal projection operator onto 𝚺_ℎ. For simplicity, we still write 𝐐_ℎ (resp., 𝐏_ℎ) for 𝐐_ℎ∣_𝐾 (resp., 𝐏_ℎ∣_𝐾). Using the standard scaling arguments (cf. [24]), we can easily obtain error estimates for the operators 𝐐_ℎ and 𝐏_ℎ.

Lemma 3.4. Let 𝐯∈𝐻^𝑠+2(𝐾)^𝑑, 𝝈∈(𝐻^𝑠+1(Ω))^𝑠_𝑑×𝑑 with 𝑠≥0. Then for all 𝐾∈𝒯_ℎ, one has ‖𝐯−𝐐_ℎ𝐯‖_0,𝐾+ℎ_𝐾|𝐯−𝐐_ℎ𝐯|_1,𝐾≲ℎ_𝐾^{min{𝑠+1,𝑘}+1}‖𝐯‖_𝑠+2,𝐾,‖𝐯−𝐐_ℎ𝐯‖_0,𝜕𝐾≲ℎ_𝐾^{min{𝑠+1,𝑘}+1/2}‖𝐯‖_𝑠+2,𝐾,‖𝝈−𝐏_ℎ𝝈‖_0,𝐾≲ℎ_𝐾^{min{𝑠+1,𝑙+1}}‖𝝈‖_𝑠+1,𝐾,ℎ_𝑒^1/2‖𝝈−𝐏_ℎ𝝈‖_0,𝑒≲ℎ_𝐾^{min{𝑠+1,𝑙+1}}‖𝝈‖_𝑠+1,𝐾. (3.23)

Theorem 3.5. Let (𝝈,𝐮)∈(𝐻^𝑠+1(Ω))^𝑠_𝑑×𝑑×𝐻^𝑠+2(Ω)^𝑑, where 𝑠 is a nonnegative integer and 𝐮 is the solution of (2.1). Let 𝐮_ℎ∈𝐕_ℎ be the solution of (3.7). Then one has ‖|𝐮−𝐮_ℎ|‖≲ℎ^{min{𝑠+1,𝑘,𝑙+1}}‖𝐮‖_𝑠+2. (3.24)

Proof. From the stability (3.21), the identity (3.13), the boundedness (3.18), and Lemmas 3.1 and 3.4, it follows that |‖𝐐_ℎ𝐮−𝐮_ℎ‖|²≲𝑎_ℎ(𝐐_ℎ𝐮−𝐮_ℎ,𝐐_ℎ𝐮−𝐮_ℎ)=𝑎_ℎ(𝐐_ℎ𝐮−𝐮,𝐐_ℎ𝐮−𝐮_ℎ)+𝑎_ℎ(𝐮−𝐮_ℎ,𝐐_ℎ𝐮−𝐮_ℎ)=𝑎_ℎ(𝐐_ℎ𝐮−𝐮,𝐐_ℎ𝐮−𝐮_ℎ)+_ℰℎ⟦𝐐_ℎ𝐮−𝐮_ℎ⟧∶{𝝈−𝐏_ℎ𝝈}𝑑𝑠+_{𝐾∈𝒯ℎ𝐾}(𝝈−𝐏_ℎ𝝈)∶𝐫(⟦𝐐_ℎ𝐮−𝐮_ℎ⟧)𝑑𝑥≲‖|𝐐_ℎ𝐮−𝐮|‖‖|𝐐_ℎ𝐮−𝐮_ℎ|‖+ℎ^{min{𝑠+1,𝑙+1}}‖𝐮‖_𝑠+2‖|𝐐_ℎ𝐮−𝐮_ℎ|‖, (3.25) that is, ‖|𝐐_ℎ𝐮−𝐮_ℎ|‖≲‖|𝐐_ℎ𝐮−𝐮|‖+ℎ^{min{𝑠+1,𝑙+1}}‖𝐮‖_𝑠+2. (3.26) Therefore, using the triangle inequality and Lemma 3.4, we further have ‖|𝐮−𝐮_ℎ|‖≤‖|𝐮−𝐐_ℎ𝐮|‖+‖|𝐐_ℎ𝐮−𝐮_ℎ|‖≲‖|𝐮−𝐐_ℎ𝐮|‖+ℎ^{min{𝑠+1,𝑙+1}}‖𝐮‖_𝑠+2,Ω≲ℎ^{min{𝑠+1,𝑘,𝑙+1}}‖𝐮‖_𝑠+2,Ω. (3.27)

Theorem 3.6. Let (𝝈,𝐮)∈(𝐻^𝑠+1(Ω))^𝑠_𝑑×𝑑×𝐻^𝑠+2(Ω)^𝑑, where 𝑠 is a nonnegative integer and 𝐮 is the solution of (2.1). Let 𝐮_ℎ∈𝐕_ℎ be the solution of (3.7). Then one has ‖𝐮−𝐮_ℎ‖_1,ℎ≲ℎ^{min{𝑠+1,𝑘,𝑙+1}}‖𝐮‖_𝑠+2. (3.28)

Proof. By the triangle inequality and Korn's inequality over 𝐕_ℎ (cf. (2.34)), we know that ‖𝐮−𝐮_ℎ‖_1,ℎ≲‖𝐮−𝐐_ℎ𝐮‖_1,ℎ+‖𝐐_ℎ𝐮−𝐮_ℎ‖_1,ℎ≲‖𝐮−𝐐_ℎ𝐮‖_1,ℎ+|‖𝐐_ℎ𝐮−𝐮_ℎ‖|≲‖𝐮−𝐐_ℎ𝐮‖_1,ℎ+|‖𝐮−𝐐_ℎ𝐮‖|+|‖𝐮−𝐮_ℎ‖|. (3.29) The desired result then follows from the above estimate combined with Lemma 3.4 and Theorem 3.5.

Theorem 3.7. Suppose that Ω is a convex bounded polygon or polyhedron. Let (𝝈,𝐮)∈(𝐻^𝑠+1(Ω))^𝑠_𝑑×𝑑×𝐻^𝑠+2(Ω)^𝑑, where 𝑠 is a nonnegative integer and 𝐮 is the solution of (2.1). Let 𝐮_ℎ∈𝐕_ℎ be the solution of (3.7). Then one has ‖𝐮−𝐮_ℎ‖₀≲ℎ^{min{1,𝑘}+min{𝑠+1,𝑘,𝑙+1}}‖𝐮‖_𝑠+2. (3.30)

Proof. The proof relies on the usual duality argument. Let (𝝈,𝐮) be the solution of the auxiliary problem 𝒜𝝈−𝜺(𝐮)=𝟎inΩ,−∇⋅𝝈=𝐮−𝐮_ℎinΩ,𝐮=𝟎on𝜕Ω. (3.31) Formally, (3.31) is problem (2.1) with 𝐟 replaced by 𝐮−𝐮_ℎ. Since Ω is a convex bounded polygon or polyhedron, we have an 𝐻²-regularity estimate (cf. [28]) ‖𝐮‖₂≲‖𝐮−𝐮_ℎ‖₀. (3.32)
Therefore, applying an elementwise integration by parts to the second equation of (3.31) with a test function 𝐮−𝐮_ℎ and using the definitions of 𝐫 (cf. (3.1)) and 𝑎_ℎ(⋅,⋅), the first equation of (3.31), and the technique to derive (2.21), we get‖𝐮−𝐮_ℎ‖²₀=−_{𝐾∈𝒯ℎ𝐾}(𝐮−𝐮_ℎ)⋅(∇⋅𝝈)𝑑𝑥=_{𝐾∈𝒯ℎ𝐾}𝝈∶𝜺(𝐮−𝐮_ℎ)𝑑𝑥−_ℰℎ𝝈∶⟦𝐮−𝐮_ℎ⟧𝑑𝑠=𝑎_ℎ(𝐮,𝐮−𝐮_ℎ)−_{𝐾∈𝒯ℎ𝐾}𝝈∶𝑟(⟦𝐮−𝐮_ℎ⟧)𝑑𝑥−_ℰℎ𝝈∶⟦𝐮−𝐮_ℎ⟧𝑑𝑠=𝑎_ℎ(𝐮−𝐮_ℎ,𝐮)−_ℰℎ{𝝈−𝐏_ℎ𝝈}∶⟦𝐮−𝐮_ℎ⟧𝑑𝑠−_{𝐾∈𝒯ℎ𝐾}(𝝈−𝐏_ℎ𝝈)∶𝐫(⟦𝐮−𝐮_ℎ⟧)𝑑𝑥, (3.33) and with the identity (3.13), we further have ‖𝐮−𝐮_ℎ‖²₀=𝑎_ℎ(𝐮−𝐮_ℎ,𝐮−𝐐_ℎ𝐮)+𝑎_ℎ(𝐮−𝐮_ℎ,𝐐_ℎ𝐮)−_ℰℎ{𝝈−𝐏_ℎ𝝈}∶⟦𝐮−𝐮_ℎ⟧𝑑𝑠−_{𝐾∈𝒯ℎ𝐾}(𝝈−𝐏_ℎ𝝈)∶𝐫(⟦𝐮−𝐮_ℎ⟧)𝑑𝑥=𝑎_ℎ(𝐮−𝐮_ℎ,𝐮−𝐐_ℎ𝐮)+_ℰℎ⟦𝐐_ℎ𝐮−𝐮⟧∶{𝝈−𝐏_ℎ𝝈}𝑑𝑠+_{𝐾∈𝒯ℎ𝐾}(𝝈−𝐏_ℎ𝝈)∶𝐫(⟦𝐐_ℎ−𝐮−𝐮⟧)𝑑𝑥_ℰℎ{𝝈−𝐏_ℎ𝝈}∶⟦𝐮−𝐮_ℎ⟧𝑑𝑠−_{𝐾∈𝒯ℎ𝐾}(𝝈−𝐏_ℎ𝝈)∶𝐫(⟦𝐮−𝐮_ℎ⟧)𝑑𝑥. (3.34) Therefore, it follows from the boundedness (3.18), the regularity (3.31), and Lemmas 3.1 and 3.4 that ‖𝐮−𝐮_ℎ‖²₀≲|‖𝐮−𝐮_ℎ‖||‖𝐮−𝐐_ℎ𝐮‖|+ℎ|‖𝐮−𝐮_ℎ‖|‖𝝈‖₁+ℎ^{min{𝑠+1,𝑙+1}}‖𝐮‖_𝑠+2(_{𝑒∈ℰℎ}ℎ_𝑒⁻¹‖⟦𝐐_ℎ𝐮−𝐮⟧‖²_0,𝑒)^1/2≲ℎ^min{1,𝑘}(|‖𝐮−𝐮_ℎ‖|+ℎ^{min{𝑠+1,𝑙+1}}‖𝐮‖_𝑠+2)‖𝐮‖₂≲ℎ^min{1,𝑘}(|‖𝐮−𝐮_ℎ‖|+ℎ^{min{𝑠+1,𝑙+1}}‖𝐮‖_𝑠+2)‖𝐮−𝐮_ℎ‖₀. (3.35) This, along with Theorem 3.5, immediately leads to ‖𝐮−𝐮_ℎ‖₀≲ℎ^min{1,𝑘}(|‖𝐮−𝐮_ℎ‖|+ℎ^{min{𝑠+1,𝑙+1}}‖𝐮‖_𝑠+2)≲ℎ^{min{1,𝑘}+min{𝑠+1,𝑘,𝑙+1}}‖𝐮‖_𝑠+2. (3.36)

We would like to end this section with some results on the optimality of our estimates derived. If we choose 𝑆₁(𝐾) and 𝑆₂(𝐾) to be 𝑃_𝑙(𝐾) and 𝑃_𝑘(𝐾), respectively, it is easy to check that the first condition of (2.5) implies 𝑘−1≤𝑙. Therefore, we can obtain the following result from the previous theorems directly.

Corollary 3.8. Let (𝝈,𝐮)∈(𝐻^𝑠+1(Ω))^𝑠_𝑑×𝑑×𝐻^𝑠+2(Ω)^𝑑, where 𝑠 is a nonnegative integer and 𝐮 is the solution of (2.1). Let 𝐮_ℎ∈𝐕_ℎ be the solution of (3.7) with 𝑆₁(𝐾)=𝑃_𝑙(𝐾) and 𝑆₂(𝐾)=𝑃_𝑘(𝐾) for all 𝐾∈𝒯_ℎ. Then one has the following optimal error estimates: ‖|𝐮−𝐮_ℎ|‖≲ℎ^{min{𝑠+1,𝑘}}‖𝐮‖_𝑠+2,‖𝐮−𝐮_ℎ‖_1,ℎ≲ℎ^{min{𝑠+1,𝑘}}‖𝐮‖_𝑠+2. (3.37) In addition, if Ω is a convex bounded polygon or polyhedron, ‖𝐮−𝐮_ℎ‖₀≲ℎ^{min{1,𝑘}+min{𝑠+1,𝑘}}‖𝐮‖_𝑠+2. (3.38)

4. Numerical Results

In this section, we report a numerical example to illustrate the theoretical results. Let Ω=(−1,1)×(−1,1), 𝜆=0.3, 𝜇=0.35, and 𝐟(𝑥₁,𝑥₂)=2𝜇(3−𝑥²₁−2𝑥²₂−2𝑥₁𝑥₂,3−2𝑥²₁−𝑥²₂−2𝑥₁𝑥₂)^𝑡+𝜆(2−2𝑥²₂−4𝑥₁𝑥₂,2−2𝑥²₁−4𝑥₁𝑥₂)^𝑡. (4.1) It is not difficult to check that the exact solution of (2.1) is𝐮(𝑥₁,𝑥₂)=((1−𝑥²₁)(1−𝑥²₂),(1−𝑥²₁)(1−𝑥²₂))^𝑡. (4.2) We use a quasiuniform triangulation 𝒯_ℎ over Ω. For any 𝐾∈𝒯_ℎ, we take 𝑆₁(𝐾)=𝑃_𝑙(𝐾) and 𝒮₂(𝐾)=𝑃_𝑘(𝐾) where 𝑘,𝑙≥0. On each edge 𝑒, 𝜂_𝑒 is set to be 1. The numerical results of the LDG method for a few choices of 𝑘 and 𝑙 are shown in Table 1 and Figures 1, 2, and 3. We observe that, when 𝑘=1, the numerical convergence rates of ‖𝐮−𝐮_ℎ‖₀, |𝐮−𝐮_ℎ|_1,ℎ, and |‖𝐮−𝐮_ℎ‖| are 𝑂(ℎ²), 𝑂(ℎ), and 𝑂(ℎ), respectively. However, there is no convergence when 𝑘=0. These phenomena agree with theoretical results in Theorems 3.5–3.7. When 𝑘=1, accuracies of the numerical results are nearly the same as those for 𝑙=0 and 𝑙=1. It goes without saying that it is more convenient to simulate for 𝑙=0 than for 𝑙=1.

Table 1 

Numerical errors of LDG method for several choices of 𝑘 and 𝑙. The theoretical convergence rates of ‖𝐮−𝐮_ℎ‖₀, |𝐮−𝐮_ℎ|_1,ℎ, and ‖|𝐮−𝐮_ℎ|‖ are 𝑂(ℎ^{𝑘+min{1,𝑘}}), 𝑂(ℎ^𝑘), and 𝑂(ℎ^𝑘), respectively.

Figure 1 

$\parallel \mathbf{u}-{\mathbf{u}}_h{\parallel }_0$ in the $\ln \hspace{0.17em}$-$\ln \hspace{0.17em}$ scale for several choices of $k$ and $l$.

Figure 2 

$|\mathbf{u}-{\mathbf{u}}_h{|}_{1,h}$ in the $\ln \hspace{0.17em}$-$\ln \hspace{0.17em}$ scale for several choices of $k$ and $l$.

Figure 3 

$|\parallel \mathbf{u}-{\mathbf{u}}_h\parallel |$ in the $\ln \hspace{0.17em}$-$\ln \hspace{0.17em}$ scale for several choices of $k$ and $l$.

Acknowledgments

The authors thank the referee for the valuable comments and suggestions, which greatly improved the presentation of the manuscript. The research of the second author is supported partially by the NNSFC (10371076) and E-Institutes of Shanghai Municipal Education Commission (E03004).