Numerical algorithm for the space-time fractional Fokker-Planck system with two internal states

TL;DR

This paper develops and analyzes a numerical scheme combining $L_1$ time discretization and finite element spatial approximation for a two-state fractional Fokker-Planck system with fractional Laplacian, providing error estimates and numerical validation.

Contribution

It introduces a novel numerical method for the two-state fractional Fokker-Planck system with fractional Laplacian, including rigorous error analysis without regularity assumptions.

Findings

The scheme achieves optimal error estimates in space and time.

Numerical experiments confirm the theoretical convergence rates.

The method effectively handles the nonlocal fractional Laplacian operator.

Abstract

The fractional Fokker-Planck system with multiple internal states is derived in [Xu and Deng, Math. Model. Nat. Phenom., $\mathbf{13}$, 10 (2018)], where the space derivative is Laplace operator. If the jump length distribution of the particles is power law instead of Gaussian, the space derivative should be replaced with fractional Laplacian. This paper focuses on solving the two state Fokker-Planck system with fractional Laplacian. We first provide a priori estimate for this system under different regularity assumptions on the initial data. Then we use $L_1$ scheme to discretize the time fractional derivatives and finite element method to approximate the fractional Laplacian operators. Furthermore, we give the error estimates for the space semidiscrete and fully discrete schemes without any assumption on regularity of solutions. Finally, the effectiveness of the designed scheme is verified by numerical experiments.

Tables

Table 1. Table 1: L 2 subscript 𝐿 2 L_{2} errors and convergence rates with s 1 = s 2 = s < 1 / 2 subscript 𝑠 1 subscript 𝑠 2 𝑠 1 2 s_{1}=s_{2}=s<1/2 and the initial condition ( a )

$s$	$1/h$	50	100	200	400	800
	$E_{1,h}$	1.293E-02	8.631E-03	5.752E-03	3.829E-03	2.546E-03
0.1		Rate	0.5834	0.5854	0.5872	0.5888
	$E_{2,h}$	9.139E-03	6.038E-03	3.986E-03	2.630E-03	1.734E-03
		Rate	0.5981	0.5991	0.5999	0.6005
	$E_{1,h}$	5.861E-03	3.486E-03	2.071E-03	1.229E-03	7.298E-04
0.25		Rate	0.7496	0.7514	0.7522	0.7523
	$E_{2,h}$	3.795E-03	2.247E-03	1.331E-03	7.890E-04	4.680E-04
		Rate	0.7562	0.7553	0.7544	0.7535
	$E_{1,h}$	2.334E-03	1.247E-03	6.681E-04	3.585E-04	1.925E-04
0.4		Rate	0.9044	0.9006	0.8982	0.8971
	$E_{2,h}$	1.468E-03	7.863E-04	4.218E-04	2.264E-04	1.216E-04
		Rate	0.9010	0.8985	0.8974	0.8975

Table 2. Table 2: L 2 subscript 𝐿 2 L_{2} errors and convergence rates with different s 1 , s 2 subscript 𝑠 1 subscript 𝑠 2 s_{1},s_{2} and the initial condition ( a )

$(s_1,s_2)$	$1/h$	50	100	200	400	800
	$E_{1,h}$	1.173E-02	7.797E-03	5.181E-03	3.441E-03	2.284E-03
(0.1,0.2)		Rate	0.5894	0.5899	0.5905	0.5913
	$E_{2,h}$	6.456E-03	3.988E-03	2.460E-03	1.516E-03	9.340E-04
		Rate	0.6949	0.6969	0.6982	0.6992
	$E_{1,h}$	4.105E-03	2.349E-03	1.345E-03	7.707E-04	4.420E-04
(0.3,0.4)		Rate	0.8051	0.8045	0.8034	0.8023
	$E_{2,h}$	1.853E-03	9.921E-04	5.320E-04	2.856E-04	1.534E-04
		Rate	0.9017	0.8989	0.8973	0.8968
	$E_{1,h}$	5.780E-04	2.754E-04	1.326E-04	6.425E-05	3.109E-05
(0.6,0.7)		Rate	1.0695	1.0547	1.0453	1.0472
	$E_{2,h}$	2.347E-04	1.092E-04	5.172E-05	2.479E-05	1.193E-05
		Rate	1.1032	1.0785	1.0613	1.0546
	$E_{1,h}$	1.143E-04	4.905E-05	2.194E-05	1.013E-05	4.798E-06
(0.8,0.9)		Rate	1.2207	1.1607	1.1149	1.0780
	$E_{2,h}$	3.297E-05	1.330E-05	5.799E-06	2.668E-06	1.269E-06
		Rate	1.3098	1.1974	1.1202	1.0714

Table 3. Table 3: L 2 subscript 𝐿 2 L_{2} errors and convergence rates with different s 1 , s 2 subscript 𝑠 1 subscript 𝑠 2 s_{1},s_{2} and the initial condition ( b )

$(s_1,s_2)$	$1/h$	50	100	200	400	800
	$E_{1,h}$	5.682E-02	4.505E-02	3.629E-02	2.967E-02	2.459E-02
(0.1,0.2)		Rate	0.3348	0.3121	0.2904	0.2709
	$E_{2,h}$	4.932E-02	3.583E-02	2.609E-02	1.906E-02	1.398E-02
		Rate	0.4612	0.4578	0.4530	0.4475
	$E_{1,h}$	9.879E-03	6.352E-03	4.101E-03	2.658E-03	1.729E-03
(0.3,0.4)		Rate	0.6371	0.6310	0.6255	0.6206
	$E_{2,h}$	8.644E-03	5.090E-03	2.990E-03	1.752E-03	1.025E-03
		Rate	0.7640	0.7677	0.7709	0.7737
	$E_{1,h}$	9.208E-04	4.679E-04	2.370E-04	1.197E-04	5.972E-05
(0.6,0.7)		Rate	0.9766	0.9813	0.9852	1.0033
	$E_{2,h}$	8.295E-04	4.080E-04	2.017E-04	1.001E-04	4.950E-05
		Rate	1.0236	1.0163	1.0111	1.0158
	$E_{1,h}$	1.290E-04	5.828E-05	2.703E-05	1.280E-05	6.170E-06
(0.8,0.9)		Rate	1.1462	1.1083	1.0784	1.0530
	$E_{2,h}$	9.645E-05	4.097E-05	1.837E-05	8.559E-06	4.089E-06
		Rate	1.2353	1.1569	1.1020	1.0657

Table 4. Table 4: L 2 subscript 𝐿 2 L_{2} errors and convergence rates with different s 1 , s 2 subscript 𝑠 1 subscript 𝑠 2 s_{1},s_{2} and the initial condition ( b ) ( ν 1 = − 0.4 subscript 𝜈 1 0.4 \nu_{1}=-0.4 , ν 2 = − 0.3 subscript 𝜈 2 0.3 \nu_{2}=-0.3 )

$(s_1,s_2)$	$1/h$	50	100	200	400	800
	$E_{1,h}$	3.524E-02	2.633E-02	1.993E-02	1.528E-02	1.184E-02
(0.1,0.2)		Rate	0.4207	0.4016	0.3837	0.3677
	$E_{2,h}$	2.688E-02	1.785E-02	1.182E-02	7.810E-03	5.151E-03
		Rate	0.5908	0.5947	0.5978	0.6005
	$E_{1,h}$	6.941E-03	4.271E-03	2.631E-03	1.622E-03	1.001E-03
(0.3,0.4)		Rate	0.7006	0.6991	0.6977	0.6965
	$E_{2,h}$	5.211E-03	2.880E-03	1.586E-03	8.700E-04	4.759E-04
		Rate	0.8556	0.8609	0.8659	0.8705

Table 5. Table 5: L 2 subscript 𝐿 2 L_{2} errors and convergence rates with different s 1 , s 2 subscript 𝑠 1 subscript 𝑠 2 s_{1},s_{2} and the initial condition ( b ) ( ν 1 = 0 subscript 𝜈 1 0 \nu_{1}=0 , ν 2 = 0.4999 subscript 𝜈 2 0.4999 \nu_{2}=0.4999 )

$(s_1,s_2)$	$1/h$	50	100	200	400	800
	$E_{1,h}$	3.630E-04	1.980E-04	1.080E-04	5.894E-05	3.217E-05
(0.4,0.1)		Rate	0.8746	0.8742	0.8739	0.8737
	$E_{2,h}$	2.591E-02	2.269E-02	1.985E-02	1.736E-02	1.517E-02
		Rate	0.1916	0.1926	0.1936	0.1947
	$E_{1,h}$	3.417E-04	1.849E-04	1.001E-04	5.411E-05	2.924E-05
(0.4,0.2)		Rate	0.8858	0.8862	0.8869	0.8882
	$E_{2,h}$	1.122E-02	8.629E-03	6.622E-03	5.072E-03	3.878E-03
		Rate	0.3784	0.3818	0.3848	0.3874
	$E_{1,h}$	8.932E-05	4.412E-05	2.196E-05	1.100E-05	5.464E-06
(0.6,0.3)		Rate	1.0175	1.0065	0.9974	1.0095
	$E_{2,h}$	4.877E-03	3.321E-03	2.252E-03	1.521E-03	1.024E-03
		Rate	0.5544	0.5606	0.5660	0.5707

Table 6. Table 6: L 2 subscript 𝐿 2 L_{2} errors and convergence rates with different α 1 subscript 𝛼 1 \alpha_{1} , α 2 subscript 𝛼 2 \alpha_{2} and the initial condition ( a )

$({\alpha }_1,{\alpha }_2)$	$1/\tau$	100	200	400	800	1600
	$E_{1,\tau }$	3.980E-02	1.957E-02	9.745E-03	4.876E-03	2.444E-03
(0.3,0.6)		Rate	1.0241	1.0059	0.9988	0.9966
	$E_{2,\tau }$	1.038E-01	5.130E-02	2.565E-02	1.288E-02	6.478E-03
		Rate	1.0173	0.9999	0.9935	0.9920
	$E_{1,\tau }$	1.662E-02	8.178E-03	4.063E-03	2.027E-03	1.013E-03
(0.4,0.7)		Rate	1.0234	1.0092	1.0031	1.0007
	$E_{2,\tau }$	4.338E-02	2.145E-02	1.070E-02	5.358E-03	2.685E-03
		Rate	1.0159	1.0031	0.9982	0.9967
	$E_{1,\tau }$	8.279E-03	4.071E-03	2.019E-03	1.006E-03	5.020E-04
(0.25,0.8)		Rate	1.0242	1.0115	1.0055	1.0026
	$E_{2,\tau }$	2.198E-02	1.086E-02	5.410E-03	2.703E-03	1.352E-03
		Rate	1.0167	1.0059	1.0013	0.9996

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}(-\Delta {)}^{s_1},-{}_0{D}_t^{1-{\alpha }_2}(-\Delta {)}^{s_2})\mathbf{G}\mathrm{in}\Omega ,t\in (0,T],\\ & \mathbf{G}(\cdot ,0)={\mathbf{G}}_0\mathrm{in}\Omega ,\\ & \mathbf{G}=0\mathrm{in}{\Omega }^c,t\in [0,T],\end{array}$$

$${}_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);$$

$$(-\Delta {)}^{s_i}u(x)=c_{n,s_i}P.V.{\int }_{{\mathbb{R}}^n}\frac{u(x)-u(y)}{|x-y{|}^{n+2s_i}}dy,s_i\in (0,1),$$

$$\{\begin{array}{rl}& \frac{\partial }{\partial t}G=-{}_0{D}_t^{1-\alpha }(-\Delta {)}^{s}G\mathrm{in}\Omega ,t\in (0,T],\\ & G(\cdot ,0)=G_0\mathrm{in}\Omega ,\\ & G=0\mathrm{in}{\Omega }^c,t\in [0,T].\end{array}$$

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

$${\Gamma }_{\theta ,\kappa }=\{r{e}^{-\mathbf{i}\theta }:r\ge \kappa \}\cup \{\kappa e^{\mathbf{i}\psi }:|\psi |\le \theta \}\cup \{r{e}^{\mathbf{i}\theta }:r\ge \kappa \},$$

$$H^s(\Omega )=\left\{w\in L^2(\Omega ):|w{|}_{H^s(\Omega )}={\left(\int {\int }_{{\Omega }^2}\frac{|w(x)-w(y){|}^2}{|x-y{|}^{n+2s}}dxdy\right)}^{1/2}<\infty \right\}$$

$$H^s(\Omega )=\{w\in H^{\lfloor s\rfloor }(\Omega ):|D^{\alpha }w{|}_{H^{\sigma }(\Omega )}<\infty \mathrm{for}\mathrm{all}\alpha s.t.|\alpha |=\lfloor s\rfloor \},$$

$${\hat{H}}^s(\Omega )=\{w\in H^s({\mathbb{R}}^n):\mathbf{supp}w\subset \bar{\Omega }\},$$

$$⟨u,w{⟩}_s:=c_{n,s}\int {\int }_{({\mathbb{R}}^n\times {\mathbb{R}}^n)\backslash ({\Omega }^c\times {\Omega }^c)}\frac{(u(x)-u(y))(w(x)-w(y))}{|x-y{|}^{n+2s}}dydx.$$

$$\{\begin{array}{rl}(-\Delta {)}^{s}u& =g\mathrm{in}\Omega ,\\ u& =0\mathrm{in}{\Omega }^c,\end{array}$$

$$|u{|}_{H^{s+\gamma }({\mathbb{R}}^n)}\le C\Vert g{\Vert }_{H^r(\Omega )},$$

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

$$({\mathbf{M}}^T{)}^{-1}=\left[\begin{array}{cc}\frac{m}{2m-1}& \frac{m-1}{2m-1}\\ \frac{m-1}{2m-1}& \frac{m}{2m-1}\end{array}\right].$$

$$\{\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}_1{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}_2{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{in}{\Omega }^c,t\in [0,T],\end{array}$$

$$z{\tilde{G}}_1+a{z}^{1-{\alpha }_1}{\tilde{G}}_1+z^{1-{\alpha }_1}{A}_1{\tilde{G}}_1=a{z}^{1-{\alpha }_2}{\tilde{G}}_2+G_{1,0},$$

$$z{\tilde{G}}_2+a{z}^{1-{\alpha }_2}{\tilde{G}}_2+z^{1-{\alpha }_2}{A}_2{\tilde{G}}_2=a{z}^{1-{\alpha }_1}{\tilde{G}}_1+G_{2,0}.$$

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

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

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

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

$$+aH(z,A_2,{\alpha }_2,{\alpha }_2-{\alpha }_1){\tilde{G}}_1),$$

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

$$+aH(z,A_1,{\alpha }_1,{\alpha }_1-{\alpha }_2){\tilde{G}}_2).$$

$$\Vert H(z,A,\alpha ,\beta ){\Vert }_{L^2(\Omega )\to L^2(\Omega )}\le C|z{|}^{\beta -\alpha },\Vert AH(z,A,\alpha ,\beta ){\Vert }_{L^2(\Omega )\to L^2(\Omega )}\le C|z{|}^{\beta },$$

$$u=(z^{\alpha }+A{)}^{-1}(-au+z^{\beta }v).$$

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

$$\Vert H(z,A,\alpha ,\beta ){\Vert }_{{\hat{H}}^{\sigma }(\Omega )\to {\hat{H}}^{\sigma }(\Omega )}\le C|z{|}^{\beta -\alpha },\Vert AH(z,A,\alpha ,\beta ){\Vert }_{{\hat{H}}^{\sigma }(\Omega )\to {\hat{H}}^{\sigma }(\Omega )}\le C|z{|}^{\beta },$$

$$\Vert AH(z,A,\alpha ,\beta ){\Vert }_{{}_{{\hat{H}}^{\tilde{\sigma }+2\mu s}(\Omega )\to {\hat{H}}^{\tilde{\sigma }}(\Omega )}}\le C|z{|}^{(\beta -\mu \alpha )},$$

$$\Vert u{\Vert }_{{\hat{H}}^{2s}(\Omega )}\le \Vert Au{\Vert }_{L^2(\Omega )}$$

$$\le C|z{|}^{\beta -\alpha }\Vert Av{\Vert }_{L^2(\Omega )}\le C|z{|}^{\beta -\alpha }\Vert v{\Vert }_{{\hat{H}}^{2s}(\Omega )},$$

$$\Vert H(z,A,\alpha ,\beta ){\Vert }_{{\hat{H}}^{2s}(\Omega )\to {\hat{H}}^{2s}(\Omega )}\le C|z{|}^{\beta -\alpha }.$$

$$\Vert H(z,A,\alpha ,\beta ){\Vert }_{{\hat{H}}^{2ms}(\Omega )\to {\hat{H}}^{2ms}(\Omega )}\le C|z{|}^{\beta -\alpha }\mathrm{for}(2m-1)s<1/2,$$

Full text

1 Introduction

Anomalous diffusion phenomena are widespread in the nature world Klafter2005 . Important progresses for modelling these phenomena have been made both microscopically by stochastic processes and macroscopically by partial differential equations (PDEs) Deng2019 . Generally, the PDEs govern the probability density function (PDF) of some particular statistical observables, say, position, functional, first exit time, etc. The fractional Fokker-Planck equation models the PDF of the position of the particles Barkai2001 ; Barkai2000 . So far, there are many numerical methods for solving FFPE, such as finite difference method, finite element method, and even the stochastic methods Deng2009 ; Deng2007 ; Heinsalu2006 ; Meerschaert2006 ; Sun2019 .

Anomalous diffusions with multiple internal states not only are often observed natural phenomena but also some challenge problems, e.g., smart animal searching for food, can be easily treated by taking them as a problem with multiple internal states. Recently, multiple-internal-state Lévy walk and CTRW with independent waiting times and jump lengths are carefully discussed and the PDEs governing the PDF of some statistical observables are derived Xu2018 ; Xu2018-2 . In the CTRW model if the distributions of the jump lengths are power law instead of Gaussian, then the corresponding PDEs involve fractional Laplacian. In this paper, we provide and analyze the numerical scheme for the following fractional Fokker-Planck system (FFPS) with two internal states Xu2018 and the appropriate boundary condition is specified Deng2019 , i.e.,

[TABLE]

where $\Omega$ denotes a bounded convex polygonal domain in $\mathbb{R}^{n}$ $(n=1,2,3)$; $\Omega^{c}$ means the complementary set of $\Omega$ in $\mathbb{R}^{n}$; $\mathbf{M}$ is the transition matrix of a Markov chain, being a $2\times 2$ invertible matrix here; $\mathbf{M}^{T}$ means the transpose of $\mathbf{M}$; $\mathbf{G}=[G_{1},G_{2}]^{T}$ denotes the solution of the system (1); $\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; $\,{}_{0}D^{1-\alpha_{i}}_{t}$ $(i=1,2)$ are the Riemann-Liouville fractional derivatives defined by Podlubny1999

[TABLE]

and $(-\Delta)^{s_{i}}$ $(i=1,2)$ are the fractional Laplacians given as

[TABLE]

where $c_{n,s_{i}}=\frac{2^{2s_{i}}s_{i}\Gamma(n/2+s_{i})}{\pi^{n/2}\Gamma(1-s_{i})}$ and ${\rm P.V.}$ denotes the principal value integral. Without loss of generality, we set $s_{1}\leq s_{2}$ in this paper.

In some sense, the system (1) can be seen as the extension of the model

[TABLE]

It is well known that Eq. (3) has a wide range of practical applications, and there are also some discussions on its regularity and numerical issues Acosta20171 ; Acosta20172 ; Acosta20173 ; Acosta20174 ; in particular, Acosta20174 provides an optimal spatial convergence rates when $s\in(1/2,1)$. Compared with (3), the solutions of the system (1) are coupled with each other and two different space fractional derivatives bring about a huge challenge on the priori estimates of the solutions. Here, we provide a priori estimate for the system (1) with $G_{1}(0),\,G_{2}(0)\in L^{2}(\Omega)$ (see Theorem 2.2) and discuss the regularity of the system (1) detailedly with $s_{1},s_{2}<1/2$ under different regularity assumptions on initial data (see Theorems 2.3 and 2.4). Then we use the finite element method to discretize the fractional Laplacians and provide error analysis for spatial semidiscrete scheme. Lastly, we use $L_{1}$ scheme to discretize the time fractional derivatives and get the first order accuracy without any assumption on the regularity of the solutions. Besides, the proof ideas used in this paper can also be applied to (3) and an optimal spatial convergence rates can be got for $s\in(0,1)$ rather than $s\in(1/2,1)$.

The paper is organized as follows. In Section 2, we first introduce the notations and then focus on the Sobolev regularity of the solutions for the system (1) under different regularity assumptions on initial data. In Section 3, we do the space discretizations by the finite element method and provide the error estimates for the semidiscrete scheme. In Section 4, we use the $L_{1}$ scheme to discretize the time fractional derivatives and provide error estimates for the fully discrete scheme. In Section 5, we confirm the theoretically predicted convergence orders by numerical examples. Finally, we conclude the paper with some discussions. Throughout this paper, $C$ denotes a generic positive constant, whose value may differ at each occurrence, and $\varepsilon>0$ is an arbitrary small constant.

2 Regularity of the solution

In this section, we focus on the regularity of the system (1).

2.1 Preliminaries

Here we make some preparations. Denote $G_{1}(t)$, $G_{2}(t)$ as the functions $G_{1}(\cdot,t)$, $G_{2}(\cdot,t)$ respectively, use the notation ‘$\tilde{~{}}$’ for taking Laplace transform, and introduce $\|\cdot\|_{X\rightarrow Y}$ as the operator norm from $X$ to $Y$, where $X$, $Y$ are Banach spaces. Furthermore, for $\kappa>0$ and $\pi/2<\theta<\pi$, we denote sectors $\Sigma_{\theta}$ and $\Sigma_{\theta,\kappa}$ as

[TABLE]

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

[TABLE]

where the circular arc is oriented counterclockwise, the two rays are oriented with an increasing imaginary part, and $\mathbf{i}$ denotes the imaginary unit. For convenience, in the following we denote $\Gamma_{\theta}=\Gamma_{\theta,0}$ and $A_{i}$ as the fractional Laplacian $(-\Delta)^{s_{i}}$ $(i=1,2)$ with homogeneous Dirichlet boundary condition.

Then we recall some fractional Sobolev spaces Acosta20171 ; Acosta20172 ; Acosta20174 ; Di2012 . Let $\Omega\subset\mathbb{R}^{n}$ $(n=1,2,3)$ be an open set and $s\in(0,1)$. Then the fractional Sobolev space $H^{s}(\Omega)$ can be defined by

[TABLE]

with the norm $\|\cdot\|_{H^{s}(\Omega)}=\|\cdot\|_{L^{2}(\Omega)}+|\cdot|_{H^{s}(\Omega)}$, which constitutes a Hilbert space. As for $s>1$ and $s\notin\mathbb{N}$, the fractional Sobolev space $H^{s}(\Omega)$ can be defined as

[TABLE]

where $\sigma=s-\lfloor s\rfloor$ and $\lfloor s\rfloor$ means the biggest integer not larger than $s$. Another space we use is composed of functions in $H^{s}(\mathbb{R}^{n})$ with support in $\bar{\Omega}$, i.e.,

[TABLE]

whose inner product can be defined as the bilinear form

[TABLE]

Remark 1

According to Acosta20174 , the norm induced by (4) is a multiple of the $H^{s}(\mathbb{R}^{n})$-seminorm, which is equivalent to the full $H^{s}(\mathbb{R}^{n})$-norm on this space because of the fractional Poincaré-type inequality Di2012 . Moreover, from Acosta20171 , $\hat{H}^{s}(\Omega)$ coincides with $H^{s}(\Omega)$ when $s<1/2$.

Next we recall the properties and elliptic regularity of the fractional Laplacian. Reference Acosta20174 claims that $(-\Delta)^{s}:H^{l}(\mathbb{R}^{n})\rightarrow H^{l-2s}(\mathbb{R}^{n})$ is a bounded and invertible operator. Besides, Ref. Grubb2015 proposes the regularity of the following problem

[TABLE]

and the main results are described as

Theorem 2.1 (Grubb2015 )

Let $\Omega\subset\mathbb{R}^{n}$ be a bounded domain with smooth boundary, $g\in H^{r}(\Omega)$ for some $r\geq-s$ and consider $u\in\hat{H}^{s}(\Omega)$ as the solution of the Dirichlet problem (5). Then, there exists a constant $C$ such that

[TABLE]

where $\gamma=\min(s+r,1/2-\epsilon)$ with $\epsilon>0$ arbitrarily small.

2.2 A priori estimate of the solution

According to the property of the transition matrix of a Markov chain Xu2018 , the matrix $\mathbf{M}$ can be denoted as

[TABLE]

and the fact that matrix $\mathbf{M}$ is invertible leads to

[TABLE]

So the system (1) can be rewritten as

[TABLE]

where $a=\frac{1-m}{2m-1}$, $m\in[0,1/2)\cup(1/2,1]$.

Taking the Laplace transforms for the first two equations of the system (6) and using the identity $\widetilde{~{}_{0}D^{\alpha}_{t}u}(z)=z^{\alpha}\tilde{u}(z)$ Podlubny1999 , we have

[TABLE]

Denote

[TABLE]

where $A$ is an operator. Then from (7) and (8) we have

[TABLE]

Thus

[TABLE]

Lemma 1

Let $A$ be the fractional Laplacian $(-\Delta)^{s}$ with homogeneous Dirichlet boundary condition. When $z\in\Sigma_{\theta,\kappa}$, $\pi/2<\theta<\pi$ and $\kappa$ is large enough, we have the estimates

[TABLE]

where $H(z,A,\alpha,\beta)$ is defined in (8).

Proof

Let $u=H(z,A,\alpha,\beta)v$. By simple calculations, we obtain

[TABLE]

Taking $L^{2}$ norm on both sides and using the resolvent estimates provided in Acosta20174 , we have

[TABLE]

which leads to the first desired estimate by taking $\kappa$ large enough and $|z|>\kappa$. Since $AH(z,A,\alpha,\beta)=z^{\beta}(I-(z^{\alpha}+a)H(z,A,\alpha,0))$, it can be easily got the second estimate.

Then we provide the resolvent estimate in $\hat{H}^{1/2+s-\epsilon}(\Omega)$.

Lemma 2

Let $A$ be the fractional Laplacian $(-\Delta)^{s}$ with homogeneous Dirichlet boundary condition and $s<1/2$. When $z\in\Sigma_{\theta,\kappa}$, $\pi/2<\theta<\pi$ and $\kappa$ is large enough, we have the estimates

[TABLE]

where $H(z,A,\alpha,\beta)$ is defined in (8) and $\sigma\in[0,1/2+s)$. Furthermore, there exists

[TABLE]

where $\tilde{\sigma}\in[0,1/2)$ and $\mu\in[0,1]$.

Proof

Assume $u=H(z,A,\alpha,\beta)v$ and $v=0$ in $\Omega^{c}$. Using Theorem 2.1 and Lemma 1, we have

[TABLE]

which leads to

[TABLE]

By induction, one can get

[TABLE]

where $m$ is a positive integer. By Lemma 1 and the interpolation property Adams2003 , there exists

[TABLE]

Taking $\sigma=1/2-s-\epsilon$ in the above equation and using (11), there is

[TABLE]

Similarly by interpolation property, one can obtain

[TABLE]

Noting that $AH(z,A,\alpha,\beta)=z^{\beta}(I-(z^{\alpha}+a)H(z,A,\alpha,0))$, the second estimate can be got. On the other hand, let $u=AH(z,A,\alpha,\beta)v$ and $v=0$ in $\Omega^{c}$. For $\tilde{\sigma}\in[0,1/2)$, we have

[TABLE]

which leads to $\|AH(z,A,\alpha,\beta)\|_{{}_{\hat{H}^{\tilde{\sigma}+2s}(\Omega)\rightarrow\hat{H}^{\tilde{\sigma}}(\Omega)}}\leq C|z|^{\beta-\alpha}$. Using the property of interpolation, we obtain

[TABLE]

Lemma 3

Let $\kappa$ satisfy the conditions given in Lemma 1 and $\Omega\subset\mathbb{R}^{n}$. Then we have the estimate

[TABLE]

Proof

By simple calculation and taking $r=|z|$, we have

[TABLE]

When $\alpha\geq-1$, using the fact $T/t>1$, we can get the desired estimate. And when $\alpha<-1$, the desired estimate can be got by taking $\kappa>1/t$.

Next, we provide the following Grönwall inequality which is similar to the one provided in Laesson1992 .

Lemma 4

Let the function $\phi(t)\geq 0$ be continuous for $0<t\leq T$. If

[TABLE]

for some constants $\{a_{k}\}_{k=0}^{N}$, $\{\alpha_{k}\}_{k=0}^{N}$, $b\geq 0$, $\beta>0$, then there is a constant $C=C(b,T,\alpha,\beta)$ such that

[TABLE]

Proof

The proof is similar to the one provided in Laesson1992 .

Then we present the priori estimates for the solutions $G_{1}$ and $G_{2}$ of (6) with nonsmooth initial value.

Theorem 2.2

Let $\gamma_{1}<\min(s_{1},1/2-\epsilon)$ and $\gamma_{2}<\min(s_{2},1/2-\epsilon)$. If $G_{1,0},\ G_{2,0}\in L^{2}(\Omega)$, then we have

[TABLE]

and

[TABLE]

Proof

By Eq. (10), Lemmas 1, 3 and taking the inverse Laplace transform for (10), we obtain

[TABLE]

According to Lemma 4 and the fact $T/t>1$, one can get the desired $L^{2}$ estimates. Similarly, acting $A_{i}$ on both sides of Eq. (10) respectively and using Lemmas 1, 3, one can obtain

[TABLE]

In view of the $L^{2}$ estimates of $G_{1}$ and $G_{2}$ proved above and $T/t>1$, the desired estimates can be got.

Lastly, we provide a detailed discussion on the regularity of the solutions when $s_{1}<1/2$.

Theorem 2.3

Assume $s_{1}\leq s_{2}<1/2$. If $G_{1,0},\ G_{2,0}\in\hat{H}^{\sigma}(\Omega)$, $\sigma<1/2$, then we have

[TABLE]

where $\gamma_{i}=\min(1/2-\epsilon,s_{i}+\sigma)$ $(i=1,2)$.

Remark 2

The proof of Theorem 2.3 is similar to the one of Theorem 2.2.

Theorem 2.4

Assume $s_{1}\leq s_{2}<1/2$, $G_{i,0}\in\hat{H}^{\sigma_{i}}(\Omega)$, and $\sigma_{i}<1/2$ $(i=1,2)$. Denote $\mu_{1}=\max(\frac{\alpha_{1}-\alpha_{2}}{\alpha_{1}}+\epsilon,0)$, $\mu_{2}=\max(\frac{\alpha_{2}-\alpha_{1}}{\alpha_{2}}+\epsilon,0)$, and $\bar{\gamma}_{i}=\min(1/2-\epsilon,s_{i}+\sigma_{i})$ $(i=1,2)$.

• •

If $\sigma_{1}+2\mu_{1}s_{1}<s_{2}+\bar{\gamma}_{2}$ and $\sigma_{2}+2\mu_{2}s_{2}<s_{1}+\bar{\gamma}_{1}$, then we have

[TABLE]

• •

Assume $\sigma_{1}>\sigma_{2}$. If $\sigma_{1}+2\mu_{1}s_{1}>s_{2}+\bar{\gamma}_{2}$ or $\sigma_{2}+2\mu_{2}s_{2}>s_{1}+\bar{\gamma}_{1}$, then we get

[TABLE]

• •

Assume $\sigma_{1}<\sigma_{2}$. If $\sigma_{1}+2\mu_{1}s_{1}>s_{2}+\bar{\gamma}_{2}$ or $\sigma_{2}+2\mu_{2}s_{2}>s_{1}+\bar{\gamma}_{1}$, then we obtain

[TABLE]

Here $\gamma_{1}=\min(1/2-\epsilon,s_{1}+\sigma_{1},s_{1}+s_{2}+\bar{\gamma}_{2}-2\mu_{1}s_{1})$ and $\gamma_{2}=\min(1/2-\epsilon,s_{2}+\sigma_{2},s_{2}+s_{1}+\bar{\gamma}_{1}-2\mu_{2}s_{2})$.

Proof

When $\sigma_{1}+2\mu_{1}s_{1}\leq s_{2}+\bar{\gamma}_{2}$ and $\sigma_{2}+2\mu_{2}s_{2}\leq s_{1}+\bar{\gamma}_{1}$, acting $A_{i}$ on both sides of (9) respectively, according to Theorem 2.1, embedding theorem, Lemmas 2, 3 and taking the inverse Laplace transform of (9), there are

[TABLE]

where $\gamma_{i}=\min(1/2-\epsilon,s_{i}+\sigma_{i})$ $(i=1,2)$. Thus by Lemma 4, embedding theorem Adams2003 and the fact $T/t>1$, the desired estimates can be got.

If $\sigma_{1}+2\mu_{1}s_{1}>s_{2}+\bar{\gamma}_{2}$ or $\sigma_{2}+2\mu_{2}s_{2}>s_{1}+\bar{\gamma}_{1}$, consider $\sigma_{1}>\sigma_{2}$ first. According to Theorem 2.3, there exists

[TABLE]

Thus

[TABLE]

where $\gamma_{1}=\min(1/2-\epsilon,s_{1}+s_{2}+\bar{\gamma}_{2}-2\mu_{1}s_{1},s_{1}+\sigma_{1})$. Similarly for $\sigma_{1}<\sigma_{2}$, one has

[TABLE]

So

[TABLE]

where $\gamma_{2}=\min(1/2-\epsilon,s_{2}+s_{1}+\bar{\gamma}_{1}-2\mu_{2}s_{2},s_{2}+\sigma_{2})$. Thus by embedding theorem Adams2003 and the fact $T/t>1$, the desired estimates are obtained.

Theorem 2.5

Assume $s_{1}<1/2\leq s_{2}$. If $G_{1,0}\in\hat{H}^{\sigma_{i}}(\Omega)$, $\sigma_{1}<1/2-s_{1}$ and $\sigma_{2}=0$, then

[TABLE]

where $\gamma_{i}=\min(1/2-\epsilon,s_{i}+\sigma_{i})$ $(i=1,2)$.

Remark 3

Combining the proofs of Theorems 2.2 and 2.4, Theorem 2.5 can be obtained.

Remark 4

For Eq. (3), one can obtain that $\|G\|_{\hat{H}^{s+1/2-\epsilon}(\Omega)}\leq Ct^{-\alpha}\|G_{0}\|_{L^{2}(\Omega)}$ for $s\in[1/2,1)$ and $\|G\|_{\hat{H}^{s+\gamma}(\Omega)}\leq Ct^{-\alpha}\|G_{0}\|_{\hat{H}^{\sigma}(\Omega)}$ for $s\in(0,1/2)$, where $\gamma=\min(1/2-\epsilon,s+\sigma)$, $\sigma>0$.

3 Space discretization and error analysis

In this section, we discretize the fractional Laplacian by the finite element method and provide the error estimates for the space semidiscrete scheme of system (6). Let $\mathcal{T}_{h}$ be a shape regular quasi-uniform partitions of the domain $\Omega$, where $h$ is the maximum diameter. Denote $X_{h}$ as the piecewise linear finite element space

[TABLE]

where $\mathcal{P}^{1}$ denotes the set of piecewise polynomials of degree $1$ over $\mathcal{T}_{h}$. Then we define the $L^{2}$-orthogonal projection $P_{h}:\ L^{2}(\Omega)\rightarrow X_{h}$ by

[TABLE]

which has the the following approximation property.

Lemma 5 (Bazhlekova2015 )

The projection $P_{h}$ satisfies

[TABLE]

Denote $(\cdot,\cdot)$ as the $L_{2}$ inner product. The semidiscrete Galerkin scheme for system (6) reads: Find $G_{1,h}\in X_{h}$ and $G_{2,h}\in X_{h}$ such that

[TABLE]

where $v_{1,h},\ v_{2,h}\in X_{h}$. As for $G_{1,h}(0)$ and $G_{2,h}(0)$, we take $G_{1,h}(0)=P_{h}G_{1,0}$, $G_{2,h}(0)=P_{h}G_{2,0}$.

Define the discrete operators $A_{i,h}$: $X_{h}\rightarrow X_{h}$ as

[TABLE]

Then (12) can be rewritten as

[TABLE]

Taking the Laplace transforms of (13), we get

[TABLE]

Next we introduce two lemmas, which will be used in the error estimate between system (6) and space semidiscrete scheme (12).

Lemma 6

For any $\phi\in\hat{H}^{s}(\Omega)$, $z\in\Sigma_{\theta,\kappa}$ with $\theta\in(\pi/2,\pi)$ and $\kappa$ being taken to be large enough to ensure $g(z)=z^{\alpha}+a\in\Sigma_{\theta}$, there exists

[TABLE]

Lemma 7

Let $v\in L^{2}(\Omega)$, $A=(-\Delta)^{s}$ with homogeneous Dirichlet boundary condition, and $z\in\Gamma_{\theta,\kappa}$. Denote $w=H(z,A,\alpha,0)v$ and $w_{h}=H(z,A_{h},\alpha,0)P_{h}v$. Then one has

[TABLE]

where

[TABLE]

and $\gamma\leq\min(s,1/2-\epsilon)$ with $\epsilon>0$ being arbitrarily small.

Remark 5

The proofs of Lemmas 6 and 7 are similar to the ones in Acosta20174 ; Bazhlekova2015 .

When $v\in\hat{H}^{\sigma}(\Omega)$, we modify the estimate in Lemma 7 as

Lemma 8

Let $A=(-\Delta)^{s}$ with homogeneous Dirichlet boundary condition, $s<1/2$, and $z\in\Gamma_{\theta,\kappa}$. Assume $v\in\hat{H}^{\sigma}(\Omega)$ with $\sigma<1/2-s$. Denote $w=H(z,A,\alpha,0)v$ and $w_{h}=H(z,A_{h},\alpha,0)P_{h}v$. Then there exists

[TABLE]

where

[TABLE]

Proof

Let $\epsilon>0$ be arbitrarily small. Here we first consider $\sigma=1/2-s-\epsilon$. Using the notation of $g(z)$ in Lemma 6 and the definitions of $w$ and $w_{h}$, there exist

[TABLE]

Thus

[TABLE]

where $e=w-w_{h}$. By Lemma 6, one has

[TABLE]

Taking $\chi=\pi_{h}w$ as the Lagrange interpolation of $w$ and using the Cauchy-Schwarz inequality, we obtain

[TABLE]

According to Lemma 6, there is

[TABLE]

Thus

[TABLE]

Therefore, $\|(g(z)+A)^{-1}\|_{L^{2}(\Omega)\rightarrow\hat{H}^{s}(\Omega)}\leq C|g(z)|^{-1/2}$. Similar to Lemma 2, there exist

[TABLE]

Using the interpolation property, we get

[TABLE]

Further using the interpolation property leads to

[TABLE]

On the other hand, using Theorem 2.1 and Lemma 2, we obtain

[TABLE]

Thus

[TABLE]

which leads to

[TABLE]

Similarly, for $\phi\in L^{2}(\Omega)$ we set

[TABLE]

By a duality argument, one has

[TABLE]

Then

[TABLE]

where we have used the fact that $|g(z)|^{1/2}\|\psi-\psi_{h}\|_{L^{2}(\Omega)}+\|\psi-\psi_{h}\|_{\hat{H}^{s}(\Omega)}\leq Ch^{s}\|\phi\|_{L^{2}(\Omega)}$Acosta20174 . Thus

[TABLE]

Combining Lemma 7 and using interpolation property, one can get the desired estimate.

For (6), we give the error estimates for the space semidiscrete scheme with nonsmooth initial values.

Theorem 3.1

Let $G_{1}$, $G_{2}$ and $G_{1,h}$, $G_{2,h}$ be the solutions of the systems (6) and (13), respectively, $G_{1,0},\ G_{2,0}\in L^{2}(\Omega)$ and $G_{1,h}(0)=P_{h}G_{1,0}$, $G_{2,h}(0)=P_{h}G_{2,0}$. Then

[TABLE]

where $\gamma_{1}\leq\min(s_{1},1/2-\epsilon)$ and $\gamma_{2}\leq\min(s_{2},1/2-\epsilon)$ with $\epsilon>0$ being arbitrarily small.

Proof

From (14), one can get

[TABLE]

Denote $e_{1}(t)=G_{1}(t)-G_{1,h}(t)$ and $e_{2}(t)=G_{2}(t)-G_{2,h}(t)$. Combining the above equation with (9) leads to

[TABLE]

We first consider $e_{1}$. For $I_{1}$, using inverse Laplace transform and Lemma 7 leads to

[TABLE]

For $I_{2}$, taking inverse Laplace transform, using Eq. (9), Lemma 7, and Theorem 2.2, we have

[TABLE]

where we have used the fact $T/t\leq 1$. As for $I_{3}$, similar to Lemma 1, one has

[TABLE]

Then the inverse Laplace transform and the $L_{2}$ stability of projection $P_{h}$ lead to

[TABLE]

Thus

[TABLE]

Similarly, there also exists

[TABLE]

Thus, the desired estimates can be obtained by Lemma 4 and the fact $T/t>1$.

Finally, combining the proof of Theorem 3.1, the priori estimate provided in Section 2, and Lemma 8, there are the following spatial error estimates for $s_{1}<1/2$. Theorems 3.2, 3.3, and 3.4, are with different assumptions on the regularities of the initial values and/or the range of $s_{2}$.

Theorem 3.2

Let $G_{1}$, $G_{2}$ and $G_{1,h}$, $G_{2,h}$ be the solutions of the systems (6) and (13), respectively. Assume $s_{1}\leq s_{2}<1/2$, $G_{1,0},\ G_{2,0}\in\hat{H}^{\sigma}(\Omega)$, $\sigma<\max(1/2-s_{1},1/2-s_{2})$ and $G_{1,h}(0)=P_{h}G_{1,0}$, $G_{2,h}(0)=P_{h}G_{2,0}$. Then

[TABLE]

where $\gamma_{1}\leq\min(s_{1}+\sigma,1/2-\epsilon)$ and $\gamma_{2}\leq\min(s_{2}+\sigma,1/2-\epsilon)$ with $\epsilon>0$ being arbitrarily small.

Theorem 3.3

Let $G_{1}$, $G_{2}$ and $G_{1,h}$, $G_{2,h}$ be the solutions of the systems (6) and (13), respectively. Assume $s_{1}\leq s_{2}<1/2$, $G_{i,0}\in\hat{H}^{\sigma_{i}}(\Omega)$, $\sigma_{i}<1/2-\epsilon$ $(i=1,2)$ and $G_{1,h}(0)=P_{h}G_{1,0}$, $G_{2,h}(0)=P_{h}G_{2,0}$. Denote $\mu_{1}=\max(\frac{\alpha_{1}-\alpha_{2}}{\alpha_{1}}+\epsilon,0)$, $\mu_{2}=\max(\frac{\alpha_{2}-\alpha_{1}}{\alpha_{2}}+\epsilon,0)$ and $\bar{\gamma}_{i}=\min(1/2-\epsilon,s_{i}+\sigma_{i})$ $(i=1,2)$.

• •

If $\sigma_{1}+2\mu_{1}s_{1}\leq 2s_{2}+\sigma_{2}$ and $\sigma_{2}+2\mu_{2}s_{2}\leq 2s_{1}+\sigma_{1}$, then

[TABLE]

• •

Assume $\sigma_{1}>\sigma_{2}$. If $\sigma_{1}+2\mu_{1}s_{1}>2s_{2}+\sigma_{2}$ or $\sigma_{2}+2\mu_{2}s_{2}>2s_{1}+\sigma_{1}$, then

[TABLE]

• •

Assume $\sigma_{1}<\sigma_{2}$. If $\sigma_{1}+2\mu_{1}s_{1}>2s_{2}+\sigma_{2}$ or $\sigma_{2}+2\mu_{2}s_{2}>2s_{1}+\sigma_{1}$, then

[TABLE]

Here $\gamma_{1}=\min(s_{1}+\sigma_{1},s_{1}+s_{2}+\bar{\gamma}_{2}-2\mu_{1}s_{1},1/2-\epsilon)$ and $\gamma_{2}=\min(s_{2}+\sigma_{2},s_{2}+s_{1}+\bar{\gamma}_{1}-2\mu_{2}s_{2},1/2-\epsilon)$ with $\epsilon>0$ being arbitrarily small.

Theorem 3.4

Let $G_{1}$, $G_{2}$ and $G_{1,h}$, $G_{2,h}$ be the solutions of the systems (6) and (13), respectively. Assume $s_{1}<1/2\leq s_{2}$, $G_{i,0}\in\hat{H}^{\sigma_{i}}(\Omega)$ $(i=1,2)$, $\sigma_{1}<1/2-s_{1}$, $\sigma_{2}=0$ and $G_{1,h}(0)=P_{h}G_{1,0}$, $G_{2,h}(0)=P_{h}G_{2,0}$. Then

[TABLE]

where $\gamma_{i}\leq\min(s_{i}+\sigma_{i},1/2-\epsilon)$, $i=1,2$ with $\epsilon>0$ arbitrary small.

Remark 6

From the numerical experiments, we find that the errors aroused by $aH(z,A_{1},\alpha_{1},\alpha_{1}-\alpha_{2})P_{h}\tilde{G}_{2}-aH(z,A_{1,h},\alpha_{1},\alpha_{1}-\alpha_{2})\tilde{G}_{2,h}$ and $aH(z,A_{2},\alpha_{2},\alpha_{2}-\alpha_{1})P_{h}\tilde{G}_{1}-aH(z,A_{2,h},\alpha_{2},\alpha_{2}-\alpha_{1})\tilde{G}_{1,h}$ in (15) have almost no effect on convergence rates.

Remark 7

As for Eq. (3), the spatial semidiscrete scheme can be written as

[TABLE]

where $v_{h}\in X_{h}$ and $G_{h}(0)=P_{h}G_{0}$. According to Lemma 8, if $s<1/2$ and $G_{0}\in\hat{H}^{\sigma}(\Omega)$, the error between $G(t)$ and $G_{h}(t)$ can be written as

[TABLE]

where $\gamma=\min(1/2-\epsilon,s+\sigma)$. And according to Lemma 7, if $s\geq 1/2$ and $G_{0}\in L^{2}(\Omega)$, the error between $G(t)$ and $G_{h}(t)$ is as follows

[TABLE]

4 Time discretization and error analysis

In this section, we use the $L_{1}$ scheme to discretize the Riemann-Liouville time fractional derivatives and perform the error analysis for the fully discrete scheme. We first introduce the notations as

[TABLE]

Lemma 9 (Nie2018 )

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

[TABLE]

According to (16), the solution of (13) in Laplace space can be reconstructed as

[TABLE]

Next, we use the Backward Euler scheme to discretize $\partial/\partial t$ and $L_{1}$ scheme to approximate $~{}_{0}D^{\alpha}_{t}$. Let the time step size $\tau=T/L$, $L\in\mathbb{N}$, $t_{i}=i\tau$, $i=0,1,\ldots,L$ and $0=t_{0}<t_{1}<\cdots<t_{L}=T$. Recall the approximation of Caputo fractional derivative by $L_{1}$ scheme (see, e.g., Jin2015 )

[TABLE]

where

[TABLE]

Using the relationship between the Caputo fractional derivative and the Riemann-Liouville fractional derivative, i.e.,

[TABLE]

we obtain

[TABLE]

where

[TABLE]

For the system (6), we have the fully discrete scheme

[TABLE]

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

To get the error estimate between (6) and (19), we introduce $Li_{p}(z)$ Flajolet1999 ; Jin2015 defined by

[TABLE]

and recall the Lemmas about $Li_{p}(z)$ .

Lemma 10 (Flajolet1999 ; Jin2015 )

For $p\neq 1,2,\cdots$, the function $Li_{p}(e^{-z})$ satisfies the singular expansion

[TABLE]

where $\varsigma(z)$ denotes the Riemann zeta function.

Lemma 11 (Flajolet1999 ; Jin2015 )

Let $|z|\leq\frac{\pi}{\sin(\theta)}$ with $\theta\in(\frac{\pi}{2},\frac{5\pi}{6})$ and $-1<p<0$. Then

[TABLE]

converges absolutely.

At the same time, we have the estimates.

Lemma 12 (Jin2015 ; Jin2016 )

Assume $z\in\Sigma_{\theta}$, $|z|\leq\frac{\pi}{\tau\sin(\theta)}$ and $\theta\in(\pi/2,\pi)$. Then there are

[TABLE]

Next, we give the error estimates of the fully discrete scheme. To get the solutions of the system (19), multiplying $\zeta^{n}$ on both sides of the first two equations in (19), summing $n$ from $1$ to $\infty$ and using (18), there exist

[TABLE]

[TABLE]

where

[TABLE]

As for $\psi^{\alpha}(\zeta)$, using the definition of $d^{\alpha}_{j}$ and $Li_{p}(z)$, we have

[TABLE]

Then, there is the following estimate.

Lemma 13

Let $z\in\Gamma_{\theta,\kappa}$, $|z\tau|\leq\frac{\pi}{sin(\theta)}$ and $\theta\in(\pi/2,5\pi/6)$. Then we have

[TABLE]

Proof

By Lemma 11, there exists

[TABLE]

Now we give the error estimates between the solutions of the systems (13) and (19).

Theorem 4.1

Let $G_{1,h}$, $G_{2,h}$ and $G^{n}_{1,h}$, $G^{n}_{2,h}$ be the solutions of the systems (13) and (19), respectively. Then

[TABLE]

Proof

We first consider the error estimates between $G^{n}_{1,h}$ and $G_{1,h}$. By (20), for small $\xi_{\tau}=e^{-\tau(\kappa+1)}$, there is

[TABLE]

Letting $\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\}\cup\{z\in\mathbb{C}:|z|=\kappa,|\arg z|\leq\theta\}$. Thus

[TABLE]

In view of (17), there exists

[TABLE]

Combining (22) and (23) leads to

[TABLE]

According to Lemma 9, there exists

[TABLE]

For $II$, similarly it has

[TABLE]

Next for $III$ and $IV$, there are

[TABLE]

and

[TABLE]

As for $III_{1}$ and $IV_{1}$, using Lemmas 9 and 12 leads to

[TABLE]

and

[TABLE]

Thus

[TABLE]

and

[TABLE]

Denote $H^{(a,b)}(z_{1},z_{2},A_{1},A_{2})$ as the $a$-th order derivative of $z_{1}$ and $b$-th order derivative of $z_{2}$. Using

[TABLE]

the mean value theorem, and the Lemmas 12 and 13, there are

[TABLE]

and

[TABLE]

Thus

[TABLE]

and

[TABLE]

By simple calculations, we have

[TABLE]

[TABLE]

In summary,

[TABLE]

Analogously, it has

[TABLE]

The proof has been completed.

5 Numerical experiments

In this section, we perform the numerical experiments to verify the effectiveness of the designed schemes. Since the exact solutions $G_{1}$ and $G_{2}$ are unknown, to get the spatial convergence rates, we calculate

[TABLE]

where $G^{n}_{1,h}$ and $G^{n}_{2,h}$ mean the numerical solutions of $G_{1}$ and $G_{2}$ at time $t_{n}$ with mesh size $h$; similarly, to obtain the temporal convergence rates, we calcuate

[TABLE]

where $G_{1,\tau}$ and $G_{2,\tau}$ are the numerical solutions of $G_{1}$ and $G_{2}$ at the fixed time $t$ with step size $\tau$. Then the spatial and temporal convergence rates can be, respectively, obtained by

[TABLE]

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)$.

Here we first give some examples to show the influence of the regularity of initial data on convergence rates.

Example 1

We take $a=2$, $\tau=1/800$, and $T=1$ to solve the system (6) with the initial condition (a), and $s_{1}=s_{2}<1/2$, $\alpha_{1}=0.4$, $\alpha_{2}=0.7$. Here $G_{1,0},G_{2,0}\in\hat{H}^{1/2-\epsilon}(\Omega)$ satisfy the conditions of Theorem 3.1. Table 1 shows that the convergence rates can be achieved as $O(h^{s+1/2-\epsilon})$, which agree with Theorems 3.2 and 3.3.

Example 2

We take $\alpha_{1}=0.4$, $\alpha_{2}=0.6$, $a=-2$, $\tau=1/800$, and $T=1$ to solve the system (6) with the initial condition (a). Table 2 shows the $L_{2}$ errors and convergence rates for different values of $s_{1},\,s_{2}$. The convergence rates are consistent with the results of Theorem 3.1 when $s_{1},s_{2}>1/2$; when $s_{1},s_{2}<1/2$, the convergence rates of $G_{2}$ are higher than the predicted ones in Theorem 3.2 (or Theorem 3.3) and the convergence rates of $G_{1}$ are the same as the predicted ones, the reason of which may be the less effect of $aH(z,A_{1},\alpha_{1},\alpha_{1}-\alpha_{2})P_{h}\tilde{G}_{2}-aH(z,A_{1,h},\alpha_{1},\alpha_{1}-\alpha_{2})\tilde{G}_{2,h}$ and $aH(z,A_{2},\alpha_{2},\alpha_{2}-\alpha_{1})P_{h}\tilde{G}_{1}-aH(z,A_{2,h},\alpha_{2},\alpha_{2}-\alpha_{1})\tilde{G}_{1,h}$ in (15) on convergence rates.

Example 3

The parameters are taken as $\alpha_{1}=0.8$, $\alpha_{2}=0.9$, $a=2$, $\tau=1/800$, and $T=1$. First, we solve the system (6) with the initial condition (b). Letting $\nu_{1}=\nu_{2}=0.4999$ leads to $G_{1,0},G_{2,0}\in L^{2}(\Omega)$. According to Table 3, the convergence rates agree with Theorem 3.1 when $s_{1},\,s_{2}>1/2$; when $s_{1},\,s_{2}<1/2$, the convergence rates of $G_{2}$ are higher than the predicted ones in Theorem 3.1 and the convergence rates of $G_{1}$ are the same as the predicted ones, the reason of which is the same as that stated in Example 2.

Then we take $\nu_{1}=-0.4$ and $\nu_{2}=-0.3$, which may lead to $G_{1,0}\in\hat{H}^{0.1}(\Omega)$ and $G_{1,0}\in\hat{H}^{0.2}(\Omega)$. Table 4 shows the convergence results and we find the convergence rates of $G_{2}$ are higher than the predicted ones in Theorem 3.1 and the convergence rates of $G_{1}$ are the same as the predicted ones, and the reason for these phenomena is the same as the one in Example 2.

Example 4

In this example, we take $\alpha_{1}=0.7$, $\alpha_{2}=0.6$, $a=0.1$, $\tau=1/50$, and $T=20$. The system (6) is solved with the initial condition (b) and we take $\nu_{1}=0$, $\nu_{2}=0.4999$, which implies $G_{1,0}\in\hat{H}^{1/2-\epsilon}(\Omega)$, $G_{2,0}\in L^{2}(\Omega)$. According to Table 5, the results for $s_{1}=0.25$ and $s_{2}=0.8$ agree with Theorem 3.4; when $s_{1},\,s_{2}<1/2$, the convergence rates of $G_{1}$ are higher than the predicted ones in Theorem 3.3 and the convergence rates of $G_{2}$ are the same as the predicted ones, the reason of which is the same as that stated in Example 2.

Finally, we verify the temporal convergence rates in the following example.

Example 5

Here we take $s_{1}=0.25$, $s_{2}=0.75$, $a=2$, and $h=1/400$ to solve the system (6) with the initial condition (a). Table 6 shows the $L_{2}$ errors and convergence rates for different $\alpha_{1},\,\alpha_{2}$, which can be used to validate the results of Theorem 4.1.

6 Conclusion

The power law distributions are widely observed in heterogeneous media, relating to the fields of physics, biology, and social science, etc. This paper focuses on the regularity and numerical methods of the two state model with fractional Laplacians, characterizing the power law properties. The priori estimates are obtained under various different regularity assumptions of initial values and/or different powers of fractional Laplacians. The designed numerical scheme is with finite element approximation for fractional Laplacians and $L_{1}$ scheme to discretize the time fractional Riemann-Liouville derivative. For the scheme, the complete error analyses are provided, and the extensive numerical experiments are performed to validate their effectiveness.

Acknowledgements.

We thank Buyang Li for the discussions.