Mixed methods for degenerate elliptic problems and application to fractional laplacian

TL;DR

This paper extends mixed finite element methods to degenerate elliptic equations with coefficients in the Muckenhoupt class, providing error analysis and applying it to fractional Laplacian problems.

Contribution

It develops a generalized error analysis for mixed finite element methods with degenerate coefficients and applies it to fractional Laplacian equations.

Findings

Extended error analysis to coefficients in Muckenhoupt class A2.

Obtained optimal error estimates for Raviart-Thomas spaces.

Applied results to fractional Laplace equation.

Abstract

We analyze the approximation by mixed finite element methods of solutions of equations of the form $-\mbox{div\,} (a\nabla u) = g$, where the coefficient $a=a(x)$ can degenerate going to cero or infinity. First, we extend the classic error analysis to this case provided that the coefficient $a$ belongs to the Muckenhoupt class $A_2$. The analysis developed applies to general mixed finite element spaces satisfying the standard commutative diagram property, whenever some stability and interpolation error estimates are valid in weighted norms. Next, we consider in detail the case of Raviart-Thomas spaces of lowest order, obtaining optimal order error estimates for general regular elements as well as for some particular anisotropic ones which are of interest in problems with boundary layers. Finally we apply the results to a problem arising in the solution of the fractional Laplace equation.

Equations

$$\left\{\begin{array}[]{rl}-\mbox{div\,}(a\nabla u)=g&\ \ \mbox{in}\ {\mathcal{D}}\\ u=0&\ \ \mbox{on}\ \Gamma_{D}\\ -a\nabla u\cdot{\boldsymbol{n}}=f&\ \ \mbox{on}\ \Gamma_{N}\end{array}\right.$$

$$[a{]}_{A_2}:=\underset{Q}{\sup }\left(\frac{1}{|Q|}{\int }_{Q}a\right)\left(\frac{1}{|Q|}{\int }_Q{a}^{-1}\right)<\infty ,$$

$$[a{]}_{A_2^s}:=\underset{R}{\sup }\left(\frac{1}{|R|}{\int }_{R}a\right)\left(\frac{1}{|R|}{\int }_R{a}^{-1}\right)<\infty .$$

$$H_a^1(S)=\left\{v\in L_a^2(S):|\nabla v|\in L_a^2(S)\right\}$$

$$\left\{\begin{array}[]{rl}{\boldsymbol{\sigma}}+a\nabla u=0&\ \ \mbox{in}\ {\mathcal{D}}\\ \mbox{div\,}{\boldsymbol{\sigma}}=g&\ \ \mbox{in}\ {\mathcal{D}}\\ u=0&\ \ \mbox{on}\ \Gamma_{D}\\ {\boldsymbol{\sigma}}\cdot{\boldsymbol{n}}=f&\ \ \mbox{on}\ \Gamma_{N}\end{array}\right.$$

$$H(\text{mbox}div,\mathcal{D})=\{\mathbf{\tau }\in L^2(\mathcal{D}{)}^n:\text{mbox}div\mathbf{\tau }\in L^2(\mathcal{D})\}$$

$$\Vert \mathbf{\tau }{\Vert }_{H(\text{mbox}div,\mathcal{D})}^2=\Vert \mathbf{\tau }{\Vert }_{L^2(\mathcal{D})}^2+\Vert \text{mbox}div\mathbf{\tau }{\Vert }_{L^2(\mathcal{D})}^2.$$

$$H_{{\Gamma }_N}(\text{mbox}div,\mathcal{D})=\{\mathbf{\tau }\in H(\text{mbox}div,\mathcal{D}):\mathbf{\tau }\cdot \mathbf{n}=0\text{mbox}on{\Gamma }_N\}.$$

$$a^{-1}\mathbf{\sigma }+\nabla u=0\text{mbox}in\mathcal{D},$$

$$\mathbf{\sigma }\cdot \mathbf{n}=f\text{mbox}on{\Gamma }_N$$

$$\{\begin{array}{rlr}{\int }_{\mathcal{D}}{a}^{-1}\mathbf{\sigma }\cdot \mathbf{\tau }-{\int }_{\mathcal{D}}u\text{mbox}div\mathbf{\tau }& =0& \forall \mathbf{\tau }\in H_{{\Gamma }_N}(\text{mbox}div,\mathcal{D})\\ {\int }_{\mathcal{D}}v\text{mbox}div\mathbf{\sigma }& ={\int }_{\mathcal{D}}gv& \forall v\in L^2(\mathcal{D})\end{array}$$

$${\mathbf{S}}_{h,N}={\mathbf{S}}_h\cap H_{{\Gamma }_N}(\text{mbox}div,\mathcal{D}),$$

$$\text{mbox}div{\mathbf{S}}_{h,N}=V_h$$

$${\int }_{\mathcal{D}}\text{mbox}div(\mathbf{\tau }-{\Pi }_h\mathbf{\tau })v=0\forall v\in V_h.$$

$$\text{mbox}div{\Pi }_h=P_h\text{mbox}div.$$

$$({\mathbf{\sigma }}_h,u_h)\in {\mathbf{S}}_h\times V_h$$

$${\mathbf{\sigma }}_h\cdot \mathbf{n}={\Pi }_h\mathbf{\sigma }\cdot \mathbf{n}\text{mbox}on{\Gamma }_N$$

$$\{\begin{array}{rlr}{\int }_{\mathcal{D}}{a}^{-1}{\mathbf{\sigma }}_h\cdot \mathbf{\tau }-{\int }_{\mathcal{D}}{u}_h\text{mbox}div\mathbf{\tau }& =0& \forall \mathbf{\tau }\in {\mathbf{S}}_{h,N},\\ {\int }_{\mathcal{D}}v\text{mbox}div{\mathbf{\sigma }}_h,& ={\int }_{\mathcal{D}}gv& \forall v\in V_h.\end{array}$$

$$\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{L_{a^{-1}}^2}\le \Vert \mathbf{\sigma }-{\Pi }_h\mathbf{\sigma }{\Vert }_{L_{a^{-1}}^2}.$$

$${\int }_{\mathcal{D}}\text{mbox}div({\Pi }_h\mathbf{\sigma }-{\mathbf{\sigma }}_h)v=0\forall v\in V_h.$$

$${\int }_{\mathcal{D}}{a}^{-1}(\mathbf{\sigma }-{\mathbf{\sigma }}_h)\cdot ({\Pi }_h\mathbf{\sigma }-{\mathbf{\sigma }}_h)=0$$

$$\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{L_{a^{-1}}^2}^2$$

$$\le \Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{L_{a^{-1}}^2}\Vert \mathbf{\sigma }-{\Pi }_h\mathbf{\sigma }{\Vert }_{L_{a^{-1}}^2},$$

$$\text{mbox}div\mathbf{\tau }=\varphi$$

$$\Vert \mathbf{\tau }{\Vert }_{H_{a^{-1}}^1(\mathcal{D})}\le C\Vert \varphi {\Vert }_{L_{a^{-1}}^2(\mathcal{D})},$$

$$\tilde{\varphi }=\{\begin{array}{rl}\varphi (x)& ,x\in \mathcal{D}\\ -\frac{{\int }_{\mathcal{D}}\varphi }{|\tilde{\mathcal{D}}\setminus \mathcal{D}|}& ,x\in \tilde{\mathcal{D}}\setminus \mathcal{D}\end{array}$$

$$\Vert \mathbf{\tau }{\Vert }_{H_{a^{-1}}^1(\tilde{\mathcal{D}})}\le C\Vert \tilde{\varphi }{\Vert }_{L_{a^{-1}}^2(\tilde{\mathcal{D}})}$$

$$\Vert {\Pi }_h\mathbf{\tau }{\Vert }_{L_{a^{-1}}^2}\le C\Vert \mathbf{\tau }{\Vert }_{H_{a^{-1}}^1}.$$

$$\Vert u-u_h{\Vert }_{L_a^2}\le \Vert u-P_{h}u{\Vert }_{L_a^2}+C\Vert \mathbf{\sigma }-{\mathbf{\sigma }}_h{\Vert }_{L_{a^{-1}}^2},$$

$$\text{mbox}div\mathbf{\tau }=(P_{h}u-u_h)a$$

$$\Vert \mathbf{\tau }{\Vert }_{H_{a^{-1}}^1}\le C\Vert (P_{h}u-u_h)a{\Vert }_{L_{a^{-1}}^2}.$$

Full text

Supported by ANPCyT under grant PICT 2014-1771, by CONICET under grant 11220130100006CO and by Universidad de Buenos Aires under grant 20020120100050BA. The first author has a fellowship from CONICET, Argentina.

1. Introduction

In this paper we analyze the approximation by mixed finite element methods of degenerate second order elliptic problems. There is a vast bibliography concerning this kind of methods (see for example the books [7, 6] and references therein). However, as far as we know, only very few papers have considered the degenerate case (we can mention [5, 22]).

Let ${\mathcal{D}}\subset{\mathbb{R}}^{n}$ be a bounded Lipschitz polytope and $a$ be a non-negative measurable function. We assume that the boundary is decomposed into two disjoint parts $\Gamma_{D}$ and $\Gamma_{N}$. Given $g\in L^{2}({\mathcal{D}})$ and $f\in L^{2}(\Gamma_{N})$ we consider the problem

[TABLE]

where ${\boldsymbol{n}}$ denotes the unit exterior normal vector. If $\Gamma_{N}=\partial{\mathcal{D}}$ we assume the usual compatibility condition $\int_{\mathcal{D}}g=\int_{\partial{\mathcal{D}}}f$.

We have written the problem in this form in order to simplify notation. However, it is easy to see that all our arguments apply to general problems where the coefficient $a$ is replaced by a matrix $A=A(x)$ satisfying $\lambda a(x)|\xi|^{2}\leq\xi^{T}A(x)\xi\leq\Lambda a(x)|\xi|^{2}$, for all $x\in{\mathcal{D}}$, where $\lambda$ and $\Lambda$ are positive constants.

We are interested in degenerate problems in the sense that the coefficient $a$ can become infinite or zero in subsets of $\overline{\mathcal{D}}$ with vanishing $n-$dimensional measure. We will assume that $a$ belongs to the Muckenhoupt class $A_{2}$, in particular $a^{-1}\in L^{1}_{loc}({\mathcal{D}})$ and, therefore, the usual mixed method is well defined.

Recall that a non-negative measurable function $a\in L^{1}_{loc}({\mathbb{R}}^{n})$ belongs to $A_{2}$ if

[TABLE]

where the supremum is taken over all cube $Q$ with faces parallel to the coordinate axes.

The class $A_{2}$ was introduced to characterize the weights for which the Hardy-Littlewood maximal operator is bounded in the associated weighted norm (See for instance [9, 23]). After that, it was used in the theory of elliptic equations (see for example the pioneering work [15]) and, more recently, in the analysis of finite element approximations [4, 24, 25].

When dealing with anisotropic estimates we will work with the more restrictive strong $A_{2}$ class, which will be denoted by $A^{s}_{2}$ and is defined by

[TABLE]

where the supremum is taken now over all $n$-dimensional rectangles with faces parallel to the coordinate axes. It is known that $a\in A^{s}_{2}$ if and only if $a$ belongs to $A_{2}$ of one variable for each variable, uniformly in the other variables (see [17, 20]).

Given a weight $a$, for any measurable set $S$ we will denote with $L_{a}^{2}(S)$ the usual Hilbert space with measure $a\,dx$. We will also work with the weighted Sobolev space

[TABLE]

with its natural norm. We will omit the domain in these notations when it is clear from the context.

Under appropriate assumptions on $a$ (particularly if $a\in A_{2}$) and the data $f$ and $g$, it is possible to prove by standard arguments that there exists a unique solution of problem (1.1) belonging to $H_{a}^{1}({\mathcal{D}})$.

Introducing the variable vector field ${\boldsymbol{\sigma}}=-a\nabla u$, problem (1.1) can be transformed into the equivalent first order system

[TABLE]

Then, mixed finite element methods are based on a weak formulation of this system and they approximate simultaneously ${\boldsymbol{\sigma}}$ and $u$. One motivation for using this type of methods is that, in many applications, the variable of physical interest is ${\boldsymbol{\sigma}}$ and, therefore, it might be more efficient to approximate it directly instead of obtaining it from a computed approximation of $u$. A typical example of this situation is the Darcy equation arising in the simulation of flows in porous media. Indeed, it is many times argued that ${\boldsymbol{\sigma}}$ is smoother than $\nabla u$. Although this is probably true in practice, it is not possible to give a mathematical foundation to this statement in general (see [16] for an interesting discussion on this subject).

As an application of our results we will consider a problem arising in the solution of the fractional Laplace equation $(-\Delta)^{s}v=f$. As we will show, in the case $0<s<1/2$, the mixed method is more convenient than the standard one in the sense that almost optimal order of convergence can be obtained with a weaker grading of the meshes.

The rest of the paper is organized as follows. In Section 2 we recall the mixed finite element method for Problem (1.1) and extend the classic error analysis to the case of degenerate problems. A fundamental tool is the existence of right inverses of the divergence in weighted norms when the weight belongs to the class $A_{2}$. The analysis given in this section can be applied to general mixed finite element spaces which satisfy the so called commutative diagram property whenever a stability property in a weighted norm for the interpolation operator is valid. Next, in Section 3, we consider the case of Raviart-Thomas elements of lowest order and prove the stability property mentioned above and error estimates in weighted norms under the regularity assumption on the family of meshes. Then, in Section 4, we continue the analysis for the Raviart-Thomas spaces of lowest order and prove some weighted interpolation error estimates, where the weights involve the distance to some part of the boundary, for anisotropic rectangular and prismatic elements which are of interest in problems with boundary layers. An important tool in this part of the analysis is the so called improved Poincaré inequality. Finally, in Section 5, we consider the approximation of the fractional Laplace equation which leads to a particular degenerate problem of the type considered in the previous sections. We show in this example how the weighted error estimates proved for anisotropic elements can be used to design a priori adapted meshes giving almost optimal order with respect to the number of degrees of freedom. We include in this section some numerical results.

2. Mixed finite element approximations

First we recall some usual notation and known results on mixed methods. The appropriate space for the vector variable is

[TABLE]

which is a Hilbert space with norm given by

[TABLE]

Moreover, since in the mixed formulation Neumann type boundary conditions are imposed in an essential way, we will work with the subspace

[TABLE]

Dividing by $a$, the first equation in (1.2) can be rewritten as

[TABLE]

and multiplying by test functions and integrating by parts, we obtain the standard weak mixed formulation of problem (1.2), namely, find ${\boldsymbol{\sigma}}\in H(\mbox{div\,},{\mathcal{D}})$ and $u\in L^{2}({\mathcal{D}})$ such that

[TABLE]

and

[TABLE]

Observe that the Dirichlet boundary condition is implicit in the weak formulation. When $\Gamma_{N}=\partial{\mathcal{D}}$, $L^{2}({\mathcal{D}})$ has to be replaced by $L^{2}_{0}({\mathcal{D}})$, the subspace of functions with vanishing mean value.

As usual, the error analysis is divided in two steps. The first one consists in proving estimates for the finite element approximation error in terms of the error for some appropriate interpolation or projection operator. This part of the analysis can be done for general mixed finite element spaces provided they satisfy the so called commutative diagram property as well as some weighted stability estimates for the appropriate projections. Therefore, we will develop this part of the error analysis for general spaces stating the necessary assumptions that afterwards have to be proved for each particular choice of approximation spaces. The second part consists in estimating the interpolation error. For simplicity, we will restrict this analysis to the lowest order Raviart-Thomas elements. Higher order elements as well as other approximation spaces could be treated similarly but this require non trivial technical modifications.

We assume that we have a family of partitions $\{{\mathcal{T}}_{h}\}$ of the domain ${\mathcal{D}}$ such that each ${\mathcal{T}}_{h}$ is consistent with the boundary conditions, i. e., the exterior boundary of an element is completely contained in $\Gamma_{D}$ or in $\Gamma_{N}$. Associated with these partitions we assume that we have finite element spaces ${\boldsymbol{S}}_{h}\subset H(\mbox{div\,},{\mathcal{D}})$, $V_{h}\subset L^{2}({\mathcal{D}})$ (or $V_{h}\subset L^{2}_{0}({\mathcal{D}})$ when $\Gamma_{N}=\partial{\mathcal{D}}$), such that, if

[TABLE]

then

[TABLE]

and there exists an operator $\Pi_{h}:{\boldsymbol{S}}\longrightarrow{\boldsymbol{S}}_{h}$, defined in an appropriate subspace ${\boldsymbol{S}}\subset H(\mbox{div\,},{\mathcal{D}})$ containing the solution ${\boldsymbol{\sigma}}$, such that, if ${\boldsymbol{\tau}}\in{\boldsymbol{S}}\cap H_{\Gamma_{N}}(\mbox{div\,},{\mathcal{D}})$ then $\Pi_{h}{\boldsymbol{\tau}}\in{\boldsymbol{S}}_{h,N}$ and, for all ${\boldsymbol{\tau}}\in{\boldsymbol{S}}$,

[TABLE]

Introducing the $L^{2}$-orthogonal projection $P_{h}:L^{2}({\mathcal{D}})\longrightarrow V_{h}$, (2.3) and (2.4) yield the commutative diagram property

[TABLE]

The mixed finite element approximation of problem (1.2) is given by

[TABLE]

such that,

[TABLE]

and

[TABLE]

Existence and uniqueness of the discrete solution and the following error estimate follow by well known arguments (see for example [7, 6]). For completeness we include the proof of the error estimate to show that the usual arguments can be adapted for degenerate problems and for the mixed boundary conditions considered here. We neglect numerical integration errors assuming that all the integrals can be computed exactly.

Lemma 2.1**.**

Assume that $a^{-1}\in L^{1}({\mathcal{D}})$ and ${\boldsymbol{\sigma}}\in L^{2}_{a^{-1}}({\mathcal{D}})^{n}$. If ${\boldsymbol{\sigma}}$ is the solution of (2.1) and (2.2)and ${\boldsymbol{\sigma}}_{h}$ that of (2.6) and (2.7), then

[TABLE]

Proof.

Subtracting the second equation in (2.7) to the second one in (2.2) and using (2.4) we obtain

[TABLE]

From (2.6) it follows that $\Pi_{h}{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\in{\boldsymbol{S}}_{h,N}$, and then, by (2.3) we conclude that $\mbox{div\,}(\Pi_{h}{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h})=0$. Moreover, taking ${{\boldsymbol{\tau}}}=\Pi_{h}{{\boldsymbol{\sigma}}}-{{\boldsymbol{\sigma}}}_{h}$ in (2.2) and (2.7), we obtain

[TABLE]

and so,

[TABLE]

and the lemma is proved.∎

To estimate the error in the approximation of the scalar variable we need a stronger assumption on the coefficient $a$. Indeed, we will prove the following result that generalizes to the weighted case the existence of continuous right inverses of the divergence.

Lemma 2.2**.**

If $a\in A_{2}$ then, given $\phi\in L^{2}_{a^{-1}}({\mathcal{D}})$ (satisfying $\int_{\mathcal{D}}\phi=0$ in the case $\Gamma_{N}=\partial{\mathcal{D}}$), there exists ${\boldsymbol{\tau}}\in H_{a^{-1}}^{1}({\mathcal{D}})^{n}\cap H_{\Gamma_{N}}(\mbox{div\,},{\mathcal{D}})$ such that

[TABLE]

and

[TABLE]

where the constant $C$ depends on ${\mathcal{D}}$ and $a$.

Proof.

In the case $\Gamma_{N}=\partial{\mathcal{D}}$ we have $\int_{\mathcal{D}}\phi=0$ and the result is known. Indeed, for domains which are star-shaped with respect to a ball it was proved in [14, Th. 3.1] and [26, Th.1.1] using Bogovskii’s solution of the divergence and the theory of singular integrals. The arguments used there can be extended for the class of John domains using the generalization of Bogovskii’s operator introduced in [3] (For more details see also [2]). A different proof was given in [10, Th. 5.2] also for the class of John domains.

Suppose now that $\Gamma_{N}\neq\partial{\mathcal{D}}$. Enlarging the domain in an appropriate way we can obtain a Lipschitz domain $\widetilde{\mathcal{D}}$ such that ${\mathcal{D}}\varsubsetneq\widetilde{\mathcal{D}}$ and $\Gamma_{N}\subset\partial\widetilde{\mathcal{D}}$. For example, we can make a smooth deformation of part of $\Gamma_{D}$.

Now, we extend $\phi$ to $\widetilde{\mathcal{D}}$ as

[TABLE]

and then, since $\int_{\widetilde{\mathcal{D}}}\widetilde{\phi}\,dx=0$, there exists ${\boldsymbol{\tau}}\in H_{a^{-1}}^{1}(\widetilde{\mathcal{D}})^{n}$, vanishing on $\partial\widetilde{\mathcal{D}}$ and satisfying

[TABLE]

It is easy to see that $\|\widetilde{\phi}\|_{L_{a^{-1}}^{2}(\widetilde{\mathcal{D}})}\leq C\|\phi\|_{L_{a^{-1}}^{2}({\mathcal{D}})}$, and, therefore, the restriction of ${\boldsymbol{\tau}}$ to ${\mathcal{D}}$ satisfies the required properties. ∎

For the next lemma we need to use the following stability result in a weighted norm:

[TABLE]

Assuming that $a\in A_{2}$, we will prove this estimate for the lowest order Raviart-Thomas spaces in a forthcoming section.

Lemma 2.3**.**

Let $({\boldsymbol{\sigma}},u)$ and $({\boldsymbol{\sigma}}_{h},u_{h})$ be the solutions of (2.1) and (2.2), and (2.6) and (2.7) respectively. If $a\in A_{2}$ and $\Pi_{h}$ satisfies (2.8) then

[TABLE]

where $C$ depends on the constant in Lemma 2.2.

Proof.

Assume first that $\Gamma_{D}\neq\emptyset$. According to Lemma 2.2 there exists ${\boldsymbol{\tau}}\in H_{a^{-1}}^{1}({\mathcal{D}})^{n}\cap H_{\Gamma_{N}}(\mbox{div\,},{\mathcal{D}})$ such that

[TABLE]

and

[TABLE]

Then,

[TABLE]

where we have used (2.4) and (2.8). Then, (2.9) follows by the triangular inequality.

Now, if $\Gamma_{N}=\partial{\mathcal{D}}$, there exists ${\boldsymbol{\tau}}\in H_{\Gamma_{N}}(\mbox{div\,},{\mathcal{D}})$ such that

[TABLE]

where $\overline{(P_{h}u-u_{h})a}$ denotes the average of $(P_{h}u-u_{h})a$, and

[TABLE]

Indeed, this follows from Lemma 2.2 and the estimate

[TABLE]

Since $\int_{\mathcal{D}}(P_{h}u-u_{h})=0$ we have

[TABLE]

The rest of the argument follows as in the previous case. ∎

Combining Lemmas 2.1 and 2.3 we obtain the following

Corollary 2.4**.**

Under the same hypotheses of Lemma 2.3 we have

[TABLE]

3. Error estimates for Raviart-Thomas elements

We now consider the approximation by the lowest order Raviart-Thomas mixed finite elements. To apply the results obtained in the previous section we have to prove error estimates for the corresponding operators $\Pi_{h}$ and $P_{h}$.

Recall that the local Raviart-Thomas space of lowest degree for a simplex $T$ is

[TABLE]

while for $R$ an $n$-dimensional rectangular element with faces parallel to the coordinate axes, is

[TABLE]

Then, the global space for the mixed approximation of the vector variable for a partition ${\mathcal{T}}_{h}$ made of any kind of elements is

[TABLE]

The associated space for the scalar variable is given by the piecewise constant functions, namely,

[TABLE]

where $V=L^{2}_{0}({\mathcal{D}})$ when $\Gamma_{N}=\partial{\mathcal{D}}$ or $V=L^{2}({\mathcal{D}})$ otherwise. Then, the projection $P_{h}$ is given by $(P_{h}v)|_{K}=P_{K}v:=\frac{1}{|K|}\int_{K}v$. The fundamental tool for the error analysis is the well known Raviart-Thomas operator $\Pi_{h}$ defined on each element $K$ as $(\Pi_{h}{\boldsymbol{\tau}})|_{K}=\Pi_{K}{\boldsymbol{\tau}}$ where

[TABLE]

for all face $F$ of $K$ where ${\boldsymbol{n}}_{F}$ denotes a unitary vector normal to $F$ (here $K$ may be a simplex or a rectangle). This operator is well defined whenever the ${\boldsymbol{\tau}}\cdot{\boldsymbol{n}}\in L^{1}(F)$ which is known to be true for any ${\boldsymbol{\tau}}\in W^{1,1}(K)^{n}$. Moreover, it is not difficult to check that (2.3), (2.4), and consequently (2.5), are satisfied.

We consider first the case of regular partitions, namely, if $h_{K}$ and $\rho_{K}$ are the diameters of $K$ and the biggest ball contained in $K$ respectively, we assume that the family of meshes $\{{\mathcal{T}}_{h}\}$ satisfy $h_{K}/\rho_{K}\leq\eta$ with a constant $\eta$ independent of $h$.

Basic tools for interpolation error estimates are the Poincaré type inequalities. Given a set $S$ and a function $v$ we will denote with $v_{S}$ the average of $v$ over $S$ (both for $S=K$ or some face of $K$). In what follows the constant $C$ will depend on the weight $[a]_{A_{2}}$, although it is possible to give an explicit bound for this dependence this is not of interest for our purposes because we will work with a fixed $a$. Let us also remark that the arguments given below can be applied to obtain analogous interpolation error estimates in $L_{a}^{p}(K)$, $1\leq p<\infty$, provided $a\in A_{p}$ (see, for example, [12] for the definition of these classes).

To simplify notation we will prove all the estimates for the weight $a$ although some of them will be used later for $a^{-1}$. Note that, from the definition of $A_{2}$, it follows immediately that $a\in A_{2}$ if and only if $a^{-1}\in A_{2}$.

In what follows we will use the following observation: under the regularity assumption it is easy to see that

[TABLE]

with $C$ depending only on $n$ and $\eta$.

Lemma 3.1**.**

For $a\in A_{2}$ there exists a constant depending only on $a$ and $\eta$, such that,

[TABLE]

and

[TABLE]

Proof.

We have

[TABLE]

where we have used the Schwarz inequality in the last step. Therefore, (3.5) follows from (3.4).

On the other hand, (3.6) is the well known weighted Poincaré inequality. It was first proved in [15] for the case of a ball and extended for very general domains in several papers (see, for example, [8, 14, 19]). The dependence of the constant on $h_{T}$ can be obtain by usual scaling arguments. ∎

Observe that, since $a\in A_{2}$ then $a^{-1}\in L^{1}(K)$ and, therefore,

[TABLE]

Consequently $H^{1}_{a}(K)\subset W^{1,1}(K)$ and, in particular, traces of functions in $H^{1}_{a}(K)$ on a face $F$ are well defined and belong to $L^{1}(F)$.

Our error estimates are based on the following generalized Poincaré inequality.

Lemma 3.2**.**

Given $a\in A_{2}$, a simplex or a rectangle $K$ and a face $F$ of $K$, there exists a constant $C$, depending only on $a$ and the regularity constant $\eta$, such that

[TABLE]

for all $v\in H^{1}_{a}(K)$.

Proof.

In view of (3.6) it is enough to estimate $\|v_{K}-v_{F}\|_{L^{2}_{a}(K)}$.

As we have mentioned above $H^{1}_{a}(K)\subset W^{1,1}(K)$, in particular $v_{F}$ is well defined. Writing now

[TABLE]

and using the classic trace theorem

[TABLE]

combined with the classic Poincaré inequality in $L^{1}(K)$, we obtain

[TABLE]

and consequently,

[TABLE]

and so, in view of (3.4), the lemma is proved. ∎

We can now prove the error estimates for the Raviart-Thomas interpolation. We will denote with $D{\boldsymbol{\sigma}}$ the differential matrix of ${\boldsymbol{\sigma}}$.

Lemma 3.3**.**

Given $a\in A_{2}$ and a simplex or a rectangle $K$ there exists a constant depending only on $a$ and the regularity constant $\eta$, such that

[TABLE]

Proof.

Consider first the case of a simplex. We choose three faces $F_{i}$ with corresponding normals ${\boldsymbol{n}}_{i}$.

From (3.3) we have $\int_{F_{i}}({\boldsymbol{\sigma}}-\Pi_{K}{\boldsymbol{\sigma}})\cdot{\boldsymbol{n}}_{i}=0$ and, therefore, using Lemma 3.2 we obtain,

[TABLE]

But, $\frac{\partial(\Pi_{K}{\boldsymbol{\sigma}})_{j}}{\partial x_{k}}=0$ for $k\neq j$, while

[TABLE]

where we have used the commutative diagram property. Then, in view of (3.5) we obtain

[TABLE]

Now observe that,

[TABLE]

with a constant depending only on $\eta$, and so, (3.7) follows from (3.8) and (3.9).

For a rectangular element we proceed in the same way. The only difference is that to prove (3.9) we use now (3.5) combined with

[TABLE]

∎

Combining the error estimates obtained above with the results of the previous section we can now state the main theorem for approximation by Raviart-Thomas of lowest order on regular families of meshes.

Theorem 3.4**.**

Let ${\mathcal{T}}_{h}$ be a family of meshes with regularity constant $\eta$ and $h=\max_{K\in{\mathcal{T}}_{h}}h_{K}$. If $({\boldsymbol{\sigma}},u)$ and $({\boldsymbol{\sigma}}_{h},u_{h})$ are the solutions of (2.1) and (2.2), and (2.6) and (2.7) respectively then, for $a\in A_{2}$, there exists a constant $C$ depending only on ${\mathcal{D}}$, $a$ and $\eta$ such that

[TABLE]

and

[TABLE]

Proof.

The error estimate for ${\boldsymbol{\sigma}}$ follows from Lemma 2.1 combined with the estimate (3.7) applied to the weight $a^{-1}$ (recall that $a\in A_{2}$ if and only if $a^{-1}\in A_{2}$).

On the other hand, observe that (3.7) implies the hypothesis (2.8) assumed in Lemma 2.3. Then, to bound the error for $u$ we apply that lemma, (3.7) again, and (3.6). ∎

4. Anisotropic error estimates

Our next goal is to prove anisotropic error estimates suitable for problems with boundary layers. For this kind of problems it is useful to have estimates involving a weighted norm on the right hand side where the weight is a power of the distance to some part of the boundary.

To present the main arguments we consider first the case of rectangular elements. Then we show how similar ideas can be applied to prismatic elements which are of interest in the application that we are going to consider in the next section, and more generally, in many problems with solutions presenting boundary layers. The case of simplex can be treated in a similar way but, as in the un-weighted case, anisotropic error estimates are valid only for some particular kind of degenerate elements (see [1]).

Proceeding as in the previous section, we need now the following weighted improved Poincaré inequality, which is well known (see, for example, [18, 11]). For $a\in A_{2}$ and $Q$ a cube,

[TABLE]

where $d$ denotes the distance to $\partial Q$. Consider an arbitrary rectangle

[TABLE]

If we replace $Q$ by $R$ in the above inequality, it is known that the constant in (4.1) blows up when the ratio between outer and inner diameter goes to infinity. However, we have the following anisotropic version if the weight belongs to the smaller class $A^{s}_{2}$ defined in the introduction. For $i=1,\cdots,n$ we define

[TABLE]

Lemma 4.1**.**

For $a\in A^{s}_{2}$,

[TABLE]

Proof.

It follows immediately from (4.1) that, if $Q$ is the unitary cube

[TABLE]

Then, (4.2) follows by standard arguments making the change of variables $x_{i}=h_{i}\hat{x}_{i}+a_{i}$ and using that, for $\hat{a}(\hat{x}):=a(x)$, $a\in A_{2}^{s}\Longrightarrow\hat{a}\in A_{2}^{s}$. ∎

Lemma 4.2**.**

For $a\in A^{s}_{2}$ and $F$ the face contained in $x_{j}=a_{j}$ we have

[TABLE]

Proof.

By a simple integration by parts in the $x_{j}$ variable we have

[TABLE]

Then,

[TABLE]

and therefore,

[TABLE]

but, multiplying and dividing by $a^{1/2}$ and using the Schwarz inequality we obtain

[TABLE]

and consequently,

[TABLE]

Therefore, (4.3) follows from (4.2). ∎

We can now prove anisotropic error estimates for the Raviart-Thomas interpolation $\Pi_{R}$. Observe that each component $(\Pi_{R}{\boldsymbol{\sigma}})_{j}$ depends only on ${\boldsymbol{\sigma}}_{j}$, and so, to simplify notation we will write simply $\Pi_{R}{\boldsymbol{\sigma}}_{j}$.

Lemma 4.3**.**

For $a\in A^{s}_{2}$ and $1\leq j\leq n$,

[TABLE]

Proof.

Since ${\boldsymbol{\sigma}}_{j}-\Pi_{R}{\boldsymbol{\sigma}}_{j}$ has vanishing mean value on the face defined by $x_{j}=a_{j}$ we obtain from (4.3),

[TABLE]

But, for $i\neq j$, $\frac{\partial(\Pi_{R}{\boldsymbol{\sigma}})_{j}}{\partial x_{i}}=0$. On the other hand from the definition of $\Pi_{R}$ we have

[TABLE]

and a simple argument using the Schwarz inequality shows that, for any $v\in L^{2}_{a}(R)$,

[TABLE]

and therefore,

[TABLE]

and the lemma is proved. ∎

Now we analyze the case of prismatic elements. For notational convenience we work in ${\mathbb{R}}^{n+1}$ and introduce the variables $(x,y)$, with $x=(x_{1},\cdots,x_{n})\in{\mathbb{R}}^{n}$ and $y\in{\mathbb{R}}$. Therefore, the class $A^{s}_{2}$ denotes now the class of weights satisfying

[TABLE]

where the supremum is taken over all $n+1$-dimensional rectangles.

We consider elements $P=K\times[y_{0},y_{1}]$ where $K$ is an $n$-dimensional simplex and $y_{j}\in{\mathbb{R}}$ for $j=0,1$.

Similar arguments than those used above for the anisotropic estimates in rectangular elements can be used in this case. To simplify notation we will prove only the particular weighted estimates that we will need for the application considered in the next section. We will denote by $h_{K}$ the diameter of $K$. The elements considered are anisotropic because no relation between $h_{K}$ and $y_{1}-y_{0}$ is required. On the other hand, for the simplices we assume the regularity condition $h_{K}/\rho_{K}\leq\eta$.

Lemma 4.4**.**

Given $a\in A^{s}_{2}$, $P=K\times[y_{0},y_{1}]$ a prismatic element, and $F_{P}$ a face of $P$ given by $F_{P}:=F\times[y_{0},y_{1}]$, where $F$ is a face of $K$, we have

[TABLE]

Proof.

Proceeding as in the proof of (4.2) we can prove the Poincaré type inequality

[TABLE]

We will denote with $dS_{F}$ and $dS_{F_{P}}$ the surface measures on $F$ and $F_{P}$ respectively. Calling $x_{F}$ the vertex of $K$ opposite to $F$ and integrating by parts we have,

[TABLE]

but, for $x\in F$, $(x-x_{F})\cdot{\boldsymbol{n}}_{F}=n|K|/|F|$, and therefore,

[TABLE]

Then, integrating in the variable $y$,

[TABLE]

and dividing this equation by $(y_{1}-y_{0})$ we obtain

[TABLE]

which, using (4.6) and proceeding as in the last part of the proof of Lemma 4.2, implies (4.5). ∎

Lemma 4.5**.**

Given $a\in A^{s}_{2}$, $P=K\times[y_{0},y_{1}]$ a prismatic element, and $F_{P}$ a face of $P$ given by $F_{P}:=K\times\{y_{j}\}$, $j=0$ or $1$, we have

[TABLE]

Proof.

It is analogous to the proof of Lemma 4.2. ∎

The local Raviart-Thomas space for $P=K\times[y_{0},y_{1}]$ is given by

[TABLE]

Given a vector field ${\boldsymbol{\sigma}}$ we define $\tilde{\boldsymbol{\sigma}}=(\sigma_{1},\cdots,\sigma_{n})$ and write ${\boldsymbol{\sigma}}=(\tilde{\boldsymbol{\sigma}},\sigma_{n+1})$. Since the normals to the top and bottom faces of $P$ are orthogonal to the other ones, the Raviart-Thomas interpolation can be written as

[TABLE]

where $\Pi_{K}$ and $\Pi_{n+1}$ depend on $\tilde{\boldsymbol{\sigma}}$ and $\sigma_{n+1}$ respectively. Indeed, they are defined by

[TABLE]

for all face $F$ of $K$ and

[TABLE]

for $j=1,2$.

Lemma 4.6**.**

For $a\in A^{s}_{2}$ and $P=K\times[y_{0},y_{1}]$, we have

[TABLE]

and

[TABLE]

where $C$ depends only on $a$ and the regularity constant $\eta$.

Proof.

Since $(\tilde{\boldsymbol{\sigma}}-\Pi_{K}\tilde{\boldsymbol{\sigma}})\cdot{\boldsymbol{n}}_{F}$ has vanishing mean value on $F_{P}=F\times[y_{0},y_{1}]$ we can apply (4.5) to obtain

[TABLE]

and using this estimate for $n$ different faces of $K$ together with the regularity assumption, we arrive at

[TABLE]

But $\frac{\partial(\Pi_{K}\tilde{\boldsymbol{\sigma}})}{\partial y}=0$ and $\frac{\partial(\Pi_{K}\tilde{\boldsymbol{\sigma}})_{i}}{\partial x_{j}}=0$ for $i\neq j$. On the other hand, $\frac{\partial(\Pi_{K}\tilde{\boldsymbol{\sigma}})_{i}}{\partial x_{i}}=\frac{\mbox{div\,}_{x}\Pi_{K}\tilde{\boldsymbol{\sigma}}}{n}$ and $\mbox{div\,}_{x}\Pi_{K}\tilde{\boldsymbol{\sigma}}=\frac{1}{|P|}\int_{P}\mbox{div\,}_{x}\tilde{\boldsymbol{\sigma}}$, and so, a simple argument using the Cauchy-Schwarz inequality yields,

[TABLE]

and puting all together we obtain (4.8).

The proof of (4.9) is analogous using now that $\sigma_{n+1}-\Pi_{n+1}\sigma_{n+1}$ has vanishing mean value on the face $K\times\{y_{0}\}$, applying (4.7), and using that $\nabla_{x}(\Pi_{n+1}\sigma_{n+1})=0$ and $\frac{\partial}{\partial y}(\Pi_{n+1}\sigma_{n+1})=\frac{1}{|P|}\int_{P}\frac{\partial\sigma_{n+1}}{\partial y}$. ∎

5. Fractional Laplacian

As an interesting application of the general results for degenerate problems we consider the fractional Laplace equation. Given $\Omega\subset{\mathbb{R}}^{n}$ and $f\in L^{2}(\Omega)$ we want to solve

[TABLE]

for $0<s<1$.

Caffarelli and Silvestre have shown that the solution of this problem can be obtained as $v(x)=u(x,0)$ where $u(x,y)$ is the solution of a degenerate elliptic problem in a cylindrical domain in $n+1$ variables, namely,

[TABLE]

with $d_{s}=2^{1-2s}\frac{\Gamma(1-s)}{\Gamma(s)}$ and $\alpha=1-2s$. To solve this equation numerically one has to approximate the domain $\mathcal{C}$ by a bounded one. With this goal we consider a problem analogous to (5.2) with $\mathcal{C}$ replaced by $\mathcal{C}_{L}=\Omega\times(0,L)$ and adding a homogeneous Dirichlet boundary condition on the upper boundary of $\mathcal{C}_{L}$, namely, we look for $u_{L}$ such that,

[TABLE]

We will use several results proved in [24], therefore, we recall some notation used in that paper. For $0<s<1$, we denote $H^{s}(\Omega)$ the fractional Sobolev space of order $s$. We define for $s\neq\frac{1}{2}$, $\mathbb{H}^{s}(\Omega):=H_{0}^{s}(\Omega)$, the closure of $C_{0}^{\infty}(\Omega)$ in $H^{s}(\Omega)$ and $\mathbb{H}^{1/2}(\Omega):=H_{00}^{1/2}(\Omega)$, the interpolation space $[H_{0}^{1}(\Omega),L^{2}(\Omega)]_{1/2}$ obtained by the K-method (for details see [21]). $\mathbb{H}^{-s}(\Omega)$ denotes the dual space of $\mathbb{H}^{s}(\Omega)$ for $s\in(0,1)$.

For our error estimates we will need some a priori bounds for the derivatives of the exact solution.

In [24] the following a priori estimates for the solution of problem (5.2) were proved,

[TABLE]

and, for $1\leq i,j\leq n$,

[TABLE]

We will use the following estimate: For $\gamma>-1$ and $v\in L^{1}(\mathcal{C}_{L})\cap L^{2}_{y^{\gamma}}(\mathcal{C}_{L})$ such that $\int_{\mathcal{C}_{L}}v=0$, there exists a constant $C$ independent of $L$ such that,

[TABLE]

This estimate can be proved using the arguments introduced in [11]. Details of the proof are given in [13, Lemma 2.2] for a square domain but the arguments apply to more general domains, in particular to the cylindrical ones considered here. That the constant C does not depend on $L$ follows from the case $L=1$ combined with a standard scaling argument.

Lemma 5.1**.**

Let $u$ be the solution of (5.2) and ${\boldsymbol{\sigma}}=(\sigma_{1},\cdots,\sigma_{n+1})=-y^{\alpha}\nabla u$. Then, for $1\leq j\leq n$ and $1\leq i\leq n+1$,

[TABLE]

and for $1\leq i\leq n$ and $\beta>1-\alpha$,

[TABLE]

Proof.

The bound for the first term in (5.7) follows immediately from (5.5). To estimate the second term observe that, from (5.2),

[TABLE]

and use (5.5).

For $1\leq i\leq n$ we have

[TABLE]

To bound the second term we use again (5.5). For the first one we observe that $\int_{\mathcal{C}_{L}}\frac{\partial u}{\partial x_{i}}=0$ because $u$ vanishes on $\partial\Omega\times(0,\infty)$, and therefore, since $\beta>1-\alpha$ we can use (5.6) with $\gamma=\alpha-2+\beta$ to obtain

[TABLE]

where we have used (5.5) for the last inequality.∎

Our goal is to approximate $u$ and ${\boldsymbol{\sigma}}=-y^{\alpha}\nabla u$ given by (5.2). Since the problem is posed in the unbounded domain $\mathcal{C}$ we need to replace it by $\mathcal{C}_{L}$ where $L$ will be chosen in terms of the mesh parameter $h$ in such a way that $L\to\infty$ when $h\to 0$.

It was shown in [24, Theorem 3.5] that for $f\in\mathbb{H}^{-s}(\Omega)$ and $L\geq 1$, if $u_{L}(x,y)$ is extended by zero for $y>L$, there exists a constant $C$ such that

[TABLE]

where $\lambda_{1}>0$ is the first eigenvalue of the Laplacian with Dirichlet boundary conditions in $\Omega$.

Moreover, using the Poincaré inequality

[TABLE]

which follows easily applying the standard Poincaré inequality in $\Omega$ for each $y$, multiplying by the weight, and integrating in $y$, we also have

[TABLE]

Now we consider the mixed finite element approximation of (5.3). We will apply the results of the previous sections for ${\mathcal{D}}=\mathcal{C}_{L}$ and $\Gamma_{N}=\Omega\times\{0\}$. However, since we want error estimates in terms of ${\boldsymbol{\sigma}}$ instead of ${\boldsymbol{\sigma}}_{L}$, to take advantage of the known a priori estimates, we need to introduce some minor modifications in the error analysis.

Given a family of meshes ${\mathcal{T}}_{h}$ made by prismatic elements as those considered in the last part of Section 4 and the associated spaces ${\boldsymbol{S}}_{h}$ and $V_{h}$ defined as in (3.1) and (3.2), the approximate solutions $u_{L,h}\in V_{h}$ and ${\boldsymbol{\sigma}}_{L,h}\in{\boldsymbol{S}}_{h}$ are given by,

[TABLE]

for every face $F$ contained in $\Omega$, and

[TABLE]

where ${\boldsymbol{S}}_{h,N}:={\boldsymbol{S}}_{h}\cap H_{\Gamma_{N}}(\mbox{div\,},{\mathcal{D}})$.

Theorem 5.2**.**

Let $u$ and $u_{L}$ be the solutions of (5.2) and (5.3) respectively, ${\boldsymbol{\sigma}}=-y^{\alpha}\nabla u$ and ${\boldsymbol{\sigma}}_{L}=-y^{\alpha}\nabla u_{L}$. If $u_{L,h}$ and ${\boldsymbol{\sigma}}_{L,h}$ are the approximate solutions given by (5.13), then

[TABLE]

and

[TABLE]

Proof.

Observing that $\Pi_{h}{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{L,h}\in{\boldsymbol{S}}_{h,N}$ and $\mbox{div\,}(\Pi_{h}{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{L,h})=0$ and proceeding as in the proof of Lemma 2.1 we obtain,

[TABLE]

Then,

[TABLE]

and therefore,

[TABLE]

which combined with a triangular inequality yields (5.14).

On the other hand, for our domain $\mathcal{C}_{L}$ the inequality from Lemma 2.3 can be written as

[TABLE]

where the constant $C$ is independent of $L$. Indeed, this follows from the proof of that lemma once we know that the constant in Lemma 2.2 is proportional to $L$, which follows from the case $L=1$ and a scaling argument.

To bound the second term in the right hand side of (5.17) we use (5.16), while for the first one we have

[TABLE]

where in the last inequality we have used the version for prisms of (4.2). To conclude the proof we observe that

[TABLE]

and, therefore, from the Poincaré inequality (5.10) we obtain

[TABLE]

∎

Next we are going to show that introducing appropriate meshes, graded in the $y$-direction, we obtain almost optimal order of convergence with respect to the number of nodes, i. e., the same order than that valid for problems with smooth solutions using uniform meshes, up to a logarithmic factor.

Given a mesh-size $h>0$, to define ${\mathcal{T}}_{h}$ we start with a quasi-uniform triangulation of $\Omega$ made of simplices of diameter less than or equal to $h$. Then, for $L\geq 1$ to be chosen below in terms of $h$, we introduce a partition of $[0,L]$ given by

[TABLE]

where $N\sim 1/h$ (we take $N=1/h$ if it is an integer or some approximation of it if not), and $\beta\in(1-\alpha,2)$ to be chosen (in the numerical experiments we have taken $\beta$ as the midpoint of this interval). Finally, the partition ${\mathcal{T}}_{h}$ of $\mathcal{C}_{L}$ is formed by the prismatic elements $P=K\times[y_{j},y_{j+1}]$, where $K$ are the elements in the partition of $\Omega$.

It follows from this definition that, for $j\geq 1$,

[TABLE]

indeed, by the mean value theorem and using that $h\sim 1/N$ we have

[TABLE]

Using the notation introduced for prismatic elements in the previous section, the Raviart-Thomas interpolation is given by $\Pi_{h}{\boldsymbol{\sigma}}=(\tilde{\Pi}_{h}\tilde{\boldsymbol{\sigma}},\Pi_{h,n+1}\sigma_{n+1})$ where $\tilde{\Pi}_{h}$ and $\Pi_{h,n+1}$ are given locally by $\Pi_{K}$ and $\Pi_{n+1}$ respectively. We recall that, since $-1<\alpha<1$, $y^{\alpha}$ and $y^{-\alpha}$ belong to $A_{2}^{s}$.

Theorem 5.3**.**

For some $\beta\in(1-\alpha,2)$, consider the family of meshes ${\mathcal{T}}_{h}$ defined above. Let $u$ be the solution of (5.2), ${\boldsymbol{\sigma}}=-y^{\alpha}\nabla u$, and $(u_{L,h},{\boldsymbol{\sigma}}_{L,h})$ be the approximation given by (5.12) and (5.13). Then, if $L=C_{1}|\log h|$ with $C_{1}\geq 4/\sqrt{\lambda_{1}}$, we have

[TABLE]

and

[TABLE]

where the constant $C$ depends on $\Omega$, $\alpha$, and $\beta$.

Proof.

From (5.14) and (5.9) we have

[TABLE]

Applying (4.8) for the elements of the form $P=K\times[0,y_{1}]$ and summing over all of them we obtain,

[TABLE]

But,

[TABLE]

where in the last inequality we have used the definition of $y_{1}$. Then,

[TABLE]

Analogously, applying now (4.9), we have

[TABLE]

and therefore, using again the definition of $y_{1}$, we obtain

[TABLE]

Consequently, combining the estimates above, we conclude

[TABLE]

Aplying now (4.8)and (4.9) for the elements of the form $P=K\times[y_{j},y_{j+1}]$, for each $j\geq 1$, and summing over these elements we obtain

[TABLE]

and using (5.19),

[TABLE]

and then, observing that

[TABLE]

summing over $j$, and combining this with (5.23), we obtain

[TABLE]

where, here and in what follows, the constant $C$ depends on $C_{\beta}$.

Applying now Lemma 5.1 and the bound (5.24) it follows from (5.22) that

[TABLE]

From the hypothesis on $C_{1}$ we have $e^{-\sqrt{\lambda_{1}}L/4}\leq h$ and, therefore, (5.20) is proved.

In view of (5.15), to finish the proof of (5.21) it is enough to show that

[TABLE]

Using (4.6) for elements of the form $K\times[0,y_{1}]$ we obtain

[TABLE]

because $2/(2-\beta)\geq 1$ and $L\geq 1$.

On the other hand, (4.6) and (5.19) yields

[TABLE]

and, therefore, taking into account (5.4), (5.25) is proved. ∎

Now we give some numerical examples showing the asymptotic behavior of the error proved in Theorem 5.3. We solve Problem (5.2) with $\Omega=(0,1)\times(0,1)$ and

[TABLE]

Recall that $0<s<1$ and $\alpha=1-2s$. In this case, assuming that $s\neq 1/2$ (i. e. $\alpha\neq 0$), the solution is given by

[TABLE]

where $K_{s}$ is a modified Bessel function of the second kind (see [24]).

We use prismatic elements given by a uniform mesh of triangles in $\Omega$ and the refinement given by (5.18) in the $y$-direction. Observe that for these meshes $h\sim(DOF)^{-1/3}$ where $DOF$ denotes the degrees of freedom. Moreover, we choose $L$ as in Theorem 5.3 with $C_{1}=1$, i. e., $L=|\log h|$.

The next graphics show the order of the errors $\|{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{L,h}\|_{L_{y^{-\alpha}}^{2}(\mathcal{C}_{L})}$ and $\|u-u_{L,h}\|_{L_{y^{\alpha}}^{2}(\mathcal{C}_{L})}$ for several values of $\alpha$.

Finally, to solve (5.1), we need to approximate $u(x,0)$ where $u$ is the solution of (5.2). We will use the approximations $u_{L,h}$ and ${\boldsymbol{\sigma}}_{L,h}$ obtained above.

Since $u_{L,h}$ is only an approximation in the $L^{2}$-norm, one cannot expect that its restriction to $y=0$ be a good approximation of $u(x,0)$. In order to obtain a better approximation we will make a local correction of $u_{L,h}$ using also the computed ${\boldsymbol{\sigma}}_{L,h}$. This correction corresponds to a first order Taylor expansion, indeed, the formula that we are going to prove in the next lemma is motivated by

[TABLE]

We will prove that in this way we obtain an approximation in $L^{2}(\Omega)$ of at least the same order than the mixed finite element approximation of (5.2).

Given $x\in\Omega$ and $0<j<N$ we introduce the jumps

[TABLE]

If $x$ is not in the interior of an element $K$ in the partition of $\Omega$ we choose arbitrary an element containing it to evaluate $u_{L,h}$ (this is irrelevant because afterwards we are going to integrate in $x$).

We will use the standard piecewise linear basis functions, namely, for $1\leq j\leq N-1$,

[TABLE]

[TABLE]

and

[TABLE]

Lemma 5.4**.**

For any $x\in\Omega$ we have

[TABLE]

Proof.

Since $u_{L,h}$ is piecewise constant one can see that

[TABLE]

Let $K$ be the element containing $x$. For $1\leq j\leq N-1$, taking the function $({\boldsymbol{0}},\tau_{j})$ supported in $K\times[y_{j-1},y_{j+1}]$ as test function in (2.7), we have

[TABLE]

and, since $\sigma_{L,h,n+1}(x,y)$ is independent of $x$ for $x\in K$, we obtain

[TABLE]

Analogously, using now $(0,\tau_{N})$ yields

[TABLE]

Therefore, replacing in (5.27) we have

[TABLE]

which immediately gives (5.26) because $\sum_{j=0}^{N}\tau_{j}\equiv 1$. ∎

To approximate the solution of (5.1) given by $v(x)=u(x,0)$ we introduce

[TABLE]

We also define $v_{L}(x)=u_{L}(x,0)$.

Lemma 5.5**.**

[TABLE]

Proof.

Since $u_{L}(x,L)=0$ and, recalling that $\frac{\partial u_{L}}{\partial y}=-y^{-\alpha}{\boldsymbol{\sigma}}_{n+1}$, we have

[TABLE]

Therefore, using (5.26) and the definition of $v_{L,h}$, we obtain

[TABLE]

and, applying the Schwarz inequality,

[TABLE]

and integrating now in $x$ we conclude the proof. ∎

We can now prove the error estimate for the approximation of the solution of the Fractional Laplacian.

Theorem 5.6**.**

Under the hypotheses of Theorem 5.3 we have

[TABLE]

where the constant is as in Theorem 5.3 an depends also on $\alpha$.

Proof.

From Lemma 5.5 and, recalling that $L=C_{1}|\log h|$, we have

[TABLE]

where the constant depends on $\alpha$. Combining this estimate with (5.20) we obtain

[TABLE]

It remains to estimate $v-v_{L}$. But, from the trace theorem given in [24, Proposition 2.5] combined with (5.11)

[TABLE]

and, from the definition of $L$ and $C_{1}$, we obtain

[TABLE]

which combined with (5.28) concludes the proof. ∎

The next graphics show the order of the error $\|v-v_{L,h}\|_{L^{2}(\Omega)}$ for Problem (5.1) with

[TABLE]

which has as exact solution

[TABLE]

Remark 5.1**.**

The order of the error for the approximation of $v$ in the $L^{2}$-norm is probably not the optimal possible. Indeed, with a more complicated postprocessing one could approximate the solution $u$ of Problem (5.2) with order almost $O(h)$ in $H^{1}_{y^{\alpha}}(\mathcal{C})$ and, by the trace theorem $\|v\|_{\mathbb{H}^{s}(\Omega)}\leq C\|v\|_{H^{1}_{y^{\alpha}}(\mathcal{C})}$ proved in [24, Proposition 2.5], one would have the same order for the approximation of $v$ in the $\mathbb{H}^{s}$-norm. Therefore, it is reasonable to expect a higher order in $L^{2}$. Let us mention also that, as far as we know, such a higher order error estimate has not been proved either for the standard method analyzed in [24]. This problem requires a different analysis and will be the object of our further research.

Acknowledgement: We thank Enrique Otárola for helpful comments.