Local discontinuous Galerkin methods for the time tempered fractional diffusion equation

TL;DR

This paper develops and analyzes local discontinuous Galerkin methods for solving time-tempered fractional diffusion equations, proving stability, convergence, and validating results with numerical experiments.

Contribution

It introduces a semi-discrete LDG scheme for tempered fractional diffusion and extends it to fully discrete schemes with error estimates, a novel approach for this class of equations.

Findings

Semi-discrete LDG scheme is unconditionally stable in L2 norm.

Optimal convergence rate of O(h^{k+1}) is achieved.

Numerical experiments confirm theoretical error estimates.

Abstract

In this article, we consider discrete schemes for a fractional diffusion equation involving a tempered fractional derivative in time. We present a semi-discrete scheme by using the local discontinuous Galerkin (LDG) discretization in the spatial variables. We prove that the semi-discrete scheme is unconditionally stable in $L^2$ norm and convergence with optimal convergence rate $\mathcal{O}(h^{k+1})$. We develop a class of fully discrete LDG schemes by combining the weighted and shifted Lubich difference operators with respect to the time variable, and establish the error estimates. Finally, numerical experiments are presented to verify the theoretical results.

Tables

Table 1. Table 1: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different λ 𝜆 \lambda at T = 1 𝑇 1 T=1 with α = 0.5 𝛼 0.5 \alpha=0.5 , τ = h 2 / 3 𝜏 superscript ℎ 2 3 \tau=h^{2/3} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =2$
$k$	$h$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$1$	1/10	1.757029e-02		8.266055e-03		4.823736e-03
	1/20	4.472764e-03	1.9739	2.091684e-03	1.9825	1.215416e-03	1.9887
	1/40	9.995738e-04	2.1618	4.874263e-04	2.1014	2.918227e-04	2.0583
	1/80	2.487744e-04	2.0065	1.215344e-04	2.0038	7.285848e-05	2.0019

Table 2. Table 2: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different λ 𝜆 \lambda at T = 1 𝑇 1 T=1 with α = 0.5 𝛼 0.5 \alpha=0.5 , τ = h 3 / 2 𝜏 superscript ℎ 3 2 \tau=h^{3/2} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =2$
$k$	$h$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$2$	1/10	7.817372e-04		3.849712e-04		2.321083e-04
	1/20	9.719738e-05	3.0077	4.805674e-05	3.0019	2.905722e-05	2.9978
	1/40	1.203610e-05	3.0135	5.976656e-06	3.0073	3.624898e-06	3.0029
	1/80	1.506415e-06	2.9982	7.477283e-07	2.9988	4.533749e-07	2.9992

Table 3. Table 3: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different λ 𝜆 \lambda at T = 1 𝑇 1 T=1 with α = 0.5 𝛼 0.5 \alpha=0.5 , τ = h 4 / 3 𝜏 superscript ℎ 4 3 \tau=h^{4/3} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =2$
$k$	$h$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$3$	1/10	3.174252e-05		1.553122e-05		9.321072e-06
	1/20	1.971099e-06	4.0093	9.698783e-07	4.0012	5.844179e-07	3.9954
	1/40	1.211423e-07	4.0242	6.009536e-08	4.0125	3.642279e-08	4.0041
	1/80	7.570034e-09	4.0003	3.756951e-09	3.9996	2.277746e-09	3.9992

Table 4. Table 4: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different λ 𝜆 \lambda at T = 1 𝑇 1 T=1 with α = 0.5 𝛼 0.5 \alpha=0.5 , τ = h 5 / 4 𝜏 superscript ℎ 5 4 \tau=h^{5/4} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =2$
$k$	$h$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$4$	1/5	4.435466e-05		2.097285e-05		1.228325e-05
	1/10	1.380329e-06	5.0060	6.796152e-07	4.9477	4.096977e-07	4.9060
	1/20	4.402038e-08	4.9707	2.160624e-08	4.9752	1.299606e-08	4.9784
	1/40	1.374369e-09	5.0013	6.805691e-10	4.9886	4.119552e-10	4.9794

Table 5. Table 5: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different λ 𝜆 \lambda at T = 1 𝑇 1 T=1 with α = 0.5 𝛼 0.5 \alpha=0.5 , τ = h 3 / 2 𝜏 superscript ℎ 3 2 \tau=h^{3/2} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =2$
$k$	$h$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$5$	1/5	6.414671e-06		3.026083e-06		1.769351e-06
	1/10	1.061861e-07	5.9167	5.229200e-08	5.8547	3.152810e-08	5.8104
	1/20	1.689493e-09	5.9739	8.353260e-10	5.9681	5.050749e-10	5.9640
	1/40	2.658552e-11	5.9898	1.320140e-11	5.9836	8.006770e-12	5.9791

Table 6. Table 6: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different α 𝛼 \alpha , λ 𝜆 \lambda at T = 1 𝑇 1 T=1 with τ = h 3 / 2 𝜏 superscript ℎ 3 2 \tau=h^{3/2} , P 5 superscript 𝑃 5 P^{5} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =2$
$\alpha$	$h$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$0.3$	1/10	3.949105e-06		1.862967e-06		1.089277e-06
	1/20	6.459789e-08	5.9339	3.181162e-08	5.8719	1.917998e-08	5.8276
	1/40	1.030441e-09	5.9702	5.094751e-10	5.9644	3.080510e-10	5.9603
$0.6$	1/10	7.681143e-06		3.623533e-06		2.118680e-06
	1/20	1.271142e-07	5.9171	6.259814e-08	5.8551	3.774191e-08	5.8109
	1/40	2.017666e-09	5.9773	9.975830e-10	5.9715	6.031825e-10	5.9674
$0.9$	1/10	1.101658e-05		5.197007e-06		3.038690e-06
	1/20	1.808903e-07	5.9284	8.908050e-08	5.8664	5.370876e-08	5.8221
	1/40	2.847323e-09	5.9894	1.407785e-09	5.9836	8.512089e-10	5.9795

Table 7. Table 7: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different λ 𝜆 \lambda with fixed h = 1 / 500 ℎ 1 500 h=1/500 and P 2 superscript 𝑃 2 P^{2} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =3$
$q$	$\tau$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$1$	1/5	6.268259e-05		3.254490e-05		8.428462e-06
	1/10	2.943336e-05	1.0906	1.563827e-05	1.0574	4.264206e-06	0.9830
	1/20	1.418827e-05	1.0528	7.640964e-06	1.0333	2.143689e-06	0.9922
	1/40	6.955695e-06	1.0284	3.773311e-06	1.0179	1.074573e-06	0.9963

Table 8. Table 8: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different λ 𝜆 \lambda with fixed h = 1 / 300 ℎ 1 300 h=1/300 and P 3 superscript 𝑃 3 P^{3} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =3$
$q$	$\tau$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$2$	1/5	4.081267e-05		2.158503e-05		5.934015e-06
	1/10	1.131338e-05	1.8510	6.164704e-06	1.8079	1.817960e-06	1.7067
	1/20	2.989321e-06	1.9201	1.657132e-06	1.8953	5.083344e-07	1.8385
	1/40	7.690827e-07	1.9586	4.302745e-07	1.9454	1.347830e-07	1.9151

Table 9. Table 9: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different λ 𝜆 \lambda with fixed h = 1 / 200 ℎ 1 200 h=1/200 and P 4 superscript 𝑃 4 P^{4} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =3$
$q$	$\tau$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$3$	1/5	2.470358e-05		1.349055e-05		4.084818e-06
	1/10	3.798010e-06	2.7014	2.168691e-06	2.6371	7.345770e-07	2.4753
	1/20	5.249281e-07	2.8551	3.072493e-07	2.8193	1.108525e-07	2.7283
	1/40	6.894422e-08	2.9286	4.087384e-08	2.9102	1.523835e-08	2.8629

Table 10. Table 10: The L 2 superscript 𝐿 2 L^{2} errors and convergence orders for different λ 𝜆 \lambda with fixed h = 1 / 100 ℎ 1 100 h=1/100 and P 5 superscript 𝑃 5 P^{5} .

		$\lambda =0.8$		$\lambda =1.5$		$\lambda =3$
$q$	$\tau$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$4$	1/5	1.842158e-05		1.074646e-05		3.850820e-06
	1/10	1.452152e-06	3.6651	9.187289e-07	3.5481	4.024233e-07	3.2584
	1/20	1.010676e-07	3.8448	6.674299e-08	3.7830	3.267246e-08	3.6226
	1/40	6.652413e-09	3.9253	4.489239e-09	3.8941	2.325680e-09	3.8124

Table 11. Table 11: The L 2 superscript 𝐿 2 L^{2} errors and orders of accuracy of problem ( 5.2 ) calculated by the fully LDG schemes ( 4.5 ) for different λ 𝜆 \lambda with α = 0.5 𝛼 0.5 \alpha=0.5 , τ = h 2 𝜏 superscript ℎ 2 \tau=h^{2} , and T = 1 𝑇 1 T=1 .

		$\lambda =0.5$		$\lambda =2.5$		$\lambda =5$
$k$	$h$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$1$	1/5	5.969030e-04		8.078203e-05		6.630993e-06
	1/10	1.490511e-04	2.0017	2.017187e-05	2.0017	1.655808e-06	2.0017
	1/20	3.725150e-05	2.0004	5.041443e-06	2.0004	4.138268e-07	2.0004

Table 12. Table 12: The L 2 superscript 𝐿 2 L^{2} errors and orders of accuracy of problem ( 5.2 ) calculated by the fully LDG schemes ( 4.5 ) for different λ 𝜆 \lambda with α = 0.5 𝛼 0.5 \alpha=0.5 , τ = h 3 𝜏 superscript ℎ 3 \tau=h^{3} , and T = 1 𝑇 1 T=1 .

		$\alpha =0.1$		$\alpha =0.5$		$\alpha =0.9$
$k$	$h$	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order	$\parallel \cdot \parallel$-error	order
$2$	1/5	2.100809e-05		1.359967e-05		6.184385e-06
	1/10	2.669664e-06	2.9762	1.723894e-06	2.9798	7.878049e-07	2.9727
	1/20	3.351008e-07	2.9940	2.163795e-07	2.9940	1.080580e-07	2.8660

Equations

$$u_t(x,t)={\kappa }_{\alpha }{}_0{D}_t^{1-\alpha ,\lambda }\left(u_{xx}(x,t)\right)-\lambda u(x,t),$$

$${}_0{D}_t^{\gamma ,\lambda }u(t)={{}_{0}D}_t^{n,\lambda }{}_0{I}_t^{n-\gamma ,\lambda }u(t),$$

$${}_0{I}_t^{\sigma ,\lambda }u(t)=\frac{1}{\Gamma (\sigma )}{\int }_0^t{e}^{-\lambda (t-s)}(t-s{)}^{\sigma -1}u(s)ds,\sigma =n-\alpha ,$$

$$_{0}D_{t}^{n,\lambda}=\bigg{(}\frac{d}{dt}+\lambda\bigg{)}^{n}=\underbrace{\bigg{(}\frac{d}{dt}+\lambda\bigg{)}\cdots\bigg{(}\frac{d}{dt}+\lambda\bigg{)}}_{n\ times}.$$

$$\frac{d}{dt}\left(e^{\lambda t}u(x,t)\right)={\kappa }_{\alpha }{}_0{D}_t^{1-\alpha }\left(e^{\lambda t}{u}_{xx}(x,t)\right).$$

$${{}_{0}^{C}D}_{t}^{\alpha}\big{(}e^{\lambda t}u(x,t)\big{)}=\kappa_{\alpha}e^{\lambda t}u_{xx}(x,t),$$

$${{}_0^{C}D}_t^{\alpha }u(t)=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{1}{(t-s{)}^{\alpha }}\frac{du(s)}{ds}ds,0<\alpha <1.$$

$${{}_{0}^{C}D}_{t}^{\alpha,\lambda}u(t)=e^{-\lambda t}~{}{{}_{0}^{C}D}_{t}^{\alpha}\big{(}e^{\lambda t}u(t)\big{)}=\frac{e^{-\lambda t}}{\Gamma(1-\alpha)}\int_{0}^{t}\frac{1}{(t-s)^{\alpha}}\frac{d(e^{\lambda s}u(s))}{ds}ds,$$

$${{}_0^{C}D}_t^{\alpha ,\lambda }u(x,t)={\kappa }_{\alpha }{u}_{xx}(x,t).$$

$$u(x,0)=u_0(x),x\in \Omega$$

$$u(x,t)=u(x+L,t),L>0,t\in (0,T].$$

$$\frac{1}{2}{{}_0^{C}D}_t^{\alpha }(u^2(t))\le u(t){{}_0^{C}D}_t^{\alpha }(u(t)),0<\alpha <1.$$

$$\Vert u(\cdot ,t){\Vert }^2+2{\kappa }_{\alpha }{{}_{0}I}_t^{\alpha ,2\lambda }\Vert u_x(\cdot ,t){\Vert }^2\le e^{-2\lambda t}\Vert u_0{\Vert }^2,$$

$${{}_0^{C}D}_t^{\gamma }v(x,t)={\kappa }_{\alpha }{v}_{xx}(x,t),$$

$$\big{(}v(\cdot,t),{{}_{0}^{C}D}_{t}^{\gamma}v(\cdot,t)\big{)}_{\Omega}=-\kappa_{\alpha}\big{(}v_{x}(x,t),v_{x}(x,t)\big{)}_{\Omega}+\kappa_{\alpha}v_{x}(x,t)v(x,t)\big{|}^{L}_{0},$$

$${{}_0^{C}D}_t^{\alpha }\Vert v(\cdot ,t){\Vert }^2+2{\kappa }_{\alpha }\Vert v_x(\cdot ,t){\Vert }^2\le 0.$$

$${{}_{0}I}_t^{\alpha }{{}_{0}D}_t^{\alpha }(w(t))=w(t)-w(0),$$

$$\Vert v(\cdot ,t){\Vert }^2+2{\kappa }_{\alpha }{{}_{0}I}_t^{\alpha }\Vert v_x(\cdot ,t){\Vert }^2\le \Vert v_0{\Vert }^2.$$

$$0=x_{1/2}<x_{3/2}<\cdots <x_{N-1/2}<x_{N+1/2}=L.$$

$$V_{h}=\big{\{}v:v\mid_{I_{j}}\in P^{k}(I_{j}),j=1,2,\cdots,N\big{\}},$$

$$\begin{array}[]{l}\displaystyle{{}_{0}^{C}\!D}_{t}^{\alpha,\lambda}u(x,t)-\kappa_{\alpha}p_{x}=0,\\ \displaystyle p-u_{x}=0.\end{array}$$

$$\begin{array}[]{l}\displaystyle\big{(}{{}_{0}^{C}\!D}_{t}^{\alpha,\lambda}u_{h},v_{h}\big{)}_{I_{j}}+\kappa_{\alpha}\left(\big{(}p_{h},(v_{h})_{x}\big{)}_{I_{j}}-(\widehat{p}_{h}v_{h}^{-})_{j+\frac{1}{2}}+(\widehat{p}_{h}v_{h}^{+})_{j-\frac{1}{2}}\right)=0,\\ \displaystyle\big{(}p_{h},w_{h}\big{)}_{I_{j}}+\big{(}u_{h},(w_{h})_{x}\big{)}_{I_{j}}-(\widehat{u}_{h}w_{h}^{-})_{j+\frac{1}{2}}+(\widehat{u}_{h}w_{h}^{+})_{j-\frac{1}{2}}=0.\end{array}$$

$${\hat{u}}_h=u_h^{-},{\hat{p}}_h=p_h^{+}.$$

$$e^{-\lambda t}\frac{1}{2}\;_{0}^{C}\!D_{t}^{\alpha,\lambda}\big{(}u(t)\big{)}^{2}\leq u(t)~{}_{0}^{C}\!D_{t}^{\alpha,\lambda}\big{(}u(t)\big{)},0<\alpha<1,\lambda>0.$$

$$\frac{1}{2}\;_{0}^{C}\!D_{t}^{\alpha}\big{(}e^{\lambda t}u(t)\big{)}^{2}\leq e^{\lambda t}u(t)~{}_{0}^{C}\!D_{t}^{\alpha}\big{(}e^{\lambda t}u(t)\big{)}.$$

$$e^{-2\lambda t}\frac{1}{2}\;_{0}^{C}\!D_{t}^{\alpha}\big{(}e^{\lambda t}u(t)\big{)}^{2}\leq u(t)~{}_{0}^{C}\!D_{t}^{\alpha,\lambda}\big{(}u(t)\big{)}.$$

$$\big{(}u_{h},\,{{}_{0}^{C}\!D}_{t}^{\alpha,\lambda}u_{h}\big{)}_{\Omega}+\kappa_{\alpha}\|p_{h}\|^{2}=0.$$

$$B(u_h,p_h;v_h,w_h)$$

$$B(u_h,p_h;u_h,{\kappa }_{\alpha }{p}_h)$$

$$\|u_{h}(\cdot,t)\|^{2}+2\kappa_{\alpha}~{}{{}_{0}I}_{t}^{\alpha,\lambda}\big{(}e^{\lambda t}\|p_{h}(\cdot,t)\|^{2}\big{)}\leq e^{-\lambda t}\|u_{h}(\cdot,0)\|^{2}.$$

Full text

fourierlargesymbols147

1 Introduction

In this paper we discuss a local discontinuous Galerkin method to solve the following time tempered fractional subdiffusion equation [1, 2]

[TABLE]

where $u(x,t)$ represents the probability density of finding a particle on $x$ at time $t$, $\kappa_{\alpha}>0$ is the diffusion coefficient, and ${}_{0}D_{t}^{1-\alpha,\lambda}(0<\alpha<1)$ denotes the Riemann-Liouville tempered fractional derivative operator. The Riemann-Liouville tempered fractional derivative of order $\gamma~{}(n-1<\gamma<n)$ is defined by [3, 4, 5]

[TABLE]

where ${{}_{0}I_{t}^{n-\alpha,\lambda}}$ denotes the Riemann-Liouville fractional tempered integral [3, 4, 5]

[TABLE]

and

[TABLE]

Tempered fractional calculus can be recognized as the generalization of fractional calculus. If we taking $\lambda=0$ in (1.2), then the tempered fractional integral and derivative operators reduce to the Riemann-Liouville fractional integral ${{}_{0}I_{t}^{\sigma}}$ and derivative ${}_{0}D_{t}^{\gamma}$ operators, respectively.

In recent years, many numerical methods such as finite difference methods [6, 7, 8, 9, 10, 11], finite element methods [12, 13, 14, 15, 16] and spectral methods [17, 18] have been developed for the numerical solutions of fractional subdiffusion and superdiffusion equations. Limited works are reported for solving the tempered fractional differential equations, when compared with a large volume of literature on numerical solutions of fractional differential equations. In this literature, Baeumera and Meerschaert [19] provide finite difference and particle tracking methods for solving the tempered fractional diffusion equation with the second order accuracy. The stability and convergence of the provided schemes are discussed. Cartea and del-Castillo-Negrete [20] derive a finite difference scheme to numerically solve a Black-Merton-Scholes model with tempered fractional derivatives. Hanert and Piret [22] presented a Chebyshev pseudo-spectral scheme to solve the space-time tempered fractional diffusion equation, and proved that the method yields an exponential convergence rate. Zayernouri et al. [23] derived an efficient Petrov-Galerkin method for solving tempered fractional ODEs by using the eigenfunctions of tempered fractional Sturm-Liouville problems. By using the weighted and shifted Gr$\ddot{\rm u}$nwald difference (WSGD) operators, Li and Deng [5] designed a series of high order numerical schemes for the space tempered fractional diffusion equations. This technique is used to solve the tempered fractional Black-Scholes equation for European double barrier option [21]. Using the properties of generalized Laguerre functions, Huang et al. [24] used Laguerre functions to approximate the substantial fractional ODEs on the half line. Li et al. [25] analysed the well-posedness and developed the Jacobi-predictor-corrector algorithm for the tempered fractional ordinary differential equation. Yu et al. [26] developed the third and fourth order quasi-compact approximations for one and two dimensional space tempered fractional diffusion equations. By using the weighted and shifted Grünwald-Letnikov formula suggested in [29], Hao et al. [27] constructed a second-order approximation for the time tempered fractional diffusion equation. By introducing fractional integral spaces, Zhao et al. [28] discussed spectral Galerkin and Petrov-Galerkin methods for tempered fractional advection problems and tempered fractional diffusion problems.

Recently, the high order and fast numerical methods for fractional differential equations draw the wide interests of the researchers [29, 30, 31, 32, 33]. The local discontinuous Galerkin (LDG) method is one of the most popular methods in this literature. The LDG method was first introduced to solve a convection-diffusion problems by Cockburn and Shu [34]. These methods have recently become increasingly popular due to their flexibility for adaptive simulations, suitability for parallel computations, applicability to problems with discontinuous solutions, and compatibility with other numerical methods. Nowadays, the LDG method has been successfully used in solving linear and nonlinear elliptic, parabolic, hyperbolic equations, and some mixed schemes. For the recent development of discontinuous Galerkin methods, see the the monograph and review articles [35, 36, 37, 38] and the references therein. More recently, some researchers pay attention to solving the fractional partial equations by the LDG method. For the time fractional differential equations, Mustapha and McLean [39] employed a piecewise-linear, discontinuous Galerkin method for the time discretization of a sub-diffusion equation. A LDG method for space discretization of a time fractional diffusion equation is discussed in Xu and Zheng’s work [40]. By using L1 time discretization, Wei et al.[41] developed an implicit fully discrete LDG finite element method for solving the time-fractional Schröinger equation, Guo et al. [42] studied a LDG method for some time fractional fourth-order differential equations. Liu et al. [43] proposed LDG method combined with a third order weighted and shifted Grünwald difference operators for a fractional subdiffusion equation. For the space fractional differential equations, based on two (or four) auxiliary variables in one dimension (or two dimensions) and the Caputo derivative as the spatial derivative, Ji and Tang [44] developed the high-order accurate Runge-Kutta LDG methods for one- and two-dimensional space-fractional diffusion equations with the variable diffusive coefficients. Deng and Hesthaven [45] have developed a LDG method for space fractional diffusion equation and given a fundamental frame to combine the LDG methods with fractional operators.

In this paper, we will develop and analyze a new class of LDG method for the model (1.1). Our new method is based on a combination of the weighted and shifted Lubich difference approaches in the time direction and a LDG method in the space direction. Stability and convergence of semi-discrete and fully discrete LDG schemes are rigorously analyzed. We show that the fully discrete scheme is unconditionally stable with convergence order of $\mathcal{O}(\tau^{q}+h^{k+1}),q=1,2,3,4,5$.

The rest of the article is organized as follows. In section 2, we first consider the initial boundary value problem of the tempered fractional diffusion equation. Then in section 3 we construct a semi-discrete LDG method for the considered equation. We perform the detailed theoretical analysis for the stability and error estimate of the semi-discrete numerical scheme in this section. In section 4, we apply the weighted and shifted Lubich difference approximation for the temporal discretization of the time fractional tempered equation. The error estimates are provided for the full-discrete LDG scheme. Finally, some numerical examples and physical simulations are presented in section 5 which confirm the theoretical statement. Some concluding remarks are given in the final section.

2 Initial boundary value problem of the tempered fractional diffusion equation

Instead of designing the numerical scheme of the equation (1.1) directly, we constructing the numerical scheme for its equivalent form. By simple calculation, the equation (1.1) can be rewritten as

[TABLE]

Performing Riemann-Liouville fractional integral ${}_{0}I_{t}^{1-\alpha}$ on both side of (1.1), we arrive at

[TABLE]

where ${{}_{0}^{C}D}_{t}^{\alpha}$ denotes the Caputo fractional derivative [3]

[TABLE]

In view of the Caputo tempered fractional derivative [4, 5]

[TABLE]

we get the following time tempered fractional diffusion equation

[TABLE]

Let $\Omega=(0,L)$ be the space domain. We now consider the initial boundary value problem of the fractional tempered diffusion equation (2.5) in the domain $(0,L)\times(0,T]$, subject to the initial condition

[TABLE]

and the boundary conditions

[TABLE]

Lemma 2.1**.**

[46]** For function $u(t)$absolutely continuous on $[0,t]$, holds the the following inequality

[TABLE]

For the initial-boundary value problem (2.5)-(2.7), we have the following stability.

Theorem 2.1**.**

The solution $u(x,t)$ of the initial-boundary value problem (2.5)-(2.7) holds the prior estimate

[TABLE]

where ${{}_{0}I}_{t}^{\alpha,2\lambda}$ denotes the Riemann-Liouville fractional tempered integral operator defined by (1.3).

Proof.

Taking $v(x,t)=e^{\lambda t}u(x,t)$ in (2.5), we have

[TABLE]

Taking inner product in equation (2.10) with the space variable, we get

[TABLE]

where $(u(x),v(x))_{\Omega}=\int_{\Omega}u(x)v(x)dx$. Denoting the $L^{2}$-norm as $\|u\|=(u(x),u(x))_{\Omega}^{1/2}$, using the boundary conditions (2.7) and using the inequality (2.8), we get

[TABLE]

By applying the fractional integral operator ${{}_{0}I}_{t}^{\alpha}$ to both sides of inequality (2.12), using the composite properties of fractional calculus [3]

[TABLE]

we obtain

[TABLE]

Taking $u(x,t)=e^{-\lambda t}v(x,t)$ in (2.13) we get the desired estimate (2.9). ∎

3 The semi-discrete LDG scheme

In this section, we present and analyze a local discontinuous Galerkin method for the equation (2.5) subjects to the initial condition (2.6) and the periodic boundary conditions (2.7). For the interval $\Omega=[0,L]$, we divide it into $N$ cells as follows

[TABLE]

We denote $I_{j}=(x_{j-1/2},x_{j+1/2}),x_{j}=(x_{j-1/2}+x_{j+1/2})/2,$ and $h_{j}=x_{j+1/2}-x_{j-1/2},h={\rm max}_{1\leq j\leq N}h_{j}.$ Furthermore, we define the mesh $\mathcal{T}=\{I_{j}=(x_{j-1/2},x_{j+1/2}),j=1,2,...,N\}$. The finite element space is defined by

[TABLE]

where $P^{k}(I_{j})$ denotes the set of all polynomials of degree at most $k$ on cell $I_{j}$. We define $u^{-}_{j+1/2},u^{+}_{j+1/2}$ represent the values of $u$ at $x_{j+1/2}$ from the left cell $I_{j}$ and the right cell $I_{j+1}$, respectively.

To define the local discontinuous Galerkin method, we rewrite (2.5) as a first-order system

[TABLE]

Now we can define the local discontinuous Galerkin method to the system (3.1). Find $u_{h},p_{h}$, for all $v_{h},w_{h}\in V_{h}$ such that

[TABLE]

The ’hat’ terms in (3.2) are the numerical probability density fluxes. We chose the alternating numerical fluxes [36, 48] as

[TABLE]

3.1 Stability analysis of the semi-discrete LDG scheme

In this section, we present the stability and convergence analysis for the semi-discrete scheme (3.2) in $L^{2}$ sense. To do so, we follow the technique used by Cockburn and Shu [34].

Lemma 3.1**.**

For function $u(t)$absolutely continuous on $[0,t]$, holds the the following inequality

[TABLE]

Proof.

Taking $v(t)=e^{\lambda t}u(t)$ in the inequality (2.8), we have

[TABLE]

which leads

[TABLE]

∎

Lemma 3.2**.**

Let $u_{h},p_{h}$ to be the solution of semi-discrete LDG scheme (3.2) with the flux $\widehat{u}_{h}$ ,$\widehat{p}_{h}$ defined in (3.3), holds

[TABLE]

Proof.

For simplicity, we denote

[TABLE]

If we take $v_{h}=u_{h},w_{h}=\kappa_{\alpha}p_{h}$ in (3.1), we have

[TABLE]

Combining the numerical flux defined by (3.3) and periodic boundary conditions (2.7) we arrive at (3.7). ∎

Corollary 3.1**.**

($L^{2}$-stability) The semi-discrete LDG scheme (3.2) with the flux choice (3.3) is $L^{2}$-stable, i.e.

[TABLE]

Proof.

Using the inequality (3.4), we can obtain

[TABLE]

Recall the composite properties of the Caputo tempered fractional derivative and the Riemann-Liouville tempered fractional integral, we have

[TABLE]

Applying the operator ${}_{0}I_{t}^{\alpha,\lambda}$ on both sides of the inequality (3.11) leads to (3.10). ∎

3.2 Convergence analysis of the semi-discrete LDG scheme

Now, we given the $L^{2}$ error estimate. In order to give more detailed error estimate, the following two special projections operators introduced in [36] will be used.

[TABLE]

[TABLE]

Lemma 3.3**.**

[34]** For projection operators $\mathcal{P^{\pm}}$, the following estimate holds

[TABLE]

where $w^{e}=\mathcal{P^{\pm}}w-w,$ $C$ is a positive constant depending $u$ and its derivatives but independent of $h$.

Lemma 3.4**.**

[47]** Let $u(t)$ be continuous and non-negative on $[0,T]$. If

[TABLE]

where $0\leq\alpha<1$. $\varphi(t)$ is nonnegative monotonic increasing continuous function on $[0,T]$, and $M$ is a positive constant, then

[TABLE]

where $E_{1-\alpha}(z)$ denotes the Mittag-Leffler function defined for all $0\leq\alpha<1$ by [3]

[TABLE]

Lemma 3.5**.**

Let $u_{h},p_{h}$ to be the solution of semi-discrete LDG scheme (3.2), and $u,p$ be the exact solution of (2.5) with initial condition (2.6) and the periodic boundary (2.7), the following error estimate holds

[TABLE]

Proof.

With the denote in (3.1), we directly get

[TABLE]

and

[TABLE]

Subtracting (3.21) from (3.20), then we obtain the error equation

[TABLE]

where we denote $e=u-u_{h},\overline{e}=p-p_{h}$. We divide the error both $e$ and $\overline{e}$ into two parts

[TABLE]

If we take $v_{h}=e_{h},w_{h}=\kappa_{\alpha}\overline{e}_{h}$ in (3.22), we get

[TABLE]

For the left side of (3.24), using the equation (3.7) in Lemma 3.2, we have

[TABLE]

Obviously, the right of (3.24) can be written as

[TABLE]

Since $(e_{h})_{x}$ and $(\overline{e}_{h})_{x}$ are polynomials of degree at most $k-1$, applying the properties (3.13) and (3.14) of the projections $\mathcal{P^{\pm}}$, we obtain

[TABLE]

In other way,

[TABLE]

By the Cauchy’s inequality, we get

[TABLE]

From the Lemma 3.3, we conclude that

[TABLE]

Using the inequality (3.4), we have

[TABLE]

Combining the composite properties (3.12) and the definition of Riemann-Liouville tempered fractional integral, we arrive at

[TABLE]

where we used the fact

[TABLE]

Furthermore, using the fractional Gronwall’s lemma 3.4, we have

[TABLE]

where $E_{\alpha}(\cdot)$ denotes the Mittag-Leffler function is defined by (3.18). ∎

4 Fully discrete LDG schemes

In this section we discrete the time in the semi-discrete scheme by virtue of high order approximation. Let $0=t_{0}<t_{1}<\cdots<t_{n}<t_{n+1}<\cdots<t_{M}=T$ be the subdivision of the time interval $[0,T]$, with the time step $\tau=t_{n+1}-t_{n}$. To achieve the high order accuracy, we employ the $q$-th order approximations given in [30] to approximate the Riemann-Liouville tempered derivative, i.e.

[TABLE]

where $R^{n}=\mathcal{O}(\tau^{q})$ and

[TABLE]

More details of ${l}_{k}^{q,\alpha}$, one can refer to [30]. Using (4.1), we find

[TABLE]

Recalling the relation of Riemann-Liouville and Caputo tempered fractional derivatives [4, 5]

[TABLE]

the weak form of the first order system (3.1) at $t_{n}$ can be rewritten as

[TABLE]

Let $u_{h}^{n},p_{h}^{n}\in V_{h}$ be the approximate solution of $u(x,t_{n}),p(x,t_{n})$, respectively. We propose the fully discrete LDG schemes as follows: Find $u_{h}^{n},p_{h}^{n}\in V_{h}$,

[TABLE]

for all $v,w\in V_{h}$, $j=1,2,...,N$. We take numerical flux to be $\widehat{u}_{h}^{n}=(u_{h}^{n})^{-},\widehat{q}_{h}^{n}=(q_{h}^{n})^{+}$ with the same choice of (3.3). In the following, we prove the stability and error estimate of the schemes (4.5) with $q=1$ in $L_{2}$ norm. For convenience, we denote ${l}_{k}^{1,\alpha}$ by $w_{k}$, where the coefficients

[TABLE]

Lemma 4.1**.**

[9]** The coefficients $w_{k}$ defined in (4.6) satisfy

[TABLE]

and

[TABLE]

Theorem 4.1**.**

The fully discrete LDG scheme (4.5) of initial-boundary problem (2.5)-(2.7) is unconditional stability and holds

[TABLE]

Proof.

Setting $v=u_{h}^{n},w=\kappa_{\alpha}\tau^{\alpha}p_{h}^{n}$ in (4.5) and summing over all elements, we obtain

[TABLE]

where

[TABLE]

and

[TABLE]

Recalling the numerical flux in (3.3), we have $\Theta(u_{h}^{n},p_{h}^{n})=0$. On the other hand, in view of the periodic boundary conditions (2.7), we get

[TABLE]

Using the Cauchy-Schwartz inequality, we arrive at

[TABLE]

Therefore, we have

[TABLE]

Next we need to prove the following estimate by mathematical induction

[TABLE]

From the inequality (4.13), we can see the inequality (4.14) holds obviously when $n=1$. Assuming

[TABLE]

then from the inequality (4.13), we obtain

[TABLE]

The proof is complete. ∎

Theorem 4.2**.**

Let $u(x,t_{n})$ be the exact solution of the problem (2.5)-(2.7), which is sufficiently smooth such that $u\in H^{m+1}$ with $0\leq m\leq k+1$. Let $u_{h}^{n}$ be the numerical solution of the fully discrete LDG scheme (4.5), then there holds the following error estimate

[TABLE]

where $C$ is a constant depending on $u,T,\alpha$ but independent of $\tau$ and $h$.

Proof.

To simplify the notation, we decompose the errors as follows:

[TABLE]

Combining (4.4) and (4.5), we have

[TABLE]

Substituting (4.16) into (4.17) and notice that $e_{u}^{0}=0$, we can get the error equation

[TABLE]

Taking $v=\mathcal{P^{-}}e_{u}^{n}$, $w=\kappa_{\alpha}\tau^{\alpha}\mathcal{P^{+}}e_{p}^{n}$ in (4.18), we get

[TABLE]

Using the properties of projections $\mathcal{P^{\pm}}$ and $w_{0}=1$, we can further get

[TABLE]

Applying the Cauchy-Schwarz inequality, we have

[TABLE]

Combining the inequality $2ab\leq a^{2}+b^{2}$ and $e^{-\lambda k\tau}\in(0,1]$, from the above inequality (4.21), we can derive

[TABLE]

Moreover, we have

[TABLE]

Next, we prove the following estimate by mathematical introduction

[TABLE]

For $n=1$, using the properties (3.15) of the projections $\mathcal{P^{\pm}}$, it can be seen that the inequality (4.24) holds obviously. Assuming

[TABLE]

Remembering $-\sum_{k=1}^{n-1}w_{k}<1$ and the properties of the projections $\mathcal{P^{\pm}}$, we have

[TABLE]

Finally, combining the triangle inequality and lemma 3.3 to have

[TABLE]

∎

5 Numerical experiments

In this section, we perform three examples to illustrate the effectiveness of our numerical schemes and confirm our theoretical results.

Example 5.1**.**

Without loss of generality, we add a force term $f(x,t)$ on the right hand side of the equation (2.5), we consider

[TABLE]

with periodic boundary conditions $u(x+\pi,t)=u(x,t)$ and initial condition $u(x,0)=sin(2\pi x)$. If we take the force term $f(x,t)$ as f(x,t)=e^{-\lambda t}\bigg{(}\frac{\Gamma(\beta+1)}{\Gamma(\beta+1-\alpha)}t^{\beta-\alpha}+4\kappa_{\alpha}\pi^{2}(t^{\beta}+1)\bigg{)}sin(2\pi x), then the exact solution of the problem (5.1) with the corresponding initial-boundary condition gives

[TABLE]

The $L^{2}$ errors and orders of the fully discrete LDG scheme (4.5) on uniform meshes are present in Table 1-Table 10. Table 1-Table 5 list the $L^{2}$ errors and orders of accuracy for schemes (4.5) with different $k$ and fixed $\alpha=0.5$. In these tests we take $\tau=h^{(k+1)/q}$. All the numerical results given in Table 1-Table 5 are consistent with the theoretical analysis which presented in theorem 4.2. Table 6 shows the errors and orders of scheme (4.5) for solving the problem (5.1) with the different parameters $\alpha$ and $\lambda$. As expected, we observe that our scheme can achieve higher order accuracy in space, as well as in time. To test the high order of scheme (4.5) in time direction, we list the errors and orders of scheme (4.5) in Table 7-Table 10. We can again clearly observe the desired orders of accuracy from these tables.

Example 5.2**.**

In this example, we examine the following homogeneous equation

[TABLE]

subjects to the boundary conditions $u(0,t)=0,u(1,t)=0$, and the initial value $u(x,0)=\sin(2\pi x).$ We can check that the exact solution of this initial-boundary value problem (5.2) is $u(x,t)=e^{-\lambda t}E_{\alpha}(-~{}4\pi^{2}t^{\alpha})\sin(2\pi x),$ where the generalized Mittag-Leffler function $E_{\alpha}(\cdot)$ defined in (3.18).

In this test, the finite element space is piecewise linear and piecewise quadratic polynomials for the second and third order schemes, respectively. The numerical results are shown in Table 11-Table 12. The evolution of numerical solutions with different $\alpha,\lambda$ at different times are given in Fig. 1.

Example 5.3**.**

In this example, we will test the dynamics behavior of the tempered fractional diffusion equation (5.2) with homogeneous Dirichlet boundary conditions on a finite domain $[-4,4]$. We take the Gaussian function

[TABLE]

as the initial condition.

The numerical results for this example are calculated by the fully discrete scheme (4.5). In the computation, we set $h=1/40,\tau=h^{2},\sigma=0.01.$ The probability density function of a diffusion particle for different values of $\alpha,\lambda$ at different times are given in Fig. 2. It can be seen that, the different parameters $\alpha,\lambda,t$ have different effect for the probability density of a particle, which is in agreement with the analytic results given in [1, 2]. The effectiveness of our numerical schemes is confirmed once again.

6 Conclusions

We have presented a numerical method for a time fractional tempered diffusion equation. The proposed method is based on a combination of the weighted and shifted Lubich difference approaches in the time direction and a LDG method in the space direction. The convergence rate of the method is proven by providing a priori error estimate, and confirmed by a series of numerical tests. It has been proved that the proposed scheme is unconditionally stable and of $q$-order convergence in time and $k+1$-order convergence in space. Some numerical experiments have been carried out to support the theoretical results.

Acknowledgments

This research was partially supported by the National Natural Science Foundation of China under Grant No.11426174, the Starting Research Fund from the Xi an university of Technology under Grant Nos. 108-211206, 2014CX022, the Natural Science Basic Research Plan in Shaanxi Province of China under Grant No.2015JQ1022, the Shaanxi science and technology research projects under Grant No.2015GY004.

Reference