---
abstract: 'An initial-boundary value problem for the time-fractional diffusion equation is discretized in space using continuous piecewise-linear finite elements on a polygonal domain with a re-entrant corner. Known error bounds for the case of a convex polygon break down because the associated Poisson equation is no longer $H^2$-regular. In particular, the method is no longer second-order accurate if quasi-uniform triangulations are used. We prove that a suitable local mesh refinement about the re-entrant corner restores second-order convergence. In this way, we generalize known results for the classical heat equation due to Chatzipantelidis, Lazarov, Thomée and Wahlbin.'
author:
- Kim Ngan Le
- William McLean
- Bishnu Lamichhane
bibliography:
- 'nonconvexrefs.bib'
title: 'Finite element approximation of a time-fractional diffusion problem in a non-convex polygonal domain[^1]'
---

Introduction
============

In a standard model of subdiffusion [@KlafterSokolov2011], each particle undergoes a continuous-time random walk with a common waiting-time distribution that obeys a power law. Consequently, the mean-square displacement of a particle is proportional to $t^\alpha$ with $0<\alpha<1$, and the macroscopic concentration $u(x,t)$ of the particles satisfies the time-fractional diffusion equation $$\label{eq: fpde}
\partial_t u-\partial_t^{1-\alpha}K\nabla^2u=f(x,t).$$ Here, $\partial_t=\partial/\partial t$ and $\nabla^2$ denotes the spatial Laplacian. The fractional time derivative is of Riemann–Liouville type: $$\partial_t^{1-\alpha}v(x,t)=\frac{\partial}{\partial t}\int_0^t
    \omega_\alpha(t-s)v(x,s)\,ds,\qquad
\omega_\alpha(t)=\frac{t^{\alpha-1}}{\Gamma(\alpha)}
\quad\text{for~$t>0$.}$$ If no sources or sinks are present, then the inhomogeneous term $f$ is identically zero. We assume for simplicity that the generalized diffusivity $K$ is a positive constant, and that the fractional PDE  holds for $x$ in a polygonal domain $\Omega\subseteq{\mathbb{R}}^2$ subject to homogeneous Dirichlet boundary conditions, with the initial condition $$\label{eq: ic}
u(x,0)=u_0(x)\quad\text{for $x\in\Omega$.}$$ In the limiting case when $\alpha\to1$, the fractional PDE  reduces to the classical heat equation that arises when the diffusing particles instead undergo Brownian motion.

Consider a spatial discretization of the preceding initial-boundary value problem using continuous piecewise-linear finite elements to obtain a semidiscrete solution $u_h$. The behaviour of $u_h$ is well understood if $\Omega$ is convex [@JinLazarovZhou2013; @McLeanThomee2010]: in this case, for general initial data $u_0\in L_2(\Omega)$ and an appropriate choice of $u_h(0)$, $$\|u_h(t)-u(t)\|\le Ct^{-\alpha}h^2\|u_0\|,\qquad0<t\le T,$$ whereas for smoother initial data $u_0\in H^2(\Omega)$, $$\|u_h(t)-u(t)\|\le Ch^2\|u_0\|_{H^2(\Omega)},\qquad0\le t\le T,$$ where  $\|\cdot\|=\|\cdot\|_{L_2(\Omega)}$. The error analysis establishing these bounds relies on the $H^2$-regularity property of the associated elliptic equation in $\Omega$, namely, that if $$\label{eq: elliptic bvp}
-K\nabla^2u=f\quad\text{in~$\Omega$,}\quad
\text{with $u=0$ on~$\partial\Omega$,}$$ then $u\in H^2(\Omega)$ with $\|u\|_{H^2(\Omega)}\le C\|f\|$.

In the present work, our aim is to study $u_h$ in the case when $\Omega$ is not convex. Since the above $H^2$-regularity breaks down, we can no longer expect $O(h^2)$ convergence if the finite element mesh is quasi-uniform. Our results generalize those of Chatzipantelidis, Lazarov, Thomée and Wahlbin [@ChatzipantelidisEtAl2006] for the heat equation (the limiting case $\alpha=1$) to the fractional-order case ($0<\alpha<1$). Our method of analysis relies on Laplace transformation, extending the approach of McLean and Thomee [@McLeanThomee2010] for the fractional order problem on a convex domain.

(0,0) – (0.4,0) arc \[radius=0.4, start angle=0, end angle=296.565\] – (0,0); (0,0) – (1,0) – (1,1) – (-1,1) – (-1,-1) – (0.5,-1) – (0,0); (-1.25,0) – (1.25,0); (0,-1.25) – (0,1.25); (0.4,-0.025) – (0.4,0.025); at (0.4,0) [$r_0$]{};

(0,0) circle \[radius=[0.015]{}\]; at (0,0) [$p_0=p_6$]{}; (1,0) circle \[radius=[0.015]{}\]; at (1,0) [$p_1$]{}; (1,1) circle \[radius=[0.015]{}\]; at (1,1) [$p_2$]{}; (-1,1) circle \[radius=[0.015]{}\]; at (-1,1) [$p_3$]{}; (-1,-1) circle \[radius=[0.015]{}\]; at (-1,-1) [$p_4$]{}; (0.5,-1) circle \[radius=[0.015]{}\]; at (0.5,-1) [$p_5$]{}; at (0.5,0) [$\Gamma_0$]{}; at (1,0.5) [$\Gamma_1$]{}; at (0,1) [$\Gamma_2$]{}; at (-1,0) [$\Gamma_3$]{}; at (-0.25,-1) [$\Gamma_4$]{}; at (0.25,-0.5) [$\Gamma_5$]{};

To focus on the essential difficulty, we assume that $\Omega$ has only a single re-entrant corner with angle $\pi/\beta$ for $1/2<\beta<1$. Without loss of generality, we assume that this corner is located at the origin and that, for some $r_0>0$, the intersection of $\Omega$ with the open disk $|x|<r_0$ is described in polar coordinates by $$\label{eq: 0 nbhd}
0<r<r_0\quad\text{and}\quad 0<\theta<\pi/\beta,$$ as illustrated in Figure \[fig: Omega\]. We denote the vertices of $\Omega$ by $p_0=(0,0)$, $p_1$, $p_2$, …, $p_J=p_0$, and the $j$th side by $$\Gamma_j=(p_j,p_{j+1})
    =\{\,(1-\sigma)p_j+\sigma p_{j+1}: 0<\sigma<1\,\}
    \quad\text{for $0\le j\le J-1$.}$$

Section \[sec: elliptic\] summarizes some key facts about the singular behaviour of the solution to the elliptic problem . In Section \[sec: FEM\], we describe a family of shape-regular triangulations ${\mathcal{T}}_h$ (indexed by the mesh parameter $h$) that depend on a local refinement parameter $\gamma\ge1$. The elements near the origin have sizes of order $h^\gamma$, so the ${\mathcal{T}}_h$ are quasi-uniform if $\gamma=1$ but become more highly refined with increasing $\gamma$. Our error bounds will be stated in terms of the quantity $$\label{eq: epsilon}
\epsilon(h,\gamma)=\begin{cases}
        
h^{\gamma\beta}/\sqrt{\gamma^{-1}-\beta},&1\le\gamma<1/\beta,\\
        h\sqrt{\log(1+h^{-1})},&\gamma=1/\beta,\\
        h/\sqrt{\beta-\gamma^{-1}},&\gamma>1/\beta,
\end{cases}$$ which ranges in size from $O(h^\beta)$ when $\gamma=1$ (the quasiuniform case) down to $O(h)$ when $\gamma>1/\beta$. We briefly review results for the finite element approximation of the elliptic problem, needed for our subsequent analysis: the error in $H^1(\Omega)$ is of order $\epsilon(h,\gamma)$, and the error in $L_2(\Omega)$ is of order $\epsilon(h,\gamma)^2$.

Section \[sec: time-dep\] gathers together some pertinant facts about the solution of the time-dependent problem  and its Laplace transform. Next, in Section \[sec: semidiscrete\], we introduce the semidiscrete finite element solution $u_h(t)$ of the time-dependent problem, and see that its stability properties mimic those of $u(t)$. In Section \[sec: error\] we study first the homogeneous equation (i.e., the case $f=0$), showing that the error in $L_2(\Omega)$ is of order $t^{-\alpha}\epsilon(h,\gamma)^2$ when $u_0\in L_2(\Omega)$. For smoother initial data, the $L_2$-error is of order $\epsilon(h,\gamma)^2$ uniformly for $0\le t\le T$. We also prove that for the inhomogeneous equation ($f\ne0$) with vanishing initial data ($u_0=0$), the error in $L_2(\Omega)$ is of order $t^{1-\alpha}\epsilon(h,\gamma)^2$. Thus, by choosing the mesh refinement parameter $\gamma>1/\beta$ we can restore second-order convergence in $L_2(\Omega)$. Section \[sec: alt bc\] outlines briefly how these results are affected by different choices of the boundary conditions. We conclude in Section \[sec: numerical\] with some numerical examples that illustrate our theoretical error bounds.

Singular behaviour in the elliptic problem {#sec: elliptic}
==========================================

In the weak formulation of the elliptic boundary-value problem  we introduce the Sobolev space $$V={\widetilde H}^1(\Omega)=H^1_0(\Omega)$$ and seek $u\in V$ satisfying $$a(u,v)={\langlef,v\rangle}\quad\text{for all $v\in V$,}$$ where $$\label{eq: a(u,v)}
a(u,v)=K\int_\Omega\nabla u\cdot\nabla v\,dx
\quad\text{and}\quad
{\langlef,v\rangle}=\int_\Omega fv\,dx.$$ Here, $f$ may belong to the dual space $V^*=H^{-1}(\Omega)$ if ${\langlef,v\rangle}$ is interpreted as the duality pairing on $V^*\times V$. Since $a(u,v)$ is bounded and coercive on ${\widetilde H}^1(\Omega)$, there exists a unique weak solution $u$, and $$\label{eq: Htilde u}
\|u\|_{{\widetilde H}^1(\Omega)}\le C\|f\|_{H^{-1}(\Omega)}.$$

To understand the difficulty created by the re-entrant corner, we separate variables in polar coordinates to construct the functions $$u_n^\pm(x)=r^{\pm n\beta}\sin(n\beta\theta)\quad
  \text{for $x=(r\cos\theta,r\sin\theta)$ and $n=1$, $2$, $3$, \dots,}$$ satisfying $$\label{eq: Laplacian in sector}
\nabla^2 u_n^\pm=0\quad
  \text{for $0<r<\infty$ and $0<\theta<\pi/\beta$,}$$ with $u_n^\pm=0$ if $\theta=0$ or $\theta=\pi/\beta$. Introducing a $C^\infty$ cutoff function $\eta$ such that $$\eta(x)=1\quad\text{for $|x|\le r_0/2$}
\quad\text{and}\quad
\eta(x)=0\quad\text{for $|x|\ge r_0$,}$$ we find that $$\eta u_n^+\in{\widetilde H}^1(\Omega)
\quad\text{but}\quad
\eta u_n^-\notin{\widetilde H}^1(\Omega)
\quad\text{for all $n\ge1$,}$$ and that $$\eta u_n^+\in H^2(\Omega)\quad\text{for all $n\ge2$,}
\quad\text{but}\quad
\eta u_1^+\notin H^2(\Omega).$$ Now consider the function $f=-K\nabla^2(\eta u_1^+)$. The choice of $\eta$ means that $f(x)=0$ for $|x|\le r_0/2$, and consequently $f$ is $C^\infty$ on $\overline{\Omega}$. Nevertheless, the (unique weak) solution of , namely $u=\eta u_1^+$, fails to belong to $H^2(\Omega)$.

Put $A=-K\nabla^2$ and $$\label{eq: V2 def}
V^2=H^2(\Omega)\cap{\widetilde H}^1(\Omega)
    =\{\,v\in H^2(\Omega):\text{$v=0$ on $\partial\Omega$}\,\}.$$

The bounded linear operator defined by the restriction $$A|_{V^2}:V^2\to L_2(\Omega)$$ is one-one and has closed range.

See Grisard [@Grisvard1992 Section 2.3].

Our task now is to identify the orthogonal complement in $L_2(\Omega)$ of the range $${\mathcal{R}}=\{\,f\in L_2(\Omega):
    \text{$f=Au$ for some $u\in V^2$}\,\}.$$ To this end, we define in the usual way the Hilbert space $$L_2(\Omega,A)=\{\,\phi\in L_2(\Omega): A\phi\in L_2(\Omega)\,\}$$ with the graph norm $\|\phi\|_{L_2(\Omega,A)}^2
=\|\phi\|^2+\|A\phi\|^2$. Let $\partial_n$ denote the outward normal derivative operator. It can be shown that the trace map $\phi\mapsto\bigl(\phi|_{\Gamma_j},\partial_n\phi|_{\Gamma_j}\bigr)$, has unique extensions from $C^1(\overline{\Omega})$ to bounded linear operators [@Grisvard1992 Theorems 1.4.2 and 1.5.2] $$H^2(\Omega)\to H^{3/2}(\Gamma_j)\times H^{1/2}(\Gamma_j)
\quad\text{and}\quad
L_2(\Omega,A)\to
{\widetilde H}^{-1/2}(\Gamma_j)\times{\widetilde H}^{-3/2}(\Gamma_j),$$ and that the second Green identity holds in the form [@Grisvard1992 Theorem 1.5.3] $$\int_\Omega\bigl[(Au)v-u(Av)\bigr]\,dx=\sum_{j=0}^{J-1}
    K\bigl[{\langleu,\partial_nv\rangle}_{\Gamma_j}
        -{\langle\partial_nu,v\rangle}_{\Gamma_j}\bigr]$$ for $u\in H^2(\Omega)$ and $v\in L_2(\Omega,A)$. Hence, $${\langleAu,\phi\rangle}={\langleu,A\phi\rangle}\quad
    \text{if $u\in V^2$, $\phi\in L_2(\Omega,A)$ and
    $\phi|_{\Gamma_j}=0$ for all $j$,}$$ implying that ${\mathcal{R}}$ is orthogonal in $L_2(\Omega)$ to the closed subspace $${\mathcal{N}}=\{\,\phi\in L_2(\Omega,A):
    \text{$A\phi=0$ in $\Omega$, 
    and $\phi|_{\Gamma_j}=0$ for every~$j$}\,\}.$$ Notice that ${\mathcal{N}}\cap{\widetilde H}^1(\Omega)=\{0\}$ because if $f=0$ then the unique weak solution of  in ${\widetilde H}^1(\Omega)$ is $u=0$.

The Hilbert space $L_2(\Omega)$ is the orthogonal direct sum of ${\mathcal{R}}$ and ${\mathcal{N}}$, and $\dim{\mathcal{N}}=1$ (assuming $\Omega$ has only a single re-entrant corner).

See [@Grisvard1992 Theorem 2.3.7].

Thus, given any $f\in L_2(\Omega)$, the unique weak solution $u\in{\widetilde H}^1(\Omega)$ of  belongs to $H^2(\Omega)$ if and only if $f\perp{\mathcal{N}}$. For general $f$, the following holds.

\[thm: singular term\] There exists $q\in{\mathcal{N}}$ (depending only on $\Omega$ and $\eta$) such that if $f\in L_2(\Omega)$ then the weak solution $u$ of  satisfies $u-{\langlef,q\rangle}\eta u_1^+\in V^2$ with $$\bigl\|u-{\langlef,q\rangle}\eta u_1^+\bigr\|_{H^2(\Omega)}\le C\|f\|.$$

Choose any nonzero $\phi\in{\mathcal{N}}$. Since $\eta u_1^+\in L_2(\Omega, A)$ but $\eta u_1^+\notin V^2$, we have ${\langleA(\eta u_1^+),\phi\rangle}\ne0$ and may therefore define $q=c\phi\in{\mathcal{N}}$ by letting $c=1/{\langleA(\eta u_1^+),\phi\rangle}$, so that ${\langleA(\eta u_1^+),q\rangle}=1$. Define $$u_1=u-{\langlef,q\rangle}\eta u_1^+\in{\widetilde H}^1(\Omega),$$ and observe that $u_1$ satisfies $Au_1=f_1$ where $f_1=f-{\langlef,q\rangle}A(\eta u_1^+)$. Since ${\langlef_1,q\rangle}=0$, the we deduce that $u_1\in H^2(\Omega)$ and $$\|u_1\|_{H^2(\Omega)}\le C\|f_1\|\le C\|f\|+C|{\langlef,q\rangle}|\le C\|f\|,$$ because $A(\eta u_1^+)\in C^\infty(\overline{\Omega})$ and $q\in L_2(\Omega)$.

Finite element approximation {#sec: FEM}
============================

Consider a family ${\mathcal{T}}_h$ of shape-regular triangulations of $\Omega$, indexed by the maximum element diameter $h$. For each element $\triangle\in{\mathcal{T}}_h$, let $$h_\triangle={\operatorname{diam}}(\triangle)\quad\text{and}\quad
r_\triangle={\operatorname{dist}}(0,\triangle),$$ and suppose that for some $\gamma\ge1$, $$\label{eq: h gamma}
chr_\triangle^{1-1/\gamma}\le h_\triangle
    \le Chr_\triangle^{1-1/\gamma}
    \quad\text{whenever $h^\gamma\le r_\triangle\le1$,}$$ with $$\label{eq: h gamma 0}
ch^\gamma\le h_\triangle\le Ch^\gamma
    \quad\text{whenever $r_\triangle\le h^\gamma$.}$$ Thus, if $\gamma=1$ then the mesh is globally quasiuniform, but for $\gamma>1$ the element diameter decreases from order $h$, when $r_\triangle\ge1$, to order $h^\gamma$, when $r_\triangle\le 
h^\gamma$. Such triangulations are widely used for elliptic problems on domains with re-entrant corners; see for instance Apel et al. [@ApelSaendigWhiteman1996 Section 3].

For each triangulation ${\mathcal{T}}_h$, we let $V_h$ denote the corresponding space of continuous piecewise-linear functions that vanish on $\partial\Omega$, so that $V_h\subseteq V={\widetilde H}^1(\Omega)$. Since the bilinear form  is bounded and coercive on $V$, there exists a unique finite element solution $u_h\in V_h$ defined by $$\label{eq: uh weak}
a(u_h,v)={\langlef,v\rangle}\quad\text{for all $v\in V_h$.}$$ This solution is stable in ${\widetilde H}^1(\Omega)$, $$\label{eq: Htilde uh}
\|u_h\|_{{\widetilde H}^1(\Omega)}\le C\|f\|_{H^{-1}(\Omega)},$$ and satisfies the quasi-optimal error bound $$\label{eq: quasioptimal}
\|u_h-u\|_{{\widetilde H}^1(\Omega)}
    \le C\min_{v\in V_h}\|v-u\|_{{\widetilde H}^1(\Omega)}.$$ Let $\Pi_h:C(\overline{\Omega})\to V_h$ denote the nodal interpolation operator, define the seminorm $$|v|_{m,\Omega}=\biggl(\sum_{j_1+j_2=m}\int_\Omega
    |\partial^jv(x)|^2\,dx\biggr)^{1/2},$$ where $\partial^j=\partial_{x_1}^{j_1}\partial_{x_2}^{j_2}$, and recall the standard interpolation error bounds [@Ciarlet2002] $$\label{eq: Pi_h}
|v-\Pi_hv|_{m,\triangle}\le Ch_\triangle^{2-m}|v|_{2,\triangle},
\qquad m\in\{0,1\}.$$ The next theorem reflects the influence of the singular behaviour of $u$ and the local mesh refinement parameter $\gamma$ on the accuracy of the approximation $u\approx\Pi_h u$.

\[thm: Pi\_h\] If $f\in L_2(\Omega)$ then the solution $u\in V$ of the elliptic problem  satisfies $$\|u-\Pi_h u\|\le Ch\epsilon(\gamma,h)\|f\|
\quad\text{and}\quad
\|u-\Pi_h u\|_{{\widetilde H}^1(\Omega)}\le C\epsilon(\gamma,h)\|f\|,$$ where $\epsilon(h,\gamma)$ is given by .

We use Theorem \[thm: singular term\] to split $u$ into singular and regular parts, $$u={u_{\operatorname{s}}}+{u_{\operatorname{r}}},\qquad {u_{\operatorname{s}}}={\langlef,q\rangle}\eta u_1^+,\qquad 
    {u_{\operatorname{r}}}\in H^2(\Omega),$$ with $\|{u_{\operatorname{r}}}\|_{H^2(\Omega)}\le C\|f\|$, leading to a corresponding decomposition of the interpolation error, $$u-\Pi_h u=({u_{\operatorname{s}}}-\Pi_h{u_{\operatorname{s}}})+({u_{\operatorname{r}}}-\Pi_h{u_{\operatorname{r}}}).$$ We see from  that $$\|{u_{\operatorname{r}}}-\Pi_h{u_{\operatorname{r}}}\|\le Ch^2|{u_{\operatorname{r}}}|_{2,\Omega}
    \le Ch^2\|f\|\le Ch\epsilon(h,\gamma)\|f\|$$ and $$|{u_{\operatorname{r}}}-\Pi_h{u_{\operatorname{r}}}|_{1,\Omega}\le Ch|{u_{\operatorname{r}}}|_{2,\Omega}
    \le Ch\|f\|\le C\epsilon(h,\gamma)\|f\|,$$ so it suffices to consider ${u_{\operatorname{s}}}-\Pi_h{u_{\operatorname{s}}}$. Note that $|\partial^j{u_{\operatorname{s}}}(x)|\le C\|f\||x|^{\beta-|j|}$ for any multi-index $j$, because $u_1^+$ is homogeneous of degree $\beta$.

We partition the triangulation into three subsets, $$\begin{gathered}
{\mathcal{T}}_h^1=\{\,\triangle\in{\mathcal{T}}_h:r_\triangle<h^\gamma\,\},\qquad
{\mathcal{T}}_h^2=\{\,\triangle\in{\mathcal{T}}_h:h^\gamma\le r_\triangle<1\,\},\\
{\mathcal{T}}_h^3=\{\,\triangle\in{\mathcal{T}}_h:r_\triangle\ge1\,\},\end{gathered}$$ and write $$|{u_{\operatorname{s}}}-\Pi_h{u_{\operatorname{s}}}|_{1,\Omega}^2=S_1+S_2+S_3
\quad\text{where}\quad
S_p=\sum_{\triangle\in{\mathcal{T}}_h^p}|{u_{\operatorname{s}}}-\Pi_h{u_{\operatorname{s}}}|_{1,\triangle}^2.$$ If $r_\triangle<h^\gamma$, then $|{u_{\operatorname{s}}}-\Pi_h{u_{\operatorname{s}}}|_{1,\triangle}\le|{u_{\operatorname{s}}}|_{1,\triangle}+
|\Pi_h{u_{\operatorname{s}}}|_{1,\triangle}$, and we estimate separately $$|{u_{\operatorname{s}}}|_{1,\triangle}^2\le C\|f\|^2\int_\triangle|x|^{2(\beta-1)}\,dx$$ and, using , $$|\Pi_h{u_{\operatorname{s}}}|_{1,\triangle}^2
    \le Ch_\triangle^{-2}|\Pi_hu_s|_{0,\triangle}^2
    \le Ch^{-2\gamma}\|f\|^2\int_\triangle|x|^{2\beta}\,dx.$$ Since $|x|\le r_\triangle+h_\triangle\le Ch^\gamma$ for $x\in\triangle$, we see that $$\label{eq: S1}
S_1\le C\|f\|^2\int_{|x|\le Ch^\gamma}|x|^{2(\beta-1)}\,dx
    +Ch^{-2\gamma}\|f\|^2\int_{|x|\le Ch^\gamma}|x|^{2\beta}\,dx
    \le Ch^{2\gamma\beta}\|f\|^2.$$ If $h^\gamma\le r_\triangle<1$, then gives $$|{u_{\operatorname{s}}}-\Pi_h{u_{\operatorname{s}}}|_{1,\triangle}^2\le Ch_\triangle^2\|f\|^2
    \int_\triangle|x|^{2(\beta-2)}\,dx,$$ and our assumption  on the mesh implies that for $x\in\triangle$, $$h_\triangle|x|^{\beta-2}
    \le Chr_\triangle^{1-1/\gamma}|x|^{\beta-2}
    =Ch\biggl(\frac{r_\triangle}{|x|}\biggr)^{1-1/\gamma}
    |x|^{\beta-1-1/\gamma}
    \le Ch|x|^{\beta-1-1/\gamma}$$ so $$\label{eq: S2}
\begin{aligned}
S_2&\le Ch^2\|f\|^2\int_{h^\gamma\le|x|\le 1+h}
        |x|^{2(\beta-1-1/\gamma)}\,dx\\
    &\le Ch^2\|f\|^2\int_{h^\gamma}^{1+h}r^{2(\beta-1/\gamma)-1}\,dr
    \le C\epsilon(h,\gamma)^2\|f\|^2.
\end{aligned}$$ In the remaining case $r_\triangle\ge1$, putting $R=\sup\{\,|x|:x\in\Omega\,\}$ we have $1\le|x|\le R$ for $x\in\triangle$, and thus $$\label{eq: S3}
S_3\le\sum_{\triangle\in{\mathcal{T}}_h^3} 
Ch_\triangle^2\|f\|^2\int_\triangle\,dx 
    \le Ch^2\|f\|^2\int_{1\le|x|\le R}\,dx\le Ch^2\|f\|^2.$$ Together, – show that $|{u_{\operatorname{s}}}-\Pi_h{u_{\operatorname{s}}}|_{1,\Omega}\le C\epsilon(h,\gamma)\|f\|$.

A similar argument shows $$\sum_{\triangle\in{\mathcal{T}}_h^1}|{u_{\operatorname{s}}}-I_h{u_{\operatorname{s}}}|_{0,\triangle}^2
    \le C\|f\|^2\int_0^{Ch^\gamma}h^{2\gamma\beta}r\,dr
    \le Ch^{2\gamma(\beta+1)}\|f\|^2$$ and $$\begin{aligned}
\sum_{\triangle\in{\mathcal{T}}_h^2\cup{\mathcal{T}}_h^3}
    |{u_{\operatorname{s}}}-I_h{u_{\operatorname{s}}}|_{0,\triangle}^2
    &\le Ch^4\|f\|^2\int_{h^\gamma}^{R}
        r^{2(\beta-1/\gamma)-1}r^{2(1-1/\gamma)}\,dr
    \le Ch^2\epsilon(h,\gamma)^2\|f\|^2.\end{aligned}$$ Hence, $\|{u_{\operatorname{s}}}-\Pi_h{u_{\operatorname{s}}}\|\le Ch\epsilon(h,\gamma)\|f\|$ and the desired bounds follow.

\[thm: elliptic error\] If $f\in L_2(\Omega)$ then the finite element solution $u_h\in V_h$ of the elliptic problem  satisfies $$\|u_h-u\|\le C\epsilon(\gamma,h)^2\|f\|
\quad\text{and}\quad
\|u_h-u\|_{{\widetilde H}^1(\Omega)}\le C\epsilon(\gamma,h)\|f\|,$$ where $\epsilon(h,\gamma)$ is given by .

The bound in ${\widetilde H}^1(\Omega)$ follows at once from  and Theorem \[thm: Pi\_h\]. The error bound in $L_2(\Omega)$ is proved via the usual duality argument. In fact, given any $\phi\in L_2(\Omega)$, the dual variational problem $$a(w,\psi)={\langlew,\phi\rangle}\quad\text{for all $w\in{\widetilde H}^1(\Omega)$,}$$ has a unique solution $\psi\in{\widetilde H}^1(\Omega)$. Since the bilinear form $a$ is symmetric, the preceding estimate for $u-\Pi_hu$ carries over, with $\phi$ playing the role of $f$, to yield $\|\psi-\Pi_h\psi\|_{{\widetilde H}^1(\Omega)}
\le C\epsilon(h,\gamma)\|\phi\|$. Thus, $$\begin{aligned}
|{\langleu_h-u,\phi\rangle}|&=|a(u_h-u,\psi)|=|a(u_h-u,\psi-\Pi_h\psi)|\\
    &\le C\|u_h-u\|_{{\widetilde H}^1(\Omega)}
    \|\psi-\Pi_h\psi\|_{{\widetilde H}^1(\Omega)}
    \le C\epsilon(h,\gamma)^2\|f\|\|\phi\|,\end{aligned}$$ implying that $\|u_h-u\|\le C\epsilon(h,\gamma)^2\|f\|$.

The time-dependent problem {#sec: time-dep}
==========================

We may view $A=-K\nabla^2$ as an unbounded operator on $L_2(\Omega)$ with domain $V^2$ given by . Since the associated bilinear form  is symmetric and coercive, and since the inclusion ${\widetilde H}^1(\Omega)\subseteq L_2(\Omega)$ is compact, there exists a complete orthonormal sequence of eigenfunctions $\phi_1$, $\phi_2$, $\phi_3$, …and corresponding real eigenvalues $\lambda_1$, $\lambda_2$, $\lambda_3$, …with $\lambda_j\to\infty$ as $j\to\infty$. Thus, $$A\phi_n=\lambda_n\phi_n\quad\text{and}\quad
{\langle\phi_m,\phi_n\rangle}=\delta_{mn}\quad
\text{for all $m$, $n\in\{1,2,3,\ldots\}$,}$$ and we may assume that $0<\lambda_1\le\lambda_2\le\lambda_3\le\cdots$. Moreover, $$(zI-A)^{-1}:L_2(\Omega)\to L_2(\Omega)$$ is a bounded linear operator for each complex number $z$ not in the spectrum ${\operatorname{spec}}(A)=\{\,\lambda_1,\lambda_2,\lambda_3,\ldots\}$, and given any $\theta_0\in(0,\pi)$ we have a resolvent estimate in the induced operator norm [@LeGiaMcLean2011 Lemma 1], $$\label{eq: resolvent}
\|(zI-A)^{-1}\|\le\frac{1+2/\lambda_1}{\sin\theta_0}\,\frac{1}{1+|z|}
    \quad\text{for $|\arg z|>\theta_0$.}$$

Define the Laplace transform $\hat f={\mathcal{L}}f$ of a suitable $f:[0,\infty)\to L_2(\Omega)$ by $$\hat f(z)=({\mathcal{L}}f)(z)={\mathcal{L}}\{f(t)\}_{t\to z}
  =\int_0^\infty e^{-zt}f(t)\,dt,$$ for $\Re z$ sufficiently large. Since ${\mathcal{L}}\{\partial_t^{1-\alpha}f\}_{t\to z}=z^{1-\alpha}\hat f(z)$, a formal calculation implies that the fractional diffusion equation  transforms to an elliptic problem (with complex coefficients) for $\hat u(z)$, $$z\hat u(z)+z^{1-\alpha}A\hat u(z)=u_0+\hat f(z),$$ and so $$\label{eq: u hat}
\hat u(z)=z^{\alpha-1}(z^\alpha I+A)^{-1}\bigl(u_0+\hat f(z)\bigr).$$ The boundary condition $u(t)=0$ on $\partial\Omega$ tranforms to give $\hat u(z)=0$ on $\partial\Omega$. Using ${\mathcal{L}}\{t^{p\alpha}/\Gamma(1+p\alpha)\}_{t\to z}=z^{-1-p\alpha}$ we find that for $\lambda>0$ and $|z|>\lambda^{-1/\alpha}$, $$z^{\alpha-1}(z^\alpha+\lambda)^{-1}
    =z^{-1}\sum_{p=0}^\infty\bigl(-\lambda z^{-\alpha}\bigr)^p
    ={\mathcal{L}}\biggl\{\sum_{p=0}^\infty
        \frac{(-\lambda t^\alpha)^p}{\Gamma(1+p\alpha)}
        \biggr\}_{t\to z}
    ={\mathcal{L}}\{E_\alpha(-\lambda t^\alpha)\},$$ where $E_\alpha$ is the Mittag–Leffler function. The inequalities $0\le E_\alpha(-t)\le1$ for $0\le t<\infty$ imply that the sum $$\label{eq: E(t)}
{\mathcal{E}}(t)v=\sum_{n=1}^\infty E_\alpha(-\lambda_nt^\alpha)
    {\langlev,\phi_n\rangle}\phi_n$$ defines a bounded linear operator ${\mathcal{E}}(t):L_2(\Omega)\to L_2(\Omega)$ satisfying $$\label{eq: ||E(t)||}
\|{\mathcal{E}}(t)v\|\le\|v\|\quad\text{for $0\le t<\infty$.}$$ Thus, for each eigenfunction $\phi_n$, $$z^{\alpha-1}(z^\alpha I+A)^{-1}\phi_n
    =z^{\alpha-1}(z^\alpha+\lambda_n)^{-1}\phi_n
    ={\mathcal{L}}\{E_\alpha(-\lambda_nt^\alpha)\phi_n\}_{t\to z},$$ and we conclude that $$\label{eq: E hat}
{\widehat{\mathcal{E}}}(z)=z^{\alpha-1}(z^\alpha I+A)^{-1}.$$ Writing  as $\hat u(z)={\widehat{\mathcal{E}}}(z)u_0+{\widehat{\mathcal{E}}}(z)\hat f(z)$ then yields a Duhamel formula, $$\label{eq: u(t) mild}
u(t)={\mathcal{E}}(t)u_0+\int_0^t{\mathcal{E}}(t-s)f(s)\,ds\quad\text{for $t>0$,}$$ that serves to define the *mild solution* of our initial-boundary value problem for . In particular, ${\mathcal{E}}(t)$ is the solution operator for the homogeneous problem ($f\equiv0$) with initial data $u_0\in L_2(\Omega)$. Also, the bound  immediately implies a stability estimate in $L_2(\Omega)$ for a general, locally integrable $f:[0,\infty)\to L_2(\Omega)$, namely $$\|u(t)\|\le\|u_0\|+\int_0^t\|f(s)\|\,ds\quad\text{for $t>0$.}$$

The semidiscrete finite element solution {#sec: semidiscrete}
========================================

Let $P_h$ denote the orthoprojector $L_2(\Omega)\to V_h$, that is, $P_hv\in V_h$ satisfies ${\langleP_hv,w\rangle}={\langlev,w\rangle}$ for all $v\in L_2(\Omega)$ and $w\in V_h$. There exists a unique linear operator $A_h:V_h\to V_h$ such that $${\langleA_hv,w\rangle}=a(v,w)\quad\text{for all $v$, $w\in V_h$,}$$ and the operator equation $A_hu_h=P_hf$ is equivalent to the variational equation  used to define the finite element solution $u_h\in V_h$ of the elliptic problem . Denote the number of degrees of freedom by $N=\dim V_h$ and equip $V_h$ with the norm induced from $L_2(\Omega)$. The finite element space $V_h$ has an orthonormal basis of eigenfunctions $\Phi_1$, $\Phi_2$, …, $\Phi_N$ with corresponding real eigenvalues $\Lambda_1$, $\Lambda_2$, …, $\Lambda_N$. Thus, $$A_h\Phi_n=\Lambda_n\Phi_n\quad\text{and}\quad
{\langle\Phi_m,\Phi_n\rangle}=\delta_{mn}\quad
\text{for $m$, $n\in\{1,2,\ldots,N\}$,}$$ and we may assume that $0<\Lambda_1\le\Lambda_2\le\cdots\le\Lambda_N$. Moreover, the resolvent $$(zI-A_h)^{-1}:V_h\to V_h$$ exists for every $z\notin{\operatorname{spec}}(A_h)=\{\Lambda_1,\Lambda_2,\ldots,\Lambda_N\}$, and we have the an estimate corresponding to : $$\label{eq: resolvent h}
\|(zI-A_h)^{-1}\|\le\frac{1+2/\Lambda_1}{\sin\theta_0}
    \,\frac{1}{1+|z|}\quad\text{for $|\arg z|>\theta_0$.}$$ Note that $\lambda_1\le\Lambda_1$ so this bound is uniform in $h$.

The first Green identity yields the variational formulation for , $$\label{eq: u(t) variational}
{\langle\partial_tu,v\rangle}+a(\partial_t^{1-\alpha}u,v)={\langlef(t),v\rangle}
    \quad\text{for all $v\in{\widetilde H}^1(\Omega)$ and $t>0$,}$$ so we define the finite element solution $u_h:[0,\infty)\to V_h$ by $$\label{eq: uh(t) variational}
{\langle\partial_tu_h,v\rangle}+a(\partial_t^{1-\alpha}u_h,v)={\langlef(t),v\rangle}
    \quad\text{for all $v\in V_h$ and $t>0$,}$$ with $u_h(0)=u_{0h}$, where $u_{0h}\in V_h$ is a suitable approximation to the initial data $u_0$. Thus, the vector of nodal values ${\mathbf{U}}(t)$ satisfies the integro-differential equation in ${\mathbb{R}}^N$, $$\label{eq: MOL}
{\mathbf{M}}\partial_t{\mathbf{U}}+{\mathbf{S}}\partial_t^{1-\alpha}
{\mathbf{U}}={\mathbf{F}}(t) ,$$ where ${\mathbf{M}}$ and ${\mathbf{S}}$ denote the $N\times N$ mass and stiffness matrices, respectively, and ${\mathbf{F}}(t)$ denotes the load vector. In the limiting case as $\alpha\to1$, when  becomes the heat equation, we see that reduces to the usual system of (stiff) ODEs arising in the method of lines.

The variational equation  is equivalent to $$\partial_tu_h+\partial_t^{1-\alpha}A_hu_h=P_hf(t)
    \quad\text{for $t>0$.}$$ Taking Laplace transforms as in Section \[sec: time-dep\], we find that $$\hat u_h(z)=z^{\alpha-1}(z^\alpha I+A_h)^{-1}
    \bigl(u_{0h}+P_h\hat f(z)\bigr)$$ and thus $$u_h(t)={\mathcal{E}}_h(t)u_{0h}+\int_0^t{\mathcal{E}}_h(t-s)P_hf(s)\,ds
    \quad\text{for $t>0$,}$$ where $${\mathcal{E}}_h(t)v=\sum_{n=1}^N E_\alpha(-\Lambda_n t^\alpha)
    {\langlev,\Phi_n\rangle} \Phi_n.$$ In the same way as  we have $$\label{eq: ||Eh(t)||}
\|{\mathcal{E}}_h(t)v\|\le\|v\|\quad\text{for $t>0$~and $v\in v_h$,}$$ implying that the finite element solution is stable in $L_2(\Omega)$, $$\label{eq: uh stable}
\|u_h(t)\|\le\|u_{0h}\|+\int_0^t\|f(s)\|\,ds.$$

For convenience, we put $$B(z)=(z^\alpha I+A)^{-1}
\quad\text{and}\quad
B_h(z)=(z^\alpha I+A_h)^{-1},$$ which satisfy the following bounds.

\[lem: B(z)\] If $|\arg z^\alpha|<\pi-\theta_0$, then

1.  $\|B(z)v\|\le C\|v\|/(1+|z|^\alpha)$ and $\|B(z)v\|_{{\widetilde H}^1(\Omega)}\le C\|v\|$ for $v\in L_2(\Omega)$;

2.  $\|B_h(z)v\|\le C\|v\|/(1+|z|^\alpha)$ and $\|B_h(z)v\|_{{\widetilde H}^1(\Omega)}\le C\|v\|$ for $v\in V_h$.

First let $v\in L_2(\Omega)$. The resolvent estimate  immediately implies the desired bounds for $w(z)=B(z)v$ in $L_2(\Omega)$. To estimate the norm of $w(z)$ in ${\widetilde H}^1(\Omega)$, observe that $Aw(z)=v-z^\alpha w(z)$ so by , $$\|w(z)\|_{{\widetilde H}^1(\Omega)}\le C\|v-z^\alpha w(z)\|_{H^{-1}(\Omega)}
    \le C\|v-z^\alpha B(z)v\|\le C\|v\|.$$ When $v\in V_h$, the estimates for $B_h(z)v$ follow in the same way from  and .

Since ${\mathcal{L}}\{{\mathcal{E}}(t)\phi_n\}_{t\to z}
=z^{\alpha-1}(z^\alpha+\lambda_n)^{-1}\phi_n=z^{\alpha-1}B(z)\phi_n$, the Laplace inversion formula implies that $$\begin{aligned}
{\mathcal{E}}(t)\phi_n&=\lim_{M\to\infty}\biggl(\frac{1}{2\pi i} 
  \int_{1-iM}^{1+iM} 
  e^{zt}z^{\alpha-1}(z^\alpha+\lambda_n)^{-1}\,dz\biggr)\phi_n
  =\frac{1}{2\pi i}\int_\Gamma e^{zt}z^{\alpha-1}B(z)\phi_n\,dz,\end{aligned}$$ for $t>0$ and for any contour $\Gamma$ that begins at $\infty$ in the third quadrant, ends at $\infty$ in the second quadrant and avoids the negative real axis. The factor $e^{zt}$ is exponentially small as $\Re z\to-\infty$, so and Lemma \[lem: B(z)\] ensure that $$\label{eq: E integral}
{\mathcal{E}}(t)v=\frac{1}{2\pi i}\int_\Gamma e^{zt}z^{\alpha-1}B(z)v\,dz
\quad\text{for $t>0$ and $v\in L_2(\Omega)$,}$$ where the integral over $\Gamma$ is absolutely convergent in $L_2(\Omega)$.

Likewise, ${\mathcal{L}}\{{\mathcal{E}}_h(t)\Phi_n\}_{t\to z}
=z^{\alpha-1}(z^\alpha+\Lambda_n)^{-1}\Phi_n=z^{\alpha-1}B_h(z)\Phi_n$ and we have a corresponding integral representation $$\label{eq: Eh integral}
{\mathcal{E}}_h(t)v=\frac{1}{2\pi i}\int_\Gamma e^{zt}z^{\alpha-1}
    B_h(z)v\,dz\quad\text{for $t>0$ and $v\in V_h$.}$$

Error bounds {#sec: error}
============

The homogeneous equation
------------------------

We now consider the error $u_h(t)-u(t)$ in the case $f\equiv0$. The main difficulty will be to estimate the difference $$\label{eq: E(t) error}
{\mathcal{E}}_h(t)P_hu_0-{\mathcal{E}}(t)u_0=\frac{1}{2\pi i}\int_\Gamma 
    e^{zt}z^{\alpha-1}G_h(z)u_0\,dz,$$ where, by  and , $G_h(z)=B_h(z)P_h-B(z)$. We begin by estimating $G_h(z)v$.

\[lem: Gh(z)\] If $v\in L_2(\Omega)$ and $|\arg z^\alpha|<\pi-\theta_0$, then $$\|G_h(z)v\|\le C\epsilon(h,\gamma)^2\|v\|
\quad\text{and}\quad
\|G_h(z)v\|_{{\widetilde H}^1(\Omega)}\le 
C(1+|z|^\alpha)\epsilon(h,\gamma)\|v\|.$$

Given $v\in L_2(\Omega)$, let $w(z)=B(z)v\in V$ so that $Aw(z)=v-z^\alpha w(z)$, and let $w_h(z)\in V_h$ be the solution of $$A_hw_h(z)=P_h[v-z^\alpha w(z)].$$ In this way, $P_hv=z^\alpha P_h w(z)+A_hw_h(z)=(z^\alpha I+A_h)w_h(z)
-z^\alpha[w_h(z)-P_hw(z)]$, and thus $B_h(z)P_hv=w_h(z)-z^\alpha B_h(z)[w_h(z)-P_hw(z)]$, implying that $$\label{eq: Gh(z)}
G_h(z)v=w_h(z)-w(z)-z^\alpha B_h(z)P_h[w_h(z)-w(z)].$$ Lemma \[lem: B(z)\] shows that $\|w(z)\|\le C\|v\|/(1+|z|^\alpha)$ so $\|v-z^\alpha w(z)\|\le C\|v\|$. By applying Theorem \[thm: elliptic error\], with $w(z)$ and $v-z^\alpha w(z)$ playing the roles of $u$ and $f$, respectively, we deduce that $$\|w_h(z)-w(z)\|\le C\epsilon(h,\gamma)^2\|v\|
\quad\text{and}\quad
\|w_h(z)-w(z)\|_{{\widetilde H}^1(\Omega)}\le C\epsilon(h,\gamma)\|v\|.$$ The result now follows from  after another application of Lemma \[lem: B(z)\].

\[thm: non-smooth u0\] Assume that $f\equiv0$ and $u_0\in L_2(\Omega)$. Then the mild solution $u(t)={\mathcal{E}}(t)u_0$ and its finite element approximation $u_h(t)={\mathcal{E}}_h(t)u_{0h}$ satisfy the error bounds $$\|u_h(t)-u(t)\|\le\|u_{0h}-P_hu_0\|
    +Ct^{-\alpha}\epsilon(h,\gamma)^2\|u_0\|$$ and $$\|u_h(t)-u(t)\|_{{\widetilde H}^1(\Omega)}
    \le Ct^{-\alpha}\|u_{0h}-P_hu_0\|
    +C(t^{-2\alpha}+t^{-\alpha})\epsilon(h,\gamma)\|u_0\|$$ for $t>0$, where $\epsilon(h,\gamma)$ is given by .

We split the error into two terms, $$u_h(t)-u(t)={\mathcal{E}}_h(t)\bigl(u_{0h}-P_hu_0\bigr)
    +\bigl[{\mathcal{E}}_h(t)P_hu_0-{\mathcal{E}}(t)u_0\bigr].$$ It follows from  that $\bigl\|{\mathcal{E}}_h(t)\bigl(u_{0h}-P_hu_0\bigr)\bigr\|
\le\|u_{0h}-P_hu_0\|$, and to estimate the second term we use the integral representation  with $\Gamma=\Gamma_+-\Gamma_-$, where $\Gamma_\pm$ is the contour $z=se^{\pm i3\pi/4}$ for $0<s<\infty$. Applying Lemma \[lem: Gh(z)\] and making an obvious substitution, we find that $$\begin{aligned}
\bigl\|{\mathcal{E}}_h(t)P_hu_0-{\mathcal{E}}(t)u_0\bigr\|&\le C\epsilon(h,\gamma)^2\|u_0\|
    \int_0^\infty e^{-st/\sqrt{2}}s^\alpha\,\frac{ds}{s}\\
    &=C\epsilon(h,\gamma)^2\|u_0\|t^{-\alpha}
    \int_0^\infty e^{-s/\sqrt{2}}s^\alpha\,\frac{ds}{s},\end{aligned}$$ which proves the first error bound of the theorem.

Choosing $\Gamma=\Gamma_+-\Gamma_-$ in the integral representation  of ${\mathcal{E}}_h(t)v$, and using Lemma \[lem: B(z)\], we have for $v\in V_h$, $$\|{\mathcal{E}}_h(t)v\|_{{\widetilde H}^1(\Omega)}\le C\int_0^\infty e^{-st/\sqrt{2}}
    s^\alpha\|v\|\,\frac{ds}{s}
    =Ct^{-\alpha}\|v\|\int_0^\infty e^{-s/\sqrt{2}}s^\alpha\,
    \frac{ds}{s}\le Ct^{-\alpha}\|v\|,$$ so in particular, when $v=u_{0h}-P_hu_0$, $$\bigl\|{\mathcal{E}}_h(t)\bigl(u_{0h}-P_hu_0\bigr)\bigr\|_{{\widetilde H}^1(\Omega)}
\le Ct^{-\alpha}\|u_{0h}-P_hu_0\|.$$ Finally, using  and Lemma \[lem: Gh(z)\] again, $$\begin{gathered}
\bigl\|{\mathcal{E}}_h(t)P_hu_0-{\mathcal{E}}(t)u_0\bigr\|_{{\widetilde H}^1(\Omega)}
    \le C\epsilon(h,\gamma)\|u_0\|
    \int_0^\infty e^{-st/\sqrt{2}}s^\alpha(1+s^\alpha)\,\frac{ds}{s}\\
    =C\epsilon(h,\gamma)\|u_0\|t^{-2\alpha}\int_0^\infty 
    e^{-s/\sqrt{2}}s^\alpha(t^\alpha+s^\alpha)\,\frac{ds}{s}
    \le C\epsilon(h,\gamma)(t^{-\alpha}+t^{-2\alpha})\|u_0\|,\end{gathered}$$ proving the second error estimate of the theorem.

When $u_0$ is sufficiently regular we obtain an error bound that is uniform in $t$. The proof uses the Ritz projector $R_h:{\widetilde H}^1(\Omega)\to V_h$, defined by $$\label{eq: Ritz R_h}
a(R_hu,v)=a(u,v)\quad\text{for all $v\in V_h$,}$$ and relies on the regularity estimate [@McLean2010 Theorem 4.4] $$\label{eq: Au_t reg}
t^\alpha\|A\partial_tu\|\le Ct^{\sigma\alpha-1}\|A^\sigma u_0\|
\quad\text{for $0<t\le T$ and $0\le\sigma\le2$.}$$

\[thm: u0 smooth\] Assume that $f\equiv0$ and $0<\delta\le1$. If $A^{1+\delta}u_0\in L_2(\Omega)$, then $$\|u_h(t)-u(t)\|\le\|u_{0h}-R_hu_0\| 
    +C\delta^{-1}\epsilon(h,\gamma)^2\|A^{1+\delta}u_0\|
    \quad\text{for $0\le t\le T$.}$$

To begin with, we permit $f\ne0$, and in the usual way, decompose the error as $$u_h(t)-u(t)=\vartheta(t)+\varrho(t),$$ where $\vartheta(t)=u_h(t)-R_hu(t)\in V_h$ and $\varrho(t)=R_hu(t)-u(t)$. Furthermore, since $\varrho(t)=\varrho(0)+\int_0^t\partial_t\varrho(s)\,ds$ it follows that $$\|\varrho(t)\|\le\|\varrho(0)\|+\int_0^t\|\partial_t\varrho(s)\|\,ds.$$ By , if $v\in V_h$ then $${\langle\partial_t\vartheta,v\rangle}+a(\partial_t^{1-\alpha}\vartheta,v)
={\langlef,v\rangle}-{\langle\partial_tR_hu,v\rangle}-a(\partial_t^{1-\alpha}R_hu,v),$$ and using the definition of the Ritz projector  followed by , we have $$a(\partial_t^{1-\alpha}R_hu,v)
=a(\partial_t^{1-\alpha}u,v)={\langlef,v\rangle}-{\langle\partial_tu,v\rangle},$$ so $$\label{eq: theta}
{\langle\partial_t\vartheta,v\rangle}+a(\partial_t^{1-\alpha}\vartheta,v)
=-{\langle\partial_t\varrho,v\rangle}.$$ In other words, $\vartheta$ is the finite element solution of the fractional diffusion problem with source term $-\partial_t\varrho(t)$. Thus, the stability estimate  gives $$\|\vartheta(t)\|\le\|\vartheta(0)\|
    +\int_0^t\|\partial_t\varrho(s)\|\,ds.$$ Theorem \[thm: elliptic error\] implies that $\|v-R_hv\|\le C\epsilon(h,\gamma)^2\|Av\|$, so $$\begin{aligned}
\|u_h(t)-u(t)\|&\le\|\vartheta(0)\|+\|\varrho(0)\|+2\int_0^t
    \|\partial_t\varrho(s)\|\,ds\\
    &\le\|u_{0h}-R_hu_0\|+C\epsilon(h,\gamma)^2\biggl(\|Au_0\|
+\int_0^t\|A\partial_tu(s)\|\,ds\biggr),\end{aligned}$$ and, assuming now that $f\equiv0$, we use  to bound the integral by $C\delta^{-1}\|A^{1+\delta}u_0\|$.

Intermediate regularity of $u_0$ ensures a milder growth of the error as $t\to0$.

Assume that $f\equiv0$ and $0<\delta\le1$. If $0<\theta<1$ and $u_{0h}=P_hu_0$, then $$\|u_h(t)-u(t)\|\le C\delta^{-\theta}t^{-\alpha(1-\theta)}
\epsilon(h,\gamma)^2\|A^{(1+\delta)\theta}u_0\|
    \quad\text{for $0<t\le T$.}$$

Since $\|u_{0h}-R_hu_0\|=\|P_h(u_0-R_hu_0)\|\le\|u_0-R_hu_0\|
    \le C\epsilon(h,\gamma)^2\|Au_0\|$, we see that $$\|u_h(t)-u(t)\|\le Ct^{-\alpha}\epsilon(h,\gamma)^2\|u_0\|$$ and $$\|u_h(t)-u(t)\|
    \le C\delta^{-1}\epsilon(h,\gamma)^2\|A^{1+\delta}u_0\|.$$ By interpolation, $\|u_h(t)-u(t)\|\le C\bigl(
t^{-\alpha}\epsilon(h,\gamma)^2\bigr)^{1-\theta}
\bigl(\delta^{-1}\epsilon(h,\gamma)^2\bigr)^\theta
\|A^{(1+\delta)\theta}u_0\|$.

The inhomogeneous equation
--------------------------

When $u_0=0$ and $f\ne0$, we readily adapt the proof of Theorem \[thm: u0 smooth\] to show the following error bound; see also McLean and Thomée [@McLeanThomee2010 Lemma 4.1]. Instead of , we now rely on the regularity result [@McLean2010 Theorem 4.1] $$\label{eq: Au reg}
\|A{\mathcal{E}}(t)v\|\le Ct^{-\alpha}\|v\|\quad\text{for $0<t\le T$.}$$

\[thm: u0=0\] If $u_0=u_{0h}=0$ then $$\|u_h(t)-u(t)\|\le Ct^{1-\alpha}\epsilon(h,\gamma)^2\biggl(
    \|f(0)\|+\int_0^t\|\partial_tf(s)\|\,ds\biggr)
    \quad\text{for $t>0$.}$$

The equation  for $\vartheta$ holds for a general $f$, and now $\vartheta(0)=\varrho(0)=0$, so stability of the finite element solution implies that $$\|u_h(t)-u(t)\|\le2\int_0^t\|\partial_t\varrho(s)\|\,ds
    \le C\epsilon(h,\gamma)^2\int_0^t\|A\partial_tu(s)\|\,ds.$$ Since $$Au(s)=A\int_0^s{\mathcal{E}}(s-\tau)f(\tau)\,d\tau
    =\int_0^sA{\mathcal{E}}(\tau)f(s-\tau)\,d\tau,$$ we have $$A\partial_tu(s)=A{\mathcal{E}}(s)f(0)
    +\int_0^sA{\mathcal{E}}(\tau)\partial_tf(s-\tau)\,d\tau,$$ so using , $$\|A\partial_tu(s)\|\le Cs^{-\alpha}\|f(0)\|
    +\int_0^s \tau^{-\alpha}\|\partial_tf(s-\tau)\|\,d\tau.$$ Thus, $$\int_0^t\|A\partial_tu(s)\|\,ds\le Ct^{1-\alpha}\|f(0)\|
+\int_0^t\int_0^s\tau^{-\alpha}\|\partial_tf(s-\tau)\|\,d\tau\,ds$$ and the double integral equals $$\int_0^t\int_0^s(s-\tau)^{-\alpha}\|\partial_tf(\tau)\|\,d\tau\,ds
    =\int_0^t\|\partial_tf(\tau)\|(t-\tau)^{1-\alpha}\,d\tau,$$ implying the desired estimate.

Alternative boundary conditions {#sec: alt bc}
===============================

Neumann boundary conditions
---------------------------

Separation of variables in polar coordinates yields the functions $$u_n^\pm=r^{\pm n\beta}\cos(n\beta\theta)
\quad\text{for $n=1$, $2$, $3$, \dots}$$ satisfying with $\partial_\theta 
u_n^\pm=0$ if $\theta=0$ or $\theta=\pi/\beta$. In addition, for $n=0$ we find $u_0^+=1$ and $u_0^-=\log r$, and can readily check that $$\label{eq: un H1}
\eta u_n^+\in H^1(\Omega)
\quad\text{but}\quad
\eta u_n^-\notin H^1(\Omega)
\quad\text{for all $n\ge0$,}$$ and that $\eta u^+_n\in H^2(\Omega)$ iff $n\ne1$. If we impose a homogeneous Neumann boundary condition $\partial_n 
u=0$ on $\partial\Omega$, then our results are essentially unchanged, but the fact that $A=-K\nabla^2$ now possesses a zero eigenvalue complicates the analysis [@MustaphaMcLean2011 Section 4].

Mixed boundary conditions
-------------------------

The functions $$u_n^\pm=r^{(n-\tfrac12)\beta}\sin(n-\tfrac12)\beta\theta
\quad\text{for $n=1$, $2$, $3$, \dots}$$ satisfy with $u_n^\pm=0$ if $\theta=0$ and $\partial_\theta u_n^\pm=0$ if $\theta=\pi/\beta$. Once again, holds, for all $n\ge1$, however now we have $$\eta u_n^+\in H^2(\Omega)\quad\text{for all $n\ge3$,}
\quad\text{but}\quad
\eta u_1^+, \eta u_2^+\notin H^2(\Omega),$$ assuming $1/2<\beta<1$. A new feature is that $\eta u_1^+\notin 
H^2(\Omega)$ also when $1\le\beta<2$, that is, for an interior angle between $\pi/2$ and $\pi$, in which case $\Omega$ is in fact convex. The proof of Theorem 3.1 must be modified by replacing $\beta$ with $\beta/2$, and replacing $\epsilon(h,\gamma)$ with $$\label{eq: epsmix}
{\epsilon_{\text{mix}}}(h,\gamma)=\begin{cases}
    h^{\gamma\beta/2}/\sqrt{\gamma^{-1}-\beta/2},&
        1\le\gamma<2/\beta,\\
    h\sqrt{\log(1+h^{-1})},&\gamma=2/\beta,\\
    h/\sqrt{\beta/2-\gamma^{-1}},&\gamma>2/\beta,
\end{cases}$$ provided the interior angles at $p_1$, $p_2$, …, $p_{J-1}$ are all less than or equal to $\pi/2$. We may then proceed as for Dirichlet boundary conditions (since all the eigenvalues of $A$ are strictly positive), with $\epsilon(h,\gamma)$ replaced by ${\epsilon_{\text{mix}}}(h,\gamma)$ in our error estimates.

Numerical experiments {#sec: numerical}
=====================

We consider two problems posed on a domain of the form $$\Omega=\{\,(r\cos\theta,r\sin\theta):
    \text{$0<r<1$ and $0<\theta<\pi/\beta$}\,\},$$ with $\beta=2/3$. Although $\Omega$ is not a polygon, the additional error in $u_h$ due to approximation of the curved part of $\partial\Omega$ is of order $h^2$ in $L_2(\Omega)$, and hence our error bounds should apply unchanged. To fix the time scale for the solutions of the fractional diffusion equation , we choose the generalized diffusivity $K$ so that the smallest eigenvalue of $A=-K\nabla^2$ equals $1$. Figure \[fig: meshes\] shows two successive meshes out of a sequence satisfying our assumptions  and for $\gamma=1/\beta=3/2$; notice that these meshes are not nested. The mesh generation code takes a specified $h_*$ and $\gamma$ and produces a triangulation with maximum element diameter $h$ equivalent to $h_*$. All source files were written in Julia 0.4 [@Julia] with some calls to Gmsh 2.10.1 [@Gmsh], and all computations performed on a desktop PC with 16GB of RAM and an Intel Core i7-4770 CPU.

![Meshes with $h_*=2^{-3}$ (left) and $h_*=2^{-4}$ (right) from a sequence satisfying  and for $\gamma=3/2$.[]{data-label="fig: meshes"}](mesh3-crop.eps "fig:") ![Meshes with $h_*=2^{-3}$ (left) and $h_*=2^{-4}$ (right) from a sequence satisfying  and for $\gamma=3/2$.[]{data-label="fig: meshes"}](mesh4-crop.eps "fig:")

For the time integration, we use a technique [@McLeanThomee2010; @WeidemanTrefethen2007] based on a quadrature approximation to the Laplace inversion formula, $$u_h(t)=\frac{1}{2\pi i}\int_\Gamma
    e^{zt}\hat u_h(z)\,dz
    =\frac{1}{2\pi i}\int_{-\infty}^\infty 
        e^{z(\xi)t}\hat u\bigl(z(\xi)\bigr)z'(\xi)
        \,d\xi,$$ where the contour $\Gamma$ has the parametric representation $$z(\xi)=\mu\bigl(1-\sin(\delta-i\xi)\bigr)
    \quad\text{for $-\infty<\xi<\infty$,}$$ with $\delta=1.1721\,0423$ and $\mu=4.4920\,7528\,M/t$ for given $t>0$ and a chosen positive integer $M$. Therefore, the contour $\Gamma$ is the left branch of an hyperbola with asymptotes $y=\pm(x-\mu)\cot\delta$ for $z=x+iy$. Putting $$z_j=z(\xi_j),\quad z'_j=z'(\xi_j),\quad
\xi_j=j\,{\Delta\xi},\quad{\Delta\xi}=\frac{1.0817\,9214}{M},$$ we define $$U_{M,h}(t)=\frac{{\Delta\xi}}{2\pi i}\sum_{j=-M}^Me^{z_jt}
    \hat u_h(z_j)z'_j\approx u_h(t).$$ To compute $\hat u_h(z_j)$ we solve the (complex) finite element equations $$z_j^\alpha{\langle\hat u_h(z_j),\chi\rangle}+a\bigl(\hat u_h(z_j),\chi\bigr)
    =z_j^{\alpha-1}{\langleu_{0h}+\hat f(z_j),\chi\rangle},\quad\chi\in V_h,$$ and since we choose real $u_{0h}$ and $f$, it follows that $\hat u_h(z_{-j})=\hat u_h(\bar z_j)=\overline{\hat u_h(z_j)}$ so the number of elliptic solves needed to evaluate $U_{M,h}(t)$ is $M+1$, not $2M+1$. An error bound for the quadrature error $\|U_{M,h}(t)-u_h(t)\|$ includes a decay factor $10.1315^{-M}$, and we observe in practice that the overall error $\|U_{M,h}(t)-u(t)\|$ is dominated by the finite element error $\|u_h(t)-u(t)\|$ for modest values of $M$. In the computations reported below we use $M=8$ to compute $U_{M,h}(t)\approx u_h(t)$, and choose $u_{0h}=P_hu_0$ for the discrete initial data.

![Behaviour of the $L_2$-error $\|u_h(t)-u(t)\|$ for Example 1 when $t=1$ — quasiuniform versus locally refined triangulations.[]{data-label="fig: locally-refined"}](Figure8_2.eps)

Example 1 {#example-1 .unnumbered}
---------

In our first example, we use $\alpha=1/2$ and choose $u_0$ and $f$ so that the solution of the initial-boundary value problem for  is $$u(x,y,t)=\bigl(1+\omega_{\alpha+1}(t)\bigr)
    r^\beta(1-r)\sin(\beta\theta).$$ In view of  and , the singular behaviour of $u$ as $r\to0$ or $t\to0$ is typical for such problems. Figure \[fig: locally-refined\] compares the behaviour of the $L_2$-error at $t=1$ for quasi-uniform ($\gamma=1$) and locally-refined ($\gamma=1/\beta=3/2$) triangulations. From Theorems \[thm: non-smooth u0\] and \[thm: u0=0\], we expect errors of order $\epsilon(h,1)^2=h^{2\gamma\beta}=h^{4/3}$ and $\epsilon(h,3/2)^2=h^2\log^2(1+h^{-1})$, respectively. The number of degrees of freedom is of order $h^{-2}$ in both cases, so in Figure \[fig: locally-refined\] we expect the corresponding error curves to be straight lines with gradients $-2/3$ and $-1$, which are in fact close to the observed values $-0.7249$ and $-0.9707$, respectively, as determined by simple linear least squares fits.

![The $L_2$-error as a function of $t$ for Example 2 with $\alpha=1/2$ and $\gamma=2/\beta$.[]{data-label="fig: t dependence"}](Figure8_3.eps)

  ---------- -------- ----------- ------- ----------- ------- ----------- -------
    $h_*$                error     rate      error     rate      error     rate
   $2^{-4}$      1957  1.465e-03           1.485e-03           1.452e-03  
   $2^{-5}$      7593  3.673e-04   1.996   3.723e-04   1.996   3.640e-04   1.997
   $2^{-6}$     29771  9.471e-05   1.955   9.597e-05   1.956   9.380e-05   1.956
   $2^{-7}$    117039  2.420e-05   1.969   2.451e-05   1.970   2.391e-05   1.972
   $2^{-8}$    466089  6.059e-06   1.998   6.119e-06   2.002   5.931e-06   2.011
  ---------- -------- ----------- ------- ----------- ------- ----------- -------

  : $L_2$-Errors and empirical convergence rates (powers of $h$) for Example 2 when $t=1$, with $\gamma=2/\beta$ and different choices of $\alpha$. (Recall that $N$ denotes the number of degrees of freedom in the finite element triangulation.)[]{data-label="tab: mixed"}

Example 2 {#example-2 .unnumbered}
---------

In our second example, we impose mixed boundary conditions: homogeneous Dirichlet for $\theta=0$ or $r=1$, and homogeneous Neumann for $\theta=\pi/\beta$. As the initial data we choose the first eigenfunction of $A=-K\nabla^2$, $$u_0(x,y)=J_{\beta/2}(\omega r)\sin(\tfrac12\beta\theta),$$ where $\omega$ is the first positive zero of the Bessel function $J_{\beta/2}$. We put $f=0$ so (recalling that our choice of $K$ means that the corresponding eigenvalue equals 1) the solution is $u(x,y,t)={\mathcal{E}}(t)u_0=E_\alpha(-t^\alpha)u_0(x,y)$, and choose $\gamma=2/\beta=3$ so that ${\epsilon_{\text{mix}}}(h,\gamma)$ is of order $h\log(1+h^{-1})$; see . Since $A^ru_0\propto u_0\in L_2(\Omega)$ for all $r>0$, we conclude from Theorem \[thm: u0 smooth\] that the $L_2$-error $\|u_h(t)-u(t)\|$ is of order $h^2\log^2(1+h^{-1})$ uniformly for $0\le t\le T$. Figure \[fig: t dependence\] confirms this behaviour in the case $\alpha=1/2$. Finally, Table \[tab: mixed\] shows that at a fixed positive time $t=1$ the $L_2$-error does not vary much with $\alpha$.

[^1]: This work was supported by the Australian Research Council grant DP140101193.
