Finite element approximations for fractional evolution problems

TL;DR

This paper develops and analyzes a finite element scheme for fractional evolution problems involving both time and space derivatives, capturing memory and long-range dispersion effects, with numerical validation in 1D and 2D.

Contribution

It introduces a novel finite element method combining convolution quadrature and linear elements for fractional evolution equations, with theoretical analysis and numerical experiments.

Findings

The scheme is well-posed and stable.

Numerical experiments confirm accuracy in 1D and 2D.

Regularity estimates support the method's convergence.

Abstract

This work introduces and analyzes a finite element scheme for evolution problems involving fractional-in-time and in-space differentiation operators up to order two. The left-sided fractional-order derivative in time we consider is employed to represent memory effects, while a nonlocal differentiation operator in space accounts for long-range dispersion processes. We discuss well-posedness and obtain regularity estimates for the evolution problems under consideration. The discrete scheme we develop is based on piecewise linear elements for the space variable and a convolution quadrature for the time component. We illustrate the method's performance with numerical experiments in one- and two-dimensional domains.

Tables

Table 1. Table 1. The L 2 ​ ( Ω ) superscript 𝐿 2 Ω L^{2}(\Omega) error at time t = 0.1 𝑡 0.1 t=0.1 with s = 0.75 𝑠 0.75 s=0.75 using a uniform mesh with size h = 1 / 5000 ℎ 1 5000 h=1/5000 . The expected convergence rate in τ 𝜏 \tau is 1 1 1 .

Example	$\alpha \setminus \tau$	0.01	0.005	0.0025	0.001	Rate (in $\tau$)
(a)	0.5	4.227e-3	2.105e-3	1.055e-3	4.493e-4	0.98
(a)	1.5	2.512e-2	1.261e-2	6.349e-3	2.602e-3	0.99
(b)	1.5	4.867e-3	2.402e-3	1.188e-3	5.362e-4	0.96

Table 2. Table 2. The L 2 ​ ( Ω ) superscript 𝐿 2 Ω L^{2}(\Omega) error at time t = 0.1 𝑡 0.1 t=0.1 with s = 0.75 𝑠 0.75 s=0.75 using τ = 1 / 5000 𝜏 1 5000 \tau=1/5000 . The expected convergence rate in h ℎ h is 1.

Example	$\alpha \setminus$ mesh size $h$	$1/250$	$1/500$	$1/1000$	$1/1500$	Rate (in $h$)
(a)	0.5	8.837e-3	3.856e-3	1.781e-3	1.198e-3	1.12
(a)	1.5	9.978e-3	4.350e-3	1.967e-3	1.252e-3	1.16
(b)	0.5	1.162e-3	5.158e-4	4.010e-3	1.640e-4	1.09
(b)	1.5	8.453e-4	3.644e-4	1.632e-4	1.035e-4	1.17

Table 3. Table 3. The L 2 ​ ( Ω ) superscript 𝐿 2 Ω L^{2}(\Omega) error at time t = 0.1 𝑡 0.1 t=0.1 with s = 0.25 𝑠 0.25 s=0.25 using τ = 1 / 5000 𝜏 1 5000 \tau=1/5000 . The expected convergence rate in h ℎ h is 0.5.

Example	$\alpha \setminus$ mesh size $h$	$1/250$	$1/500$	$1/1000$	$1/1500$	Rate (in $h$)
(a)	0.5	8.571e-2	4.999e-2	2.924e-2	2.139e-2	0.77
(a)	1.5	1.125e-1	6.596e-2	3.859e-2	2.818e-2	0.77
(b)	0.5	1.171e-2	6.845e-3	4.010e-3	2.937e-3	0.77
(b)	1.5	1.154e-2	6.811e-3	4.014e-3	2.943e-3	0.76

Table 4. Table 4. The L 2 ​ ( Ω ) superscript 𝐿 2 Ω L^{2}(\Omega) error at time t = 0.02 𝑡 0.02 t=0.02 with s ∈ { 0.25 , 0.75 } 𝑠 0.25 0.75 s\in\{0.25,0.75\} and α = 0.8 𝛼 0.8 \alpha=0.8 , using τ = 1 / 5000 𝜏 1 5000 \tau=1/5000 for example (c). The observed rates are 0.77 0.77 0.77 and 1.38 1.38 1.38 , respectively.

Mesh size $h$	$s=0.25$	$s=0.75$
0.1	1.790e-1	5.673e-2
0.05	1.102e-1	2.342e-2
0.03	7.077e-2	1.054e-2
0.02	5.206e-2	6.255e-3

Equations

$$\text{leftidx}{}^C{\partial }_t^{\alpha }u+(-\Delta {)}^{s}u=f\text{mbox}in\Omega \times (0,T).$$

$$\leftidx{{}^{C}}\partial_{t}^{\alpha}u(x,t)=\left\{\begin{array}[]{rl}\frac{1}{\Gamma(k-\alpha)}\int_{0}^{t}\frac{1}{(t-r)^{\alpha-k+1}}\frac{\partial^{k}u}{\partial r^{k}}(x,r)\,dr&\mbox{ if }k-1<\alpha<k,\ k\in{\mathbb{N}},\\ \frac{\partial^{k}u}{\partial t^{k}}u(x,t)&\mbox{ if }\alpha=k\in{\mathbb{N}},\end{array}\right.$$

$$(-\Delta {)}^{s}u(x)=C(n,s)\text{mbox}P.V.{\int }_{{\mathbb{R}}^n}\frac{u(x)-u(y)}{|x-y{|}^{n+2s}}dy.$$

$$\partial_{t}^{\alpha}u(x,t)=\left\{\begin{array}[]{rl}\frac{1}{\Gamma(k-\alpha)}\frac{\partial^{k}}{\partial t^{k}}\int_{0}^{t}\frac{1}{(t-r)^{\alpha-k+1}}u(x,r)\,dr&\mbox{ if }k-1<\alpha<k,\ k\in{\mathbb{N}},\\ \frac{\partial^{k}u}{\partial t^{k}}u(x,t)&\mbox{ if }\alpha=k\in{\mathbb{N}}.\end{array}\right.$$

$$G_{\alpha ,s}(x,t)=t^{\frac{-\alpha }{2s}}{\Phi }_{\alpha ,s}\left(\frac{x}{t^{\frac{\alpha }{2s}}}\right),$$

$$\left\{\begin{array}[]{rl}u=0&\mbox{ in }\Omega^{c}\times[0,T),\\ u(\cdot,0)=v&\mbox{ in }\Omega,\end{array}\right.$$

$${\partial }_{t}u(\cdot ,0)=b\text{mbox}in\Omega .$$

$$H^s(\Omega )=\left\{w\in L^2(\Omega ):|w{|}_{H^s(\Omega )}:={\left({\iint }_{{\Omega }^2}\frac{|w(x)-w(y){|}^2}{|x-y{|}^{n+2s}}dxdy\right)}^{\frac{1}{2}}<\infty \right\}.$$

$$H^s(\Omega )=\left\{w\in H^m(\Omega ):|D^{\alpha }w{|}_{H^{\sigma }(\Omega )}<\infty \text{ for all }\alpha \text{ s.t }|\alpha |=m\right\}.$$

$${\tilde{H}}^s(\Omega )=\left\{w\in H^s({\mathbb{R}}^n):\text{ supp }w\subset \bar{\Omega }\right\}.$$

$$⟨u,w{⟩}_s:=C(n,s){\iint }_{({\mathbb{R}}^n\times {\mathbb{R}}^n)\setminus ({\Omega }^c\times {\Omega }^c)}\frac{(u(x)-u(y))(w(x)-w(y))}{|x-y{|}^{n+2s}}dydx$$

$$H^{-s}(\Omega )={\left({\tilde{H}}^s(\Omega )\right)}^{\prime }.$$

$${\int }_0^T\left(\text{leftidx}{}^C{\partial }_t^{\alpha }u(\cdot ,t)w(\cdot ,t)\right)dt+{\int }_0^T⟨u(\cdot ,t),w(\cdot ,t){⟩}_{s}dt={\int }_0^T\left(f(\cdot ,t),w(\cdot ,t)\right)dt$$

$$\left\{\begin{array}[]{rl}(-\Delta)^{s}u=g&\text{in }\Omega,\\ u=0&\text{in }\Omega^{c}.\end{array}\right.$$

$$|u{|}_{H^{s+\gamma }({\mathbb{R}}^n)}\le C\Vert g{\Vert }_{H^r(\Omega )}.$$

$$E_{\alpha ,\mu }(z)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\mu )}.$$

$$\text{leftidx}{}^C{\partial }_t^{\alpha }{E}_{\alpha ,1}(-\lambda t^{\alpha })=-\lambda E_{\alpha ,1}(-\lambda t^{\alpha }).$$

$$\frac{d}{dt}{E}_{\alpha ,1}(-\lambda t^{\alpha })=-\lambda t^{\alpha -1}{E}_{\alpha ,\alpha }(-\lambda t^{\alpha }),$$

$$\frac{d}{dt}\left(t{E}_{\alpha ,2}(-\lambda t^{\alpha })\right)=E_{\alpha ,1}(-\lambda t^{\alpha }).$$

$$|E_{\alpha ,\mu }(z)|\le \frac{C}{1+|z|}\text{mbox}for\delta <|\arg (z)|\le \pi .$$

$$\left\{\begin{array}[]{rl}(-\Delta)^{s}u=\lambda\,u&\mbox{in }\Omega,\\ u=0&\mbox{in }\Omega^{c}.\\ \end{array}\right.$$

$$0<{\lambda }_1<{\lambda }_2\le \dots ,{\lambda }_k\to \infty \text{mbox}ask\to \infty ,$$

$$u(x,t)=\sum _{k=1}^{\infty }{u}_k(t){\varphi }_k(x).$$

$$\left\{\begin{array}[]{rl}\leftidx{{}^{C}}\partial_{t}^{\alpha}u_{k}+\lambda_{k}u_{k}=f_{k},&\\ u_{k}(0)=v_{k},&\\ u^{\prime}_{k}(0)=b_{k},&\text{if }\alpha\in(1,2].\end{array}\right.$$

$$F_k(t)w={\int }_0^t(t-r{)}^{\alpha -1}{E}_{\alpha ,\alpha }(-{\lambda }_k(t-r{)}^{\alpha })w(r)dr,$$

$$u_k(t)$$

$$u_k(t)={\int }_0^t{e}^{-{\lambda }_k(t-r)}{f}_k(r)dr+v_k{e}^{-{\lambda }_{k}t},$$

$$u_k(t)=\frac{1}{\sqrt{{\lambda }_k}}{\int }_0^t\sin (\sqrt{{\lambda }_k}(t-r))f_k(r)dr+v_k\cos (\sqrt{{\lambda }_k}t)+b_k\frac{\sin (\sqrt{{\lambda }_k}t)}{\sqrt{{\lambda }_k}}.$$

$$\Vert u{\Vert }_{C([0,T];L^2(\Omega ))}\le C\Vert v{\Vert }_{L^2(\Omega )},$$

$$\Vert u(\cdot ,t){\Vert }_{H^{s+\gamma }(\Omega )}+\Vert \text{leftidx}{}^C{\partial }_t^{\alpha }u(\cdot ,t){\Vert }_{L^2(\Omega )}\le C{t}^{-\alpha }\Vert v{\Vert }_{L^2(\Omega )}.$$

Full text

1. Introduction

Since the introduction of Continuous Time Random Walks (CTRW) by Montroll and Weiss [39], anomalous diffusion phenomena has been an active area of research among the scientific community. The CTRW assign a joint space-time distribution to individual particle motions: when the tails of these distributions are heavy enough, non-Fickian dispersion results for all time and space scales. A heavy-tailed jump (waiting time) distribution implies the absence of a characteristic space (time) scale.

The equivalence between these heavy-tailed motions and transport equations that use fractional-order derivatives has been shown by several authors; see, for example [22]. Space nonlocality is a direct consequence of the existence of arbitrarily large jumps in space, whereas time nonlocality is due to the history dependence introduced in the dynamics by the presence of anomalously large waiting times.

The evidence of anomalous diffusion phenomena has been thoroughly reported in physical and social environments, such as plasma turbulence [14, 15], hidrology [8, 9, 42], finance [34], and human travel [13] and predator search [46] patterns. Models of transport dynamics in complex systems taking into account this non-Fickian behavior have been proposed accordingly. Also, evolution processes intermediate between diffusion and wave propagation have been shown to govern the propagation of stress waves in viscoelastic materials [18, 33].

Integer-order differentiation operators are local, because the derivative of a function at a given point depends only on the values of the function in a neighborhood of it. In contrast, fractional-order derivatives are nonlocal, integro-differential operators. A left-sided fractional-order derivative in time may be employed to represent memory effects, while a nonlocal differentiaton operator in space accounts for long-range dispersion processes.

We now describe the problems we are going to consider in this work. Let $\Omega\subset{\mathbb{R}^{n}}$ be a domain with smooth enough boundary, $\alpha\in(0,2]$, $s\in(0,1)$ and a forcing term $f\colon\Omega\times(0,T)\to{\mathbb{R}}$. We aim to solve the fractional differential equation

[TABLE]

Here, $\leftidx{{}^{C}}\partial_{t}^{\alpha}$ denotes the Caputo derivative, given by

[TABLE]

while $(-\Delta)^{s}$ is the fractional Laplace operator, defined as

[TABLE]

Above, $C(n,s)=\frac{2^{2s}s\Gamma(s+\frac{n}{2})}{\pi^{n/2}\Gamma(1-s)}$ is a normalization constant.

Closely related to the Caputo derivative, the Riemann-Liouville fractional derivative is needed in the sequel. Let us recall here its definition,

[TABLE]

For $0<\alpha\leq 1$, problem (1.1) is usually called a fractional diffusion equation. On the other hand, for $1<\alpha\leq 2$ it is sometimes called a fractional diffusion-wave equation. Analyzing scaling and similarity properties of the Green function $G_{\alpha,s}$ associated to the operator $\leftidx{{}^{C}}\partial_{t}^{\alpha}+(-\Delta)^{s}$, in [32] it is shown that

[TABLE]

for a certain one-variable function $\Phi_{\alpha,s}$. Notice that in case $\alpha=s$, although the CTRW associated to equation (1.1) has the same scaling properties as Brownian motion, the lack of finite moments makes the diffusion process anomalous. On the other hand, the term fractional wave equation has been utilized to refer to the problem with $1<\alpha=2s<2$, since for this choice of the parameters some features of the standard wave equation are preserved [31]. For example, the maximum, gravity center and mass center of the fundamental solution such possess constant propagation velocity.

Let $v,b\in L^{2}(\Omega)$ be given data. We complement problem (1.1) with the initial and boundary value conditions

[TABLE]

and, additionally for $1<\alpha\leq 2$, we require that

[TABLE]

A variational formulation for the time-fractional problem involving a Caputo derivative of order $\alpha\in(0,1)$ has been recently studied by Karkulik [27]. In that work, the author shows that the Caputo derivative is a linear and bounded operator on a time-fractional Sobolev-Bochner space, considers a variational formulation based exclusively on Sobolev regularity and proves that if the initial condition belongs to $H^{1-1/\alpha+\varepsilon}(\Omega)$ for some $\varepsilon>0$, then the time-fractional problem is well posed.

As fractional-order differential operators involve singular kernels, numerical approximations of equations involving them is a delicate task. Moreover, their nonlocal character calls for the design of efficient numerical schemes for the discretization. Several strategies have been shown to succeed in the one-dimensional context. For example, finite element schemes were treated in [17]; finite difference methods have been employed to discretize either the fractional-order space diffusion term [37, 38] and time derivative [28]; pseudospectral methods were considered in [22]. A comparison between these approaches may be found in [23]. Also, a convolution quadrature rule has been proposed in [29, 30], and in turn, has been utilized for the approximation of fractional time derivatives in [26].

Nevertheless, analysis and implementation of two-dimensional schemes is scarce. In [3], the Dirichlet problem for the fractional Laplace operator was analyzed and a simple finite element implementation for the two-dimensional problem was introduced. Later on, the authors presented the code employed for such implementation in [2]. The space discretization in the current work is based on such a code. Recently, implementations based on finite elements [5, 6] and integral representation formulas [11] have been proposed as well.

It is noteworthy that the fractional Laplace operator defined by (1.2) does not coincide with the operator considered, for example, in [10, 41, 47]. Indeed, the spatial operator considered in those works is a power of the Laplacian in the spectral sense.

Our work does not include the case $s=1$, which corresponds to a local-in-space process, as it is already covered by other authors’ work. For the range $0<\alpha\leq 1$, reference [25] develops a semidiscrete Galerkin method and studies the error both for smooth and non-smooth data. Naturally, the local-in-space case is also covered by the previously mentioned works [10, 41, 47] regarding spectral fractional powers of the Laplacian. For the full range of time derivatives we are considering in this work, [36] deals with an alternative formulation of (1.1) and a method based on the Laplace transform is developed, while in [40] an approach via discontinuous Galerkin discretization in time is introduced.

Organization of the paper

Preliminary concepts regarding fractional Sobolev spaces, elliptic regularity for the fractional Laplacian and the Mittag-Leffler function are discussed in Section 2. Moreover, that section also deals with the well-posedness of (1.1) with conditions (1.3), (1.4). Afterwards, in Section 3 we take advantage of the representation of solutions to derive regularity estimates for the fractional evolution problems under consideration. A numerical scheme, based on standard Galerkin finite element approximations for the space variable and a convolution quadrature for the time component, is proposed in Section 4, and error bounds for this scheme are presented in Section 5. In Section 6 we present some numerical examples that illustrate the accuracy of our convergence estimates and qualitative behavior of solutions to fractional evolution problems. Finally, Appendices A and B include details on the convolution quadrature rule utilized for the time discretization and on the error analysis of the semi-discrete scheme, respectively.

2. Preliminaries

In this section we set the basic notation and present some preliminary concepts necessary for the analysis of the fractional evolution problems under consideration. We recall elliptic regularity results for the fractional Laplacian and some important properties of the Mittag-Leffler function. These properties are then utilized to derive a representation formula for solutions that allows to prove the well-posedness of problem (1.1).

2.1. Sobolev spaces and fractional Laplace operator

Let $\Omega\subset{\mathbb{R}^{n}}$ be an open set and $s\in(0,1)$. The fractional Sobolev space $H^{s}(\Omega)$ is defined by

[TABLE]

This set, furnished with the norm $\|\cdot\|_{H^{s}(\Omega)}=\|\cdot\|_{L^{2}(\Omega)}+|\cdot|_{H^{s}(\Omega)},$ constitutes a Hilbert space. If $s>1$ and it is not an integer, considering the decomposition $s=m+\sigma$, where $m\in\mathbb{N}$ and $\sigma\in(0,1)$, allows to define $H^{s}(\Omega)$ by means of

[TABLE]

Another important space of interest for the problems under consideration is that of functions in $H^{s}({\mathbb{R}^{n}})$ supported within $\overline{\Omega}$,

[TABLE]

The bilinear form

[TABLE]

constitutes an inner product on $\widetilde{H}^{s}(\Omega)$. The norm it induces, which is just a multiple of the $H^{s}({\mathbb{R}^{n}})$-seminorm, is equivalent to the full $H^{s}({\mathbb{R}^{n}})$-norm on this space, because a Poincaré-type inequality holds in it. See, for example, [3] for details.

Remark 2.1*.*

We emphasize that functions in $\widetilde{H}^{s}(\Omega)$ are defined in the whole space ${\mathbb{R}^{n}}$. In particular, this allows to consider the action of the fractional Laplacian over this set.

Finally, Sobolev spaces of negative order are defined by duality, using $L^{2}(\Omega)$ as pivot space. Of interest in the problems we are considering is the space

[TABLE]

Recalling definition (2.1) and denoting by $\left(\cdot,\cdot\right)$ the inner product in $L^{2}(\Omega)$, the weak formulation of problem (1.1) reads: find $u\in L^{2}(0,T;\widetilde{H}^{s}(\Omega))$ such that

[TABLE]

for all $w\in L^{2}(0,T;\widetilde{H}^{s}(\Omega)).$

2.2. Elliptic regularity

We recall regularity results for the homogeneous problem

[TABLE]

Even though the fractional Laplacian is an operator of order $2s$ in ${\mathbb{R}^{n}}$, in the sense that $(-\Delta)^{s}\colon H^{\ell}({\mathbb{R}^{n}})\to H^{\ell-2s}({\mathbb{R}^{n}})$ is bounded and invertible, theory is much more delicate for problems posed on bounded domains.

Grubb [20] provides regularity estimates for solutions of (2.2) in the setting of Hörmander $\mu-$spaces. We express these results in terms of standard Sobolev spaces, and refer to [20] for further details.

Proposition 2.2**.**

Let $\Omega\subset{\mathbb{R}^{n}}$ be a bounded domain with smooth boundary, $g\in H^{r}(\Omega)$ for some $r\geq-s$ and consider $u\in\widetilde{H}^{s}(\Omega)$, the solution of the Dirichlet problem (2.2). Then, there exists a constant $C(n,s)$ such that

[TABLE]

Here, $\gamma=\min\{s,1/2-\varepsilon\}$, with $\varepsilon>0$ arbitrary small.

Remark 2.3*.*

Observe that, in general, it is not true that solutions of (2.2) have $2s$ more derivatives than the right-hand side function $g$. No matter how regular $g$ is, the solution of (2.2) is not expected to be more regular than $H^{s+\gamma}({\mathbb{R}^{n}})$. In spite of this, the singular behavior of solutions can be localized at the boundary and described more appropriately in weighted spaces [3, 4].

Remark 2.4*.*

In view of Proposition 2.2, given $v\in\widetilde{H}^{s}(\Omega)$, assuming that it satisfies $(-\Delta)^{s}v\in L^{2}(\Omega)$ is weaker than assuming that $v\in H^{2s}({\mathbb{R}^{n}})$. Hypotheses of this type on the initial/boundary conditions are employed throughout this work.

2.3. Mittag-Leffler function

Let $\alpha>0$ and $\mu\in{\mathbb{R}}$, then the Mittag-Leffler function $E_{\alpha,\mu}\colon\mathbb{C}\to\mathbb{C}$ is defined by

[TABLE]

This is a complex function that depends on two parameters; in particular, it generalizes the exponentials, in view of the identity $E_{1,1}(z)=e^{z}$ for all $z\in\mathbb{C}$. The following properties of this family of functions are useful to derive the regularity estimates we present in Section 3.

Lemma 2.5** (cf. [16, Theorem 4.3]).**

If $\alpha,\lambda>0$, then

[TABLE]

Moreover, the following identities hold for integer-order differentiation:

[TABLE]

and

[TABLE]

Lemma 2.6** (cf. [43, Theorem 1.6]).**

Let $0<\alpha<2$ and $\mu\in{\mathbb{R}}$ satisfying $\pi\alpha/2<\delta<\min\{\pi,\pi\gamma\},$ then there exists a constant $C(\alpha,\mu,\delta)>0$ such that

[TABLE]

2.4. Solution representation

Let $\{(\phi_{k},\lambda_{k})\}_{k=1}^{\infty}$ denote the solutions of the fractional eigenvalue problem

[TABLE]

It is well-known that the fractional Laplacian has a sequence of eigenvalues

[TABLE]

and that the eigenfunctions’ set $\{\phi_{k}\}_{k=1}^{\infty}$ may be taken to constitute an orthonormal basis of $L^{2}(\Omega)$.

Remark 2.7*.*

Unlike the classical Laplacian, eigenfunctions of the fractional Laplacian are in general non-smooth [21, 44]. Indeed, considering a smooth function $d$ that behaves like $\text{dist}(x,\partial\Omega)$ near to $\partial\Omega$, all eigenfunctions $\phi_{k}$ belong to the space $d^{s}C^{2s(-\varepsilon)}(\Omega)$ (the $\varepsilon$ is active only if $s=1/2$) and $\frac{\phi_{k}}{d^{s}}$ does not vanish near $\partial\Omega$. The best Sobolev regularity guaranteed for solutions of (2.3) is $\phi_{k}\in H^{s+1/2-\varepsilon}({\mathbb{R}^{n}})$ for $\varepsilon>0$ (see [12]).

The reduced Sobolev regularity of eigenfunctions precludes the possibility of solutions to equation (1.1) being smooth, even for $\alpha=1$. This is in stark contrast with the case of the classical Laplacian. However, solutions of diffusion equations with memory –local in space but fractional in time– are known to be less regular than their classical counterparts [35]. The effect of fractional differentiation in time is that high-frequency modes are less strongly damped than in classical diffusion, and the time derivatives of the solution are unbounded as $t\to 0$. In the next section we shall show that solutions of (1.1) present the same behavior. See also the numerical experiments in Subsection 6.3.

We write solutions of (1.1) by means of separation of variables,

[TABLE]

Then, for every $k\geq 1$ it must hold that

[TABLE]

where $f_{k}=\left(f,\phi_{k}\right)$, $v_{k}=\left(v,\phi_{k}\right),$ and $b_{k}=\left(b,\phi_{k}\right).$ Existence and uniqueness of solutions to (2.5) follow from standard theory for fractional-order differential equations [16, Theorem 7.2]. Moreover, solutions of (2.5) may be represented as the superposition of the respective solution of the problem with forcing term and initial data equal to zero and the solution of the problem with vanishing forcing term. Namely, defining

[TABLE]

the solution of (2.5) may be written as

[TABLE]

For the particular value of $\alpha=1$ the above expression yields the well-known formula

[TABLE]

usually derived by the method of variation of parameters. Also, considering $\alpha=2$, in virtue of the identities $E_{2,1}(z)=\cosh(\sqrt{z})$ and $E_{2,2}(z)=\frac{\sinh(\sqrt{z})}{\sqrt{z}}$, expression (2.7) becomes

[TABLE]

Summing the solutions for every eigenmode, we obtain the following result.

Theorem 2.8**.**

Let $\Omega$ be a bounded, smooth domain, $s\in(0,1)$ and $\alpha\in(0,2]$. Assume that $f\in L^{\infty}(0,T;L^{2}(\Omega))$, $v\in L^{2}(\Omega)$ and $b\in L^{2}(\Omega)$ are given. Then, problem (1.1), with boundary conditions (1.3) (and (1.4) if $\alpha\in(1,2]$) admits a solution $u\in L^{2}(0,T;\widetilde{H}^{s}(\Omega))$ that can be represented by (2.4). The modes $u_{k}$ are given, accordingly, by (2.6) if $\alpha\in(0,1]$ and (2.7) if $\alpha\in(1,2]$.

3. Regularity of solutions

In this section we state some regularity results for solutions of the problems under consideration. We split the estimates according to whether the initial values or the forcing term are null. For the sake of brevity we omit the proofs; these can be obtained following the arguments of [45], and are based on expansion (2.4) and utilizing Lemmas 2.5 and 2.6 to bound the corresponding coefficients. We also recall also that throughout this paper we are assuming that $\Omega$ is a domain with smooth boundary, so that Proposition 2.2 holds. According to that proposition, we fix the notation $\gamma:=\min\{s,1/2-\varepsilon\}$, with $\varepsilon>0$ arbitrarily small.

Theorem 3.1**.**

Let $0<\alpha\leq 1$ and suppose that $f\equiv 0$. Let $u$, given by (2.4), be the solution of (1.1) with initial and boundary conditions according to (1.3).

1. (a)

If $v\in L^{2}(\Omega)$, then $u\in C([0,T];L^{2}(\Omega))\cap C((0,T];\widetilde{H}^{s}(\Omega)\cap H^{s+\gamma}(\Omega))$ and $\leftidx{{}^{C}}\partial_{t}^{\alpha}u\in C((0,T];L^{2}(\Omega))$. Moreover, there exists a constant $C>0$ such that

[TABLE] 2. (b)

Assume that $v\in\widetilde{H}^{s}(\Omega)$. Then, $u\in L^{2}(0,T;\widetilde{H}^{s}(\Omega)\cap H^{s+\gamma}(\Omega))$, $\leftidx{{}^{C}}\partial_{t}^{\alpha}u\in L^{2}(\Omega\times(0,T))$, and the following estimate holds:

[TABLE] 3. (c)

Furthermore, if $v\in\widetilde{H}^{s}(\Omega)$ is such that $(-\Delta)^{s}v\in L^{2}(\Omega)$, then $u\in C([0,T];\widetilde{H}^{s}(\Omega)\cap H^{s+\gamma}(\Omega))$, $\leftidx{{}^{C}}\partial_{t}^{\alpha}u\in C([0,T];L^{2}(\Omega))$ and the bound

[TABLE]

is satisfied for some $C>0$ independent of $v$.

Regularity estimates for the fractional diffusion problem with a non-homogeneous right-hand side function $f$ are also attainable.

Theorem 3.2**.**

Let $0<\alpha\leq 1$ and $v\equiv 0$. Consider $u$, given by (2.4), the solution of (1.1) with homogeneous initial and boundary conditions. If $f\in L^{\infty}(0,T;L^{2}(\Omega))$, then $u\in L^{2}(0,T;\widetilde{H}^{s}(\Omega)\cap H^{s+\gamma}(\Omega))$, $\leftidx{{}^{C}}\partial_{t}^{\alpha}u\in L^{2}(\Omega\times(0,T))$ and

[TABLE]

Estimates for the fractional diffusion-wave case are obtained similarly.

Theorem 3.3**.**

Let $1<\alpha\leq 2$ and suppose that $f\equiv 0$. Let $u$, given by (2.4), be the solution of (1.1) with initial/boundary conditions (1.3) and (1.4).

1. (a)

Assume that $v\in L^{2}(\Omega)$ and $b\in L^{2}(\Omega)$. Then, $u\in C([0,T];L^{2}(\Omega))\cap C((0,T];\widetilde{H}^{s}(\Omega)\cap H^{s+\gamma}(\Omega))$ and $\leftidx{{}^{C}}\partial_{t}^{\alpha}u\in C((0,T];L^{2}(\Omega))$. Moreover, there exists a constant $C>0$ such that

[TABLE] 2. (b)

If $v\in\widetilde{H}^{s}(\Omega)$ and $b\in L^{2}(\Omega)$, then $\partial_{t}u\in C([0,T];H^{-s}(\Omega))$, and

[TABLE] 3. (c)

Moreover, if $v\in\widetilde{H}^{s}(\Omega)$ is such that $(-\Delta)^{s}v\in L^{2}(\Omega)$ and $b\in\widetilde{H}^{s}(\Omega)$, then $u\in C([0,T];\widetilde{H}^{s}(\Omega)\cap H^{s+\gamma}(\Omega))\cap C^{1}([0,T];L^{2}(\Omega))$, $\leftidx{{}^{C}}\partial_{t}^{\alpha}u\in C([0,T];L^{2}(\Omega))$, and the following estimates hold:

[TABLE]

Finally, estimates for problems with non-null forcing term and $\alpha\in(1,2]$ have the following form.

Theorem 3.4**.**

Let $1<\alpha\leq 2$, $v\equiv 0$ and $b\equiv 0$. Consider $u$, given by (2.4), be the solution of (1.1) with homogeneous initial and boundary conditions. If $f\in C([0,T];L^{2}(\Omega))$ is such that $(-\Delta)^{s}f\in L^{2}(\Omega\times(0,T))$, then $u\in C([0,T];\widetilde{H}^{s}(\Omega)\cap H^{s+\gamma}(\Omega))$, $\leftidx{{}^{C}}\partial_{t}^{\alpha}u\in C([0,T];L^{2}(\Omega))$ and

[TABLE]

4. Numerical approximations

In this section we devise a discrete scheme to approximate (1.1). To this end, standard Galerkin finite elements are utilized in the spatial discretization (following [3]) and a convolution quadrature is used for the time variable (following [26]).

4.1. Semi discrete scheme

For an appropriate treatment, it is convenient to derive the numerical scheme in two steps. In first place we discretize in space, and afterwards in time. We follow the ideas developed in [26], taking advantage of the fact that, from the theoretical point of view, minor changes are required to handle with the fractional Laplacian instead of its classical counterpart.

Indeed, let $\mathcal{T}_{h}$ be a shape regular and quasi-uniform admissible simplicial mesh on $\Omega$, and let $X_{h}\subset\widetilde{H}^{s}(\Omega)$ be the piecewise linear finite element space associated with $\mathcal{T}_{h}$, namely,

[TABLE]

The semidiscrete problem reads: find $u_{h}\colon[0,T]\to X_{h}$ such that

[TABLE]

Here, $v_{h}=P_{h}v$, $b_{h}=P_{h}b$, and $P_{h}$ denotes the $L^{2}(\Omega)$ projection on $X_{h}$.

Observe that, defining the discrete fractional Laplacian $A_{h}:X_{h}\rightarrow X_{h}$ as the unique operator that satisfies

[TABLE]

and considering $f_{h}:=P_{h}f$, we may rewrite (4.1) as

[TABLE]

4.2. Fully discrete scheme

At this point, a suitable discretization of the Caputo differentiation operator is required to obtain a fully discrete scheme. To this end, we employ the convolution quadrature technique described by Lubich in [29, 30], which allows us to derive discrete estimations of an integral which involve singular kernels.

Upon dividing $[0,T]$ uniformly with a time step size $\tau=T/N$, and letting $t=n\tau$ ($n\in\{1,\ldots,N\}$), by means of the convolution quadrature rule we are able to estimate the Riemann-Liouville operator of a function $g$ by

[TABLE]

where the weights $\{\omega_{j}\}_{j\in\mathbb{N}_{0}}$ are obtained as the coefficients of the power series

[TABLE]

For the reader’s convenience we give an overview of the main ideas of this technique in Appendix A and refer the reader, for example, to [29, 30] for further details.

We are now able to suitably discretize the Caputo differentiation operator. To this end, we need to reformulate (1.1) using the Riemann-Liouville derivative instead of the Caputo one. It is well-known that these two operators are related by (see, for example, [16, Theorem 3.1])

[TABLE]

Thus, we rewrite (4.2) for the fractional diffusion case as

[TABLE]

and for the fractional diffusion-wave case as

[TABLE]

Replacing the Riemann-Liouville derivative by its discrete version given by (4.3), and that we will denote by $\overline{\partial_{\tau}}^{\alpha}$, we formulate the fully discrete problem as: find $U^{n}_{h}\in X_{h}$, with $n=\{1,\ldots,N\}$, such that

[TABLE]

or

[TABLE]

for fractional diffusion and fractional diffusion-wave problems respectively, where $F^{n}_{h}=P_{h}f(t_{n})$.

In order to obtain a better error estimation in the diffusion-wave case, it is necessary to replace $F^{n}_{h}$ with a corrected term $G^{n}_{h}:=\overline{\partial_{\tau}}\partial^{-1}_{t}f_{h}(t_{n})$. See [26] for further details.

For the sake of the reader’s convenience, we conclude this section by giving the vectorial form of the fully discrete scheme. Let $\{\varphi_{i}\}_{i=1,\ldots,\mathcal{N}}$ be the Lagrange nodal basis that generates $X_{h}$. Let $U^{n}$, $F^{n}$ and $G^{n}\in\mathbb{R}^{\mathcal{N}}$, $n=0,\ldots,N$ be such that $U^{n}_{h}=\sum_{i=1}^{\mathcal{N}}U^{n}_{i}\varphi_{i}$, $F^{n}_{h}=\sum_{i=1}^{\mathcal{N}}F^{n}_{i}\varphi_{i}$ and $G^{n}_{h}=\sum_{i=1}^{\mathcal{N}}G^{n}_{i}\varphi_{i}$, where $U_{h}^{n}$ denotes the solution of the fully discrete problem. Then we formulate (4.4) and (4.5), respectively, in the following vectorial equations:

[TABLE]

and

[TABLE]

Above, $M,K\in\mathbb{R}^{\mathcal{N}\times\mathcal{N}}$ are the mass and stiffness matrices, respectively. Namely, $M_{i,j}=(\varphi_{i},\varphi_{j})$ and $K_{i,j}=\langle\varphi_{i},\varphi_{j}\rangle_{s}$.

The computation and assembly of the stiffness matrix in dimension greater than one is not a trivial task. However, this problem for two-dimensional domains has been tackled in [2], where the authors provide a short MATLAB implementation.

There are several options to compute the coefficients $\{\omega_{j}\}_{j\in\mathbb{N}_{0}}$. Recalling that

[TABLE]

Fast Fourier Transform can be used for an efficient computation of $\{\omega_{j}\}_{j\in\mathbb{N}_{0}}$ (see [43, Section 7.5]) . Alternatively, a useful recursive expression is given also in [43, formula (7.23)]:

[TABLE]

For the numerical experiments we exhibit in Section 6 we have taken advantage of this identity.

5. Error bounds

This section shows error estimates for the numerical scheme discussed in Section 4. The derivation of the error bounds can be carried out following the guidelines from [26]. Most of these results can be extrapolated to our case. Details of the proofs in those cases where the generalization is not direct are provided in Appendix B.

5.1. Error bounds for the semi-discrete scheme

Here we focus only on the diffusion-wave case ($1<\alpha<2$), where the generalization of the error bounds becomes more laborious. In this context, the error bounds are given by the generalization of two theorems. The first one, in the spirit of [26, Theorem 3.2], provides an error estimation for the homogeneous case.

Theorem 5.1**.**

Let $1<\alpha<2$, $u$ be the solution of (1.1) with $v\in\widetilde{H}^{q}(\Omega)$, $b\in\widetilde{H}^{r}(\Omega)$ $q,r\in[0,2s]$, and $f=0$; and let $u_{h}$ be the solution of (4.1) with $v_{h}=P_{h}v$, $b_{h}=P_{h}b$, and $f_{h}=0$. Writing $e_{h}(t)=u(t)-u_{h}(t)$, there exists a positive constant $C=C(s,n)$ such that

[TABLE]

To complete the error estimate for the semi-discrete scheme, it still remains to analyze the case $v=0,b=0$ and $f\not=0$. A proper generalization of [24, Theorem 3.2] can be carried out following the guidelines outlined in that work.

Theorem 5.2**.**

Let $1<\alpha<2$, $f\in L^{\infty}(0,T;L^{2}(\Omega))$, and let $u$ and $u_{h}$ be the solutions of (1.1) and (4.1) respectively, with $f_{h}=P_{h}f$, and all the initial data equal to zero. Then, there exists a positive constant $C=C(s,n)$ such that

[TABLE]

5.2. Error bounds for the fully-discrete scheme

Considering all the theory displayed up to this point, error estimates for the fully-discrete scheme can be derived in the same way as in [26]. We refer the reader to that work for the details.

Theorem 5.3**.**

Let $u$ be the solution of problem (1.1) with $v\in\widetilde{H}^{q}(\Omega)$, $b\in\widetilde{H}^{r}(\Omega)$ $q,r\in[0,2s]$, and $f=0$; and let $U^{n}_{h}$ be the solution of (4.4) or (4.5) with $v_{h}=P_{h}v$, $b_{h}=P_{h}b$, and $F^{n}_{h}=0$. Then, there exists a positive constant $C=C(s,n)$ such that

• •

If $0<\alpha<1,$ then

[TABLE]

• •

If $1<\alpha<2,$ then

[TABLE]

Remark 5.4*.*

In the previous theorem –and in Theorem 5.1 as well – we wrote the orders of convergence in term of various Sobolev norms of the data. For clarity, hypotheses in theorems 3.1 and 3.3 just involved either $L^{2}$ or $H^{s}$ norms of the data. For instance, assuming that $v\in\widetilde{H}^{s}(\Omega)$ is such that $(-\Delta)^{s}v\in L^{2}(\Omega)$ and $b\in\widetilde{H}^{s}(\Omega)$, the conclusions of Theorem 5.3 read

[TABLE]

We emphasize that, as stated in Remark 2.4, the identity $\|(-\Delta)^{s}v\|_{L^{2}(\Omega)}\leq C|v|_{H^{2s}({\mathbb{R}^{n}})}$ holds for all $v\in\widetilde{H}^{2s}(\Omega)$.

Finally, we state the order of convergence of the fully-discrete scheme for the problems with a non-null source term.

Theorem 5.5**.**

Let $u$ be the solution of (1.1) with homogeneous initial data and with $f\in L^{\infty}(0,T;L^{2}(\Omega))$; and let $U^{n}_{h}$ be the solution of (4.4) or (4.5) with $f_{h}=P_{h}f$. Then, there exists a positive constant $C=C(s,n)$ such that

• •

For $0<\alpha<1$, if $\int^{t}_{0}(t-s)^{\alpha-1}\|f^{\prime}(s)\|_{L^{2}(\Omega)}ds<\infty$ for $t\in(0,T]$, then

[TABLE]

• •

If $1<\alpha<2,$ then

[TABLE]

6. Numerical experiments

This section exhibits the results of numerical tests for discretizations of problems posed in one- and two-dimensional domains. Numerical solutions of (1.1) were obtained by applying the scheme described in Section 4. The experiments in two-dimensional geometries were carried out with a code based on the one presented in [2].

6.1. Explicit Solutions

In [2] it is shown how some families of non-trivial solutions for the fractional Poisson problem can be constructed. For the sake of brevity, we refer the reader to that work for details. Here we summarize these results in order to be applied to the evolution equation in the cases in which $\Omega$ corresponds to a) $(-1,1)\subset{\mathbb{R}}$ and, more generally, b) $B(0,1)\subset{\mathbb{R}}^{n}$.

Define for $n\geq 1$, the function $\omega^{s}:{\mathbb{R}}^{n}\to{\mathbb{R}},$

[TABLE]

Then,

[TABLE]

solves

[TABLE]

with $f(x)=\mu_{s}^{k}g_{k}^{(s)}(x)$, where in case a)

[TABLE]

and in case b)

[TABLE]

Above, $C_{k}^{(s+1/2)}$ and $P_{k}^{(s,\,n/2-1)}$ denote a Gegenbauer and a Jacobi polynomial [1], respectively.

Next, let $h(t)$ be a function such that $\leftidx{{}^{C}}\partial_{t}^{\alpha}h(t)$ can be easily computed. By means of separation of variables we can construct explicit solutions of the fractional evolution problem of the form

[TABLE]

6.2. Orders of convergence

In order to confirm the predicted convergence rate, we show the results we obtained in three example problems:

1. (a)

$u(x,t)=E_{\alpha,1}(-t^{\alpha})\cdot\omega^{s}(x)C^{(s)}_{3}(x)$, $\Omega=(-1,1)$; 2. (b)

$u(x,t)=\sin(t)\cdot\omega^{s}(x)C^{(s)}_{3}(x)$, $\Omega=(-1,1)$; 3. (c)

$u(x,t)=E_{\alpha,1}(-t^{\alpha})\cdot\omega^{s}(x)P^{(s,0)}_{k}(2|x|^{2}-1)$, $\Omega=B(0,1)\subset{\mathbb{R}}^{2}.$

For examples (a) and (b) we examine the time and spatial convergence over a fixed time $t=0.1$. A fixed small time step is taking to see the spatial convergence and viceversa. Our results are summarized in tables 1, 2, 3 and 4.

The experimental orders of convergence (e.o.c.) are in agreement with the theory in the case $s>1/2$ while our numerical examples exhibit e.o.c. (in space) higher than those predicted if $s<1/2$ (see Tables 3 and 4). This behavior seems to be due to the fact that the extra regularity of the data present in our examples can not be exploited in our arguments; the actual solutions are more regular than what is predicted by Theorems 3.1 and 3.3.

6.3. Qualitative aspects

Finally, we present experiments that illustrate some qualitative effects of the fractional derivatives. In Figure 1, we fix $s=0.5$, and show the evolution in time for different values of the parameter $\alpha$, ranging from fractional diffusion to fractional diffusion-wave. Memory effects are present for $\alpha=0.5$, while the solution oscillates for $\alpha=1.5$.

Figure 2, in turn, displays the effect of moving the parameters $\alpha$ and $s$ for a fixed time. It can be seen that increasing the spatial differentiation order $s$ leads to a faster spreading of the initial condition. Apparent differences can be noticed among the three different problems with $\alpha=2s$ exhibited along the diagonal of the figure.

Our last example, in Figure 3, exhibits the persistence of a singularity along the time due to the memory induced by the fractional in-time derivative. In that experiment, we have set $\alpha=0.99$ and $s=0.9$. Notice that the solution vanishes to [math] as time increases, but even though the differentiation parameters are both close to $1$, which corresponds to the classical heat equation, the singular behavior of the initial condition persists in time.

Acknowledgments

The authors thank Prof. Bangti Jin for pointing out reference [26].

Appendix A Convolution Quadrature Rule

The aim of this appendix is to describe a numerical approximation technique for convolutions that plays an important role in the assemblage of the numerical scheme we propose.

Dividing $[0,T]$ uniformly with a time step size $\tau=T/N$, and letting $t=n\tau$ ($n\in\{1,\ldots,N\}$), we seek for a numerical approximation of the convolution integral

[TABLE]

by means of a finite sum $\sum_{j=0}^{n}\omega_{j}g(t-j\tau).$ The weights $\{\omega_{j}\}_{j\in\mathbb{N}_{0}}$ are obtained as the coefficients of the power series

[TABLE]

where $K$ denotes the Laplace transform of the kernel $k$, and $\delta(\xi)$ is the quotient of the generating polynomials of a linear multistep method.

To obtain the weights $\omega_{j}$, suppose that we extend the kernel $k$ by zero over $r\leq 0$ and that for all $r>0$ it satisfes

[TABLE]

for some $c,\mu>0$. Then, the inversion formula

[TABLE]

holds, where $\Gamma$ is a contour lying in the sector of analyticity of $K$, parallel to its boundary and oriented with an increasing imaginary part. Furthermore, defining $\Sigma_{\theta}:=\{z\in\mathbb{C}:|\text{arg}(z)|\leq\theta\}$, $\theta\in(\pi/2,\pi)$, it holds that $K$ is analytic in $\Sigma_{\theta}$ and satisfies

[TABLE]

This condition is in turn equivalent to (A.2).

Replacing (A.3) in (A.1) and switching the order of integration gives

[TABLE]

Since the inner integral in the right-hand side is the solution of the ordinary differential equation $y^{\prime}=zy+g$, with $y(0)=0$, we can obtain a numerical estimation by using some multistep method. For simplicity, suppose we utilize Backward Euler discretization (BE), that gives the scheme

[TABLE]

Multiplying by $\xi^{n}$ both sides of the equality, and summing over $n$, we obtain

[TABLE]

where $\mathbf{y}(\xi):=\sum^{\infty}_{n=0}y_{n}\xi^{n}$ and $\mathbf{g}(\xi):=\sum^{\infty}_{n=0}g_{n}\xi^{n}$. Defining $\delta(\xi):=(1-\xi)$, from (A.6) we deduce

[TABLE]

Thus, the numerical approximation of $y$ at time $n\tau$ is given by the $n$-th coefficient of the power series $\left(\frac{\delta(\xi)}{\tau}-z\right)^{-1}\mathbf{g}(\xi)$.

In order to obtain the desired numerical approximation of (A.1) we utilize the former expression, fix $\xi$ and integrate in $z$ the right hand side in (A.5). Using Cauchy’s integral formula gives

[TABLE]

Therefore, the numerical approximation of (A.1) at $t=n\tau$ is given by the $n$-th coefficient of the power series $K\left(\frac{\delta(\xi)}{\tau}\right)\cdot\mathbf{g}(\xi)$. Finally, noticing that the coefficients of the series are the Cauchy product of the sequences $\{\omega_{n}\}_{n\in\mathbb{N}_{0}}$ and $\{g(n\tau)\}_{n\in\mathbb{N}_{0}}$, where $\{\omega_{n}\}$ are the coefficients of the power series expansion of $K\left(\frac{\delta(\xi)}{\tau}\right)$, we obtain an expression for the weights in (4.3).

Given a complex valued function $K$, analytic in $\Sigma_{\theta}$ and satisfying (A.4), we use the transfer function notation for (A.1),

[TABLE]

where $k$ is given by (A.3), and the notation

[TABLE]

for the discrete approximation.

Next, we generalize the definition of the Convolution Quadrature Rule to operators that satisfy (A.4) with a negative value of $\mu$. Indeed, let $m$ be a positive integer such that $\mu+m>0$, setting $\tilde{K}(z):=z^{-m}K(z)$ we define

[TABLE]

with $\tilde{k}$ the kernel associated with $\tilde{K}$. All the results and estimates that are achieved in the former case are still true upon this generalization (see [29, Section 5]). This is convenient because we are interested in the particular case of $K(z)=z^{\alpha}$, that delivers

[TABLE]

with $m$ a positive integer such that $m-1\leq\alpha<m$. Considering this, we set the notation $K(\partial_{t}):=K(z)$ and $K(\overline{\partial_{\tau}}):=K(\delta(\xi)/\tau)$.

Appendix B Error analysis of the semi-discrete scheme

We first set an integral representation of $u$ for the homogeneous case $f=0$. Define the sector $\Sigma_{\theta}:=\{z\in\mathbb{C}\text{ such that }|\arg(u)|<\theta\}$, then $u(t):[0,T]\rightarrow L^{2}(\Omega)$ can be analytically extended to $\Sigma_{\pi/2}$ (see [45, Theorem 2.3]). Applying the Laplace transform in (1.1) we obtain

[TABLE]

where $A$ is the fractional Laplacian with homogeneous Dirichlet conditions. Therefore, via the Laplace inversion formula, we write the integral representation

[TABLE]

where $\Gamma_{\theta,\delta}=\{z\in\mathbb{C}:|z|=\delta,|\arg(z)|\leq\theta\}\cup\{z\in\mathbb{C}:z=re^{\pm i\theta},r\geq\delta\}$.

If we choose $\theta$ such that $\pi/2<\theta<\min\{\pi,\pi/\alpha\}$, then $z^{\alpha}\in\Sigma_{\theta^{\prime}}$ with $\theta^{\prime}=\alpha\theta$ for all $z\in\Sigma_{\theta}.$ Considering $\theta$ in this mode, there exists a constant $C$ only depending on $\theta$ and $\alpha$ such that

[TABLE]

As in (B.1), we can write an analogous expression for $u_{h}$,

[TABLE]

The following technical result can be proved analogously to [19, Lemma 7.1].

Lemma B.1**.**

Let $\varphi\in\widetilde{H}^{s}(\Omega)$, and $z\in\Sigma_{\theta}$ with $\pi/2<\theta<\min\{\pi,\pi/\alpha\}$. Then there exists a positive constant $c(\theta)$ such that

[TABLE]

The next lemma sets an error estimate between $(z^{\alpha}I+A)^{-1}f$ and its discrete approximation $(z^{\alpha}I+A_{h})^{-1}P_{h}f$, analogous to [7, Lemma 3.4].

Lemma B.2**.**

Let $f\in L^{2}(\Omega)$, $z\in\Sigma_{\theta}$, $\omega:=(z^{\alpha}I+A)^{-1}f$, $\omega_{h}:=(z^{\alpha}I+A_{h})^{-1}P_{h}f$. Then there exists a positive constant $C(s,n,\theta)$ such that

[TABLE]

As before, $\gamma=\min\{s,1/2-\varepsilon\}$, with $\varepsilon>0$ arbitrary small.

Proof.

We consider first the case $s\geq 1/2$. By definition of $\omega$ and $\omega_{h}$, it holds that

[TABLE]

If we set $e_{h}:=\omega-\omega_{h}$ and subtract these two expressions, we derive

[TABLE]

Applying Lemma B.1 and this identity, we arrive to

[TABLE]

Taking $\varphi=\Pi_{h}\omega$ in the former expression, where $\Pi_{h}$ is a suitable quasi-interpolation operator (see, for example, [3, Section 4.1]), we deduce

[TABLE]

where we have used the fact that $s\geq 1/2$ implies that $h^{s}\leq h^{1/2-\varepsilon}$.

On the other hand, if we choose $\varphi=\omega$ in Lemma B.1, we obtain

[TABLE]

Consequently,

[TABLE]

From Proposition 2.2, we know that for the case $z=0$ the estimate $|\omega|_{H^{s+1/2-\varepsilon}({\mathbb{R}^{n}})}\leq\|f\|_{L^{2}(\Omega)}$ holds. Utilizing this estimate with $-z^{\alpha}\omega+f$ instead of $f$, we obtain

[TABLE]

where in the last inequality we used (B.4). Combining this with (B.3), we derive

[TABLE]

This implies that

[TABLE]

and gives the bound for $|e_{h}|_{H^{s}({\mathbb{R}^{n}})}$. Next, we aim to estimate $\|e_{h}\|^{2}_{L^{2}(\Omega)}$. For this purpose, we proceed via the following duality argument. Given $\varphi\in L^{2}(\Omega)$, define

[TABLE]

Thus, we write

[TABLE]

We aim to bound the supremum in the identity above. Resorting to (B.2) and the Cauchy-Schwarz inequality, we bound

[TABLE]

Finally, applying (B.5) we arrive at

[TABLE]

from where we can derive the desired inequality.

The analysis of the case $s\leq 1/2$ can be carried out in analogously. Indeed, using that $h^{s}\geq h^{1/2-\varepsilon}$ we obtain, instead of (B.3), the inequality

[TABLE]

and proceeding as before we arrive at the desired estimate. ∎

At this point, we are able to give a sketch of the proof of Theorem 5.1.

Proof of Theorem 5.1.

The proof outlined in [26, Theorem 3.2] can be reproduced using Lemma B.2 instead of Lemma 3.1 from that work and approximation properties of the elliptic projection along with the estimate $\|(-\Delta)^{s}w\|_{L^{2}(\Omega)}\leq C|w|_{H^{2s}({\mathbb{R}^{n}})}.$ The details are therefore omitted. ∎

Finally, we consider the problem with zero initial conditions and a non-zero source term.

Proof of Theorem 5.2.

Using the approximation properties of the elliptic projection together with an appropriate generalization of [24, Lemma 2.3] (in the same spirit of Lemma B.2), the proof given in that work can be adapted to our purpose without major difficulties. ∎