Singularly perturbed reaction-diffusion problems as first order systems

We consider a singularly perturbed reaction diffusion problem as a first order two-by-two system. Using piecewise discontinuous polynomials for the first component and $H_{div}$-conforming elements for the second component we provide a convergence analysis on layer adapted meshes and an optimal convergence order in a balanced norm that is comparable with a balanced $H^2$-norm for the second order formulation.


Introduction
Consider the singularly perturbed reaction diffusion problem, given in Ω = (0, 1) 2 by where 0 < ε ≪ 1, c ∈ W 1,∞ , c ∞ ≥ c ≥ c 0 > 0 and u = 0 on ∂Ω. We rewrite the problem, using u = −ε grad • u, into a first order system where grad • denotes the gradient in H 1 0 (Ω) and div its adjoint, the divergence. This formulation is also called a mixed formulation. For its weak formulation let ·, · denote the L 2 -scalar product over Ω. Then (2) becomes with V = (v, v) ∈ L 2 (Ω) × H div (Ω) cu, v + ε div u, v + u, v + ε grad • u, v = f, v * Institute of Scientific Computing, Technische Universität Dresden, Germany. e-mail: sebastian.franz@tu-dresden.de which can also be written for U = (u, u) ∈ L 2 (Ω) × H div (Ω) as ( This is the weak form we will discretise and analyse. Singularly perturbed reaction diffusion problems were analysed in many papers, see e.g. [2,18] The associated norm to (1) is the ε-weighted H 1 -norm, also called energy norm. Unfortunately, that norm is not strong enough to see the boundary layers. For the boundary x = 0 the corresponding layer function is of the type e −x/ε . Here it holds e −x/ε L 2 (Ω) + ε grad e −x/ε L 2 (Ω) Therefore over the last years convergence in a balanced norm, where the boundary layers do not vanish for ε → 0, was considered, see [1,9,11,17]. For the lowest order Raviart-Thomas elements on a Shishkin mesh the system (3) was also considered in [15,Section 5] and analysed in a balanced H 1 -comparable norm.
In this paper we prove optimal convergence orders in a stronger balanced H 2 -comparable norm for a variety of H div -conforming elements on general layer-adapted meshes. The paper is organised as follows. In Section 2 we define the numerical method and recall results for a solution decomposition and interpolation errors. In Section 3 we provide the convergence analysis and in the final Section 4 some numerical examples illustrating our theoretical results are given. Notation: We denote vector valued functions with a bold font. L p (D) with the norm · L p (D) is the classical Lebesque space of function integrable to the power p over a domain D ⊂ R 2 and W ℓ,p (D) the corresponding Sobolev space for derivatives up to order ℓ. Furthermore, we write A B if there exists a generic constant C > 0 such that A ≤ C ·B.

Numerical method and interpolation errors
In order to define our numerical method, we need discrete spaces defined over an appropriate mesh. A basic tool for defining this mesh is the knowledge of a solution decomposition, especially the structure of layers.
Assumption 2.1. The solution u of (1) can be written as where s is the smooth part, w i are boundary layers and w ij are corner layers (both counted counterclockwise). To be more precise, for any given degree k it holds for 0 ≤ i, j ≤ k + 2, and analogously for the other boundary layers and corner layers.
Remark 2.2. Using above solution decomposition for u we derive a similar decomposition for the solution U of (3), as U = (u, u) and u = −ε grad u. Such assumptions on a solution decomposition are very common in the analysis of singularly perturbed problems. They hold true under compatibility and regularity conditions on the data, see e.g. [6,10,12].
We follow [16] and construct an S-type mesh using the information of the solution decomposition. First, we define a transition point λ, such that a typical boundary layer function is small enough: for a constant σ > 0 specified later. We additionally assume as otherwise ε is large enough to facilitate a standard numerical analysis. Now 0 = x 0 < x 1 < · · · < x N = 1 are given by where φ is a mesh-generating function with the properties • φ is monotonically increasing, • φ(0) = 0 and φ(1/2) = ln N, The first three conditions are given in [16], while the last one allows the mesh-widths inside the boundary layers to be bounded from below, see also [8]. Related to φ we define the mesh characterising function ψ by Several S-type meshes are given in [16] fulfilling above properties. We only provide the definitions of the two mostly used. For the Shishkin mesh we have φ(t) = 2t ln N, ψ(t) = N −2t , max |ψ ′ | = 2 ln N and the Bakhvalov-S-mesh In addition to the mesh generating and characterising functions also max |ψ ′ | := max is given, that enters all the error estimates on S-type meshes.
The two-dimensional mesh T N is then defined by all cells and also the simpler bound Let us denote two subdomains of Ω per layer function, exemplarily given for w 1 by where Q k (K) is the space of polynomials with degree up to k in each variable on the cell K of T N . For the discretisation of H div with D k (K) we can use • the Raviart-Thomas space introduced by Raviart and Thomas in [14] on triangular meshes, see also [4] for rectangular meshes, where Q p,q (K) is the space of polynomials with degree p in x and degree q in y on the cell K or see [5], where P k (K) is the space of polynomials of total degree k on the cell K and curl w = (∂ y w, −∂ x w).
Then the discrete method reads: Note that the solution U of (3) does also fulfill (5) and we therefore have Galerkin orthogonality

Numerical analysis
Let us start with an interpolation operator into U N given by its two components. The first one I 1 will be a weighted local L 2 -projection, defined on any K ⊂ T N by Note that we have similar to the standard L 2 -projection anisotropic interpolation error estimates, i.e. for any 0 ≤ ℓ ≤ k + 1 it holds on a cell K with dimension h x × h y if v ∈ H ℓ (K). These estimates can be shown with the help of a pointwise interpolation operator and the L 2 -stability of the L 2 -projection, see [17], or directly by using the theory of Apel, [2]. Note that the L 2 -projection is also L ∞ -stable, see also [13]. The second operator I 2 utilises the classical interpolation operator J on D k . It is defined on each cell K for It holds the anisotropic interpolation error estimates for J , see [7,19], if v ∈ H k+1 (K). Note that for RT k an even sharper result involving only pure derivatives of v holds, see [7]. In addition, we have also anisotropic interpolation error estimates for the L 2 -norm of the divergence, see [7], if v is such that div v ∈ H k+1 (K) for RT k and div v ∈ H k (K) for BDM k .
If we only want to prove convergence in the L 2 -norm of U = (u, −ε grad u) the interpolation operator J is enough. For a stronger convergence result we define a more sophisticated operator, following ideas from [11]. Recalling the decomposition of u, we have for u = −ε grad u the decomposition with the obvious definition of the bold font letters. Now the operator I 2 is defined piecewise on each cell K ⊂ T N with id ∈ {1, 2, 3, 4, 12, 23, 34, 41} Using Γ id := ∂Ω id \ Γ we define the remaining operators on each K ⊂ Ω id \ Ω * id using the same definition as for J with the exception of the first condition in (8) and (9). This one is replaced by the two conditions For our analysis let us define a norm that is associated with B(·, ·). Here it holds that is equivalent to coercivity in the energy norm of the weak formulation of (1). But we can actually use the stronger norm . This norm is equivalent to the weighted H 2 -norm u L 2 (Ω) +ε grad u L 2 (Ω) + ε 2 ∆u L 2 (Ω) , which is stronger than the energy norm and, unfortunately, also not balanced. To repair this weakness we also introduce a balanced version of this norm The remainder of this section is devoted to proving optimal uniform convergence orders in the balanced norm. Of course, convergence in the unbalanced norm then follows.
Proof. By (13) we already have and together with δ ≤ c 0 In addition it holds Let us split the error U − U N into an interpolation error and a discrete error Using above inf-sup inequality and the Galerkin orthogonality (6) we arrive at and we are left with estimating B((η, η), V ) for any V ∈ U N . Here it holds using (5) due to I 1 being the weighted L 2 -projection. Note that in the case of constant c, the last term would also vanish due to div v| K ∈ Q k (K).
Proof. Using the solution decomposition (12) and the anisotropic interpolation error estimate (10) we obtain For the boundary layer terms we use the special structure of I 2 and estimate differently on the subdomains of Ω. We show the procedure for w 1 , the estimates of the other terms follow similarly. In Ω \ Ω 1 the interpolant is zero and we get In Ω 1 we have where P 1 :=Ĵ 1 − J . For the first term we obtain using (10) and (4) J due to σ > k +1. In the remaining ply of elements the operator P 1 in the Raviart-Thomas case is given by F P 1 w 1 · n · q = 0 for all faces F ⊂ ∂K \ Γ 1 , ∀q ∈ P k (F ), F P 1 w 1 · n · q = F w 1 · n · q for all faces F ⊂ ∂K ∩ Γ 1 , ∀q ∈ P k (F ), and similarly for the other finite elements. Thus P 1 w 1 depends only on w 1 · n| Γ 1 . With P 1 on Γ 1 being defined by weighted integrals, we have (15) Applying the same techniques to the other boundary and corner layer terms, and collecting the result finishes the first part of the proof. For the divergence we can apply the same techniques with the difference of applying (11) instead of (10). We obtain for σ > k + 1 and The last term to estimate is the error on the ply of elements in Ω 1 \Ω * 1 . A closer inspection of P 1 w 1 reveals (P 1 w 1 ) 2 = 0. Thus, an inverse inequality followed by (15) where h N/4 ≥ h min ≥ εN −1 holds due to the assumptions on φ. The analysis for the other terms of the decomposition follows the same lines. For D(K) = BDM k (K) the same analysis can be done, only replacing the convergence orders by k for σ ≥ k + 1/2.
Proof. Let c K := 1 meas(K) K c ≥ c 0 be a piecewise constant approximation of c. Now due to the weighted L 2 -projection and div v| K ∈ Q k (K). It holds Thus we obtain, using h x and h y as abbreviations for the dimensions of K, an inverse inequality and (7) for any v ∈ H k+1 (K) and similarly for the y-derivative of the second component. Let us start with the smooth part s of the solution decomposition and denote the coarse part of Ω by Ω c := Ω \ 4 i=1 Ω i and the union of corners by Ω cor := Ω 12 ∪ Ω 23 ∪ Ω 34 ∪ Ω 41 . Then we obtain due to h min ≥ ε −1 N and the condition on hε. In the final estimate on Ω cor we have used For ∂ y v 2 holds a similar estimate due to symmetry. Next we look at the boundary layer term w 1 . We obtain in Ω 1 using again the condition on εh. Now for the y-derivative it holds In the remainder of the domain we apply the L ∞ -stability and get The estimation of the other boundary layer terms and of the corner layer terms is similar. Combining all the individual results proves the assertion. Proof. Using the inf-sup estimate (14) and the previous lemmas we obtain for D(K) = RT k (K) which are the three components of |||ξ||| bal . Similar results follow for D(K) = BDM k (K) with the additional (h + N −1 max |ψ ′ |)(ln N) 1/2 1.
Proof. With the triangle inequality |||U − U N ||| bal ≤ |||(ξ, ξ)||| bal + |||(η, η)||| bal and the previous lemmas it only remains to estimate η L 2 (Ω) , which can be done using the local anisotropic interpolation error estimates (7) by standard techniques Corollary 3.6. Under the same conditions as the previous theorem we also have for D(K) = RT k (K) in the unbalanced norm Remark 3.7. On a Shishkin mesh we have h i = εN −1 ln N inside the boundary domain. Thus the condition on εh becomes On a Bakhvalov S-mesh we have h ∼ ε and the condition becomes Remark 3.8. The same analysis can also be conducted for the Arnold-Boffi-Falk element see [3], using div D(K) = Q k+1 (K) \ span{x k+1 y k+1 } as discrete space for the first component. Anisotropic interpolation error estimates are given in [7]. Although η L 2 (Ω) and div η L 2 (Ω) can be estimated with order k + 2, we obtain only convergence rates of order k + 1 due to η L 2 (Ω) ε 1/2 (h + N −1 max |ψ ′ |) k+1 .  and an exact solution is prescribed, see [1,17] for c = 1. The solution has only boundary layers at x = 0 and y = 0, and a corner layer at (0, 0). Therefore, we modify our mesh accordingly. For our experiments we will always use Bakhvalov-S-meshes. All computations were done in SOFE, a finite-element framework in Matlab and Octave, see github.com/SOFE-Developers/SOFE. Let us start the numerical investigation by looking at the dependence on ε. For that we fix N = 16 and use RT 1 -elements, and vary ε ∈ {10 −3 , 10 −4 , 10 −5 , 10 −6 }. We obtain the numbers in Table 1. We observe u − u h L 2 (Ω) to be independent of ε, as expected, while u − u h L 2 (Ω) ε 1/2 and div(u − u h ) L 2 (Ω) ε −1/2 , also both as expected. Consequently, |||U − U h ||| stays independent of ε, due to the dominating effect of u − u h L 2 (Ω) , and the larger balanced norm |||U − U h ||| bal is independent too due to the correct weighting of the other two norms. Now let us come to the convergence orders. For that purpose we fix ε = 10 −4 and vary for different values of k the number N of cells per dimension. We start with Raviart-Thomas elements and obtain the results of Table 2. Here along with the computed errors also the estimated rates of convergence are given and they are close to the expected rates of k + 1 for the balanced norm. In the case of Brezzi-Douglas-Marini elements we get Table 3. As expected we only see rates of k in the balanced version with slightly better results for the lowest order case. The reason for this behaviour lies in the components of the balanced norms, where the faster converging ones dominate for smaller values of N the balanced norm. A closer look reveals u − u h L 2 (Ω) and div(u − u h ) L 2 (Ω) only to be convergent with order 1, see Table 4.