Fast algorithms for convolution quadrature of Riemann-Liouville fractional derivative

TL;DR

This paper develops fast, memory-efficient algorithms for computing Riemann-Liouville fractional derivatives using convolution quadrature, enabling efficient solutions to fractional Fokker-Planck equations with improved accuracy and reduced computational costs.

Contribution

The paper introduces novel fast algorithms for fractional derivatives that do not assume solution regularity, significantly reducing computational time and memory for fractional PDEs.

Findings

Algorithms achieve first- and second-order accuracy in time.

Numerical examples confirm convergence and efficiency.

Memory and computational costs are substantially decreased.

Abstract

Recently, the numerical schemes of the Fokker-Planck equations describing anomalous diffusion with two internal states have been proposed in [Nie, Sun and Deng, arXiv: 1811.04723], which use convolution quadrature to approximate the Riemann-Liouville fractional derivative; and the schemes need huge storage and computational cost because of the non-locality of fractional derivative and the large scale of the system. This paper first provides the fast algorithms for computing the Riemann-Liouville derivative based on convolution quadrature with the generating function given by the backward Euler and second-order backward difference methods; the algorithms don't require the assumption of the regularity of the solution in time, while the computation time and the total memory requirement are greatly reduced. Then we apply the fast algorithms to solve the homogeneous fractional Fokker-Planck equations with two internal states for nonsmooth data and get the first- and second-order accuracy in time. Lastly, numerical examples are presented to verify the convergence and the effectiveness of the fast algorithms.

Tables

Table 1. Table 1. L 2 subscript 𝐿 2 L_{2} error, convergence rates, and CPU time at t = 1 𝑡 1 t=1 with h = 1 / 256 ℎ 1 256 h=1/256

$({\alpha }_1,{\alpha }_2)$		$1/\tau$	100	200	400	800	1600
		$E_{1,\tau }$	8.434E-06	4.141E-06	2.001E-06	9.335E-07	4.000E-07
			Rate	1.0263	1.0489	1.1003	1.2228
	BE	$E_{2,\tau }$	1.347E-04	6.609E-05	3.193E-05	1.489E-05	6.379E-06
			Rate	1.0276	1.0496	1.1007	1.2230
(0.3,0.6)		CPU time(s)	1.34	3.19	7.08	22.33	78.53
		$E_{1,\tau }$	8.434E-06	4.141E-06	2.001E-06	9.335E-07	4.000E-07
			Rate	1.0263	1.0489	1.1003	1.2228
	FBE	$E_{2,\tau }$	1.347E-04	6.609E-05	3.193E-05	1.489E-05	6.379E-06
			Rate	1.0277	1.0496	1.1007	1.2229
		CPU time(s)	1.39	2.92	5.78	11.97	23.95
		$E_{1,\tau }$	9.318E-06	4.574E-06	2.211E-06	1.031E-06	4.417E-07
			Rate	1.0266	1.0490	1.1004	1.2228
	BE	$E_{2,\tau }$	1.390E-04	6.811E-05	3.289E-05	1.533E-05	6.568E-06
			Rate	1.0290	1.0503	1.1010	1.2231
(0.4,0.7)		CPU time(s)	1.48	3.14	7.27	22.42	79.66
		$E_{1,\tau }$	9.318E-06	4.574E-06	2.210E-06	1.031E-06	4.417E-07
			Rate	1.0266	1.0490	1.1004	1.2228
	FBE	$E_{2,\tau }$	1.390E-04	6.811E-05	3.289E-05	1.533E-05	6.568E-06
			Rate	1.0290	1.0503	1.1010	1.2231
		CPU time(s)	1.47	2.98	5.72	12.33	24.61

Table 2. Table 2. L 2 subscript 𝐿 2 L_{2} error, convergence rates, and CPU time at t = 1 𝑡 1 t=1 with h = 1 / 256 ℎ 1 256 h=1/256

$({\alpha }_1,{\alpha }_2)$		$1/\tau$	100	200	400	800	1600
		$E_{1,\tau }$	1.916E-05	9.411E-06	4.550E-06	2.123E-06	9.095E-07
			Rate	1.0254	1.0485	1.1001	1.2227
	BE	$E_{2,\tau }$	3.467E-05	1.699E-05	8.204E-06	3.824E-06	1.638E-06
			Rate	1.0290	1.0502	1.1010	1.2231
(0.3,0.7)		CPU time(s)	1.36	3.16	7.13	22.09	78.36
		$E_{1,\tau }$	1.916E-05	9.411E-06	4.550E-06	2.123E-06	9.095E-07
			Rate	1.0254	1.0485	1.1001	1.2227
	FBE	$E_{2,\tau }$	3.467E-05	1.699E-05	8.203E-06	3.824E-06	1.638E-06
			Rate	1.0290	1.0502	1.1010	1.2231
		CPU time(s)	1.41	2.84	5.89	11.94	24.02
		$E_{1,\tau }$	2.510E-05	1.233E-05	5.959E-06	2.779E-06	1.191E-06
			Rate	1.0260	1.0487	1.1003	1.2228
	BE	$E_{2,\tau }$	3.364E-05	1.650E-05	7.971E-06	3.717E-06	1.592E-06
			Rate	1.0276	1.0496	1.1007	1.2230
(0.4,0.6)		CPU time(s)	1.41	3.08	7.16	22.17	79.25
		$E_{1,\tau }$	2.510E-05	1.233E-05	5.959E-06	2.779E-06	1.191E-06
			Rate	1.0260	1.0487	1.1002	1.2227
	FBE	$E_{2,\tau }$	3.364E-05	1.650E-05	7.971E-06	3.717E-06	1.592E-06
			Rate	1.0276	1.0496	1.1006	1.2229
		CPU time(s)	1.41	2.92	6.03	12.05	24.41

Table 3. Table 3. L 2 subscript 𝐿 2 L_{2} error and convergence rates at t = 1 𝑡 1 t=1 with h = 1 / 1024 ℎ 1 1024 h=1/1024

$({\alpha }_1,{\alpha }_2)$		$1/\tau$	20	40	80	160	320
		$E_{1,\tau }$	9.445E-06	2.274E-06	5.521E-07	1.304E-07	2.598E-08
			Rate	2.0546	2.0419	2.0825	2.3273
	SBD	$E_{2,\tau }$	3.976E-05	9.551E-06	2.317E-06	5.467E-07	1.090E-07
(0.3,0.4)			Rate	2.0577	2.0435	2.0831	2.3271
		$E_{1,\tau }$	9.463E-06	2.274E-06	5.519E-07	1.303E-07	2.611E-08
			Rate	2.0572	2.0427	2.0829	2.3190
	FSBD	$E_{2,\tau }$	3.981E-05	9.552E-06	2.316E-06	5.465E-07	1.095E-07
			Rate	2.0593	2.0442	2.0831	2.3188
		$E_{1,\tau }$	3.189E-05	7.594E-06	1.835E-06	4.343E-07	8.803E-08
			Rate	2.0703	2.0487	2.0792	2.3028
	SBD	$E_{2,\tau }$	9.081E-05	2.158E-05	5.225E-06	1.255E-06	2.675E-07
(0.7,0.8)			Rate	2.0733	2.0461	2.0581	2.2296
		$E_{1,\tau }$	3.193E-05	7.594E-06	1.836E-06	4.349E-07	8.947E-08
			Rate	2.0720	2.0487	2.0774	2.2812
	FSBD	$E_{2,\tau }$	9.089E-05	2.158E-05	5.228E-06	1.263E-06	2.863E-07
			Rate	2.0744	2.0454	2.0499	2.1406

Table 4. Table 4. L 2 subscript 𝐿 2 L_{2} error and convergence rates at t = 1 𝑡 1 t=1 with h = 1 / 1024 ℎ 1 1024 h=1/1024

$({\alpha }_1,{\alpha }_2)$		$1/\tau$	20	40	80	160	320
		$E_{1,\tau }$	3.262E-06	7.860E-07	1.910E-07	4.511E-08	8.992E-09
			Rate	2.0534	2.0412	2.0817	2.3269
	SBD	$E_{2,\tau }$	5.178E-06	1.246E-06	3.024E-07	7.139E-08	1.422E-08
(0.2,0.4)			Rate	2.0555	2.0424	2.0827	2.3275
		$E_{1,\tau }$	3.228E-06	7.859E-07	1.909E-07	4.507E-08	9.028E-09
			Rate	2.0382	2.0415	2.0825	2.3197
	FSBD	$E_{2,\tau }$	5.152E-06	1.246E-06	3.023E-07	7.134E-08	1.429E-08
			Rate	2.0482	2.0429	2.0832	2.3193
		$E_{1,\tau }$	1.141E-05	2.729E-06	6.607E-07	1.561E-07	3.136E-08
			Rate	2.0636	2.0461	2.0814	2.3155
	SBD	$E_{2,\tau }$	1.543E-05	3.669E-06	8.887E-07	2.137E-07	4.577E-08
(0.6,0.8)			Rate	2.0727	2.0455	2.0561	2.2231
		$E_{1,\tau }$	1.142E-05	2.729E-06	6.606E-07	1.562E-07	3.160E-08
			Rate	2.0655	2.0464	2.0804	2.3053
	FSBD	$E_{2,\tau }$	1.545E-05	3.669E-06	8.892E-07	2.151E-07	4.918E-08
			Rate	2.0740	2.0448	2.0473	2.1290

Table 5. Table 5. L 2 superscript 𝐿 2 L^{2} error and CPU time at t=10

	$({\alpha }_1,{\alpha }_2)$	(0.3,0.8)	(0.4,0.7)	(0.5,0.6)
	$E_{1,\tau }$	2.628E-10	4.019E-10	7.399E-10
SBD	$E_{2,\tau }$	7.571E-09	3.471E-09	1.561E-09
	CPU time(s)	870.13	883.70	882.48
	$E_{1,\tau }$	3.543E-10	4.862E-10	1.008E-09
FSBD	$E_{2,\tau }$	1.321E-08	5.697E-09	2.375E-09
	CPU time(s)	209.44	210.83	211.36

Equations

$$\{\begin{array}{rl}& {\mathbf{M}}^T\frac{\partial }{\partial t}\mathbf{G}=({\mathbf{M}}^T-\mathbf{I})\mathrm{diag}({}_0{D}_t^{1-{\alpha }_1},{}_0{D}_t^{1-{\alpha }_2})\mathbf{G}\\ & +{\mathbf{M}}^T\mathrm{diag}({}_0{D}_t^{1-{\alpha }_1},{}_0{D}_t^{1-{\alpha }_2})\Delta \mathbf{G}\mathrm{in}\Omega ,t\in [0,T],\\ & \mathbf{G}(\cdot ,0)={\mathbf{G}}_0\mathrm{in}\Omega ,\\ & \mathbf{G}=0\mathrm{on}\partial \Omega ,t\in [0,T],\end{array}$$

$$\mathbf{M}=\left[\begin{array}{cc}m& 1-m\\ 1-m& m\end{array}\right],m\in (0,1/2)\bigcup (1/2,1);$$

$${}_0{D}_t^{1-{\alpha }_i}G=\frac{1}{\Gamma ({\alpha }_i)}\frac{\partial }{\partial t}{\int }_0^t(t-\xi {)}^{{\alpha }_i-1}G(\xi )d\xi ,{\alpha }_i\in (0,1),i=1,2.$$

$$\Vert \vartheta {\Vert }_{{\dot{H}}^r(\Omega )}^2=\sum _{j=1}^{\infty }{\lambda }_j^r(\vartheta ,{\phi }_j{)}^2,$$

$${\Sigma }_{\theta }=\{z\in \mathbb{C}\setminus \{0\},|\arg z|\le \theta \},$$

$${\Sigma }_{\theta ,\kappa }=\{z\in \mathbb{C}:|z|>\kappa ,|\arg z|\le \theta \},$$

$${\Gamma }_{\theta ,\kappa }=\{z\in \mathbb{C}:|z|=\kappa ,|\arg z|\le \theta \}\bigcup \{z\in \mathbb{C}:z=r{e}^{\pm \mathbf{i}\theta }:r\ge \kappa \},$$

$$\{\begin{array}{rl}& \frac{\partial G_1}{\partial t}+a{}_0{D}_t^{1-{\alpha }_1}{G}_1+{}_0{D}_t^{1-{\alpha }_1}A{G}_1=a{}_0{D}_t^{1-{\alpha }_2}{G}_2\mathrm{in}\Omega ,t\in [0,T],\\ & \frac{\partial G_2}{\partial t}+a{}_0{D}_t^{1-{\alpha }_2}{G}_2+{}_0{D}_t^{1-{\alpha }_2}A{G}_2=a{}_0{D}_t^{1-{\alpha }_1}{G}_1\mathrm{in}\Omega ,t\in [0,T],\\ & \mathbf{G}(\cdot ,0)={\mathbf{G}}_0\mathrm{in}\Omega ,\\ & \mathbf{G}=0\mathrm{on}\partial \Omega ,t\in [0,T],\end{array}$$

$${\tilde{G}}_1=$$

$${\tilde{G}}_2=$$

$$H(z)={\left((z^{{\alpha }_1}+a+A)(z^{{\alpha }_2}+a+A)-a^2\right)}^{-1},$$

$$H_{{\alpha }_1}(z)=H(z)(z^{{\alpha }_2}+a+A),H_{{\alpha }_2}(z)=H(z)(z^{{\alpha }_1}+a+A).$$

$$\Vert {\left(z^{{\alpha }_1}+a+A\right)}^{-1}\Vert \le C|z{|}^{-{\alpha }_1},\Vert {\left(z^{{\alpha }_2}+a+A\right)}^{-1}\Vert \le C|z{|}^{-{\alpha }_2},$$

$$\Vert H(z)\Vert \le C|z{|}^{-{\alpha }_1-{\alpha }_2},\Vert H_{{\alpha }_1}(z)\Vert \le C|z{|}^{-{\alpha }_1},\Vert H_{{\alpha }_2}(z)\Vert \le C|z{|}^{-{\alpha }_2},$$

$$\Vert AH(z)\Vert \le C\min \left(|z{|}^{-{\alpha }_1},|z{|}^{-{\alpha }_2}\right),\Vert A{H}_{{\alpha }_1}(z)\Vert \le C,\Vert A{H}_{{\alpha }_2}(z)\Vert \le C,$$

$$\Vert (a+A)H(z)\Vert \le C\min \left(|z{|}^{-{\alpha }_1},|z{|}^{-{\alpha }_2}\right),\Vert (a+A)H_{{\alpha }_1}(z)\Vert \le C,\Vert (a+A)H_{{\alpha }_2}(z)\Vert \le C.$$

$$(a+A)H(z)$$

$$=(z^{{\alpha }_2}+a+A)H(z)-z^{{\alpha }_2}H(z).$$

$$(z^{{\alpha }_2}+a+A)u-a^2(z^{{\alpha }_1}+a+A{)}^{-1}u=v,$$

$$u={\left(z^{{\alpha }_2}+a+A\right)}^{-1}v+a^2{\left((z^{{\alpha }_1}+a+A)(z^{{\alpha }_2}+a+A)\right)}^{-1}u.$$

$$\Vert u{\Vert }_{L^2(\Omega )}\le C|z{|}^{-{\alpha }_2}\Vert v{\Vert }_{L^2(\Omega )}+C{a}^2|z{|}^{-{\alpha }_1-{\alpha }_2}\Vert u{\Vert }_{L^2(\Omega )}.$$

$$\Vert (z^{{\alpha }_1}+a+A)H(z)\Vert \le C|z{|}^{-{\alpha }_2}.$$

$$\Vert (z^{{\alpha }_2}+a+A)H(z)\Vert \le C|z{|}^{-{\alpha }_1}.$$

$$\Vert z^{{\alpha }_1}H(z)\Vert \le C|z{|}^{-{\alpha }_2},\Vert z^{{\alpha }_2}H(z)\Vert \le C|z{|}^{-{\alpha }_1}.$$

$$\frac{1}{t^{\beta }}=\frac{1}{\Gamma (\beta )}{\int }_0^{\infty }{e}^{-ts}{s}^{\beta -1}ds.$$

$${}_0{D}_t^{\alpha }G(t)=$$

$$=$$

$$z^{\alpha }=\frac{1}{\Gamma (1-\alpha )\Gamma (\alpha )}{\int }_0^{\infty }{s}^{\alpha -1}\frac{z}{s+z}ds.$$

$${\left(\frac{1-\zeta }{\tau }\right)}^{\alpha }=\frac{1}{\Gamma (1-\alpha )\Gamma (\alpha )}{\int }_0^{\infty }{s}^{\alpha -1}\left(1-\frac{\tau s}{\tau s+1-\zeta }\right)ds.$$

$${\left(\frac{1-\zeta }{\tau }\right)}^{\alpha }=\sum _{i=0}^{\infty }{d}_{1,i}^{\alpha }{\zeta }^i$$

$$\frac{\tau s}{\tau s+1-\zeta }=\frac{\tau s}{1+\tau s}\sum _{i=0}^{\infty }{\left(\frac{\zeta }{1+\tau s}\right)}^i,$$

$$d_{1,0}^{\alpha }=\frac{1}{\Gamma (1-\alpha )\Gamma (\alpha )}{\int }_0^{\infty }{s}^{\alpha -1}\left(1-\frac{\tau s}{\tau s+1}\right)ds,$$

$$d_{1,i}^{\alpha }=-\frac{1}{\Gamma (1-\alpha )\Gamma (\alpha )}{\int }_0^{\infty }{s}^{\alpha -1}\frac{\tau s}{(\tau s+1{)}^{i+1}}ds,i\ge 1.$$

Full text

1. Introduction

Nowadays, it is widely recognized that anomalous diffusions are ubiquitous in the natural world, which also naturally become an interdisciplinary research topic. Microscopically, various stochastic processes are introduced, including Lévy process, Lévy walk, Lévy flights, continuous time random walks with power law waiting times and/or jump lengths, while macroscopically the diverse partial differential equations (PDEs) governing the probability density functions (PDFs) of a variety of statistical observables, e.g., position and functionals, are derived [4]. A lot of efforts are made for numerically solving these PDEs [1, 3, 5, 6, 14, 13, 15], and generally they are nonlocal, which urges people to develop the fast algorithm to deal with the challenges of huge storage and computational complexity [7, 23, 24].

More recently, the anomalous diffusions with multiple internal states are carefully studied, and the corresponding macroscopic PDEs are built [21, 22]. Then, [19] provides a numerical scheme and does the numerical analyses for the homogeneous fractional Fokker-Planck equations with two internal states [21], i.e.,

[TABLE]

where $\Omega$ denotes a bounded convex polygonal domain in $\mathbb{R}^{d}$ $(d=1,2,3)$; $\mathbf{M}$ is the transition matrix of a Markov chain, being a $2\times 2$ invertible matrix here, namely, we set

[TABLE]

$\mathbf{G}=[G_{1},G_{2}]^{T}$ denotes the solution of the system (1.1) and $\mathbf{G}_{0}=[G_{1}(0),G_{2}(0)]^{T}$ is the initial value; $\mathbf{I}$ is an identity matrix; ‘diag’ denotes a diagonal matrix formed from its vector argument, and $~{}_{0}D^{1-\alpha_{i}}_{t}$, $i=1,2$ are the Riemann-Liouville fractional derivatives defined by [20]

[TABLE]

The convolution quadrature [16, 17, 18] is a popular strategy to approximate (1.3), since it doesn’t require the assumption of regularity of the solution in time and it can achieve high order accuracy after the suitable modification [8, 9, 10, 11, 12]. The backward Euler (BE) convolution quadrature is used to solve the system (1.1) with first-order accuracy in [19]; it’s worth pointing out that the time step size $\tau$ must be very small to ensure the stability and convergence of the algorithm when $m$ is less than but close to $0.5$, which results in huge storage cost and computational complexity. So fast algorithm is expected to be developed.

In this paper, besides the first-order approximation given in [19], we also discuss the second-order approximation of (1.3) designed by convolution quadrature with generating function given by second-order backward difference (SBD). We modify the k-th order $(k=1,2)$ approximations of $~{}_{0}D^{\alpha}_{t}G(t_{n})$ based on convolution quadrature to speed up the calculations, that is, we use the sum of geometric sequences to approximate the weights generated by BE and SBD convolution quadrature. According to the property of the geometric sequences, the computation can be performed iteratively, which greatly reduce the computational complexity and the storage cost (for the details, refer to Section 3). Afterwards, we apply the designed fast algorithms to solve the system (1.1) and get the first- and second-order accuracy in time. Compared with the existing fast algorithms for fractional derivatives, the ones provided in this paper have the advantage of weakening the requirement of the regularity of the solution in time.

This paper is organized as follows. In Section 2, we introduce some needed notations and lemmas. In Section 3, we develop the fast algorithms based on convolution quadrature for Riemann-Liouville fractional derivatives, i.e., fast BE and SBD discretizations. In Section 4, we use the fast algorithms to solve the homogeneous fractional Fokker-Planck equations (1.1) with the first- and second-order accuracy, respectively, in time. Section 5 shows the effectiveness of the fast algorithms by numerical experiments.

2. Preliminaries

Let’s begin with some needed notations. Throughout the paper, $C$ denotes a generic positive constant, whose value may differ at each occurrence. We denote $G_{1}(t)$, $G_{2}(t)$ as the functions $G_{1}(\cdot,t)$, $G_{2}(\cdot,t)$ respectively, introduce $\|\cdot\|$ as the operator norm from $L^{2}(\Omega)$ to $L^{2}(\Omega)$, and denote $A=-\Delta$ in the following. For any $r\geq 0$, denote the space $\dot{H}^{r}(\Omega)=\{\vartheta\in L^{2}(\Omega):A^{\frac{r}{2}}\vartheta\in L^{2}(\Omega)\}$ with the norm [2]

[TABLE]

where ${(\lambda_{j},\varphi_{j})}$ are the eigenvalues ordered non-decreasingly and the corresponding eigenfunctions normalized in the $L^{2}(\Omega)$ norm of $-\Delta$ on the domain $\Omega$ with a zero Dirichlet boundary condition. Thus $\dot{H}^{0}(\Omega)=L^{2}(\Omega)$, $\dot{H}^{1}(\Omega)=H^{1}_{0}(\Omega)$, and $\dot{H}^{2}(\Omega)=H^{2}(\Omega)\bigcap H^{1}_{0}(\Omega)$.

After that, for $\kappa>0$ and $\pi/2<\theta<\pi$, we define sectors $\Sigma_{\theta}$ and $\Sigma_{\theta,\kappa}$ in the complex plane $\mathbb{C}$ as

[TABLE]

and the contour $\Gamma_{\theta,\kappa}$ is defined by

[TABLE]

oriented with an increasing imaginary part, where $\mathbf{i}$ denotes the imaginary unit and $\mathbf{i}^{2}=-1$. According to the results in [19], the system (1.1) can be rewritten as

[TABLE]

where $a=\frac{1-m}{2m-1}$ and $m$ is defined in (1.2); and the system (2.1) has the solution of the form

[TABLE]

where ‘$\tilde{~{}}$’ stands for taking Laplace transform,

[TABLE]

and

[TABLE]

Then we provide some estimates related to (2.2) and (2.3), which will be used in the error estimates.

Lemma 2.1** ([19]).**

When $z\in\Sigma_{\theta,\kappa}$, $\pi/2<\theta<\pi$, and $\kappa>\max\left(2|a|^{1/\alpha_{1}},2|a|^{1/\alpha_{2}}\right)$, we have the estimates

[TABLE]

where $H(z)$, $H_{\alpha_{1}}(z)$ and $H_{\alpha_{2}}(z)$ are defined in (2.2) and (2.3), respectively.

Lemma 2.2**.**

When $z\in\Sigma_{\theta,\kappa}$, $\pi/2<\theta<\pi$, and $\kappa>\max\left(2|a|^{1/\alpha_{1}},2|a|^{1/\alpha_{2}}\right)$, there are the estimates of $H(z)$, $H_{\alpha_{1}}(z)$ and $H_{\alpha_{2}}(z)$ in (2.2) and (2.3),

[TABLE]

Proof.

First, consider the estimate of $\|(a+A)H(z)\|$. Obviously, there exist the equalities

[TABLE]

To estimate $(a+A)H(z)$, one can estimate $(z^{\alpha_{1}}+a+A)H(z)$, $z^{\alpha_{1}}H(z)$, $(z^{\alpha_{2}}+a+A)H(z)$, and $z^{\alpha_{2}}H(z)$. As for $(z^{\alpha_{1}}+a+A)H(z)$, let

[TABLE]

which results in

[TABLE]

Performing $L_{2}$ norm on both sides of the above equality and using Lemma 2.1, we have

[TABLE]

Taking $\kappa$ sufficiently large leads to

[TABLE]

Similarly, we also have

[TABLE]

From Lemma 2.1, there exist

[TABLE]

Thus, we obtain the estimate of $\|(a+A)H(z)\|$. The estimates of $\|(a+A)H_{\alpha_{1}}(z)\|$ and $\|(a+A)H_{\alpha_{2}}(z)\|$ can be similarly obtained. ∎

3. Fast evaluation of the Riemann-Liouville fractional derivative

In this section, we provide the fast BE and fast SBD approximations based on the convolution quadrature of the Riemann-Liouville fractional derivative. Suppose that $N$ is the total number of time steps, and the time step size $\tau=T/N$ and $t_{n}=n\tau$, $1\leq n\leq N$.

Let’s start from the integral representation of the power function.

Lemma 3.1** ([7]).**

For any $\beta>0$, there is

[TABLE]

By using the property of convolution and Lemma 3.1, Eq. (1.3) can be rewritten as

[TABLE]

Taking the Laplace transform on the both sides of (3.1), we obtain

[TABLE]

Following the classical convolution quadrature, we only need to take $z=\delta(\zeta)$ on the left side of (3.2) to get the discretization of $~{}_{0}D^{\alpha}_{t}G(t)$, where $\delta(\zeta)$ is the generating function given by BE or SBD methods, i.e., $\delta(\zeta)=(1-\zeta)/\tau$ or $\delta(\zeta)=((1-\zeta)+(1-\zeta)^{2}/2)/\tau$. Here we get the integral representation of the weights generated by convolution quadrature according to (3.2) to speed up the evaluation.

3.1. Fast BE discretization

Firstly we take $z=\delta(\zeta)=(1-\zeta)/\tau$ in (3.2) and get

[TABLE]

Setting

[TABLE]

and using the fact

[TABLE]

we obtain from (3.3)

[TABLE]

[TABLE]

For $i\geq 2$, simple calculation (see Appendix A) leads to

[TABLE]

where

[TABLE]

Then the Riemann-Liouville fractional derivative can be discretized as

[TABLE]

where

[TABLE]

$\{w^{\alpha}_{j}\}_{j=1}^{N_{p}}$ denote the integration weights, $\{s^{\alpha}_{j}\}_{j=1}^{N_{p}}$ signify the integration points, $N_{p}$ is the number of integration points, and $\{\epsilon^{\alpha}_{1,i}\}_{i=2}^{\infty}$ indicate the errors caused by integral approximation. Obviously, $\epsilon^{\alpha}_{1,0}=\epsilon^{\alpha}_{1,1}=0$. To make $|\epsilon^{\alpha}_{1,i}|$ small enough, we can use the Gauss-Jacobi rule to generate $w^{\alpha}_{j}$ and $s^{\alpha}_{j}$, and for $2<i\leq 2N_{p}+1$, $d^{\alpha}_{1,i}$ can be exactly approximated, i.e., $\epsilon^{\alpha}_{1,i}=0$.

To get a fast evaluation for (3.8), we rewrite it as

[TABLE]

where

[TABLE]

and we call it history part. According to Eq. (3.7), it’s easy to know that $\{\bar{d}_{1,i}(s)\}_{i=2}^{\infty}$ is a geometric sequence, so we can get

[TABLE]

Thus we can get $\mathcal{G}^{1}_{hist,j}(t_{n})$ from $\mathcal{G}^{1}_{hist,j}(t_{n-1})$ and $G(t_{n-2})$ instead of calculating the sum of $\bar{d}_{1,n-i}(s^{\alpha}_{j})G(t_{i})$. So the computation time is reduced from $\mathcal{O}(N^{2})$ to $\mathcal{O}(NN_{p})$ and the total memory requirement is cut down from $\mathcal{O}(N)$ to $\mathcal{O}(N_{p})$, where $N$ stands for the total number of time steps and $N_{p}$ is the number of integration points.

3.2. Fast SBD discretization

In this subsection, we take $z=\delta(\zeta)=((1-\zeta)+(1-\zeta)^{2}/2)/\tau$ in Eq. (3.2), which leads to

[TABLE]

Set

[TABLE]

By simple calculation, we have

[TABLE]

where $\sigma$ is the solution of $\sigma^{2}+2\sigma+2\tau s=0$. Then from (3.11) we obtain

[TABLE]

[TABLE]

For $i\geq 3$, by simple calculations and Jordan’s Lemma (see Appendix B), we have

[TABLE]

where

[TABLE]

Thus $d^{\alpha}_{2,i}$ can be approximated as

[TABLE]

where $\{w^{\alpha}_{1,j}\}_{j=1}^{N_{p,1}}$, $\{w^{\alpha}_{2,j}\}_{j=1}^{N_{p,2}}$ denote the integration weights, $\{s^{\alpha}_{1,j}\}_{j=1}^{N_{p,1}}$, $\{s^{\alpha}_{2,j}\}_{j=1}^{N_{p,2}}$ are the integral points, $N_{p,1}$, $N_{p,2}$ signify the number of integration points, and $\{\epsilon^{\alpha}_{2,i}\}_{i=3}^{\infty}$ indicate the errors caused by the integration approximation. For convenience, we set $\epsilon^{\alpha}_{2,0}=\epsilon^{\alpha}_{2,1}=\epsilon^{\alpha}_{2,2}=0$. To make $|\epsilon^{\alpha}_{2,i}|$ small enough, the Gauss-Jacobi rule can be used to obtain $w^{\alpha}_{1,j}$, $w^{\alpha}_{2,j}$ and $s^{\alpha}_{1,j}$, $s^{\alpha}_{2,j}$.

Remark 3.2*.*

According to (3.12), the weights $d^{\alpha}_{2,i}~{}(i=0,1,\cdots)$ are real numbers, but we find that $\bar{d}^{\alpha}_{2,i}(s)$ for $s<-\frac{1}{3}$ and $w^{\alpha}_{1,j}$ are complex numbers, which are caused by two terms $(1+3s)^{\alpha}$ and $(-1)^{-\alpha}$, respectively. So, to reduce the computation time but without losing precision, we use the real parts of $\bar{d}^{\alpha}_{2,i}(s)$ and $w^{\alpha}_{1,j}$ to accomplish the simulation.

Then, the Riemann-Liouville fractional derivative can be discretized as

[TABLE]

where $N_{s}$ is a parameter that ensures the accuracy of the discretization.

Remark 3.3*.*

Here the reason that we introduce the parameter $N_{s}$ is that $d_{2,i}^{\alpha}$ can’t be approximated effectively by the Gauss-Jacobi rule when $i$ is small, so we start from the $N_{s}$-th term to approximate $d_{2,i}^{\alpha}$, i.e., we still use the first $N_{s}$ weights $d_{2,i}^{\alpha}$ generated by convolution quadrature in the fast SBD discretization and $\epsilon^{\alpha}_{2,i}=0$ for $i<N_{s}$. The detailed discussions on the value of $N_{s}$ will be presented in the numerical experiments.

To get the fast evaluation of (3.17), we rewrite it as

[TABLE]

where

[TABLE]

and we also call the two terms history parts. Using the property of geometrical sequences $\{\bar{d}^{1}_{2,i}(s)\}_{i=2}^{\infty}$ and $\{\bar{d}^{2}_{2,i}(s)\}_{i=2}^{\infty}$ defined in (3.15), we obtain

[TABLE]

Thus we can get $\mathcal{G}^{2,1}_{hist,j}(t_{n})$ and $\mathcal{G}^{2,2}_{hist,j}(t_{n})$ from $\mathcal{G}^{2,1}_{hist,j}(t_{n-1})$, $\mathcal{G}^{2,2}_{hist,j}(t_{n-1})$ and $G(t_{n-N_{s}})$ instead of calculating the sum of $\bar{d}^{1}_{2,n-i}(s^{\alpha}_{1,j})G(t_{i})$ and the sum of $\bar{d}^{2}_{2,n-i}(s^{\alpha}_{2,j})G(t_{i})$. So the computation time is reduced from $\mathcal{O}(N^{2})$ to $\mathcal{O}(N(N_{p,1}+N_{p,2}+N_{s}))$ and the total memory requirement is reduced from $\mathcal{O}(N)$ to $\mathcal{O}(N_{p,1}+N_{p,2}+N_{s})$, where $N$ is the total number of the time steps, $N_{p,1}$ and $N_{p,2}$ are the number of the integral points, and $N_{s}$ is a parameter that ensures the accuracy of the approximation.

4. Error analysis

Now, we apply the fast BE and SBD algorithms developed in Section 3 to solve the system of fractional partial differential equations (2.1) and give error analyses of the fast BE and SBD schemes, respectively.

4.1. Error estimates for the fast BE scheme

According to [19], we have the following BE scheme

[TABLE]

where $\bar{G}^{n}_{1}$, $\bar{G}^{n}_{2}$ are the numerical solutions of $G_{1}$, $G_{2}$ at time $t_{n}$. According to (3.9), we modify the system (4.1) as the fast BE scheme, i.e.,

[TABLE]

where $G^{n}_{1}$, $G^{n}_{2}$ are the numerical solutions of $G_{1}$, $G_{2}$ at time $t_{n}$. From Eq. (3.8), the system (4.2) can be rewritten as

[TABLE]

Remark 4.1*.*

Following [18], we omit $G(t_{0})$ in (3.8) when we construct the BE scheme of system (2.1), which is helpful in the approximation of the Riemann-Liouville fractional derivative in the system (4.3). Accordingly, the history part defined in (3.10) is modified as

[TABLE]

Lemma 4.2** ([12]).**

For any $\theta\in(\pi/2,\pi)$, where $\theta=\arg(z)$, there exist positive constants $c_{1}$, $c_{2}$ such that

[TABLE]

where $\delta(\zeta)=(1-\zeta)/\tau$ or $\delta(\zeta)=((1-\zeta)+(1-\zeta)^{2}/2)/\tau$.

According to [19], we have the following estimates between (2.1) and (4.1).

Theorem 4.3**.**

Let $G_{1}$, $G_{2}$ and $\bar{G}^{n}_{1}$, $\bar{G}^{n}_{2}$ be, respectively, the solutions of the systems (2.1) and (4.1). Then

[TABLE]

Now we provide the regularity of solutions for the system (4.1).

Theorem 4.4**.**

Let $\bar{G}^{n}_{1}$, $\bar{G}^{n}_{2}$ be the solutions of the system (4.1). Then we have

[TABLE]

for $\nu=0,1$; furthermore, if $G_{1}(0),G_{2}(0)\in H^{1}_{0}(\Omega)\bigcap H^{2}(\Omega)$, there exist

[TABLE]

Proof.

Here we take $\kappa\geq 1/t$ for given $t$ and ensure that $\kappa$ is large enough to satisfy the conditions in Lemmas 2.1 and 2.2. To get the solutions of the system (4.1), we multiply by $\zeta^{n}$ and sum from $1$ to $\infty$ for both sides of the first two equations in (4.1) and get

[TABLE]

According to (3.4), we have

[TABLE]

which result in, after simple calculations,

[TABLE]

where $H$, $H_{\alpha_{1}}$, and $H_{\alpha_{2}}$ are defined by (2.2) and (2.3). By (4.4), for $\xi_{\tau}=e^{-\tau(\kappa+1)}$, there is

[TABLE]

Taking $\zeta=e^{-z\tau}$ leads to

[TABLE]

where $\Gamma^{\tau}=\{z=\kappa+1+\mathbf{i}y:y\in\mathbb{R}~{}{\rm and}~{}|y|\leq\pi/\tau\}$. Next, we deform the contour $\Gamma^{\tau}$ to $\Gamma^{\tau}_{\theta,\kappa}=\{z\in\mathbb{C}:\kappa\leq|z|\leq\frac{\pi}{\tau\sin(\theta)},|\arg z|=\theta\}\bigcup\{z\in\mathbb{C}:|z|=\kappa,|\arg z|\leq\theta\}$. Thus

[TABLE]

Combining Lemmas 2.1 and 4.2, we obtain

[TABLE]

where $\Re(z)$ denotes the real part of $z$. Similarly, we have

[TABLE]

Multiplying the operator $A$ on both sides of (4.5), we have

[TABLE]

When $G_{1}(0),G_{2}(0)\in L^{2}(\Omega)$, according to Lemmas 2.1 and 4.2, there is

[TABLE]

When $G_{1}(0),G_{2}(0)\in\dot{H}^{2}(\Omega)$, we obtain

[TABLE]

Similarly, when $G_{1}(0),G_{2}(0)\in L^{2}(\Omega)$, we have

[TABLE]

When $G_{1}(0),G_{2}(0)\in\dot{H}^{2}(\Omega)$, there exists

[TABLE]

∎

Theorem 4.5**.**

Let $\bar{G}^{n}_{1}$, $\bar{G}^{n}_{2}$ and $G^{n}_{1}$, $G^{n}_{2}$ be, respectively, the solutions of the systems (4.1) and (4.3). Then

[TABLE]

where $\epsilon_{i}=\max(|\epsilon^{1-\alpha_{1}}_{1,i}|,|\epsilon^{1-\alpha_{2}}_{1,i}|)$.

Proof.

In this proof, we take $\kappa\geq 1/t$ for given $t$ and ensure that $\kappa$ is large enough to satisfy the conditions in Lemmas 2.1 and 2.2. Subtracting (4.3) from (4.1) and denoting $e^{n}_{1}=\bar{G}^{n}_{1}-G^{n}_{1}$, $e^{n}_{2}=\bar{G}^{n}_{2}-G^{n}_{2}$, we have

[TABLE]

Multiplying $\zeta^{n}$ and summing from $1$ to $\infty$ for the both sides of equations in (4.6) lead to

[TABLE]

Introducing

[TABLE]

and using (3.4), we have

[TABLE]

Making simple calculations leads to

[TABLE]

Denote

[TABLE]

Consider the estimate of $E^{n}_{0,1,1}$. Taking $\xi_{\tau}=e^{-\tau(\kappa+1)}$ results in

[TABLE]

which leads to

[TABLE]

Letting $\zeta=e^{-z\tau}$, we obtain

[TABLE]

where $\Gamma^{\tau}=\{z=\kappa+1+\mathbf{i}y:y\in\mathbb{R}~{}{\rm and}~{}|y|\leq\pi/\tau\}$. Next we deform the contour $\Gamma^{\tau}$ to $\Gamma^{\tau}_{\theta,\kappa}=\{z\in\mathbb{C}:\kappa\leq|z|\leq\frac{\pi}{\tau\sin(\theta)},|\arg z|=\theta\}\bigcup\{z\in\mathbb{C}:|z|=\kappa,|\arg z|\leq\theta\}$. Then

[TABLE]

Using Lemmas 2.1 and 4.2, we have

[TABLE]

Similarly, we have

[TABLE]

Thus

[TABLE]

Combining Grönwall’s inequality and Theorem 4.4 leads to

[TABLE]

∎

Combining Theorem 4.3 and Theorem 4.5, we get the following error estimates for the fast BE scheme.

Theorem 4.6**.**

Let $G_{1}$, $G_{2}$ and $G^{n}_{1}$, $G^{n}_{2}$ be the solutions of the systems (2.1) and (4.3), respectively. Then we have the estimates, if $G_{1}(0)$, $G_{2}(0)\in L^{2}(\Omega)$,

[TABLE]

where $\epsilon_{i}=\max(|\epsilon^{1-\alpha_{1}}_{1,i}|,|\epsilon^{1-\alpha_{2}}_{1,i}|)$.

4.2. Error estimates for the fast SBD scheme

Here, we first provide the SBD scheme of (2.1) by SBD convolution quadrature and give the error estimate of the SBD scheme. Then we present the error estimate of the fast SBD scheme. According to [18, 11], to keep the accuracy of the scheme, one needs to modify the discretization (3.17). Namely, denoting $\bar{\partial}^{\alpha}_{\tau}$ as the discretization of ${}_{0}D^{\alpha}_{t}$ and letting $\partial^{-1}_{t}$ be the integration on time, from (2.1) we obtain

[TABLE]

According to the fact $(0,3/2,1,1,\ldots)=\bar{\partial}_{\tau}\partial^{-1}_{t}1$ [11] and setting $\bar{G}^{-1}_{1}=G_{1}(0)$, $\bar{G}^{-1}_{2}=G_{2}(0)$, the SBD scheme of (2.1) can be written as

[TABLE]

where $\bar{G}^{n}_{1}$, $\bar{G}^{n}_{2}$ are the numerical solutions of $G_{1}$, $G_{2}$ at $t_{n}$. Similarly using Eq. (3.17), and noting that $\epsilon^{\alpha}_{2,0}=\epsilon^{\alpha}_{2,1}=\epsilon^{\alpha}_{2,2}=0$, we obtain the fast SBD scheme

[TABLE]

where $G^{n}_{1}$, $G^{n}_{2}$ are the numerical solutions of $G_{1}$, $G_{2}$ at time $t_{n}$.

Remark 4.7*.*

To keep the accuracy of SBD scheme (4.7), the discretization of $\,{}_{0}D^{\alpha}_{t}G(t)$ should be modified and correspondingly the definition of history parts in (3.18) should be changed as

[TABLE]

Next, we provide the error estimates between (2.1) and (4.7).

Theorem 4.8**.**

Let $G_{1}$, $G_{2}$ and $\bar{G}^{n}_{1}$, $\bar{G}^{n}_{2}$ be, respectively, the solutions of the systems (2.1) and (4.7). Then

[TABLE]

Proof.

Here we take $\kappa\geq 1/t$ for given $t$ and ensure that $\kappa$ is large enough to satisfy the conditions in Lemmas 2.1 and 2.2. Introduce $V_{1}(t)=G_{1}(t)-G_{1}(0)$, $V_{2}(t)=G_{2}(t)-G_{2}(0)$, $\bar{V}^{n}_{1}=\bar{G}^{n}_{1}-G_{1}(0)$, $\bar{V}^{n}_{2}=\bar{G}^{n}_{2}-G_{2}(0)$; and denote $\mathbf{V}=[V_{1},V_{2}]^{T}$. Thus (2.1) can be written as

[TABLE]

Therefore we get the solutions of system (4.9) as

[TABLE]

where ‘$\tilde{~{}}$’ stands for taking Laplace transform. Let $\bar{V}^{-1}_{1}=\bar{V}^{-1}_{2}=0$, and then (4.7) can be rewritten as

[TABLE]

Multiplying $\zeta^{n}$ and summing from $1$ to $\infty$ for both sides of the first two formulas in (4.11), we obtain

[TABLE]

Using Eq. (3.12) and the facts

[TABLE]

and $\delta(\zeta)=\frac{(1-\zeta)+(1-\zeta)^{2}/2}{\tau}$, we obtain

[TABLE]

Thus, after simple calculations, we have

[TABLE]

To get error estimate between $V_{1}(t_{n})$ and $\bar{V}^{n}_{1}$, we need to get the error estimates between $V_{1,i}(t_{n})$ and $\bar{V}^{n}_{1,i}$ $(i=1,2,3,4)$. Then we consider the error estimate between $V_{1,1}(t_{n})$ and $\bar{V}^{n}_{1,1}$. Using Eq. (4.12) and denoting $\mu(\zeta)$ as $\mu(\zeta)=\tau\delta(\zeta)\nu(\zeta)=\frac{\zeta(3-\zeta)^{2}}{4}$, when $\xi_{\tau}=e^{-(\kappa+1)\tau}$, we have

[TABLE]

Taking $\zeta=e^{-z\tau}$, we obtain

[TABLE]

where $\Gamma^{\tau}=\{z=\kappa+1+\mathbf{i}y:y\in\mathbb{R}~{}{\rm and}~{}|y|\leq\pi/\tau\}$. Next we deform the contour $\Gamma^{\tau}$ to $\Gamma^{\tau}_{\theta,\kappa}=\{z\in\mathbb{C}:\kappa\leq|z|\leq\frac{\pi}{\tau\sin(\theta)},|\arg z|=\theta\}\bigcup\{z\in\mathbb{C}:|z|=\kappa,|\arg z|\leq\theta\}$. Thus

[TABLE]

Taking the inverse Laplace transform for (4.10), and combining (4.13), we obtain

[TABLE]

According to Lemma 2.2, we have

[TABLE]

Using $\|\frac{d}{dz}H_{\alpha_{1}}(z)z^{-1}(-(a+A))\|\leq C|z|^{-2}$, $\delta(e^{-z\tau})=z+\mathcal{O}(\tau^{2}z^{3})$ and the mean value theorem, we obtain

[TABLE]

Combining the fact $\mu(e^{-z\tau})=1+O(z^{2}\tau^{2})$ [18, 12] and Lemmas 2.2 and 4.2 leads to

[TABLE]

Thus

[TABLE]

Similarly, we can get

[TABLE]

Using the fact $T/t>1$, we have

[TABLE]

Similarly, there is

[TABLE]

Thus, the proof is completed. ∎

Then we have the following regularity estimates of the solutions of system (4.7).

Theorem 4.9**.**

Let $\bar{G}^{n}_{1}$, $\bar{G}^{n}_{2}$ be the solutions of the systems (4.7). Then

[TABLE]

Proof.

The proof is similar to the proof of Theorem 4.4. ∎

Next we provide the error estimate between $\bar{G}^{n}_{1}$, $\bar{G}^{n}_{2}$ and $G^{n}_{1}$, $G^{n}_{2}$, which are the solutions of the systems (4.7) and (4.8), respectively.

Theorem 4.10**.**

Let $\bar{G}^{n}_{1}$, $\bar{G}^{n}_{2}$ and $G^{n}_{1}$, $G^{n}_{2}$ be, respectively, the solutions of the systems (4.7) and (4.8). Then

[TABLE]

where $\epsilon_{i}=\max(|\epsilon^{1-\alpha_{1}}_{2,i}|,|\epsilon^{1-\alpha_{2}}_{2,i}|)$.

Proof.

In the proof, we take $\kappa\geq 1/t$ for given $t$ and ensure $\kappa$ is large enough to satisfy the conditions in Lemmas 2.1 and 2.2. Subtracting (4.8) from (4.7), denoting $e^{n}_{1}=\bar{G}^{n}_{1}-G^{n}_{1}$, $e^{n}_{2}=\bar{G}^{n}_{2}-G^{n}_{2}$, $e^{-1}_{1}=e^{-1}_{2}=0$ and using $e^{0}_{1}=e^{0}_{2}=0$, we have

[TABLE]

By simple calculation, we obtain

[TABLE]

Multiplying $\zeta^{n}$ and summing from $1$ to $\infty$ for the both sides of equations in (4.14) lead to

[TABLE]

Introducing

[TABLE]

and using Eq. (3.12), we obtain

[TABLE]

Thus

[TABLE]

Denote

[TABLE]

Here we first consider the estimate of $E^{n}_{0,1,1}$. Taking $\xi_{\tau}=e^{-\tau(\kappa+1)}$, we get

[TABLE]

which leads to

[TABLE]

Letting $\zeta=e^{-z\tau}$, there is

[TABLE]

where $\Gamma^{\tau}=\{z=\kappa+1+\mathbf{i}y:y\in\mathbb{R}~{}{\rm and}~{}|y|\leq\pi/\tau\}$. Next, we deform the contour $\Gamma^{\tau}$ to $\Gamma^{\tau}_{\theta,\kappa}=\{z\in\mathbb{C}:\kappa\leq|z|\leq\frac{\pi}{\tau\sin(\theta)},|\arg z|=\theta\}\bigcup\{z\in\mathbb{C}:|z|=\kappa,|\arg z|\leq\theta\}$. Thus

[TABLE]

Using Lemmas 2.1 and 4.2, we have

[TABLE]

Similarly, we have

[TABLE]

Thus

[TABLE]

Combining Grönwall inequality and Theorem 4.9 leads to

[TABLE]

∎

Combining Theorems 4.8 and 4.10, we get the following error estimates for the fast SBD scheme.

Theorem 4.11**.**

Let $G_{1}$, $G_{2}$ and $G^{n}_{1}$, $G^{n}_{2}$ be the solutions of the systems (2.1) and (4.8), respectively. Then we have the estimates, if $G_{1}(0)$, $G_{2}(0)\in L^{2}(\Omega)$,

[TABLE]

where $\epsilon_{i}=\max(|\epsilon^{1-\alpha_{1}}_{2,i}|,|\epsilon^{1-\alpha_{2}}_{2,i}|)$.

5. Numerical experiments

In this section, we verify the effectiveness of the fast algorithms by comparing with the classical convolution quadrature schemes. Here, we denote $h$ as the mesh size and define

[TABLE]

where $G^{N}_{1,\tau}$ and $G^{N}_{2,\tau}$, respectively, signify the numerical solutions of $G_{1}$ and $G_{2}$ at the fixed time $t_{N}$ with time step size $\tau$. The temporal convergence rates can be calculated by

[TABLE]

In the numerical experiments, the following two groups of initial values are used:

1. (a)

[TABLE] 2. (b)

[TABLE]

where $\chi_{(a,b)}$ denotes the characteristic function on $(a,b)$.

5.1. Performance of fast BE scheme

We, respectively, use the BE scheme and fast BE scheme to solve the system (2.1) with $a=2$. Use the numerical solution with $\tau=1/3200$ and $h=1/256$ as the ‘exact’ solution. Tables 1 and 2 give the $L_{2}$ errors, convergence rates and CPU time for solving the system (2.1) with the initial values (a) and (b) for different $\alpha_{1}$ and $\alpha_{2}$ respectively, which show that the fast BE scheme has the same convergence rates as BE scheme and it takes much less CPU time when $\tau$ is small.

5.2. Performance of fast SBD scheme

Here we first discuss the choice of $N_{s}$. According to (3.16), we need to make $|\epsilon^{\alpha}_{2,i}|$ small enough such that the approximation of $d^{\alpha}_{2,i}$ effective, and there are two ways to achieve it: the increase of the number of the integral points and the selection of a suitable parameter $N_{s}$. For small $N_{s}$, one needs to use a large number of integral points to ensure the accuracy of the approximation of $d^{\alpha}_{2,i}$ because of the low convergence rates of the Gauss-Jacobi rule, which increases the computation time and storage cost tremendously. Figure 1 shows the change of the error $|\epsilon^{\alpha}_{2,i}|$ as $i$ increases when $T=1$, $\tau=1/1000$, $\alpha=0.3$, and $N_{p,1}+N_{p,2}=62$. It is found that when $i<15$, $d^{\alpha}_{2,i}$ can’t be well approximated, i.e., $|\epsilon^{\alpha}_{2,i}|$ is large, so we start from the $15$-th term to approximate $d^{\alpha}_{2,i}$, that is, we take $N_{s}=15$ to ensure the accuracy of the approximation. The same thing happens in Figure 2 when $T=1$, $\tau=1/1000$, $\alpha=0.8$, and $N_{p,1}+N_{p,2}=82$, hence we take $N_{s}=17$ in the fast SBD scheme. In a word, choosing a suitable value for $N_{s}$ not only ensures the accuracy of the scheme, but also saves the computation time.

Then we solve the system (2.1) by the SBD scheme and fast SBD scheme, respectively, with $a=-1$. Use the numerical solution with $\tau=1/640$ and $h=1/1024$ as the ‘exact’ solution. Tables 3 and 4 provide the $L_{2}$ errors and convergence rates for solving the system (2.1), respectively, with the initial values (a) and (b) for different $\alpha_{1}$ and $\alpha_{2}$, which show that the fast SBD scheme has the same convergence rates as SBD scheme.

Finally, we use the SBD scheme and fast SBD scheme to solve the system (2.1) with $a=-1$ to verify the efficiency of the latter scheme. Use the numerical solution with $\tau=1/4000$ and $h=1/1024$ as the ‘exact’ solution. Table 5 shows the $L_{2}$ error and CPU time when $T=10$, $\tau=1/2000$ and $h=1/1024$ for different $\alpha_{1}$ and $\alpha_{2}$. The $L^{2}$ errors of SBD scheme are close to the errors of fast SBD scheme but the CPU times of fast SBD scheme are much less than SBD scheme for different $\alpha_{1}$ and $\alpha_{2}$, which show that the fast algorithm can greatly reduce computation time.

Conclusion

The fast algorithms based on BE and SBD convolution quadratures are developed to solve the homogeneous fractional Fokker-Planck equations with two internal states, and they, respectively, have the first- and second-order convergence rates. One of the advantages of the provided fast algorithms is that the assumption of the regularity of the solution in time is not required. The effectiveness of the fast algorithms is verfied by numerical experiments.

Acknowledgements

This work was supported by the National Natural Science Foundation of China under grant no. 11671182, and the Fundamental Research Funds for the Central Universities under grant no. lzujbky-2018-ot03.

Appendix A Derivation of (3.6)

For $i\geq 2$, letting $\tau s=t$ in Eq. (3.5), we have

[TABLE]

Then taking $\eta=1/t$ leads to

[TABLE]

Lastly, we take $\eta=(s+1)/2$ and get

[TABLE]

Appendix B Derivation of (3.14)

For $i\geq 3$, taking $\tau s=t$ in (3.13) results in

[TABLE]

and $\sigma^{2}+2\sigma+2\tau s=\sigma^{2}+2\sigma+2t=0$, which lead to

[TABLE]

with $\Gamma_{1}=\{z=-1+\sqrt{1-2t},~{}0\leq t\leq\infty\}$. It is easy to find that for any $n\geq 3$ and $-\pi<\theta_{1},\theta_{2}<\pi$ when $0<\alpha<1$,

[TABLE]

and according to Jordan’s Lemma, we have

[TABLE]

Taking $\sigma=1/t$ in the above equation results in

[TABLE]

Further letting $t=(s+1)/4$ leads to

[TABLE]

Similarly, for the second part of Eq. (3.13), we can rewrite it as

[TABLE]