Implementation of high-order,
discontinuous Galerkin time stepping
for fractional diffusion problems

The discontinuous Galerkin (dG) method provides an effective numerical procedure for the time integration of diffusion problems. In the mid-1980s, Eriksson, Johnson and Thomée [2] provided the first detailed error analysis, which has been subsequently extended and refined by numerous authors [[]and references therein]MakridakisNochetto2006,SchmutzWihler2019,SchotzauSchwab2001. The dG method has also proved effective for time stepping of fractional diffusion problems [9, 12] of the form

$$\partial_{t}u+\partial_{t}^{1-\alpha}Au=f(t)\quad\text{for $0<t\leq T$, with $u(0)=u_{0}$.}$$

(1)

Here, $A$ is a linear, second-order, elliptic partial differential operator over a spatial domain $\Omega$, subject to a homogeneous Dirichlet boundary condition $u=0$ on $\partial\Omega$. (Our notation suppresses the dependence of $u$ and $f$ on the spatial variables.) The fractional diffusion exponent is assumed to satisfy $0<\alpha<1$ (the sub-diffusive case), and the fractional time derivative is understood in the Riemann–Liouville sense: for $t>0$ and $\mu>0$,

$$\partial_{t}^{\mu}v=\frac{\partial}{\partial t}\int_{0}^{t}\omega_{\mu}(t-s)v(s)\,ds\quad\text{where}\quad\omega_{\mu}(t)=\frac{t^{\mu-1}}{\Gamma(\mu)}.$$

The partial integro-differential equation (1) arises in a variety of physical models [3, 11] of diffusing particles whose behaviour is described by a continuous-time random walk for which the waiting-time distribution is a power law that decays like $1/t^{1+\alpha}$. The expected waiting time is therefore infinite, and the mean-square displacement is proportional to $t^{\alpha}$. Standard Brownian motion is recovered in the limit as $\alpha\to 1$, when (1) reduces to the classical diffusion equation.

Our main concern in the present work is with the practical implementation of dG time stepping for (1), and in particular with the accurate evaluation of certain coefficients $H^{n,n-\ell}_{ij}$ used during the $n$th step. Section 2 introduces the dG method for the fractional ODE case of (1), in which the operator $A$ is replaced by a scalar $\lambda>0$. We will see in the simplest, lowest-order scheme, when the dG solution is piecewise-constant in time, that

$$H^{n,0}_{11}=\int_{t_{n-1}}^{t_{n}}\frac{d}{dt}\biggl{(}\int_{t_{n-1}}^{t}\omega_{\alpha}(t-s)\,ds\biggr{)}\,dt$$

and

$$H^{n,n-\ell}_{11}=\int_{t_{n-1}}^{t_{n}}\frac{d}{dt}\biggl{(}\int_{t_{\ell-1}}^{t_{\ell}}\omega_{\alpha}(t-s)\,ds\biggr{)}\,dt\quad\text{for $1\leq\ell\leq n-1$,}$$

where $0=t_{0}<t_{1}<t_{2}<\cdots$ are the discrete time levels. We easily verify that $H^{n,0}_{11}=\omega_{\alpha+1}(k_{n})=k_{n}^{\alpha}/\Gamma(\alpha+1)$, for a step-size $k_{n}=t_{n}-t_{n-1}$, and

	$\displaystyle H^{n,n-\ell}_{11}$	$\displaystyle=\omega_{\alpha+1}(t_{n}-t_{\ell-1})-\omega_{\alpha+1}(t_{n}-t_{\ell})$		(2)
		$\displaystyle\qquad{}-\omega_{\alpha+1}(t_{n-1}-t_{\ell-1})+\omega_{\alpha+1}(t_{n-1}-t_{\ell}),$

but for higher-order schemes the coefficients become progressively more complicated. Although the $H^{n,n-\ell}_{ij}$ can always be evaluated via repeated integration by parts, the resulting expressions are likely to suffer from roundoff when evaluated in floating-point arithmetic if $n-\ell$ is large. Consider just the lowest order case (2) with uniform time steps $t_{n}=nk$, so that

$$H^{n,n-\ell}_{11}=k^{\alpha}\bigl{[}\omega_{\alpha+1}(n-\ell+1)-2\omega_{\alpha+1}(n-\ell)+\omega_{\alpha+1}(n-\ell-1)\bigr{]}.$$

Since the factor in square brackets is a second-difference of $\omega_{\alpha+1}$, its magnitude decays like $(n-\ell)^{\alpha-2}$ as $n-\ell$ increases, but the individual terms grow like $(n-\ell)^{\alpha}$.

We are therefore led to evaluate the coefficients $H^{n,n-\ell}_{ij}$ via quadratures with positive weights. No special techniques are needed for $\ell\leq n-2$, but when $\ell=n$ or $n-1$ we must deal with weakly singular integrands. In Section 3, we show how certain substitutions reduce the problem to dealing with integrands that are either smooth, or are products of smooth functions and standard Jacobi weight functions. Similar substitutions, known as Duffy transformations [1], have long been used to compute singular integrals arising in the boundary element method.

Section 4 introduces a spatial discretisation for the fractional PDE (1) and describes the structure of the linear system that must be solved at each time step. In Section 5, we specialise the expressions for the coefficients by choosing Legendre polynomials as the shape functions employed in the dG time stepping.

Section 6 describes a post-processing technique that, when applied to the dG solution $U$, produces a more accurate approximate solution $\widehat{U}$, known as the reconstruction [6] of $U$. If $U$ is a piecewise polynomial of degree at most $r-1$, then $\widehat{U}$ is a piecewise polynomial of degree at most $r$. For a classical diffusion problem, both $U$ and $\widehat{U}$ are known to be quasi-optimal, that is, accurate of order $k^{r}$ and $k^{r+1}$, respectively. Thus, it is natural to ask what happens in the fractional-order case, and we investigate this question in numerical experiments reported in Section 7.

Our central concern is present already in the zero-dimensional case when we replace the elliptic operator $A$ with a scalar $\lambda\geq 0$, so that the solution $u(t)$ is a real-valued function satisfying the fractional ODE

$$u^{\prime}+\lambda\partial_{t}^{1-\alpha}u=f(t)\quad\text{for $0<t\leq T$, with $u(0)=u_{0}$.}$$

(3)

For the time discretisation, we introduce a grid

$$0=t_{0}<t_{1}<t_{2}<\cdots<t_{N}=T,$$

and form the vector $\boldsymbol{t}=(t_{0},t_{1},\ldots,t_{N})$. Let $k_{n}=t_{n}-t_{n-1}$ denote the length of the $n$th (open) subinterval $I_{n}=(t_{n-1},t_{n})$. We form the disjoint union

$$I=I_{1}\cup I_{2}\cup\cdots\cup I_{N},$$

and for any function $v:I\to\mathbb{R}$ write

$$v^{n}_{+}=\lim_{\epsilon\downarrow 0}v(t+\epsilon),\qquad v^{n}_{-}=\lim_{\epsilon\downarrow 0}v(t-\epsilon),\qquad\llbracket v\rrbracket^{n}=v^{n}_{+}-v^{n}_{-},$$

provided the one-sided limits exist.

Given a vector $\boldsymbol{r}=(r_{1},r_{2},\ldots,r_{N})$ of integers $r_{n}\geq 0$, the trial space $\mathcal{X}=\mathcal{X}(\boldsymbol{t},\boldsymbol{r})$ consists of the functions $X:I\to\mathbb{R}$ such that $X|_{I_{n}}\in\mathbb{P}_{r_{n}-1}$ for $1\leq n\leq N$. Here, $\mathbb{P}_{m}$ denotes the space of polynomials of degree at most $m\geq 0$, with real coefficients. The dG solution $U\in\mathcal{X}$ of (3) is then defined by [9, 12]

$$\llbracket U\rrbracket^{n-1}X^{n-1}_{+}+\int_{I_{n}}(U^{\prime}+\lambda\partial_{t}^{1-\alpha}U)X\,dt=\int_{I_{n}}fX\,dt$$

(4)

for $X\in\mathbb{P}_{r_{n}-1}$ and $1\leq n\leq N$, where, in the case $n=1$, we set $U^{0}_{-}=u_{0}$ so that $\llbracket U\rrbracket^{0}=U^{0}_{+}-U^{0}_{-}=U^{0}_{+}-u_{0}$. (The monograph of Thomée [15, Chapter 12] is a standard reference providing a general introduction to dG time stepping for classical diffusion problems.)

To compute $U$, we choose for each $n$ a basis $\psi_{n1}$, $\psi_{n2}$, …, $\psi_{nr_{n}}$ for $\mathbb{P}_{r_{n}-1}$ and write

$$U(t)=\sum_{j=1}^{r_{n}}U^{nj}\psi_{nj}(t)\quad\text{for $t\in I_{n}$.}$$

(5)

When $X=\psi_{ni}$, we find that

$$U^{n-1}_{+}X^{n-1}_{+}+\int_{I_{n}}U^{\prime}X\,dt=\sum_{j=1}^{r_{n}}G^{n}_{ij}U^{nj}\quad\text{and}\quad U^{n-1}_{-}X^{n-1}_{+}=\sum_{j=1}^{r_{n-1}}K^{n,n-1}_{ij}U^{n-1,j},$$

with coefficients given by

$$G^{n}_{ij}=\psi_{nj}(t_{n-1})\psi_{ni}(t_{n-1})+\int_{I_{n}}\psi_{nj}^{\prime}\psi_{ni}\,dt$$

and

$$K^{n,n-1}_{ij}=\psi_{n-1,j}(t_{n-1})\psi_{ni}(t_{n-1}).$$

Owing to the convolutional structure of the fractional derivative, it is convenient to introduce the notation

$$\bar{\ell}=n-\ell$$

and define, if $t\in I_{n}$,

$$\rho^{n\bar{\ell}}_{j}(t)=\rho^{n,n-\ell}_{j}(t)=\int_{I_{\ell}}\omega_{\alpha}(t-s)\psi_{\ell j}(s)\,ds\quad\text{for $1\leq\ell\leq n-1$,}$$

with

$$\rho^{n\bar{n}}_{j}(t)=\rho^{n0}_{j}(t)=\int_{t_{n-1}}^{t}\omega_{\alpha}(t-s)\psi_{nj}(s)\,ds.$$

We find that

$$\partial_{t}^{1-\alpha}U=\sum_{\ell=1}^{n}\sum_{j=1}^{r_{\ell}}U^{\ell j}(\rho^{n\bar{\ell}}_{j})^{\prime}(t)\quad\text{for $t\in I_{n}$,}$$

and thus

$$\int_{I_{n}}(\partial_{t}^{1-\alpha}U)X\,dt=\sum_{\ell=1}^{n}\sum_{j=1}^{r_{\ell}}H^{n\bar{\ell}}_{ij}U^{\ell j}\quad\text{where}\quad H^{n\bar{\ell}}_{ij}=H^{n,n-\ell}_{ij}=\int_{I_{n}}(\rho^{n\ell}_{j})^{\prime}\psi_{ni}\,dt.$$

Hence, putting

$$F^{ni}=\int_{I_{n}}f\psi_{ni}\,dt,$$

the dG method (4) requires

$$\sum_{j=1}^{r_{n}}\bigl{(}G^{n}_{ij}+\lambda H^{n0}_{ij}\bigr{)}U^{nj}=F^{ni}-\sum_{\ell=1}^{n-1}\sum_{j=1}^{r_{\ell}}\lambda H^{n,n-\ell}_{ij}U^{\ell j}\\ +\begin{cases}\psi_{1i}(0)u_{0},&n=1,\\ \sum_{j=1}^{r_{n-1}}K^{n,n-1}_{ij}U^{n-1,j},&2\leq n\leq N.\end{cases}$$

(6)

At the $n$th time step, this $r_{n}\times r_{n}$ linear system must be solved to determine $U^{n1}$, $U^{n2}$,…, $U^{nr_{n}}$ and hence $U(t)$ for $t\in I_{n}$.

To compute $G^{n}_{ij}$, $H^{n\ell}_{ij}$ and $K^{n,n-1}_{ij}$ it is convenient to map each closed subinterval $\bar{I}_{n}=[t_{n-1},t_{n}]$ to the reference element $[-1,1]$. We therefore define the affine function $\mathsf{t}_{n}:[-1,1]\to\bar{I}_{n}$ by

$$\mathsf{t}_{n}(\tau)=\tfrac{1}{2}\bigl{[}(1-\tau)t_{n-1}+(1+\tau)t_{n}\bigr{]}\quad\text{for $-1\leq\tau\leq 1$,}$$

and let

$$\Psi_{nj}(\tau)=\psi_{nj}(t)\quad\text{for $t=\mathsf{t}_{n}(\tau)$ and $-1\leq\tau\leq 1$.}$$

In this way,

$$G^{n}_{ij}=\Psi_{nj}(-1)\Psi_{ni}(-1)+\int_{-1}^{1}\Psi_{nj}^{\prime}(\tau)\Psi_{ni}(\tau)\,d\tau$$

(9)

and

$$K^{n,n-1}_{ij}=\Psi_{n-1,j}(+1)\Psi_{ni}(-1).$$

(10)

Both of these coefficients are readily computed; the remainder of this section is devoted to $H^{n\bar{\ell}}_{ij}$. The formulae in the next lemma allow us to compute $H^{n0}_{ij}$ to machine precision via Gauss–Legendre and Gauss–Jacobi quadrature.