Numerical Approximation of Space-time Fractional Parabolic Equations

TL;DR

This paper introduces a novel numerical scheme for space-time fractional parabolic equations, utilizing Dunford-Taylor integrals, finite element methods, and sinc quadrature to achieve optimal space discretization and efficient fully discrete solutions.

Contribution

It develops a new numerical approach combining finite element approximation with Dunford-Taylor integrals and sinc quadrature for space-time fractional parabolic equations, avoiding traditional time-stepping.

Findings

Optimal order space error up to logarithm of 1/h

Fully discrete method free of time stepping

Effective approximation of convolution via sinc quadrature

Abstract

In this paper, we develop a numerical scheme for the space-time fractional parabolic equation, i.e., an equation involving a fractional time derivative and a fractional spatial operator. Both the initial value problem and the non-homogeneous forcing problem (with zero initial data) are considered. The solution operator $E(t)$ for the initial value problem can be written as a Dunford-Taylor integral involving the Mittag-Leffler function $e_{\alpha,1}$ and the resolvent of the underlying (non-fractional) spatial operator over an appropriate integration path in the complex plane. Here $\alpha$ denotes the order of the fractional time derivative. The solution for the non-homogeneous problem can be written as a convolution involving an operator $W(t)$ and the forcing function $F(t)$. We develop and analyze semi-discrete methods based on finite element approximation to the underlying (non-fractional) spatial operator in terms of analogous Dunford-Taylor integrals applied to the discrete operator. The space error is of optimal order up to a logarithm of $1/h$. The fully discrete method for the initial value problem is developed from the semi-discrete approximation by applying an exponentially convergent sinc quadrature technique to approximate the Dunford-Taylor integral of the discrete operator and is free of any time stepping. To approximate the convolution appearing in the semi-discrete approximation to the non-homogeneous problem, we apply a pseudo midpoint quadrature. This involves the average of $W_h(s)$, (the semi-discrete approximation to $W(s)$) over the quadrature interval. This average can also be written as a Dunford-Taylor integral. We first analyze the error between this quadrature and the semi-discrete approximation. To develop a fully discrete method, we then introduce sinc quadrature approximations to the Dunford-Taylor integrals for computing the averages.

Equations

$$\{\begin{array}{rl}{\partial }_t^{\gamma }u+L^{\beta }u& =f,\text{in }(0,\mathsf{T}]\times \Omega ,\\ u& =0,\text{on }(0,\mathsf{T}]\times \partial \Omega ,\\ u& =v,\text{on }\{0\}\times \Omega .\end{array}$$

$${\partial }_t^{\gamma }u(t):=\frac{1}{\Gamma (1-\gamma )}{\int }_0^t\frac{1}{(t-r{)}^{\gamma }}\frac{\partial u(r)}{\partial r}dr.$$

$${\partial }_t^{\gamma }u={}^R{\partial }_t^{\gamma }(u-u(0))$$

$$L^{\beta }v:=\sum _{j=1}^{\infty }{\lambda }_j^{\beta }(v,{\psi }_j){\psi }_j,$$

$$\{\begin{array}{rl}⟨{\partial }_t^{\gamma }u,\varphi ⟩+A(u,\varphi )& =(f,\varphi ),\text{for all }\varphi \in D(L^{\beta /2})\text{ and for a.e. }t\in (0,\mathsf{T}],\\ u(0)& =v.\end{array}$$

$$u(t):=u(t,\cdot )=E(t)v+{\int }_0^{t}W(r)f(t-r)dr.$$

$$E(t)w:=e_{\gamma ,1}(-t^{\gamma }{L}^{\beta })w=\sum _{j=1}^{\infty }{e}_{\gamma ,1}(-t^{\gamma }{\lambda }_j^{\beta })(w,{\psi }_j){\psi }_j$$

$$W(t)w:=t^{\gamma -1}{e}_{\gamma ,\gamma }(-t^{\gamma }{L}^{\beta })w=\sum _{j=1}^{\infty }{t}^{\gamma -1}{e}_{\gamma ,\gamma }(-t^{\gamma }{\lambda }_j^{\beta })(w,{\psi }_j){\psi }_j,$$

$$\int_{t_{j-1}}^{t_{j}}W_{h}({r}){\pi_{h}}f(t-{r})\,d{r}\approx\bigg{(}\int_{t_{j-1}}^{t_{j}}W_{h}({r})\,d{r}\bigg{)}\ {\pi_{h}}f(t-{t_{j-\frac{1}{2}}}),$$

$$\|u\|:=\|u\|_{L^{2}}:=\bigg{(}\int_{\Omega}|u|^{2}\,dx\bigg{)}^{1/2}\quad\text{and}\quad\|u\|_{H^{1}}:=(\|u\|^{2}_{L^{2}}+\||\nabla u|\|_{L^{2}}^{2})^{1/2}.$$

$$(v,w):={\int }_{\Omega }v(x)\overline{w}(x)dx.$$

$$\Vert F{\Vert }_{H^{-1}}:=\underset{\theta \in H_0^1}{\sup }\frac{⟨F,\theta ⟩}{|\theta {|}_{H^1}},$$

$$\Vert A{\Vert }_{B_1\to B_2}=\underset{v\in B_1,v\ne 0}{\sup }\frac{\Vert Av{\Vert }_{B_2}}{\Vert v{\Vert }_{B_1}}$$

$$c_0|v{|}_{H^1}^2\le d(v,v);|d(v,w)|\le c_1|v{|}_{H^1}|w{|}_{H^1},for\; allv,w\in H_0^1.$$

$$d(w,\theta )=⟨f,\theta ⟩,for\; all\theta \in H_0^1.$$

$$v=\sum _{j=1}^{\infty }(v,{\psi }_j){\psi }_j$$

$${\dot{H}}^s:=\left\{v\in L^2\text{s.t.}\sum _{j=1}^{\infty }{\lambda }_j^s|(v,{\psi }_j){|}^2<\infty \right\}.$$

$$\|v\|_{{\dot{H}}^{s}}:=\bigg{(}\sum_{j=1}^{\infty}\lambda_{j}^{s}|(v,\psi_{j})|^{2}\bigg{)}^{1/2}.$$

$${\dot{H}}^{-s}=\left\{F\in H^{-1}\ \text{s.t}\ \|F\|_{{\dot{H}}^{-s}}:=\bigg{(}\sum_{j=1}^{\infty}\lambda_{j}^{-s}|\langle F,\psi_{j}\rangle|^{2}\bigg{)}^{1/2}<\infty\right\}.$$

$$L^{\beta }v:=\sum _{j=1}^{\infty }{\lambda }_j^{\beta }(v,{\psi }_j){\psi }_j,\forall v\in D(L^{\beta }):={\dot{H}}^{2\beta }.$$

$$A(v,w):=(L^{\beta /2}v,L^{\beta /2}w)=\sum _{j=1}^{\infty }{\lambda }_j^{\beta }(v,{\psi }_j)\overline{(w,{\psi }_j)},$$

$${\mathbb{H}}^s:=\{\begin{array}{rlr}& [H_0^1,H_0^1\cap H^2{]}_{s-1,2},& 1\le s\le 2,\\ & [L^2,H_0^1{]}_{s,2},& 0\le s\le 1,\\ & [H^{-1},L^2{]}_{s+1,2},& -1\le s\le 0.\end{array}$$

$$d(u,v)={\int }_{\Omega }a(x)\nabla u\cdot \nabla vdx$$

$$0<c_0\le a(x)\le c_1\text{for a.e.}x\in \Omega .$$

$$e_{\gamma ,\mu }(z):=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (k\gamma +\mu )},z\in \mathbb{C}.$$

$${\partial }_t^{\gamma }u+\lambda u=0,$$

$${\partial }_t{e}_{\gamma ,1}(-t^{\gamma }{\lambda }^{\beta })={\lambda }^{\beta }{t}^{\gamma -1}{e}_{\gamma ,\gamma }(-t^{\gamma }{\lambda }^{\beta })$$

$${\partial }_t{e}_{\gamma ,\gamma }(-t^{\gamma }{\lambda }^{\beta })={\lambda }^{\beta }{t}^{\gamma -1}\left((\gamma -1)e_{\gamma ,2\gamma }(-t^{\gamma }{\lambda }^{\beta })-e_{\gamma ,2\gamma -1}(-t^{\gamma }{\lambda }^{\beta })\right).$$

$$|e_{\gamma ,\mu }(z)|\le \frac{C}{1+|z|},for\zeta \le |\arg (z)|\le \pi .$$

$$u(t)=E(t)v+{\int }_0^{t}W(s)f(t-s)ds,$$

Full text

We develop and analyze semi-discrete methods based on finite element approximation to the underlying (non-fractional) spatial operator in terms of analogous Dunford-Taylor integrals applied to the discrete operator. The space error is of optimal order up to a logarithm of $1/h$. The fully discrete method for the initial value problem is developed from the semi-discrete approximation by applying an exponentially convergent sinc quadrature technique to approximate the Dunford-Taylor integral of the discrete operator and is free of any time stepping.

To approximate the convolution appearing in the semi-discrete approximation to the non-homogeneous problem, we apply a pseudo midpoint quadrature. This involves the average of $W_{h}(s)$, (the semi-discrete approximation to $W(s)$) over the quadrature interval. This average can also be written as a Dunford-Taylor integral. We first analyze the error between this quadrature and the semi-discrete approximation. To develop a fully discrete method, we then introduce sinc quadrature approximations to the Dunford-Taylor integrals for computing the averages.

1. Introduction

In this paper, we investigate the numerical approximation to the following time dependent problem: given a bounded Lipschitz polygonal domain $\Omega$, a final time ${\mathsf{T}}>0$, an initial value $v\in L^{2}(\Omega)$ (a complex valued Sobolev space) and a forcing function $f\in L^{\infty}(0,{\mathsf{T}};L^{2}(\Omega))$, we seek $u:[0,{\mathsf{T}}]\times\Omega\rightarrow{\mathbb{R}}$ satisfying

[TABLE]

Here the fractional derivative in time $\partial^{\gamma}_{t}$ with ${\gamma}\in(0,1)$ is defined by the left-sided Caputo fractional derivative of order ${\gamma}$,

[TABLE]

Note that (2) holds for smooth $u$ and extends by continuity to a bounded operator on $H^{\gamma}(0,{\mathsf{T}})\cap C[0,{\mathsf{T}}]$ satisfying

[TABLE]

where ${}^{R}\partial^{\gamma}_{t}$ denotes the Riemann-Liouville fractional derivative. The differential operator $L$ appearing in (1) is an unbounded operator associated with a Hermitian, coercive and sesquilinear form $d(\cdot,\cdot)$ on $H^{1}_{0}(\Omega)\times H^{1}_{0}(\Omega)$. For ${\beta}\in(0,1)$, the fractional differential operator $L^{\beta}$ is defined by the following eigenfunction expansion

[TABLE]

where $(\cdot,\cdot)$ denotes the $L^{2}(\Omega)$ inner product and $\{\psi_{j}\}$ is an $L^{2}(\Omega)$-orthonormal basis of eigenfunctions of $L$ with eigenvalues $\{\lambda_{j}\}$. The above definition is valid for $v\in D(L^{\beta})$, where $D(L^{\beta})$ denotes the functions $v\in L^{2}(\Omega)$ such that $L^{\beta}v\in L^{2}(\Omega)$. A weak formulation of (1) reads: find $u\in L^{2}(0,{\mathsf{T}};D(L^{{\beta}/2}))\cap C([0,{\mathsf{T}}];L^{2}(\Omega))$ and $\partial_{t}^{\gamma}u\in L^{2}(0,{\mathsf{T}};D(L^{-{\beta}/2}))$ satisfying

[TABLE]

Here the bilinear form $A(u,\phi):=(L^{{\beta}/2}u,L^{{\beta}/2}\phi)$ and $\langle\cdot,\cdot\rangle$ denotes the duality pairing between $D(L^{-\beta/2})$ and $D(L^{\beta/2})$. As a consequence of [17, Theorem 6], the above problem has a unique solution, which can be explicitly written as

[TABLE]

Here, for $w\in L^{2}(\Omega)$,

[TABLE]

and

[TABLE]

with $e_{\gamma,\mu}(z)$ denoting the Mittag-Leffler function (see the defintion (13)). We also refer to Theorem 2.1 and 2.2 of [19] for a detailed proof of the above formula when ${\beta}=1$, noting that the argument is similar for any ${\beta}\in(0,1)$.

A major difficulty in approximation solutions of (4) involves time stepping in the presence of the fractional time derivative. The L1 time stepping method was developed in [14] and applied for the case ${\beta}=1$. Letting ${\tau}$ be the time step, it was shown in [14] that the L1 scheme gives the rate of convergence $O({\tau}^{2-\gamma})$ provided that the solution is twice continuously differentiable in time. For the homogeneous problem ($f=0$), the L1 scheme is guaranteed to yield first order convergence assuming the initial data $v$ is in $L^{2}(\Omega)$ (see [10]). See also [11] and the reference therein for other time discretization methods and error analyses. We also refer to [13] for the backward time stepping scheme for the case ${\gamma}=1$.

The numerical approximation to the solution (5) has been studied recently in [17]. The main difficulty is to discretize the fractional differential operators $\partial_{t}^{\gamma}$ and $L^{\beta}$ simultaneously. In [16], the factional-in-space operator $L^{\beta}$ was approximated as a Dirichlet-to-Neumann mapping via a Caffarelli-Silvestre extension problem [8] on $\Omega\times(0,\infty)$. In [17], Nochetto et. al. analyze an L1 time stepping scheme for (4) in the context of the Caffarelli-Silvestre extension problem and obtain a rate of convergence in time of $O({\tau}^{\theta})$ with $\theta\in(0,1/2)$ (see Theorem 3.11 in [17]).

The goal of the paper is to approximate the solution of (4) directly based on the solution formula (5). Our approximation technique and its numerical analysis relies on the Dunford-Taylor integral representation of the solution formula (5). Such a numerical method has been developed for the classical parabolic problem [3, 13] (i.e. the case ${\gamma}=1$) and the stationary problem [4]; see also [5] when the differential operator $L$ is regularly accretive [12].

The outline of the remainder of the paper is as follows. Section 2 provides some notation and preliminaries related to (1). In Section 3, we review some classical results from the finite element discretization and provide a key result (Theorem 3.3) instrumental to derive error estimates for semi-discrete schemes. In Section 4, we study the semi-discrete approximation $E_{h}(t)v:=e_{{\gamma},1}(-t^{\gamma}L_{h}^{\beta}){\pi_{h}}v$ to $E(t)v$. Here $L_{h}$ is the Galerkin finite element approximation of $L$ in the continuous piecewise linear finite element space $\mathbb{V}_{h}$ and ${\pi_{h}}$ denote the $L^{2}$ projection onto $\mathbb{V}_{h}$. We subsequently apply a sinc quadrature scheme to the Dunford-Taylor integral representation of the semi-discrete solution. For the sinc approximation, we choose the hyperbolic contour $z(y)=b(\cosh(y)+i\sinh(y))$ for $y\in{\mathbb{R}}$, with $b\in(0,\lambda_{1}/\sqrt{2})$. Here $\lambda_{1}$ denotes the smallest eigenvalue of $L$. Theorem 3.3 directly gives an error estimate for the semi-discrete approximation in fractional Sobolev spaces of order $s$, with $s\in[0,1]$. As expected, the rate of convergence depends on the smoothness of the solution which, in term, depends on the smoothness of the initial data and the regularity pickup associated with the spatial exponent $\beta$. Theorem 4.3 proves that for a quadrature of $2N+1$ points with quadrature spacing $k=cN^{-1/2}$ and $c$ depending on ${\beta}$, the sinc quadrature error is bounded by $Ct^{-{\gamma}}\exp(-c\sqrt{N})$, where the constant $C$ is independent of $t$ and $N$. In Section 5, we focus on the approximation scheme for the non-homogeneous forcing problem. The approximation in time is based on a pseudo-midpoint quadrature applied to the convolution in (5), i.e., given a partition $\{t_{j}\}$ on $[0,t]$,

[TABLE]

where $W_{h}(t)$ is the semi-discrete approximation to $W(t)$. Assuming that the forcing function $f$ is in $H^{2}(0,t;L^{2})$, We show in Theorem 5.3 that the error in the approximation (8) in time is $O({\mathcal{N}}^{-2})$ under a geometric partition refined towards $t=0$ (with $C({\gamma}){\mathcal{N}}\log_{2}{\mathcal{N}}$ subintervals). We then apply an exponentially convergent sinc quadrature scheme to approximate the Dunford-Taylor integral representation of the discrete operator $\int_{t_{j-1}}^{t_{j}}W_{h}({r})\,d{r}$. Theorem 5.5 shows that the sinc quadrature leads to an additional error which is $O(\log({\mathcal{N}})\exp(-\sqrt{cN}))$. Some technical proofs are given in Appendices A and B.

Throughout this paper, $c$ and $C$ denote generic constants. We shall sometimes explicity indicate their dependence when appropriate.

2. Notation and Preliminaries

2.1. Notation

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded polygonal domain with Lipschitz boundary. Denote by $L^{2}(\Omega)$ and $H^{1}(\Omega)$ (or in short $L^{2}$ and $H^{1}$) the standard Sobolev spaces of complex valued functions equipped with the norms

[TABLE]

The $L^{2}$ scalar product is denoted $(\cdot,\cdot)$:

[TABLE]

We also denote by $H^{1}_{0}:=H^{1}_{0}(\Omega)\subset H^{1}(\Omega)$ the closed subspace of $H^{1}$ consisting of functions with vanishing traces. Thanks to the Poincaré inequality, we will use the semi-norm $|\cdot|_{H^{1}}:=\||\nabla(\cdot)|\|$ as the norm on $H^{1}_{0}$. The dual space of $H^{1}_{0}$ is denoted $H^{-1}:=H^{-1}(\Omega)$ and is equipped with the dual norm:

[TABLE]

where $\langle\cdot,\cdot\rangle$ stands for the duality pairing between $H^{-1}$ and $H^{1}_{0}$.

The norm of an operator $A:B_{1}\rightarrow B_{2}$ between two Banach spaces $(B_{1},\|.\|_{B_{1}})$ and $(B_{2},\|.\|_{B_{2}})$ is given by

[TABLE]

and in short $\|A\|$ when $B_{1}=B_{2}=L_{2}$.

2.2. The Unbounded Operator $L$

Let us assume that $d(\cdot,\cdot)$ is a Hermitian, coercive and sesquilinear form on ${H^{1}_{0}}\times{H^{1}_{0}}$. We denote by $c_{0}$ and $c_{1}$ the two positive constants such that

[TABLE]

Furthermore, we let $T:H^{-1}\rightarrow H^{1}_{0}$ be the solution operator, i.e. for $f\in H^{-1}$, $Tf:=w\in H^{1}_{0}$, where $w$ is the unique solution (thanks to Lax-Milgram lemma) of

[TABLE]

Following Section 2 of [12], see also Section 2.3 in [5], we denote $L$ to be the inverse of $T|_{L^{2}}$ and define $D(L):=\text{Range}(T|_{L^{2}})$.

2.3. The Dotted Spaces

The operator $T$ is compact and symmetric on $L^{2}$. Fredholm theory guarantees the existence of an $L^{2}$-orthonormal basis of eigenfunctions $\{\psi_{j}\}_{j=1}^{\infty}$ with non-increasing real eigenvalues $\mu_{1}\geq\mu_{2}\geq\mu_{3}\geq...>0$. For every positive integer $j$, $\psi_{j}$ is also an eigenfunction of $L$ with corresponding eigenvalue $\lambda_{j}=1/\mu_{j}$. The decay of the coefficients $(v,\psi)$ in the representation

[TABLE]

characterizes the dotted spaces ${\dot{H}}^{s}$. Indeed, for $s\geq 0$, we set

[TABLE]

On ${\dot{H}}^{s}$, we consider the natural norm

[TABLE]

We also denote by ${\dot{H}}^{-s}$ the dual space of ${\dot{H}}^{s}$ for $s\in[0,1]$. It is known that (see for instance [5])

[TABLE]

Note that, we identify $L^{2}$ functions with $H^{-1}$ functionals by $\langle F,\cdot\rangle:=(f,\cdot)\in H^{-1}$ for $f\in L^{2}$.

2.4. Fractional Powers of Elliptic Operators

Let $L$ be defined from a Hermitian, coercive and sesquilinear form on ${H^{1}_{0}}\times{H^{1}_{0}}$ as described in Section 2.2. For ${\beta}\in(0,1)$, the fractional power of $L$ is given by

[TABLE]

In addition, we define the associated sesquilinear form $A:{\dot{H}}^{\beta}\times{\dot{H}}^{\beta}\rightarrow\mathbb{C}$ by

[TABLE]

which satisfies $A(v,v)=\|v\|_{{\dot{H}}^{{\beta}}}^{2}$.

2.5. Intermediate Spaces and the Regularity Assumption

As we saw above, the dotted spaces relies on the eigenfunction decomposition of a compact operator. These are natural spaces to consider fractional powers of operators but are less adequate to describe standard smoothness properties. The latter are better characterized by the intermediate spaces ${\mathbb{H}}^{s}$ defined for $s\in[-1,2]$ by real interpolation

[TABLE]

In order to link the two set of functional spaces introduced above, we assume the following elliptic regularity condition:

Assumption 2.1** (Elliptic Regularity).**

There exists $\alpha\in(0,1]$ so that

1. (a)

$T$* is a bounded map of ${\mathbb{H}}^{-1+\alpha}$ into ${\mathbb{H}}^{1+\alpha}$;* 2. (b)

$L$* is a bounded operator from ${\mathbb{H}}^{1+\alpha}$ to ${\mathbb{H}}^{-1+\alpha}$.*

Under the above assumption we have the following equivalence property:

Proposition 2.1** (Equivalence, Proposition 4.1 in [4]).**

Suppose that Assumption 2.1 holds for $\alpha\in(0,1]$. Then the spaces ${\mathbb{H}}^{s}$ and $\dot{H}^{s}$ coincide for $s\in[-1,1+\alpha]$ with equivalent norms.

Notice that Assumption 2.1 is quite standard and holds for a large class of sesquilinear forms $d(\cdot,\cdot)$. An important example is the diffusion process given by

[TABLE]

defined on $H^{1}_{0}\times H^{1}_{0}$, where $a\in L^{\infty}(\Omega)$ satisfies

[TABLE]

The $\alpha$ in Assumption  2.1 is related to the domain $\Omega$ and the smoothness of the coefficients. For example, if $\Omega$ is convex and $a$ is smooth, Assumption  2.1 holds for any $\alpha$ in $(0,1]$. In contrast, for the two dimensional L-shaped domain and smooth $a$, Assumption  2.1 only holds for $\alpha\in(0,2/3)$.

2.6. The Mittag-Leffler Function

The Mittag-Leffler functions are instrumental to represent the solution of fractional time evolution, see (6) and (7). We briefly introduce them together with their properties used in our argumentation. We refer to Section 1.8 in [20] for more details.

For ${\gamma}>0$ and $\mu\in{\mathbb{R}}$, the two-parameter Mittag-Leffler function ${e_{{\gamma},\mu}}(z)$ is defined by

[TABLE]

These functions are entire functions (analytic in ${\mathbb{C}}$). We note that [20, equation (3.1.42)] (see also [9]) implies that $u(t)={e_{{\gamma},1}}(-\lambda t^{\gamma})$ for $t,\lambda>0$ satisfies

[TABLE]

i.e., is a solution of the scalar homogeneous version of the first equation of (1). For this reason, the function ${e_{{\gamma},1}}(-\lambda t^{\gamma})$ will play a major role in our analysis. We also note that

[TABLE]

and

[TABLE]

Recall that $\partial^{\gamma}_{t}$ always denotes the left-sided Caputo fractional derivative (2).

Another critical property for our study is their decay when $|z|\to\infty$ in a positive sector: For $0<{\gamma}<1$, $\mu\in{\mathbb{R}}$ and $\frac{{\gamma}\pi}{2}<\zeta<{\gamma}\pi$, there is a constant $C$ only depending on ${\gamma},\mu,\zeta$ so that

[TABLE]

2.7. Solution via superposition

The solution $u$ of (4) is the superposition of two solutions: the homogeneous solution $f=0$ and the non-homogeneous solution $v=0$,

[TABLE]

where $E(t)$ is defined by (6) and $W(t)$ by (7). Following [19], we have that $u\in C^{0}([0,T];L^{2})$ and in particular $u(0)=v$.

We discuss the approximation of each term in the decomposition separately. For the homogeneous problem ($f=0$), we use the Dunford-Taylor integral representation of $u(t)=E(t)v$,

[TABLE]

Here $R_{z}(L):=(zI-L)^{-1}$ and $z^{{\beta}}:=e^{{{\beta}}\ln z}$ with the logarithm defined with branch cut along the negative real axis. Given $r_{0}\in(0,\lambda_{1})$, the contour ${\mathcal{C}}$ consists of three segments (see Figure 1):

[TABLE]

We use an analogous representation for $W(s)$, namely,

[TABLE]

The justification of (18) and (20) are a consequence of (16) and standard Dunford-Taylor integral techniques, see, [21, 2] for additional details.

3. Finite Element Approximations

3.1. Subdivisions and Finite Element Spaces

Let $\{\mathcal{T}_{h}\}_{h>0}$ be a sequence of globally shape regular and quasi-uniform conforming subdivisions of $\Omega$ made of simplexes, i.e. there are positive constants $\rho$ and $c$ independent of $h$ such that if for $\uptau\in\mathcal{T}_{h}$, $h_{\uptau}$ denotes the diameter of $\uptau$ and $r_{\uptau}$ denotes the radius of the largest ball which can be inscribed in $\uptau$, then

[TABLE]

Fix $h>0$ and denote by $\mathbb{V}_{h}\subset{H^{1}_{0}}$ the space of continuous piecewise linear finite element functions with respect to $\mathcal{T}_{h}$ and by ${M_{h}}$ the dimension of $\mathbb{V}_{h}$.

The $L^{2}$ projection onto $\mathbb{V}_{h}$ is denoted by ${\pi_{h}}:L^{2}\rightarrow\mathbb{V}_{h}$ and satisfies

[TABLE]

For $s\in[0,1]$ and $\sigma>0$ satisfying $s+\sigma\leq 2$, Lemma 5.1 in [5] guarantees the existence of a constant $c(s,\sigma)$ independent of $h$ such that

[TABLE]

In addition, for any $s\in[0,1]$, there exists a constant $c$ such that

[TABLE]

The case $s=0$ follows from the definition of the $L^{2}$-projection, the case $s=1$ is treated in [1, 6] and the general case follows by interpolation.

3.2. Discrete Operators

The finite element analogues of the operators $T$ and $L$ given in Section 2.2 are defined as follows: For $F\in H^{-1}$, $T_{h}:H^{-1}\rightarrow\mathbb{V}_{h}$ is defined by

[TABLE]

and $L_{h}:\mathbb{V}_{h}\rightarrow\mathbb{V}_{h}$ is given by

[TABLE]

so that $T_{h}|_{\mathbb{V}_{h}}=L_{h}^{-1}$.

We now recall the following finite element error estimates.

Proposition 3.1** (Lemma 6.1 in [5]).**

Let Assumption 2.1 (a) holds for some $\alpha\in(0,1]$. Let $s\in[0,\frac{1}{2}]$ and set $\alpha^{*}:=\frac{1}{2}(\alpha+\min(\alpha,1-2s))$. There is a constant $C$ independent of $h$ such that

[TABLE]

Similar to the operator $T$, $T_{h}|_{\mathbb{V}_{h}}$ has positive eigenvalues $\{\mu_{j,h}\}_{j=1}^{M_{h}}$ with corresponding $L^{2}$-orthonormal eigenfunctions $\{\psi_{j,h}\}_{j=1}^{M_{h}}$. The eigenvalues of $L_{h}$ are denoted as $\lambda_{j,h}:=\mu_{j,h}^{-1}$ for $j=1,2,\ldots,{M_{h}}$. Then the discrete fractional operator $L_{h}^{\beta}\ :\ \mathbb{V}_{h}\rightarrow\mathbb{V}_{h}$ is then given by

[TABLE]

Its associated sesquilinear form reads

[TABLE]

for all $v_{h},w_{h}\in\mathbb{V}_{h}$.

For any $s\in[0,1]$, the dotted spaces described in Section 2.3 also have discrete counterparts $\dot{H}_{h}^{s}$, which are characterized by their norms

[TABLE]

On $\mathbb{V}_{h}$, the two dotted norms are equivalent: For $s\in[0,1]$, there exists a constant $c$ independent of $h$ such that for all $v_{h}\in\mathbb{V}_{h}$,

[TABLE]

(see Appendix A.2 in [7]). From the property $\max_{j}{\lambda_{j,h}}\leq ch^{-2}$ (cf. [7, equation (2.8)]), we also deduce an inverse inequality in discrete dotted spaces: For $s,\sigma\geq 0$, we have

[TABLE]

3.3. The Semi-discrete Scheme in Space

We propose a Galerkin finite element method for the space discretization of (5). This is to find $u_{h}(t)\in\mathbb{V}_{h}$ satisfying

[TABLE]

where the bilinear form $A_{h}(\cdot,\cdot)$ is defined by (26) and $\pi_{h}$ is the $L^{2}$-projection onto $\mathbb{V}_{h}$. Similarly to the continuous case (see discussion in Section 2.7), the solution of the above discrete problem is given by

[TABLE]

where

[TABLE]

and ${\mathcal{C}}$ is as in (19).

3.4. A semi-discrete estimate

The purpose of this section (Theorem 3.3) is to obtain estimates for

[TABLE]

which, in view of representations (17) and (31), will be instrumental to derive error estimates for the space discretization.

The following lemma assesses the discrepancy between the resolvant $R_{z}(L)=(z-L)^{-1}$ and its finite element approximation. Its somewhat technical proof is postponed to Appendix A.

Lemma 3.2** (Space Discretization of Resolvant).**

Assume that Assumption 2.1 holds for some $\alpha\in(0,1]$. Let $s\in[0,\frac{1}{2}]$ and $\delta\in[0,(1+\alpha)/2]$. Then, there exists a positive constant $C$ independent of $h$ such that for all $\tilde{\alpha}$ with $2\tilde{\alpha}\in(0,\alpha+\min(\alpha,1-2s)]$, $z\in{\mathcal{C}}$ and $v\in{\dot{H}}^{2\delta}$

[TABLE]

We are now in a position to prove the error estimate for the semi-discrete approximation in space. Before doing so, for $s\in[0,1/2]$ and $0<\epsilon\ll 1$, we set

[TABLE]

We assume that

[TABLE]

The assumption (36) is sufficient to guarantee that the solution $e_{\gamma,\mu}(-t^{\gamma}L^{\beta})$ is in ${\dot{H}}^{2s+\epsilon}$ and we have the following theorem.

Theorem 3.3** (Space Discretization of $e_{\gamma,\mu}(-t^{\gamma}L^{\beta})$).**

Let $0<\gamma<1$, $s\in[0,1/2]$, $\mu\in\mathbb{R}$ and $\alpha^{*}$ be as in (35). Assume that Assumption 2.1 holds for $\alpha\in(0,1]$, and that $\delta$ satisfies (36). Then there exists a constant $C$ such that

[TABLE]

where

[TABLE]

Proof.

Without loss of generality we assume that $2\delta\leq 1+{\alpha^{*}}$ as the case $2\delta>1+{\alpha^{*}}$ follows from the continuous embedding

[TABLE]

Also, we use the notation $E^{\gamma,\mu}(t):=e_{\gamma,\mu}(-t^{\gamma}L^{\beta})$, $E^{\gamma,\mu}_{h}(t):=e_{\gamma,\mu}(-t^{\gamma}L_{h}^{\beta})$ and decompose the error in two terms:

[TABLE]

$\boxed{1}$ For the first term on the right hand side above, we note that the assumptions on the parameters imply that ${\alpha^{*}}+s\leq(\alpha+1)/2\leq 1$ and so the approximation property (23) of $\pi_{h}$ yields

[TABLE]

We estimate $\|E^{\gamma,\mu}(t)v\|_{\dot{H}^{2({\alpha^{*}}+s)}}$ by expanding $v$ in Fourier series with respect to the eigenfunctions of $L$ (see Section 2.3) and denote by $c_{j}:=(v,{\psi}_{j})$ the Fourier coefficient of $v$ so that

[TABLE]

Two cases need to be considered:

Case 1: $\delta\geq{\alpha^{*}}+s$. Here, the regularity of the initial condition is large enough to directly use the bound $|{e_{{\gamma},\mu}}(-t^{\gamma}\lambda_{j}^{\beta})|\leq C$ deduced from (16) to get

[TABLE]

Case 2: $\delta<{\alpha^{*}}+s$. In this case, we need to rely on the parabolic regularity for $t>0$. We apply (16) again and obtain

[TABLE]

Noting that $0<{\alpha^{*}}+s-\delta<{\beta}$, a Young’s inequality implies

[TABLE]

Whence,

[TABLE]

Returning to (39) after gathering the estimates obtained for the two different cases, we obtain

[TABLE]

$\boxed{2}$ We return to (38) and estimate now $\|\pi_{h}(E(t)-E_{h}(t)\pi_{h})v\|_{{\dot{H}}^{2s}}$. This time we use the integral representations and the resolvant approximation (Lemma 3.2) to get

[TABLE]

Furthermore, the decay estimate (16) of the Mittag-Leffler function evaluated at $-t^{\gamma}z^{\beta}$ for $z\in{\mathcal{C}}$ yields

[TABLE]

$\boxed{3}$ To prove

[TABLE]

it remains to show that

[TABLE]

This is done separately on each part of the contour ${\mathcal{C}}$, see (19). On ${\mathcal{C}}_{2}$, $|z|=r_{0}$ so that we directly have

[TABLE]

On ${\mathcal{C}}_{1}\cup{\mathcal{C}}_{3}$, we use the parametrization $z(r)=re^{{\pm}i\pi/4}$ to write

[TABLE]

When $\delta>{\alpha^{*}}+s$, we have enough decay to directly obtain

[TABLE]

When $\delta\leq{\alpha^{*}}+s$, we perform the change of variable $y:=t^{\gamma}|z|^{\beta}$ and obtain

[TABLE]

Thus,

[TABLE]

$\boxed{4}$ Gathering the estimates for each part of the contour yields (43) and thus (42), which, combined with (40), yields the desired result. ∎

4. Approximation of the Homogeneous Problem

This section presents and analyzes the proposed approximation algorithm in the case $f=0$. We note that the bound for the finite element approximation for the space discretization error is contained in Theorem 3.3. In this section, we define a sinc quadrature approximation to $E_{h}(t)$ and analyze the resulting quadrature error.

4.1. The Sinc Quadrature Approximation

We discuss the approximation of the contour integral in

[TABLE]

The first step involves replacing the contour ${\mathcal{C}}$ by one more suitable for application of the sinc quadrature technique. For $y\in{\mathbb{C}}$, we set

[TABLE]

and, for $0<b<\lambda_{1}/\sqrt{2}$, consider the hyperbolic contour ${{{\mathcal{C}}^{\prime}}}:=\{z(y)\ :\ y\in{\mathbb{R}}\}$. Using this contour, we have

[TABLE]

Given a positive integer $N$ and a quadrature spacing $k>0$, we set $y_{j}:=jk$ for $j=-N,...,N$ and define the sinc quadrature approximation of $e_{\gamma,1}(-t^{\gamma}L_{h}^{\beta})g_{h}$ by

[TABLE]

4.2. Quadrature Error

We now discuss the quadrature error. Expanding $(E_{h}(t)-Q_{h,k}^{N}(t))g_{h}$ in term of the discrete eigenfunction $\{\psi_{j,h}\}_{j=1}^{M_{h}}$ (see Section 3.2), for $s>0$ we have

[TABLE]

where

[TABLE]

and

[TABLE]

The function $g_{\lambda}(y,t)$ is well defined for $t>0$, $\lambda\geq\lambda_{1}$, $y\in{\mathbb{C}}$ with $z(y)\neq\lambda$ and $z(y)$ not on the branch cut for the logarithm.

Following [15], we show that when $k=c/\sqrt{N}$ for some constant $c$, the quantity ${\mathcal{E}}(\lambda,t)\to 0$ when $k\to 0$ uniformly with respect to $\lambda\geq\lambda_{1}$. Moreover, the convergence rate is $O(\exp{(-c\sqrt{N})})$. We then use this estimate in (47) to deduce exponential rate of convergence for the sinc quadrature scheme (46).

This program requires additional notations and we start with the class of functions $S(B_{d})$.

Definition 4.1**.**

Given $d>0$, we define the space $S(B_{d})$ to be the set of functions $f$ defined on ${\mathbb{R}}$ having the following properties:

1. (i)

$f$* extends to an analytic function in the infinite strip*

[TABLE]

and is continuous on $\overline{B_{d}}$. 2. (ii)

There exists a constant $C$ independent of $y\in{\mathbb{R}}$ such that

[TABLE] 3. (iii)

We have

[TABLE]

Note that condition $(ii)$ is more restrictive than actually needed (see Definition 2.12 in [15]) but sufficient for our considerations. In addition, For $f\in S(B_{d})$, Theorem 2.20 in [15] provides the error estimate for the quadrature approximation to $\int_{\mathbb{R}}f(x)\,dx$ using an infinite number of equally spaced quadrature points with spacing $k>0$:

[TABLE]

The lemma below is proved in Appendix B and is the first step in estimating the sinc quadrature error.

Lemma 4.1**.**

Let $\lambda\geq\lambda_{1}$ and $t>0$. The function $w\mapsto g_{\lambda}(w,t)$ belongs to $S(B_{d})$ for $0<d<\pi/4$. Moreover, there exists a constant $C$ only depending on ${\beta}$, $d$ and $b$ such that

[TABLE]

The above lemma together with the quadrature estimate (50) leads to exponential decay for ${\mathcal{E}}(\lambda,t)$ as provided in the following lemma.

Lemma 4.2**.**

Let $0<d<\pi/4$. There exists a constant $C$ only depending on $d$, $b$, $\beta$ and $\lambda_{1}$ such that for $k<1$, $N>0$, $t>0$ and $\lambda\geq\lambda_{1}$,

[TABLE]

Proof.

In order to derived the desired estimate, we write

[TABLE]

Lemma 4.1 guarantees that $g_{\lambda}(.,t)\in S(B_{d})$ and so in view of (50), we obtain

[TABLE]

where $C$ is the constant in (51). For the truncation term, we use (83) (in the appendix) to write

[TABLE]

where $C$ is a constant only depending on $d$, $b$ and $\lambda_{1}$. Next we bound the infinite sum by the integral and arrive at

[TABLE]

where now the constant depends on $\beta$ as well. Gathering the above estimates completes the proof. ∎

Remark 4.1** (Choice of $k$ and $N$).**

The optimal combination of $k$ and $N$ is obtained by balancing the two exponentials on the right hand side of (52). Hence, we select $k$ and $N$ such that $\pi d/k={\beta}Nk$, i.e. $k=\sqrt{\frac{\pi d}{{\beta}N}}$, and the estimate on ${\mathcal{E}}(\lambda,t)$ becomes

[TABLE]

Estimates on the difference between $E_{h}(t)$ defined by (31) and $Q_{h,k}^{N}$ defined by (46) follow from (53) and (47) as stated in the following theorem.

Theorem 4.3**.**

Let $s\in[0,1/2]$, $d\in(0,\pi/4)$, and let $N$ be a positive integer. Set $k=\sqrt{\frac{\pi d}{{\beta}N}}$. Then there exists a constant $C$ independent of $k$, $N$, $t$ and $h$ such that for every $g_{h}\in{\dot{H}}_{h}^{2s}$

[TABLE]

4.3. The Total Error

The discrete approximation after space and quadrature discretization is

[TABLE]

Gathering the space and quadrature error estimates, we obtain the final estimate for the approximation of the homogeneous problem.

Theorem 4.4** (Total error).**

Assume that the conditions of Theorem 3.3 and Theorem 4.3 hold. Then there exists a constant $C$ independent of $h$, $t$ and $N$ such that

[TABLE]

provided the initial condition $v$ is in ${\mathbb{H}}^{2s}\cap{\mathbb{H}}^{2\delta}$. Here $D(t)$ is the constant given by (37).

Proof.

We use the decomposition

[TABLE]

and invoke Theorem 3.3 with $\mu=1$ and Lemma 4.3 with $g_{h}=\pi_{h}v$ to arrive at

[TABLE]

The equivalence of norms (28) together with stability of the $L^{2}$ projection (24) and the equivalence property between the dotted spaces and interpolation spaces (12) (see Proposition 2.1) yield the desired result. ∎

Remark 4.2** (Implementation).**

Denote $\widetilde{U}(t)$ the vector of coefficients of $u_{h}^{N}(t)$ with respect to the finite element local basis functions and $\widetilde{V}$ the vector of inner product between $v$ and local basis functions. Let $\widetilde{A}$ and $\widetilde{M}$ be the stiffness and mass matrices. Then

[TABLE]

Remark 4.3** (Complexity of the Implementation).**

We take advantage of the exponential decay of the sinc quadrature by setting $N=c({\alpha^{*}}\ln(1/h))^{2}$ so that

[TABLE]

Hence, computing $u_{h}^{N}(t)$ for a fixed $t$ requires $O(\log(1/h)^{2})$ complex finite element system solves.

4.4. Numerical Illustration

In this section, we provide numerical illustrations of the rate of convergence predicted by Theorem 3.3 and Lemma 4.3.

Space Discretization Error

In order to illustrate the space discretization error, we start with a one dimensional problem and use a spectral decomposition to compute the exact solution without resorting to quadrature. Set $\Omega=(0,1)$ , $Lu:=-u^{\prime\prime}$. We chose the initial condition to be $v\equiv 1$ or, using the eigenvalues $\lambda_{\ell}=\pi^{2}\ell^{2}$ and associated eigenfunctions $\psi_{\ell}(x)=\sqrt{2}\sin(\pi\ell x)$,

[TABLE]

The number of term used before the truncation is chosen large enough not to influence the space discretization ($50000$). With these notations, the exact solution for $\gamma=1/2$ and $0<\beta<1$ is approximated by

[TABLE]

For the space discretization, we consider a sequence of uniform meshes with mesh sizes $h_{j}=2^{-j}$, where $j=1,2,\dots$ and denote by $\{\varphi_{k,h}\}_{k=1,\ldots,M_{h_{j}}}$ the continuous piecewise linear finite element basis of $\mathbb{V}_{h}$. The eigenvalues of $L_{h_{j}}$ corresponds to the eigenvalues of $M_{h_{j}}^{-1}S_{h_{j}}$, where $M_{h_{j}}$ and $S_{h_{j}}$ are the mass and stiffness matrices and are given by

[TABLE]

The associated eigenfunctions to $L_{h}$ are

[TABLE]

Similar to (56), we use the discrete spectral representation below of $u_{h_{j}}(t)$ for our computation

[TABLE]

with

[TABLE]

Note that $\alpha$ in Assumption 2.1 is 1, $v\in{\dot{H}}^{1/2-\epsilon}$ for any $\epsilon>0$ so that $\delta=1/4-\epsilon$. The error will be computed in $L^{2}$ and $H^{1}$, i.e. $s=0$ and $s=1/2$. For the latter we need $\beta>1/4$. The predicted convergence rates (Theorem 3.3) are

[TABLE]

for every $\epsilon>0$, i.e.

[TABLE]

We use the MATLAB code [18] to evaluate $e_{\gamma,1}(z)$ for any $z\in\mathbb{C}$ and fix $t=0.5$. In Figure 2, we report the errors

[TABLE]

for $j=3,4,5,6,7$ and different values of ${\beta}$. The observed rate of convergence

[TABLE]

are also reported in this figure and match the rates predicted by Theorem 3.3.

Effect of the Sinc Quadrature

We examine the error between the semi-discrete approximation and its sinc quadrature approximation. To this end and in order to factor out the space discretization, it suffices to observe ${\mathcal{E}}({\lambda},t)$ defined by (48) for all $\lambda\geq\lambda_{1}$. Here we fix $t=0.5$ and approximate $\|{\mathcal{E}}(.,t)\|_{L^{\infty}(\lambda_{1},\infty)}$ with $\lambda_{1}=10$ using the method discussed in Section 5.2 in [3]. For the hyperbolic contour $z(y)$ in (45), we choose $b=1$ so that $b\in(0,\lambda_{1}/\sqrt{2})$. Following Remark 4.1, we fix the number of quadrature points to $2N+1$ and balance the two source of errors by setting $k=\sqrt{{\pi d}/{({\beta}N)}}$ with $d={\pi}/{8}$. According to (53), we have

[TABLE]

The left graph of Figure 3 illustrates the exponential decay of $\|{\mathcal{E}}(\lambda,t)\|_{L^{\infty}(10,\infty)}$ as $N$ increases for ${\gamma}=0.5$ and ${\beta}=0.3,0.5,0.7$. We also report (right) the singular behavior of $\|{\mathcal{E}}(.,t)\|_{L^{\infty}(10,\infty)}$ in time for $N=100$, ${\beta}=0.5$ and ${\gamma}=0.3,0.5,0.7$.

A Two Dimensional Problem

We now focus our attention to the total error in a two dimensional problem. Let $\Omega=(0,1)^{2}$, $L=-\Delta$ and the initial condition be the eigenfunction of $L$ given by

[TABLE]

The exact solution is then given by

[TABLE]

The space discretizations are subordinate to a sequence of uniform subdivisions made of triangles with the mesh size $h_{j}=2^{-j}\sqrt{2}$. For the quadrature, we chose $N=400$ and set $k=\sqrt{{\pi^{2}}/{(8{\beta}N)}}$ for the quadrature error not to affect the space discretization error. Since $\lambda_{1}=\pi^{2}$, we again set $b=1$ in (45). We fix $t=0.5$, $\gamma=0.5$ and report in Figure 4, the quantities $\|u(t)-u_{h_{j}}^{N}(t)\|$ for $j=3,4,5,6,7,8$ and different ${\beta}$. As announced in Theorem 4.4, a second order rate of convergence is observed.

5. Approximation of the Non-homogeneous Problem

We now turn our attention to the non-homogeneous problem, i.e. $f\neq 0$ and $v=0$ in (1), for which the solution reads

[TABLE]

5.1. The Semi-discrete Scheme

According to (31), the finite element approximation of (57) is given by

[TABLE]

As in the homogeneous case, the finite element approximation error is derived from Lemma 3.3 and we have the following lemma.

Lemma 5.1** (Space Discretization for the non-homogeneous problem).**

Assume that Assumption 2.1 holds for $\alpha\in(0,1]$. Let $\gamma\in(0,1)$, $s\in[0,\frac{1}{2}]$ and let $\alpha^{*}$ and $\delta$ be as in (35) and (36), respectively. There exists a constant $C$ such that

[TABLE]

where

[TABLE]

Proof.

Applying Theorem 3.3 gives

[TABLE]

where $D(t)$ is given by (37). The conclusion follow from $\int_{0}^{t}r^{\gamma-1}D(r)\,dr=\widetilde{D}(t)$. ∎

5.2. Time Discretization via Numerical Integration

Given a final time ${\mathsf{T}}$, we discuss first a numerical approximation of the integral

[TABLE]

For simplicity, we set

[TABLE]

so that the above integral becomes

[TABLE]

For a positive integer ${\mathcal{M}}$, let $0=t_{0}<t_{1}<...<t_{{\mathcal{M}}}={\mathsf{T}}$ be a partition of the time interval $[0,{\mathsf{T}}]$. On each subinterval we set $t_{j-\frac{1}{2}}=\frac{1}{2}(t_{j}+t_{j-1})$ and propose the pseudo-midpoint approximation

[TABLE]

where to achieve the last step, we used relation (14).

Before going further, we note that numerical methods based on (60) cannot perform optimally when using a uniform decomposition of the time interval because $W_{h}(t)$ is singular at $t=0$. Hence, the performance of algorithms based on uniform partitions are bound to the error on the first interval $(0,t_{1})$. Measuring in the ${\dot{H}}^{2s}$-norm for $s\in[0,1/2]$, we have

[TABLE]

To overcome this deterioration, we propose a geometric refinement of the partition near $t_{0}=0$ which depends on two positive integers ${\mathcal{M}}$ and ${\mathcal{N}}$ (see also Section 3.1 of [4]). We first set

[TABLE]

We decompose further all but the first interval

[TABLE]

onto ${\mathcal{N}}$ subintervals

[TABLE]

where, for $l=0,...,{\mathcal{N}}$,

[TABLE]

As in (60), we approximate

[TABLE]

on each subinterval $I_{j,l}:=[t_{j,l-1},t_{j,l}]$ by

[TABLE]

Here $t_{j,l-1/2}:=\frac{1}{2}(t_{j,l-1}+t_{j,l})$. We use the bar symbol to denote average quantities over the interval $[t_{j,l-1},t_{j,l}]$, e.g.,

[TABLE]

and

[TABLE]

The approximate solution after time integration is thus given by

[TABLE]

We start by assessing the local integration error

[TABLE]

Lemma 5.2** (Local Approximation).**

Let $\gamma\in(0,1)$ and $s\in[0,1/2]$. Let $j\geq 2$ and assume that $g(t)=f({\mathsf{T}}-t)$ belongs to $H^{2}(t_{j-1},t_{j};{\dot{H}}^{2s})$. There exists a constant $C$ independent of $h$, and $\tau_{j}$ such that on every interval $I_{j}=[t_{j-1},t_{j}]$, we have

[TABLE]

Proof.

We use the following decomposition on each sub-interval:

[TABLE]

$\boxed{1}$ We estimate $E_{1}$

[TABLE]

We now bound $\|{\overline{W}}_{j,l}\pi_{h}\|_{{\dot{H}}^{2s}\to{\dot{H}}^{2s}}$ and $\|{\overline{g}}_{j,l}-g({t}_{j,l-\frac{1}{2}})\|_{{\dot{H}}^{2s}}$ separately. For the latter, we expand $g(\eta)$ at $\eta=t_{j,l-\frac{1}{2}}$ to get

[TABLE]

where $g_{t}$ and $g_{tt}$ denote the first and second partial derivative in time of $g$. As a consequence, taking advantage of $t_{j,l-\frac{1}{2}}$ being the midpoint of the interval $I_{j,l}$, we obtain

[TABLE]

and so using a Cauchy-Schwarz inequality

[TABLE]

In order to bound $\|{\overline{W}}_{j,l}\pi_{h}\|_{{\dot{H}}^{2s}\to{\dot{H}}^{2s}}$, we note that from the definition of the discrete dotted spaces ${\dot{H}}^{2s}_{h}$ (see (27)), we have

[TABLE]

Therefore, from the expression of $W_{h}(t)$ in (58), the equivalence of norms (28) and the stability estimate (24) for $\pi_{h}$, we derive that

[TABLE]

Estimates (66) and (67) into (65) give the final bound for $E_{1}$

[TABLE]

$\boxed{2}$ We estimate $E_{2}$

[TABLE]

In this case as well, we need to estimate two terms separately, namely $\|(W_{h}({r})-{\overline{W}}_{j,l})\pi_{h}\|_{{\dot{H}}^{2s}\to{\dot{H}}^{2s}}$ and $\|g({r})-g({t}_{j,l-\frac{1}{2}})\|_{{\dot{H}}^{2s}}$. For the latter, we write

[TABLE]

Next, we bound $\|(W_{h}({r})-{\overline{W}}_{j,l})\pi_{h}\|_{{\dot{H}}^{2s}\to{\dot{H}}^{2s}}$. As before, it suffices to estimate $\|W_{h}({r})-{\overline{W}}_{j,l}\|_{{\dot{H}}^{2s}_{h}\to{\dot{H}}^{2s}_{h}}$. To achieve this, we use the eigenfunctions $\{\psi_{i,h}\}_{i=1}^{M_{h}}$ of $L_{h}$. By (15),

[TABLE]

This and (16) with $z=-r^{\gamma}\lambda_{i,h}^{\beta}$ imply that for $r\in I_{j,l}$,

[TABLE]

where the constant in the above inequality is independent of $j$, $l$ and $h$. Whence,

[TABLE]

and

[TABLE]

The above estimate and (70) in (69) yield the final bound on $E_{2}$

[TABLE]

[math] Summing up the contribution from each subinterval and using a Cauchy-Schwarz inequality, yields the desired result. ∎

Remark 5.1** (Uniform time-stepping).**

In the case of uniform time-stepping, i.e. ${\mathcal{N}}=0$ and $t_{j}=j\tau$, $\tau={\mathsf{T}}/{\mathcal{M}}$, we derive from the estimate provided in Lemma 5.2 and the first interval estimate (61) that the quadrature error behaves asymptotically like $\tau^{1+\gamma}$. We do not pursue this further but rather investigate errors coming from the geometric partition.

Theorem 5.3** (Time Discretization of Non-Homogeneous Problem).**

Let $\gamma\in(0,1)$, $s\in[0,1/2]$, ${\mathsf{T}}\geq{\mathsf{T}}_{0}>0$. Let ${\mathcal{N}}$ be a positive integer and

[TABLE]

Assume that $f$ is in $H^{2}(0,{\mathsf{T}};{\dot{H}}^{2s})$ and let $u_{h}^{{\mathcal{N}}}({\mathsf{T}}):=u_{h}^{{\mathcal{N}},{\mathcal{M}}}$ be defined by (64) and let $u_{h}({\mathsf{T}})$ be the semi-discrete in space solution (58). Then there exists a constant $C$ independent of ${\mathcal{N}}$, $h$ and ${\mathsf{T}}$ satisfying

[TABLE]

Proof.

Using the definitions of $u_{h}({\mathsf{T}})$ and $u_{h}^{{\mathcal{N}}}({\mathsf{T}})$, we write

[TABLE]

For the first term, we note that (16) immediately implies that $\|e_{\gamma,\gamma}(-r^{\gamma}L_{h}^{\beta})\|_{{\dot{H}}^{2s}\to{\dot{H}}^{2s}}\leq C$. The stability of the $L^{2}$ projection (24) and (72) give

[TABLE]

For the second term, we apply Lemma 5.2 on each interval $I_{j}$, $j=1,\ldots,{\mathcal{M}}-1$ to get

[TABLE]

where we use the fact that $C^{-1}t_{j}\leq t_{j,l}\leq Ct_{j}$ for some constant $C$ independent of ${\mathcal{N}}$ and ${\mathcal{M}}$. Hence, a Cauchy-Schwarz inequality and the definitions of $t_{j}$ and $\tau_{j}$ yield

[TABLE]

This, together with the estimate for the first interval, implies

[TABLE]

To conclude, we observe that

[TABLE]

and that the embedding $H^{1}(0,{\mathsf{T}})\subset L^{\infty}(0,{\mathsf{T}})$ is continuous with norm independent of ${\mathsf{T}}\geq{\mathsf{T}}_{0}$. ∎

5.3. A Sinc Approximation of the Contour Integral

In view of (63), one remaining problem is to compute

[TABLE]

for $t>0$, $\tau>0$ and $g_{h}\in\mathbb{V}_{h}$. We proceed as in the homogeneous case discussed in Section 4.1.

Let $N$ be a positive integer and let $k>0$ be a quadrature spacing. For $t,{\tau}>0$ and $g_{h}\in\mathbb{V}_{h}$, we propose the following sinc approximation of $H_{h}(t,\tau)$:

[TABLE]

where $z(y)$ for $y\in{\mathbb{R}}$ is the hyperbolic contour (45). With this, the computable approximation of the solution to the non-homogeneous problem becomes

[TABLE]

We start with the approximation of $H_{h}(t,\tau)$ by $Q_{h,k}^{N}(t,\tau)$.

Lemma 5.4**.**

Let $t,{\tau}>0$, $s\in[0,1/2]$ and $d\in(0,\pi/4)$. There exists a constant $C$ only depending on $d$, $b$, $\beta$, $\lambda_{1}$ such that for any $g_{h}\in\mathbb{V}_{h}$,

[TABLE]

Proof.

For $y\in B_{d}$, define

[TABLE]

and note that

[TABLE]

Here we applied (16) replacing $z$ with $-z(y)^{\beta}s^{\gamma}$ so that

[TABLE]

Hence, the desired estimate follows upon proceeding as in the proofs of Lemmas 4.1 and 4.2. ∎

We are now in a position to prove the error estimate for the sinc quadrature on the non-homogeneous problem.

Lemma 5.5**.**

Let ${\mathsf{T}}>0$, $s\in[0,1/2]$ and assume that $f\in L^{\infty}(0,{\mathsf{T}};{\dot{H}}^{2s})$. Let $N$ be a positive integer, $d\in(0,\pi/4)$ and set $k=\sqrt{\frac{\pi d}{{\beta}N}}$. Let $u_{h}^{{\mathcal{N}},{\mathcal{M}}}$ be as in (64) and let $u_{h,k}^{{\mathcal{N}},N}$ be as in (74). There exists a constant $C$ independent of $h$, ${\mathsf{T}}$, $k$, $N$, ${\mathcal{N}}$, ${\mathcal{M}}$ satisfying

[TABLE]

Proof.

Note that both $u^{{\mathcal{N}},{\mathcal{M}}}_{h}$ and $u^{{\mathcal{N}},{\mathcal{M}}}_{h,k}$ are approximations starting at $t_{1}$ (the first interval $I_{0}=[0,t_{1}]$ is skipped). Hence, applying Lemma 5.4 on each interval $I_{j,l}$ (i.e. with $\tau=\tau_{j}$, $t=t_{j,l}$ and $g_{h}=\pi_{h}f({\mathsf{T}}-t_{j,l-\frac{1}{2}})$) for $j=1,...,{\mathcal{M}}-1$ and $l=0,...,{\mathcal{N}}$, yields

[TABLE]

where we have used the definition (62) of $t_{j,l}$ to guarantee that $C^{-1}t_{j}\leq t_{j,l}\leq Ct_{j}$ as well as the definition of $\tau_{j}=2^{-({\mathcal{M}}-j)}{\mathsf{T}}/{\mathcal{N}}$. This is the desired result. ∎

5.4. Total Error

We summarize this section by the following total error estimate for the fully discrete approximation (74) to the solution of the non-homogeneous problem. Since $k=k(N)$ and ${\mathcal{M}}={\mathcal{M}}({\mathcal{N}})$, we denote by $u_{h}^{{\mathcal{N}},N}({\mathsf{T}})$ the fully discrete solution (74).

Theorem 5.6** (Total Error).**

Assume that Assumption 2.1 holds for $\alpha\in(0,1]$. Furthermore, let $\gamma\in(0,1)$, $\delta\geq 0$, $s\in[0,\min(1/2,\delta)]$ and let $\alpha^{*}$ be as in (35). Let ${\mathsf{T}}>0$, ${\mathcal{N}}$ a positive integer and ${\mathcal{M}}=\left\lceil\frac{2\log_{2}{{\mathcal{N}}}}{{\gamma}}\right\rceil$. Let $N$ be a positive integer, $d\in(0,\pi/4)$ and set $k=\sqrt{\frac{\pi d}{{\beta}N}}$. There exists a constant C independent of $h$, ${\mathcal{N}}$, ${\mathsf{T}}$ and $N$ such that for every $f\in H^{2}(0,T;{\mathbb{H}}^{2\delta})$ we have

[TABLE]

where $\tilde{D}(T)$ is given by (59).

Proof.

This is in essence Lemmas 5.1, 5.3 and 5.5 together with the equivalence property between the dotted spaces and interpolation spaces (12) (see Proposition 2.1).

∎

Remark 5.2** (Choice of ${\mathcal{N}}$ and $N$).**

In practice, we balance the three error terms in Theorem 5.6 by setting

[TABLE]

for some positive constants $c_{1}$ and $c_{2}$ so that the total error behaves like $h^{2{\alpha^{*}}}$. We note that the number of the finite element systems that need to be solved for the non-homogeneous problem is the same as for the homogeneous problem, i.e. $O(\ln(1/h)^{2})$ complex systems (see the numerical illustration below).

5.5. Numerical illustration

To minimize the number of system solves in the computation of (74), we rewrite

[TABLE]

where

[TABLE]

To implement the above we proceed as follows:

Compute the inner product vectors, i.e., the integral of $f(t-{t}_{j,l-1/2})$ against the finite element basis vectors, for all $(j,l)$. 2. 2)

For each, $n$:

1. a)

compute the sums in (75) but replacing $\pi_{h}f(T-t_{j,l-1/2})$ by the corresponding inner product vector, and 2. b)

compute $z(y_{n})^{-\beta}z^{\prime}(y_{n})(z(y_{n})I-L_{h})^{-1}{\mathcal{H}}_{n}$ by inversion of the corresponding stiffness matrix applied to the vector of Part a). 3. 3)

Sum up all contribution and multiply the result by $\frac{k}{2\pi i}$.

We illustrate the error behavior in time on a two dimensional problem with domain $\Omega=(0,1)^{2}$ and $L=-\Delta$ with homogeneous Dirichlet boundary conditions. We set ${\beta}=0.5$ and consider the exact soltuion $u(t,x_{1},x_{2})=t^{3}\sin{(\pi x_{1})}\sin{(\pi x_{2})}$ which vanishes at $t=0$. This corresponds to

[TABLE]

We partition $\Omega$ using uniform triangles with the mesh size $h=2^{-5}\sqrt{2}$ and use $N=400$ for the sinc quadrature parameter. We also set $b=1$ in the hyperbolic contour (45). In Figure 5 (left), we report $\|u(0.5)-Q_{h}^{{\mathcal{N}},N}(0.5)\|$ for ${\mathcal{N}}=2,4,8,16,32$ and different values of ${\gamma}$. In each cases, as predicted by Theorem 5.3, the rate of convergence ${\mathcal{N}}^{-2}$ is observed. For comparison, the approximation based on a uniform partition is also provided. In this case, the error decay behaves like $\tau^{1+\gamma}$ (see Remark 5.1).

Appendix A Proof of Lemma 3.2

The following lemma proved in [3] (see, Lemma 3.1 of [3]) and is instrumental in the proof of Lemma 3.2.

Lemma A.1**.**

There is a positive constant $C$ only depending on $s\in[0,1]$ such that

[TABLE]

The same inequality holds $\mathbb{V}_{h}$, i.e. with $T$ replaced by $T_{h}$ and $f\in\mathbb{V}_{h}$.

Proof of Lemma 3.2.

Noting that

[TABLE]

and

[TABLE]

we obtain

[TABLE]

where for the last step we used the definition of $T_{h}$ to deduce that $\pi_{h}(zT_{h}-I)^{-1}=(zT_{h}-I)^{-1}\pi_{h}$. We have left to prove:

[TABLE]

for a constant $C$ is independent of $h$ and $z$ and where

[TABLE]

To show this, we write

[TABLE]

where $\gamma:=2\alpha^{*}-\alpha$. We estimate the three terms on the right hand side above separately.

We start with $\mathrm{III}$ and use the definition of the dotted spaces (see Section 2.3) to write

[TABLE]

Applying Lemma A.1 (recall that $\delta\in[0,(1+\alpha)/2]$ and $\alpha\in[0,1]$ so that $(1+\alpha)/2-\delta\in[0,1]$), we obtain

[TABLE]

where $C$ is the constant in (76).

To estimate $\mathrm{I}$, we start with the equivalence of norms (28) so that

[TABLE]

Whence, the stability of the $L^{2}$ projection (24) together with the equivalence property between dotted spaces and interpolation spaces (Proposition 2.1) as well as the definition of the discrete dotted space norm (27) lead to

[TABLE]

We recall that $\alpha\in(0,1]$ and ${\gamma}=2\tilde{\alpha}-\alpha$ so that $(1+{\gamma})/2+s\in(0,1]$. Hence, Lemma A.1 ensures the following estimate:

[TABLE]

For the remaining term, Proposition 3.1 with $2s=1-{\gamma}$ gives

[TABLE]

Combining the above estimate with (79) and (81) yields (77) and completes the proof. ∎

Appendix B Sinc quadrature Lemma.

The results of the next lemma are contained in the proof of Theorem 4.1 of [3].

Lemma B.1**.**

Let $0<d<\pi/4$ and $\lambda>\lambda_{1}$. let $z(y)$ be defined by (45) and $B_{d}=\left\{z\in\mathbb{C}:\ \mathfrak{Im}(z)<d\right\}.$ The following assertions hold.

1. (a)

There exists a constant $C>0$ only depending on $\lambda_{1}$, $b$ and $d$ such that

[TABLE] 2. (b)

There exists a constant $C>0$ only depending on $\lambda_{1}$, $b$ and $d$ such that

[TABLE] 3. (c)

There is a constant $C>0$ only depending on $b$, $d$ and $\beta$ such that

[TABLE]

Proof of the Lemma 4.1.

From the expression (4.13) of $\mathfrak{Re}(z(y))$ in [3], we deduce that $\mathfrak{Re}(z(y))$ is strictly positive for $y\in\bar{B}_{d}=\{w\in\mathbb{C}:\ \mathfrak{Im}(w)\leq d\}$. It follows from this and Part (a) of Lemma B.1 that Condition (i) of Definition 4.1 holds for $g_{\lambda}(\cdot,t)$ for $\lambda\geq\lambda_{1}$ and $t>0$.

We now give a proof of (ii) and (iii) of Definition 4.1 simultaneously. Note that Part (b) in Lemma B.1 together with (16) imply that for $y\in\overline{B}_{d}$,

[TABLE]

Furthermore, the estimate on $\mathfrak{Re}(z(y)^{\beta})$ in Part (c) of Lemma B.1 yields

[TABLE]

This guarantees that

[TABLE]

and

[TABLE]

which yield (ii), (iii) and the bound on $N(B_{d})$. ∎