$hp$-Finite Elements for Fractional Diffusion

TL;DR

This paper introduces an $hp$-finite element numerical scheme for efficiently solving boundary value problems involving the spectral fractional Laplacian, reducing computational complexity and improving convergence.

Contribution

It presents a novel $hp$-finite element discretization approach on a truncated semi-infinite cylinder for spectral fractional Laplacian problems, enhancing efficiency and accuracy.

Findings

Significant reduction in degrees of freedom needed for computation

Slightly improved convergence properties over linear finite element methods

Validated effectiveness through numerical experiments

Abstract

The purpose of this work is to introduce and analyze a numerical scheme to efficiently solve boundary value problems involving the spectral fractional Laplacian. The approach is based on a reformulation of the problem posed on a semi-infinite cylinder in one more spatial dimension. After a suitable truncation of this cylinder, the resulting problem is discretized with linear finite elements in the original domain and with $hp$-finite elements in the extended direction. The proposed approach yields a drastic reduction of the computational complexity in terms of degrees of freedom and even has slightly improved convergence properties compared to a discretization using linear finite elements for both the original domain and the extended direction. The performance of the method is illustrated by numerical experiments.

Equations

$$(-\Delta {)}^s\mathfrak{u}$$

$$\mathfrak{u}$$

$$\mathrm{div}(y^{1-2s}\nabla u)$$

$$u$$

$${\partial }_{{\nu }^{1-2s}}u$$

$${\int }_{\Omega }\nabla {\phi }_k\cdot \nabla vdx={\lambda }_k{\int }_{\Omega }{\phi }_{k}vdx\forall v\in H_0^1(\Omega ).$$

$$\text{mathds}{H}^s(\Omega )=\left\{\mathfrak{v}\in L^2(\Omega )\right\}\Vert \mathfrak{v}{\Vert }_{\text{mathds}{H}^s(\Omega )}^2=\sum _{k=1}^{\infty }{\lambda }_k^s{\mathfrak{v}}_k^2<\infty \text{ with }{\mathfrak{v}}_k={\int }_{\Omega }\mathfrak{v}{\phi }_{k}dx.$$

$$(-\Delta {)}^s\mathfrak{u}=\sum _{k=1}^{\infty }{\lambda }_k^s{\mathfrak{u}}_k{\phi }_k\in \text{mathds}{H}^{-s}(\Omega )\text{with}{\mathfrak{u}}_k={\int }_{\Omega }\mathfrak{u}{\phi }_{k}dx.$$

$$\sum _{k=1}^{\infty }{\lambda }_k^s{\mathfrak{u}}_k{\mathfrak{v}}_k={\int }_{\Omega }\mathfrak{fv}dx\forall \mathfrak{v}\in \text{mathds}{H}^s(\Omega )\text{with}{\mathfrak{v}}_k={\int }_{\Omega }\mathfrak{v}{\phi }_{k}dx.$$

$$\mathfrak{u}=\sum _{k=1}^{\infty }{\mathfrak{u}}_k{\phi }_k\text{with}{\mathfrak{u}}_k={\lambda }_k^{-s}{\mathfrak{f}}_k\text{and}{\mathfrak{f}}_k={\int }_{\Omega }\mathfrak{f}{\phi }_{k}dx.$$

$$\text{mathds}{H}^s(\Omega )=\{\begin{array}{ll}{H}^s(\Omega ),& \text{for }0<s<\frac{1}{2},\\ H_{00}^{\frac{1}{2}}(\Omega ),& \text{for }s=\frac{1}{2},\\ H_0^s(\Omega ),& \text{for }\frac{1}{2}<s<1.\end{array}$$

$${\mathring{H}}_L^1(C,y^{\alpha })=\left\{v\in H^1(y^{\alpha },C)\right\}v=0\text{ on }{\partial }_{L}C.$$

$${\mathring{H}}_L^1(C_Y,y^{\alpha })=\left\{v\in H^1(y^{\alpha },C_Y)\right\}v=0\text{ on }{\partial }_L{C}_Y\cup (\Omega \times \{Y\}).$$

$${\mathrm{tr}}_{\Omega }{\mathring{H}}_L^1(C,y^{\alpha })=\text{mathds}{H}^s(\Omega )\text{and}\Vert {\mathrm{tr}}_{\Omega }v{\Vert }_{\text{mathds}{H}^s(\Omega )}\le c\Vert v{\Vert }_{{\mathring{H}}_L^1(C,y^{\alpha })}.$$

$${\int }_C{y}^{\alpha }\nabla u\cdot \nabla vd(x,y)=d_s⟨\mathfrak{f},{\mathrm{tr}}_{\Omega }v{⟩}_{\text{mathds}{H}^{-s}(\Omega ),\text{mathds}{H}^s(\Omega )}\forall v\in {\mathring{H}}_L^1(C,y^{\alpha })$$

$$\mathrm{div}(y^{\alpha }\nabla u)$$

$$u$$

$${\partial }_{{\nu }^{\alpha }}u$$

$$\Vert y^{\frac{\alpha }{2}}{\nabla }_x^{2}u{\Vert }_{L^2(C)}^2+\Vert y^{\frac{\alpha }{2}}{\partial }_y{\nabla }_{x}u{\Vert }_{L^2(C)}^2\le c\Vert \mathfrak{f}{\Vert }_{\text{mathds}{H}^{1-s}(\Omega )}^2.$$

$$\Vert y^{\frac{\alpha }{2}}{\Delta }_{x}u{\Vert }_{L^2(C)}^2+\Vert y^{\frac{\alpha }{2}}{\partial }_y{\nabla }_{x}u{\Vert }_{L^2(C)}^2=d_s\Vert \mathfrak{f}{\Vert }_{\text{mathds}{H}^{1-s}(\Omega )}^2$$

$${\psi }_s(z)=c_s{z}^s{K}_s(z)\text{with}{c}_s=\frac{2^{1-s}}{\Gamma (s)}.$$

$$u(x,y)=\sum _{k=1}^{\infty }{\mathfrak{u}}_k{\phi }_k(x){\psi }_{s,k}(y),$$

$$c_{\frac{1}{2}}{z}^{\frac{1}{2}}{K}_{\frac{1}{2}}(z)=e^{-z}={\psi }_{\frac{1}{2}}(z).$$

$$|y^n{\psi }_{s,k}^{(n)}(y)|\le c{8}^{n}n!{\lambda }_k^{s-\frac{r}{2}}{y}^{2s-r}.$$

$${\psi }_{s,k}^{(n)}(y)=\frac{d^n}{d{y}^n}\psi (\sqrt{{\lambda }_k}y)=(\sqrt{{\lambda }_k}{)}^n{\psi }^{(n)}(\sqrt{{\lambda }_k}y).$$

$$|y^n{\psi }_{s,k}^{(n)}(y)|=|(\sqrt{{\lambda }_k}y{)}^n{\psi }^{(n)}(\sqrt{{\lambda }_k}y)|\le c{8}^{n}n!(\sqrt{{\lambda }_k}y{)}^{2s-r}=c{8}^{n}n!{\lambda }_k^{s-\frac{r}{2}}{y}^{2s-r}.\text{qed}$$

$$|y^{r_1}{\psi }_{s,k}(y)|\le c_1{\lambda }_k^{-\frac{r_1}{2}}{e}^{-\frac{\sqrt{{\lambda }_k}}{2}y}\text{and}|y^{r_2}{\psi }_{s,k}^{\prime }(y)|\le c_2{\lambda }_k^{-\frac{r_2}{2}}{e}^{-\frac{\sqrt{{\lambda }_k}}{2}y}.$$

$$(\sqrt{{\lambda }_k}{)}^{s_0}{e}^{\sqrt{{\lambda }_k}}{K}_s(\sqrt{{\lambda }_k})\le (\sqrt{{\lambda }_1}{)}^{s_0}{e}^{\sqrt{{\lambda }_1}}{K}_s(\sqrt{{\lambda }_1})$$

$$|y^{r_1}{\psi }_{s,k}(y)|={\lambda }_k^{-\frac{r_1}{2}}|(\sqrt{{\lambda }_k}y{)}^{r_1}{\psi }_s(\sqrt{{\lambda }_k}y)|.$$

$$\Vert y^{\frac{\alpha }{2}}\nabla u{\Vert }_{L^2(C\setminus C_Y)}\le c{e}^{-\frac{\sqrt{{\lambda }_1}}{2}Y}\Vert \mathfrak{f}{\Vert }_{\text{mathds}{H}^{-s}(\Omega )}.$$

$$\Vert y^{\frac{\alpha }{2}}\nabla u{\Vert }_{L^2(C\setminus C_Y)}^2$$

$$\displaystyle=\sum_{k=1}^{\infty}\mathfrak{u}_{k}^{2}\int_{Y}^{\infty}\bigl{\{}\lambda_{k}\lvert y^{\frac{\alpha}{2}}\psi_{s,k}(y)\rvert^{2}+\lvert y^{\frac{\alpha}{2}}\psi_{s,k}^{\prime}(y)\rvert^{2}\bigr{\}}\,dy$$

$$\displaystyle\leq c\sum_{k=1}^{\infty}\mathfrak{u}_{k}^{2}\int_{Y}^{\infty}\bigl{\{}\lambda_{k}^{s+\frac{1}{2}}e^{-\sqrt{\lambda_{k}}y}+\lambda_{k}^{s-\frac{1}{2}}e^{-\sqrt{\lambda_{k}}y}\bigr{\}}\,dy,$$

$$\lVert y^{\frac{\alpha}{2}}\nabla u\rVert_{L^{2}(C\setminus C_{Y})}^{2}\leq c\sum_{k=1}^{\infty}\mathfrak{u}_{k}^{2}\bigl{\{}\lambda_{k}^{s}e^{-\sqrt{\lambda_{k}}Y}+\lambda_{k}^{s-1}e^{-\sqrt{\lambda_{k}}Y}\bigr{\}}\leq ce^{-\sqrt{\lambda_{1}}Y}\sum_{k=1}^{\infty}\lambda_{k}^{s}\mathfrak{u}_{k}^{2}=ce^{-\sqrt{\lambda_{1}}Y}\lVert\mathfrak{f}\rVert^{2}_{\mathds{H}^{-s}(\Omega)},$$

Full text

\minisec

AMS subject classification 35S15, 65R20, 65N12, 65N30

1 Introduction

In this work, we are concerned with boundary value problems involving the fractional Laplacian, being the prototype of a nonlocal operator. To be more specific: Let $\Omega\subset\mathds{R}^{d}$ for $d\in\set{1,2,3}$ be a bounded, convex, polygonal or polyhedral domain. We are interested in the solution of the boundary value problem

[TABLE]

where $(-\Delta)^{s}$ denotes the spectral fractional Laplacian of order $s\in(0,1)$ defined by the eigenvalues and eigenfunctions of the standard Laplacian, precisely introduced in Section 2. The main purpose of this paper is to introduce and analyze a numerical scheme to efficiently solve problem (1.1).

Our approach is based on the following equivalent reformulation of problem (1.1) posed on the semi-infinite cylinder $\set{(x,y)\in\Omega\times(0,\infty)}\subset\mathds{R}^{d+1}$: Let $u$ be the weak solution of the extended problem

[TABLE]

see Section 2 for more details. Then, the trace $\mathfrak{u}=u(\cdot,0)$ is the solution of the fractional boundary value problem (1.1).

In contrast to the nonlocal problem (1.1), the extended problem is localized. However, a direct application of a finite element method to the extended problem is not feasible because of the semi-infinite domain. As remedy, the exponential decay of $u$ in direction $y$ towards infinity (see Proposition 2.9) can be employed such that a truncation of the semi-infinite cylinder to $\Omega\times(0,Y)$ becomes reasonable. The extended problem posed on the truncated cylinder can be discretized using finite elements. However, due to the degenerate/singular nature of the extended problem, anisotropic meshes are favorable in order to obtain an optimally convergent numerical scheme. Moreover, the height $Y$ of the truncated cylinder needs to be chosen dependent on the mesh parameter to ensure the aforementioned convergence. This approach was already pursued in [27] using a discretization with first degree tensor product finite elements on graded meshes in the extended direction, see also [14, 15, 26, 28] for related results. If $h_{\Omega}$ denotes the mesh parameter and $\mathcal{N}_{\Omega}$ the number of degrees of freedom in $\Omega$ then the approach from [27] yields a discretization error of order $O(h_{\Omega}\lvert\ln h_{\Omega}\rvert^{s})$ in the corresponding energy norm associated with (1.1) while solving problems with $O(\mathcal{N}_{\Omega}^{1+1/d})$ degrees of freedom.

In this work, we introduce and analyze a discretization of the truncated problem with linear finite elements in the original domain $\Omega$ and with $hp$-finite elements on a geometric mesh in the extended direction. This drastically reduces the computational complexity to $O(\mathcal{N}_{\Omega}(\ln\mathcal{N}_{\Omega})^{2})$ degrees of freedom and even yields a slightly better convergence rate of order $O(h_{\Omega})$. Especially, when $\mathcal{N}_{\Omega}$ is large, the difference between the factor $(\ln\mathcal{N}_{\Omega})^{2}$ and $\mathcal{N}_{\Omega}^{1/d}$ becomes clearly perceptible. For instance, in our numerical experiments we could reduce the number of degrees of freedom by a factor of about $111$ to obtain an error of less than $9\cdot 10^{-3}$ in the case $s=0.8$, see Section 5 for more details. We also notice that our approach and results are not limited to the spectral fractional Laplacian. They naturally extend by only minor modifications to fractional powers of general second order elliptic operators.

Let us briefly give an overview on other numerical approaches from the literature to solve boundary value problems involving the fractional Laplacian: Due to the spectral definition of the operator, it seems to be natural to compute an approximating, discrete spectral decomposition of the standard Laplacian in order to get an approximation of the solution of (1.1), see [21, 22, 31]. However, this may result in solving a large number of discrete eigenvalue problems. Another approach to determine an approximation to the solution of problem (1.1) is analyzed in [8], see also [6, 7, 9] for related results. In that reference, $(-\Delta)^{-s}$ is represented in terms of Bochner integrals involving $(I-t^{2}\Delta)^{-1}$ for $t\in(0,\infty)$. Subsequently, different quadrature formulas to approximate this integral are analyzed which require multiple evaluations of $(I-t_{i}^{2}\Delta_{h})^{-1}$ with $t_{i}$ being a quadrature point and $-\Delta_{h}$ denoting a finite element discretization to $-\Delta$. Numerical approaches for the integral definition of the fractional Laplacian, which is not equivalent to the spectral definition considered in the present paper, can be found in [2, 3, 4, 5, 10, 16, 18, 19, 20].

This paper is organized as follows: In Section 2, we state the definition of the fractional Laplacian, formulate the extended problem in detail, and introduce the functional framework needed for the subsequent error analysis. Moreover, in this section, we are concerned with several properties of the solution of the extended problem such as a series representation and corresponding regularity results. The discrete, extended problem posed on the truncated cylinder is formulated at the beginning of Section 3. In the extended direction, we distinguish between graded meshes and $h$-FEM, and geometric meshes and $hp$-FEM, see Sections 3.1 and 3.2. The error analysis is given in Section 4. Thereby, in Section 4.1, we mainly recover the results of [27]. The reason for doing this is twofold. First, we are able to slightly improve the mesh grading condition used in [27]. However, the main reason to analyze the $h$-FEM on graded meshes before developing the analysis for the $hp$-method considered in Section 4.2 is, that the techniques we use are almost identical for both cases, but the details are simpler for $h$-FEM. Implementation aspects and numerical experiments, which underline the efficiency of our approach, are presented in Section 5. In the appendix, we collect different results for special functions defined by the modified Bessel functions of second kind. These are especially needed in Section 2 for the discussion of the solution of the extended problem.

Finally, we notice that, in the following, $c$ denotes a generic constant which will always be independent of the mesh parameter $h_{\Omega}$ when we analyze the discretization error.

2 Continuous Problem

Let $-\Delta$ be the $L^{2}(\Omega)$ realization of the Laplacian with homogeneous Dirichlet boundary conditions. It is well-known that $-\Delta$ has a compact resolvent and its eigenvalues form a non-decreasing sequence $0<\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{k}\leq\cdots$ satisfying $\lim_{k\to\infty}\lambda_{k}=\infty$. We denote by $\varphi_{k}$ the orthonormal $L^{2}(\Omega)$ eigenfunctions associated with $\lambda_{k}$ fulfilling

[TABLE]

For any $s\geq 0$, we introduce the fractional order Sobolev space

[TABLE]

Moreover, we denote by $\mathds{H}^{-s}(\Omega)$ the dual space of $\mathds{H}^{s}(\Omega)$. Then, the spectral fractional Laplacian is defined for $s\in(0,1)$ on the space $\mathds{H}^{s}(\Omega)$ as the limit

[TABLE]

Due to the Cauchy-criterion the limit exists for any $\mathfrak{u}\in\mathds{H}^{s}(\Omega)$. Thus, problem (1.1) has to be understood as: Given $\mathfrak{f}\in\mathds{H}^{s}(\Omega)$, find $\mathfrak{u}\in\mathds{H}^{s}(\Omega)$ such that

[TABLE]

Proposition 2.1**.**

For any $\mathfrak{f}\in\mathds{H}^{-s}(\Omega)$, problem (1.1) admits a unique solution $\mathfrak{u}\in\mathds{H}^{s}(\Omega)$ fulfilling $\lVert\mathfrak{u}\rVert_{\mathds{H}^{s}(\Omega)}=\lVert\mathfrak{f}\rVert_{\mathds{H}^{-s}(\Omega)}$. Moreover, there is the series representation

[TABLE]

Proof.

The existence of an unique solution and the equality of the norms is a consequence of the Riesz representation theorem. The series representation of $\mathfrak{u}$ is obtained by testing (2.1) with $\varphi_{m}\in\mathds{H}^{s}(\Omega)$ and using the orthogonality of the eigenfunctions. ∎

Remark 2.2*.*

Due to the definition of the fractional Laplacian and the previous result, we observe that problem (1.1) is already meaningful without additionally imposing the homogeneous Dirichlet boundary conditions since these are already included in the definition of the operator. Moreover, we notice that the regularity of $\mathfrak{u}$ can be described in classical fractional Sobolev spaces as well, since

[TABLE]

For more details we refer to, e.g., [27].

Problem (1.1) can equivalently be posed on a semi-infinite cylinder. In $\mathds{R}^{d}$, this is due to Caffarelli and Silvestre [12]. The restriction to bounded domains $\Omega$ was considered by Stinga and Torrea in [30], see also [11, 13]. This kind of extension is the basis for the computational approaches in the subsequent sections.

In order to state the extended problem, we first introduce the required notation. We denote by $C=\Omega\times(0,\infty)$ the aforementioned semi-infinite cylinder and by $\partial_{L}C=\partial\Omega\times[0,\infty)$ its lateral boundary. We also need to define a truncated cylinder: for $Y>0$, the truncated cylinder is given by $C_{Y}=\Omega\times(0,Y)$ with its lateral boundary $\partial_{L}C_{Y}=\partial\Omega\times[0,Y]$. As $C$ and $C_{Y}$ are objects in $\mathds{R}^{d+1}$, we use $y$ to denote the extended variable, such that a vector $(x,y)\in\mathds{R}^{d+1}$ admits the representation $(x,y)=(x_{1},x_{2},\dots,x_{d},y)$. Similarly, the gradient in $\mathds{R}^{d+1}$ has the representation $\nabla=(\nabla_{x},\partial_{y})=(\partial_{x_{1}},\partial_{x_{2}},\dots,\partial_{x_{d}},\partial_{y})$.

Next we introduce weighted Sobolev spaces with a weight function $y^{\alpha}$ for $\alpha\in(-1,1)$. In this regard, let $D\subset\mathds{R}^{d}\times(0,\infty)$ be an open set, such as $C$ or $C_{Y}$. Then, we define the weighted space $L^{2}(D,y^{\alpha})$ as the space of all measurable functions on $D$ with finite norm $\lVert v\rVert_{L^{2}(D,y^{\alpha})}=\lVert y^{\frac{\alpha}{2}}v\rVert_{L^{2}(D)}$. Similarly, the space $H^{1}(D,y^{\alpha})$ denotes the space of all functions $v\in L^{2}(D,y^{\alpha})$ whose weak derivatives of first order belong to $L^{2}(D,y^{\alpha})$.

To study the extended problems, we introduce the space

[TABLE]

The space $\mathring{H}^{1}_{L}(C_{Y},y^{\alpha})$ is defined analogously, but endowed with zero Dirichlet boundary conditions also on $\Omega\times\set{Y}$:

[TABLE]

For $v\in\mathring{H}^{1}_{L}(C,y^{\alpha})$, we denote by $\operatorname{tr}_{\Omega}v$ the trace of $v$ onto $\Omega\times\set{0}$, i.e., $\operatorname{tr}_{\Omega}v=v(\cdot,0)$.

Proposition 2.3**.**

For $\alpha=1-2s$, it holds

[TABLE]

Proof.

See [11, Proposition 1.8] for $s=\frac{1}{2}$ and [13, Proposition 2.1] for $s\neq\frac{1}{2}$. ∎

Now, we are able to state the extended problem: Given $\mathfrak{f}\in\mathds{H}^{-s}(\Omega)$, find $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ such that

[TABLE]

with $\alpha=1-2s$ and $d_{s}=2^{\alpha}\frac{\Gamma(1-s)}{\Gamma(s)}$. That is, the function $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ is a weak solution of

[TABLE]

where we have set $\partial_{\nu^{\alpha}}u(x,0)=\lim_{y\to 0}y^{\alpha}\partial_{y}u(x,y)$. Note that subsequently, the parameter $\alpha$ will always be equal to $1-2s$ for the considered $s\in(0,1)$.

In the remainder of this section, we discuss several properties of the solution $u$ to (2.2). Due to the following proposition, it is reasonable to determine the solution $u$ of (2.2) in order to get the solution $\mathfrak{u}$ of (1.1).

Proposition 2.4**.**

For $\mathfrak{f}\in\mathds{H}^{-s}(\Omega)$, the extended problem (2.2) admits a unique solution $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$. Furthermore, $\mathfrak{u}=\operatorname{tr}_{\Omega}u\in\mathds{H}^{s}(\Omega)$ solves (1.1).

Proof.

See [13, Lemma 2.1]. ∎

We have the following regularity result.

Proposition 2.5**.**

Let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2) and $\mathfrak{f}\in\mathds{H}^{1-s}(\Omega)$. Then, it holds

[TABLE]

Proof.

[27, Theorem 2.7] yields

[TABLE]

with $d_{s}$ given above. Convexity of $\Omega$ then implies the assertion. ∎

For a series expansion of the solution to (2.2), we introduce the function $\psi_{s}\colon[0,\infty)\to\mathds{R}$ given by

[TABLE]

Here, $K_{s}(z)$ denotes the modified Bessel function of second kind, see, e.g., [1, Section 9.6] and $\Gamma(s)$ denotes the gamma function.

Proposition 2.6**.**

For $\mathfrak{f}\in\mathds{H}^{-s}(\Omega)$, let $\mathfrak{u}\in\mathds{H}^{s}(\Omega)$ be the solution of (1.1) and $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Then, it holds

[TABLE]

where $\mathfrak{u}_{k}=\int_{\Omega}\mathfrak{u}\varphi_{k}\,dx$ and $\psi_{s,k}(y)=\psi_{s}(\sqrt{\lambda_{k}}y)$.

Proof.

For $s\neq\frac{1}{2}$, see [13, Proposition 2.1]. For $s=\frac{1}{2}$, the result was proved with $\psi_{\frac{1}{2}}(z)=e^{-z}$ in [11, Proposition 2.2]. However, by [1, 9.6.23] and the relation $\Gamma(\frac{1}{2})=\sqrt{\pi}$, it holds

[TABLE]

Hence, $\psi_{s}(z)=c_{s}z^{s}K_{s}(z)$ holds also in the case $s=\frac{1}{2}$. ∎

We state an estimate of the derivatives of $\psi_{s,k}$ which we will use later. The result is based on properties of $\psi_{s}$ which are analyzed in the appendix.

Corollary 2.7**.**

Let $r\in[0,1]$. There exists a constant $c>0$ depending only on $s$, such that for any $y>0$ and $n\in\mathds{N}$ it holds

[TABLE]

Proof.

Simple calculations yield

[TABLE]

Consequently, we obtain by means of Lemma A.3

[TABLE]

Next, we state a result about the exponential decay of $\psi_{s,k}$ and its derivative. It is based on corresponding results for $\psi_{s}$ and its derivative, proved in the appendix.

Corollary 2.8**.**

Let $y\geq 1$, $r_{1}\geq\min(s,\frac{1}{2})-s$ and $r_{2}\geq\min(1-s,\frac{1}{2})-s$. Then, there exists a constant $c_{1}$ only depending on $r_{1}$, $s$, and $\lambda_{1}$ and a constant $c_{2}$ only depending on $r_{2}$, $s$, and $\lambda_{1}$ such that

[TABLE]

Proof.

Let $s_{0}=\min(s,\frac{1}{2})$. According to [25, Theorem 5], we get for $z\geq 0$ that $z^{s_{0}}e^{z}K_{s}(z)$ is a decreasing function. Consequently, we obtain

[TABLE]

since the sequence $(\lambda_{k})_{k\in\mathds{N}}$ is non-decreasing. Moreover, due to the definition of $\psi_{s,k}$, we deduce

[TABLE]

Thus, by setting $a=\sqrt{\lambda_{k}}$ and $z=\sqrt{\lambda_{k}}y\geq a$ in Lemma A.4 (a), we obtain the validity of the first inequality of the assertion by means of (2.4). The second inequality can be deduced in the same manner employing Lemma A.4 (b). ∎

As already mentioned, for computational reasons, the semi-infinite cylinder $C$ will be truncated to $C_{Y}=\Omega\times(0,Y)$ for some $Y>0$ later on. Because of this, the behavior of $\nabla u$ for $y\to\infty$ will play a role. It can be estimated as follows:

Proposition 2.9**.**

For $\mathfrak{f}\in\mathds{H}^{-s}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Then, there exists a constant $c>0$ such that for every $Y\geq 1$, it holds

[TABLE]

Proof.

The result can be found in [27, Proposition 3.1]. For the sake of completeness, we state a (slightly different) proof here. According to Proposition 2.6, we obtain by using the definition of the eigenfunctions $\varphi_{k}$ and its orthogonality

[TABLE]

where we employed Corollary 2.8 in the last step with $r_{1}=r_{2}=\frac{\alpha}{2}=\frac{1}{2}-s$. Calculating the integral and using that $(\lambda_{k})_{k\in\mathds{N}}$ is a non-decreasing sequence, yields

[TABLE]

where we used Proposition 2.1 in the last step. ∎

3 Discretization

Let $\mathcal{T}_{\Omega}$ be a conforming and quasi-uniform triangulation of $\Omega$ which is admissible in the sense of Ciarlet. For each $\mathcal{T}_{\Omega}=\set{K}$, let $K\subset\mathds{R}^{d}$ be an element that is isoparametrically equivalent either to the unit cube or to the unit simplex in $\mathds{R}^{d}$. We introduce the global mesh parameter $h_{\Omega}$ with respect to the triangulation of $\Omega$ by $h_{\Omega}=\max_{K\in\mathcal{T}_{\Omega}}\operatorname{diam}K$. We always assume that $h_{\Omega}\leq\frac{1}{2}$. On $\mathcal{T}_{\Omega}$, we define a finite element space $V_{h}$ as

[TABLE]

In case that $K$ is a simplex then $\mathds{P}_{1}(K)=\mathcal{P}_{1}(K)$, the set of polynomials on the element $K$ of degree at most $1$. If $K$ is a cube then $\mathds{P}_{1}(K)$ equals $\mathcal{Q}_{1}(K)$, the set of polynomials on $K$ of degree at most 1 in each variable. The number of degrees of freedom in $V_{h}$ is denoted by $\mathcal{N}_{\Omega}$. It holds $\mathcal{N}_{\Omega}=O(h_{\Omega}^{-d})$.

Furthermore, let $\mathcal{I}_{Y}=\set{I_{m}}$ be a triangulation of the interval $(0,Y)$ in the sense that $[0,Y]=\bigcup_{m=1}^{M}I_{m}$ with $I_{m}=[y_{m-1},y_{m}]$ and $M\in\mathds{N}$ exactly specified below in the Sections 3.1 and 3.2. Moreover, let $h_{m}=\lvert I_{m}\rvert$. Next, we introduce a polynomial degree vector $p=(p_{1},p_{2},\dots,p_{M})\in\mathds{N}^{M}$ which will assign to each element $I_{m}\in\mathcal{I}_{Y}$ a maximal polynomial degree $p_{m}$. It will be exactly specified in the Sections 3.1 and 3.2. On $\mathcal{I}_{Y}$, we define the finite element space $V_{M}$ as

[TABLE]

where $\mathcal{P}_{p_{m}}(I_{m})$ denotes the space of polynomials up to degree $p_{m}$ on $I_{m}$.

Now, the triangulations $\mathcal{T}_{\Omega,Y}$ of the cylinder $C_{Y}$ are constructed as tensor product triangulations by means of $\mathcal{T}_{\Omega}$ and $\mathcal{I}_{Y}$, i.e., $\mathcal{T}_{\Omega,Y}=\set{T}$ with $T=K\times I_{m}$ for $K\in\mathcal{T}_{\Omega}$ and $I_{m}\in\mathcal{I}_{Y}$. By means of the previous considerations, we define the finite element space $V_{h,M}$ posed on the tensor product mesh $\mathcal{T}_{\Omega,Y}$ by

[TABLE]

Note, that each function $v_{h}\in V_{h,M}$ vanishes on the lateral boundary of $C_{Y}$ and on its top. As a consequence, the extension by zero of $v_{h}$ to the semi-infinite cylinder $C$ belongs to $\mathring{H}^{1}_{L}(C,y^{\alpha})$. Without further mention, we consider this type of extension for each $v_{h}\in V_{h,M}$ whenever needed.

With the just introduced notation, we define approximations to the solution $u$ of (2.2) as follows: Find $u_{h}\in V_{h,M}$ satisfying

[TABLE]

where we recall that $\alpha=1-2s$ and $d_{s}=2^{\alpha}\tfrac{\Gamma(1-s)}{\Gamma(s)}$. Note that $\operatorname{tr}_{\Omega}u_{h}=u_{h}(\cdot,0)$ will be used as an approximation of $\mathfrak{u}$.

We distinguish two possible types of discretization in the artificial $y$ direction, which will be defined in the following two sections.

3.1 Graded meshes and $h$-FEM

In this section, let $M\in\mathds{N}$ to be determined later. We set

[TABLE]

where $\mu\in(0,1]$ represents the grading parameter in direction $y$. Hence, for $\mu<1$ the triangulation $\mathcal{T}_{\Omega,Y}$ is anisotropic. Moreover, due to the choice of the polynomial degree vector $p$, the discrete space $V_{h,M}$ consists of globally continuous and piecewise multilinear functions on $C_{Y}$.

We start with a result regarding the diameter $h_{m}$ of the elements $I_{m}$.

Lemma 3.1**.**

It holds

[TABLE]

Proof.

The first equality is obvious due to the definition of $y_{1}$. For the second, we observe that there exists a $m_{*}\in(m-1,m)$ such that

[TABLE]

due to the mean value theorem. Since $\frac{1}{2}m\leq m-1\leq m$ for $m\geq 2$, the result follows. ∎

Since we consider a discretization with linear polynomials in direction $y$, the number of degrees of freedom $\mathcal{N}_{Y}$ in $V_{M}$ is proportional to $M$, i.e., it holds $\mathcal{N}_{Y}=O(M)$.

3.2 Geometric meshes and $hp$-FEM

For the second possibility of discretization presented here, the mesh in direction $y$ is chosen as a geometric mesh. That is, for chosen $\sigma\in(0,1)$ and $M\in\mathds{N}$ to be determined later, the nodes $y_{0},y_{1},\dots,y_{M}$ are given by

[TABLE]

Further, we define $p=(p_{1},p_{2},\dots,p_{M})\in\mathds{N}^{M}$ to be a linear degree vector with slope $\beta>0$. That is, there exists a constant $c\geq 1$ such that for all $m=1,2,\dots,M$ the following relation is fulfilled:

[TABLE]

We start with collecting some basic results for the discretization considered in this subsection.

Lemma 3.2**.**

For a geometric mesh given by (3.2), there holds for $h_{m}=\lvert I_{m}\rvert$

[TABLE]

Proof.

The first equation of the assertion is obvious due to the definition of the geometric mesh. For $m=2,3,\dots,M$, we obtain

[TABLE]

As a consequence, we get

[TABLE]

Lemma 3.3**.**

For a geometric mesh given by (3.2) and a linear degree vector $p\in\mathds{N}^{M}$ in the sense of (3.3), it holds $p_{1}=1$ and there is a constant $c\geq 1$ such that

[TABLE]

for all $m=2,3,\dots,M$.

Proof.

Since $\ln\frac{h_{1}}{h_{1}}=0$, there obviously holds $p_{1}=1$. According to Lemma 3.2, we get for $m=2,3\dots,M-1$

[TABLE]

which shows the second assertion. ∎

Lemma 3.4**.**

Let $\mathcal{I}_{Y}$ be a geometric triangulation of $(0,Y)$ given by (3.2) and let $p\in\mathds{N}^{M}$ be a linear degree vector as defined in (3.3). Then, for the number of degrees of freedom $\mathcal{N}_{Y}$ in $V_{M}$, it holds

[TABLE]

Proof.

On each interval $I_{m}$ we locally have $p_{m}+1$ degrees of freedom. Moreover, two neighboring intervals always share one degree of freedom. Thus, we have

[TABLE]

which proves the assertion having in mind Lemma 3.3. ∎

4 Error Estimates

We start with providing a general error estimate between the solutions $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ of (2.2) and $u_{h}\in V_{h,M}$ of (3.1) in weighted norms.

Lemma 4.1**.**

Let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2) and $u_{h}\in V_{h,M}$ be the solution of (3.1). Then, it holds

[TABLE]

Proof.

Testing (2.2) with $v_{h}\in V_{h,M}$ and then subtracting (3.1) from this equation yields

[TABLE]

Based on this, we deduce for all $v_{h}\in V_{h,M}$

[TABLE]

Dividing by $\lVert y^{\frac{\alpha}{2}}\nabla(u-u_{h})\rVert_{L^{2}(C)}$ ends the proof. ∎

Whereas the last term in the estimate of Lemma 4.1 can be treated by means of Proposition 2.9, we have to estimate the first term. To this end, we introduce the following approximation operators separately for the $x$ and $y$ variables:

By $\pi_{x}\colon L^{2}(\Omega)\to V_{h}$, we denote the $L^{2}$ projection with respect to the $x$ variable on $\Omega$. For the interpolation with respect to the $y$ direction, we consider each interval $I_{m}$ separately. For fixed $m\in\set{1,2,\dots,M}$, let $q\in\mathds{N}$ and $y_{m-1}=x_{0},x_{1},\ldots,x_{q}=y_{m}$ be the Gauss-Lobatto points in $I_{m}$ and let $l_{i,q}$ denote the corresponding Lagrange polynomials of order $q$. Then, we define the Gauss-Lobatto interpolant $i_{q}\colon C^{0}(\overline{I_{m}})\to\mathcal{P}_{q}(I_{m})$ by

[TABLE]

Further, we define an interpolant $\tilde{i}_{q}$ which admits $(\tilde{i}_{q}\xi)(y_{m})=0$ given by

[TABLE]

Based on this, we define the interpolation $i_{y}^{p}\colon C^{0}((0,Y])\to V_{M}$ for a linear degree vector $p=(p_{1},p_{2},\dots,p_{M})\in\mathds{N}^{M}$ by

[TABLE]

In particular, it holds that $i_{y}^{p}\xi$ is constant on $I_{1}$ and $(i_{y}^{p}\xi)(Y)=0$. For any function $v\in\mathring{H}^{1}_{L}(C,y^{\alpha})$, we set

[TABLE]

Moreover, for the solution $u$ of (2.2), the application of $i_{y}^{p}$ is defined as

[TABLE]

which is well-defined because $\psi_{s,k}\in C^{0}(\mathds{R}_{+})$. Thus, by construction, we have $\pi_{x}i_{y}^{p}u\in V_{h,M}$.

For the first term on the right-hand side of the estimate in Lemma 4.1, we have the following result.

Lemma 4.2**.**

Let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Then, it holds

[TABLE]

Proof.

First, we set $v_{h}=\pi_{x}i_{y}^{p}u\in V_{h,M}$. Then, by introducing $\pi_{x}u$ as an intermediate function, we deduce

[TABLE]

According to the definition of $\pi_{x}$, we have $\partial_{y}\pi_{x}(u-i_{y}^{p}u)=\pi_{x}\partial_{y}(u-i_{y}^{p}u)$ almost everywhere in $C_{Y}$. As a consequence, we deduce by well known stability estimates for the $L^{2}$ projection $\pi_{x}$

[TABLE]

which shows the assertion. ∎

Lemma 4.3**.**

For $\mathfrak{f}\in\mathds{H}^{1-s}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Then, it holds

[TABLE]

Proof.

Analogously to the foregoing proof, by classical estimates for the $L^{2}$ projection $\pi_{x}$, we obtain

[TABLE]

Then, Proposition 2.5 yields the assertion. ∎

Next, we are concerned with estimates for the second term in Lemma 4.2. To this end, we consider the interpolation error on each subinterval $I_{m}\in\mathcal{I}_{Y}$. Employing the decomposition from Proposition 2.6 and the definition of the eigenfunctions $\varphi_{k}$, we obtain

[TABLE]

By using this identity in the following two subsections, we will estimate the terms

[TABLE]

for the two types of triangulations $\mathcal{I}_{Y}$ and polynomial spaces $V_{M}$ introduced in the Section 3.1 and 3.2.

Thereby, in Section 4.1, we will mainly recover the results of [27]. The reason for doing this is twofold. First, we are able to slightly improve the grading condition from $\mu<\frac{2}{3}s$ (in our notation) of [27, Section 5.2] to $\mu<s$. However, the main reason to analyze the $h$-FEM on graded meshes before developing the analysis for the considered $hp$-method is, that the techniques we use are almost identical for both cases, but the details are simpler for $h$-FEM, of course.

Later, in Section 4.2, we will analyze the $hp$-method introduced in Section 3.2, which yields a slightly improved rate of convergence ($h_{\Omega}$ vs. $h_{\Omega}\lvert\ln h_{\Omega}\rvert^{s}$) compared to $h$-FEM but a drastic reduction of the computational complexity in terms of degrees of freedom from $O(\mathcal{N}_{\Omega}^{1+1/d})$ to $O(\mathcal{N}_{\Omega}(\ln\mathcal{N}_{\Omega})^{2})$.

Note that in the following estimates, we will track the dependence on $Y$ explicitly since as a last step $Y$ will be chosen $h$-dependent.

4.1 Graded meshes and $h$-FEM

As announced, we are concerned with estimates for (4.4) for the discretization defined in Section 3.1. For simplicity, in this subsection, we will write $i_{y}u$ for $i_{y}^{p}u$, since we have $p=(1,1,\dots,1)$ here.

Lemma 4.4** (Estimate on $I_{1}$).**

For $\mathfrak{f}\in L^{2}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2) and let $M\geq h_{\Omega}^{-1}$. Then, it holds

[TABLE]

Proof.

First, we observe that the interpolant $i_{y}\psi_{s,k}$ is constant on $I_{1}$. By its definition (4.2), it holds (i_{y}\psi_{s,k})\bigr{\rvert}_{I_{1}}=\psi_{s,k}(y_{1}). Integration by parts and noting that y^{\alpha+1}(\psi_{s,k}-\psi_{s,k}(y_{1}))^{2}\bigr{\rvert}_{0}^{y_{1}}=0 yields

[TABLE]

Then, dividing by $\lVert y^{\frac{\alpha}{2}}(\psi_{s,k}-i_{y}\psi_{s,k})\rVert_{L^{2}(I_{1})}$ implies

[TABLE]

with $c=\frac{1}{1-s}$. By means of Corollary 2.7 with $n=1$ and $r=1$, we obtain

[TABLE]

Hence, we get by Lemma 3.1 together with the assumption on $M$ that

[TABLE]

In a similar fashion, we obtain by Corollary 2.7 with $n=1$ and $r=0$ the relation

[TABLE]

where we again used Lemma 3.1 in the last step.

The previous estimates together with (4.4) and Proposition 2.1 yield

[TABLE]

which implies the assertion. ∎

Lemma 4.5** (Estimates on $I_{m}$ for $2\leq m\leq M-1$).**

For $\mathfrak{f}\in L^{2}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Moreover, let $M\geq h_{\Omega}^{-1}$ and $\mu\neq s$. Then, it holds

[TABLE]

Proof.

For $m\geq 2$, we have that $y_{m-1}\leq y_{m}\leq 2^{\frac{1}{\mu}}y_{m-1}$. It holds $i_{y}\psi_{s,k}=i_{1}\psi_{s,k}$ such that we conclude with standard estimates for the linear Lagrange interpolant $i_{1}$, Lemma 3.1, and the assumption on $M$ that

[TABLE]

By using Corollary 2.7 with $n=1$ and $r=1$, this implies

[TABLE]

Similarly, using Corollary 2.7 with $n=1$ and $r=0$, we obtain for the term involving the derivative

[TABLE]

The previous estimates in combination with (4.4) yield

[TABLE]

Finally, applying Proposition 2.1, we get

[TABLE]

which states the assertion. ∎

Lemma 4.6** (Estimate on $I_{M}$).**

For $\mathfrak{f}\in L^{2}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Moreover, let $2h_{\Omega}^{-1}\geq M\geq h_{\Omega}^{-1}$, $\mu\neq s$, and

[TABLE]

Then, it holds

[TABLE]

Proof.

We recall that $Y=y_{M}$ and $i_{y}\psi_{s,k}=\tilde{i}_{1}\psi_{s,k}$ on $I_{M}$. We introduce the Lagrange interpolation $i_{1}$ on $I_{M}$ as an intermediate function such that

[TABLE]

where we used that $\lVert l_{1,1}\rVert_{L^{\infty}(I_{M})}=\lVert\frac{y-y_{M-1}}{h_{M}}\rVert_{L^{\infty}(I_{M})}=1$ in the last step. As in the proof of Lemma 4.5, we deduce

[TABLE]

Since $Y\geq 1$ by assumption, we obtain using Corollary 2.8 with $r_{1}=\frac{\alpha+1}{2}=1-s$ together with the monotonicity of $e^{-\sqrt{\lambda_{k}}y}$

[TABLE]

where we notice that $(\lambda_{k})_{k\in\mathds{N}}$ is a non-decreasing sequence. Combining the previous results yields

[TABLE]

Similarly, we deduce

[TABLE]

where we used that $\lVert l_{1,1}^{\prime}\rVert_{L^{\infty}(I_{M})}=\lVert h_{M}^{-1}\rVert_{L^{\infty}(I_{M})}=h_{M}^{-1}$ in the last step. The first term can again be estimated as in the proof Lemma 4.5 such that

[TABLE]

Employing Corollary 2.8 with $r_{1}=\frac{\alpha}{2}=\frac{1}{2}-s$ together with the monotonicity of $e^{-\sqrt{\lambda_{k}}y}$ yields for $Y\geq 1$

[TABLE]

where we used once again that the sequence $(\lambda_{k})_{k\in\mathds{N}}$ is non-decreasing. Due to the previous results, we arrive at

[TABLE]

By combining (4.4), (4.5) and (4.6), we obtain

[TABLE]

where we used $\sum_{k=1}^{\infty}\lambda_{k}^{2s}\mathfrak{u}_{k}^{2}=\sum_{k=1}^{\infty}\mathfrak{f}_{k}^{2}$ (see Proposition 2.1) and the definition of $\lVert\mathfrak{f}\rVert^{2}_{L^{2}(\Omega)}$. According to Lemma 3.1, there holds $h_{M}^{-1}\leq cY^{\mu-1}Y^{-\mu}M=cY^{-1}M$. Since

[TABLE]

by assumption, we obtain

[TABLE]

which ends the proof. ∎

Corollary 4.7**.**

For $\mathfrak{f}\in L^{2}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Moreover, let $2h_{\Omega}^{-1}\geq M\geq h_{\Omega}^{-1}$, $\mu<s$, and

[TABLE]

Then, it holds

[TABLE]

Proof.

By the Lemmas 4.4, 4.5, and 4.6, we obtain

[TABLE]

where we have used $\mu<s$ and the upper bound on $Y$. ∎

Now, we are able to state the main result for this subsection analyzing the $h$-FEM on graded meshes.

Theorem 4.8**.**

For $\mathfrak{f}\in\mathds{H}^{1-s}(\Omega)$, let $\mathfrak{u}\in\mathds{H}^{s}(\Omega)$ and $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solutions of (1.1) and (2.2), respectively, and let $u_{h}\in V_{h,M}$ be the solution of (3.1). Moreover, let $2h_{\Omega}^{-1}\geq M\geq h_{\Omega}^{-1}$, $\mu<s$, and

[TABLE]

Then, it holds

[TABLE]

Proof.

The first inequality of the assertion is due to Propositions 2.3 and 2.4. Using the Lemmas 4.1 and 4.2, we get

[TABLE]

The three terms on the right-hand side are estimated in Lemma 4.3, Corollary 4.7, and Proposition 2.9. Hence, we get

[TABLE]

Then, the lower bound on $Y$ yields $e^{-\frac{\sqrt{\lambda_{1}}}{2}Y}\leq ch_{\Omega}^{\frac{3}{2}}\leq ch_{\Omega}$, which implies the assertion. ∎

Theorem 4.9**.**

The total number degrees of freedom $\mathcal{N}_{\Omega,Y}$ in $V_{h,M}$ to achieve the order of convergence given in Theorem 4.8 behaves like

[TABLE]

where $d$ denotes the dimension of $\Omega$.

Proof.

For the number of degrees of freedom $\mathcal{N}_{\Omega,Y}$ of the discretization considered in this section, it holds $\mathcal{N}_{\Omega,Y}=\mathcal{N}_{\Omega}\mathcal{N}_{Y}=\mathcal{N}_{\Omega}M$. Then, the assertion follows from $M=O(h_{\Omega}^{-1})=O(\mathcal{N}_{\Omega}^{1/d})$. ∎

4.2 Geometric meshes and $hp$-FEM

In this section, we derive discretization error estimates for the $hp$-method described in Section 3.2, which results in a slightly improved rate of convergence of $h_{\Omega}$ compared to the previous subsection. However, we will have a drastically reduced computational complexity in terms of the number of degrees of freedom. To this end, we do neither fix the number of elements $M$ in direction $y$ nor the slope $\beta$ of the linear degree vector $p$ yet. These will be set below. As before, we start with estimates for $\lVert y^{\frac{\alpha}{2}}\nabla(u-i_{y}^{p}u)\rVert_{L^{2}(\Omega\times I_{m})}^{2}$ based on (4.4).

Lemma 4.10** (Estimate on $I_{1}$).**

For $\mathfrak{f}\in L^{2}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2) and let

[TABLE]

for some $\varepsilon\geq 0$. Then, it holds

[TABLE]

Proof.

Notice that $i_{y}^{p}\psi_{s,k}=\psi_{s,k}(y_{1})$ on the first interval $I_{1}$ as in the previous section. Thus, as in the proof of Lemma 4.4 but using

[TABLE]

from Lemma 3.2 and the assumption on $M$, we get

[TABLE]

Hence, we obtain

[TABLE]

As in the proof of Lemma 4.4, this yields the assertion. ∎

In order to derive estimates on $I_{m}$ for $2\leq m\leq M-1$, we recall the following result which is a direct consequence of [24, Lemma 3.2.6].

Proposition 4.11**.**

Let $w$ be analytic on $\hat{I}=(0,1)$ and satisfy for some $c_{w},\delta>0$ the estimate

[TABLE]

Then, there are constants $c,b>0$ depending only on $\delta$ such that the Gauss-Lobatto interpolant $i_{q}w$ of degree $q\in\mathds{N}$ on $\hat{I}$ satisfies

[TABLE]

Lemma 4.12** (Estimates on $I_{m}$ for $2\leq m\leq M-1$).**

For $\mathfrak{f}\in L^{2}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Moreover, let $p\in\mathds{N}^{M}$ be a linear degree vector as in (3.3) with some $\beta>0$ and let

[TABLE]

for some $\varepsilon\geq 0$, where $b>0$ is a constant depending on $\sigma$ only. Then, it holds

[TABLE]

Proof.

For $m\geq 2$ it holds $y_{m-1}\leq y_{m}\leq\sigma^{-1}y_{m-1}$. By transforming to the reference element $\hat{I}=(0,1)$, we obtain for (i_{y}^{p}\psi_{s,k})\bigr{\rvert}_{I_{m}}=i_{p_{m}}\psi_{s,k}\bigr{\rvert}_{I_{m}} that

[TABLE]

By means of Corollary 2.7, applied with $r=1$, and Lemma 3.2, we have

[TABLE]

Moreover, due to well known series representations of $K_{s}$ from [1, 9.6.2 and 9.6.10], we directly conclude that $\hat{\psi}_{s,k}$ is analytic on $\hat{I}$. Hence, Proposition 4.11 implies with $\delta=8(\sigma^{-1}-1)$ that

[TABLE]

Then, we get by Lemma 3.2

[TABLE]

Analogously, we obtain

[TABLE]

By means of Corollary 2.7 with $r=0$ and Lemma 3.2, we have

[TABLE]

Consequently, Proposition 4.11 yields with $\delta=8(\sigma^{-1}-1)$ that

[TABLE]

This implies

[TABLE]

Collecting the previous results yields

[TABLE]

Relation (3.3) implies

[TABLE]

Thus, we deduce

[TABLE]

Let us now distinguish two cases:

• •

$s\leq\beta b$: Since $h_{m}\geq h_{1}$, we then obtain

[TABLE]

As before, the relation $M\geq\frac{(1+\varepsilon)\lvert\ln h_{\Omega}\rvert}{s\lvert\ln\sigma\rvert}$ together with Lemma 3.2 implies $h_{1}\leq ch_{\Omega}^{\frac{1+\varepsilon}{s}}Y$. Hence, we get

[TABLE]

• •

$s>\beta b$: With $h_{m}\leq Y$, we get from (4.7) that

[TABLE]

Similarly as before, the relation $M\geq\frac{(1+\varepsilon)\lvert\ln h_{\Omega}\rvert}{\beta b\lvert\ln\sigma\rvert}$ together with Lemma 3.2 implies $h_{1}\leq ch_{\Omega}^{\frac{1+\varepsilon}{\beta b}}Y$. Thus, we get

[TABLE]

The previous results in combination with (4.4) imply

[TABLE]

Finally, applying Proposition 2.1 yields the assertion. ∎

Lemma 4.13** (Estimate on $I_{M}$).**

For $\mathfrak{f}\in L^{2}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Moreover, let $p\in\mathds{N}^{M}$ be a linear degree vector as in (3.3) with some $\beta>0$ and let

[TABLE]

for some $\varepsilon\geq 0$, where $b>0$ is a constant depending on $\sigma$ only. Then, it holds

[TABLE]

Proof.

We proceed as in the proof of Lemma 4.6 and recall that $Y=y_{M}=\sigma^{-1}y_{M-1}$ and $h_{M}=(1-\sigma)Y=(\sigma^{-1}-1)y_{M-1}$. Moreover, according to [17], the Lagrange basis functions $l_{i,q}$ of order $q\in\mathds{N}$ on $I_{M}$ have the property

[TABLE]

As a consequence, we obtain by means of an inverse inequality (see, e.g., [24, Lemma 3.2.2])

[TABLE]

Noting that $i_{y}^{p}\psi_{s,k}=\tilde{i}_{p_{M}}\psi_{s,k}$ on $I_{M}$, we introduce the Gauss-Lobatto interpolant $i_{p_{M}}\psi_{s,k}$ on $I_{M}$ as an intermediate function such that

[TABLE]

where we used (4.8) in the last step. As in the proof of Lemma 4.12, we deduce for $M\geq\frac{(1+\varepsilon)\lvert\ln h_{\Omega}\rvert}{\min(s,\beta b)\lvert\ln\sigma\rvert}$ that

[TABLE]

Since $Y\geq 1$ by assumption, we obtain using Corollary 2.8 with $r_{1}=\frac{\alpha+1}{2}=1-s$ together with the monotonicity of $e^{-\sqrt{\lambda_{k}}y}$

[TABLE]

where we notice that $(\lambda_{k})_{k\in\mathds{N}}$ is a non-decreasing sequence. Combining the previous results yields

[TABLE]

Similarly, we deduce by means of (4.9)

[TABLE]

Using $M\geq\frac{(1+\varepsilon)\lvert\ln h_{\Omega}\rvert}{\min(s,\beta b)\lvert\ln\sigma\rvert}$, the first term can again be estimated as in the proof Lemma 4.12 such that

[TABLE]

Employing Corollary 2.8 with $r_{1}=\frac{\alpha+1}{2}+2=3-s$ together with the monotonicity of $e^{-\sqrt{\lambda_{k}}y}$ yields for $Y\geq 1$

[TABLE]

where we used once again that the sequence $(\lambda_{k})_{k\in\mathds{N}}$ is non-decreasing. Due to the previous results, we arrive at

[TABLE]

By combining (4.4), (4.10) and (4.11), we obtain as in the proof of Lemma 4.6

[TABLE]

According to Lemma 3.3, there holds $p_{M}\leq cM$. Since

[TABLE]

by assumption, we obtain

[TABLE]

which ends the proof. ∎

Corollary 4.14**.**

For $\mathfrak{f}\in L^{2}(\Omega)$, let $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solution of (2.2). Moreover, let $p\in\mathds{N}^{M}$ be a linear degree vector as in (3.3) with some $\beta>0$ and let

[TABLE]

for some $\varepsilon>0$, where $b>0$ is a constant depending on $\sigma$ only. Then, it holds

[TABLE]

Proof.

By the Lemmas 4.10, 4.12, and 4.13, we obtain

[TABLE]

where we have used the upper bounds on $M$ and $Y$ and the boundedness of $h_{\Omega}^{2\varepsilon}\lvert\ln h_{\Omega}\rvert^{2s+1}$ for $\varepsilon>0$. ∎

Now, we are able to state the main result for this subsection analyzing the $hp$-FEM on geometric meshes.

Theorem 4.15**.**

For $\mathfrak{f}\in\mathds{H}^{1-s}(\Omega)$, let $\mathfrak{u}\in\mathds{H}^{s}(\Omega)$ and $u\in\mathring{H}^{1}_{L}(C,y^{\alpha})$ be the solutions of (1.1) and (2.2), respectively, and let $u_{h}\in V_{h,M}$ be the solution of (3.1). Moreover, let $p\in\mathds{N}^{M}$ be a linear degree vector as in (3.3) with some $\beta>0$ and let

[TABLE]

for some $\varepsilon>0$, where $b>0$ is a constant depending on $\sigma$ only. Then, it holds

[TABLE]

Proof.

The first inequality of the assertion is due to Propositions 2.3 and 2.4. Using Lemmas 4.1 and 4.2, we get

[TABLE]

The three terms on the right-hand side are estimated in Lemma 4.3, Corollary 4.14, and Proposition 2.9. Hence, we get

[TABLE]

Then, the lower bound on $Y$ yields $e^{-\frac{\sqrt{\lambda_{1}}}{2}Y}\leq ch_{\Omega}$, which implies the assertion. ∎

Theorem 4.16**.**

The total number degrees of freedom $\mathcal{N}_{\Omega,Y}$ in $V_{h,M}$ to achieve the order of convergence given in Theorem 4.15 behaves like

[TABLE]

Proof.

As a direct consequence of Lemma 3.4, we obtain that the number of degrees of freedom $\mathcal{N}_{\Omega,Y}$ of the discretization considered in this section fulfills $\mathcal{N}_{\Omega,Y}=\mathcal{N}_{\Omega}\mathcal{N}_{Y}=O(\mathcal{N}_{\Omega}M^{2})$. Then, the assertion follows from $M=O(\lvert\ln h_{\Omega}\rvert)=O(\ln\mathcal{N}_{\Omega})$. ∎

5 Numerical Experiments

5.1 Implementation

For the discretization with respect to $y$ in $h$-FEM and $hp$-FEM, we use hierarchical Lobatto polynomials (see, e.g., [29]) as local shape functions on $\hat{I}=(0,1)$ which then are transformed onto each interval $I_{m}$. Without the weight, i.e., for $s=\frac{1}{2}$, this would result in a very sparse structure of the local stiffness matrix, since those shape functions are orthogonal. For $s\neq\frac{1}{2}$ the latter does not hold; nevertheless the global matrix is of course still sparse.

Let $\eta_{i}$ for $i=1,2,\dots,\mathcal{N}_{\Omega}$ be the ansatz functions for the discretization of $\Omega$ and $\tau_{j}$ for $j=1,2,\dots,\mathcal{N}_{Y}$ ansatz functions for the discretization of $(0,Y)$. On the cylinder $C_{Y}$ we then use ansatz functions of the form $\varphi_{i,j}(x,y)=\eta_{i}(x)\tau_{j}(y)$ with $i=1,2,\dots,\mathcal{N}_{\Omega}$ and $j=1,\dots,\mathcal{N}_{Y}$.

Due to this special structure, the system matrix $S\in\mathds{R}^{\mathcal{N}_{\Omega,Y}\times\mathcal{N}_{\Omega,Y}}$ for solving (3.1) can be expressed by means of the Kronecker product as

[TABLE]

Here, $A\in\mathds{R}^{\mathcal{N}_{\Omega}\times\mathcal{N}_{\Omega}}$ denotes matrices arising from discretization of $\Omega$ and $B\in\mathds{R}^{\mathcal{N}_{Y}\times\mathcal{N}_{Y}}$ denotes matrices arising from discretization of $(0,Y)$ given as

[TABLE]

We observe that one can assemble the matrices $A$ and $B$ completely independent from each other. This is advantageous since the weight $y^{\alpha}$ only affects the $B$ matrices, while the $A$ matrices are standard FEM matrices, which can be computed by any FEM software.

Using the special structure of $S$, one can implement a memory efficient solvers for the algebraic systems without ever fully assembling $S$. This will be the topic of a forthcoming paper.

5.2 Numerical Results

We take the following configuration from [27, Section 6.1]. For $\Omega=(0,1)^{2}\subset\mathds{R}^{2}$ the eigenfunctions of the Dirichlet-Laplacian are known to be $\varphi_{k,l}(x)=\sin(k\pi x_{1})\sin(l\pi x_{2})$ with corresponding eigenvalues $\lambda_{k,l}=\pi^{2}(k^{2}+l^{2})$ for $k,l=1,2,\dots$. For the right-hand side $\mathfrak{f}(x)=\lambda_{1,1}^{s}\varphi_{1,1}(x)=(2\pi^{2})^{s}\sin(\pi x_{1})\sin(\pi x_{2})$, the solution $\mathfrak{u}$ of (1.1) and $u$ of (2.2) are then given by

[TABLE]

For the discretization by means of $h$-FEM (cf. the Section 3.1 and 4.1), we choose the parameters

[TABLE]

whereas for the discretization by means of $hp$-FEM (cf. the Section 3.2 and 4.2), we choose the following parameters:

[TABLE]

The orders of convergence stated by the Theorems 4.8 and 4.15 in terms of $h_{\Omega}$ are confirmed by the results of the numerical computations given in Figure 1. Note, that the error $\lVert y^{\frac{\alpha}{2}}\nabla(u-u_{h})\rVert_{L^{2}(C)}$ is evaluated by means of the identity

[TABLE]

which holds due to the Galerkin orthogonality (4.1).

In Figure 2, we depict the errors for both types of discretizations over the total numbers of degrees of freedom $\mathcal{N}_{\Omega,Y}$. Thereby, the slower growth of $\mathcal{N}_{\Omega,Y}$ for the $hp$-discretization given by Theorem 4.16 in comparison to Theorem 4.9 clearly leads to a drastic reduction of the number of degrees of freedom compared to $h$-FEM on a graded mesh. For instance the number of degrees of freedom to achieve an error of less than $9\cdot 10^{-3}$ in the case $s=0.8$ reduces from $\mathcal{N}_{\Omega,Y}=1\,072\,692\,225$ for $h$-FEM to $\mathcal{N}_{\Omega,Y}=9\,661\,477$ for $hp$-FEM, which is a factor of about $111$.

Appendix A Estimates for $\psi_{s}$ and its derivatives

We begin with a representation of the derivatives of the expression $z^{s}K_{s}(z)$, where $K_{s}$ are the modified Bessel functions of second kind. It which will be used in the sequel to derive estimates for the derivatives of $\psi_{s}$.

Lemma A.1**.**

The derivatives of $z^{s}K_{s}(z)$ of order $n\in\mathds{N}_{0}$ can be calculated as

[TABLE]

where the coefficients $a_{m}^{n}$ are given by

[TABLE]

Proof.

We prove this assertion by induction. To this end, we first collect some basic results for the modified Bessel function of second kind. In [1, 9.6.28], we find for all $\nu\in\mathds{R}$

[TABLE]

As a consequence, there holds

[TABLE]

Using the latter result, we get the following formula for $m,l\in\mathds{N}_{0}$:

[TABLE]

By means of this, we obtain for $m,n\in\mathds{N}_{0}$ with $m\leq n$ by setting $l=n-m\geq 0$

[TABLE]

These elementary results build the basis for the induction: The hypothesis (A.1) clearly holds for $n=0$.

Assuming that (A.1) holds for some $n\in\mathds{N}_{0}$, we deduce

[TABLE]

Employing (A.6), we continue with

[TABLE]

It remains to show that

[TABLE]

The first and third equation are obvious due to (A.2) and (A.4). Thus, we only elaborate on the second. We distinguish three cases for $1\leq m\leq n$:

• •

m\geq\bigl{\lfloor}\frac{n+1}{2}\bigr{\rfloor}+2: Again due to (A.4), we have $a_{m}^{n}=a_{m-1}^{n}=0$, since m,m-1>\bigl{\lfloor}\frac{n}{2}\bigr{\rfloor}. Hence, it holds

[TABLE]

• •

m=\bigl{\lfloor}\frac{n+1}{2}\bigr{\rfloor}+1: Here, it holds m>\bigl{\lfloor}\frac{n}{2}\bigr{\rfloor} and we already know that $a_{m}^{n}=0$. Moreover, in case that $n$ is even, we deduce $m=\frac{n}{2}+1$ and

[TABLE]

If $n$ is odd, we obtain $a_{m-1}^{n}=0$ since m-1=\frac{n+1}{2}>\bigl{\lfloor}\frac{n}{2}\bigr{\rfloor}. As a consequence, we get

[TABLE]

• •

m\leq\bigl{\lfloor}\frac{n+1}{2}\bigr{\rfloor}: Here, we again distinguish between $n$ even and $n$ odd. In the first case, we have \bigl{\lfloor}\frac{n+1}{2}\bigr{\rfloor}=\bigl{\lfloor}\frac{n}{2}\bigr{\rfloor}. Hence, it holds m-1,m\leq\bigl{\lfloor}\frac{n}{2}\bigr{\rfloor} and by means of (A.3), we get

[TABLE]

In case that $n$ is odd, we have that \bigl{\lfloor}\frac{n+1}{2}\bigr{\rfloor}=\bigl{\lfloor}\frac{n}{2}\bigr{\rfloor}+1. Thus, for m\leq\bigl{\lfloor}\frac{n+1}{2}\bigr{\rfloor}-1, we can reuse the calculations from before. If m=\bigl{\lfloor}\frac{n+1}{2}\bigr{\rfloor}=\frac{n+1}{2}, we have $a^{n}_{m}=0$ such that

[TABLE]

This ends the proof. ∎

We next analyze $\psi_{s}$ defined in (2.3) and its derivatives with respect to its boundedness properties.

Lemma A.2**.**

For $z\geq 0$, it holds

[TABLE]

Proof.

Since $z^{\nu}K_{\nu}(z)>0$ for all $z>0$ and $\nu>-1$, see, e.g., [1, 9.6.1], and due to (A.5), the function $\psi_{s}$ is positive and monotone decreasing such that $\psi_{s}(z)\leq\psi_{s}(0)$ for all $z\geq 0$.

In [1, 9.6.9] one can find for $\nu>0$ the following behavior of the modified Bessel function of the second kind for $z\rightarrow 0$:

[TABLE]

As a consequence, we obtain by (2.3)

[TABLE]

which yields together with the foregoing observations the assertion. ∎

Lemma A.3**.**

Let $r\in[0,1]$. There exists a constant $c>0$ depending only on $s$, such that for any $z>0$ and $n\in\mathds{N}$ it holds

[TABLE]

Proof.

In order to deduce the bounds for the derivatives of $\psi_{s}$, we continue with collecting some auxiliary results. As before in the proof of Lemma A.2, we have that $z^{r}z^{\nu}K_{\nu}(z)$ is positive for all $z>0$ and $\nu>-1$. Let \nu_{0}=\min\bigl{(}\nu,\frac{1}{2}\bigr{)}. From [25, Theorem 5], we obtain for $z>0$ that $z^{\nu_{0}}e^{z}K_{\nu}(z)$ is a decreasing function for all $\nu\geq 0$. To employ this, we consider the product

[TABLE]

and note that

[TABLE]

Hence, $z^{\nu+r}K_{\nu}(z)$ admits its maximum in the interval $[0,\nu+r-\nu_{0}]$. Due to Lemma A.2, we consequently get by (2.3)

[TABLE]

Next, we are concerned with the bounds for the derivatives of $\psi_{s}$. Employing Lemma A.1 and the relation $K_{\nu}(z)=K_{-\nu}(z)$ for $\nu\in\mathds{R}$, see [1, 9.6.6], we obtain by (2.3)

[TABLE]

where the coefficients are given by Lemma A.1. Let $\nu(m,n)=n-m-s$. We observe that $\nu(m,n)>0$ since m\leq\bigl{\lfloor}\frac{n}{2}\bigr{\rfloor}. Consequently, (A.8) yields

[TABLE]

with \nu_{0}(m,n)=\min\bigl{(}\nu(m,n),\frac{1}{2}\bigr{)}. Since $n-m\geq 1$, it holds

[TABLE]

Further, using

[TABLE]

from [23, estimate (8)], which holds for all $l\in\mathds{N}_{0}$ and $0<\rho<1$ , we get by choosing $l=n-m-1$ and $\rho=1-s$

[TABLE]

with a constant $c$ depending only on $s$. Using (A.9) and (A.10), we get

[TABLE]

Estimating each summand separately, yields by means of Lemma A.1

[TABLE]

For the last step, notice that

[TABLE]

Finally, (A.11) and (A.12) yield the assertion since \bigl{\lfloor}\frac{n}{2}\bigr{\rfloor}+1\leq 2^{n}. ∎

Finally, we state a result about the exponential decay of $\psi_{s}$ and its derivative.

Lemma A.4**.**

The following assertions hold:

1. (a)

Let $s_{0}=\min(s,\frac{1}{2})$. Moreover, let $z\geq a>0$ and $r\geq s_{0}-s$. Then there exists a constant $c$ only depending on $r$ and $s$ such that

[TABLE] 2. (b)

Let $s_{0}^{\prime}=\min(1-s,\frac{1}{2})$. Moreover, let $z\geq a>0$ and $r\geq s_{0}^{\prime}-s$. Then there exists a constant $c$ only depending on $r$ and $s$ such that

[TABLE]

Proof.

We start as in the proof of Lemma A.3. According to [25, Theorem 5], we get for $z\geq 0$ that $z^{s_{0}}e^{z}K_{s}(z)$ is a decreasing function. Consequently, having in mind the definition of $\psi_{s}$ and that $z^{\nu}K_{\nu}(z)>0$ for all $z\geq 0$ and $\nu>-1$, see the proof of Lemma A.2, we obtain

[TABLE]

This is already the desired result for $r=s_{0}-s$ noticing that $e^{-\frac{z}{2}}<1=c$ for $z\geq 0$. For $r>s_{0}-s$, we observe that

[TABLE]

where we used that

[TABLE]

Combining the previous results yields the first inequality of the assertion. Next, we deduce by means of the definition of $\psi_{s}$ and (A.5)

[TABLE]

such that the second inequality of the assertion follows from the first one noting that $r+2s-1\geq s_{0}^{\prime}-s+2s-1=\min(1-s,\frac{1}{2})-(1-s)$ by the assumption on $r$. ∎

Acknowledgement

The authors acknowledge Christian Kahle’s and Martin Stoll’s support in the efficient solution of the arising linear systems.