Two – Grid hp – DGFEM for Second Order Quasilinear Elliptic PDEs Based on an Incomplete Newton Iteration

In this paper we propose a class of so-called two-grid hp-version discontinuous Galerkin finite element methods for the numerical solution of a second-order quasilinear elliptic boundary value problem based on the application of a single step of a nonlinear Newton solver. We present both the a priori and a posteriori error analysis of this two-grid hp–version DGFEM as well as performing numerical experiments to validate the bounds.


Introduction
In our recent articles [4,5] we have considered a class of two-grid finite element methods for strongly monotone partial differential equations.Here, the underlying problem is first approximated on a coarse finite element space; the resulting coarse solution is then used to linearise the underlying problem on a finer finite element space, so that only a linear system of equations is solved on this richer space.In this paper we consider an alternative two-grid interior penalty (IP) discontinuous Galerkin finite element method (DGFEM), based on employing a single step of a Newton solver on the finer space, cf.[1], [9, Section 5.2], for the numerical solution of the following quasilinear elliptic boundary value problem: where Ω is a bounded polygonal domain in R 2 , with boundary Γ and f ∈ L 2 (Ω).
We assume that µ ∈ C 2 ( Ω × [0, ∞)) satisfies the condition: there exists positive constants m µ and M µ such that the following monotonicity property is satisfied: For ease of notation we write µ(t) instead of µ(x, t).The outline of this article is as follows.In Section 2 we state the proposed two-grid IP DGFEM.In Sections 3 and 4 we consider the a priori and a posteriori error analysis, respectively, of the two-grid IP DGFEM.Finally, in Section 5 we present some numerical results to validate the theoretical error bounds.

Two-Grid hp-Version IP DGFEM
We consider shape-regular meshes T h that partition Ω ⊂ R 2 into open disjoint elements κ such that Ω = κ∈T h κ.By h κ we denote the element diameter of κ ∈ T h , h = max κ∈T h h κ , and n κ signifies the unit outward normal vector to κ.We allow the meshes T h to be 1-irregular ; further, we suppose that T h is of bounded local variation, i.e., there exists a constant ρ 1 ≥ 1, independent of the element sizes, such that ρ −1 1 ≤ hκ /h κ ≤ ρ 1 , for any pair of elements κ, κ ∈ T h which share a common edge e = ∂κ ∩ ∂κ .To each κ ∈ T h we assign a polynomial degree p κ ≥ 1 and define the degree vector p = {p κ : κ ∈ T h }.We suppose that p is also of bounded local variation, i.e., there exists a constant ρ 2 ≥ 1, independent of the element sizes and p, such that, for any pair of neighbouring elements κ, κ ∈ T h , ρ −1 2 ≤ pκ /p κ ≤ ρ 2 .With this notation, we introduce the finite element space Here, for p ≥ 0, P p (κ) denotes the space of polynomials of degree at most p on κ, while Q p (κ) is the space of polynomials of degree at most p in each variable on κ.
For the mesh T h , we write E I h to denote the set of all interior edges of the partition T h of Ω, E B h the set of all boundary edges of T h , and set Let v and q be scalar-and vector-valued functions, respectively, which are sufficiently smooth inside each element κ ∈ T h .Given two adjacent elements, κ + , κ − ∈ T h which share a common edge e ∈ E I h , i.e., e = ∂κ + ∩ ∂κ − , we write v ± and q ± to denote the traces of the functions v and q, respectively, on the edge e, taken from the interior of κ ± , respectively.With this notation, the averages of v and q at x ∈ e are given by { {v} } = 1 /2(v + + v − ) and { {q} } = 1 /2(q + + q − ), respectively.Similarly, the jumps of v and q at x ∈ e are given by [[v] respectively, where n κ ± denotes the unit outward normal vector on ∂κ ± , respectively.On a boundary edge e ∈ E B h , we set { {v} } = v, { {q} } = q, [[v]] = vn and [[q]] = q • n, with n denoting the unit outward normal vector on the boundary Γ.For e ∈ E h , we define h e to be the length of the edge; moreover, we set p e = max(p κ , p κ ), if e = ∂κ ∩ ∂κ ∈ E I h , and p e = p κ , if e = ∂κ ∩ Γ ∈ E B h .2.1.Standard IP DGFEM discretisation.Given a fine mesh partition T h of Ω, with the corresponding polynomial degree vector p, the standard IP DGFEM is defined as follows: find Here, θ ∈ [−1, 1], ∇ h is the element-wise gradient operator and σ h,p = γp 2 e /h e , where γ > 0 is a sufficiently large constant.We define the energy norm on V (T h , p): Lemma 2.1 (See [6]).The semilinear form A h,p (•, •) is strongly monotone in the sense that, there exists γ min > 0, such that for any γ ≥ γ min (2.2) , where C m is a positive constant, independent of the discretisation parameters.
2.2.Two-grid IP DGFEM discretisation.We now introduce a two-grid IP DGFEM based on employing a single step of the Newton iteration on the fine mesh.To this end, we consider two partitions T h and T H of Ω, with granularity h and H, respectively.We assume that T h and T H are nested in that sense that for any element κ h ∈ T h there exists an element κ H ∈ T H such that κ h ⊆ κ H . Moreover for each mesh, T h and T H , we have a corresponding polynomial degree vector p = {p κ : κ ∈ T h } and P = {p κ : κ ∈ T H }, respectively, where given an element κ h ∈ T h and an element κ H ∈ T H , such that κ h ⊆ κ H , the polynomial degree vectors satisfy the condition that p κ h ≥ p κ H . Thereby, the finite element spaces V (T h , p) and V (T H , P ) satisfy the following the condition: Using this notation we introduce the hp-version two-grid IP DGFEM discretisation of (1.1) based on a single Newton iteration step, cf.[1], [9, Section 5.2]: (1) Compute the coarse grid approximation u H,P ∈ V (T H , P ) such that (2) Determine the fine grid solution for all v h,p ∈ V (T h , p).Here, A h,p [u](φ, v) denotes the Fréchet derivative of u → A h,p (u, v), for fixed v, evaluated at u; thereby, given φ we have Remark 2.2.For simplicity of presentation, throughout the rest of this article we shall only consider the incomplete IP variation of the DGFEM, i.e., when θ = 0. Lemma 2.3.Under the assumptions on µ, the following inequality holds: Proof.Setting w 1 = u + tv and Taking the limit as t → 0, we deduce the statement of the Lemma.

A Priori Error Analysis
For simplicity of presentation, in this section we assume that the mesh is quasiuniform with mesh size h and that p is uniform over the mesh, i.e., p ≡ p.
3.1.Auxiliary Results.We first state the following auxiliary results.Lemma 3.2.For a function v ∈ V (T h , p) we have the inverse inequality where C is a positive constant, independent of the discretisation parameters.
Proof.Given κ ∈ T h , employing standard inverse inequalities, see [8], gives Summing over κ ∈ T h , employing the inequality √ a i 2 , a i ≥ 0, i = 1, . . ., n, and taking the fourth root of both sides, completes the proof.
where the remainder Q satisfies DG ∇φ DG , and C is a positive constant, independent of the discretisation parameters.
Proof.We follow the proof outlined by [9, Lemma 3.1]; to this end, setting ξ(t) = v + t(w − v) and η(t) = A h,p (ξ(t), φ), we note that the first equation follows from the identity Here . Secondly, term T 2 is bounded in an analogous fashion as follows: Term T 3 is bounded via the inverse trace inequality, see [8], and Lemma 3.2: We can bound T 4 in an analogous manner as follows: Combining these bounds for terms T 1 , T 2 , T 3 and T 4 completes the proof.
2 , and u h,p ∈ V (T h , p) be the IP DGFEM defined by (2.1), we have that where C is a positive constant, independent of the discretisation parameters.
Proof.Writing P u to denote the projection of u onto the finite element space V (T h , p) defined in [2], we have that u − P u H q (Ω) ≤ C h 2−q p 2−q u H 2 (Ω) and ∇(u − P u ) L ∞ (Ω) ≤ C u H 2 (Ω) for all q ≤ 2. Exploiting these bounds, standard inverse inequalities, [8], and the a priori bound for the IP DGFEM, [6], gives 2 , the quantities u H 2 (Ω) and ∇u L ∞ (Ω) are both bounded uniformly by a constant; this then completes the proof.
3.2.Proof of Theorem 3.1.We now exploit the above results to prove Theorem 3.1.For the first bound (3.1), we employ Lemma 2.3, (2.1), (2.4) and (3.3); thereby, with φ = u h,p − u 2G , we deduce that Hence, from Lemma 3.3 we get that Applying Lemma 3.4, noting that p 3/2 ≥ P 3/2 ≥ 1, and the a priori bound for the standard IP DGFEM, cf.[6, Theorem 3.3], gives Noting that h ≤ H and p ≥ P completes the proof of the first bound (3.1).To prove the second inequality (3.2), we first employ the triangle inequality Thereby, applying the a priori error bound for the standard IP DGFEM, together with the bound (3.1), completes the proof of Theorem 3.1.

A Posteriori Error Analysis
Here, we state an a posteriori error bound for the two-grid IP DGFEM.
Theorem 4.1.Let u ∈ H 1 0 (Ω) be the analytical solution of (1.1), u H,P ∈ V (T H , P ) and u 2G ∈ V (T h , p) the numerical approximations obtained from (2.3) and (2.4), respectively; then the following hp-a posteriori error bound holds with a constant C > 0, which is independent of h, H, p and P .Here, for κ ∈ T h , and Π κ,pκ denotes the (elementwise) L 2 -projection onto V (T h , p).
Proof.The proof of this error bound follows in an analogous manner to the a posteriori proof presented in [5], cf. also [7].For details, we refer to [3].

Numerical Experiments
In this section we perform numerical experiments to validate the a priori error bound, Theorem 3.1 and demonstrate the performance of the a posteriori error bound, Theorem 4.1; here, we set γ = 10 and θ = 0. Throughout this section, we let Ω be the unit square (0, 1) 2 ⊂ R 2 and define the nonlinear coefficient as µ(x, |∇u|) = 2+(1+|∇u|) −1 .We select the right-hand forcing function f so that the analytical solution to (1.1) is given by u(x, y) = x(1 − x)y(1 − y)(1 − 2y)e −20(2x−1) 2 .We first validate the bound given in Theorem 3.1; to this end we first solve the standard IP DGFEM on a 256 × 256 uniform mesh of quadrilaterals to compute u h,p for a fixed constant polynomial degree p = 1, 2, 3. We then compute the solution u 2G to (2.3)-(2.4),for p = 1, 2, 3, on a fixed fine 256 × 256 mesh, while performing uniform h-refinement of the coarse mesh, starting from a 4 × 4 mesh with polynomial degree P = p. Figure 1 shows the convergence rate of the error between u h,p and u 2G , measured in the DG norm, compared to the size of the coarse mesh.Here, we observe that u h,p − u 2G DG tends to zero at the optimal rate O(H 2P ), for each fixed P , cf.Theorem 3.1.

Adaptive
Refinement using Theorem 4.1.For this experiment we use the two-grid mesh adaptation algorithm from [5], with the local error indicators η κ and local two-grid error indicators ξ κ from Theorem 4.1, to automatically refine the coarse and fine meshes employing both h-and hp-adaptive mesh refinement.Figure 2 shows u − u 2G DG compared to the third root of the degrees of freedom, as well as the effectivity indices of the error estimator.As can be seen for both hand hp-adaptive refinement, the effectivity indices are roughly constant, indicating that the error bound overestimates the error by a roughly constant factor.For reference purposes, we also calculate the standard IP DGFEM solution u h,p , using , µ ∇u (|•|) and µ ∇u (|•|) denote the first and second derivatives of µ(|•|), respectively.First consider T 1 : given that µ ∈ C 2 ( Ω × [0, ∞)), Lemma 3.2 gives

Figure 1 .Figure 2 .
Figure 1.Convergence of error between u 2G and u h,p .