WEAK GALERKIN FINITE ELEMENT METHODS COMBINED WITH CRANK-NICOLSON SCHEME FOR PARABOLIC INTERFACE PROBLEMS

This article is devoted to the a priori error estimates of the fully discrete Crank-Nicolson approximation for the linear parabolic interface problem via weak Galerkin finite element methods (WG-FEM). All the finite element functions are discontinuous for which the usual gradient operator is implemented as distributions in properly defined spaces. Optimal order error estimates in both L∞(H1) and L∞(L2) norms are established for lowest order WG finite element space (Pk(K), Pk−1(∂K), [ Pk−1(K) ]2 ). Finally, we give numerical examples to verify the theoretical results.


Introduction
Let Ω be a convex polygonal domain in R 2 with boundary ∂Ω and Here η is the outward pointing unit normal to Ω 1 and [v] denotes the jump of a quantity v across the interface Γ i.e., [v] The coefficient function β(x) is assumed to be positive and piecewise constant across Γ, i.e., β(x) = β k for x ∈ Ω k , k = 1, 2. Across the interface Γ, the source function f : Ω × (0, T ] → R can be singular. We assume that f is sufficiently smooth locally. Jump functions ψ : Γ × (0, T ] → R and ϕ : Γ × (0, T ] → R are given, and T < ∞. Medium heterogeneities makes the construction of stable and accurate numerical schemes for differential equations more challenging. The past few decades have witnessed intensive research activity in interface problems via finite element algorithms. For recent literature, one may refer to [2,[5][6][7][8] and the references therein. The objective of the present work is to propose and analyze weak Galerkin finite element method (WG-FEM ) for parabolic interface problems. The WG-FEM introduced in [9] refers to the numerical techniques for partial differential equations where the differential operators are approximated by weak forms. In [10], a WG-FEM was developed for the second order elliptic equation in mixed form. The resulting WG mixed finite element schemes turned out to be applicable for general finite element partitions consisting of shape regular polytopes, and the stabilization idea opened a new door for weak Galerkin method. Recently, WG-FEM have been applied to interface problems [4,5,8]. The WG algorithm in [8] allows the use of finite element partitions consisting of general polytopal meshes and assume that grid line exactly follows the actual interface. Optimal order error estimate in H 1 norm is established for WG finite element space (P k (K), P k (∂K), P k−1 (K) 2 ). Here, K is any polygonal domain with boundary ∂K and P k (K) (P k (∂K)) is a set of polynomials on K (∂K) with degree no more than k ≥ 1. Then the work of [8] has been extended to elliptic and parabolic interface problems in [4,5] for lowest order WG finite element . The time discretization in [5] is based on the backward Euler approximation. In this paper, we study WG-FEM with Crank-Nicolson scheme for solving parabolic interface problems. Optimal order error estimates in L ∞ (L 2 ) norm is established for the fully discrete scheme. WG-FEM with second-order accuracy in time for parabolic problems without interface can be found in [11].

Preliminaries and weak Galerkin discretization
Throughout the work, we will follow the usual notation for Sobolev spaces and norms (cf.  Let T h be a partition of the domain Ω consisting of polygons in two dimension satisfying a set of conditions specified in [8,10]. Denote by E h the set of all edges in T h and let E 0 h = E h \∂Ω be the set of all interior edges. Let Γ h be the subset of E h of all edges on Γ. For every element K ∈ T h , we denote by h K its diameter and (2.1) Note that T 1 contains all elements in Ω 1 and non-interface elements in Ω 2 .
Let K be any polygonal domain with interior K 0 and boundary ∂K. A weak function on the region K refers to a pair of scalar-valued functions where n is the outward normal to ∂K. The usual L 2 -inner product can be written locally on each element as follows For each element K ∈ T h , denote by Q 0 the usual L 2 projection operator from L 2 (K) onto P k (K) and by Q b the L 2 projection from L 2 (e) onto P k−1 (e) for any e ∈ E h . We shall combine Q 0 with Q b by writing Q h = {Q 0 , Q b }. We recall following crucial approximation properties for local projections Q h . Lemma 2.1 (Lemma 3.4, [8]). Let T h be a finite element partition of Ω satisfying the shape regularity assumption as specified in [10].
Let e ⊂ Γ be shared by two elements K 1 ⊂ Ω 1 and K 2 ⊂ Ω 2 . We define four forms Here, ⟨·, ·⟩ e denotes the L 2 inner product on e ∈ E h . The continuous-time weak Galerkin finite element approximation to (1.1)-(1.3) can be obtained by seeking The bilinear map a(·, ·) on V 0 h is given by where the stabilizer s(·, ·) is defined as The finite element space V 0 h is a normed linear space with respect to a triple-bar norm given by (cf. [5]) We, now, turn our attention to discrete time weak Galerkin procedures. We first divide the interval [0, T ] into M equally-spaced subintervals by the following points 0 = t 0 < t 1 < · · · < t M = T, with t n = nτ , τ = T /M be the time step. For a smooth function ξ on [0, T ], define ξ n = ξ(t n ) and∂ h be the fully discrete approximation of u at t = t n which we shall define through the following scheme:

Error analysis for the fully discrete scheme
This section deals with the error analysis for the fully discrete scheme. We derive optimal order error bounds in H 1 -norm and L 2 -norm. The basic idea applied is to use elliptic projection. Let X * be the collection of all v ∈ L 2 (Ω) with the property Further, in view of (3.1), this definition may be expressed by saying that R h v is the weak Galerkin finite element solution of the elliptic interface problem with exact solution v ∈ H 1 0 (Ω) (cf. [4,8])

2)
[ The error estimate for R h , as shows in the following lemma, should be applied.

Lemma 3.1 ( [4, 8]). Let R h be defined by (3.1). Assume that the exact solution of problem (3.2) is so regular that
For fully discrete error estimates, we now split the errors at t = t n as follows We denote our error as Using ρ and θ, error e n can be further separated as where w n =∂R h u n −û n t . For simplicity of the exposition, we write w n = w n 1 + w n 2 , where w n 1 =∂R h u n −∂u n and w n 2 =∂u n −û n t . Now, setting v =θ n in (3.5), we have (∂θ n ,θ n ) + a(θ n ,θ n ) = −(w n ,θ n ).
Since a(θ n ,θ n ) ≥ 0, we have The term w j 1 can be expressed as Summing over j from j = 1 to j = n in (3.10), we obtain Combining (3.9), (3.11) and (3.7), and further using the fact that An application of Lemma 3.1 for ρ n yields Again, it is easy to verify that Thus, we have Combining estimates (3.12) and (3.13) along with Lemma 2.1, we obtain following optimal L ∞ (L 2 ) norm error estimate (∥u∥ H 1 (H k+1 (Ωi)) + ∥u ttt ∥ L 2 (L 2 (Ωi)) ) .

Numerical Example
We present in this section numerical results to validate the theoretical estimates presented in Section 3. For our numerical experiment, we use lowest order weak Galerkin space (P 1 (K), P 0 (∂K), P 0 (K) 2 ) based on uniform triangulations of Ω i , i = 1, 2. The nodes of the triangulations of Ω 1 and Ω 2 coincide on the interface Γ. Note that for each iteration, the spatial mesh size becomes half of the previous mesh size. We choose the uniform time step τ = 1 10 h. Example 4.1. We consider the two dimensional domain Ω = (−1, 1) × (−1, 1) and the interface is taken to be the circle x 2 + y 2 = 1 4 . We select the appearing in (1.1)-(1.3) setting exact solution as with β 1 = 10 −4 and β 2 = 1.
For the L ∞ (0, T ; L 2 (Ω)) error with τ = O(h), we observe its experimental order of convergence (EOC). For each run i, EOC of a given sequence of L ∞ (0, T ; L 2 (Ω)) errors e(i) defined on a sequence of meshes of size h(i) by EOC(e(i)) = log e(i + 1)/e(i) The convergence behavior of the fully discrete weak Galerkin solutions at final time T = 1 with respect to the L 2 -norm and H 1 -norm are also depicted in Tables 1-2. It is clear from these tables that we have achieved optimal order of convergence in both the norms, which confirm the theoretical prediction as proved in Theorems 3.1-3.2.  Example 4.2. In our second numerical example, we consider the square domain Ω = (−1, 1)×(−1, 1) and the interface is taken to be the ellipse {(x, y) : 4x 2 +16y 2 = r 2 = 1}. We select the data in (3.14) such that the exact solution u is given by  Tables 3-4 represent the numerical solution errors and convergence rates in L 2 and H 1 norms, respectively. In both cases, errors are calculated at time t = 1 and clearly demonstrates the second order of convergence in L 2 norm and first order of convergence in H 1 norm.