Superconvergence analysis of DG-FEM based on the polynomial preserving recovery for Helmholtz equation with high wave number
SSUPERCONVERGENCE ANALYSIS OF DG-FEM BASED ON THEPOLYNOMIAL PRESERVING RECOVERY FORHELMHOLTZ EQUATION WITH HIGH WAVE NUMBER
YU DU ∗∗ AND
ZHIMIN ZHANG †† Abstract.
We study superconvergence property of the linear discontinuous Galerkin finiteelement method with the polynomial preserving recovery (PPR) and Richardson extrapolation forthe two dimensional Helmholtz equation. The error estimate with explicit dependence on the wavenumber k , the penalty parameter µ and the mesh condition parameter α is derived. First, weprove that under the assumption k ( kh ) ≤ C ( h is the mesh size) and certain mesh condition,the estimate between the finite element solution and the linear interpolation of the exact solutionis superconvergent under the (cid:107)|·|(cid:107) -seminorm. Second, we prove a superconvergence result for therecovered gradient by PPR. Furthermore, we estimate the error between the finite element gradientand recovered gradient, which motivate us to define the a posteriori error estimator. Finally, Somenumerical examples are provided to confirm the theoretical results of superconvergence analysis. Alltheoretical findings are verified by numerical tests. Key words.
Helmholtz equation, large wave number, pollution errors, superconvergence, poly-nomial preserving recovery, discontinuous Galerkin finite element methods
AMS subject classifications.
1. Introduction.
This paper is devoted to superconvergence analysis of thelinear Galerkin finite element method for the following Helmholtz problem: − ∆ u − k u = f in Ω , (1.1) ∂u∂n + i ku = g on Γ , (1.2)where Ω ∈ R is a bounded polygon with boundary Γ := ∂ Ω, i = √− n denotes the unit outward normal to Γ. The above Helmholtzproblem is an approximation of the acoustic scattering problem (with time dependence e i ωt ) and k is known as the wave number. The Robin boundary condition (1.2) isknown as the first order approximation of the radiation condition (cf. [16]). TheHelmholtz problem (1.1)–(1.2) also arises in applications as a consequence of frequencydomain treatment of attenuated scalar waves (cf. [13]).It is well known that the finite element method of fixed order for the Helmholtzproblem (1.1)–(1.2) at high frequencies ( k (cid:29)
1) is subject to the effect of pollution: theratio of the error of the finite element solution to the error of the best approximationfrom the finite element space cannot be uniformly bounded with respect to k [2, 5, 4,12, 19, 20, 21]. More precisely, the linear finite element method for a 2-D Helmholtzproblem satisfies the following error estimate under the mesh constraint k ( kh ) ≤ C ∗ Beijing computational science research center, Beijing, 100193, China. [email protected],[email protected] . This research work is supported by a Tianhe–2JK computing time award at theBeijing Computational Science Research Center (CSRC). The research of this author was supportedin part by the China Postdoctoral Science Foundation under grant 2016M591053 and the NationalNatural Science Foundation of China under grants 11601026 † Beijing Computational Science Research Center, Beijing, 100193 and Department of Mathemat-ics, Wayne State University, Detroit, MI 48202. [email protected], [email protected] . Theresearch of this author was supported in part by the National Natural Science Foundation of Chinaunder grants 11471031, 91430216, U1530401, and the U.S. National Science Foundation throughgrant DMS–1419040. 1 a r X i v : . [ m a t h . NA ] D ec Y. Du & Z. Zhang [37, 14]: (cid:107)∇ ( u − u h ) (cid:107) L (Ω) ≤ C kh + C k ( kh ) . (1.3)Here u h is the linear finite element solution, h is the mesh size and C i , i = 1 , k and h . It is easy to see that the order of the firstterm on the right hand side of (1.3) is the same to that of the interpolation errorin H -seminorm and it can dominate the error bound only if k ( kh ) is small enough.However, the second term on the right-hand side of (1.3) dominates the estimateunder other mesh conditions. For example, kh is fixed and k is large enough. Theterm C k ( kh ) is called the pollution error of the finite element solution.Considerable efforts have been made in analysis of different numerical methodsfor the Helmholtz problem with large wave number in the literature. The readersare referred to [3, 13, 28] for asymptotic error estimates of general DG methods and[20, 21] for pre-asymptotic error estimates of a one-dimensional problem discretized onequidistant grid. For more pre-asymptotic error estimates, Please refer to [24, 25] and[8, 37] for classical finite element methods as well as interior penalty finite elementmethods. For other methods solving the Helmholtz problems, such as the interiorpenalty discontinuous Galerkin method or the source transfer domain decompositionmethod, one can read [23, 17, 18, 36, 15, 10].In this work, we investigate the superconvergence property of the linear discon-tinuous Galerkin (DG) finite element method when being post-processed by our poly-nomial preserving recovery (PPR) for the Helmholtz problem. PPR was proposedby Zhang and Naga [35] in 2004 and has been successfully applied to finite elementmethods. COMSOL Multiphysics adopted PPR as a post-processing tool since 2008[1]. [ ? ] has applied the technique to the Helmholtz problem and prove its superconver-gence property. In this paper, we generalize the technique over the DG finite elementspace and prove its superconvergence property. To learn more about PPR, readersare referred to [33, 32, 26, 29]. Some theoretical results about recovery techniquesand recovery-type error estimators can be found in [6, 22, 34, 30, 31].Our purpose of this paper is to prove the superconvergence error estimates forthe linear discontinuous Galerkin finite element method and analyze the influenceof the PPR technique on the pollution error. Note the superconvergence error es-timates depend on the triangulation, the penalty parameter and the wave numberunder certain mesh condition. In order to prove the estimates, we first assume somemesh constraints, called “Condition α ”, and then redefine the PPR method on thediscontinuous Galerkin finite element space. Finally, all the estimates motivate usto combine the PPR technique and the Richardson extrapolation to reduce the errorfurther and define the a posterior error estimator.The remainder of this paper is organized as follows: some notations, DG-FEMand the mesh constraints are introduced in section 2. In section 3, we prove the super-convergence between the interpolant and the finite element solution to the problemwith Robin boundary (1.1)–(1.2). In section 4, we redefine the PPR technique overthe linear DG finite element space and prove the superconvergence property of G h inthe Sobolev space H and show the most important result, that is the error estimateof G h u h . Then we try to give the reason for the effect of G h to the pollution error insection 5. Finally, we simulate a model problem by the linear DG-FEM, PPR methodand the Richardson extrapolation in section 6. It is shown that the recovered gradientcan be improved by the Richardson extrapolation further and the a posterior errorestimator based on the PPR and Richardson extrapolation is exact asymptotically. uperconvergence analysis of linear FEM x Ω ∈ Ω and a positive constant c Ω depending only on Ω such that( x − x Ω ) · n ≥ c Ω ∀ x ∈ Γ .
2. Preliminaries.
Throughout this paper, we assume that for any node point z ∈ N h , there exists at least one interior edge e ∈ E Ih having z .To introduce the method and simplify the analysis, we introduce some notationfirst. The standard Sobolev and Hilbert space, norm, and inner product notation areadopted (cf. [7, 11]). In particular, ( · , · ) Q and (cid:104)· , ·(cid:105) Σ for Σ = ∂Q denote the L -innerproduct on complex-valued L ( Q ) and L (Σ) spaces, respectively.Let T h be a regular triangulation of the domain Ω. For any τ ∈ T h , we denote by h τ its diameter and by | τ | its area. Let h = max τ ∈T h h τ . Assume that h τ (cid:104) h .Let V h be the approximation space of piecewise linear polynomials, that is, V h := (cid:8) v h ∈ L (Ω) : v h | τ ∈ P ( τ ) ∀ τ ∈ T h (cid:9) , where P ( τ ) denotes the set of all polynomials on τ with degree ≤ E h be the set of all edges of T h and N h be the set of all nodal points. Denote allthe boundary edges by E Bh := { e ∈ E h : e ⊂ Γ } and the interior edges by E Ih := E h \E Bh .For each edge e ∈ E h , define h e := diam( e ). For e = ∂τ ∩ τ (cid:48) ∈ E Ih , let n e be a unitnormal vector to e . We assume that the normal vector n e is oriented from τ (cid:48) to τ and define [ v ] | e := v | τ (cid:48) − v | τ , { v } := 12 ( v | τ (cid:48) + v | τ ) . We define the space E := (cid:81) τ ∈T h H ( τ ) and introduce the sesquilinear form a h ( · , · )on E × E as follows: a h ( u, v ) := (cid:88) K ∈T h ( ∇ u, ∇ v ) K − (cid:88) e ∈E Ih (cid:18)(cid:28)(cid:26) ∂u∂ n e (cid:27) , [ v ] (cid:29) e + (cid:28) [ u ] , (cid:26) ∂v∂ n e (cid:27)(cid:29) e (cid:19) + J ( u, v ) J ( u, v ) := (cid:88) e ∈E Ih ρ ,e h µe (cid:104) [ u ] , [ v ] (cid:105) e , The linear DG method is defined as follows: find u h ∈ V h such that a h ( u h , v h ) − k ( u h , v h ) + i k (cid:104) u h , v h (cid:105) = ( f, v h ) + (cid:104) g, v h (cid:105) ∀ v h ∈ V h . (2.1) Remark 2.1 (a) Note that the method is the standard symmetric DG method(cf. [27]). So we have the proposition (cf. proposition 2.9 in [27]): u ∈ H (Ω) is thesolution to (1.1)–(1.2) if and only if u satisfies the general DG variational formulation a h ( u, v ) − k ( u, v ) + i k (cid:104) u, v (cid:105) = ( f, v ) + (cid:104) g, v (cid:105) ∀ v ∈ E. (b) A similar DG method, called interior penalty discontinuous Galerkin method(IPDG), was introduced and analyzed by Feng, Wu and so on for the Helmholtz Y. Du & Z. Zhang problem (1.1)–(1.2). The reader is referred to [17, 18, 36, 15] for both asymptotic andpreasymptotic error estimates.The following two norms will be used in the forthcoming sections : (cid:107) v (cid:107) ,h := (cid:32) (cid:88) K ∈T h (cid:107) v (cid:107) L ( K ) + J ( v, v ) (cid:33) / ∀ v ∈ E, (cid:107)| v |(cid:107) ,h := (cid:16) (cid:107) v (cid:107) ,h + k (cid:107) v (cid:107) (cid:17) / ∀ v ∈ E. We denote by |·| H ( ω ) the broken semi- H norm (cid:16)(cid:80) τ ∈ ω |·| H ( τ ) (cid:17) / where ω ⊂ T h .Throughout the paper C denotes a generic positive constant which is independentof h, k, f, g and the penalty parameters. We use the shorthand notation A (cid:46) B for A ≤ CB and assume k (cid:29) u to the problem (1.1)–(1.2) is H -regular over Ω and thedata g is H -regular over Γ. Denote by C u,g = (cid:88) j =1 k − ( j − (cid:107) u (cid:107) j + (cid:88) j =1 k − j | g | H j (Γ) . (2.2)The function C u,g could be treated as a constant in this paper since (cid:107) u (cid:107) j isbounded by max( k , k j − ) when u is the solution to the Helmholtz problem (1.1)–(1.2). The reader is referred to [23, 24, 25] for the estimates of u .Before estimating the errors, we state the coercivity and continuity properties forthe sesquilinear form a h ( · , · ). Since they easily follow from the difinitions 2.10–2.11in [27], the details are omitted. Lemma 2.1.
For any v, w ∈ E , the sesquilinear form a h ( · , · ) satisfies | a h ( v, w ) | , | a h ( w, v ) | (cid:46) (cid:107) w (cid:107) ,h (cid:107) v (cid:107) ,h . In addition, there exists constant ρ such that if ρ ≤ ρ ,e , (cid:107) v h (cid:107) ,h (cid:46) a h ( v h , v h ) ∀ v h ∈ V h , (cid:107)| v h |(cid:107) ,h (cid:46) a h ( v h , v h ) + k ( v h , v h ) ∀ v h ∈ V h . The following lemma shows the preasymptotic error estimates for the solution u to(1.1)–(1.2). The results can be derived by arguments same to those in [15, 36] andwe omit the details to save space. Lemma 2.2.
Assume that u is the solution the problem (1.1) – (1.2) and u h is thediscrete solution of the scheme (2.1) . Then there exists constant C independent of k and h , such that if k ( kh ) ≤ C then the following estimates hold: (cid:107) u − u h (cid:107) ,h (cid:46) (cid:0) kh + k ( kh ) (cid:1) C u,g , (2.3) k (cid:107) u − u h (cid:107) (cid:46) (cid:0) ( kh ) + k ( kh ) (cid:1) C u,g . (2.4)We begin with some definitions regarding meshes. For an interior edge e ∈ E Ih ,we denote Ω e = τ e ∪ τ (cid:48) e , a patch formed by the two elements τ e and τ (cid:48) e sharing e ,see Figures 2.1-2.2. For any edge e ∈ E h and an element τ with e ⊂ τ , θ e denotes uperconvergence analysis of linear FEM (cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:64)(cid:64)(cid:64)(cid:64)(cid:64)(cid:64)(cid:64)(cid:64) eτ τ (cid:48) e + 1 e − θ e θ (cid:48) e (cid:0)(cid:0)(cid:9) (cid:0)(cid:0)(cid:18) n (cid:48) e n e Fig. 2.1 . Notation in the patch Ω e . (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0) Γ e e − e + 1 (cid:64)(cid:64)(cid:73) θ e e e − e + 1 (cid:64)(cid:64)(cid:73) θ e e e − e + 1 (cid:64)(cid:64)(cid:73) θ e (cid:64)(cid:64)(cid:64)(cid:64)(cid:64)(cid:64) Γ (cid:12)(cid:12)(cid:12)(cid:12) (cid:76)(cid:76)(cid:76)(cid:76) e e + 1 e − (cid:0)(cid:0)(cid:9) θ e (cid:8)(cid:8)(cid:8) e e − e + 1 (cid:64)(cid:64)(cid:82) θ e Fig. 2.2 . Notation in the boundary elements. the angle opposite of the edge e in τ , t e denotes the unit tangent vector of e withcounterclockwise orientation and n e , the unit outward normal vector of e , h e , h e +1 ,and h e − denote the lengths of the three edges of τ , respectively. Here the subscript e + 1 or e − (cid:48) is added for the corresponding quantitiesin τ (cid:48) with t e = − t (cid:48) e and n e = − n (cid:48) e due to the orientation.For any e ∈ E Ih (cf. Figure 2.1), we say that Ω e is an ε approximate parallelogramif the lengths of any two opposite edges differ by at most ε , that is, (cid:12)(cid:12) h e − − h (cid:48) e − (cid:12)(cid:12) + (cid:12)(cid:12) h e +1 − h (cid:48) e +1 (cid:12)(cid:12) ≤ ε. For any e ∈ E Bh (cf. Figure 2.2), we say that τ e is an ε approximate isoscelestriangle if the lengths of its two edges e − e + 1 differ by at most ε , that is, | h e +1 − h e − | ≤ ε. Definition 2.3.
The triangulation T h is said to satisfy mesh condition α if thereexists a constant α ≥ such that(a) the patch Ω e is an O ( h α ) approximate parallelogram for any interior edge e ∈ E Ih ;(b) the triangle τ e is an O ( h α ) approximate isosceles triangle for any boundaryedge e ∈ E Bh ;Remark α ) in Definition 2.3 is often used to prove the su-perconvergence property for problems with the Dirichlet boundary condition [9, 29].Note that this restriction is technique and just for theoretical purpose. In fact, one su-perconvergence results still can be obtained under general meshes which do not satisfythe condition, such as Delaunay triangulation and Chevron pattern triangulation. Y. Du & Z. Zhang
3. Superconvergence between the discontinuous finite element solutionand linear interpolant.
First we introduce a quadratic interpolant ψ Q = Π Q ψ of ψ based on nodal values and moment conditions on edges,(Π Q φ )( z ) = φ ( z ) , (cid:90) e Π Q φ = (cid:90) e φ ∀ z ∈ N h , e ∈ E h . (3.1)The following fundamental identity for v h ∈ P ( τ ) has been proved in [9]: (cid:90) τ ∇ ( φ − φ I ) · ∇ v h = (cid:88) e ∈ ∂τ (cid:18) β e (cid:90) e ∂ φ Q ∂ t e ∂v h ∂ t e + γ e (cid:90) e ∂ φ Q ∂ t e ∂ n e ∂v h ∂ t e (cid:19) (3.2)where β e = 112 cot θ e ( h e +1 − h e − ) , γ e = 13 cot θ e | τ | , (3.3)and φ I ∈ P ( τ ) is the linear interpolant of φ on τ . The following lemma can be easilyobtained [9, 29]. Lemma 3.1.
Let m e denote t e or n e . Assume that T h satisfies the mesh condition α , then we have the following estimates:(a) For any interior edge e ∈ E Ih , | β e | + | β (cid:48) e | (cid:46) h , | γ e | + | γ (cid:48) e | (cid:46) h ;(3.4) | β e − β (cid:48) e | (cid:46) h α , | γ e − γ (cid:48) e | (cid:46) h α . (3.5) (b) For two adjacent edges e , e ∈ E Bh , that is e ∩ e (cid:54) = ∅ , | β e | + | β e | (cid:46) h α , | γ e | + | γ e | (cid:46) h (3.6) | γ e − γ e | (cid:46) h α (3.7) (c) For any edge e ∈ E h , e ⊂ ∂τ e , (cid:90) e ∂ φ∂ t e ∂ m e ∂v h ∂ t e (cid:46) ( (cid:107) φ (cid:107) H ( τ e ) + h − (cid:107) φ (cid:107) H ( τ e ) ) | v h | H ( τ e ) ;(3.8) (cid:90) e ∂ ( φ − φ Q ) ∂ t e ∂ m e ∂v h ∂ t e (cid:46) | φ | H ( τ e ) | v h | H ( τ e ) . (3.9) Proof . The inequalities (3.4)–(3.6) follow from the mesh condition α . From thecondition (a) and (b) in Definition 2.3, we have for any e , e ∈ E Bh satisfying e ∩ e (cid:54) = ∅ (cf. Figure 2.2), (cid:12)(cid:12) h e − h e +1 cos θ e − h e − h e +1 cos θ e (cid:12)(cid:12) h (cid:46) h α , which implies (3.7).Finally, the inequalities (3.8) and (3.9) follow from the trace theorem. Lemma 3.2.
Assume that T h satisfies the α approximation ondition . Then forany v h ∈ V h , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) K ∈T h (cid:90) K ∇ ( u − u I ) · ∇ v h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:46) (cid:16) ( kh ) + kh α + kh µ/ (cid:17) (cid:107)| v h |(cid:107) ,h C u,g . (3.10) uperconvergence analysis of linear FEM Here u I is the linear interpolant of u on Ω .Proof . From (3.2), we have (cid:90) Ω ∇ ( u − u I ) · ∇ v h = (cid:88) τ ∈T h (cid:88) e ⊂ ∂τ (cid:18) β e (cid:90) e ∂ u Q ∂ t e ∂v h ∂ t e + γ e (cid:90) e ∂ u Q ∂ t e ∂ n e ∂v h ∂ t e (cid:19) = I + I , where I = (cid:88) e ∈E Ih (cid:20)(cid:90) e ∂ u∂ t e (cid:18) β e ∂v h ∂ t e | K − β (cid:48) e ∂v h ∂ t e | K (cid:19) + (cid:90) e ∂ u∂ t e ∂ n e (cid:18) γ e ∂v h ∂ t e | K − γ (cid:48) e ∂v h ∂ t e | K (cid:19) + β e (cid:90) e ∂ ( u Q − u ) ∂ t e ∂v h ∂ t e + γ e (cid:90) e ∂ ( u Q − u ) ∂ t e ∂ n e ∂v h ∂ t e + β (cid:48) e (cid:90) e ∂ ( u − u Q ) ∂ t e ∂v h ∂ t e + γ (cid:48) e (cid:90) e ∂ ( u − u Q ) ∂ t e ∂ n e ∂v h ∂ t e (cid:21) , = I , + I , + I , ,I = (cid:88) e ∈E Bh (cid:20) β e (cid:90) e ∂ u∂ t e ∂v h ∂ t e + γ e (cid:90) e ∂ u∂ t e ∂ n e ∂v h ∂ t e + β e (cid:90) e ∂ ( u Q − u ) ∂ t e ∂v h ∂ t e + γ e (cid:90) e ∂ ( u Q − u ) ∂ t e ∂ n e ∂v h ∂ t e (cid:21) . First, I , and I , can be estimated by Lemma 3.1 and H¨older’s inequality: | I , + I , | (cid:46) (cid:88) e ∈E Ih (cid:16) ( h α + h ) (cid:107) u (cid:107) H ( τ e ) + h α (cid:107) u (cid:107) H ( τ e ) (cid:17) | v h | H ( τ e ) (cid:46) (cid:0) ( h α + h ) (cid:107) u (cid:107) + h α (cid:107) u (cid:107) (cid:1) | v h | H (cid:46) (cid:0) ( kh ) + kh α (cid:1) | v h | H C u,g . From Lemma 3.1 and the inverse inequality, I , = (cid:88) e ∈E Ih (cid:20) ( β e − β (cid:48) e ) (cid:90) e ∂ u∂ t e ∂v h ∂ t e | K + ( γ e − γ (cid:48) e ) (cid:90) e ∂ u∂ t e ∂ n e ∂v h ∂ t e | K + β (cid:48) e (cid:90) e ∂ u∂ t e (cid:20) ∂v h ∂ t e (cid:21) + γ e (cid:90) e ∂ u∂ t e ∂ n e (cid:20) ∂v h ∂ t e (cid:21)(cid:21) (cid:46) (cid:88) e ∈E Ih (cid:20) (cid:16) h α (cid:107) u (cid:107) H ( τ e ) + h α (cid:107) u (cid:107) H ( τ e ) (cid:17) | v h | H ( τ e ) + (cid:16) h / (cid:107) u (cid:107) H ( τ e ) + h / (cid:107) u (cid:107) H ( τ e ) (cid:17) (cid:18) (cid:88) e ∈E Ih (cid:13)(cid:13)(cid:13)(cid:13)(cid:20) ∂v h ∂ t e (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) L ( e ) (cid:19) / (cid:21) (cid:46) (cid:16) ( kh ) + kh α + kh µ/ (cid:17) (cid:107) v h (cid:107) ,h C u,g , where we have used the ineqaulities (cid:18) (cid:88) e ∈E Ih (cid:13)(cid:13)(cid:13)(cid:13)(cid:20) ∂v h ∂ t e (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) L ( e ) (cid:19) / = (cid:18) (cid:88) e ∈E Ih (cid:13)(cid:13)(cid:13)(cid:13) ∂ [ v h ] ∂ t e (cid:13)(cid:13)(cid:13)(cid:13) L ( e ) (cid:19) / (cid:46) (cid:18) (cid:88) e ∈E Ih h − e (cid:107) [ v h ] (cid:107) L ( e ) (cid:19) / (cid:46) h µ/ − / J ( v h , v h ) / ≤ h µ/ − / (cid:107) v h (cid:107) ,h . Y. Du & Z. Zhang
By combining the inequalities above, we get | I | (cid:46) (cid:16) ( kh ) + kh α + kh µ/ (cid:17) (cid:107) v h (cid:107) ,h C u,g . (3.11)Then we turn to the estimate of I . From (3.6) and (3.9), (cid:88) e ∈E Bh (cid:20) β e (cid:90) e ∂ u∂t e ∂v h ∂t e + β e (cid:90) e ∂ ( u Q − u ) ∂t e ∂v h ∂t e + γ e (cid:90) e ∂ ( u Q − u ) ∂t e ∂n e ∂v h ∂t e (cid:21) (3.12) (cid:46) (cid:88) e ∈E Bh (cid:0) h α (cid:107) u (cid:107) H ( τ e ) + ( h α + h ) (cid:107) u (cid:107) H ( τ e ) (cid:1) | v h | H ( τ e ) (cid:46) (cid:0) kh α + ( kh ) (cid:1) | v h | H C u,g . Therefore, we only need to estimate the remaining terms of I . Denote by z i thenodes on Γ. Denote by e i , e i +1 ∈ E Bh sharing z i with counterclockwise orientation (cf.Figure 2.2). Denote by [ · ] z i = ( ·| e i +1 )( z i ) − ( ·| e i )( z i ) and {·} z i = (cid:0) ·| e i ( z i )+ ·| e i +1 ( z i ) (cid:1) / (cid:88) e ∈E Bh γ e (cid:90) e ∂ u∂t e ∂n e ∂v h ∂t e = (cid:88) z i ∈ Γ (cid:84) N h [ γ e v h ] z i ∂ u∂t e ∂n e ( z i ) − (cid:88) e ∈E Bh γ e (cid:90) e ∂ u∂t e ∂n e v h , (3.13) = I , + I , + I , , where I , = (cid:88) z i ∈ Γ (cid:84) N h [ γ e ] z i { v h } z i ∂ u∂t e ∂n e ( z i ) ,I , = (cid:88) z i ∈ Γ (cid:84) N h { γ e } z i [ v h ] z i ∂ u∂t e ∂n e ( z i ) ,I , = − (cid:88) e ∈E Bh γ e (cid:90) e ∂ u∂t e ∂n e v h . Suppose that w ∈ H ([ a, b ]) and denote by h ab = b − a , then we have w ( b ) = (cid:90) ba (cid:18) x − ab − a w ( x ) (cid:19) (cid:48) d x = 1 b − a (cid:90) ba w + 2 (cid:90) ba x − ab − a ww (cid:48) (3.14) ≤ h ab (cid:107) w (cid:107) L ([ a,b ]) + 2 | w | H ([ a,b ]) (cid:107) w (cid:107) L ([ a,b ]) , which implies | I , | ≤ (cid:88) z i ∈ Γ (cid:84) N h (cid:12)(cid:12) [ γ e ] z i (cid:12)(cid:12) (cid:18) h e i (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n e i (cid:12)(cid:12)(cid:12)(cid:12) H ( e i ) + 2 (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n e i (cid:12)(cid:12)(cid:12)(cid:12) H ( e i ) (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n e i (cid:12)(cid:12)(cid:12)(cid:12) H ( e i ) (cid:19) / (3.15) · (cid:18) h e i (cid:107) v h (cid:107) L ( e i ∩ e i +1 ) + 2 | v h | H ( e i ) (cid:107) v h (cid:107) L ( e i ∩ e i +1 ) (cid:19) / (cid:46) max z i ∈ Γ (cid:84) N h (cid:12)(cid:12) [ γ e ] z i (cid:12)(cid:12) (cid:18) h (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n (cid:12)(cid:12)(cid:12)(cid:12) H (Γ) + (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n (cid:12)(cid:12)(cid:12)(cid:12) H (Γ) (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n (cid:12)(cid:12)(cid:12)(cid:12) H (Γ) (cid:19) / uperconvergence analysis of linear FEM · (cid:18) h (cid:107) v h (cid:107) L (Γ) + | v h | H (Γ) (cid:107) v h (cid:107) L (Γ) (cid:19) / (cid:46) max z i ∈ Γ (cid:84) N h (cid:12)(cid:12) [ γ e ] z i (cid:12)(cid:12) h (cid:18)(cid:0) | g | H (Γ) + k | u | H (Γ) (cid:1) + h (cid:0) | g | H (Γ) + k | u | H (Γ) (cid:1) · (cid:0) | g | H (Γ) + k | u | H (Γ) (cid:1)(cid:19) / · (cid:18) (cid:107) v h (cid:107) L (Γ) + h | v h | H (Γ) (cid:107) v h (cid:107) L (Γ) (cid:19) / (cid:46) h α k / (cid:18) (cid:107) v h (cid:107) (cid:107) v h (cid:107) H + h / (cid:107) v h (cid:107) / H (cid:107) v h (cid:107) / (cid:19) / C u,g (cid:46) kh α (cid:107)| v h |(cid:107) ,h C u,g , where we have used (3.7).For any z i ∈ N Bh , let τ i, , τ i, , · · · , τ i,n i be n i triangles clockwise around z i (cf.3.1). If n i = 1 for some z i ∈ N Bh , it is easy to see that [ γ e v h ] z i in (3.13) is equalto zero. Thus we assume that n i ≥ z i ∈ N Bh for simplicity of presentation.Denote by e ij = τ i,j ∩ τ i,j +1 for j = 1 , , · · · , n i −
1. Clearly, e ij are in E Ih and we have (cid:12)(cid:12) [ v h ] z i (cid:12)(cid:12) ≤ n i − (cid:88) j =1 (cid:12)(cid:12)(cid:12) [ v h ] e ij (cid:12)(cid:12)(cid:12) ≤ √ n i − (cid:18) n i − (cid:88) j =1 (cid:12)(cid:12)(cid:12) [ v h ] e ij (cid:12)(cid:12)(cid:12) (cid:19) / . It is well known that n i can be bounded by a constant independent of the mesh size h for any regular triangulation. Therefore, we have | I , | ≤ (cid:88) z i ∈N Bh (cid:12)(cid:12) { γ e } z i (cid:12)(cid:12) (cid:18) h e i (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n e i (cid:12)(cid:12)(cid:12)(cid:12) H ( e i ) + 2 (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n e i (cid:12)(cid:12)(cid:12)(cid:12) H ( e i ) (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n e i (cid:12)(cid:12)(cid:12)(cid:12) H ( e i ) (cid:19) / (3.16) · (cid:18) n i − (cid:88) j =1 h e ij (cid:107) [ v h ] (cid:107) L ( e ij ) + 2 | [ v h ] | H ( e ij ) (cid:107) [ v h ] (cid:107) L ( e ij ) (cid:19) / (cid:46) max z i ∈N Bh (cid:12)(cid:12) { γ e } z i (cid:12)(cid:12) (cid:18) h (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n (cid:12)(cid:12)(cid:12)(cid:12) H (Γ) + (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n (cid:12)(cid:12)(cid:12)(cid:12) H (Γ) (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n (cid:12)(cid:12)(cid:12)(cid:12) H (Γ) (cid:19) / · (cid:18) h µ ρ (cid:88) z i ∈N Bh n i − (cid:88) j =1 ρ ,e ij h µe ij (cid:107) [ v h ] (cid:107) L ( e ij ) (cid:19) / (cid:46) h ( µ − / max z i ∈N Bh (cid:12)(cid:12) { γ e } z i (cid:12)(cid:12) (cid:18)(cid:0) | g | H (Γ) + k | u | H (Γ) (cid:1) + h (cid:0) | g | H (Γ) + k | u | H (Γ) (cid:1) · (cid:0) | g | H (Γ) + k | u | H (Γ) (cid:1)(cid:19) / J ( v h , v h ) / (cid:46) k / h (3+ µ ) / (cid:107)| v h |(cid:107) ,h C u,g where we have used (3.6).On the other hand, | I , | (cid:46) h (cid:12)(cid:12)(cid:12)(cid:12) ∂u∂n (cid:12)(cid:12)(cid:12)(cid:12) H (Γ) (cid:107) v h (cid:107) L (Γ) (3.17) (cid:46) h (cid:0) | g | H (Γ) + k | u | H (Γ) (cid:1) · (cid:107) v h (cid:107) / (cid:107) v h (cid:107) / H Y. Du & Z. Zhang (cid:46) k / h (cid:107)| v h |(cid:107) ,h C u,g . From (3.13)–(3.17), (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) e ∈E Bh γ e (cid:90) e ∂ u∂t e ∂n e ∂v h ∂t e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:46) (cid:0) ( kh ) + kh α + k / h (3+ µ ) / (cid:1) (cid:107)| v h |(cid:107) C u,g . (3.18)Then the estimate of I can be obtained from (3.12) and (3.18), | I | (cid:46) (cid:0) ( kh ) + kh α + k / h (3+ µ ) / (cid:1) (cid:107)| v h |(cid:107) ,h C u,g . (3.19)Finally, by combining (3.11) and (3.19), we prove the lemma. Theorem 3.3.
Assume that T h satisfies the mesh condition α . u h is the DGfinite element solution of the scheme (2.1) and u I is the linear interpolation of thesolution u to (1.1) – (1.2) . There exists a constant C independent of k and h , suchthat if k ( kh ) ≤ C , we have (cid:107)| u h − u I |(cid:107) ,h (cid:46) (cid:16) kh α + kh µ/ + k ( kh ) (cid:17) C u,g . (3.20) Proof . For simplicity of presentation, we denote v h = u h − u I . By Lemma 2.1and the Galerkin orthogonality, we have (cid:107)| u h − u I |(cid:107) ,h (cid:46) a h ( u h − u I , v h ) + k ( u h − u I , v h )(3.21)= (cid:60) (cid:0) a h ( u h − u I , v h ) − k ( u h − u I , v h ) + i k (cid:104) u h − u I , v h (cid:105) + 2 k ( u h − u I , v h ) (cid:1) = (cid:60) (cid:0) a h ( u − u I , v h ) − k ( u − u I , v h ) + i k (cid:104) u − u I , v h (cid:105) + 2 k ( u h − u I , v h ) (cid:1) = I + I + I + I . It is well known that | I | ≤ k (cid:107) u − u I (cid:107) k (cid:107) v h (cid:107) (cid:46) kh (cid:107) u (cid:107) (cid:107)| v h |(cid:107) ,h . (3.22)By the trace inequality, we have | I | ≤ k (cid:107) u − u I (cid:107) L ( ∂ Ω) (cid:107) v h (cid:107) L ( ∂ Ω) (3.23) (cid:46) kh (cid:107) u (cid:107) H ( ∂ Ω) (cid:107) v h (cid:107) L ( ∂ Ω) (cid:46) k / h (cid:107) u (cid:107) / (cid:107) u (cid:107) / (cid:107)| v h |(cid:107) ,h (cid:46) ( kh ) (cid:107)| v h |(cid:107) ,h C u,g . From Lemma 2.2, we know that there exists a constant C independent on k and h ,such that if k ( kh ) ≤ C , the following inequality holds, k (cid:107) u − u h (cid:107) L (Ω) (cid:46) (cid:0) ( kh ) + k ( kh ) (cid:1) C u,g , which implies | I | ≤ k (cid:107) u h − u I (cid:107) · k (cid:107) v h (cid:107) (3.24) ≤ k (cid:107) u − u h (cid:107) + k (cid:107) u − u I (cid:107) ) · k (cid:107) v h (cid:107) (cid:46) (cid:0) ( kh ) + k ( kh ) (cid:1) (cid:107)| v h |(cid:107) ,h C u,g . uperconvergence analysis of linear FEM (cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8) z τ z, τ z, τ z, τ z, τ z, τ z, Fig. 3.1 . n z triangles sharing z with n z = 6 . Next, we estimate I . From the definition of a h ( · , · ) and the fact that u I is continuousin Ω, we have | I | ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) τ ∈T h (cid:90) τ ∇ ( u − u I ) · ∇ ¯ v h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) e ∈E Ih (cid:28)(cid:26) ∂ ( u − u I ) ∂ n e (cid:27) , [ v h ] (cid:29) e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = I , + I , . Then I , ≤ (cid:88) e ∈E Ih (cid:13)(cid:13)(cid:13)(cid:13)(cid:26) ∂ ( u − u I ) ∂ n e (cid:27)(cid:13)(cid:13)(cid:13)(cid:13) L ( e ) (cid:107) [ v h ] (cid:107) L ( e ) (3.25) (cid:46) h / (cid:88) e ∈E Ih (cid:107) u (cid:107) H ( τ e ) (cid:107) [ v h ] (cid:107) L ( e ) (cid:46) h µ/ (cid:107) u (cid:107) J ( v h , v h ) / (cid:46) kh µ/ (cid:107)| v h |(cid:107) ,h C u,g . From Lemma 3.2 and (3.25), | I | (cid:46) (cid:16) ( kh ) + kh α + kh µ/ (cid:17) (cid:107)| v h |(cid:107) ,h C u,g . (3.26)Therefore, by combining the equations (3.21)–(3.24) and (3.26), we have if k ( kh ) ≤ C , the following estimate holds: (cid:107)| u h − u I |(cid:107) ,h (cid:46) (cid:16) ( kh ) + kh α + kh µ/ (cid:17) (cid:107)| v h |(cid:107) ,h C u,g . This completes the proof.In this paper, for any node point z ∈ N h , let n z be the number of trianglesassociated with z and let τ z, , τ z, , · · · , τ z,n z be the elements having z such that theyare counterclockwise around z as shown in Figure 3.1. τ z,j means τ z, mod( j,n z ) for theinteger j > n z . Lemma 3.4.
Assume that T h satisfies the mesh condition α and u h is the DGfinite element solution of the scheme (2.1) . For any z ∈ N h and τ ∈ T h having thenode point z , let u h ( z, τ ) be u h | τ ( z ) . There exists a constant C independent of k and h , such that if k ( kh ) ≤ C , we have (cid:32) (cid:88) z ∈N h n z (cid:88) i =1 | u h ( z, τ i +1 ) − u h ( z, τ i ) | (cid:33) / (cid:46) h µ/ (cid:16) kh α + kh µ/ + k ( kh ) (cid:17) C u,g . Y. Du & Z. Zhang
Proof . For any e ∈ E Ih , let z e , z (cid:48) e be its two endpoints and τ e , τ (cid:48) e be two elementssharing e . Then from (3.14) and the discrete inverse inequalities, we have (cid:18) (cid:88) z ∈N h n z (cid:88) i =1 | u h ( z, τ i +1 ) − u h ( z, τ i ) | (cid:19) / = (cid:18) (cid:88) e ∈E Ih | u h ( z e , τ e ) − u h ( z e , τ (cid:48) e ) | + | u h ( z (cid:48) e , τ e ) − u h ( z (cid:48) e , τ (cid:48) e ) | (cid:19) / (cid:46) (cid:18) (cid:88) e ∈E Ih h e (cid:107) [ u h ] (cid:107) L ( e ) + (cid:107) [ u h ] (cid:107) L ( e ) (cid:107) [ u h ] (cid:107) H ( e ) (cid:19) / (cid:46) (cid:18) (cid:88) e ∈E Ih h e (cid:107) [ u h ] (cid:107) L ( e ) (cid:19) / . From the fact that the linear interpolation u I of the solution u to (1.1)–(1.2) is con-tinuous , that is [ u I ] e = 0 ∀ e ∈ E IH , and Theorem 3.3, (cid:18) (cid:88) e ∈E Ih h e (cid:107) [ u h ] (cid:107) L ( e ) (cid:19) / = (cid:18) (cid:88) e ∈E Ih h e (cid:107) [ u h − u I ] (cid:107) L ( e ) (cid:19) / (cid:46) h µ/ J ( u h − u I , u h − u I ) / (cid:46) h µ/ (cid:107)| u h − u I |(cid:107) ,h (cid:46) h µ/ (cid:16) ( kh ) + kh α + kh µ/ (cid:17) C u,g . This completes the proof.
4. The gradient recovery operator G h and its superconvergence. In thissection, we define a polynomial preserving recovery method for the discontinuous finiteelement space and derive the superconvergent error estimate.We first recall a gradient recovery operator developed in 2004 for the continuousfinite element methods, which is called polynomial preserving recovery (PPR). Let ˜ V h be the approximation space of continuous piecewise linear polynomials over T h andlet ˜ G h : C (Ω) (cid:55)→ ˜ V h × ˜ V h be the gradient recovery operator. Given a node z ∈ N h ,we select n ≥ z j ∈ N h , j = 1 , , · · · , n , in an element patch ω z containing z ( z is one of z j ) and fit a polynomial of degree 2, in the least squaressense, with values of w ∈ C (Ω) at those sampling points. First, we find p ∈ P ( ω z )for some w ∈ C (Ω) such that n (cid:88) j =1 ( p − w ) ( z j ) = min q ∈ P n (cid:88) j =1 ( q − w ) ( z j ) . (4.1)Here P ( ω z ) is the well-known piecewise quadratic polynomial space defined on ω z .The recovery gradient at z is then defined as˜ G h w ( z ) = ( ∇ p )( z ) . (4.2)We define a PPR operator G h over the DG finite element space. For any u h ∈ V h ,we define a continuous piecewise linear polynomial ˜ u h ∈ ˜ V h by˜ u h ( z ) = λ u h ( z, τ z, ) + λ u h ( z, τ z, ) + · · · + λ n z u h ( z, τ z,n z ) ∀ z ∈ N h , uperconvergence analysis of linear FEM λ j , j = 1 , , · · · , n z ,are non-negative numbers satisfying λ + λ + · · · + λ n z = 1.We define the gradient recovery operator G h from V h to ˜ V h × ˜ V h by G h u h := ˜ G h ˜ u h . (4.3)Clearly, for any v h ∈ ˜ V h , G h v h is completely equal to ˜ G h v h . Lemma 4.1.
For any node point z ∈ N h , let τ z, , τ z, , · · · , τ z,n z be the n z elements counterclockwise around z as shown in Figure 3.1. The following inequalityholds for any u h ∈ V h (cid:107) G h u h (cid:107) (cid:46) | u h | H ( T h ) + (cid:18) (cid:88) z ∈N h n z − (cid:88) j =1 | u h ( z, τ z,j +1 ) − u h ( z, τ z,j ) | (cid:19) / . (4.4) Proof . We have the property (cid:13)(cid:13)(cid:13) ˜ G h ˜ v h (cid:13)(cid:13)(cid:13) (cid:46) (cid:107)∇ ˜ v h (cid:107) for ˜ v h ∈ ˜ V h (cf. [26]), whichimplies (cid:107) G h u h (cid:107) = (cid:13)(cid:13)(cid:13) ˜ G h ˜ u h (cid:13)(cid:13)(cid:13) (cid:46) (cid:107)∇ ˜ u h (cid:107) (cid:46) | ˜ u h − u h | H ( T h ) + | u h | H ( T h ) . For any τ ∈ T h , let φ τ, , φ τ, and φ τ, be its node bases and let z τ, , z τ, , z τ, ∈ N h be its three vertices satisfying φ τ,i ( z τ,j ) = δ ( i − j ). We have | ˜ u h − u h | H ( T h ) = (cid:18) (cid:88) τ ∈T h (cid:107)∇ (˜ u h − u h ) (cid:107) L ( τ ) (cid:19) / = (cid:18) (cid:88) τ ∈T h (cid:90) τ (cid:12)(cid:12) (cid:88) j =1 (˜ u h ( z τ,j ) − u h ( z τ,j , τ )) ∂ x φ τ,j (cid:12)(cid:12) + (cid:90) τ (cid:12)(cid:12) (cid:88) j =1 (˜ u h ( z τ,j ) − u h ( z τ,j , τ )) ∂ y φ τ,j (cid:12)(cid:12) (cid:19) / ≤ (cid:18) (cid:88) τ ∈T h (cid:0) (cid:88) j =1 | ˜ u h ( z τ,j ) − u h ( z τ,j , τ ) | (cid:1) · (cid:0) (cid:88) j =1 (cid:90) τ | ∂ x φ τ,j | + | ∂ y φ τ,j | (cid:1)(cid:19) / . Since T h is a uniform regular triangulation, it is well known that (cid:80) j =1 (cid:82) τ | ∂ x φ τ,j | + | ∂ y φ τ,j | can be bounded by some constant C independent of the mesh size h and thetriangle τ for any τ ∈ T h . Therefore, we have | ˜ u h − u h | H ( T h ) (cid:46) (cid:18) (cid:88) τ ∈T h (cid:0) (cid:88) j =1 | ˜ u h ( z τ,j ) − u h ( z τ,j , τ ) | (cid:1)(cid:19) / (cid:46) (cid:18) (cid:88) z ∈N h n z (cid:88) j =1 | ˜ u h ( z ) − u h ( z, τ z,j ) | (cid:19) / . Y. Du & Z. Zhang
From the definition of ˜ u h , we have for any z ∈ N h , n z (cid:88) j =1 | ˜ u h ( z ) − u h ( z, τ z,j ) | = n z (cid:88) j =1 | (cid:88) i (cid:54) = j λ i ( u h ( z, τ z,i ) − u h ( z, τ z,j )) | ≤ n z (cid:88) j =1 (cid:0) (cid:88) i (cid:54) = j λ i | u h ( z, τ z,i ) − u h ( z, τ z,j ) | (cid:1) ≤ n z (cid:88) j =1 ( n z − (cid:88) i (cid:54) = j λ i | u h ( z, τ z,i ) − u h ( z, τ z,j ) | ≤ ( n z − n z − (cid:88) j =1 n z (cid:88) i = j +1 ( λ i + λ j ) | u h ( z, τ z,i ) − u h ( z, τ z,j ) | ≤ ( n z − n z − (cid:88) j =1 n z (cid:88) i = j +1 ( λ i + λ j ) (cid:0) i − (cid:88) t = j | u h ( z, τ z,t +1 ) − u h ( z, τ z,t ) | (cid:1) ≤ ( n z − n z − (cid:88) j =1 n z (cid:88) i = j +1 j − (cid:88) t = i ( λ i + λ j ) | u h ( z, τ z,t +1 ) − u h ( z, τ z,t ) | ≤ ( n z − n z − (cid:88) t =1 (cid:0) ( n z − t ) t (cid:88) i =1 λ i + t n z (cid:88) i = t +1 λ i (cid:1) | u h ( z, τ z,t +1 ) − u h ( z, τ z,t ) | ≤ ( n z − n z − (cid:88) t =1 | u h ( z, τ z,t +1 ) − u h ( z, τ z,t ) | . Since n z can be bounded by a constant independent of h and the vertex z , we completethe proof. Lemma 4.2.
For any element τ ∈ T h and any function φ ∈ H (˜ τ ) , (cid:107) G h φ I − ∇ φ (cid:107) L ( τ ) (cid:46) h (cid:107) φ (cid:107) H (˜ τ ) , (4.5) where ˜ τ = (cid:83) { ω z : z ∈ N h ∩ τ } and φ I is the linear interpolant of φ .Proof . The proof is completed by the fact that G h u I = ˜ G h u I and Lemma 4.1 in??. Since u I is continuous, that is˜ u I ( z ) − u I ( z, τ z,j ) = 0 ∀ z ∈ N h and j = 1 , , · · · , n z , we have (cid:107) G h u h − ∇ u (cid:107) ≤ (cid:107) G h ( u h − u I ) (cid:107) + (cid:107) G h u I − ∇ u (cid:107) (4.6) (cid:46) | u h − u I | H ( T h ) + (cid:107) G h u I − ∇ u (cid:107) + (cid:18) (cid:88) z ∈N h n z (cid:88) j =1 | u h ( z ) − u h ( z, τ z,j ) | (cid:19) / . Then by combining Lemmas 3.2–4.2 and the inequality (4.6), we have the followingtheorem which is our main result in the paper.
Theorem 4.3.
Let u and u h be the solutions to (1.1) – (1.2) and the discretesolution, respectively. Assume that T h satisfies the mesh condition α . Then thereexists a constant C independent of k and h such that if k ( kh ) ≤ C , (cid:107) G h u h − ∇ u (cid:107) (cid:46) (cid:16) kh α + kh µ/ + k ( kh ) (cid:17) C u,g . (4.7) uperconvergence analysis of linear FEM
5. The influence of the operator G h to the pollution error. In this section,we estimate the error between G h u h and ∇ u h , which motivate us to combine theRichardson extrapolation and the ppr technique to reduce the numerical errors, anddefine the a posterior estimator in Section 6.First, we define an elliptic projection from V to V h : find u + h ∈ V h such that a h ( u + h , v h ) + i k (cid:10) u + h , v h (cid:11) = a h ( u, v h ) + i k (cid:104) u, v h (cid:105) ∀ v h ∈ V h . (5.1)In other words, the elliptic projection u + h of u is the finite element approximation tothe solution of the following (complex-valued) Poisson problem: − ∆ u = F in Ω , (5.2) ∂u∂n + i ku = g on Γ , (5.3)for some given function F which are determined by u . This kind of elliptic projectionis often used to study some properties, such as stability and convergence, of the FEMfor the Helmholtz problem. Readers are referred to [37, 36, 15, 14]. Lemma 5.1.
Assume that u is H -regular. u + h is its elliptic projection defined by (5.1) . There hold the following estimates: (cid:13)(cid:13) u − u + h (cid:13)(cid:13) ,h (cid:46) inf v h ∈ V h (cid:107)| u − v h |(cid:107) ,h (cid:13)(cid:13) u − u + h (cid:13)(cid:13) (cid:46) h inf v h ∈ V h (cid:107)| u − v h |(cid:107) ,h . Proof . From Lemma 3.5 in [36], we know that (cid:13)(cid:13) u − u + h (cid:13)(cid:13) ,h (cid:46) inf v h ∈ V h (cid:0) (cid:107) u − v h (cid:107) + k (cid:107) u − v h (cid:107) L (Γ) (cid:1) / , (cid:13)(cid:13) u − u + h (cid:13)(cid:13) (cid:46) h inf v h ∈ V h (cid:0) (cid:107) u − v h (cid:107) + k (cid:107) u − v h (cid:107) L (Γ) (cid:1) / . Then the estimates follow from k (cid:107) u − v h (cid:107) L (Γ) (cid:46) k (cid:107) u − v h (cid:107) (cid:107) u − v h (cid:107) H (Ω) (cid:46) k (cid:107) u − v h (cid:107) + (cid:107) u − v h (cid:107) H (Ω) (cid:46) (cid:107)| u − v h |(cid:107) ,h . Lemma 5.2.
Assume that u is the exact solution to (1.1) – (1.2) and u + h is itselliptic projection defined by (5.1) . Assume that T h satisfies the mesh condition α .We have (cid:13)(cid:13)(cid:12)(cid:12) ∇ u + h − ∇ u I (cid:12)(cid:12)(cid:13)(cid:13) ,h (cid:46) (cid:0) kh α + kh µ/ + ( kh ) (cid:1) C u,g . (5.4) Proof . Denote v h = u + h − u I . By the Galerkin orthogonality, (cid:13)(cid:13)(cid:12)(cid:12) ∇ u + h − ∇ u I (cid:12)(cid:12)(cid:13)(cid:13) ,h (cid:46) (cid:60) (cid:0) a h ( u + h − u I , v h ) + i k (cid:10) u + h − u I , v h (cid:11) (cid:1) + k ( u + h − u I , v h ) (cid:46) (cid:60) (cid:0) a h ( u − u I , v h ) + i k (cid:104) u − u I , v h (cid:105) (cid:1) + k ( u + h − u I , v h ) (cid:46) | a h ( u − u I , v h ) | + | k (cid:104) u − u I , v h (cid:105)| + k (cid:13)(cid:13) u + h − u I (cid:13)(cid:13) · k (cid:107) v h (cid:107) . Y. Du & Z. Zhang
By some arguments same to those in 3.3, it is obtained that | a h ( u − u I , v h ) | (cid:46) (cid:16) ( kh ) + kh α + kh µ/ (cid:17) (cid:107)| v h |(cid:107) ,h C u,g and | k (cid:104) u − u I , v h (cid:105)| (cid:46) ( kh ) (cid:107)| v h |(cid:107) ,h C u,g . From Lemma 5.1, we have k (cid:13)(cid:13) u + h − u I (cid:13)(cid:13) · k (cid:107) v h (cid:107) ≤ (cid:0) k (cid:13)(cid:13) u + h − u (cid:13)(cid:13) + k (cid:107) u − u I (cid:107) (cid:1) (cid:107)| v h |(cid:107) ,h (cid:46) ( kh ) (cid:0) √ kh (cid:1) (cid:107)| v h |(cid:107) ,h C u,g . By combining the inequalities above, we complete the proof.
Theorem 5.3.
Let u and u + h be the solution to (1.1) – (1.2) and the elliptic pro-jection defined by (5.1) , respectively. Assume that the mesh condition α is satisfied.Then the following error estimate holds: (cid:13)(cid:13) G h u + h − ∇ u (cid:13)(cid:13) (cid:46) (cid:0) kh α + kh µ/ + ( kh ) (cid:1) C u,g . (5.5) Theorem 5.4.
Let u and u + h be the solution to (1.1) – (1.2) and the elliptic pro-jection defined by (5.1) , respectively. Assume that the mesh condition α is satisfied.We have (cid:18) (cid:88) τ ∈T h (cid:107) G h u h − ∇ u h (cid:107) L ( τ ) (cid:19) / (cid:46) (cid:0) kh + k ( kh ) (cid:1) C u,g . (5.6) Proof . Denote by θ h = u h − u + h . u h can be written as u h = u + h + θ h , where u + h isthe elliptic projection of u defined by (5.1). We have (cid:18) (cid:88) τ ∈T h (cid:107) G h u h − ∇ u h (cid:107) L ( τ ) (cid:19) / (5.7) = (cid:18) (cid:88) τ ∈T h (cid:13)(cid:13) G h ( u + h + θ h ) − ∇ ( u + h + θ h ) (cid:13)(cid:13) L ( τ ) (cid:19) / ≤ (cid:18) (cid:88) τ ∈T h (cid:13)(cid:13) G h u + h − ∇ u + h (cid:13)(cid:13) L ( τ ) (cid:19) / + (cid:18) (cid:88) τ ∈T h (cid:107) G h θ h − ∇ θ h (cid:107) L ( τ ) (cid:19) / . From Lemma 5.1 and Theorem 5.3, (cid:18) (cid:88) τ ∈T h (cid:13)(cid:13) G h u + h − ∇ u + h (cid:13)(cid:13) L ( τ ) (cid:19) / (5.8) ≤ (cid:13)(cid:13) G h u + h − ∇ u (cid:13)(cid:13) + (cid:12)(cid:12) u − u + h (cid:12)(cid:12) H (cid:46) (cid:0) kh α + kh µ/ + ( kh ) (cid:1) C u,g + khC u,g (cid:46) ( kh + ( kh ) ) C u,g . uperconvergence analysis of linear FEM θ h satisfies the following equation: a h ( θ h , v h ) + i k (cid:104) θ h , v h (cid:105) = − k ( u − u h , v h ) . (5.9)Therefor, θ h can be understood as the numerical approximation to the following Pois-son problem with Robin boundary: − ∆ θ = − k ( u − u h ) in Ω ,∂θ∂n + i kθ = 0 on Γ . Therefore, (cid:18) (cid:88) τ ∈T h (cid:107) G h θ h − ∇ θ h (cid:107) L ( τ ) (cid:19) / ≤ (cid:107) G h θ h − ∇ θ (cid:107) + | θ − θ h | H (5.10) (cid:46) h (cid:107) θ (cid:107) (cid:46) k h (cid:107) u − u h (cid:107) (cid:46) k h ( kh + k h ) C u,g (cid:46) (( kh ) + k ( kh ) ) C u,g . The proof is completed by combining (5.7)–(5.10).
6. Numerical Tests.
In this section, some numerical tests are implemented inorder to demonstrate our theoretical results. We simulate the Helmholtz problem(1.1)–(1.2) where the source data f and g is so chosen that the exact solution is u = cos( kr ) r − cos k + i sin kk (cid:0) J ( k ) + i J ( k ) (cid:1) J ( kr )in polar coordinates, where J ( z ) is a Bessel function of the first kind, and Ω =[0 . , . × [0 . , . γ = 5 and we set λ = 1 , λ = 0 , · · · , λ n z = 0in the definition of G h .We first simulate the problem over the regular pattern uniform triangulation anddenote by T N the triangulation consisting of 2 N triangles of size h which is equivalentto 1 /N .From Theorem 3.3, there exists the estimate that if k ( kh ) ≤ C , (cid:107) u h − u I (cid:107) ,h (cid:46) kh µ/ + k ( kh ) . The left graphs in Figure 6.1–Figure 6.3 show the numerical errors (cid:107) u h − u I (cid:107) ,h withpenalty parameters µ = 0 , , k = 5 , ,
50 and 100, respectively. The rightgraphs in Figure 6.1–Figure 6.3 show the convergence orders of the errors (cid:107) u h − u I (cid:107) ,h shown in the left graphs, respectively. As we expected, (cid:107) u h − u I (cid:107) ,h decays at therates of O ( h ), O ( h / ) and O ( h ) for the small wave numbers k = 5 ,
10, respectively.However, we can see that for the large wave number k = 50 , (cid:107) u h − u I (cid:107) ,h does notconverge at first, then begins to decay at the rates which are greater than O ( h µ/ )when N is large enough, which implies the existence of the constraint k ( kh ) ≤ C and the so-called pollution error k ( kh ) .Figure 6.4–Figure ?? show the numerical errors (cid:107) G h u h − ∇ u (cid:107) in left graphs andthe convergence order in right graphs for k = 5 , , ,
100 with penalty parameters µ = 0 , ,
2, respectively. Clearly, the recovered gradients super-converge at the rategreater than O ( h µ/ ). Therefore, whether the estimate (4.7) is sharp with respectto µ is still open. The constraint k ( kh ) ≤ C and the so-called pollution error k ( kh ) can also be observed.8 Y. Du & Z. Zhang N -3 -2 -1 k u h − u I k , h k=5k=10k=50k=100 N -1.5-1-0.500.511.52 c o n v e r g e n ce o r d e r o f k u h − u I k , h k=5k=10k=50k=100 Fig. 6.1 . (cid:107) u h − u I (cid:107) ,h (left) and the convergence order of (cid:107) u h − u I (cid:107) ,h (right) for k =5 , , , , where u h is the numerical solution over the regular pattern uniform triangulation T N with µ = 0 . N -4 -3 -2 -1 k u h − u I k , h k=5k=10k=50k=100 N -1.5-1-0.500.511.52 c o n v e r g e n ce o r d e r o f k u h − u I k , h k=5k=10k=50k=100 Fig. 6.2 . (cid:107) u h − u I (cid:107) ,h (left) and the convergence order of (cid:107) u h − u I (cid:107) ,h (right) for k =5 , , , , where u h is the numerical solution over the regular pattern uniform triangulation T N with µ = 1 / . N -5 -4 -3 -2 -1 k u h − u I k , h k=5k=10k=50k=100 N -1.5-1-0.500.511.52 c o n v e r g e n ce o r d e r o f k u h − u I k , h k=5k=10k=50k=100 Fig. 6.3 . (cid:107) u h − u I (cid:107) ,h (left) and the convergence order of (cid:107) u h − u I (cid:107) ,h (right) for k =5 , , , , where u h is the numerical solution over the regular pattern uniform triangulation T N with µ = 1 . uperconvergence analysis of linear FEM N -5 -4 -3 -2 -1 k G h u h − ∇ u k k=5k=10k=50k=100 N -0.500.511.52 c o n v e r g e n ce o r d e r o f k G h u h − ∇ u k k=5k=10k=50k=100 Fig. 6.4 . (cid:107) G h u h − ∇ u (cid:107) (left) and the convergence order of (cid:107) G h u h − ∇ u (cid:107) (right) for k =5 , , , , where u h is the numerical solution over the regular pattern uniform triangulation T N with µ = 0 . N -6 -5 -4 -3 -2 -1 k G h u h − ∇ u k k=5k=10k=50k=100 N -0.500.511.52 c o n v e r g e n ce o r d e r o f k G h u h − ∇ u k k=5k=10k=50k=100 Fig. 6.5 . (cid:107) G h u h − ∇ u (cid:107) (left) and the convergence order of (cid:107) G h u h − ∇ u (cid:107) (right) for k =5 , , , , where u h is the numerical solution over the regular pattern uniform triangulation T N with µ = 1 / . m k=10 k=50 E E E E E E Table 6.1
The numerical errors E := | u − u h | H ( T h ) , E := (cid:107)∇ u − G h u h (cid:107) and E := (cid:107)∇ u − RG H u h (cid:107) with µ = 0 over T m ( m = 4 , , , . . . , ) for k = 10 , . Y. Du & Z. Zhang m k=10 k=50 E E E E E E Table 6.2
The numerical errors E := | u − u h | H ( T h ) , E := (cid:107)∇ u − G h u h (cid:107) and E := (cid:107)∇ u − RG H u h (cid:107) with µ = 1 over T m ( m = 4 , , , . . . , ) for k = 10 , . m k=10 k=50 E E E E E E Table 6.3
The numerical errors E := | u − u h | H ( T h ) , E := (cid:107)∇ u − G h u h (cid:107) and E := (cid:107)∇ u − RG H u h (cid:107) with µ = 2 over T m ( m = 4 , , , . . . , ) for k = 10 , . m k=10 k=60 k=120 E η h E η h E η h Table 6.4
The numerical errors E := | u − u h | H ( T h ) and η h with µ = 0 over T m ( m = 4 , , , . . . , )for k = 10 , , . uperconvergence analysis of linear FEM
21m k=10 k=60 k=120 E η h E η h E η h Table 6.5
The numerical errors E := | u − u h | H ( T h ) and η h with µ = 1 over T m ( m = 4 , , , . . . , )for k = 10 , , . REFERENCES[1]
C. AB. , COMSOL MultiPhysics User’s Guide , 3.5a ed., 2008.[2]
M. Ainsworth , Discrete dispersion relation for hp-version finite element approximation athigh wave number , SIAM J. Numer. Anal., 42 (2004), pp. 553–575.[3]
A. Aziz and R. Kellogg , A scattering problem for the Helmholtz equation , in Advances inComputer Methods for Partial Differential Equations-III, vol. 1, 1979, pp. 93–95.[4]
I. Babuˇska, F. Ihlenburg, E. Paik, and S. Sauter , A generalized finite element methodfor solving the Helmholtz equation in two dimensions with minimal pollution , Comput.Methods Appl. Mech. Engrg., 128 (1995), pp. 325–359.[5]
I. Babuˇska and S. Sauter , Is the pollution effect of the FEM avoidable for the Helmholtzequation considering high wave numbers? , SIAM Rev., 42 (2000), pp. 451–484.[6]
R. E. Bank and J. Xu , Asymptotically exact a posteriori error estimators, Part I: Grid withsuperconvergence , SIAM J. Numer. Anal., 41 (2003), pp. 2294–2312.[7]
S. Brenner and L. Scott , The mathematical theory of finite element methods , Springer, NewYork, third ed., 2008.[8]
E. Burman, H. Wu, and L. Zhu , Continuous interior penalty finite element method forHelmholtz equation with high wave number: One dimensional analysis , arXiv:1211.1424.[9]
L. Chen and J. Xu , Topics on adaptive finite element methods, in Adaptive Computations:Theory and Algorithms, T. Tang and J. Xu, eds. , Science Press, Beijing, 2007.[10]
Z. Chen and X. Xiang , A source transfer domain decomposition method for helmholtz equa-tions in unbounded domain , SIAM J. Numer. Anal., 51 (2013), pp. 2331–2356.[11]
P. G. Ciarlet , The finite element method for elliptic problems , North-Holland Pub. Co., NewYork, 1978.[12]
A. Deraemaeker, I. Babuˇska, and P. Bouillard , Dispersion and pollution of the FEMsolution for the Helmholtz equation in one, two and three dimensions , Internat. J. Numer.Methods Engrg., 46 (1999), pp. 471–499.[13]
J. Douglas Jr, J. Santos, and D. Sheen , Approximation of scalar waves in the space-frequency domain , Math. Models Methods Appl. Sci., 4 (1994), pp. 509–531.[14]
Y. Du and H. Wu , Preasymptotic error analysis of higher order fem and cip-fem for Helmholtzequation with high wave number , SIAM J. Numer. Anal., 53 (2015), pp. 782–804.[15]
Y. Du and L. Zhu , Preasymptotic error analysis of high order interior penalty discontinuousGalerkin methods for the Helmholtz equation with high wave number , J. Sci. Comput.,Accepted, (2015).[16]
B. Engquist and A. Majda , Radiation boundary conditions for acoustic and elastic wavecalculations , Comm. Pure Appl. Math., 32 (1979), pp. 313–357.[17]
X. Feng and H. Wu , Discontinuous Galerkin methods for the Helmholtz equation with largewave numbers , SIAM J. Numer. Anal., 47 (2009), pp. 2872–2896.[18] , hp -discontinuous Galerkin methods for the Helmholtz equation with large wave number ,Math. Comp., 80 (2011), pp. 1997–2024. Y. Du & Z. Zhang[19]
I. Harari , Reducing spurious dispersion, anisotropy and reflection in finite element analysisof time-harmonic acoustics , Comput. Meth. Appl. Mech. Engrg., 140 (1997), pp. 39–58.[20]
F. Ihlenburg and I. Babuˇska , Finite element solution of the Helmholtz equation with highwave number. I. The h -version of the FEM , Comput. Math. Appl., 30 (1995), pp. 9–37.[21] , Finite element solution of the Helmholtz equation with high wave number. II. The h - p version of the FEM , SIAM J. Numer. Anal., 34 (1997), pp. 315–358.[22] A. M. Lakhany, I. Marek, and J. R. Whiteman , Superconvergence results on mildly struc-tured triangulations , Comput. Methods Appl. Mech. Engrg., 189 (2000), pp. 1–75.[23]
J. Melenk, A. Parsania, and S. Sauter , General DG-methods for highly indefinite Helmholtzproblems , Journal of Scientific Computing, 57 (2013), pp. 536–581.[24]
J. M. Melenk and S. Sauter , Convergence analysis for finite element discretizations ofthe Helmholtz equation with Dirichlet-to-Neumann boundary conditions , Math. Comp.,79 (2010), pp. 1871–1914.[25] ,
Wavenumber explicit convergence analysis for Galerkin discretizations of the Helmholtzequation , SIAM J. Numer. Anal., 49 (2011), pp. 1210–1243.[26]
A. Naga and Z. Zhang , A posteriori error estimates based on the polynomial preservingrecovery , SIAM J. Numer. Anal., 42 (2004), pp. 1780–1800.[27]
B. Rivi`ere , Discontinous Galerkin Methods for Solving Elliptic and Parabolic Equations: the-ory and implementation , Philadelphia, PA : SIAM, Society for Industrial and AppliedMathematics, 2008.[28]
A. Schatz , An observation concerning Ritz–Galerkin methods with indefinite bilinear forms ,Math. Comp., 28 (1974), pp. 959–962.[29]
H. Wu and Z. Zhang , Can we have superconvergent gradient recovery under adaptive meshes? ,SIAM J. Numer. Anal., 45 (2007), pp. 1701–1722.[30]
J. Xu and Z. Zhang , Analysis of recovery type a posteriori error estimators for mildly struc-tured grids , Math. Comp., 73 (2003), pp. 1139–1152.[31]
N. Yan and A. Zhou , Gradient recovery type a posteriori error estimates for finite elementapproximations on irregular meshes , Comput. Methods Appl. Mech. Engrg., 190 (2001),pp. 4289–4299.[32]
Z. Zhang , Polynomial preserving gradient recovery and a posteriori estimate for bilinear ele-ment on irregular quadrilaterals , Internat. J. Numer. Anal. Model., 1 (2004), pp. 1–24.[33] ,
Polynomial preserving recovery for anisotropic and irregular grids , J. Comput. Math.,22 (2004), pp. 331–340.[34]
Z. Zhang and B. Li , Analysis of a class of superconvergence patch recovery techniques forlinear and bilinear finite elements , Numer. Methods Partial Differential Equations, 15(1999), pp. 151–167.[35]
Z. Zhang and A. Naga , A new finite element gradient recovery method: Superconvergenceproperty , SIAM J. Sci. Comput., 26 (2005), pp. 1192–1213.[36]
L. Zhu and Y. Du , Pre-asymptotic error analysis of hp -interior penalty discontinuous Galerkinmethods for the Helmholtz equation with large wave number , Comput. Math. Appl., 70(2015), pp. 917–933.[37] L. Zhu and H. Wu , Pre-asymptotic error analysis of CIP-FEM and FEM for Helmholtzequation with high wave number. Part II: hp versionversion