Using hierarchical matrices in the solution of the time-fractional heat equation by multigrid waveform relaxation

TL;DR

This paper introduces an efficient multigrid waveform relaxation method utilizing hierarchical matrices for solving the time-fractional heat equation on non-uniform grids, achieving optimal complexity and good convergence.

Contribution

It proposes a novel space-time multigrid approach with hierarchical matrix approximation to efficiently solve non-uniform grid discretizations of the time-fractional heat equation.

Findings

Method achieves ${ m O}(k N M \log(M))$ complexity.

Numerical experiments confirm efficiency and convergence.

Applicable to 1D and 2D problems.

Abstract

This work deals with the efficient numerical solution of the time-fractional heat equation discretized on non-uniform temporal meshes. Non-uniform grids are essential to capture the singularities of "typical" solutions of time-fractional problems. We propose an efficient space-time multigrid method based on the waveform relaxation technique, which accounts for the nonlocal character of the fractional differential operator. To maintain an optimal complexity, which can be obtained for the case of uniform grids, we approximate the coefficient matrix corresponding to the temporal discretization by its hierarchical matrix (${\cal H}$-matrix) representation. In particular, the proposed method has a computational cost of ${\cal O}(k N M \log(M))$, where $M$ is the number of time steps, $N$ is the number of spatial grid points, and $k$ is a parameter which controls the accuracy of the ${\cal H}$-matrix approximation. The efficiency and the good convergence of the algorithm, which can be theoretically justified by a semi-algebraic mode analysis, are demonstrated through numerical experiments in both one- and two-dimensional spaces.

Figures

Tables

Table 1. Table 1: Maximum errors E M = ‖ u − u N , M ‖ subscript 𝐸 𝑀 norm 𝑢 subscript 𝑢 𝑁 𝑀 E_{M}=\|u-u_{N,M}\| obtained for three different values of δ 𝛿 \delta , with N = 1024 𝑁 1024 N=1024 spatial grid-points and different numbers of time steps M 𝑀 M , and corresponding reduction orders log 2 ⁡ ( E M / E 2 ​ M ) subscript 2 subscript 𝐸 𝑀 subscript 𝐸 2 𝑀 \log_{2}(E_{M}/E_{2M}) .

$\delta$		$M=32$	$M=64$	$M=128$	$M=256$	$M=512$
$0.4$	$E_M=\Vert u-u_{N,M}\Vert$	$1.9E-3$	$7.0E-4$	$2.4E-4$	$8.5E-5$	$2.9E-5$
	${\log }_2\left(E_M/E_{2M}\right)$	$1.44$	$1.50$	$1.53$	$1.55$
$0.6$	$E_M=\Vert u-u_{N,M}\Vert$	$3.3E-3$	$1.4E-3$	$5.5E-4$	$2.1E-4$	$8.3E-5$
	${\log }_2\left(E_M/E_{2M}\right)$	$1.23$	$1.35$	$1.35$	$1.36$
$0.8$	$E_M=\Vert u-u_{N,M}\Vert$	$5.0E-3$	$2.4E-3$	$1.1E-3$	$5.0E-4$	$2.2E-4$
	${\log }_2\left(E_M/E_{2M}\right)$	$1.05$	$1.12$	$1.13$	$1.14$

Table 2. Table 2: Comparison between SAMA predictions and experimentally computed convergence factors (in brackets) with one smoothing step, for four values of δ 𝛿 \delta and different space-time grid sizes.

$\delta$	$64\times 64$	$128\times 128$	$256\times 256$
$0.2$	0.118 (0.118)	0.114 (0.116)	0.114 (0.114)
$0.4$	0.123 (0.124)	0.122 (0.123)	0.115 (0.117)
$0.6$	0.106 (0.107)	0.088 (0.089)	0.067 (0.068)
$0.8$	0.066 (0.067)	0.043 (0.045)	0.026 (0.028)

Table 3. Table 3: Number of V ​ ( 0 , 1 ) − limit-from 𝑉 0 1 V(0,1)- WRMG iterations necessary to reduce the initial residual in a factor of 10 − 10 superscript 10 10 10^{-10} for different fractional orders δ 𝛿 \delta and for different grid-sizes. The corresponding average convergence factors (in brackets) are also included.

$\delta$	$128\times 128$	$256\times 256$	$512\times 512$	$1024\times 1024$	$2048\times 2048$
0.2	11 (0.11)	11 (0.11)	11 (0.11)	11 (0.11)	11 (0.11)
0.4	11 (0.11)	11 (0.11)	11 (0.11)	10 (0.09)	9 (0.08)
0.6	10 (0.09)	9 (0.06)	8 (0.05)	7 (0.04)	7 (0.04)
0.8	8 (0.04)	7 (0.04)	7 (0.04)	7 (0.04)	7 (0.04)

Table 4. Table 4: Number of V ​ ( 0 , 1 ) − limit-from 𝑉 0 1 V(0,1)- WRMG iterations necessary to reduce the initial residual in a factor of 10 − 10 superscript 10 10 10^{-10} for different fractional orders δ 𝛿 \delta and for different values of k 𝑘 k (fixed grid-size 512 × 512 512 512 512\times 512 ). The corresponding average convergence factors (in brackets) are also included.

$\delta$	$k=5$	$k=10$	$k=15$	$k=20$	$k=25$	$k=30$
0.2	11 (0.11)	11 (0.11)	11 (0.11)	11 (0.11)	11 (0.11)	11 (0.11)
0.4	11 (0.11)	11 (0.11)	11 (0.11)	11 (0.11)	10 (0.11)	11 (0.11)
0.6	8 (0.05)	8 (0.05)	8 (0.05)	8 (0.05)	8 (0.05)	8 (0.05)
0.8	7 (0.04)	7 (0.04)	7 (0.04)	7 (0.04)	7 (0.04)	7 (0.04)

Table 5. Table 5: Number of V ​ ( 1 , 1 ) 𝑉 1 1 V(1,1) -WRMG iterations necessary to reduce the initial residual in a factor of 10 − 10 superscript 10 10 10^{-10} , together with the corresponding average convergence factors (in brackets), for different fractional orders δ 𝛿 \delta and for different grid-sizes.

$\delta$	$32\times 32\times 32$	$64\times 64\times 64$	$128\times 128\times 128$	$256\times 256\times 256$
0.2	9 (0.07)	9 (0.07)	9 (0.07)	9 (0.07)
0.4	9 (0.07)	9 (0.08)	9 (0.08)	9 (0.08)
0.6	9 (0.07)	9 (0.08)	9 (0.08)	9 (0.08)
0.8	9 (0.07)	9 (0.08)	9 (0.08)	9 (0.08)

Table 6. Table 6: Number of V ​ ( 1 , 1 ) 𝑉 1 1 V(1,1) -WRMG iterations necessary to reduce the initial residual in a factor of 10 − 10 superscript 10 10 10^{-10} , together with the corresponding average convergence factors (in brackets), for different fractional orders δ 𝛿 \delta and for different rank k 𝑘 k (fixed grid-size 64 × 64 × 64 64 64 64 64\times 64\times 64 ).

$\delta$	$k=5$	$k=10$	$k=15$	$k=20$	$k=25$	$k=30$
0.2	9 (0.07)	9 (0.07)	9 (0.07)	9 (0.07)	9 (0.07)	9 (0.07)
0.4	9 (0.07)	9 (0.08)	9 (0.07)	9 (0.08)	9 (0.08)	9 (0.07)
0.6	9 (0.08)	9 (0.08)	9 (0.08)	9 (0.08)	9 (0.08)	9 (0.08)
0.8	9 (0.08)	9 (0.08)	9 (0.08)	9 (0.08)	9 (0.08)	9 (0.08)

Equations

$$D_t^{\delta }u-\frac{{\partial }^{2}u}{\partial x^2}=f(x,t),0<x<L,t>0,$$

$$u(0,t)=0,u(L,t)=0,t>0,$$

$$u(x,0)=g(x),0\le x\le L.$$

$$D_t^{\delta }u(x,t):=\left[J^{1-\delta }\left(\frac{\partial u}{\partial t}\right)\right](x,t),0\le x\le L,t>0,$$

$$\left(J^{1-\delta }u\right)(x,t):=\left[\frac{1}{\Gamma (1-\delta )}{\int }_0^t(t-s{)}^{-\delta }u(x,s)ds\right],0\le x\le L,t>0,$$

$$G_h=\left\{x_n=nh,n=0,1,\dots ,N+1\right\},$$

$$D_t^{\delta }{u}_h(t)+A_h{u}_h(t)=f_h(t),u_h(0)=g_h,t>0,$$

$$D_t^{\delta }u(t)=f(t),t\in (0,T],$$

$${\int }_0^T{D}_t^{\delta }u(t){\psi }_m(t)\mathrm{d}t={\int }_0^{T}f(t){\psi }_m(t)\mathrm{d}t.$$

$$\sum _{j=1}^M\left(\frac{1}{\Gamma (1-\delta )}{\int }_0^T\left({\int }_0^t(t-s{)}^{-\delta }{\phi }_j^{\prime }(s)\mathrm{d}s\right){\psi }_m(t)\mathrm{d}t\right)u_j=f(t_m).$$

$$R_{m,j}=\frac{1}{\Gamma (1-\delta )}{\int }_0^T\left({\int }_0^t(t-s{)}^{-\delta }{\phi }_j^{\prime }(s)\mathrm{d}s\right){\psi }_m(t)\mathrm{d}t.$$

$$R_{m,j}=\frac{1}{\Gamma (1-\delta )}{\int }_0^{t_m}(t_m-s{)}^{-\delta }{\phi }_j^{\prime }(s)\mathrm{d}s=0,$$

$$d_{m,0}$$

$$d_{m,1}$$

$$d_{m,k}$$

$$D_M^{\delta }{u}_m=\frac{1}{\Gamma (2-\delta )}\left(d_{m,1}{u}_m+\sum _{k=1}^{m-1}(d_{m,k+1}-d_{m,k})u_{m-k}\right),$$

$$D_M^{\delta }{u}_{n,m}-\frac{u_{n+1,m}-2u_{n,m}+u_{n-1,m}}{h^2}=f(x_n,t_m),1\le n\le N,1\le m\le M,$$

$$u_{0,m}=0,u_{N+1,m}=0,0<m\le M,$$

$$u_{n,0}=g(x_n),0\le n\le N+1.$$

$$\tilde{\mathcal{K}}(t,s):=\sum _{\nu =0}^{k-1}{p}_{\nu }(t)q_{\nu }(s).$$

$$\mathrm{diam}(I_t)\le \mathrm{dist}(I_s,I_t).$$

$${\partial }_t^{\nu }\left[(t-s{)}^{-\delta }\right]=(-1{)}^{\nu }\prod _{l=1}^{\nu }(\delta +l-1)(t-s{)}^{-(\delta +\nu )},\nu =1,2,\dots ,$$

$$\mathcal{K}(t,s)\approx \sum _{\nu =0}^{k-1}\frac{1}{\nu !}\left(\prod _{l=1}^{\nu }(1-\delta -l)\right)(t_0-s{)}^{-(\delta +\nu )}(t-t_0{)}^{\nu }=:\tilde{\mathcal{K}}(t,s).$$

$$p_{\nu }(t)$$

$$q_{\nu }(s)$$

$$R_{m,j}\approx {\tilde{R}}_{m,j}$$

$${\tilde{R}}_{m,j}=\frac{1}{\Gamma (1-\delta )}\sum _{\nu =0}^{k-1}\left({\int }_0^T{p}_{\nu }(t){\psi }_m(t)\mathrm{d}t\right)\left({\int }_0^t{q}_{\nu }(s){\phi }_j^{\prime }(s)\mathrm{d}s\right).$$

$$R{|}_{\vartheta \times \sigma }\approx \tilde{R}{|}_{\vartheta \times \sigma }=A{B}^T,A\in {\mathbb{R}}^{|\vartheta |\times k},B\in {\mathbb{R}}^{|\sigma |\times k},$$

$$A_{m,\nu }:=\frac{1}{\Gamma (1-\delta )}{\int }_0^T{p}_{\nu }(t){\psi }_m(t)\mathrm{d}t=\frac{1}{\Gamma (1-\delta )}(t_m-t_0{)}^{\nu },$$

$$B_{j,\nu }:={\int }_0^t{q}_{\nu }(s){\phi }_j^{\prime }(s)\mathrm{d}s=\frac{c_{\nu }}{1-\delta -\nu }\left(\frac{1}{{\tau }_{j-1}}(m_{j-1}-m_j)-\frac{1}{{\tau }_j}(m_j-m_{j+1})\right),$$

$$|R_{m,j}-{\tilde{R}}_{m,j}|=\mathcal{O}(3^{-k}),$$

$$|A_{m,\nu }|\le \frac{1}{\Gamma (1-\delta )}|a-t_0{|}^{\nu },\nu \ge 0,$$

$$(t_0-t_{j-1}{)}^{1-\delta -\nu }$$

$$(t_0-t_{j+1}{)}^{1-\delta -\nu }$$

$$B_{j,\nu }=\frac{c_{\nu }(-\delta -\nu )}{2}\left({\tau }_{j-1}(t_0-{\xi }_j{)}^{-(1+\delta +\nu )}+{\tau }_j(t_0-{\xi }_{j+1}{)}^{-(1+\delta +\nu )}\right).$$

Full text

††journal: Journal of Computational Physics

1 Introduction

The design of efficient numerical methods for differential equations involving fractional derivatives has become a challenging topic recently, due to its wide range of applications in many different fields [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. The nonlocal character of the fractional differential operator usually results in a dense coefficient matrix when numerical methods are applied, leading to higher memory requirements as well as a considerable increase of the solution time comparing with solving integer differential equations. Usually, if uniform meshes are employed, the coefficient matrix has a Toeplitz-like structure, and efficient solvers have been proposed to reduce the computational complexity of traditional Gaussian elimination type methods. For time-fractional differential equations, alternating direction implicit schemes (ADI) with a computational complexity of ${\cal O}(NM^{2})$, where $N$ is the number of spatial grid-points and $M$ the number of time steps, were proposed in [12]. Also an approximate inversion method with a computational cost of ${\cal O}(NM\log(M))$ has been proposed in [13], where the authors approximate the coefficient matrix by a block $\varepsilon$-circulant matrix, which can be block diagonalized by FFT and, in order to solve the resulting complex block system, the authors use a multigrid method. A parallel-in-time method based on the parareal algorithm [14] has been proposed for solving time-fractional differential equations [15]. This method consists of an iterative predictor-corrector procedure combining an inexpensive but inaccurate solver with an expensive but accurate solver. Also, an efficient algorithm for the evaluation of the Caputo fractional derivative, based on the sum-of-exponentials technique, is presented in [16] for the solution of fractional diffusion equations. Recently, an efficient, robust and parallel-in-time multigrid method based on the waveform relaxation approach has been proposed in [17]. By exploiting the Toeplitz-like structure of the coefficient matrix, the computational complexity of the method is ${\cal O}(NM\log(M))$ with a storage requirement of ${\cal O}(NM)$. Although some works have studied the solution of fractional diffusion equations on non-uniform temporal meshes, see for example [18], most existing efficient solvers can only be applied when a uniform grid is used, or they can only achieve their best performance when that is the case. The latter happens for the fast solver presented in [17], which can be applied for nonuniform meshes but with a significant increase of the computational complexity, due to the loss of the Toeplitz structure of the matrix. In the case of space-fractional PDEs, recently, a fast solver based on a geometric multigrid method for nonuniform grids has been proposed in [19]. The key in this work is to use hierarchical matrices to approximate the dense stiffness matrices. Our aim here is to combine the two approaches proposed in [17, 19] and efficiently solve, both in terms of CPU time and memory requirements, the time-fractional heat equation on nonuniform grids.

${\cal H}$-matrices [20, 21] consist of powerful data-sparse approximations of dense matrices, providing a significant reduction of the storage requirement from ${\cal O}(n^{2})$ to ${\cal O}(nk\log(n))$ units of storage ($n$ is the matrix size), where $k$ is a parameter that controls the accuracy of the approximation. Moreover, the matrix-vector multiplication in ${\cal H}$-matrix format can be done in ${\cal O}(kn\log(n))$ operations. Therefore, on non-uniform grids, the approximation of the dense matrices arising from the time discretization of the fractional partial differential equations (PDEs) by the ${\cal H}$-matrices representation is the key for maintaining an optimal computational complexity of the multigrid waveform relaxation method [22, 23, 24].

As in standard multigrid methods (see [25, 26, 27]), the multigrid waveform relaxation method accelerates the convergence of the waveform relaxation method, which is a continuous-in-time iterative algorithm for solving large systems of ordinary differential equations (ODEs), by introducing a hierarchy of coarser levels. This method combines the very fast multigrid convergence with the high parallel efficiency of the waveform relaxation. In practice, it uses a red-black zebra-in-time line relaxation together with a coarse-grid correction procedure based on coarsening only in the spatial dimension. In this work, we develop an efficient and robust multigrid waveform relaxation method based on the ${\cal H}$-matrix representation of the discretization of the time-fractional heat equation on a nonuniform temporal grid. The good convergence properties of the algorithm will be theoretically justified by applying a semi-algebraic mode analysis (SAMA) [28]. This analysis is essentially a generalization of the classical local Fourier analysis (LFA) or local mode analysis [29, 30, 26, 27, 31] and combines the standard LFA with an algebraic computation that accounts for the non-local character of the fractional differential operators. In particular, the semi-algebraic analysis developed in [17] for the time-fractional heat equation is used in this work. In addition, the classical disadvantage of waveform methods regarding the requirement of extra storage for unknowns is not a drawback here anymore since the storage of the solutions in previous steps is needed anyway for the time-fractional PDEs.

The remainder of this work is structured as follows. In Section 2, we introduce the time-fractional model problem and its discretization in the general framework of a non-uniform temporal grid. Section 3 is devoted to present the hierarchical matrix representation of the dense matrix corresponding to the time-discretization. In Section 4, the multigrid waveform relaxation method is described, and its computational cost is estimated based on the computations in ${\cal H}$-matrix framework. In order to illustrate the good behavior of the multigrid waveform relaxation method for solving the time-fractional diffusion problem, in Section 5, numerical experiments in both one- and two-dimensional spaces are considered. In addition, some results of the semi-algebraic mode analysis are also presented in this section to theoretically confirm the good convergence results obtained numerically. Finally, some conclusions are drawn in Section 6.

2 Model problem and discretization

We consider the time-fractional heat equation, arising by replacing the first-order time derivative with the Caputo derivative of order $\delta$, where $0<\delta<1$. In the literature, this model is also known as fractional sub-diffusion equation, which is a subclass of anomalous diffusive problems [32, 33]. In this section, we restrict ourselves to the one-dimensional case for the sake of simplicity and formulate the model problem as the following initial-boundary value problem,

[TABLE]

Here $D_{t}^{\delta}$ denotes the Caputo fractional derivative [34, 35], defined as follows

[TABLE]

where $J^{1-\delta}$ is the Riemann-Liouville fractional integral operator given by

[TABLE]

with $\Gamma$ being the Gamma function [36].

In order to discretize model problem (1)-(3) we consider a uniform mesh in space

[TABLE]

where $h=\displaystyle\frac{L}{N+1}$ and $N+1$ is the number of subdivisions in the spatial domain, and a non-uniform grid in time, $G_{\tau}$ given by $0=t_{0}<t_{1}<\cdots<t_{M-1}<t_{M}=T$ with $T$ being the final time, $M$ representing the number of subdivisions for temporal discretization and the time step size is $\tau_{m}=t_{m+1}-t_{m},\ m=1,\ldots,M-1$. Then, the whole grid is given by $G_{h,\tau}=G_{h}\times G_{\tau}$. For the sake of simplicity, we use a uniform spatial grid in the presentation. Notice, however, that the proposed method can be straightforwardly applied on non-uniform grids.

The diffusion term in (1) is approximated by standard spatial discretization schemes such as finite difference or finite element methods. Here, we use a standard second order finite difference scheme, yielding the following semi-discrete problem

[TABLE]

where $u_{h}$ and $f_{h}$ are functions at time $t$ defined on the discrete spatial mesh $G_{h}$, and $A_{h}$ is the discrete space approximation. For the time discretization, we use a Petrov-Galerkin approach. In order to describe such an approximation, we consider the following simplified problem,

[TABLE]

with initial condition $u(0)=0$. As standard in the Galerkin finite element framework, we consider the finite dimensional space ${\mathcal{V}}:=\hbox{span}\{\varphi_{1},\ldots,\varphi_{M}\}$, where $\varphi_{i}$ are the standard piecewise linear basis functions defined on $G_{\tau}$. For the test functions we consider the following Dirac’s delta functions $\psi_{m}(t)=\delta(t-t_{m})$. Next, equation (5) is multiplied by $\psi_{m}(t)$ and integrated over interval $(0,T)$, which gives

[TABLE]

Since $u(t)\approx u^{M}(t)=\sum_{j=1}^{M}u_{j}\varphi_{j}(t)$, by using the definition of the Caputo fractional derivative, we obtain that

[TABLE]

This leads to a linear system of equations $R\bm{u}=\bm{f}$, where $\bm{u}=(u_{1},u_{2},\cdots,u_{M})^{T}$, $\bm{f}=(f(t_{1}),f(t_{2}),\cdots,f(t_{M}))^{T}$, and the entries of the coefficient matrix are

[TABLE]

It is easy to see that when $j>m$

[TABLE]

and, therefore, $R$ is a dense lower triangular matrix whose entries are given by $R_{m,m-k}=d_{m,k+1}-d_{m,k}$ for $k=0,\ldots,m-1$, where

[TABLE]

Note that the obtained discretization of the Caputo derivative is

[TABLE]

which corresponds to the generalization of the well-known L1 scheme [37, 38, 39] for non-uniform grids [35].

Remark 2.1**.**

By considering different choices of the test or trial functions, we can obtain different discretizations for the time-fractional Caputo derivative. For example, by using piecewise quadratic basis functions for trial functions, we can obtain a temporal discretization with higher accuracy. There are many high order schemes for the time fractional problems, see e.g., [40, 41, 42], and it will be interesting to see if we can reconstruct some of the existing high order schemes by choosing different test and trial functions. In this work, however, we focus on the L1 scheme and the study of high order schemes will be the subject of our future work.

Finally, by denoting as $u_{n,m}$ the nodal approximation to the solution at each grid point $(x_{n},t_{m})$, we approximate (1)-(3) by the following discrete problem,

[TABLE]

In the case of a temporal uniform grid with time-step $\tau$, it is proved rigorously that this scheme has a rate of convergence of ${\cal O}(h^{2}+\tau^{\delta})$ for “typical” solutions of the time-fractional heat equation, that is, solutions presenting a boundary layer at the initial time. To improve the poor convergence when the fractional order $\delta$ is very small, in [35] the authors proposed to use a graded mesh in time for which the resulting scheme has a convergence order of ${\cal O}(h^{2}+M^{-(2-\delta)})$. Such a grid is defined by $t_{m}=T(m/M)^{r}$ with $m=0,\ldots,M$, where $r\geq(2-\delta)/\delta$. In particular, in this work we choose $r=(2-\delta)/\delta$, as suggested in [35].

Notice that matrix $R$ is a dense lower triangular matrix. For the case of a uniform grid, it has Toeplitz structure, which can be exploited to develop efficient algorithms. For example, a geometric multigrid (GMG) method with a computational cost of ${\cal O}(NM\log(M))$ was proposed in [17]. For the case of a non-uniform grid, $R$ does not have Toeplitz structure anymore and although the GMG method [17] can still be applied, this gives rise to a significant increase of the computational complexity. Our aim is to approximate the matrix $R$ by $\widetilde{R}$ which can be stored in a data-sparse format and then design an efficient multigrid solver with a computational cost of ${\cal O}(kNM\log(M))$ where $k$ is a parameter that controls the accuracy of such approximation.

3 Discretization based on ${\cal H}$-Matrices representation

In this section we aim to approximate matrix $R$ by $\widetilde{R}$ based on the ${\cal H}$-matrices framework. According to (6), the entries of $R$ are defined using the kernel ${\cal K}(t,s)=(t-s)^{-\delta}$. Similar to typical kernel functions, singularities only occur when $t=s$. The kernel function is smooth everywhere else and decays when $|t-s|\rightarrow\infty$. This implies that the entries of $R$ decay to [math] when they are far away from the diagonal, and then they usually can be replaced by low-rank approximations. Such approximation can be done by replacing the original kernel by a degenerate (separable) kernel,

[TABLE]

This is a partial sum of $k$ terms which usually is obtained by truncation of certain infinite sum of $\mathcal{K}$. Moreover, each term is a product of two functions, one only depends on $t$ and the other one only depends on $s$. This approximation is appropriate, however, only when the truncation error can be bounded uniformly, i.e., $t$ and $s$ should be sufficiently far away from each other. More concretely, the feasibility of such approximation is characterized by the following condition.

Definition 3.1** (Admissibility condition).**

Let $I_{t}:=[a,b],\ I_{s}:=[c,d]\subset[0,T]$ be two intervals such that $d<a$. We say that $I_{t}\times I_{s}$ satisfies the admissibility condition if the following holds,

[TABLE]

We define the set of indices $\vartheta\times\sigma=\{(m,j)\ \ni\ (\mathop{\mathrm{supp}}\psi_{m},\mathop{\mathrm{supp}}\varphi_{j})\subset I_{t}\times I_{s}\}$. We say that $\vartheta\times\sigma$ is admissible if $I_{t}\times I_{s}$ satisfies the admissibility condition.

Remark 3.1**.**

More general admissibility conditions could be considered as shown in [20, 21]. Here, we use this simple choice to demonstrate the idea and make the presentation easy to follow.

Next, we approximate the kernel ${\cal K}(t,s)$ by its truncated Taylor expansion in a set $I_{t}\times I_{s}$ satisfying the admissibility condition. Taking into account that

[TABLE]

the truncated Taylor expansion of the kernel ${\cal K}(t,s)$ at the midpoint of interval $I_{t}$, that is $t_{0}=\frac{a+b}{2}$, is given by

[TABLE]

In this way, we have obtained a degenerate kernel $\widetilde{\cal K}(t,s)$ in $I_{t}\times I_{s}$ (10) where

[TABLE]

By using the degenerate kernel, we approximate the matrix entries $R_{m,j}$ in (6), $\forall(m,j)\in\vartheta\times\sigma$, by the following

[TABLE]

We observe that the double integral can be separated into a product of two single integrals as follows,

[TABLE]

From this expression, submatrix $R|_{\vartheta\times\sigma}$ can be approximated by $\widetilde{R}|_{\vartheta\times\sigma}$, which is a ${\bm{R}}k$-matrix [43], i.e., (see Figure 1)

[TABLE]

where $|\vartheta|$ and $|\sigma|$ denote the cardinality of sets $\vartheta$ and $\sigma$, respectively. The entries of the matrices $A$ and $B$ are given by

[TABLE]

where $m_{i}=(t_{0}-t_{i})^{1-\delta-\nu}$ and $\displaystyle c_{\nu}=\frac{1}{\nu!}\left(\prod_{l=1}^{\nu}(1-\delta-l)\right)$.

Notice that submatrix $R|_{\vartheta\times\sigma}$ is approximated by a low-rank representation with at most rank $k$, which only needs $k(|\vartheta|+|\sigma|)$ number of elements to store. More importantly, this approximation in an admissible block gives rise to a uniformly bounded element-wise truncation error, as stated in the following theorem.

Theorem 3.1**.**

Let $I_{t}:=[a,b],\ I_{s}:=[c,d]$ be two intervals on $[0,T]$ such that $d<a$. Assume that $I_{t}\times I_{s}$ satisfies the admissibility condition (11). Then, for all set of indices $(m,j)\in\vartheta\times\sigma$, we have the following estimate of the approximation error,

[TABLE]

where $k$ is the number of terms considered in the truncated Taylor expansion of the kernel ${\cal K}(t,s)$.

Proof.

From (12) it is easy to see that

[TABLE]

where $t_{0}=(a+b)/2$. Taking into account the following Taylor expansions,

[TABLE]

from (13), we obtain that

[TABLE]

By using that $|c_{\nu}|\leq 1$ and that $\xi_{j},\xi_{j+1}\in[c,d]$, we have the following bound

[TABLE]

For each $\nu\geq 0$, by using the bounds in (14) and (15), the product of $A_{m,\nu}$ and $B_{j,\nu}$ is given as follows,

[TABLE]

Denoting $\displaystyle r=\frac{|a-t_{0}|}{|a-t_{0}|+|a-d|}$, we can bound the element-wise error in the following way,

[TABLE]

The sum of the series appearing in (16) is given by,

[TABLE]

Taking into account that $\displaystyle r=\frac{{\rm diam}(I_{t})}{{\rm diam}(I_{t})+2{\rm dist}(I_{t},I_{s})}$ and that $\displaystyle\frac{{\rm diam}(I_{t})}{{\rm dist}(I_{t},I_{s})}\leq 1$ due to the admissibility condition (11), we have

[TABLE]

Note that the dominating term is $3^{-k}$, and thus we can conclude that

[TABLE]

which gives us an upper bound of the element-wise truncation error. ∎

Remark 3.2**.**

Using the $\mathcal{H}$-matrix representation means we are solving the perturbed linear system $\widetilde{R}\widetilde{\bm{u}}=\bm{f}$, $\widetilde{\bm{u}}=(\widetilde{u}_{1},\widetilde{u}_{2},\cdots,\widetilde{u}_{M})^{T}$. By the standard perturbation theory of solving linear systems of equations and the element-wise error estimate $|R_{m,j}-\widetilde{R}_{m,j}|={\cal O}(3^{-k})$ from Theorem 3.1, we can easily get that $\|\bm{u}-\widetilde{\bm{u}}\|_{\infty}=\mathcal{O}(3^{-k})$, where $\|\cdot\|_{\infty}$ denotes the standard $\ell^{\infty}$-norm. By triangular inequality, we naturally have

[TABLE]

Therefore, if we choose $k$ large enough, the error introduced by the $\mathcal{H}$-matrix representation is negligible and the error estimate remains the same, i.e., $\mathcal{O}(\tau^{\delta})$ on a temporal uniform grid and $\mathcal{O}(M^{-(2-\delta)})$ on a graded mesh in time. For example, we choose $k=20$ in our experiments, since $3^{-20}\approx 3\times 10^{-10}$ is quite small, we can see that the convergence order of the numerical solution obtained by the $\mathcal{H}$-matrix representation remains the same numerically.

The approximation of the admissible blocks of matrix $R$ gives its ${\cal H}$-matrix representation $\widetilde{R}$, in which those admissible submatrices are stored in a ${\bm{R}}k$-matrix representation and the rest submatrices are stored in a full-matrix format. In practice, such an ${\cal H}$-matrix representation is usually constructed by using tree data structures. For the details of the construction we refer the reader to the books [20, 21]. Here, we only present an intuitive explanation of the algorithm for constructing the ${\cal H}$-matrix representation. We start from the original coefficient matrix, by dividing it into four submatrices, as shown in Figure 2 (a). The upper right block is a zero block because of the structure of $R$ and it is colored in light blue. Then, we check if the other subblocks satisfy the admissibility condition (11). At this level, none of them satisfies the admissibility condition and each of them is recursively divided into four blocks, yielding to the structure shown in Figure 2 (b). The blocks above the diagonal are zero again and we check the admissibility condition for the rest of the blocks. For this level, three blocks are admissible and they are stored as low-rank approximation (colored in purple in the picture). The rest of the blocks (colored in pink) are split again recursively. In this way, and by repeating the same checking and coloring procedure, we obtain the structure in Figure 2 (c). This process is recursively carried out until we get the resulting ${\cal H}$-matrix representation, in which the non-zero subblocks satisfying the admissibility condition are stored in a low-rank matrix representation and the rest of the blocks are stored in a full-matrix representation.

We would like to emphasize that the ${\cal H}$-matrix representation is computed only for the matrix corresponding to the time-discretization, i.e., matrix $R$, which is a dense matrix. Combining this representation $\widetilde{R}$ with the spatial discretization matrix $A_{h}$ in (4), we obtain the coefficient matrix of the fully-discrete system, denoted here by $A_{h,\tau}$.

4 Multigrid waveform relaxation method

In this section, we introduce the proposed multigrid waveform relaxation method for solving (7)-(9) based on the ${\cal H}$-matrix representation, which gives an optimal algorithm even on non-uniform grids.

Waveform relaxation methods, also known as dynamic iteration methods, are continuous-in-time iterative algorithms for numerically solving large systems of ordinary differential equations. Different from standard iterative techniques, waveform relaxation iterates functions in time instead of scalar values. They can also be applied to solve time dependent PDEs by replacing the spatial derivatives by the spatial discrete approximations, obtaining semi-discretizations of the problems. In this way, the PDEs are transformed into a large set of ordinary differential equations, which can then be solved by an iterative algorithm.

For the time-fractional heat equation, after semi-discretization, we obtain the system of ordinary differential equations (4). Next step is to solve (4) by an iterative method. In particular, in this work, we consider a red-black Gauss-Seidel iteration (denoted by $S_{h}$) which consists of a two-stage procedure, i.e., the updates are performed first on the even points and then on the grid-points with odd numbering. In addition, to accelerate the convergence of the red-black Gauss-Seidel waveform relaxation, a coarse-grid correction procedure based on a coarsening strategy only in the spatial dimension is performed. This results in the so-called linear multigrid waveform relaxation algorithm [24]. If standard inter-grid transfer operators, as full-weighting restriction and linear interpolation, are considered, the algorithm of the multigrid waveform relaxation (WRMG) is given in Algorithm 1, where $k$ denotes the iteration number and $t$ is the time variable.

In the practical implementation of Algorithm 1, a discrete-time algorithm should be used. Therefore, after discretizing in time by replacing the differential operator $D_{t}^{\delta}$ by $D_{M}^{\delta}$, the previous algorithm can be interpreted as a space-time multigrid method with coarsening only in space. In order to reduce the high computational cost of solving $D_{M}^{\delta}$, we use the ${\cal H}$-matrix representation of operator $D_{M}^{\delta}$. Therefore, the algorithm is applied to the coefficient matrix ${\cal A}_{h,\tau}$ which is obtained by the spatial discretization matrix $A_{h}$ and the ${\cal H}$-matrix representation $\widetilde{R}$ in time. Finally, the whole multigrid waveform relaxation combines a zebra-in-time line relaxation with a standard semi-coarsening strategy only in the spatial dimension as shown in Algorithm 2, where $k$ denotes the iteration number.

4.1 Computational cost of the algorithm

From Algorithm 2, it is clear that the most time-consuming parts of the multigrid waveform relaxation method are the calculation of the residual and the relaxation step. These components require matrix-vector multiplications and the solution of dense lower triangular systems, respectively. A standard implementation of these parts would give rise to a computational cost of at least ${\cal O}(NM^{2})$, whereas the remaining components of the algorithm can be performed with a computational cost proportional to the number of unknowns. As we are going to see next, we can reduce the computational cost of the algorithm to ${\cal O}(kNM\log(M))$ thanks to the use of the hierarchical matrices.

In the calculation of the defect, a matrix-vector multiplication is needed. The calculations corresponding to the spatial discretization can be performed with a computational cost of ${\cal O}(NM)$. For each spatial grid-point, however, the matrix-vector multiplication $\widetilde{R}x$ for a given vector $x$ is required. This can be carried out by using the standard matrix-vector multiplication based on the ${\cal H}$-matrix format (see [21]). We take advantage of the lower triangular structure of matrix $\widetilde{R}$ to avoid the calculations corresponding to the zero upper triangular part, and use a slightly modified matrix-vector multiplication algorithm, which is given in Algorithm 3. Basically, an extra if statement is added to handle the case of zero blocks.

It is well-known that the computational complexity for the matrix-vector multiplication in ${\cal H}$-matrix format is ${\cal O}(kM\log(M))$. Since this is the computational cost for each spatial grid-point, the whole matrix-vector product in the calculation of the residual requires ${\cal O}(kNM\log(M))$ operations.

In the relaxation step, for each spatial grid-point we require the solution of a lower triangular system of $M$ equations, i.e., $(\widetilde{R}+2/h^{2}I_{M})x=b$. This is done by applying the standard forward substitution method in the ${\cal H}$-matrix format (Algorithm 4) to $(\widetilde{R}+2/h^{2}I_{M})$. It is well-known that such a method has a computational cost of ${\cal O}(kM\log(M))$, see e.g., [20, 21].

Notice that in Algorithm 4, we only need to consider the two cases of the matrix being stored in a full-matrix format or in a $(2\times 2)$-block form because the diagonal blocks of $\widetilde{R}$ are always not admissible and stored in a dense matrix format.

Overall, the computational cost of an iteration of one $V$-cycle multigrid waveform relaxation method is ${\cal O}(kNM\log(M))$. Since the convergence rate of such an algorithm usually is independent of the number of unknowns, as we will see, only few cycles are needed to reach the desired accuracy, and the total computational cost for solving the time-fractional heat equation on graded meshes is ${\cal O}(kNM\log(M))$.

Remark 4.1**.**

All the parallelization techniques for multigrid methods can be used for the implementation of the multigrid waveform algorithm on parallel computers. Moreover, the proposed algorithm can be parallelized in the time direction by using the multigrid waveform relaxation with cyclic reduction [22].

5 Numerical results

This section is devoted to illustrate the good performance of the proposed multigrid waveform relaxation method for the solution of the time-fractional heat equation when graded meshes are considered. We also present some results obtained by a semi-algebraic mode analyisis, which usually provides accurate predictions of the performance of the multigrid methods, and indeed, it can be made rigorous if appropriate boundary conditions are considered. In particular, here we consider the semi-algebraic mode analysis introduced in [28], which was already applied to study the convergence of the WRMG algorithm for the time-fractional heat equation in the case of uniform temporal grids in [17]. This analysis is based on an exponential Fourier mode analysis or local Fourier analysis technique only in space and an exact analytical approach in time. This is the key that allows us to apply this analysis when non-uniform grids in time are considered. For the details of this analysis we refer the readers to [17]. The only modification that we need to do when using the graded meshes is to consider the new coefficients of the matrix corresponding to the time discretization. Notice that the analysis is applied to the original discretization rather than to the ${\cal H}$-matrix representation that we use in the implementation. As we proved in Theorem 3.1, the element-wise difference between both matrices is ${\cal O}(3^{-k})$ and, therefore, by choosing sufficiently large $k$, the results of the analysis by using the original matrix are reliable to predict the practical results obtained by using the ${\cal H}$-matrix representation (see Remark 3.2). We will see in the following experiments that this analysis gives rise to very accurate predictions.

We consider two numerical experiments: one test problem which is one dimensional in space, and a second example which is two dimensional in space. We choose rank $k=20$ in our experiments except in the cases where we vary $k$. We will consider $V$-cycles for all the experiments since their convergence rates are similar to those obtained by $W$-cycles and, therefore, they provide a more efficient multigrid method in practice.

All numerical computations were carried out using MATLAB on a MacBook Pro with a Core i5 2.7 GHz and 8 GB RAM, running OS X 10.10 (Yosemite).

5.1 One-dimensional time-fractional heat equation

In the first numerical experiment we consider a problem whose solution is smooth away from the initial time but has a certain singular behavior at $t=0$ where it presents an initial layer. These are reasonably general and realistic hypotheses on the behavior of the solution of the considered problems near the initial time. In particular, we consider problem (1)-(3) defined on $[0,\pi]\times[0,1]$, with a zero right-hand side ($f(x,t)=0$) and an initial condition $g(x)=\sin\,x$. The solution of this problem is $u(x,t)=E_{\delta}(-t^{\delta})\sin\,x$, where $E_{\delta}:{\mathbb{R}}\rightarrow{\mathbb{R}}$ is given by

[TABLE]

(see [44, 35]). In Figure 3 (a), we can observe the singularity of the analytical solution for near the initial time, where an initial layer appears. The picture corresponds to the fractional order $\delta=0.6$, and following the rule given in Section 2 to construct the optimal graded mesh, in Figure 3 (b) such a grid is displayed.

As stated in [35], for “typical” solutions of (1)-(3), a rate of convergence of ${\mathcal{O}}(h^{2}+\tau^{\delta})$ is obtained when the discrete scheme on uniform grids is considered, whereas the convergence order improves to ${\mathcal{O}}(h^{2}+M^{-(2-\delta)})$ with the use of graded meshes. This can be seen in Figure 4, where for a fractional order $\delta=0.4$ and both uniform and graded meshes, we display the maximum errors between the analytical and the numerical solution for various numbers of time-steps $M$ and by using a sufficiently fine spatial grid ($N=2048$). It is observed that the convergence order with the uniform grid is $0.4$, whereas it increases to $1.6$ when the graded mesh is used.

For different values of $\delta$, in Table 1 we display the maximum errors $E_{M}=\|u-u_{N,M}\|$ obtained when using $N=1024$ and different values of $M$, together with the corresponding reduction orders computed by $\log_{2}(E_{M}/E_{2M})$. It can be seen that the reduction orders asymptotically match with the expected convergence rates.

Next, we study the convergence of the proposed multigrid waveform relaxation method on graded meshes. For this purpose, we perform a semi-algebraic mode analysis which provides very accurate predictions of the asymptotic convergence factors of the multigrid method. In Table 2, we show the two-grid convergence factors with one smoothing step provided by the analysis for four values of the fractional order $\delta$ and different values of $N$ and $M$. The convergence factors obtained by using a multilevel $W$-cycle in numerical experiments are displayed (in brackets) in this table as well, showing a good agreement between the predictions and the real convergence rates. This demonstrates that the analysis is a very useful tool for studying the convergence of the proposed multigrid waveform relaxation method.

Since we look for an efficient solver that is robust with respect to the number of time steps ($M$), in Figure 5, we show the two-grid convergence factors predicted by the analysis for different values of $M$ as well as different fractional orders $\delta$, for a fixed value of $N=128$. Indeed, the convergence rates of the proposed multigrid method are bounded below $0.2$ for any values of the parameters $\delta$ and $M$ which demonstrates its robustness.

In order to have a robust multigrid solver, the convergence rate of the method should also be independent of the number of spatial grid-points, $N$. We perform the semi-algebraic analysis to study the performance for different values of $N$, and a fixed number of time-steps $M=256$. The results are presented in Figure 6, where the two-grid convergence factors predicted by SAMA are shown for different values of $\delta$ and for a wide range of values of $N=2^{l}$ with $l=1,\ldots,20$.

We can observe that in all cases the obtained convergence factors are bounded below $0.2$, providing a robust solver for the considered time-fractional problem for any value of the fractional order $\delta$.

In practice, $V$-cycle multigrid schemes are often preferred comparing with $W$-cycles because of their lower computational cost. Here we use $V$-cycles in the following numerical experiments since, as we see, they provide a convergence that is not far from that of the two-grid or $W$-cycles.

We consider a $V$-cycle with no pre-smoothing and only one post-smoothing step. We choose a $V(0,1)$ instead of a $V(1,0)-$cycle since the first approach provides better convergence factors. In Table 3, we display the number of WRMG iterations necessary to reduce the initial residual in a factor of $10^{-10}$ for different values of the fractional order $\delta$ and different grid-sizes varying from $128\times 128$ to $2048\times 2048$ doubling the mesh-size in both spatial and temporal dimensions. We also show the obtained average convergence factors. We can conclude that the convergence of the considered WRMG is very robust independently of the spatial discretization parameter and of the use of the graded mesh.

Next, we investigate the effect of the rank $k$ on the convergence behavior of the V-cycle multigrid. Here, we fix the grid-size to be $512\times 512$ and vary the fractional order $\delta$ and rank $k$. The results are shown in Table 4. As previously, we display the number of WRMG iterations necessary to reduce the initial residual in a factor of $10^{-10}$, together with the corresponding average convergence factors. It is clear that $k$ does not affect the performance of the V-cycle multigrid, more precisely, the number of iterations stays the same and the convergence factors also remain constant. This is expected because the rank $k$ only affects the approximation of the $\mathcal{H}$-matrix representation.

Taking into account the excellent convergence rates obtained and that, as we previously commented, the total computational cost is ${\cal O}(kNM\log(M))$ thanks to the use of the ${\cal H}$-matrix representation, we provide a very efficient solver for the time-fractional heat equation on graded meshes. In Figure 7, we show the CPU time of the proposed WRMG method for different fractional orders. We observe that, for all cases, the computational complexity is optimal, which confirms our expectation.

5.2 Two-dimensional time-fractional heat equation

The aim of this second numerical experiment is to show that the proposed strategy can be extended to problems with higher spatial dimensions. In particular, here we consider a two-dimensional time-fractional diffusion model problem defined on the spatial domain $\Omega=(0,\pi)\times(0,\pi)$ given by

[TABLE]

where the right-hand side $f$ is defined as

[TABLE]

It can be easily seen that the analytic solution is

[TABLE]

A semi-algebraic analysis is performed analogously to the one-dimensional spatial case. The only difference is that a standard two-dimensional local Fourier analysis is used now in the spatial domain, combined again with an exact analytical approach in time. First of all, in Figure 8, the asymptotic convergence factors obtained by using a multilevel $W$-cycle are compared to the two-grid convergence rates predicted by the semi-algebraic analysis. This comparison is done by choosing $\delta=0.4$, $M=128$ time-steps, and for a range of values of $N$ from $N=8\times 8$ to $N=512\times 512$. Accurate correspondence can be observed between the real and the predicted values and such comparisons are similar for different values of the fractional order, which shows that the analysis provides a very useful tool for the study of the convergence of the proposed method.

The semi-algebraic analysis is used now to demonstrate the robustness of the multigrid waveform relaxation method with respect to the fractional order $\delta$. To this end, in Figure 9, the two-grid convergence factors provided by SAMA are shown for different values of $\delta$ and different numbers of spatial grid-points $N$, for a fixed number of time-steps $M=128$. A very satisfactory convergence is observed in all cases, making the multigrid waveform relaxation method a good choice for an efficient solution of the time-fractional two-dimensional heat equation.

In order to show the efficiency and the robustness of the proposed method in the case of two spatial dimensions, we consider a multigrid $V$-cycle with one pre- and one post-smoothing steps. As for the one-dimensional case, in Table 5, we show a performance of the method independent of the grid-size and robust with respect to the value of the fractional order $\delta$. In particular, we show the number of WRMG iterations that are needed to reduce the initial maximum residual by a factor of $10^{-10}$ for different grid-sizes varying from $32\times 32\times 32$ to $256\times 256\times 256$ and for different values of the fractional order $\delta$. The average convergence factors are displayed as well.

In Table 6, we study the influence of the rank $k$ in the WRMG performance. We fix the grid-size to be $64\times 64\times 64$ and vary the fractional order $\delta$ and the rank $k$. We can see from Table 6 that $k$ does not affect the performance of the V-cycle multigrid, since the number of iterations stays the same, as well as the convergence factors remain constant. This demonstrates the robustness of our $\mathcal{H}$-matrix representation.

Again, we obtain an excellent convergence that together with the computational cost of ${\cal O}(kNM\log(M))$, makes the proposed WRMG method a very efficient solver also for the two-dimensional time-fractional heat equation on graded meshes. The CPU time of the proposed WRMG method is shown in Figure 10 for different fractional orders, where we can confirm that the computational complexity is optimal for all cases, as expected.

6 Conclusions

In this work, we have proposed a fast solver for the time-fractional heat equation, targeting in particular “typical” non-smooth solutions that arise even when right-hand sides are smooth. The use of uniform grids in these cases may lead to a very poor convergence of the numerical solution, which can be enhanced by considering graded meshes. The algorithm is based on the multigrid acceleration of the waveform relaxation method, which provides an optimal complexity only for uniform grids. For non-uniform grids in time, however, the computational cost increases and is not optimal any more. In order to overcome this drawback and keep the optimality of the overall computational complexity of the method, the hierarchical matrices framework is considered. In this way, a computational cost of ${\cal O}(kNM\log(M))$, where $M$ is the number of time steps and $N$ is the number of spatial grid points, is obtained. Numerical experiments with one- and two-spatial dimensions are presented to illustrate the efficiency of the method and its robustness with respect to the fractional orders. Within these tests, a semi-algebraic mode analysis is used to theoretically justify the good convergence rates provided by the multigrid waveform relaxation algorithm.

One direction for future work is to study how to generalize our fast solver to high order schemes. This involves rewriting the high order schemes using the finite element framework. If this is doable, we can replace the original kernel by a separable kernel and construct the corresponding $\mathcal{H}$-matrix representations. Otherwise, our method cannot be directly applied and we need to investigate other ways to construct the low rank approximations, for example, construct them algebraically. Another possible research direction is the generalization of the proposed approach to time-fractional nonlinear problems. We plan to use the $\mathcal{H}$-matrix representation to approximate the time-fractional derivatives and nonlinear multigrid schemes to handle the nonlinearity, for example, using the full approximation scheme [29].

Acknowledgements

Francisco J. Gaspar has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 705402, POROSOS. The work of Carmen Rodrigo is supported in part by the Spanish project FEDER /MCYT MTM2016-75139-R and the Diputación General de Aragón (Grupo de Referencia APEDIF, ref. E24_17R). The work of Hu is partially supported by NSF grant DMS-1620063.