Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions

Abstract

An initial-boundary value problem with a Caputo time derivative of fractional order $\alpha \in (0,1)$ is considered, solutions of which typically exhibit a singular behaviour at an initial time. For this problem, we give a simple framework for the analysis of the error of L1-type discretizations on graded and uniform temporal meshes in the $L_\infty$ and $L_2$ norms. This framework is employed in the analysis of both finite difference and finite element spatial discretiztions. Our theoretical findings are illustrated by numerical experiments.

1. Introduction

The purpose of this paper is to give a simple framework for the analysis of the error in the $L_\infty (\Omega )$ and $L_2(\Omega )$ norms for L1-type discretizations of the fractional-order parabolic problem

$$\begin{equation} \begin{array}{l} D_t^{\alpha }u+\mathcal{L}u=f(x,t)\quad \mbox{for}\;\;(x,t)\in \Omega \times (0,T],\\[5.69046pt] u(x,t)=0\quad \mbox{for}\;\;(x,t)\in \partial \Omega \times (0,T],\qquad u(x,0)=u_0(x)\quad \mbox{for}\;\;x\in \Omega . \end{array} {\tag{1.1}} \end{equation}$$

This problem is posed in a bounded Lipschitz domain $\Omega \subset \mathbb{R}^d$ (where $d\in \{1,2,3\}$). The operator $D_t^\alpha$, for some $\alpha \in (0,1)$, is the Caputo fractional derivative in time defined 2 by

$$\begin{equation} D_t^{\alpha } u(\cdot ,t) := \frac{1}{\Gamma (1-\alpha )} \int _{0}^t(t-s)^{-\alpha }\, \partial _s u(\cdot , s)\, ds \qquad \text{for }\ 0<t \le T, {\tag{1.2}} \end{equation}$$

where $\Gamma (\cdot )$ is the Gamma function, and $\partial _s$ denotes the partial derivative in $s$. The spatial operator $\mathcal{L}$ is a linear second-order elliptic operator:

$$\begin{equation} \mathcal{L}u := \sum _{k=1}^d \Bigl \{-\partial _{x_k}\!(a_k(x)\,\partial _{x_k}\!u) + b_k(x)\, \partial _{x_k}\!u \Bigr \}+c(x)\, u, {\tag{1.3}} \end{equation}$$

with sufficiently smooth coefficients $\{a_k\}$, $\{b_k\}$ and $c$ in $C(\bar{\Omega })$, for which we assume that $a_k>0$ in $\bar{\Omega }$, and also either $c\ge 0$ or $c-\frac{1}{2}\sum _{k=1}^d\partial _{x_k}\!b_k\ge 0$. All our results also apply to the case $\mathcal{L}=\mathcal{L}(t)$, while some remain valid for a more general uniformly-elliptic $\mathcal{L}$ (i.e., with mixed second-order derivatives); see Remark 3.3.

Throughout the paper, it will be assumed that there exists a unique solution of this problem in $C(\bar{\Omega }\times [0,T])$ such that $|\partial _t^l u(\cdot ,t)|\lesssim 1+t^{\alpha -l}$ for $l=0,1,2$ (the notation $\lesssim$ is rigourously defined in the final paragraph of this section). This is a realistic assumption, satisfied by typical solutions of problem 1.1, in contrast to a stronger assumption $|\partial ^l u(\cdot ,t)|\lesssim 1$ frequently made in the literature (see, e.g., references in 8, Table 1.1). Indeed, 21, Theorem 2.1 shows that if a solution $u$ of 1.1 is less singular than we assume (in the sense that $|\partial _t^l u(\cdot ,t)|\lesssim 1+t^{\gamma -l}$ for $l=0,1,2$ with any $\gamma >\alpha$), then the initial condition $u_0$ is uniquely defined by the other data of the problem, which is clearly too restrictive. At the same time, our results can be easily applied to the case of $u$ having no singularities or exhibiting a somewhat different singular behaviour at $t=0$.

We consider L1-type schemes for problem 1.1, which employ the discetization of $D^\alpha _tu$ defined, for $m=1,\ldots ,M$, by

$$\begin{equation} \delta _t^{\alpha } U^m := \frac{1}{\Gamma (1-\alpha )} \sum _{j=1}^m \delta _t U^j\!\int _{t_{j-1}}^{t_j}\!\!(t_m-s)^{-\alpha }\, ds, \qquad \delta _t U^j:=\frac{U^j-U^{j-1}}{t_j-t_{j-1}}, {\tag{1.4}} \end{equation}$$

when associated with the temporal mesh $0=t_0<t_1<\ldots <t_M=T$ on $[0,T]$. Similarly to 22, our main interest will be in graded temporal meshes as they offer an efficient way of computing reliable numerical approximations of solutions that are singular at $t=0$. We shall also consider uniform temporal meshes, as although the latter have lower convergence rates near $t=0$, they have been shown to be first-order accurate for $t\gtrsim 1$ 59.

Novelty

We present a new framework for the estimation of the error whenever an L1 scheme is used on graded or uniform temporal meshes. This framework is simple, applies to both finite difference and finite element spatial discretizations, and works for error estimation in both $L_2(\Omega )$ and $L_\infty (\Omega )$ norms. It easily extends to general elliptic operators $\mathcal{L}=\mathcal{L}(t)$, as well as quasi-uniform and quasi-graded temporal meshes. Naturally, it yields versions of some previously-known error bounds as particular cases. It is also used here to establish entirely new results.

Graded meshes for problems of type 1.1 for the case $d=1$ were recently considered in 22, where maximum norm error bounds are obtained for finite difference discretizations. In comparison, our analysis deals with temporal-discretization errors on graded meshes in an entirely different and substantially more concise way. To be more precise, we use more intuitive integral representations of the temporal truncation errors; see Lemma 2.3. Once error bounds on graded meshes are established for a paradigm problem without spatial derivatives, they seamlessly extend to finite difference and finite element spatial discretizations of 1.1 for any $d\ge 1$. Our results on graded meshes are new for finite element discretizations, as well as for finite difference discretizations for $d>1$.

The convergence behaviour of the L1 method on uniform temporal meshes is well understood. In particular, for finite element spatial discretizations, the errors in the $L_2(\Omega )$ norm have been estimated in 9 using Laplace transform techniques (for $\mathcal{L}=-\triangle$ and $f=0$). For finite difference discretizations for $d=1$, a similar error bound in the maximum norm was established in 5. Within our theoretical framework, we easily get versions of error bounds of 9 and 5. Furthermore, we give error bounds for finite element discretizations in the $L_\infty (\Omega )$ norm on uniform temporal meshes, which appear to be entirely new. (Some error bounds in the $L_\infty$ norm for linear-finite-element spatial semidiscretizations are given in 12.)

Our approach to uniform meshes is very similar to the case of graded meshes. The main difference is in that now we employ a more subtle stability property of the discrete fractional-derivative operator $\delta _t^\alpha$ from 5, a version of which is also given in 11; see Lemma 2.1*. Additionally, we give a considerably shorter and more intuitive proof of this stability result. This new proof relies on a simple barrier function, and may be of independent interest; see Appendix A.

Outline

We start by presenting, in §2, a paradigm for the temporal-error analysis using a simplest example without spatial derivatives. This error analysis is extended in §3 to temporal semidiscretizations of 1.1. Full discretizations that employ finite differences and finite elements are, respectively, addressed in §4 and §5. Finally, the assumptions on the derivatives of the exact solution are discussed in §6, and our theoretical findings are illustrated by numerical experiments in §7.

Notation

We write $a\simeq b$ when $a \lesssim b$ and $a \gtrsim b$, and $a \lesssim b$ when $a \le Cb$ with a generic constant $C$ depending on $\Omega$, $T$, $u_0$, and $f$, but not on the total numbers of degrees of freedom in space or time. Also, for $1 \le p \le \infty$, and $k \ge 0$, we shall use the standard norms in the spaces $L_p(\Omega )$ and the related Sobolev spaces $W_p^k(\Omega )$, while $H^1_0(\Omega )$ is the standard space of functions in $W_2^1(\Omega )$ vanishing on $\partial \Omega$.

2. Paradigm for the temporal-discretization error analysis

2.1. Graded temporal mesh

Throughout the paper, we shall frequently consider the graded temporal mesh $\{t_j=T(j/M)^r\}_{j=0}^M$ with some $r\ge 1$ (while $r=1$ generates a uniform mesh). For this mesh, a calculation shows that

$$\begin{equation} \tau _j:=t_j-t_{j-1}\simeq M^{-1}\,t_j^{1-1/r}\qquad \mbox{for  }j=1,\ldots ,M. {\tag{2.1}} \end{equation}$$

This follows from $\tau _1=t_1\simeq M^{-r}$ for $j=1$, and $t_j\le 2^{r} t_{j-1}$ for $j\ge 2$.

Note that all results of the paper immediately apply to a quasi-graded mesh defined by $\{t_j=T(\xi _j)^r\}_{j=0}^M$, where $\{\xi _j\}_{j=0}^M$ is a quasi-uniform mesh on $[0,1]$.

2.2. Stability properties of the discrete fractional operator $\delta _t^\alpha$

The definition 1.4 of $\delta _t^\alpha$ can be rewritten as

$$\begin{align} \delta _t^\alpha V^m &= \underbrace{{\kappa _{m,m}}}_{{}>0} V^m-\sum _{j=1}^m \underbrace{(\kappa _{m,j}-\kappa _{m,j-1})}_{{}>0}V^{j-1}, {\tag{2.2a}}\\ \kappa _{m,j}&:= \frac{\tau _j^{-1}}{\Gamma (1-\alpha )} \int _{t_{j-1}}^{t_{j}}\!\! (t_m-s)^{-\alpha }\, ds \quad \mbox{for}\;\; j=1,\ldots ,m,\quad \;\; \kappa _{m,0}:=0. {\tag{2.2b}} \end{align}$$Here $\kappa _{m,j}$ for $j\ge 1$ is the average of the function $\{\Gamma (1-\alpha )\}^{-1}(t_m-s)^{-\alpha }$ on the interval $s\in (t_{j-1},t_j)$, so $\kappa _{m,j-1}\le \kappa _{m,j}$ for all admissible $j$ and $m$.

Lemma 2.1.

(i)

For any $\{V^j\}_{j=0}^M$ on an arbitrary mesh $\{t_j\}_{j=0}^M$, one has$$\begin{equation*} |V^m-V^0|\lesssim \max _{j=1,\ldots ,m}\bigl \{t_j^{\alpha }\, |\delta _t^\alpha V^j|\bigr \} \qquad \mbox{for}\;\; m=1,\ldots , M. \end{equation*}$$

(ii)

If $V^0=0$ and $\delta _t^\alpha |V^j|\le |F^j|$ for $j=1,\ldots ,M$, then$$\begin{equation*} |V^m|\lesssim \max _{j=1,\ldots ,m}\bigl \{t_j^{\alpha }\, |F^j|\bigr \} \qquad \mbox{for}\;\; m=1,\ldots , M. \end{equation*}$$

Proof.

(i)

Let $W^j:=V^j-V^0$; then $W^0=0$, while $\delta _t^\alpha W^j=\delta _t^\alpha V^j=:F^j$, so we need to prove that $|W^m|\lesssim \max _{j\le m}\{t_j^{\alpha }\,|F^j|\}$. Let $\max _{j\le m}|W^j|=|W^n|$ for some $1\le n\le m$. Then, by 2.2a combined with $W^0=0$, one gets$$\begin{equation} \underbrace{{\kappa _{n,n}}}_{{}>0} |W^n|- \sum _{j=2}^n \underbrace{(\kappa _{n,j}-\kappa _{n,j-1})}_{{}>0}|W^n|\le |F^n| \quad \Rightarrow \quad |W^n|\le \kappa _{n,1}^{-1}\,|F^n|. {\tag{2.3}} \end{equation}$$

Next, recalling 2.2b, and also using $(t_n-s)^{-\alpha }\ge t_n^{-\alpha }$ on $(0,t_1)$, one concludes that $\kappa _{n,1}\gtrsim t_n^{-\alpha }$.

(ii)

Let $W^0=0$ and $\delta _t^\alpha W^j= |F^j|$ for $j=1,\ldots ,M$. Then $0\le |V^m|\le W^m$ (as $\delta _t^\alpha$ is associated with an $M$-matrix), while $|W^m|\lesssim \max _{j=1,\ldots ,m}\bigl \{t_j^{\alpha }\, |F^j|\bigr \}$ by the result of part (i). The desired assertion follows.

■

To deal with uniform temporal meshes, we employ a more subtle stability result.

Lemma 2.1* (5).

Let $r=1$ and $\tau :=TM^{-1}$. Given $\gamma \in (0,\alpha ]$, if $V^0=0$ and $|\delta _t^\alpha V^j|\lesssim \tau ^\gamma t_j^{-\gamma -1}$ for $j=1,\ldots ,M$, then $|V^j|\lesssim t_j^{\alpha -1}$ for $j=1,\ldots , M$.

Proof.

The desired assertion follows from 5, Lemma 3 with $\beta =1+\gamma$; see also 11, Theorem 3.3 for a similar result. We give an alternative (substantially shorter) proof in Appendix A.

■

The next lemma will be useful when dealing with Ritz projections while estimating the errors of finite element discretizations in §5.

Lemma 2.2.

Let $\{V^j\}_{j=0}^M\in \mathbb{R}^{M+1}$ and $\{\lambda ^j\}_{j=1}^M\in \mathbb{R}^{M}$, and let $\bar{\lambda }=\bar{\lambda }(t)$ be a piecewise-constant left-continuous function defined by $\bar{\lambda }(t)=\lambda ^j$ for $t\in (t_{j-1},t_j]$, $j=1,\ldots ,M$. Then, with the notation $J^{1-\alpha }v(t):=\{\Gamma (1-\alpha )\}^{-1}\!\! \int _{0}^t(t-s)^{-\alpha } v( s)\, ds$,

$$\begin{equation} \delta _t^\alpha V^j \le J^{1-\alpha }\bar{\lambda }(t_j)\quad \forall \,j\ge 1 \quad \;\;\Rightarrow \quad \;\; V^m-V^0\le \sum _{j=1}^m\tau _j\, \lambda ^j\quad \forall \,m\ge 0.{\tag{2.4}} \end{equation}$$

Proof.

Let $\Lambda ^j:=V^0+\int _0^{t_{j}}\bar{\lambda }\,dt$ so that $\lambda ^j=\delta _t\Lambda ^j$. Now, $J^{1-\alpha }\bar{\lambda }(t_j)=\delta _t^\alpha \Lambda ^j$, so we get $M$ equations $\delta _t^\alpha V^j \le \delta ^\alpha _t\Lambda ^j$ for $j=1,\ldots ,M$. Augmenting these equations by $V^0=\Lambda ^0$, we get the matrix relation $A\vec{V}\le A\vec{\Lambda }$ for the column vectors $\vec{V}:=\{V^j\}_{j=0}^M$ and $\vec{\Lambda }:=\{\Lambda ^j\}_{j=0}^M$ with an inverse-monotone $(M+1)\times (M+1)$ matrix $A$. (The latter follows from $A$ being diagonally dominant, with the entries $A_{ij}\le 0$ for $i\neq j$ in view of 2.2a.) Consequently, $\vec{V}\le \vec{\Lambda }$, which immediately yields the desired assertion.

■

2.3. Error estimation for a simplest example (without spatial derivatives)

It is convenient to illustrate our approach to the estimation of the temporal-discretization error using a very simple example. Consider a fractional-derivative problem without spatial derivatives together with its discretization:

$$\begin{align} D_t^\alpha u(t)&=f(t)&&\hspace{-1.6cm}\mbox{for}\;\;t\in (0,T],&&\hspace{-0.8cm} u(0)=u_0, \tag{2.5a}\\ \delta _t^\alpha U^j&=f(t_j)&&\hspace{-1.6cm}\mbox{for}\;\;j=1,\ldots ,M,&&\hspace{-0.8cm} U^0=u_0. \tag{2.5b} \end{align}$$

Throughout this subsection, with slight abuse of notation, $\partial _t$ will be used for $\frac{d}{dt}$, while $\delta _tu(t_j):=\tau _j^{-1}[u(t_j)-u(t_{j-1})]$ (similarly to $\delta _t$ in 1.4).

Lemma 2.3.

Let $\{t_j=T(j/M)^r\}_{j=0}^M$ for some $r\ge 1$. Then for $u$ and $U^j$ that satisfy 2.5, one has

$$\begin{equation*} |u(t_m)-U^m| \lesssim \max _{j=1,\ldots ,m}\psi ^j, \end{equation*}$$

where $m=1,\ldots ,M$, and

$$\begin{align} \psi ^1&:=\tau _1^\alpha \sup _{s\in (0,t_1)}\!\!\bigl (s^{1-\alpha }|\delta _tu(t_1)-\partial _s u(s)|\bigr ), {\tag{2.6a}}\\ \psi ^j&:=\tau _j^{2-\alpha } \,t_j^\alpha \sup _{s\in (t_{j-1},t_j)}\!\!\!|\partial _s^2u(s)|\qquad \mbox{for}\;\;j\ge 2. {\tag{2.6b}} \end{align}$$

Proof.

Using the standard piecewise-linear Lagrange interpolant $u^I$ of $u$, let

$$\begin{equation*} \chi :=u-u^I \quad \Rightarrow \quad |\chi (s)|\le \underbrace{\tau _j(t_j-s)}_{{}\le \tau _j^2}\,\sup _{s\in (t_{j-1},t_j)}\!\!\!|\partial _s^2u|\quad \mbox{for }s\in [t_{j-1},t_j]. \end{equation*}$$

As $\chi$ will appear in the truncation error, it is useful to note that, in view of 2.6b,

$$\begin{equation} |\chi (s)|\le \psi ^j\, \tau _j^\alpha \,t_j^{-\alpha }\min \{1,{(t_j-s)/\tau _j}\}\quad \;\mbox{for}\;\; s\in (t_{j-1},t_j),\;j\ge 2. {\tag{2.7a}} \end{equation}$$On $(0,t_1)$, one has $\chi '(s)=\partial _s u(s)-\delta _tu(t_1)$, which, combined with 2.6a, yields