---
abstract: |
    We consider the initial/boundary value problem for a diffusion equation involving multiple time-fractional derivatives on a bounded convex polyhedral domain. We analyze a space semidiscrete scheme based on the standard Galerkin finite element method using continuous piecewise linear functions. Nearly optimal error estimates for both cases of initial data and inhomogeneous term are derived, which cover both smooth and nonsmooth data. Further we develop a fully discrete scheme based on a finite difference discretization of the time-fractional derivatives, and discuss its stability and error estimate. Extensive numerical experiments for one and two-dimension problems confirm the convergence rates of the theoretical results.\
    [**Keywords**]{}: multi-term time-fractional diffusion equation, finite element method, error estimate, semidiscrete scheme, Caputo derivative
address:
- 'Department of Mathematics, University of California, Riverside, 900 University Ave., Riverside, CA 92521, USA (bangti.jin@gmail.com)'
- |
    Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, USA\
    (lazarov@math.tamu.edu)
- 'Graduate School of Mathematical Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153, Japan (ykliu@ms.u-tokyo.ac.jp)'
- |
    Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, USA\
    (zzhou@math.tamu.edu)
author:
- Bangti Jin
- Raytcho Lazarov
- Yikan Liu
- Zhi Zhou
bibliography:
- 'frac.bib'
date: 'started August, 2013; today is '
title: 'The Galerkin Finite Element Method for A Multi-term Time-Fractional Diffusion equation'
---

\#1[H\^[\#1]{}()]{} \#1\#2[\#1 \_[H\^[\#2]{}()]{}]{}

introduction {#sec:intro}
============

We consider the following initial/boundary value problem for a multi-term time fractional diffusion equation in $u(x,t)$: $$\begin{aligned}
{3}\label{eqn:goveq}
   \PD u-\Delta u&= f,&&\quad \text{in  } \Omega&&\quad T \ge t > 0,\notag\\
   u&=0,&&\quad\text{on}\  \partial\Omega&&\quad T \ge t > 0,\\
    u(0)&=v,&&\quad\text{in  }\Omega,&&\notag\end{aligned}$$ where $\Omega$ denotes a bounded convex polygonal domain in $\mathbb R^d\,(d=1,2,3)$ with a boundary $\partial\Omega$, $f$ is the source term, and the initial data $v$ is a given function on $\Omega$ and $T>0$ is a fixed value. Here the differential operator $\PD$ is defined by $$\PD=\partial_t^{{\alpha}}+\sum_{i=1}^m b_i \partial_t^{{\alpha}_i},$$ where $0< {\alpha}_m  \le ...\le {\alpha}_1<{\alpha}<1$ are the orders of the fractional derivatives, $b_i > 0$, $i=1,2,...,m$, with the left-sided Caputo fractional derivative $\partial_t^\beta u$ being defined by (cf. [@KilbasSrivastavaTrujillo:2006 pp.91]) $$\label{eqn:McT}
   \partial_t^\beta u(t)= \frac{1}{\Gamma(1-\beta)} \int_0^t(t-\tau)^{-\beta}\frac{d}{d\tau}u(\tau) d\tau,$$ where $\Gamma(\cdot)$ denotes the Gamma function.

In the case of $m=0$, the model reduces to its single-term counterpart $$\label{eqn:single}
  \partial_t^\alpha u -\Delta u = f\quad \mbox{ in } \Omega\times(0,T].$$ This model has been studied extensively from different aspects due to its extraordinary capability of modeling anomalous diffusion phenomena in highly heterogeneous aquifers and complex viscoelastic materials [@AdamsGelhar:1992; @Nigmatulin:1986]. It is the fractional analogue of the classical diffusion equation: with $\alpha=1$, it recovers the latter, and thus inherits some of its analytical properties. However, it differs considerably from the latter in the sense that, due to the presence of the nonlocal fractional derivative term, it has limited smoothing property in space and slow asymptotic decay in time [@SakamotoYamamoto:2011], which in turn also impacts related numerical analysis [@JinLazarovZhou:2013] and inverse problems [@JinRundell:2012; @SakamotoYamamoto:2011].

The model was developed to improve the modeling accuracy of the single-term model for describing anomalous diffusion. For example, in [@SchumerBensonMeerschaertBaeumer:2003], a two-term fractional-order diffusion model was proposed for the total concentration in solute transport, in order to distinguish explicitly the mobile and immobile status of the solute using fractional dynamics. The kinetic equation with two fractional derivatives of different orders appears also quite naturally when describing subdiffusive motion in velocity fields [@MetzlerKlafterSokolov:1998]; see also [@KellyMcGoughMeerschaert:2008] for discussions on the model for wave-type phenomena.

There are very few mathematical studies on the model . Luchko [@Luchko:2011] established a maximum principle for problem , and constructed a generalized solution for the case $f\equiv0$ using the multinomial Mittag-Leffler function. Jiang et al [@JiangLiuTurnerBurrage:2012b] derived analytical solutions for the diffusion equation with fractional derivatives in both time and space. Li and Yamamoto [@LiYamamoto:2013] established existence, uniqueness, and the Hölder regularity of the solution using a fixed point argument for problem with variable coefficients $\{b_i\}$. Very recently, Li et al [@LiLiuYamamoto:2013] showed the uniqueness and continuous dependence of the solution on the initial value $v$ and the source term $f$, by exploiting refined properties of the multinomial Mittag-Leffler function.

The applications of the model motivate the design and analysis of numerical schemes that have optimal (with respect to data regularity) convergence rates. Such schemes are especially valuable for problems where the solution has low regularity. The case $m=0$, i.e., the single-term model , has been extensively studied, and stability and error estimates were provided; see [@LinXu:2007; @ZhangSunWu:2011] for the finite difference method, [@LiXu:2009; @ZayernouriKardiadakis:2014] for the spectral method, [@McLeanThomee:2010; @Mustapha:2011; @MustaphaMcLean:2013; @JinLazarovZhou:2013; @JinLazarovPasciakZhou:2013; @JinLazarovPasciakZhou:2013a] for the finite element method, and [@BrunnerLingYamamoto:2010; @FuChenYang:2013] for meshfree methods based on radial basis functions, to name a few. In particular, in [@JinLazarovPasciakZhou:2013a; @JinLazarovPasciakZhou:2013; @JinLazarovZhou:2013], the authors established almost optimal error estimates with respect to the regularity of the initial data $v$ and the right hand side $f$ for a semidiscrete Galerkin scheme. These studies include the interesting case of very weak data, i.e., $ v \in \dH {q}$ and $f\in L^\infty(0,T;\dH q)$ for $-1 < q \le 0$.

Numerical methods for the general multi-term case for an ordinary differential equation were considered in [@Katsikadelis:2009; @ElSayedElKallaZiada:2010]. In [@ZhaoXiaoXu:2013], a scheme based on the finite element method in space and a specialized finite difference method in time was proposed for , and error estimates were derived. We also refer to [@LiuMeerschaert:2013] for a numerical scheme based on a fractional predictor-corrector method for the multi-term time fractional wave-diffusion equation. The error analysis in these works is done under the assumption that the solution is sufficiently smooth and therefore it excludes the case of low regularity solutions. This is the main goal of the present study. However, the derivation of optimal with respect to the regularity error estimates requires additional analysis of the properties of problem , e.g., stability, asymptotic behavior for $t \to 0^+$. Relevant results of this type have recently been obtained in [@LiLiuYamamoto:2013], which, however, are not enough for the analysis of the semidiscrete Galerkin scheme, and hence in Section \[sec:prelim\], we make the necessary extensions.

Now we describe the semidiscrete Galerkin scheme. Let ${\{\mathcal{T}_h\}}_{0<h<1}$ be a family of shape regular and quasi-uniform partitions of the domain $\Omega$ into $d$-simplexes, called finite elements, with a maximum diameter $h$. The approximate solution $u_h$ is sought in the finite element space $X_h$ of continuous piecewise linear functions over the triangulation $\mathcal{T}_h $ $$X_h =\left\{\chi\in H^1_0(\Omega): \ \chi ~~\mbox{is a linear function over}  ~~\tau,   \,\,\,\,\forall \tau \in \mathcal{T}_h\right\}.$$ The semidiscrete Galerkin FEM for problem is: find $ u_h (t)\in X_h$ such that $$\label{eqn:fem}
   {( \PD u_{h},\chi)}+ a(u_h,\chi)= {(f, \chi)}, \quad \forall \chi\in X_h,\ T \ge t >0, \quad u_h(0)=v_h,$$ where $a(u,w)=(\nabla u, \nabla w) ~~ \text{for}\ u, \, w\in H_0^1(\Omega)$, and $v_h \in X_h$ is an approximation of the initial data $v$. The choice of $v_h$ will depend on the smoothness of the initial data $v$. We shall study the convergence of the semidiscrete scheme for the case of initial data $ v \in \dH q$, $-1<q\leq 2$, and right hand side $f\in L^\infty(0,T;\dH q)$, $-1<q<1$. The case of nonsmooth data, i.e., $-1<q\leq 0$, is very common in inverse problems and optimal control [@JinRundell:2012; @SakamotoYamamoto:2011]; see also [@XieZou:2005; @JinLu:2012; @CasasClasonKunisch:2013; @CasasZuazua:2013] for the parabolic counterpart. The goal of this work is to develop a numerical scheme based on the finite element approximation for the model , and provide a complete error analysis. We derive error estimates optimal with respect to the data regularity for the semidiscrete scheme, and a convergence rate $O(h^2+\tau^{2
-\alpha})$ for the fully discrete scheme in case of a smooth solution. Specifically, our essential contributions are as follows. First, we obtain an improved regularity result for the inhomogeneous problem, by allowing less regular source term, cf. Theorem \[thm:regeps2\]. This is achieved by first establishing a new result, i.e., the complete monotonicity of the multinomial Mittag-Leffler function, cf. Lemma \[lem:MMLcm\]. Second, we derive nearly optimal error estimates for a semidiscrete Galerkin scheme for both homogeneous and inhomogeneous problems, cf. Theorems \[thm:SG-smooth\]-\[thm:gal:l2\], which cover both smooth and nonsmooth data. Third, we develop a fully discrete scheme based on a finite difference method in time, and establish its stability and error estimates, cf. Theorem \[thm:estfull\]. We note that the derived error estimate for the fully discrete scheme holds only for smooth solutions.

The rest of the paper is organized as follows. In Section \[sec:prelim\], we recall the solution theory for the model for both homogeneous and inhomogeneous problems, using properties of the multinomial Mittag-Leffler function. The readers not interested in the analysis may proceed directly to Section \[sec:galerkin\]. Almost optimal error estimates for their Galerkin finite element approximations are given in Section \[sec:galerkin\]. Then a fully discrete scheme based on a finite difference approximation of the Caputo fractional derivatives is given in Section \[sec:fulldis\], and an error analysis is also provided. Finally, extensive numerical experiments are presented to illustrate the accuracy and efficiency of the Galerkin scheme, and to verify the convergence theory. Throughout, we denote by $C$ a generic constant, which may differ at different occurrences, but always independent of the mesh size $h$ and time step size $\tau$.

Solution theory {#sec:prelim}
===============

In this part, we recall the solution theory for problem . We shall describe the solution representation using the multinomial Mittag-Leffler function, and derive optimal solution regularity for the homogeneous and inhomogeneous problems.

Multinomial Mittag-leffler function {#ssec:MultiML}
-----------------------------------

First we recall the multinomial Mittag-Leffler function, introduced in [@HadidLuchko:1996]. For $0<\beta<2$, $0<\beta_i<1$ and $z_i\in \mathbb C$, $i=1,...,m$, the multinomial Mittag-Leffler function $E_{(\beta_1,...,\beta_m),\beta}
(z_1,...,z_m)$ is defined by $$    E_{(\beta_1,...,\beta_m),\beta}(z_1,...,z_m)=\sum_{k=0}^{\infty}
    \sum_{\substack{l_1+...+l_m=k\\l_1\ge0,...,l_m\ge0}} (k;l_1,...,l_m)
     \frac{\prod_{i=1}^m z_i^{l_i}}{\Gamma(\beta+\Sigma_{i=1}^m \beta_i l_i)},$$ where the notation $(k;l_1,...,l_m)$ denotes the multinomial coefficient, i.e., $$(k;l_1,...,l_m)=\frac{k!}{l_1!...l_m!}\quad \mbox{with } k = \sum_{i=1}^ml_i.$$

It generalizes the exponential function $e^z$: with $m=1$ and $\beta=\beta_1=1$, it reproduces the exponential function $e^z$. It appears in the solution representation of problem , cf. below. We shall need the following two important lemmas on the function $E_{(\beta_1,...,\beta_m),\beta}(z_1,...,z_m)$, recently obtained in [@LiLiuYamamoto:2013 Section 2.1].

\[lem:MLbound\] Let $0<\beta<2$, $0<\beta_i<1$, $\beta_1>\max\{\beta_2,...,\beta_m\}$ and $\frac{\beta_1\pi}{2}<\mu<\beta_1\pi$. Assume that there is $K >0$ such that $-K\le z_i<0$, $i=2,...,m$. Then there exists a constant $C=C(\beta_1,...,\beta_m,\beta,K,\mu)>0$ such that $$    E_{(\beta_1,...,\beta_m),\beta}(z_1,...,z_m)\le \frac{C}{1+|z_1|},
        \quad \quad\quad \mu\leq|\mathrm{arg}(z_1)|\leq \pi.$$

\[lem:multiMLprop\] Let $0<\beta<2$, $0<\beta_i<1$ and $z_i \in \mathbb C$, $i=1,...,m$. Then we have $$\frac{1}{\Gamma(\beta_0)}+ \sum_{i=1}^{m}z_iE_{(\beta_1,...,\beta_m),{\beta_0+\beta_i}}(z_1,...,z_m)
    =E_{(\beta_1,...,\beta_m),\beta_0}(z_1,...,z_m).$$

Solution Representation {#ssec:represent}
-----------------------

For $s\ge-1$, we denote by $\dH s\subset H^{-1}(\Omega)$ the Hilbert space induced by the norm: $$\|v\|_{\dH s}^2=\sum_{j=1}^{\infty}{\lambda}_j^s \langle v,{\varphi}_j \rangle^2$$ with $\{{\lambda}_j\}_{j=1}^\infty$ and $\{{\varphi}_j\}_{j=1}^\infty$ being respectively the eigenvalues and the $L^2(\Omega)$-orthonormal eigenfunctions of the Laplace operator $-\Delta$ on the domain $\Omega$ with a homogeneous Dirichlet boundary condition. Then $\{{\varphi}_j\}_{j=1}^\infty$ and $\{{\lambda}_j^{1/2}
{\varphi}_j\}_{j=1}^\infty$, form an orthonormal basis in $L^2(\Omega)$ and $H^{-1}(\Omega)$, respectively. Further, $\|v\|_{\dH 0}=\|v\|_{L^2(\Omega)}=(v,v)^{1/2}$ is the norm in $L^2(\Omega)$ and $\|v\|_{\dH {-1}}
= \|v\|_{H^{-1}(\Omega)}$ is the norm in $H^{-1}(\Omega)$. It is easy to verify that $\|v\|_{\dH 1}= \|\nabla v\|_{L^2(\Omega)}$ is also the norm in $H_0^1(\Omega)$ and $\|v\|_{\dH 2}=\|\Delta v\|_{L^2(\Omega)}$ is equivalent to the norm in $H^2(\Omega)\cap H^1_0(\Omega)$ [@Thomee:2006 Lemma 3.1]. Note that $\dH s$, $s\ge -1$ form a Hilbert scale of interpolation spaces. Hence, we denote $\|\cdot\|_{H^s(\Omega)}$ to be the norm on the interpolation scale between $H^1_0(\Omega)$ and $L^2(\Omega)$ for $s\in
[0,1]$ and $\|\cdot\|_{H^{s}(\Omega)}$ to be the norm on the interpolation scale between $L^2(\Omega)$ and $H^{-1}(\Omega)$ for $s\in [-1,0]$. Then, $\| \cdot \|_{H^s(\Omega)}$ and $\|\cdot\|_{\dH s}$ are equivalent for $s\in [-1,1]$. Further, for a Banach space $B$, we define the space $$L^r(0,T;B) = \{u(t)\in B \mbox{ for a.e. } t\in (0,T) \mbox{ and } \|u\|_{L^r(0,T;B)}<\infty\},$$ for any $r\geq 1$, and the norm $\|\cdot\|_{L^r(0,T;B)}$ is defined by $$\|u\|_{L^r(0,T;B)} = \left\{\begin{aligned}\left(\int_0^T\|u(t)\|_B^rdt\right)^{1/r}, &\quad r\in [1,\infty),\\
     \displaystyle  {\mathrm{esssup}_{t\in(0,T)}}\|u(t)\|_B, &\quad r= \infty.
       \end{aligned}\right.$$

Upon denoting $\vecal=({\alpha},{\alpha}-{\alpha}_1,...,{\alpha}-{\alpha}_m)$, we introduce the following solution operator $$\label{eqn:opE}
    E(t)v=\sum_{j=1}^{\infty} \left(1-{\lambda}_j t^{{\alpha}}
    E_{\vecal,1+{\alpha}}(-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m})\right) (v,{\varphi}_j){\varphi}_j.$$ This operator is motivated by a separation of variable [@LuchkoGorenflo:1999; @Luchko:2011]. Then for problem with a homogeneous right hand side, i.e., $f\equiv0$, we have $u(x,t)=E(t)v$. However, the representation is not always very convenient for analyzing its smoothing property. We derive an alternative representation of the solution operator $E$ using Lemma \[lem:multiMLprop\]: $$\label{eqn:opE2}
\begin{split}
    E(t)v            = &\sum_{j=1}^{\infty}
    E_{\vecal,1} (-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m}) (v,{\varphi}_j){\varphi}_j\\
      &\ \ + \sum_{i=1}^m b_it^{{\alpha}-{\alpha}_i}
       \sum_{j=1}^{\infty} E_{\vecal,1+{\alpha}-{\alpha}_i}
       (-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m}) (v,{\varphi}_j){\varphi}_j.
\end{split}$$

Besides, we define the following operator $\bar{E}$ for $\chi\in L^2(\Omega)$ by $$\label{eqn:Ebar}
    \bar E(t)\chi=\sum_{j=1}^{\infty}  t^{{\alpha}-1}
    E_{\vecal,{\alpha}}(-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m}) (\chi,{\varphi}_j){\varphi}_j.$$ The operators $E(t) $ and $ {\bar E}(t) $ can be used to represent the solution $u$ of as: $$\label{eqn:solrep}
   u(t)=E(t)v + \int_0^t  {\bar E}(t-s) f(s) ds.$$

The operator $\bar{E}$ has the following smoothing property.

\[lem:barE\] For any $t>0$ and $\chi\in \dH q$, $q\in(-1,2]$, there holds for $0\le p-q \le 2$ $$\|\bar E(t) \chi \|_{\dH p} \le Ct^{-1+{\alpha}(1+(q-p)/2)}\|\chi\|_{\dH q}.$$

The definition of the operator $\bar{E}$ in and Lemma \[lem:MLbound\] yield $$\begin{split}
\|\bar E(t) \chi \|_{\dH p}^2
      &=t^{-2+(2+q-p){\alpha}}\sum_{j=1}^{\infty} ({\lambda}_j t^{{\alpha}})^{p-q}
      | E_{\vecal,{\alpha}}(-{\lambda}_j t^{\alpha},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m})|^2 {\lambda}_j^q |(\chi,{\varphi}_j)|^2\\
      &\le C t^{-2+(2+q-p){\alpha}}\sum_{j=1}^{\infty} \frac{({\lambda}_j t^{{\alpha}})^{p-q}}{(1+{\lambda}_jt^{{\alpha}})^2} {\lambda}_j^q |(\chi,{\varphi}_j)|^2\\
      &\le C t^{-2+(2+q-p){\alpha}}\sum_{j=1}^{\infty} {\lambda}_j^q |(\chi,{\varphi}_j)|^2 \le C t^{-2+(2+q-p){\alpha}}\|\chi\|_{\dH q},
\end{split}$$ where the last line follows by the inequality $\sup_{j\in \mathbb{N}}
\frac{ ({\lambda}_j t^{\alpha})^{p-q}}{(1+{\lambda}_j t^{\alpha})^2}\le C$, for $0\leq p-q\leq 2$.

Solution regularity
-------------------

First we recall known regularity results. In [@LiYamamoto:2013], Li and Yamamoto showed that in the case of variable coefficients $\{b_i(x)\}$, there exists a unique mild solution $u\in C((0,T];
\dH{\gamma})\cap C([0,T];L^2(\Omega))$ and $u\in C([0,T];\dH {\gamma}) \cap L^{\infty}(0,T;\dH 2)$ when $ v \in L^2(\Omega)$, $f=0$ and $v =0$, $f\in L^{\infty}(0,T]; L^2(\Omega))$, respectively, with $\gamma\in[0,1)$. These results were recently refined in [@LiLiuYamamoto:2013] for the case of constant coefficients, i.e., problem . In particular, it was shown that for $v\in \dH q$, $0\leq q
\leq 1$, and $f=0$, $u\in L^{1/(1-q/2)}(0,T;H^2(\Omega)\cap H_0^1(\Omega))$; and for $v=0$ and $f\in L^r(0,T;
\dH q)$, $0\leq q\leq 2$, $r\geq1$, $u\in L^r(0,T;\dH {q+2-\gamma})$ for some $\gamma\in(0,1]$. Here we follow the approach in [@LiLiuYamamoto:2013], and extend these results to a slightly more general setting of $v\in \dH q$, $-1<q\leq 2$, and $f\in L^2(0,T;\dH q)$, $-1<q\le 1$. The nonsmooth case, i.e., $-1<q\leq 0$, arises commonly in related inverse problems and optimal control problems.

We shall derive the solution regularity to the homogeneous problem, i.e., $f\equiv 0$, and the inhomogeneous problem, i.e., $v\equiv 0$, separately. These results will be essential for the error analysis of the space semidiscrete Galerkin scheme in Section \[sec:galerkin\]. First we consider the homogeneous problem with initial data $v\in \dH q$, $-1<q\leq 2$.

\[thm:homogreg\] Let $u(t)=E(t)v$ be the solution to problem with $f\equiv 0$ and $v\in \dH q$, $q\in(-1,2]$. Then there holds $$ \|\PD^\ell u(t) \|_{\dH p}
 \le Ct^{-{\alpha}(\ell+(p-q)/2)}\|v\|_{\dH q},\quad t>0,$$ where for $\ell=0$, $ 0 \le p-q\le 2$ and for $\ell=1$, $-2 \le p-q \le 0$.

We show that represents indeed the weak solution to problem with $f\equiv 0$ and further it satisfies the desired estimate. We first discuss the case $\ell=0$. By Lemma \[lem:MLbound\] and we have for $0\le p-q \le 2$ $$\begin{split}
    \| E(t)v \|_{\dH p}^2 & = \displaystyle
       \sum_{j=1}^{\infty}\lambda_j^p \Big (E_{\vecal,1} (-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m}) \Big .\\
     & \Big . \quad\quad +\sum_{i=1}^mb_it^{{\alpha}-{\alpha}_i}
       E_{\vecal,1+{\alpha}-{\alpha}_i}
       (-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m}) \Big )^2 (v,{\varphi}_j)^2\\
              &\le C t^{-(p-q){\alpha}}\sum_{j=1}^{\infty}\frac{({\lambda}_jt^{{\alpha}})^{p-q}}{(1+{\lambda}_jt^{\alpha})^2}
         {\lambda}_j^q |(v,{\varphi}_j)|^2\le C t^{-(p-q){\alpha}}  \| v \|_{\dH q}^2,
\end{split}$$ where the last line follows from the inequality $\sup_{j\in\mathbb N}\frac{ ({\lambda}_jt^{\alpha})^{p-q}
}{(1+{\lambda}_jt^{\alpha})^2}\le C$ for $0\le p-q\le 2$. The estimate for the case $\ell=1$ follows from the identity $\|\PD E(t)v \|_{\dH p} = \| E(t)v \|_{\dH {p+2}}$. It remains to show that satisfies also the initial condition in the sense that $\lim_{t\rightarrow 0^+}
\| E(t)v-v  \|_{\dH q}=0$. By identity and Lemma \[lem:MLbound\], we deduce $$\begin{split}
    \| E(t)v-v  \|_{\dH q}^2
       &=\sum_{j=1}^{\infty} {\lambda}_j^{2} t^{2{\alpha}}
       \bigg|E_{\vecal,1+{\alpha}}(-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m})\bigg|^2{\lambda}_j^q| (v,{\varphi}_j)|^2\\
       & \leq C\|v\|_{\dH q}^2 <\infty.
\end{split}$$ Using Lemma \[lem:multiMLprop\], we rewrite the term ${\lambda}_j t^{{\alpha}}E_{\vecal,1+{\alpha}}
(-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m})$ as $$\begin{aligned}
      {\lambda}_j t^{{\alpha}} E_{\vecal,1+{\alpha}}&(-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m})\\
    = & (1-E_{\vecal,1} (-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m})) \\
    & \ \ -  \sum_{i=1}^m b_it^{{\alpha}-{\alpha}_i}E_{\vecal,1+{\alpha}-{\alpha}_i}
       (-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m}).
  \end{aligned}$$ Upon noting the identity $\lim_{t\to 0^+}(1-E_{\vecal,1} (-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,
-b_mt^{{\alpha}-{\alpha}_m}))=0$, and the boundedness of $E_{\vecal,1+{\alpha}-{\alpha}_i}(-{\lambda}_jt^{{\alpha}},
-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m})$ from Lemma \[lem:MLbound\], we deduce that for all $j$ $$\lim_{t\to0^+}{\lambda}_j t^{{\alpha}} E_{\vecal,1+{\alpha}}(-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m})=0.$$ Hence, the desired assertion follows by Lebesgue’s dominated convergence theorem.

Now we turn to the inhomogeneous problem with a nonsmooth right hand side, i.e., $f\in L^{\infty}(0,T;\dH q)$, $-1<q\leq 1$, and a zero initial data $v\equiv0$.

\[thm:inhomogreg\] For $f\in L^{\infty}(0,T;\dH q)$, $-1<q\le1$, and $v\equiv 0$, the solution $u$ belongs to $L^{\infty}(0,T;\dH {q+2-{\epsilon}})$ for any $\epsilon>0$ and satisfies $$\label{eq:regeps}
    \|u(\cdot,t)\|_{\dH {q+2-\epsilon}} \leq C\epsilon^{-1}t^{{\epsilon}{\alpha}/2}\|f\|_{L^\infty(0,t;\dH q)}.$$ Hence, it is a solution to problem with a homogeneous initial data $v=0$.

By construction, it satisfies the governing equation. By Lemma \[lem:barE\], we have $$\begin{aligned}
    \|u(\cdot,t)\|_{\dH {q+2-{\epsilon}}} & = \|\int_0^t \bar{E}(t-s)f(s)ds\|_{\dH {q+2-{\epsilon}}}\\
      & \leq \int_0^t \|\bar{E}(t-s)f(s)\|_{\dH {q+2-{\epsilon}}} ds \\
      & \leq C\int_0^t (t-s)^{{\epsilon}{\alpha}/2-1} \|f(s)\|_{\dH q}ds\\
      &\leq C\epsilon^{-1}t^{{\epsilon}{\alpha}/2}\|f\|_{L^\infty(0,t;\dH q)}
  \end{aligned}$$ which shows the desired estimate. Further, it satisfies the initial condition $u(0)=0$, i.e., for any ${\epsilon}>0$, $\lim_{t\to 0^+}\|u(\cdot,t)
\|_{\dH {q+2-{\epsilon}}}=0$, and thus it is indeed a solution of .

Next we extend Theorem \[thm:inhomogreg\] to allow less regular right hand sides $f\in L^2(0,T;\dH q)$, $-1<q\le 1$. Then the function $u(x,t)$ satisfies also the differential equation as an element in the space $L^2(0,T;\dot H^{q+2}(\Omega))$. However, it may not satisfy the homogeneous initial condition $u(x,0)=0$. In Remark \[weakest\_sol\] below, we argue that the weakest class of source term that produces a legitimate weak solution of is $f \in
L^r(0,T;\dH q)$ with $r>1/{\alpha}$ and $-1 < q \le 1$. Obviously, for $1/2<{\alpha}<1$, it does give a solution $u(x,t) \in L^2(0,T;\dH {q+2})$. To this end, we introduce the shorthand notation $$\bar E_{\vecal}^j(t)  = t^{{\alpha}-1} E_{\vecal,{\alpha}}(-{\lambda}_jt^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m}).$$

The function $\bar E_{\vecal}^j(t)$ is completely monotone; see Appendix \[app:MMLcm\] for the technical proof.

\[lem:MMLcm\] The function $ \bar E_{\vecal}^j(t)$ for $j\in \mathbb{N}$ has the following properties: $$\begin{aligned}
     \bar E_{\vecal}^j(t)\ \mbox { is completely monotone} \quad \mbox{and}\quad
    \int_0^T | \bar E_{\vecal}^j(t)|\,dt<\frac{1}{{\lambda}_j}. $$

\[thm:regeps2\] For $f\in L^2(0,T;\dH q)$, $-1<q\le1$, the representation belongs to $L^{2}(0,T;\dH {q+2})$ and satisfies the [*a priori*]{} estimate $$\label{eq:regeps2}
    \|u\|_{L^2(0,t;\dH {q+2})} + \| \PD u \|_{L^2(0,t;\dH {q})} \leq C\|f\|_{L^2(0,t;\dH q)}.$$

By Young’s inequality for the convolution $\|k\ast f\|_{L^p}\le\|k\|_{L_1}\|f\|_{L^p}$, $k\in L^1$, $f\in L^p$, $p\ge 1$, and Lemma \[lem:MMLcm\], we deduce $$\begin{aligned}
      \left\|\int_0^t \bar E_{\vecal}^n (t-\tau)f_n(\tau)\,d\tau\right\|_{L^2(0,T)}^2&\le
      \left(\int_0^T|\bar E_{\vecal}^n (t)|\,dt\right)^2  \left(\int_0^T |f_n(t)|^2\,dt\right)
      \le \frac{1}{{\lambda}_n^2} \int_0^T |f_n(t)|^2\,dt.
   \end{aligned}$$ Hence, $$\begin{aligned}
    \| u\|_{L^2(0,T;\dot H^{q+2}(\Omega))}^2&\le \sum_{n=1}^{\infty}{\lambda}_n^{q+2}\left\|\int_0^t
    \bar E_{\vecal}^n (t-\tau)f_n(\tau)\,d\tau\right\|_{L^2(0,T)}^2\\
    &\le \sum_{n=1}^{\infty}\lambda_n^q\int_0^T |f_n(t)|^2\,dt= \|f\|_{L^2(0,T;\dH q)}^2.
  \end{aligned}$$ The estimate on $\| \PD u \|_{L^2(0,t;\dH {q})}$ follows analogously. This completes the proof.

\[weakest\_sol\] The condition $f\in  L^\infty(0,T;\dH q)$ in Theorem \[thm:inhomogreg\] can be weakened to $f\in L^r(0,T;\dH q)$ with $r>1/\alpha$. This follows from Lemma \[lem:barE\] and Hölder’s inequality with $r'$, $1/r' + 1/r=1$ $$\begin{aligned}
    \|u(\cdot,t)\|_{\dH q} & \leq \int_0^t \|\bar{E}(t-s)f(s)\|_{\dH {q}} ds
       \leq C\int_0^t (t-s)^{\alpha-1} \|f(s)\|_{\dH q}ds \\
     &  \leq C\left(\frac{t^{1+r'(\alpha-1)}}{1+r'(\alpha-1)}\right)^{1/r'}\|f\|_{L^r(0,t;\dot H^q(\Omega))},
  \end{aligned}$$ where $1+r'(\alpha-1)>0$ by the condition $r>1/\alpha$. It follows from this that the initial condition $u(\cdot,0)=0$ holds in the following sense: $\lim_{t\to0^+}\|u(\cdot,t)\|_{\dH q}=0$. Hence for any $\alpha\in (1/2,1)$ the representation remains a legitimate solution under the weaker condition $f\in L^2(0,T;\dH q)$.

Error Estimates for Semidiscrete Galerkin Scheme {#sec:galerkin}
================================================

Now we derive and analyze a space semidiscrete Galerkin finite element scheme. First we describe the semidiscrete scheme, and then derive almost optimal error estimates for the homogeneous and inhomogeneous problems separately. In the analysis we essentially use the technique developed in [@JinLazarovZhou:2013] and improved in [@JinLazarovPasciakZhou:2013; @JinLazarovPasciakZhou:2013a].

Semidiscrete scheme
-------------------

To describe the scheme, we need the $L^2(\Omega)$ projection $P_h:L^2(\Omega)\to X_h$ and Ritz projection $R_h:H^1_0(\Omega)\to X_h$, respectively, defined by $$\begin{aligned}
    (P_h \psi,\chi) & =(\psi,\chi) \quad\forall \chi\in X_h,\\
    (\nabla R_h \psi,\nabla\chi) & =(\nabla \psi,\nabla\chi) \quad \forall \chi\in X_h.
  \end{aligned}$$

The operators $R_h$ and $P_h$ satisfy the following approximation property.

\[lem:prh-bound\] For any $\psi\in \dH q$, $q=1,2$, the operator $R_h$ satisfies: $$\begin{aligned}
   \|R_h \psi-\psi\|_{L^2(\Omega)}+h\|\nabla(R_h \psi-\psi)\|_{L^2(\Omega)}\le Ch^q\| \psi\|_{\dot H^q(\Omega)}.\end{aligned}$$ Further, for $s\in [0,1]$ we have $$\begin{aligned}
     \|(I-P_h)\psi \|_{H^s(\Omega)} &\le Ch^{2-s} \|\psi\|_{\dH 2}\quad \forall \psi\in
       H^2(\Omega)\cap H^1_0(\Omega),\\
     \|(I-P_h)\psi \|_{H^s(\Omega)} &\le Ch^{1-s} \|\psi\|_{\dH 1}\quad \forall \psi\in H^1_0(\Omega).
  \end{aligned}$$

Now we can describe the semidiscrete Galerkin scheme. Upon introducing the discrete Laplacian $\Delta_h: X_h\to X_h$ defined by $$-(\Delta_h\psi,\chi)=(\nabla\psi,\nabla\chi)\quad\forall\psi,\,\chi\in X_h,$$ and $f_h= P_h f$, we may write the spatially discrete problem as $$\label{eqn:fem-operator}
   \PD u_{h}(t)-\Delta_h u_h(t) =f_h(t) \for t\ge0 \quad \mbox{with} \quad  u_h(0)=v_h,$$ where $v_h\in X_h$ is an approximation to the initial data $v$. Next we give a solution representation of using the eigenvalues and eigenfunctions $\{{\lambda}^h_j\}_{j=1}^{N}$ and $\{{\varphi}_j^h\}_{j=1}^{N}$ of the discrete Laplacian $-\Delta_h$. First we introduce the operators $E_h$ and $\bar{E}_h$, the discrete analogues of and , for $t>0 $, defined respectively by $$\label{eqn:Eh}
  \begin{split}
    E_h(t)v_h&=\sum_{j=1}^{N}E_{\vecal,1} (-{\lambda}_j^ht^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m}) (v,{\varphi}_j^h){\varphi}_j^h\\
      & + \sum_{i=1}^m b_it^{{\alpha}-{\alpha}_i} \sum_{j=1}^{N} E_{\vecal,1+{\alpha}-{\alpha}_i}
       (-{\lambda}_j^ht^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m}) (v,{\varphi}_j^h){\varphi}_j^h,\\
\end{split}$$ and $$\label{eqn:Ehbar}
 \Etilh(t) f_h = \sum_{j=1}^N t^{{\alpha}-1} E_{\vecal,{\alpha}}
 (-{\lambda}_j^ht^{{\alpha}},-b_1t^{{\alpha}-{\alpha}_1},...,-b_mt^{{\alpha}-{\alpha}_m})\,(f_h,{\varphi}^h_j) \, {\varphi}_j^h.$$ Then the solution $u_h$ of the discrete problem can be expressed by: $$\label{eqn:Duhamelh}
     u_h(x,t)= E_h(t) v_h + \int_0^t \Etilh(t-s) f_h(s)\,ds.$$

On the finite element space $X_h$, we introduce the discrete norm $\tribar \cdot\tribar_{\dH p}$ defined by $$\tribar \psi\tribar_{\dH p}^2 =
      \sum_{j=1}^N({\lambda}_j^h)^p(\psi,{\varphi}_j^h)^2\quad \psi\in X_h.$$ The norm $\tribar \cdot\tribar_{\dH p}$ is well defined for all real $p$. Clearly, $\tribar \psi\tribar_{\dH 1}=\|\psi\|_{\dH 1}$ and $\tribar \psi\tribar_{\dH 0}=\|\psi\|_{L^2(\Omega)}$ for any $\psi\in X_h$. Further, the following inverse inequality holds [@JinLazarovZhou:2013]: if the mesh $\mathcal{T}_h$ is quasi-uniform, then for any $l>s$ $$\label{eqn:inverse}
  \tribar \psi\tribar_{\dH l}\le Ch^{s-l}\tribar \psi\tribar_{\dH {s}}\quad \forall \psi\in X_h.$$

\[lem:regh\] Assume that the mesh $\mathcal{T}_h$ is quasi-uniform. Then for any $v_h \in X_h$ the function $u_h(t)=E_h(t)v_h$ satisfies $$\tribar \PD^\ell u_h(t) \tribar_{\dH p}
        \le Ct^{ -{\alpha}(\ell + (p-q)/2)}\tribar v_h\tribar_{\dH q}, \quad t >0,$$ where for $\ell=0$, $0\leq p-q\leq 2$ and for $\ell=1$, $p \le q \le p+2$.

Upon noting $\tribar \PD E_h(t)v_h \tribar_{\dH p}=\tribar  E_h(t)v_h \tribar_{\dH {p+2}}$, it suffices to show the case $\ell=0$. Using the representation and Lemma \[lem:MLbound\], we have for $0\leq p-q\leq 2$ $$\begin{split}
    \tribar E_h(t)v_h \tribar_{\dH p}^2 &\le C \sum_{j=1}^{N}\frac{({\lambda}_j^h)^p}{(1+{\lambda}_j^h t^{\alpha})^2} |(v_h,{\varphi}_j^h)|^2\\
       &\le C \sup_{1 \le j\le N}\frac{ ({\lambda}_j^h t^{\alpha})^{p-q}}{(1+{\lambda}_j^h t^{\alpha})^2}
       t^{-(p-q){\alpha}}\sum_{j=1}^N\frac{({\lambda}_j^h t^{{\alpha}})^{p-q}}{(1+{\lambda}_j^h t^{\alpha})^2}({\lambda}_j^h)^q |(v_h,{\varphi}_j^h)|^2\\
       &\le C t^{-(p-q){\alpha}}  \tribar v_h \tribar_{\dH q}^2,
\end{split}$$ where the last inequality follows from $\sup_{1 \le j\le N}\frac{
({\lambda}_j^h t^{\alpha})^{p-q}}{(1+{\lambda}_j^h t^{\alpha})^2}\le C$ for $0\leq p-q\leq 2$.

The next result is a discrete analogue to Lemma \[lem:barE\].

\[lem:reg-f\] Let $\Etilh$ be defined by and $ \chi \in X_h$. Then for all $t >0$ $$\tribar \Etilh(t) \chi \tribar_{\dH p} \le \left \{
\begin{array}{ll}
  Ct^{ -1 + {\alpha}(1 + (q -p)/2)}\tribar \chi \tribar_{\dH q}, & \quad 0\le p-q\le 2, \\[1.3ex]
  Ct^{ -1 + \alpha }\tribar \chi \tribar_{\dH q},  & \quad p< q.
\end{array} \right.$$

The proof for the case $0\le p-q\le2$ is similar to Lemma \[lem:barE\]. The other assertion follows from the fact that $\{\lambda_j^h\}_{j=1}^N$ are bounded from zero independent of $h$.

Error estimates for the homogeneous problem
-------------------------------------------

To derive error estimates, first we consider the case of smooth initial data, i.e., $v \in \dH 2$. To this end, we split the error $u_h(t)-u(t)$ into two terms: $$u_h-u= (u_h-R_hu)+(R_hu-u):={\vartheta}+ \varrho.$$ By Lemma \[lem:prh-bound\] and Theorem \[thm:homogreg\], we have for any $t>0$ $$\label{eqn:rho-bound}
 \| \varrho(t) \|_{L^2(\Omega)} + h \|\nabla\varrho(t)\|_{L^2(\Omega)}\le Ch^2 t^{-(1-q/2){\alpha}}\|v\|_{\dH q} \quad v\in \dH q, q=1,2.$$ So it suffices to get proper estimates for ${\vartheta}(t) $, which is given below.

\[lem:vth-smooth\] The function ${\vartheta}(t):=u_h(t)-R_hu(t)$ satisfies for $p=0,1$ $$\|{\vartheta}(t)\|_{\dH p} \leq Ch^{2-p}\|v\|_{\dH 2} .$$

Using the identity $\Delta_hR_h=P_h\Delta$, we note that ${\vartheta}$ satisfies $$\PD {\vartheta}(t) -\Delta_h{\vartheta}(t) =-P_h \PD \varrho(t) \for t > 0,$$ with ${\vartheta}(0)=0$. By the representation , $${\vartheta}(t)=-\int_0^t \Etilh(t-s) P_h \PD \varrho(s)\,ds.$$ Then by Lemmas \[lem:reg-f\] and \[lem:prh-bound\], and Theorem \[thm:homogreg\] we have for $p=0,1$ $$\begin{split}
  \|{\vartheta}(t)\|_{\dH p} &\le \int_0^t \| \Etilh(t-s) P_h \PD \varrho(s) \|_{\dH p} \,ds  \\
    & \le C \int_0^t (t-s)^{(1-p/2){\alpha}-1} \|\PD \varrho(s) \|_{L^2(\Omega)}\,ds \\
    & \le C h^{2-p}\int_0^t (t-s)^{(1-p/2){\alpha}-1} \|\PD u(s) \|_{\dH {2-p}}\,ds \\
    & \le C h^{2-p} \int_0^t (t-s)^{(1-p/2){\alpha}-1} s^{-(1-p/2){\alpha}} \,ds \|v\|_{\dH 2} \le Ch^{2-p}\|v\|_{\dH 2},
\end{split}$$ which is the desired result.

Using , Lemma \[lem:vth-smooth\] and the triangle inequality, we arrive at our first estimate, which is formulated in the following Theorem:

\[thm:SG-smooth\] Let $v\in \dH 2$ and $f\equiv 0$, and $u$ and $u_h$ be the solutions of and with $v_h=R_hv$, respectively. Then $$\|u_h(t)- u(t)\|_{L^2(\Omega)} + h\|\nabla (u_h(t) - u(t))\|_{L^2(\Omega)} \le Ch^2 \|v\|_{\dH 2}.$$

Now we turn to the nonsmooth case, i.e., $v\in \dH q$ with $-1< q \leq  1$. Since the Ritz projection $R_h$ is not well-defined for nonsmooth data, we use instead the $L^2
(\Omega)$-projection $v_h=P_hv$ and split the error $u_h-u$ into: $$\label{eqn:splitnonsmooth}
  u_h-u=(u_h-P_hu)+(P_hu-u):={{\widetilde \vartheta}}+ {{\widetilde \varrho}}.$$ By Lemma \[lem:prh-bound\] and Theorem \[thm:homogreg\] we have for $-1\leq q\leq 1$ $$  \| {{\widetilde \varrho}}(t) \|_{L^2(\Omega)} + h \|\nabla{{\widetilde \varrho}}(t) \|_{L^2(\Omega)}
  \leq Ch^2\|u(t)\|_{\dH 2}
   \leq Ch^2 t^{-{\alpha}(1-q/2)}\|v\|_{\dH q}. $$ Thus, we only need to estimate the term ${{\widetilde \vartheta}}(t)$, which is stated in the following lemma.

\[lem:vtht\] Let ${{\widetilde \vartheta}}(t)=u_h(t)-P_hu(t)$. Then for $p=0,1$, $-1< q\le1$, there holds (with $\ell_h=|\ln h|$) $$\|{{\widetilde \vartheta}}(t)\|_{\dH p}\leq Ch^{\min(q,0)+2-p}\ell_h
    t^{-{\alpha}\left(1-\max(q/2,0)\right)} \| v \|_{\dH q}.$$

Obviously, $ P_h \PD {{\widetilde \varrho}}= \PD P_h(P_hu-u)=0$ and using the identity $\Delta_hR_h=P_h\Delta$, we get the following problem for ${{\widetilde \vartheta}}$: $$\label{eqn:vtht2}
 \PD {{\widetilde \vartheta}}(t) -\Delta_h {{\widetilde \vartheta}}(t) =
- \Delta_h (R_h u - P_h u)(t), \quad t>0, \quad {{\widetilde \vartheta}}(0)=0.$$ Using , ${{\widetilde \vartheta}}(t)$ can be represented by $$\label{eqn:vtht}
  {{\widetilde \vartheta}}(t) = - \int_0^t\Etilh(t-s)\Delta_h(R_hu-P_hu)(s)\,ds.$$ Let $A=\Etilh(t-s)\Delta_h(R_hu-P_hu)(s)$. Then by Lemma \[lem:regh\], there holds for $p=0,1$: $$\begin{split}
 \| A \|_{\dH p}
   &\le C (t-s)^{{\epsilon}{\alpha}/2-1}\tribar \Delta_h(R_hu-P_hu)(s) \tribar_{\dH {p-2+{\epsilon}}} \\
   &\le C (t-s)^{{\epsilon}{\alpha}/2-1}\tribar(R_hu-P_hu)(s) \tribar_{\dH {p+{\epsilon}}}.\\
\end{split}$$ Then by , Theorem \[thm:homogreg\], Lemma \[lem:prh-bound\] we have for $p=0,1$ and $-1\le q\le 1$ $$\begin{split}
 \| A \|_{\dH p}
    &\le Ch^{\min(q,0)+2-p-{\epsilon}}(t-s)^{{\epsilon}{\alpha}/2-1}\| u(s) \|_{\dH {\min(q,0)+2}}\\
    &\le Ch^{\min(q,0)+2-p-{\epsilon}}(t-s)^{{\epsilon}{\alpha}/2-1}
    s^{-\left(1-\max(q/2,0)\right){\alpha}}\| v \|_{\dH {q}}.
\end{split}$$ Then plugging the estimate into yields $$\begin{split}
 \| {{\widetilde \vartheta}}\|_{\dH p}
    &\le Ch^{\min(q,0)+2-p-{\epsilon}} \int_0^t (t-s)^{{\epsilon}{\alpha}/2-1}
    s^{-\left(1-\max(q/2,0)\right){\alpha}} \, ds \| v \|_{\dH {q}}\\
    &\le C\epsilon^{-1}h^{\min(q,0)+2-p-{\epsilon}} t^{-{\alpha}\left(1-\max(q/2,0)\right)}\| v \|_{\dH q}.
\end{split}$$ Now with the choice ${\epsilon}=1/\ell_h$, we obtain the desired estimate.

Now the triangle inequality yields an error estimate for nonsmooth initial data.

\[thm:SG-nonsmooth\] Let $f\equiv 0$, $u$ and $u_h$ be the solutions of with $v\in \dH q$, $-1<q\le 1$, and with $v_h=P_hv$, respectively. Then with $\ell_h=|\ln h|$, there holds $$\| u_h(t) - u(t) \|_{L^2(\Omega)} +   h  \|\nabla(u_h(t) - u(t))\|_{L^2(\Omega)}
  \le Ch^{\min(q,0)+2} \, \ell_h \,t^{-{\alpha}(1-\max(q/2,0))}\|v\|_{\dH q}.$$

Error estimates for the inhomogeneous problem
---------------------------------------------

Now we derive error estimates for the semidiscrete Galerkin approximations of the inhomogeneous problem with $f\in L^{\infty}(0,T;\dH q)$, $-1<q\leq 0$, and $v\equiv0$, in both $L^2$ and $L^\infty$-norm in time. To this end, we appeal again to the splitting . By Theorem \[thm:inhomogreg\] and Lemma \[lem:prh-bound\], the following estimate holds for ${{\widetilde \varrho}}$: $$\| {{\widetilde \varrho}}(t) \|_{L^2(\Omega)} + h \|\nabla{{\widetilde \varrho}}(t)\|_{L^2(\Omega)} \le Ch^{2+q-\epsilon} \|u(t)\|_{\dH {2+q-\epsilon}}
          \le C\epsilon^{-1}h^{2+q-\epsilon}\|f\|_{L^\infty(0,t;\dH q)}.$$ Now the choice $\ell_h= |\ln h|,\, \epsilon=1/\ell_h$, yields $$\label{eqn:Ph-bound}
 \| {{\widetilde \varrho}}(t) \|_{L^2(\Omega)} + h \|\nabla{{\widetilde \varrho}}(t)\|_{L^2(\Omega)}\le C\ell_hh^{2+q} \|f\|_{L^\infty(0,t;\dH q)}.$$

Thus, it suffices to bound the term ${{\widetilde \vartheta}}$; see the lemma below.

\[lem:vtht-f\] Let ${{\widetilde \vartheta}}(t)$ be defined by , and $f\in L^\infty(0,T;\dH q)$, $-1<q\leq 0$. Then with $\ell_h=|\ln h|$, there holds $$\|{{\widetilde \vartheta}}(t)\|_{L^2(\Omega)}+h\|\nabla{{\widetilde \vartheta}}(t)\|_{L^2(\Omega)}\leq Ch^{2+q}\ell_h^2\|f\|_{L^\infty(0,t;\dH q)}.$$

By and Lemma \[lem:reg-f\], we deduce that for $p=0,1$ $$\begin{aligned}
     \|{{\widetilde \vartheta}}(t)\|_{\dH p} &\leq \int_0^t \|\Etilh(t-s)\Delta_h(R_h u-P_hu)(s)\|_{\dH p} ds\\
      &\leq C\int_0^t (t-s)^{\epsilon{\alpha}/2 - 1} \tribar\Delta_h(R_h u -P_hu)(s)\tribar_{\dH {p-2+\epsilon}}ds\\
      & \leq C\int_0^t(t-s)^{\epsilon{\alpha}/2-1} \tribar R_h u(s) -P_hu(s)\tribar_{\dH {p+\epsilon}}ds.
    \end{aligned}$$ Further, using and Lemma \[lem:prh-bound\], we deduce for $p=0,1$ $$\begin{aligned}
     \|{{\widetilde \vartheta}}(t)\|_{\dH p} & \leq C h^{-\epsilon} \int_0^t(t-s)^{\epsilon{\alpha}/2-1} \| R_h u(s) -P_hu(s)\|_{\dH p}ds \\
       &  \leq C h^{2+q-p-2\epsilon} \int_0^t(t-s)^{\epsilon{\alpha}/2-1} \| u(s) \|_{\dH {2+q-\epsilon}}ds.
   \end{aligned}$$ Now by and the choice ${\epsilon}=1/\ell_h$ we get for $p=0,1$: $$\begin{aligned}
    \|{{\widetilde \vartheta}}(t)\|_{\dH p} & \leq C \epsilon^{-1} h^{2+q-p-2\epsilon}  \|f\|_{L^\infty(0,t;\dH q)}
    \int_0^t(t-s)^{\epsilon{\alpha}/2-1}t^{\epsilon \alpha/2}ds\\
  &\leq C \epsilon^{-2} h^{2+q-p-2\epsilon}\|f\|_{L^\infty(0,t;\dH q)}
   \le Ch^{2+q-p}\ell_h^2\|f\|_{L^\infty(0,t;\dH q)}.
   \end{aligned}$$ This completes the proof of the lemma.

An inspection of the proof of Lemma \[lem:vtht-f\] indicates that for $0<q <1$, one can get rid of one factor $\ell_h$. Now we can state an error estimate in $L^\infty$-norm in time.

\[thm:gal:linf\] Let $v\equiv0$, $f\in L^\infty(0,T;\dH q)$, $-1<q\leq0$, and $u$ and $u_h$ be the solutions of and with $f_h=P_hf$, respectively. Then with $\ell_h =| \ln h|$ and $t > 0$, there holds $$\| u_h(t) - u(t) \|_{L^2(\Omega)} +
  h\|\nabla(u_h(t) - u(t))\|_{L^2(\Omega)} \le Ch^{2+q} \ell_h^{2} \|f\|_{L^\infty(0,t;\dH q)}.$$

Last, we derive an error estimate in $L^2$-norm in time. To this end, we need a discrete analogue of Theorem \[thm:regeps2\], which follows from the identical proof.

\[lem:reg-d-l2\] Let $u_h$ be the solution of with $v_h=0$. Then for arbitrary $p>-1$ $$\int_0^T\normh{ \PD u_h(t)}{p}^2 +\normh{ u_h(t) }{p+2}^2 \, dt \le  \int_0^T \normh {f_h(t)}{p}^2 dt.$$

\[thm:gal:l2\] Let $v\equiv 0$, $f\in L^\infty(0,T;\dH q)$, $-1<q\leq 0$, and $u$ and $u_h$ be the solutions of and with $f_h=P_hf$, respectively. Then $$\| u_h - u \|_{L^2(0,T;L^2(\Omega))} + h\|\nabla(u_h- u)\|_{L^2(0,T;L^2(\Omega))} \le Ch^{2+q}  \|f\|_{L^2(0,T;\dH q)}.$$

We use the splitting . By Theorem \[thm:regeps2\] and Lemma \[lem:prh-bound\] $$\begin{aligned}
    \| {{\widetilde \varrho}}\|_{L^2(0,T;L^2(\Omega))} + h  \|\nabla {{\widetilde \varrho}}\|_{L^2(0,T;L^2(\Omega))} &\le Ch^{2+q}  \|u\|_{L^2(0,T;\dH {2+q})}\\
                  & \le Ch^{2+q}  \|f\|_{L^2(0,T;\dH q)}.
\end{aligned}$$ By , and Lemmas \[lem:reg-d-l2\] and \[lem:prh-bound\], we have for $p=0,\,1$: $$\begin{split}
   \int_0^T \|{{\widetilde \vartheta}}(t)\|^2_{\dH p} dt & \le C\int_0^T \normh{\Delta_h (R_h u - P_h u)(t) }{p-2}^2 dt \\
                    & \le C\int_0^T \normh{(R_h u - P_h u)(t)}{p}^2 dt\\
                    & \le C h^{4+2q-2p} \| u(t) \|_{L^2(0,T;\dH {2+q})}^2\\
                    & \le C h^{4+2q-2p}\| f(t) \|_{L^2(0,T;\dH q)}^2.
 \end{split}$$ Combing the preceding two estimates yields the desired assertion.

A Fully Discrete Scheme {#sec:fulldis}
=======================

Now we describe a fully discrete scheme for problem based on the finite difference method introduced in [@LinXu:2007]. To discretize the time-fractional derivatives, we divide the interval $[0,T]$ uniformly with a time step size $\tau=T/K$, $K\in\mathbb{N}$. We use the following discretization: $$\label{eqn:difference}
\begin{aligned}
    \frac{\pa^{\alpha}u(x,t_{n+1})}{\partial t^{\alpha}} &= \frac{1}{\Gamma(1-{\alpha})}
    \sum_{j=0}^{n} \int_{t_j}^{t_{j+1}}(t_{n+1}-s)^{-{\alpha}}\frac{\partial u(x,s)}{\partial s}\,ds \\
    &=  \frac{1}{\Gamma(1-{\alpha})}
    \sum_{j=0}^{n} \frac{u(x,t_{j+1})-u(x,t_j)}{\tau} \int_{t_j}^{t_{j+1}}(t_{n+1}-s)^{-{\alpha}}\,ds +r_{{\alpha},\tau}^{n+1}\\
    &= \frac{1}{\Gamma(2-{\alpha})} \sum_{j=0}^{n} d_{{\alpha},j} \frac{u(x,t_{n+1-j})-u(x,t_{n-j})}{\tau^{\alpha}}+ r_{{\alpha},\tau}^{n+1},
\end{aligned}$$ where $d_{{\alpha},j}=(j+1)^{1-{\alpha}}-j^{1-{\alpha}}$ with $j=0,1,2,...,n$ and $r_{{\alpha},\tau}^{n+1}$ denotes the local truncation error, which is given by $$\label{eqn:residue}
\begin{split}
   | r_{{\alpha},\tau}^{n+1} | \le C \max_{0\le t \le T} |u_{tt}(x,t)|
   \sum_{j=1}^{n} \int_{t_j}^{t_{j+1}}\frac{2s-t_j-t_{j+1}}{({t_{n+1}-s})^{{\alpha}}}\,ds +O(\tau^2).
\end{split}$$ Lin and Xu [@LinXu:2007 Lemma 3.1] showed that the truncation error $r_{\alpha,\tau}^{n+1}$ can be bounded by $$\label{eqn:residueerror}
   | r_{{\alpha},\tau}^{n+1} | \le C \max_{0\le t \le T} \left |u_{tt}(x,t) \right | \tau^{2-{\alpha}}.$$ Then the multi-term fractional derivative $\PD u(t)$ at $t=t_{n+1}$ in can be discretized by $$\label{eqn:difference2}
    \PD u(t_{n+1}) = P_\tau (\bar \pa_t) u(t_{n+1})+R_\tau^{n+1},$$ where the discrete differential operator $P_\tau(\bar\pa_t)$ is defined by $$\label{eqn:discdiff}
{\small
   P_\tau (\bar \pa_t) u(t_{n+1}) := \frac{1}{\Gamma(2-{\alpha})} \sum_{j=0}^{n}  P_j
   \frac{u(x,t_{n+1-j})-u(x,t_{n-j})}{\tau^{\alpha}},}$$ where the coefficients $\{P_j\}$ are defined by $$P_j=d_{{\alpha},j} + \sum_{i=1}^m\frac{\Gamma(2-{\alpha})b_id_{{\alpha}_i,j}\tau^{{\alpha}-{\alpha}_i}}{\Gamma(2-{\alpha}_i)},
\quad j\in \mathbb{N}.$$ Then by the local truncation error $R_\tau^{n+1}$ of the approximation $P_\tau(\bar\partial_t)u(t_{n+1})$ is bounded by $$\label{eqn:residue2}
\begin{split}
   | R_\tau^{n+1} | \le C \max_{0\le t \le T} \left |u_{tt}(x,t) \right |    \left(\tau^{2-{\alpha}} + \sum_{i=1}^m b_i \tau^{2-{\alpha}_i}\right)\le C \tau^{2-{\alpha}}\max_{0\le t \le T}
    \left |u_{tt}(x,t) \right |. \end{split}$$ By the monotonicity and convergence of $\{d_{{\alpha},j}\}$ [@LinXu:2007 equation (13)], we know that $$\label{eqn:Pmono}
P_0 > P_1 >...>0 \quad \text{and} \quad P_j \rightarrow 0 \quad \text{for} \quad j \rightarrow \infty.$$

Now we arrive at the following fully discrete scheme: find $U^{n+1}\in X_h$ such that $$\label{eqn:fully1}
    (P_\tau(\bar \pa_t)U^{n+1},\chi) + (\nabla U^{n+1},\nabla \chi) = (F^{n+1}, \chi)\quad \forall \chi\in X_h,$$ where $F^{n+1}=f(x,t_{n+1})$. Upon setting $\gamma=\Gamma(2-{\alpha})\tau^{{\alpha}}$, the fully discrete scheme is equivalent to finding $U^{n+1}\in X_h$ such that for all $\chi \in X_h$ $$\label{eqn:fully}
    P_0(U^{n+1},\chi) + \gamma(\nabla U^{n+1}, \nabla \chi)
    =\sum_{j=0}^{n-1} (P_j-P_{j+1})(U^{n-j},\chi)+P_n (U^0 ,\chi) + \gamma(F^{n+1},\chi).$$

The next result gives the stability of the fully discrete scheme.

\[lem:fullystab\] The fully discrete scheme is unconditionally stable, i.e., for all $n\in \mathbb{N}$ $$\label{eqn:fullystab}
\| U^n \|_{L^2(\Omega)} \le  \| U^0 \|_{L^2(\Omega)}+  c\max_{1\le j\le n} \|  F^{j}  \|_{L^2(\Omega)}.$$ where the constant $c$ depends only on ${\alpha}$ and $T$.

The case $n=1$ is trivial. Then the proof proceeds by mathematical induction. By noting the monotone decreasing property of the sequence $\{P_n\}$ from and choosing $\chi=U^{n+1}$ in , we deduce $$\begin{split}
   P_0 \| U^{n+1}\|_{L^2(\Omega)} &\le \sum_{j=0}^{n-1} (P_j-P_{j+1})  \| U^{n-j} \|_{L^2(\Omega)}
   +  P_n \|  U^0 \|_{L^2(\Omega)}
   +  \gamma \|  F^{n+1}  \|_{L^2(\Omega)} \\
   &\le \sum_{j=0}^{n-1} (P_j-P_{j+1})  \| U^{n-j} \|_{L^2(\Omega)}
   +  P_n \|  U^0 \|_{L^2(\Omega)}
   +  \gamma \max_{1\le j\le n+1} \|  F^{j}  \|_{L^2(\Omega)} \\
   &\le P_0 \| U^0 \|_{L^2(\Omega)}+ \left(c(P_0-P_n)+\gamma\right)\max_{1\le j\le n+1} \|  F^{j}  \|_{L^2(\Omega)} \\
\end{split}$$ Using the monotonicity of $\{P_n\}$ again gives $$c(P_0-P_n)+\gamma \le cP_0-(c P_{N} - \gamma).$$ It suffices to choose a constant $c$ such that $c P_N - \gamma>0$. By taking $\tau=T/N$, we get $$\begin{split}
P_N &= (N+1)^{1-{\alpha}}-N^{1-{\alpha}}=((T+\tau)^{1-{\alpha}}-T^{1-{\alpha}})\tau^{{\alpha}-1}
    \leq (1-{\alpha})T^{-{\alpha}}\tau^{{\alpha}} \end{split}$$ upon noting the concavity of the function $g(\tau)=(T+\tau)^{1-\alpha}$. Then by choosing $c=\Gamma(2-{\alpha})T^{\alpha}/(1-{\alpha})$ we obtain $$\begin{split}
   P_0 \| U^{n+1}\|_{L^2(\Omega)}
   \le P_0 \| U^0 \|_{L^2(\Omega)}+ cP_0\max_{1\le j\le n+1} \|  F^{j}  \|_{L^2(\Omega)}. \\
\end{split}$$ The desired result follows by dividing both sides by $P_0$.

Next we state an error estimate for the fully discrete scheme. In order to analyze the temporal discretization error, we assume the solution is sufficiently smooth.

\[thm:estfull\] Let the solution $u$ be sufficiently smooth, and $\{U^n\}\subset X_h$ be the solution of the fully discrete scheme with $U^0$ such that $$\| U^0-v\|_{L^2(\Omega)} \le Ch^2\| v\|_{\dH 2}.$$ Then there holds $$\begin{split}
\| u(t_n)-U^n\|_{L^2(\Omega)} \le  C \bigg(h^2 (\|v\|_{H^2(\Omega)}  +\| f \|_{L^\infty(0,T;L^2(\Omega))} & +
\max_{0< t\le t_n} \| u_t \|_{\dH 2} )\\
&  +\tau^{2-{\alpha}} \max_{0< t\le t_n}\| u_{tt}(t)\|_{L^2(\Omega)} \bigg).
\end{split}$$

We split the error $e^n=u(t_n)-U^n$ into $$e^n = (u(t_n)-R_h u(t_n))+(R_h u(t_n)- U^n)=: \varrho^n+\theta^n.$$ The term $\varrho^n$ can be bounded by $$\| u(t_n)-R_h u(t_n) \|_{L^2(\Omega)}\le Ch^2\| u(t_n)\|_{\dH 2}
\le Ch^2(\|v\|_{H^2(\Omega)}+\| f \|_{L^\infty(0,T;L^2(\Omega))}).$$ It suffices to bound the term $\theta^n$. By comparing and , we have the error equation $$\label{eqn:erroreq}
  (P_\tau(\bar \pa_t)\theta^n,\chi) + (\nabla\theta^n,\nabla\chi)=(\omega^n,\chi),$$ where the right hand side $\omega^n$ is given by $$\omega^n = R_h P_\tau(\bar \pa_t)u(t_n)-\PD u(t_n)= -P_\tau(\pa_t)\varrho(t_n)-R_\tau^n
  :=\omega_1^n+\omega_2^n,$$ where the truncation error $R_\tau^n$ is defined in . Using the identity $$\varrho(x,t_{j+1})-\varrho(x,t_{j})=\int_{t_{j}}^{t_{j+1}} \varrho_t(x,t) \,dt,$$ we can bound the term $\omega_1^n$ by $$\begin{split}
    \| \omega_1^n\|_{L^2(\Omega)}
    &\le C \bigg|\hspace{-0.6mm}\bigg| \sum_{j=0}^{n-1} \frac{\varrho(t_{j+1})-\varrho(t_j)}{\tau} \int_{t_j}^{t_{j+1}}(t_{n}-s)^{-{\alpha}}+\sum_{i=1}^m b_i (t_{n}-s)^{-{\alpha}_i}\,ds\bigg|\hspace{-0.6mm}\bigg|_{L^2(\Omega)}\\
    &\leq C\sum_{j=0}^{n-1} \tau^{-1} \int^{t_{j+1}}_{t_{j}} \| \varrho_t(t) \|_{L^2(\Omega)} \,dt\int_{t_j}^{t_{j+1}}(t_{n}-s)^{-{\alpha}}+\sum_{i=1}^m b_i (t_{n}-s)^{-{\alpha}_i}\,ds\\
    &\leq C h^2\max_{0< t\le t_n} \| u_t \|_{\dH 2}
    \left(  \int_0^{t_n}(t_{n}-s)^{-{\alpha}}+\sum_{i=1}^m b_i (t_{n}-s)^{-{\alpha}_i}\,ds\right)\\
    &\le C h^2\max_{0< t\le t_n} \| u_t \|_{\dH 2}.
\end{split}$$ Meanwhile, the second term $\omega_2^n$ can be bounded using . Then by the stability from Lemma \[lem:fullystab\] for the error equation , we obtain $$\begin{aligned}
  \| \theta^n  \|_{L^2(\Omega)} & \le C\bigg( \| \theta^0  \|_{L^2(\Omega)} + \max_{1\le j\le n} \| \omega_1^j \|_{L^2(\Omega)}
  + \max_{1\le j\le n} \| \omega_2^j \|_{L^2(\Omega)}\bigg)\\
      & \le C \bigg(h^2\| v \|_{\dH 2}
     + h^2 \max_{0< t\le t_n} \| u_t \|_{\dH 2} +
     \tau^{2-{\alpha}} \max_{0< t\le t_n} \|u_{tt}(t)\|_{L^2(\Omega)} \bigg).
     \end{aligned}$$

The error estimate in Theorem \[thm:estfull\] holds only if the solution $u$ is sufficiently smooth. There seems no known error estimate expressed in terms of the initial data (and right hand side) only for fully discrete schemes for nonsmooth initial data even for the single-term time-fractional diffusion equation with a Caputo fractional derivative.

Numerical Experiments {#sec:numeric}
=====================

In this part we present one- and two-dimensional numerical experiments to verify the error estimates in Sections \[sec:galerkin\] and \[sec:fulldis\]. We shall discuss the cases of a homogeneous problem and an inhomogeneous problem separately.

The case of a smooth solution
-----------------------------

Here we consider the following one-dimensional problem on the unit interval $\Omega=(0,1)$ with $0<\beta<{\alpha}<1$ $$\label{1Dnum}
\begin{split}
       \partial_t^{\alpha}u(x,t) +  \partial_t^\beta u(x,t) -\partial_{xx}^2 u(x,t)&=f,  \quad 0< x <1 \quad 0\le t\le T,\\
       u(0,t)=u(1,t)&=0, \quad  0\le t \le T, \\
       u(x,0)&=v(x), \quad  0\le x \le 1.
\end{split}$$ In order to verify the estimate in Theorem \[thm:estfull\], we first check the case that the solution $u$ is sufficiently smooth. To this end, we set initial data $v$ to $v(x)=x(1-x)$ and the source term $f$ to $f(x,t)=(2t^{2-{\alpha}}/{\Gamma(3-{\alpha})}+2t^{2-\beta}
/{\Gamma(3-\beta)})(-x^2+x)+2(1+t^2)$. Then the exact solution $u$ is given by $u(x,t)=
(1+t^2)(-x^2+x)$, which is very smooth.

In our computation, we divide the unit interval $\Omega$ into $M$ equally spaced subintervals, with a mesh size $h = 1/N$. Similarly, we fix the time step size at $\tau=1/K$. Here we choose $N$ large enough so that the space discretization error is negligible, and the time discretization error dominates. We measure the accuracy of the numerical approximation $U^n$ by the normalized $L^2$ error $\| U^n-u(t_n) \|_{L^2(\Omega)}/
\| v \|_{L^2(\Omega)}$. In Table \[tab:smooth-sol-time\], we show the temporal convergence rates, indicated in the column `rate` (the number in bracket is the theoretical rate), for three different $\alpha$ values, which fully confirm the theoretical result, cf. also Fig. \[fig:time\_error\] for the plot of the convergence rates.

     ${\alpha}$        $\tau$     $1/10$    $1/20$    $1/40$    $1/80$    $1/160$            rate
  ----------------- ------------ --------- --------- --------- --------- --------- -------------------------
   ${\alpha}=0.25$   $L^2$-norm   5.58e-4   1.73e-4   5.25e-5   1.51e-5   3.90e-6   $\approx 1.81$ ($1.75$)
   ${\alpha}=0.5$    $L^2$-norm   1.45e-3   5.11e-4   1.78e-4   6.17e-5   2.08e-5   $\approx 1.55$ ($1.50$)
   ${\alpha}=0.95$   $L^2$-norm   7.92e-3   3.79e-3   1.82e-3   8.73e-4   4.20e-4   $\approx 1.06$ ($1.05$)

  : Numerical results for the case with a smooth solution at $t=1$ with $\beta=0.2$ and ${\alpha}=0.25, 0.5, 0.95$, discretized on a uniform mesh with $h= 2^{-10}$ and $\tau =0.2\times2^{-k}$.[]{data-label="tab:smooth-sol-time"}

Homogeneous problems
--------------------

In this part we present numerical results to illustrate the spatial convergence rates in Section \[sec:galerkin\]. We performed numerical tests on the following three different initial data:

-   Smooth data: $v(x)=\sin(2\pi x)$ which belongs to $H^2(\Omega)\cap H^1_0(\Omega)$.

-   Nonsmooth data: $v(x)=\chi_{(0,1/2]}$ which lies in the space $\dH {{\epsilon}}$ for any ${\epsilon}\in [0,1/2)$.

-   Very weak data: $v(x)=\delta_{1/2}(x)$, a Dirac $\delta_{1/2}(x)$-function concentrated at $x=1/2$, which belongs to the space $\dH {- {\epsilon}}$ for any ${\epsilon}\in (1/2,1]$.

In order to check the convergence rate of the semidiscrete scheme, we discretize the fractional derivatives with a small time step $\tau$ so that the temporal discretization error is negligible. In view of the possibly singular behavior as $t \rightarrow 0$, we set the time step $\tau$ to $\tau=t/(5\times 10^4)$, with $t$ being the terminal time. For each example, we measure the error $e(t)=u(t)-u_h(t)$ by the normalized errors $\| e(t) \|_{L^2(\Omega)}/ \| v \|_{L^2(\Omega)}$ and $\| \nabla e(t) \|_{L^2(\Omega)}/ \| v \|_{L^2(\Omega)}$. The normalization enables us to observe the behavior of the error with respect to time in case of nonsmooth initial data.

### Numerical results for example (2a): smooth initial data

The numerical results show $O(h^2)$ and $O(h)$ convergence rates for the $L^2$- and $H^1$-norms of the error, respectively, for all three different $\alpha$ values, cf. Fig. \[fig:spaceerror:smooth\]. As the value of $\alpha$ increases from $0.25$ to $0.95$, the error at $t=1$ decreases accordingly, which resembles that for the single-term time-fractional diffusion equation [@JinLazarovZhou:2013].

### Numerical results for example (2b): nonsmooth initial data

For nonsmooth initial data, we are particularly interested in errors for $t$ close to zero, and thus we also present the errors at $t=0.01$ and $t=0.001$; see Table \[tab:nonsmooth1\_initial\]. The numerical results fully confirm the theoretically predicted rates for nonsmooth initial data. Further, in Table \[tab:check\_singular\] we show the $L^2$-norm of the error for fixed $h = 2^{-6}$ and $t\rightarrow 0$. We observe that the error deteriorates as $t\rightarrow0$. Upon noting $v\in \dH{1/2-{\epsilon}}$, it follows from Theorem \[thm:SG-nonsmooth\] that the error grows like $O(t^{-3{\alpha}/4})$, which agrees well with the results in Table \[tab:check\_singular\].

      $t$         $k$         $3$       $4$       $5$       $6$       $7$              rate
  ----------- ------------ --------- --------- --------- --------- --------- -------------------------
     $t=1$     $L^2$-norm   1.86e-3   4.64e-4   1.16e-4   2.87e-5   6.88e-6   $\approx 2.02$ ($2.00$)
               $H^1$-norm   4.89e-2   2.44e-2   1.22e-2   6.07e-3   2.96e-3   $\approx 1.01$ ($1.00$)
   $t=0.01$    $L^2$-norm   8.04e-3   2.00e-3   5.01e-4   1.24e-4   2.98e-5   $\approx 2.03$ ($2.00$)
               $H^1$-norm   2.31e-2   1.16e-1   5.79e-2   2.88e-2   1.40e-2   $\approx 1.01$ ($1.00$)
   $t=0.001$   $L^2$-norm   1.65e-2   4.14e-3   1.03e-3   2.56e-4   6.18e-4   $\approx 2.01$ ($2.00$)
               $H^1$-norm   5.15e-1   2.58e-1   1.29e-2   6.41e-2   3.13e-2   $\approx 1.01$ ($1.00$)

  : Numerical results for the nonsmooth case (2b) with ${\alpha}=0.5$ and $\beta=0.2$ at $t=1, 0.01, 0.001$, discretized on a uniform mesh with $h = 2^{-k}$ and $\tau =t/(5\times10^4)$. []{data-label="tab:nonsmooth1_initial"}

     $t$       1e-3      1e-4      1e-5      1e-6      1e-7      1e-8              rate
  ---------- --------- --------- --------- --------- --------- --------- -------------------------
   Case(2b)   2.56e-4   5.39e-4   1.15e-3   2.91e-3   6.77e-3   1.55e-2   $\approx-0.37 (-0.37) $

  : $L^2$-error with ${\alpha}=0.5$ and $h=2^{-6}$ for $t \to 0$ for nonsmooth initial data (2b).[]{data-label="tab:check_singular"}

### Numerical results for example (2c): very weak initial data

The numerical results show a superconvergence with a rate of $O(h^2)$ in the $L^2$-norm and $O(h)$ in the $H^1$-norm, cf. Fig. \[fig:weak\](a). This is attributed to the fact that in one-dimension the solution with the Dirac $\delta$-function as the initial data is smooth from both sides of the support point and the finite element spaces $X_h$ have good approximation property. When the singularity point $x=1/2$ is not aligned with the grid, Fig. \[fig:weak\](b) indicates an $O(h^{3/2})$ and $O(h^{1/2})$ convergence rate for the $L^2$- and $H^1$-norm of the error, respectively, which agrees with our theory.

Inhomogeneous problems
----------------------

Now we consider the inhomogeneous problem with $v\equiv0$ on the unit interval $\Omega=(0,1)$ and test the following two examples:

1.  Nonsmooth data: $f(x,t)=(\chi_{[1/2,1]}(t)+1)\chi_{[0,{1/2}]}(x)$. The jump at $x=1/2$ leads to $f(t,\cdot) \notin \dot H^1(\Omega)$; nonetheless, for any $\epsilon >0$, $f \in L^\infty(0,T;\dH {{1/2}-\epsilon})$.

2.  Very weak data: $f(x,t)=(\chi_{[1/2,1]}(t)+1)\delta_{1/2}(x)$ where $f$ involves a Dirac $\delta_{1/2}(x)$-function concentrated at $x=0.5$.

### Numerical results for example (3a)

Since the errors are bounded independently of the time, cf. Theorem \[thm:gal:linf\], we only present the errors in $L^\infty$ in time, i.e., $\|e(t)\|_{L^2(\Omega)}$ and $\| \nabla e(t)
\|_{L^2(\Omega)}$. In Table \[tab:nonsmooth1\_right\], we present the $L^2$- and $H^1$-error at $t=1$, $0.01$, and $0.001$. The numerical results agree well with our theoretical predictions, i.e., $O(h^2)$ and $O(h)$ convergence rates for the $L^2$- and $H^1$-norms of the error, respectively.

      $t$         $k$         $3$       $4$       $5$       $6$       $7$              rate
  ----------- ------------ --------- --------- --------- --------- --------- -------------------------
     $t=1$     $L^2$-norm   1.76e-3   4.40e-4   1.10e-4   2.71e-5   6.53e-6   $\approx 2.01$ ($2.00$)
               $H^1$-norm   4.72e-2   2.36e-2   1.18e-2   5.86e-3   2.86e-3   $\approx 1.01$ ($1.00$)
   $t=0.01$    $L^2$-norm   6.34e-4   1.59e-4   3.96e-5   9.82e-6   2.38e-6   $\approx 2.01$ ($2.00$)
               $H^1$-norm   1.89e-2   9.46e-3   4.72e-3   2.35e-3   1.15e-3   $\approx 1.01$ ($1.00$)
   $t=0.001$   $L^2$-norm   4.55e-4   1.15e-4   2.88e-5   1.15e-6   1.73e-6   $\approx 2.02$ ($2.00$)
               $H^1$-norm   1.45e-2   7.31e-3   3.66e-3   1.82e-3   8.88e-4   $\approx 1.01$ ($1.00$)

  : Numerical results for example (3a) with ${\alpha}=0.5$ and $\beta=0.2$ at $t=1, 0.01, 0.001$, discretized on a uniform mesh $h = 2^{-k}$ and $\tau=t/(5\times 10^{4})$. []{data-label="tab:nonsmooth1_right"}

### Numerical results for example (3b)

In Table \[tab:weak\_right\_G\] we show convergence rates at three different times, i.e., $t=1$, $0.01$, and $0.001$. Here the mesh size $h$ is chosen to be $h=1/(2^k+1)$, and thus the support of the Dirac $\delta$-function does not align with the grid. The results indicate an $O(h^{1/2})$ and $O(h^{3/2})$ convergence rate for the $H^1$- and $L^2$-norm of the error, respectively, which agrees well with the theoretical prediction. If the Dirac $\delta$-function is supported at a grid point, both $L^2$- and $H^1$-norms of the error exhibit a superconvergence of one half order, cf. Table \[tab:weak\_right\_d\]. This, however, theoretically remains to be established.

      $t$         $k$         $3$       $4$       $5$       $6$       $7$              rate
  ----------- ------------ --------- --------- --------- --------- --------- -------------------------
    $t=0.1$    $L^2$-norm   1.02e-2   4.01e-3   1.49e-3   5.35e-4   1.82e-4   $\approx 1.49$ ($1.50$)
               $H^1$-norm   3.24e-1   2.35e-1   1.65e-1   1.11e-1   6.94e-2   $\approx 0.50$ ($0.50$)
   $t=0.01$    $L^2$-norm   4.66e-3   1.91e-3   7.29e-4   2.64e-4   9.02e-5   $\approx 1.45$ ($1.50$)
               $H^1$-norm   1.54e-1   1.14e-1   8.16e-2   5.54e-2   3.47e-2   $\approx 0.55$ ($0.50$)
   $t=0.001$   $L^2$-norm   4.30e-3   1.83e-3   7.12e-4   2.61e-4   8.97e-5   $\approx 1.45$ ($1.50$)
               $H^1$-norm   1.47e-1   1.11e-1   8.05e-2   5.50e-2   3.45e-2   $\approx 0.55$ ($0.50$)

  : Numerical results for example (3b) with ${\alpha}=0.5$ and $\beta=0.2$ at $t=0.1, 0.01, 0.001$, discretized on a uniform mesh $h = 1/(2^{k}+1)$ and $\tau=t/(5\times 10^{4})$. []{data-label="tab:weak_right_d"}

      $t$         $k$         $3$       $4$       $5$       $6$       $7$              rate
  ----------- ------------ --------- --------- --------- --------- --------- -------------------------
     $t=1$     $L^2$-norm   5.35e-4   1.34e-4   3.35e-5   8.31e-6   2.01e-6   $\approx 2.01$ ($1.50$)
               $H^1$-norm   1.49e-2   7.48e-3   3.74e-3   1.86e-3   9.07e-4   $\approx 1.01$ ($0.50$)
   $t=0.01$    $L^2$-norm   6.67e-4   1.67e-4   4.17e-5   1.04e-5   2.52e-6   $\approx 2.03$ ($1.50$)
               $H^1$-norm   2.56e-2   1.29e-2   6.44e-3   3.20e-3   1.56e-3   $\approx 1.02$ ($0.50$)
   $t=0.001$   $L^2$-norm   8.19e-4   2.08e-4   5.22e-5   1.30e-5   3.19e-6   $\approx 2.02$ ($1.50$)
               $H^1$-norm   3.96e-2   2.00e-2   1.00e-3   4.98e-3   2.45e-3   $\approx 1.01$ ($0.50$)

  : Numerical results for example (3b) with ${\alpha}=0.5$ and $\beta=0.2$ at $t=1, 0.01, 0.001$, discretized on a uniform mesh with $h= 2^{-k}$ and $\tau=t/(5\times 10^{4})$. []{data-label="tab:weak_right_G"}

Examples in two-dimension
-------------------------

In this part, we present three two-dimensional examples on the unit square $\Omega=(0,1)^2$.

1.  Nonsmooth initial data: $v=\chi_{(0,1/2)\times(0,1)}$ and $f\equiv0$.

2.  Very weak initial data: $v=\delta_\Gamma$ with $\Gamma$ being the boundary of the square $[1/4,3/4]^2$ and $\langle \delta_\Gamma,\phi\rangle=\int_\Gamma\phi(s)\,ds$. By Hölder’s inequality and the continuity of the trace operator from $\dot H^{{1/2}+\epsilon}(\Omega)$ to $L^2(\Gamma)$ [@AdamsFournier:2003], we deduce $\delta_\Gamma \in H^{-1/2-\epsilon}(\Omega)$.

3.  Nonsmooth right hand side: $f(x,t)=(\chi_{[1/20,1/10]}(t)+1)\chi_{(0,1/2)\times(0,1)}(x)$ and $v\equiv0$.

To discretize the problem, we divide each direction into $N=2^k$ equally spaced subintervals, with a mesh size $h=1/N$ so that the domain $[0,1]^2$ is divided into $N^2$ small squares. We get a symmetric mesh by connecting the diagonal of each small square.

The numerical results for example (4a) are shown in Table \[tab:nonsmooth2D\], which agree well with Theorem \[thm:SG-nonsmooth\], with a rate $O(h^2)$ and $O(h)$, respectively, for the $L^2$- and $H^1$-norm of the error. Interestingly, for example (4b), both the $L^2$-norm and $H^1$-norm of the error exhibit super-convergence, cf. Table \[tab:weak2D\_G\]. The numerical results for example (4c) confirm the theoretical results; see Table \[tab:2Dnonsmooth\_right\]. The solution profiles for examples (4b) and (4c) at $t=0.1$ are shown in Fig. \[fig:2D\_solution\], from which the nonsmooth region of the solution can be clearly observed.

      $t$         $k$         $3$       $4$       $5$       $6$       $7$              rate
  ----------- ------------ --------- --------- --------- --------- --------- -------------------------
    $t=0.1$    $L^2$-norm   5.25e-3   1.35e-3   3.38e-4   8.24e-5   1.98e-5   $\approx 2.06$ ($2.00$)
               $H^1$-norm   9.10e-2   4.53e-2   2.25e-2   1.09e-2   4.99e-3   $\approx 1.04$ ($1.00$)
   $t=0.01$    $L^2$-norm   1.25e-2   3.23e-3   8.09e-4   1.97e-4   4.65e-5   $\approx 2.05$ ($2.00$)
               $H^1$-norm   2.18e-1   1.08e-1   5.35e-2   2.62e-2   1.27e-2   $\approx 1.05$ ($1.00$)
   $t=0.001$   $L^2$-norm   3.02e-2   7.84e-3   1.97e-3   4.81e-4   1.16e-4   $\approx 2.03$ ($2.00$)
               $H^1$-norm   5.30e-1   2.64e-1   1.31e-1   6.38e-2   3.14e-2   $\approx 1.04$ ($1.00$)

  : Numerical results for (4a) with ${\alpha}=0.5$ and $\beta=0.2$ at $t=0.1, 0.01, 0.001$, discretized on a uniform mesh, $h = 2^{-k}$ and $\tau=t/10^{4}$. []{data-label="tab:nonsmooth2D"}

      $t$         $k$         $3$       $4$       $5$       $6$       $7$              rate
  ----------- ------------ --------- --------- --------- --------- --------- -------------------------
    $t=0.1$    $L^2$-norm   1.18e-2   3.18e-3   8.41e-4   2.18e-4   5.41e-5   $\approx 1.92$ ($1.50$)
               $H^1$-norm   2.25e-1   1.13e-1   6.60e-2   3.40e-2   1.66e-2   $\approx 0.92$ ($0.50$)
   $t=0.01$    $L^2$-norm   2.82e-2   7.62e-3   2.28e-3   5.26e-4   1.25e-4   $\approx 1.95$ ($2.00$)
               $H^1$-norm   5.66e-1   3.09e-1   1.65e-1   8.52e-2   4.19e-2   $\approx 0.94$ ($1.00$)
   $t=0.001$   $L^2$-norm   6.65e-2   1.83e-3   4.98e-3   1.33e-3   3.30e-4   $\approx 1.91$ ($2.00$)
               $H^1$-norm   1.66e0    8.93e-1   4.75e-1   2.43e-1   1.21e-1   $\approx 0.95$ ($1.00$)

  : Numerical results for example (4b) with ${\alpha}=0.5$ and $\beta=0.2$ at $t=0.1, 0.01, 0.001$ for a uniform mesh with $h = 2^{-k}$ and $\tau=t/10^{4}$. []{data-label="tab:weak2D_G"}

      $t$         $k$         $3$       $4$       $5$       $6$       $7$              rate
  ----------- ------------ --------- --------- --------- --------- --------- -------------------------
    $t=0.1$    $L^2$-norm   2.28e-3   5.86e-4   1.47e-4   3.58e-5   7.91e-6   $\approx 2.07$ ($2.00$)
               $H^1$-norm   3.97e-2   1.97e-2   9.77e-3   4.76e-3   2.13e-3   $\approx 1.06$ ($1.00$)
   $t=0.01$    $L^2$-norm   1.06e-3   2.73e-4   6.86e-5   1.67e-6   3.70e-6   $\approx 2.06$ ($2.00$)
               $H^1$-norm   1.85e-2   9.18e-3   4.56e-3   2.22e-3   9.94e-3   $\approx 1.06$ ($1.00$)
   $t=0.001$   $L^2$-norm   8.66e-4   2.28e-4   5.75e-5   1.40e-6   3.11e-6   $\approx 2.04$ ($2.00$)
               $H^1$-norm   1.56e-2   7.82e-3   3.88e-3   1.90e-3   8.47e-4   $\approx 1.05$ ($1.00$)

  : Numerical results for example (4c) with ${\alpha}=0.5$ and $\beta=0.2$ at $t=0.1, 0.01, 0.001$ for a uniform mesh with $h= 2^{-k}$ and $\tau=t/10^{4}$. []{data-label="tab:2Dnonsmooth_right"}

Concluding remarks
==================

In this work, we have developed a simple numerical scheme based on the Galerkin finite element method for a multi-term time fractional diffusion equation which involves multiple Caputo fractional derivatives in time. A complete error analysis of the space semidiscrete Galerkin scheme is provided. The theory covers the practically very important case of nonsmooth initial data and right hand side. The analysis relies essentially on some new regularity results of the multi-term time fractional diffusion equation. Further, we have developed a fully discrete scheme based on a finite difference discretization of the Caputo fractional derivatives. The stability and error estimate of the fully discrete scheme were established, provided that the solution is smooth. The extensive numerical experiments in one- and two-dimension fully confirmed our convergence analysis: the empirical convergence rates agree well with the theoretical predictions for both smooth and nonsmooth data.

Acknowledgements {#acknowledgements .unnumbered}
================

The research of B. Jin has been supported by US NSF Grant DMS-1319052, R. Lazarov was supported in parts by US NSF Grant DMS-1016525 and also by Award No. KUS-C1-016-04, made by King Abdullah University of Science and Technology (KAUST), and Y. Liu was supported by the Program for Leading Graduate Schools, MEXT, Japan.

Proof of Lemma \[lem:MMLcm\] {#app:MMLcm}
============================

First, we define an auxiliary function $v_j(t)$ by $$\hat v_j(z)=\mathcal{L}(v_j)
  = \frac{z^{{\alpha}-1}+\sum_{k=1}^m b_k z^{{\alpha}_k-1}}{z^{\alpha}+\sum_{k=1}^m b_k z^{{\alpha}_k}+{\lambda}_j}.$$ Now by the property of the Laplace transform $f(0^+)=\lim_{z\to \infty} z \widehat{f}(z)$, we obtain $v_j(0+)=1$. The function $\bar E_{\vecal}^j(t)$ is the inverse Laplace integral of $\widehat {\bar E_{\vecal}^j} =(z^{\alpha}+\sum_{k=1}^m b_k z^{{\alpha}_k}+{\lambda}_j)^{-1}$, i.e. $$\label{L1}
  \bar E_{\vecal}^j(t)=\frac{1}{2\pi\mathrm{i}} \int_{Br}e^{zt}\frac{1}{{z^{\alpha}+\sum_{k=1}^m b_k z^{{\alpha}_k}+{\lambda}_j}}\,dz,$$ where $Br=\{z;\ \text{Re}~~z=\sigma,\ \sigma>0\}$ is the Bromwich path. The function $\widehat{\bar E_{\vecal}^j}(z)$ has a branch point $0$, so we cut off the negative part of the real axis. Note that the function $z^{\alpha}+\sum_{k=1}^m b_k z^{{\alpha}_k}+{\lambda}_j$ has no zero in the main sheet of the Riemann surface including its boundaries on the cut. Indeed, if $z=\varrho e^{i\theta}$, with $\rho>0$, $\theta\in(-\pi,\pi)$, then $$\Im\left\{z^{\alpha}+\sum_{k=1}^m b_k z^{{\alpha}_k}+{\lambda}_j\right\}=\varrho^{\alpha}\sin{\alpha}\theta+\sum_{k=1}^m b_k \varrho^{{\alpha}_k}\sin{\alpha}_k\theta\neq 0,\ \ \forall \theta\neq 0,$$ since $\sin \alpha \theta$ and $\sin\alpha_k\theta$ have the same sign for any $\theta\in(-\pi,\pi)$ and $b_k>0$. Hence, $ \bar E_{\vecal}^j(t)$ can be found by bending the Bromwich path into the Hankel path $Ha({\epsilon})$, which starts from $-\infty$ along the lower side of the negative real axis, encircles the disc $|s|={\epsilon}$ counterclockwise and ends at $-\infty$ along the upper side of the negative real axis. Then by taking ${\epsilon}\to 0$, we obtain $$  \bar E_{\vecal}^j(t)= \int_0^\infty e^{-rz}K_n(r)\,dr,$$ where $$  K_n(r)=-\frac{1}{\pi}\Im\left\{\left.\frac{1}
  {z^{\alpha}+\sum_{k=1}^m b_k z^{{\alpha}_k}+{\lambda}_j}\right|_{z=r e^{i\pi}}\right\}.$$ It is easy to check $$K_n(r)=\frac{1}{\pi} \frac{r^{{\alpha}}\sin{\alpha}\pi+\sum_{k=1}^m b_kr^{{\alpha}_k}\sin{\alpha}_k\pi}
 {(r^{{\alpha}}\cos{\alpha}\pi+\sum_{k=1}^m b_kr^{{\alpha}_k}\cos{\alpha}_k\pi+{\lambda}_j)^2
 +(r^{{\alpha}}\sin{\alpha}\pi+\sum_{k=1}^m b_kr^{{\alpha}_k}\sin{\alpha}_k\pi)^2}
$$ which is greater than zero for all $r>0$. Therefore, $\bar E_{\vecal}^j(t)$ is completely monotone. A similar argument shows that $v_j(t)$ is also completely monotone. Consequently, $$\int_0^T |\bar E_{\vecal}^j(t)|\,dt=\int_0^T \bar E_{\vecal}^j(t)\,dt=
-\frac{1}{{\lambda}_j}\int_0^T v'_j(t)\,dt=\frac{1}{{\lambda}_n}(1-v_j(T))<\frac{1}{{\lambda}_n},$$ which concludes the proof of the lemma.
