Discontinuous Galerkin methods and their adaptivity for the tempered fractional (convection) diffusion equations

TL;DR

This paper develops and analyzes adaptive discontinuous Galerkin methods for solving tempered fractional convection-diffusion equations, demonstrating their stability, convergence, and effectiveness through theoretical analysis and numerical experiments.

Contribution

It introduces a novel adaptive DG framework with a posteriori error estimates for tempered fractional PDEs, applicable to general fractional operators.

Findings

Proved stability and convergence of the proposed DG schemes.

Designed effective local error indicators for adaptivity.

Validated the methods through extensive numerical experiments.

Abstract

This paper focuses on the adaptive discontinuous Galerkin (DG) methods for the tempered fractional (convection) diffusion equations. The DG schemes with interior penalty for the diffusion term and numerical flux for the convection term are used to solve the equations, and the detailed stability and convergence analyses are provided. Based on the derived posteriori error estimates, the local error indicator is designed. The theoretical results and the effectiveness of the adaptive DG methods are respectively verified and displayed by the extensive numerical experiments. The strategy of designing adaptive schemes presented in this paper works for the general PDEs with fractional operators.

Tables

Table 1. Table 1. Numerical errors ( L 2 subscript 𝐿 2 L_{2} ) and orders of convergence on unstructured meshes for Example 2.19 .

	K	68	211		436		702
N	$(\alpha ,\beta )$	error	error	order	error	order	error	order
1	(0.2,0.2)	4.46e-2	1.52e-2	1.90	7.60e-3	1.91	5.20e-3	1.59
	(0.5,0.5)	4.23e-2	1.42e-2	1.93	7.00e-3	1.95	4.80e-3	1.58
	(0.7,0.2)	4.29e-2	1.45e-2	1.92	7.20e-3	1.93	4.90e-3	1.62
	K	68	211		436		702
2	(0.2,0.2)	1.01e-2	2.10e-3	2.77	7.75e-4	2.75	3.80e-4	2.99
	(0.5,0.5)	9.70e-3	2.00e-3	2.79	7.19e-4	2.82	3.51e-4	3.01
	(0.7,0.2)	9.80e-3	2.10e-3	2.72	7.33e-4	2.90	3.57e-4	3.02
	K	68	211		436		702
3	(0.2,0.2)	3.00e-3	3.94e-4	3.59	1.02e-4	3.72	3.76e-5	4.21
	(0.5,0.5)	3.00e-3	3.87e-4	3.62	1.01e-4	3.71	3.70e-5	4.20
	(0.7,0.2)	3.00e-3	3.89e-4	3.61	1.01e-4	3.71	3.71e-5	4.21

Table 2. Table 2. Numerical errors ( L 2 subscript 𝐿 2 L_{2} ) and orders of convergence on unstructured meshes for Example 2.20 .

	K	68	211		436		702
N	$(\alpha ,\beta )$	error	error	order	error	order	error	order
1	(0.5,0.5)	3.33e-2	9.94e-3	2.13	4.66e-3	2.09	2.86e-3	2.05
	(1.5,1.5)	3.09e-2	9.30e-3	2.12	4.50e-3	2.00	2.80e-3	1.99
	(1.5,0.5)	3.20e-2	9.60e-3	2.12	4.50e-3	2.09	2.80e-3	1.99
	K	68	211		436		702
2	(0.5,0.5)	1.50e-3	3.09e-4	2.79	1.08e-4	2.90	5.40e-5	2.91
	(1.5,1.5)	9.73e-4	2.00e-4	2.74	7.36e-5	2.75	3.49e-5	3.13
	(1.5,0.5)	1.03e-3	1.90e-4	2.99	6.15e-5	3.11	2.65e-5	3.53
	K	68	211		436		702
3	(0.5,0.5)	3.72e-4	4.75e-5	3.81	1.24e-5	3.70	4.62e-6	4.15
	(1.5,1.5)	3.29e-4	4.31e-5	3.59	1.14e-5	3.66	4.22e-6	4.17
	(1.5,0.5)	3.51e-4	4.52e-5	3.62	1.18e-5	3.70	4.40e-6	4.14

Table 3. Table 3. Example 3.8 : Uniform refinement.

K	${\Vert u-u_k\Vert }_{L^2}$	${\Vert u-u_k\Vert }_E$
32	0.3164	1.1565
72	0.1908	0.9010
128	0.1257	0.6794
200	0.0935	0.5264
288	0.0755	0.4357
392	0.0637	0.3803
512	0.0528	0.3362

Table 4. Table 4. Example 3.8 : Refinement based on the energy-norm indicator.

k	K	${\Vert u-u_k\Vert }_{L^2}$	${\Vert u-u_k\Vert }_E$	$\eta$	$I_{eff}$
1	8	0.6752	1.2608	3.5025	2.7781
2	15	0.3930	1.0819	2.0108	1.8586
3	28	0.2614	0.7705	1.6311	2.1168
4	42	0.1652	0.6573	1.1697	1.7796
5	86	0.1169	0.4383	0.9193	2.0973
6	152	0.0872	0.4074	0.6599	1.6198
7	201	0.0684	0.3023	0.5550	1.8361

Table 5. Table 5. Example 3.8 : Refinement based on the weighted indicator.

k	K	${\Vert u-u_k\Vert }_{L^2}$	${\Vert u-u_k\Vert }_E$	$\eta$	$I_{eff}$
1	8	0.6752	1.2608	1.0679	0.8470
2	15	0.3930	1.0819	0.7608	0.7032
3	28	0.2614	0.7705	0.3620	0.4698
4	36	0.1650	0.6408	0.2558	0.3992
5	63	0.1242	0.4852	0.2260	0.4658
6	110	0.0933	0.4174	0.1603	0.3840
7	166	0.0551	0.2766	0.1146	0.4143

Equations

$$\left\{\begin{array}[]{ll}\partial_{t}u+\textbf{b}\cdot\nabla u-\kappa_{1}\nabla_{x}^{\alpha,\lambda}u-\kappa_{2}\nabla_{y}^{\beta,\lambda}u=f,&(\textbf{x},t)\in\Omega\times J,\\ u(\textbf{x},0)=u_{0}(\textbf{x}),&\textbf{x}\in\Omega,\\ u(\textbf{x},t)=0,&(\textbf{x},t)\in\mathbb{R}^{2}\backslash\Omega\times J,\end{array}\right.$$

$${}_{-\infty }{I}_x^{\alpha ,\lambda }u(x)={{\mathrm{e}}^{-\lambda x}}_{-\infty }{I}_x^{\alpha }[{\mathrm{e}}^{\lambda x}u(x)]=\frac{1}{\Gamma (\alpha )}{\int }_{-\infty }^x(x-\xi {)}^{\alpha -1}{\mathrm{e}}^{-\lambda (x-\xi )}u(\xi )d\xi ,$$

$${}_x{I}_{\infty }^{\alpha ,\lambda }u(x)={{\mathrm{e}}^{\lambda x}}_x{I}_{\infty }^{\alpha }[{\mathrm{e}}^{-\lambda x}u(x)]=\frac{1}{\Gamma (\alpha )}{\int }_x^{\infty }(\xi -x{)}^{\alpha -1}{\mathrm{e}}^{-\lambda (\xi -x)}u(\xi )d\xi .$$

$${}_{-\infty }{D}_x^{\alpha ,\lambda }u(x)={{\mathrm{e}}^{-\lambda x}}_{-\infty }{D}_x^{\alpha }[{\mathrm{e}}^{\lambda x}u(x)]={(\lambda +D{)}^n}_{-\infty }{I}_x^{n-\alpha ,\lambda }u(x),$$

$${}_x{D}_{\infty }^{\alpha ,\lambda }u(x)={{\mathrm{e}}^{\lambda x}}_x{D}_{\infty }^{\alpha }[{\mathrm{e}}^{-\lambda x}u(x)]={(\lambda -D{)}^n}_x{I}_{\infty }^{n-\alpha ,\lambda }u(x).$$

$${}_{-}{}_{\infty }^C{D}_x^{\alpha ,\lambda }u(x)={\mathrm{e}}^{-\lambda x}{}_{-}{}_{\infty }^C{D}_x^{\alpha }[{\mathrm{e}}^{\lambda x}u(x)]={}_{-\infty }{I}_x^{n-\alpha ,\lambda }(\lambda +D{)}^{n}u(x),$$

$${}_x^C{D}_{\infty }^{\alpha ,\lambda }u(x)={\mathrm{e}}^{\lambda x}{}_x^C{D}_{\infty }^{\alpha }[{\mathrm{e}}^{-\lambda x}u(x)]={}_x{I}_{\infty }^{n-\alpha ,\lambda }(\lambda -D{)}^{n}u(x).$$

$${\nabla }_x^{\alpha ,\lambda }u(\text{x},t)=-{\kappa }_{\alpha }{[}_{-\infty }{\nabla }_x^{\alpha ,\lambda }+{}_x{\nabla }_{\infty }^{\alpha ,\lambda }]u(\text{x},t),$$

$${\nabla }_y^{\beta ,\lambda }u(\text{x},t)=-{\kappa }_{\beta }{[}_{-\infty }{\nabla }_y^{\beta ,\lambda }+{}_y{\nabla }_{\infty }^{\beta ,\lambda }]u(\text{x},t),$$

$${}_{-\infty }{\nabla }_x^{\alpha ,\lambda }u(\text{x},t)={}_{-\infty }{D}_x^{\alpha ,\lambda }u(\text{x},t)-{\lambda }^{\alpha }u(\text{x},t)\mathrm{for}\alpha \in (0,1),$$

$${}_{-\infty }{\nabla }_x^{\alpha ,\lambda }u(\text{x},t)={}_{-\infty }{D}_x^{\alpha ,\lambda }u(\text{x},t)-{\lambda }^{\alpha }u(\text{x},t)-\alpha {\lambda }^{\alpha -1}\frac{\partial u(\text{x},t)}{\partial x}\mathrm{for}\alpha \in (1,2),$$

$${}_x{\nabla }_{\infty }^{\alpha ,\lambda }u(\text{x},t)={}_x{D}_{\infty }^{\alpha ,\lambda }u(\text{x},t)-{\lambda }^{\alpha }u(\text{x},t)\mathrm{for}\alpha \in (0,1),$$

$${}_x{\nabla }_{\infty }^{\alpha ,\lambda }u(\text{x},t)={}_x{D}_{\infty }^{\alpha ,\lambda }u(\text{x},t)-{\lambda }^{\alpha }u(\text{x},t)+\alpha {\lambda }^{\alpha -1}\frac{\partial u(\text{x},t)}{\partial x}\mathrm{for}\alpha \in (1,2),$$

$${}_a{D}_x^{\mu ,\lambda }f(x)={}_a^C{D}_x^{\mu ,\lambda }f(x)+\sum _{j=0}^{m-1}{D}^{j,\lambda }f(x){|}_{x=a}\frac{(x-a{)}^{j-\mu }{\mathrm{e}}^{-\lambda (x-a)}}{\Gamma (j-\mu +1)},$$

$${}_a{I}_x^{\mu ,\lambda }[{}_a{D}_x^{\mu ,\lambda }f(x)]=f(x)-\sum _{j=1}^m{}_a{D}_x^{\mu -j,\lambda }f(x){|}_{x=a}\frac{(x-a{)}^{\mu -j}{\mathrm{e}}^{-\lambda (x-a)}}{\Gamma (\mu -j+1)},$$

$${}_a{D}_x^{\nu ,\lambda }[{}_a{D}_x^{\mu ,\lambda }f(x)]={}_a{D}_x^{\nu +\mu ,\lambda }f(x)-\sum _{j=1}^m{}_a{D}_x^{\mu -j,\lambda }f(x){|}_{x=a}\frac{(x-a{)}^{-\nu -j}{\mathrm{e}}^{-\lambda (x-a)}}{\Gamma (-\nu -j+1)},$$

$$f^{(j)}(a)=f^{(j)}(b)=0\mathrm{for}j=0,1,\cdots ,m-1.$$

$${}_a{D}_x^{\mu -j,\lambda }f(x){|}_{x=a}={}_x{D}_b^{\mu -j,\lambda }f(x){|}_{x=b}=0\mathrm{for}j=1,2,\cdots ,m,$$

$${}_a{D}_x^{\mu ,\lambda }f(x)={}_a^C{D}_x^{\mu ,\lambda }f(x),$$

$${}_a{I}_x^{\mu ,\lambda }[{}_a{D}_x^{\mu ,\lambda }f(x)]=f(x),$$

$${}_a{D}_x^{\nu ,\lambda }[{}_a{D}_x^{\mu ,\lambda }f(x)]={}_a{D}_x^{\nu +\mu ,\lambda }f(x).$$

$${\int }_a^b{}_a{I}_x^{\alpha ,\lambda }u(x)\cdot v(x)dx={\int }_a^{b}u(x)\cdot {}_x{I}_b^{\alpha ,\lambda }v(x)dx.$$

$${\int }_a^b{}_a{D}_x^{\alpha ,\lambda }u(x)\cdot v(x)dx={\int }_a^{b}u(x)\cdot {}_x{D}_b^{\alpha ,\lambda }v(x)dx.$$

$$\mathcal{F}{[}_{-\infty }{I}_x^{\mu ,\lambda }u(x)](\omega )=(\lambda +i\omega {)}^{-\mu }\hat{u}(\omega ),$$

$$\mathcal{F}{[}_x{I}_{\infty }^{\mu ,\lambda }u(x)](\omega )=(\lambda -i\omega {)}^{-\mu }\hat{u}(\omega ).$$

$$\mathcal{F}{[}_{-\infty }{D}_x^{\mu ,\lambda }u(x)](\omega )=(\lambda +i\omega {)}^{\mu }\hat{u}(\omega ),$$

$$\mathcal{F}{[}_x{D}_{\infty }^{\mu ,\lambda }u(x)](\omega )=(\lambda -i\omega {)}^{\mu }\hat{u}(\omega ),$$

$$(p+q{)}^{\alpha }\sim p^{\alpha }+q^{\alpha }.$$

$$2^{\alpha -1}(p^{\alpha }+q^{\alpha })\le (p+q{)}^{\alpha }\le (p^{\alpha }+q^{\alpha })\mathrm{for}0<\alpha \le 1$$

$$(p^{\alpha }+q^{\alpha })\le (p+q{)}^{\alpha }\le 2^{\alpha -1}(p^{\alpha }+q^{\alpha })\mathrm{for}\alpha >1.$$

Full text

1. Introduction

Fractional calculus [6] is a popular mathematical tool for modeling anomalous diffusions [22], being ubiquitous in nature. Microscopically, anomalous diffusion can be described by continuous time random walk (CTRW), governed by the waiting time and jump length; generally the first moment of the waiting time and/or the second moment of the jump length diverge(s). Sometimes, it is better to temper the broad distribution(s) of the waiting time and/or the jump length [2, 12, 18], because of the boundedness of physical space or the finite lifespan of the biological particles or the slow transition of different diffusion types. Based on the tempered CTRW, the partial differential equations (PDEs) characterizing the evolution of the functional distribution of the trajectories of the particles are derived [33], which reduce to the PDEs describing the distribution of the positions of the particles if taking the parameter $p$ over there as [math], called tempered fractional PDEs; here, we discuss their (adaptive) discontinuous Galerkin (DG) methods.

There are already some progresses for numerically solving (tempered) fractional PDEs by variational methods [12, 15, 17, 20, 23, 26, 32, 34, 36]. Ervin and Roop [15] firstly present the variational formulation for the fractional advection dispersion equation. The DG methods are particularly extended to fractional problems with their majority of characteristics [5, 10, 16, 27, 31, 35], naturally being formulated for any order of accuracy in any element, being flexible in choosing element sizes in any place, suitable for adaptivity, being local and easy to invert for mass matrix, leading to an explicit formulation for time dependent problems, etc. Cockburn and Mustapha [11] provide a hybridizable DG method for fractional diffusion problems; Mclean and Mustapha [23] discuss the superconvergence of the DG method for the fractional diffusion and wave equations; Xu and Hesthaven [34], and Wang et al [32], respectively, consider DG and hybridized DG methods for the fractional convection-diffusion equations; Zayernouri and Karniadakis [36] design discontinuous spectral element methods for the time and space fractional differential equations. Du et al [13] give a convergent adaptive finite element algorithm for nonlocal diffusion and peridynamic models. It seems that there are not works for digging out the potential advantages of DG methods in adaptivity for fractional problems, by deriving posteriori error estimates and providing the local error indicators.

The model we consider in this paper is the two dimensional space tempered fractional differential equation with absorbing boundary conditions [14], i.e.,

[TABLE]

where $\alpha,\beta\in(0,2)\backslash\{1\}$, $\lambda>0$, and $\kappa_{1},\kappa_{2}>0$ in the domain $\Omega=[a,b]\times[c,d]$ and $J=[0,T]$; the function $f\in L^{2}(J;L^{2}(\Omega))$ is a source term; the convection coefficient b is a given divergence-free velocity field, i.e., $\nabla\cdot\mathbf{b}=0$, being supposed to satisfy $\textbf{b}\in L^{\infty}(J;W^{1,\infty}(\Omega)^{2})$, and the initial function $u_{0}\in L^{2}(\Omega)$. As for the discussion of the adaptivity of the fractional problems, we start from the steady state version of (1.1) with $\textbf{b}=0$. The first part of the paper focuses on designing the DG scheme of Eq. (1.1) with genuinely unstructured grids, and offering explicit theoretical analyses. Being different from [26], which constructs the LDG scheme by rewriting the fractional equation as a first order system, we adopt the primal DG methods, namely interior penalty (IP) method, still keeping the advantages over the classical continuous Galerkin method in facilitating $hp$-adaptivity and yielding block diagonal mass matrices in time-dependent problems. Generally, the non-ignorable drawback of the IP method is to specify sufficient large penalty parameter for guaranteeing numerical stability, which degrades the performance of the iterative solver of the linear system [30]. Fortunately, for the (tempered) fractional equations, this weak point disappears, since the schemes are stable for any value of the penalty parameter, say, simply taking as $1$. For the convection term, the upwind flux [10, 27] is used for ensuring numerical stability.

Mesh adaption is the basic technique of balancing the computational cost and accuracy, which introduces extra points near the singularities or the high gradient part of the solution to be computed. The key ingredient of adaptivity is a posteriori error estimators [5, 8], which are computable quantities depending on the computed solution and data. We derive a posteriori error estimators for fractional operators and obtain the local error indicators, being used to dynamically and locally refine or coarsen meshes. To show the effectiveness of the local error indicators, we adaptively solve the fractional differential equations with singularities, including both the steady state and time dependent ones. For the steady state equations, two schemes are presented. One is based on energy norm, while another one is based on dual weighted residual. It is observed that the latter performs better than the former. For the time dependent one, both the space mesh and time-step size are adapted simultaneously.

The outline of this paper is as follows. Section 2 is composed of five subsections. The first subsection reviews the definitions and properties of tempered fractional calculus. The notations and the variational formulations of DG schemes are, respectively, proposed in the second and third subsections. In the fourth subsection, we perform the stability analysis and error estimates for the two dimensional tempered fractional convection-diffusion equations. The numerical results are provided in the last subsection. Section 3 is discussing adaptivity, composed of two parts, which are, respectively, for the stationary equation and the evolution equation. We conclude the paper with some remarks in the last section.

2. DG for tempered fractional convection-diffusion equation

In this section, we design the DG scheme for the tempered fractional convection-diffusion equation (1.1), provide the detailed stability and convergence proof, and numerically verify the theoretical results.

2.1. Tempered fractional operators

We firstly introduce some preliminary definitions of tempered fractional calculus [12, 18].

Definition 2.1**.**

For any $\alpha>0,\lambda>0$, the left and right tempered Riemann-Liouville fractional integrals of function $u(x)$ defined on $\mathbb{R}$ are given by

[TABLE]

and

[TABLE]

Definition 2.2**.**

For any $\alpha>0,n-1<\alpha<n,n\in\mathbb{N}^{+},\lambda>0,$ the left and right tempered Riemann-Liouville fractional derivatives of function $u(x)$ defined on $\mathbb{R}$ are given as

[TABLE]

and

[TABLE]

Definition 2.3**.**

For any $\alpha>0,n-1<\alpha<n,n\in\mathbb{N}^{+},\lambda>0,$ the left and right tempered Caputo fractional derivatives of function $u(x)$ defined on $\mathbb{R}$ are described by

[TABLE]

and

[TABLE]

Definition 2.4**.**

The Riesz tempered fractional derivatives $\nabla_{x}^{\alpha,\lambda}$ and $\nabla_{y}^{\beta,\lambda}$ with $\alpha,\beta\in(0,2)\backslash\{1\}$, $\lambda>0$, are respectively defined as

[TABLE]

and

[TABLE]

where $\kappa_{\alpha}=\frac{1}{2\cos(\alpha\pi/2)}$ and $\kappa_{\beta}=\frac{1}{2\cos(\beta\pi/2)}$. The left and right tempered Riemann-Liouville fractional derivatives are defined by

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

respectively. The definition of $\nabla_{y}^{\beta,\lambda}$ is similar.

Remark 2.5*.*

For $\lambda=0$, the tempered Riemann-Liouville fractional derivative reduces to the Riemann-Liouville fractional derivative and the tempered Riesz fractional derivative to the Riesz fractional derivative. The fractional substantial derivative [7, 17] has similar definition to ${}_{-\infty}D_{x}^{\alpha,\lambda}$, but their physical backgrounds are totally different.

Because of the absorbing boundary condition of Eq. (1.1), and Definitions 2.1 and 2.2, we have ${}_{-\infty}I_{x}^{\alpha,\lambda}u(x)={}_{a}I_{x}^{\alpha,\lambda}u(x)$ and ${}_{x}I_{\infty}^{\alpha,\lambda}u(x)={}_{x}I_{b}^{\alpha,\lambda}u(x)$; it is similar for the tempered fractional derivative. For convenience, we use the latter notations in the following. First we introduce some properties of the tempered fractional calculus and the tempered fractional spaces. Suppose that the function $f(x)$ is $(m-1)$-times continuously differentiable in the interval $[a,b]$ and that its $m$-times derivative is integrable in $[a,b]$. Then for any $m-1<\mu,\nu<m,\lambda>0$,

[TABLE]

[TABLE]

[TABLE]

where $D^{j,\lambda}$ denotes $(\lambda+D)^{j}$ when $j$ is an integer. It can be noted that the difference between tempered Riemann-Liouville fractional derivatives and tempered Caputo fractional derivatives is the sum of the values of the derivatives of $f(x)$ at boundary. So, we state the following condition.

Condition A: For $m-1<\mu<m$, the function $f(x)$ satisfies

[TABLE]

From the discussion of ([25], p. 75-77), we know that if $f(x)$ is sufficiently smooth, then the conditions

[TABLE]

are equivalent to Condition A. Therefore, under Condition A, the formulae (2.13), (2.14) and (2.15) become

[TABLE]

[TABLE]

[TABLE]

Lemma 2.6** (Adjoint property).**

For any $\alpha>0,\lambda>0,$ the left and right tempered Riemann-Liouville fractional integral operators are adjoint for any functions $u(x),v(x)\in L^{2}([a,b]),$ i.e.,

[TABLE]

Lemma 2.7** (Adjoint property).**

For any $\alpha>0,\lambda>0,$ the left and right tempered Riemann-Liouville fractional derivative operators are adjoint for functions $u(x)$ and $v(x)$ under Condition A, i.e.,

[TABLE]

Lemma 2.8**.**

For $u(x)\in L^{2}(\mathbb{R})$ and $\mu,\lambda>0$, it holds that

[TABLE]

[TABLE]

If $u(x)\in C_{0}^{\infty}(\mathbb{R})$ further, then

[TABLE]

[TABLE]

where $\hat{u}(\omega)=\int_{-\infty}^{\infty}~{}\mathrm{e}^{-i\omega x}u(x)dx$.

Next, we discuss the tempered fractional Sobolev space. Denote $\Omega=[a,b]\times[c,d]$ as a finite domain; $A\lesssim B$ means that $A$ can be bounded by a multiple of $B$, which is independent of $B$; and $A\sim B$ means that $A\lesssim B\lesssim A$. The definitions and properties for the norms of the left and right tempered fractional derivatives are similar, so we mainly focus on the left one.

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

Let $p,q>0$ and $\alpha>0$. Then

[TABLE]

More specifically,

[TABLE]

and

[TABLE]

Let $\alpha,\beta>0$ and $u(x,y)\in L^{2}(\mathbb{R}^{2})$. Define the semi-norms

[TABLE]

[TABLE]

and norms

[TABLE]

[TABLE]

Next, we define the norms for the functions in $H^{\mu,\lambda}(\mathbb{R}^{2})$ in terms of the Fourier transform, i.e.,

[TABLE]

[TABLE]

[TABLE]

where $\omega_{1}$ and $\omega_{2}$ are two components of $\omega$, and $\hat{u}$ is the Fourier transform of $u(x,y)$. By Plancherel’s theorem, it can be proved that the spaces $J_{L,x}^{\mu,\lambda}$ and $H_{x}^{\mu,\lambda}$ are equal with equivalent semi-norms and norms while the spaces $J_{L,y}^{\mu,\lambda}$ and $H_{y}^{\mu,\lambda}$ are also equal with equivalent semi-norms and norms.

Noticing that $|\omega|^{2}=\omega_{1}^{2}+\omega_{2}^{2}$, together with (2.20), we have

[TABLE]

Therefore,

[TABLE]

which gives an equivalent form of the semi-norm $H^{\mu,\lambda}(\mathbb{R}^{2})$, being useful in the error estimate. For any $\mu>0$, the Sobolev space $H^{\mu}(\mathbb{R}^{2})$ is defined with the semi-norm

[TABLE]

and the norm

[TABLE]

In the following, we use $H^{\mu}(\Omega)$ to denote the space of functions on $\Omega$ that admit extensions to $H^{\mu}(\mathbb{R}^{2})$, equipped with the quotient norm $\|u\|_{H^{\mu}(\Omega)}:=\inf\limits_{\tilde{u}}\|\tilde{u}\|_{H^{\mu}(\mathbb{R}^{2})}$, where the infimum extends over all possible $\tilde{u}\in H^{\mu}(\mathbb{R}^{2})$ such that $\tilde{u}=u$ on $\Omega$ (in the sense of distributions).

Taking notice of (2.20) again, it holds that

[TABLE]

Therefore,

[TABLE]

Then we have the following lemmas from [12] with extension of the space $H_{0}^{\mu}(\Omega)$ to two dimension.

Lemma 2.10**.**

Let $0<\mu_{1}<\mu_{2}$, and $\mu_{1},\mu_{2}\neq n-\frac{1}{2}\,(n\in\mathbb{N}^{+})$. If $u\in H_{0}^{\mu_{2}}(\Omega)$, then

[TABLE]

Lemma 2.11**.**

If $u\in H_{0}^{\mu}(\Omega)$ with $\mu\in(0,1),\mu\neq\frac{1}{2},$ for all $\lambda\geq 0$, then

[TABLE]

and if $0<\mu_{1}<\mu$ and $\mu_{1}\neq\frac{1}{2}$, it holds that

[TABLE]

Lemma 2.12**.**

Let $0<\alpha<1$, and $\lambda>0$. If $f(x,y)\in H_{0}^{\alpha}(\Omega)$, then

[TABLE]

and

[TABLE]

Remark 2.13*.*

Lemma 2.12 shows the coercivity of the tempered fractional operator since $\cos(\pi\alpha/2)>0$ when $0<\alpha<1$. If $1<\alpha<2$, the tempered fractional operator is also coercive; for the details of the proof, see [12].

2.2. Notations for DG methods

Denote $\Omega_{h}$ as a conforming subdivision of $\Omega\subset\mathbb{R}^{2}$ with non-overlapping triangles. For an element $T\in\Omega_{h}$, $h_{T}$ denotes its diameter, and $h=\max_{T\in\Omega_{h}}h_{T}$. The family of meshes $\Omega_{h}$ is assumed to be shape-regular, i.e., there exists a constant $c$ such that for all $T\in\Omega_{h}$, $\frac{h_{T}}{\rho_{T}}\leq c$, where $\rho_{T}$ denotes the radius of the largest inscribed ball in $T$. Denote $\Gamma$ as the union of the boundaries of the elements $T$ of $\Omega_{h}$, and $\Gamma_{i}$ the set of interior faces of $\Omega_{h}$, i.e., the set of faces that are not included in the boundary $\partial\Omega$. Let $\Gamma_{b}$ be the faces that are included in $\partial\Omega$. Then $\Gamma=\Gamma_{i}\cup\Gamma_{b}$. Associated with the mesh $\Omega_{h}$, we define the broken Sobolev spaces

[TABLE]

and

[TABLE]

equipped with the broken Sobolev norm

[TABLE]

where

[TABLE]

The $H^{s,\lambda}(T)$ norm is defined through (2.21). And it also can be defined by the right derivatives.

We define the DG finite element space as follows:

[TABLE]

where $P_{N}(T)$ denotes the set of polynomials of degree less than or equal to $N$. The global solution can be approximated as

[TABLE]

and the local solution $u(\mathbf{x},t)$ can be expressed by

[TABLE]

where $l_{j}(\mathbf{x})$ denotes the two-dimensional multivariate Lagrange interpolation basis function, and $N_{p}=(N+1)(N+2)/2$ is degree of freedom of one element.

Next, we introduce some notations to manipulate numerical fluxes. If two elements $T_{e}^{1}$ and $T_{e}^{2}$ are neighbors and share one common side $e$, there are two traces of the function $v$ along $e$. We assume that the normal vector $\mathbf{n}_{e}$ is oriented from $T_{e}^{1}$ to $T_{e}^{2}$, and denote

[TABLE]

The definition of jump and average to sides that belong to the boundary $\partial\Omega$ is:

[TABLE]

We also introduce two bilinear forms $J_{0},\,J_{1}:H^{s,\lambda}(\Omega_{h})\times H^{s,\lambda}(\Omega_{h})\rightarrow\mathbb{R}$ that penalize the jump of the function values and the jump of the normal derivative values:

[TABLE]

Remark 2.14*.*

Reference [5] presents numerical analysis for the advection-reaction equation in the case of continuous, piecewise affine approximations (using the continuous interior penalty method (CIP)) and discontinuous, piecewise affine or constant approximations (DG). In both of the two cases, stability is obtained by using interior penalty; for CIP, the gradient jumps over element faces are penalized, and for DG, the jumps of the solution.

2.3. Variational formulation

From the definition of tempered fractional derivatives, we rewrite Eq (1.1) as a more clear form

[TABLE]

Here, we mainly discuss the case $\alpha,\beta\in(0,1)$. Then we formulate the weak variational formulation. For $f\in H^{-s,\lambda}(\Omega)$, find $u\in V:=H^{s,\lambda}_{0}(\Omega)\bigcap H^{1}(\Omega)$, where $s=\max\{\alpha/2,\beta/2\}$, such that $\forall v\in V$,

[TABLE]

Denote two bilinear forms as $a_{x}^{\alpha}(u,v)=({}_{a}D_{x}^{\frac{\alpha}{2},\lambda}u,{}_{x}D_{b}^{\frac{\alpha}{2},\lambda}v)+({}_{a}D_{x}^{\frac{\alpha}{2},\lambda}v,{}_{x}D_{b}^{\frac{\alpha}{2},\lambda}u),$ and $a_{y}^{\beta}(u,v)=({}_{c}D_{y}^{\frac{\beta}{2},\lambda}u,{}_{y}D_{d}^{\frac{\beta}{2},\lambda}v)+({}_{c}D_{y}^{\frac{\beta}{2},\lambda}v,{}_{y}D_{d}^{\frac{\beta}{2},\lambda}u),$ and define

[TABLE]

and

[TABLE]

Then the variational formulation becomes

[TABLE]

where $\kappa=2\lambda^{\alpha}\kappa_{1}\kappa_{\alpha}+2\lambda^{\beta}\kappa_{2}\kappa_{\beta}>0$. So the DG bilinear form should be defined as

[TABLE]

and

[TABLE]

where $\hat{u}$ corresponding to the element boundary terms from integration by parts is the so-called numerical flux. It is a single valued function defined on the faces and should be designed based on different guiding principles for different PDEs to guarantee stability and optimal order of convergence. Here we choose an upwind flux. Denote the upwind value of a function $v$ by $v^{\mathrm{up}}$. We recall that $\mathbf{n_{e}}$ is a unit normal vector pointing from $T_{e}^{1}$ to $T_{e}^{2}$:

[TABLE]

Then the DG variational formulation can be stated as: find $u\in H^{s,\lambda}(\Omega_{h})\bigcap H^{1}(\Omega_{h})$, such that $\forall v\in H^{s,\lambda}(\Omega_{h})\bigcap H^{1}(\Omega_{h})$ there exists

[TABLE]

We discrete the time derivative with backward Euler. Let $N_{T}$ be a positive integer and $\Delta t=T/N_{T}$ denote the time step. We also use the notations:

[TABLE]

Define the orthogonal projection operators, $P_{h}:L_{2}(\Omega_{h})\,\rightarrow\,V_{h}$, i.e., for each element $T$,

[TABLE]

The fully discrete DG scheme is as follows: find $u_{h}^{n+1}\in V_{h}$, such that $\forall v\in V_{h}$, there exists

[TABLE]

with known $u_{h}^{n}$, and if $n=0,~{}u_{h}^{0}=P_{h}u_{0}(x)$.

2.4. Stability analysis and error estimates

Lemma 2.15** (Discrete Grönwall inequality [27]).**

Let $\Delta t,B,C>0$ and $(a_{n})_{n}$, $(b_{n})_{n}$, $(c_{n})_{n}$ be sequences of nonnegative numbers satisfying

[TABLE]

Then, if $C\Delta t<1$,

[TABLE]

Lemma 2.16**.**

For any function $f(x),g(x,t)\in L^{2}(\Omega)$, if $|f(x)|\leq\int_{t_{1}}^{t_{2}}g(x,t)dt~{}\forall x\in\Omega$, then

[TABLE]

Proof.

This inequality can be easily obtained by using Höld’s inequality. Here we omit it. ∎

In the following, $C$ denotes a generic constant independent of $h$ and $\Delta t$, which takes different values in different occurrences. First, we prove the coercivity and continuity of the bilinear form $a_{h}(u,v)$. Define the energy norm on $V_{h}$ as

[TABLE]

where

[TABLE]

Then from (2.35) and Lemma 2.12, we have

[TABLE]

where $\gamma=\min\{\kappa_{1},\kappa_{2},1\}$.

On the other hand, from the definition of $a_{h}(u,v)$, and Cauchy-Schwarz’s inequality, we have

[TABLE]

Next we deal with another bilinear form $b_{h}(u,v)$. By an upwind flux,

[TABLE]

Then

[TABLE]

Now we examine the stability property of the scheme (2.40). Taking $u_{h}^{0}=P_{h}u_{0}$ leads to the following results.

Theorem 2.17**.**

For absorbing boundary conditions, the fully discrete scheme (2.40) is unconditionally stable, and there exists a positive constant $C$ independent of $h$, such that for all $m>0$,

[TABLE]

Proof.

Taking $v=u_{h}^{n+1}$ in (2.40), using (2.42), (2.45), and Cauchy-Schwarz’s inequality,

[TABLE]

Next, we observe that:

[TABLE]

Therefore, by using Young’s inequality,

[TABLE]

Multiplying (2.47) by $2\Delta t$ and summing from $n=0$ to $n=m-1$ lead to

[TABLE]

Then the desired result is obtained by using Lemma 2.15. ∎

Assuming that the solution of Eq. (1.1) is sufficiently regular, we have the following error estimates, showing that the numerical $L^{2}$ error is $\mathcal{O}(h^{k+\frac{1}{2}}+\Delta t)$.

Theorem 2.18**.**

Let $u^{n}$ be the exact solution of (2.38), $u_{h}^{n}$ the numerical solution of the fully discrete scheme (2.40). Then

[TABLE]

Proof.

As usual, we denote the error $e^{n}=u^{n}-u_{h}^{n}$ by two parts $\rho=u^{n}-P_{h}u^{n}$ and $\theta=u_{h}^{n}-P_{h}u^{n}$. From (2.38) and (2.40), we have

[TABLE]

Noting that $e^{n}=\rho^{n}-\theta^{n}$, we get

[TABLE]

Taking $v=\theta^{n+1}$, similar to the proof of stability, we obtain

[TABLE]

where $T_{1}=(\frac{\rho^{n+1}-\rho^{n}}{\Delta t},\theta^{n+1})$, $T_{2}=b_{h}(\rho^{n+1},\theta^{n+1})$, $T_{3}=a_{h}(\rho^{n+1},\theta^{n+1})$, $T_{4}=\kappa(\rho^{n+1},\theta^{n+1})$, and $T_{5}=(\frac{\partial}{\partial t}u^{n+1}-\frac{u^{n+1}-u^{n}}{\Delta t},\theta^{n+1})$.

Since

[TABLE]

by Lemma 2.16, we have

[TABLE]

Hence, with Höld’s, Young’s inequalities, we obtain

[TABLE]

From the definition of the projection $P_{h}$ (2.39) and trace inequalities, we have

[TABLE]

From the continuity of $a(u,v)$ (2.43), we obtain

[TABLE]

In the last inequality, using embedding theorem and trace theorem, three terms in the energy norm $\|\rho^{n+1}\|_{E(\Omega_{h})}$ both can be bounded by $|\rho^{n+1}|_{H^{1/2}(\Omega_{h})}$. Then, for a continuous interpolation function $\Pi u$ of $u$,

[TABLE]

which implies $\|\rho^{n+1}\|_{E(\Omega_{h})}^{2}\leq Ch^{2k+1}$.

Similarly, for the fourth term, we have

[TABLE]

The Taylor expansion with integral remainder has the form

[TABLE]

Thus,

[TABLE]

Then we get

[TABLE]

Substituting $T_{i},\>i=1,\cdots,5$ into (2.49), we have

[TABLE]

where $\epsilon_{2}$ and $\epsilon_{3}$ are chosen as sufficiently small numbers such that $\epsilon_{2}+\epsilon_{3}\leq\gamma$.

From the definition of energy norm (2.41), we know

[TABLE]

Then

[TABLE]

Multiplying (2.50) by $2\Delta t$, summing over $n$ from [math] to $N_{T}-1$, and using the discrete Grönwall inequality with $\theta^{0}=0$, we get

[TABLE]

By the triangle inequality, we obtain the desired result. ∎

2.5. Numerical experiment

In this section, we offer the numerical performance of the proposed schemes for two examples to validate the preceding theoretical analysis. We use the backward Euler discretization to solve the method-of-line fractional PDE, i.e., the classical ODE system. We take the time steps $\Delta t$ to be $h^{N+1}$, where $N$ denotes the order of polynomial of finite element space. As to the spatial approximation, we adopt the interpolation bases [16].

We first introduce the local and global vector and matrix notations,

[TABLE]

[TABLE]

and $\mathbf{u}_{h}^{n}$ denotes the value of $\mathbf{u}_{h}$ at time $t_{n}$. Let $f_{j,k}=(f,l_{j}^{k}(\mathbf{x}))_{T^{k}}$. Similarly, denote

[TABLE]

[TABLE]

and let $\mathbf{F}_{h}^{n}$ be the value of $\mathbf{f}_{h}$ at time $t_{n}$. Then, we define the local mass matrix $M^{k}$ and the local spatial stiffness matrix $S_{x}^{k},S_{y}^{k}$ at element $T^{k}$ as

[TABLE]

[TABLE]

[TABLE]

It is a little bit complex to build the tempered fractional spatial stiffness matrix, since tempered fractional operators are nonlocal and we need all the information of the related elements in $x$ direction or $y$ direction when generating any stiffness matrix of an element. We take the method of [26] and get the global tempered fractional spatial stiffness matrices ${}_{l}G_{x},{}_{l}G_{y},{}_{r}G_{x},{}_{r}G_{y}$, where $`l/r$’ denote left/right tempered fractional derivative and $`x/y$’ denote the $x$ or $y$ direction.

With the above notations, we rewrite the global fully discrete form (2.40) as

[TABLE]

where $M$, $S_{x}$, $S_{y}$ are global mass and stiffness matrices, and their non-zero diagonal blocks are constructed by $M^{k}$, $S_{x}^{k}$, and $S_{y}^{k}$ respectively.

Example 2.19**.**

Consider the problem

[TABLE]

where $\mathbf{b}=(0.5,0.5),\,\kappa_{1}=0.1,\,\kappa_{2}=0.2,\,\lambda=2,\,T=1$ and $\alpha,\beta\in(0,1)$ on the computational domain $\Omega=(0,2)\times(0,2)$. Its exact solution is $u=\mathrm{e}^{-t}x^{2}(2-x)^{2}y^{2}(2-y)^{2}$ with appropriate initial and boundary conditions.

Example 2.20**.**

Consider the same problem Eq. (2.51), but the parameters are taken as $\mathbf{b}=(0.5,0.5),\,\kappa_{1}=0.1,\,\kappa_{2}=0.2,\,\lambda=0.2,\,T=1$ and $\alpha\,,\beta\in(0,2)$ on the computational domain $\Omega=(0,2)\times(0,2)$. The exact solution is $u=\mathrm{e}^{-t}\sin{\frac{\pi}{2}x}\sin{\frac{\pi}{2}y}$.

For the numerical experiments, in order to validate the stability and the convergence of the preceding scheme, the order of convergence is calculated by

[TABLE]

Table 1 and Table 2 list the $L^{2}$ errors and convergence orders for different parameters $(\alpha,\beta)$ in different DG finite element space $P_{N}$, where $N$ denotes the degree of polynomial in two variables, and $K$ the total number of triangle elements. It can be seen that the convergence order $N+1/2$ is consistent with the theoretical prediction, even for the case of $(\alpha,\beta)\in(0,2)$.

3. Adaptive DG algorithm

This section focuses on the adaptive DG scheme for the fractional diffusion equations. We derive posteriori error estimates, and design the local error indicators. The numerical experiments are performed to show the performances of the adaptive schemes.

3.1. Stationary Equation

First we consider the simple stationary equation:

[TABLE]

A posteriori error estimators are an essential ingredient of adaptivity, which are computable quantities depending on the computed solution and data that provide information about the quality of approximation and may thus be used to make judicious mesh modifications. The ultimate purpose is to construct the estimator of meshes that will eventually be equivalent to the exact error. The usual method of constructing the estimator deals with error estimation in global norms like the ‘energy norm’ or the ‘$L^{2}$ norm’, which is the first scheme in the following.

3.1.1. Scheme 1 – energy norm [4]

For $T\in\Omega_{h}$, we define the operators $D_{x}^{\alpha}={}_{a}D_{x}^{\alpha}+{}_{x}D_{b}^{\alpha}$, $D_{y}^{\alpha}={}_{c}D_{y}^{\alpha}+{}_{y}D_{d}^{\alpha}$, and the local error estimator $\eta_{T}$ by

[TABLE]

where $R:=f-(D_{x}^{\alpha}u_{h}+D_{y}^{\alpha}u_{h})$. We will prove that $\eta_{T}$ bounds the exact error by inequalities in both directions, where the constants in these inequalities depend only on properties of the triangulation. The upper estimate shows that $\eta_{T}$ can be used as a reliable stopping criterion for the algorithm, while the lower estimate suggests that refinement based on $\eta_{T}$ will be efficient.

Lemma 3.1** (upper bound).**

There exists a constant $C_{1}$, depending only on the minimum angle of $\Omega_{h}$ and the ratio of the biggest diam to the smallest diam of the elements, such that

[TABLE]

where the energy norm $\|v\|_{E(\Omega_{h})}$ is defined in (2.41).

Proof.

The weak form reads as follows: find $u\in V:=H_{0}^{\alpha/2}(\Omega_{h})$, such that

[TABLE]

where the symmetry bilinear form

[TABLE]

The notations $a_{x}^{\alpha}(u,v)$ and $a_{y}^{\alpha}(u,v)$ are defined in (2.35).

Let $u_{h}\in V_{h}$ (defined in (2.30)) be the numerical solution, and $e_{h}:=u-u_{h}$. Then $e_{h}$ satisfies the residual equation

[TABLE]

Now we resort to the interpolation operator $\Pi:V\rightarrow V_{h}$, which satisfies the approximation property, for any $v\in V$,

[TABLE]

Define

[TABLE]

where $c_{e}=1/2$ if $e\in\Gamma_{i}$ and $c_{e}=1$ if $e\in\Gamma_{b}$. Then

[TABLE]

and thus

[TABLE]

Taking $v=e_{h}\in H^{\alpha/2}(\Omega_{h})$, we have

[TABLE]

∎

Lemma 3.2** (lower bound).**

There exists constant $C_{2}$ depending only on the minimum angle of $\Omega_{h}$, such that

[TABLE]

where $Q_{h}$ is the $L^{2}$ orthogonal projection onto the space of (discontinuous) piecewise polynomials of degree $k-1$.

Proof.

Let $\zeta$ be the bubble function in element $T$, vanishing outside $T$. Then for any polynomial $\phi$ [4],

[TABLE]

Taking $v=\zeta Q_{h}R\in V_{h}$ in (3.4), we have

[TABLE]

and

[TABLE]

Thus,

[TABLE]

where the inverse inequality is used in the last inequality; the proof of the inverse inequality is given in Appendix A. Then using (3.9)

[TABLE]

and hence

[TABLE]

Next, we consider the jump term. Similarly, let $\zeta$ be the bubble function of one edge $e\in T$, vanishing outside $T$. Taking $v=\zeta[u_{h}]\in V_{h}$ in (3.4),

[TABLE]

Then using inverse inequality and trace inequality lead to

[TABLE]

and

[TABLE]

Therefore, for small $h$,

[TABLE]

and thus

[TABLE]

The proof is completed by combining all these estimates. ∎

The goal of adaptive methods is the generation of a mesh which is adapted to the problem such that a given criterion, like a tolerance for the estimated error between exact and discrete solution, is fulfilled by the discrete solution on this mesh. An optimal mesh should be as coarse as possible while meeting the criterion, in order to save computational time and memory requirements. A global refinement of the mesh would lead to the best error reduction, but the amount of new unknowns might be much larger than needed to reduce the error below the given tolerance. We use the Marking Strategy C in [24] to deal with the error estimator and the oscillation simultaneously. More strategies can be found in [29], like Maximum strategy and Equidistribution strategy.

First, we define

[TABLE]

Marking Strategy C: Given a parameter $0<\theta_{1},\theta_{2}<1$, construct a subset $\hat{\Omega}_{h}$ of $\Omega_{h}$ such that

[TABLE]

Enlarge $\hat{\Omega}_{h}$ such that

[TABLE]

Then we will refine the mesh $\Omega_{h}$ by Marking Strategy C, and show that the error $e_{h}$ converges to zero in the energy norm at each refinement. Before this, we include the subscripts $j$ to signify the symbol corresponding to the $j$th refinement. We have

[TABLE]

and

[TABLE]

The following two lemmas can be proved similarly as in [4]. The main difference is the definition of the energy norm. Here we omit the details.

Lemma 3.3** (error reduction).**

Suppose an element $T\in\Omega_{j}$ contains a node of $\Omega_{j+1}$ in its interior. Then we have

[TABLE]

where the positive constant $C$ depends only on the minimum angle of $\Omega_{j}$.

Lemma 3.4** (oscillation reduction).**

Let $T\in\Omega_{j}$ be subdivided into elements in $\Omega_{j+1}$ such that

[TABLE]

where the positive constant $\gamma\leq 1$. There exist constants $0<\rho_{1}<1$ and $0<\rho_{2}$, depending only on $\gamma$, such that

[TABLE]

Combining Lemma 3.3 and Lemma 3.4 with Galerkin orthogonality

[TABLE]

we can easily get the following convergence proposition, which shows that $u_{j}$ converges to $u$ in the energy norm.

Corollary 3.5**.**

Let $\{u_{j}\}_{j\geq 1}$ be the sequence of finite element solutions generated by Marking Strategy C. There exist positive constants $\chi$ and $\xi$ such that $0<\xi<1$ and

[TABLE]

3.1.2. Scheme 2 – Dual Weighted Residual method (DWR)

Different from the scheme above, the traditional approach to adaptivity aiming at estimating the error with respect to the generic energy norm of the problem, or the global $L^{2}$ norm, the Dual Weighted Residual method (DWR) [3] for goal-oriented error estimation aims at economical computation of arbitrary quantities of physical interest. This is typically required in the design cycles of technical applications. ‘Goal-oriented’ adaptivity is designed to achieve these tasks with minimal cost.

When solving the fractional problems, the DWR method is significantly better than the traditional approach since the fractional operator and energy norm is nonlocal. In detail, for left Riemann-Liouville fractional operator, the numerical solution on one element $T$ is affected by all the elements on the left. Therefore, the lower bound is not local absolutely and the over refinement may occur. While DWR multiplies every local error indicator on one element by a weight, consisting of the dual solution. It has the feature of a ‘generalized’ Green function $G(T,T^{\prime})$, which describes the dependence of the target error quantity $J(e_{h})$ concentrated at some element $T$; we select it as $\|e_{h}\|_{E(T)}$ in the following, on local properties of the data, i.e., the error estimator on element $T^{\prime}$.

Following the general concept of the DWR method, let $z\in V$ be the solution of the associated dual problem

[TABLE]

and $z_{h}\in V_{h}$ be discontinuous finite element approximation defined by

[TABLE]

where the bilinear form $a(\cdot,\cdot)$ is defined as

[TABLE]

Using this construction together with Galerkin orghogonality, we obtain

[TABLE]

Thus, we define the local error indicator

[TABLE]

Taking $J(\varphi)$ such that $J(e_{h})=\|e_{h}\|_{E(\Omega_{h})}$, then we have the global upper bound in energy norm,

[TABLE]

If we take a rough estimate to $J(e_{h})$ by $\|z-z_{h}\|_{L^{2}(\Omega)}\leq Ch^{\alpha/2}|z|_{{H^{\alpha/2}}(\Omega_{h})}$ and $\|[z-z_{h}]\|_{L^{2}(\partial T)}\leq C\|z\|_{L^{2}(\partial T)}$, then we have, similar to (3.6), from (3.22), with a priori analysis in forms of bounds for $z$,

[TABLE]

Here, we get a global posteriori error estimate based on energy norm that is consistent with Scheme 1. In this sense, a posteriori error estimate based on DWR is more meticulous. In order to evaluate the posteriori error representation (3.24), we need information about the discontinuous dual solution $z$. Since in practice, $z$ is not explicitly known, such information has to be obtained either through a priori analysis in form of bounds for $z$ in certain Sobolev norms or through computation by solving the dual problem numerically.

Here, we approximate $z$ by a high-order method. We take $V_{h}$ as the linear discontinuous finite element and solve the dual problem by using quadratic discontinuous finite element on the current mesh yielding an approximation $z_{h}^{(2)}\in V_{h}^{(2)}$ to $z$, and $z_{h}$ can be got by linear interpolation of $z$. This yields the approximate local error indicator

[TABLE]

Remark 3.6*.*

In addition to the above points that a posteriori error estimate based on DWR is better than a global posteriori error estimate based on energy norm, the former has advantages when the derivative on $x-$direction and on $y-$direction is different, i.e., $\alpha\neq\beta$. When discussing a posteriori error estimate based on energy norm, we take $\alpha=\beta$ in Eq (3.1) for convenience. Actually, if $\alpha<\beta$, the indicator is not easy to select. We must take the exponent of $h_{T}$ in the indicator $\eta_{T}:=h_{T}^{\alpha/2}\|R\|_{L^{2}(T)}$ to be $\alpha/2$ to guarantee the upper bound and take it to be $\beta/2$ to guarantee the lower bound. A posteriori error estimate based on DWR avoids this problem, and it is still effective for complex problems.

3.1.3. Numerical experiment

In this section we will present some numerical experiments using the two schemes given above. We compare them with uniformly refinement and with each other.

Example 3.7**.**

Consider the 1D fractional equation on the domain $\Omega=[0,2]$ with $\alpha=0.8$,

[TABLE]

and its dual problem is

[TABLE]

The source term $f$ is chosen such that the exact solution is $u=\left(1-(x-1)^{2}\right)^{\gamma}$, $\gamma=0.7$, which has poor regularity near the boundary, and we use the discontinuous piecewise linear function for approximation.

This is a boundary layer problem, i.e., the solution has less regularity around the endpoints of domain $[0,2]$, where the mesh should be finer. We initially divide the interval $[0,2]$ into 8 cells uniformly, then refine the mesh based on the energy-norm indicator in Figure 1 (left) and the weighted-norm indicator in Figure 1 (right). The numerical solution approximates the exact solution both very well obviously. And from an intuitive point of view, the right one is better than the left one. Besides, we compare the two kinds of refinements with uniform refinement, and show their convergence rates. A reference line is provided which shows the optimal convergence rate $N^{-2}$. As shown in Figure 2, the uniform refinement is the worst while the refinement based on the weighted indicator is the best, almost achieving the optimal convergence rate.

Example 3.8**.**

Consider the 2D fractional equation on the domain by $\Omega:=[0,2]\times[0,2]$ with $\alpha=0.2$, $\beta=0.8$,

[TABLE]

and its dual problem is itself. The source term $f$ is chosen such that the exact solution writes

[TABLE]

Being the same as the case of 1D, we compare the three kinds of refinements: uniform refinement, the refinement based on energy-norm indicator, and weighted indicator. This exact solution has less regularity inside the domain, i.e., interior layer. Figure 3 is the shape of the exact solution. The adaptive mesh will be finer in the steep region. The experiment datum are provided in Table 3, Table 4, and Table 5. We denote the total error indicator by $\eta$, and the effectiveness index $I_{eff}=\eta/\|e\|_{E(\Omega_{h})}$. If the effectiveness index $I_{eff}$ remains roughly constant in different meshes, the indicator $\eta$ approximates the true error $e_{h}$ well. By comparing these tables, adaptive refinement is significantly better than uniform refinement. The weighted indicator is slightly better than the energy-norm indicator since the $L^{2}$ norm error and energy-norm error of the former is smaller than the latter, and the effectiveness index $I_{eff}$ of the former changes more moderately than the latter. In Figure 4, two kinds of adaptive meshes are presented, and both refine the steep region. The right one (based on the weighted indicator) is slightly better than the left one (based on the energy norm indicator).

3.2. Evolution Equation

In this section, we consider the time dependent tempered fractional equation:

[TABLE]

Different from the stationary equation, the mesh here is adapted to the solution in every time step using a posteriori error indicator, i.e., the adaptive algorithm solving the evolution equation at the $n-$th time step reads as

[TABLE]

Here the refinement/coarsening procedure includes both the mesh and time-step size modifications. In this paper, we propose the following algorithm, similar to [8], to modify the time-step size $\tau_{n}$ and mesh $\Omega_{h}^{n}$ starting from the initial time-step size $\tau_{n,0}=\tau_{n-1}$ and initial mesh $\Omega_{h}^{n,0}=\Omega_{h}^{n-1}$:

1. Refine the time-step size $\tau_{n,0}$ to the final time-step size $\tau_{n}$ such that the associated time error indicators are less than the prescribed tolerances.

2. Refine/Coarsen the mesh $\Omega_{h}^{n,0}$ to the final mesh $\Omega_{h}^{n}$ such that the associated space error indicators are less than the prescribed tolerances.

3. Enlarge the initial time-step size $\tau_{n+1,0}$ for next time step if the current time error indicator is much less than the tolerances.

3.2.1. A posteriori error analyses

While the previous section is concentrated on the spatial mesh refinement based on the local error indicator, this section describes the process of evolution of the solution and the time-step size based on the time error indicator. Then we introduce the time and space local error indicator, and prove global upper bound and local lower bound. For convenience, we take $\mathbf{b}=\mathbf{0}$, $\alpha=\beta$, and let $V_{h}^{n}$ indicate the usual space of linear discontinuous finite element in $\Omega^{n}$ and piecewise constant discontinuous in $\tau_{n}$ [3]. Let $u\in V$ be the solution of

[TABLE]

and $u_{h}^{n}$ be the fully discrete discontinuous finite element approximation defined by

[TABLE]

where $\partial_{t}u_{h}^{n}=\frac{u_{h}^{n}-u_{h}^{n-1}}{\tau_{n}}$ and $\bar{f}^{n}=\frac{1}{\tau_{n}}\int_{t^{n-1}}^{t^{n}}f(\mathbf{x},t)dt$.

Theorem 3.9** (upper bound).**

For any integer $1\leq m\leq N$, there exists a positive constant $C$ depending only on the minimum angle of meshes $\Omega_{h}^{n},\ n=1,2,\cdots,m,$ such that the posteriori error estimate

[TABLE]

holds, where the time error indicator and space error indicator are given by

[TABLE]

[TABLE]

Proof.

From (3.27) and (3.28), for a.e. $t\in(t^{n-1},t^{n}]$, and for any $v\in V$, we have

[TABLE]

where $R^{n}=\bar{f}^{n}-\frac{u_{h}^{n}-u_{h}^{n-1}}{\tau_{n}}-\kappa_{1}\kappa_{\alpha}D_{x}^{\alpha,\lambda}u_{h}^{n}-\kappa_{2}\kappa_{\alpha}D_{y}^{\alpha,\lambda}+\kappa u_{h}^{n}$ and $\Pi^{n}$ is the interpolation operator.

Taking $v=u-u_{h}^{n}:=e^{n}$, similar to (3.6), we have

[TABLE]

Integrating the above formula in time from $t^{n-1}$ to $t^{n}$ and summing over $n$ from $1$ to $m$, with the discrete Grönwall inequality, we complete the proof. ∎

Next, we prove the lower bound to ensure over-refinement will not occur based on our space error indicator. Let $u_{\ast}^{n}\in V$ be the solution of auxiliary problem [8].

[TABLE]

Note that for fixed time-step size $\tau_{n}$, by adapting the mesh $\Omega_{h}^{n}$, we are essentially controlling the error between $u_{h}^{n}$ and $u_{\ast}^{n}$, not between $u_{h}^{n}$ and the exact solution $u$. Based on this observation, we have the following analyses.

Theorem 3.10** (lower bound).**

There exist constants $C_{2},C_{3}>0$ depending only on the minimum angle of $\Omega_{h}^{n}$ such that for any $T\in\Omega_{h}^{n}$, the following estimate holds:

[TABLE]

where $\hat{C}_{n}=\underset{T\in\Omega_{h}^{n}}{\max}(h_{T}^{\alpha}/\tau_{n})$ and $\mathrm{osc}^{n}(T)^{2}=h_{T}^{\alpha}\|Q^{n}R^{n}-R^{n}\|_{L^{2}(T)}^{2}$.

Proof.

Let $p$ be a node of $\Omega_{h}^{n}$ interior to $T$ and $\xi$ be a piecewise linear function with respect to $\Omega_{h}^{n}$ that equals $1$ at $p$ and [math] at all the other nodes. Then $v=\xi Q^{n}R^{n}$ belongs to $V_{h}^{n}$ and vanishes outside $T$. Thus

[TABLE]

Combining with (3.29), we find

[TABLE]

Therefore, we have

[TABLE]

Thus,

[TABLE]

The term $\|[u_{h}]\|_{L^{2}(\partial T)}^{2}$ can be bounded by the right hand side similarly.

Combining the above estimates implies (3.30). ∎

3.2.2. Numerical experiment

In the last part, we mainly observe the performance of spatial adaptive mesh, and give the specific error of different cases on experiments. Here we present some numerical experiments using the algorithm of evolution equations given above. We focus on observing the evolution of exact solution by taking it with different regularity at different time, which will indicate the effectiveness of the indicator.

Example 3.11**.**

Consider the equation (3.26) on the domain $\Omega=[0,2]\times[0,2]$, with $\mathbf{b}=\mathbf{0},\kappa_{1}=0.1,\kappa_{2}=0.2,\lambda=0.2$, and $\alpha=\beta=0.8$. The source term $f$ is chosen such that the exact solution writes

[TABLE]

For this example, the solution has less regularity near the point $(t,t)$, where we expect to observe more refinement. The adaptation process yields the meshes shown in Figure 5, which is in consistent with the theoretical result completely.

4. Conclusion

We discuss the DG methods for the two dimensional time dependent space tempered fractional convection-diffusion equations, especially their adaptivity; and the provided strategy for designing adaptive schemes works for the general PDEs with fractional operators. DG methods are superior to finite element method in many ways, but the stability of the discrete schemes needs more attention. The interior boundaries are usually connected by two ways: interior penalty and numerical flux. We use the former to deal with the diffusion term and the latter to the convection term. The stability and convergence analyses are explicitly provided. The theoretical results are confirmed by numerical experiments. For the adaptivity of the DG methods, we consider two schemes of the stationary problem, i.e., a posteriori error estimator based on the traditional energy norm and another one based on the dual weighted residual. The numerical experiments confirm the advantage of the DWR method. Finally, we consider the fractional evolution problem. A posteriori error estimate is provided and the indicator is designed; its effectiveness is displayed by numerical simulations.

5. Acknowledgements

This work was supported by the National Natural Science Foundation of China under Grant No. 11671182.

Appendix A The inverse inequation

Lemma A.1**.**

Let $\alpha>0$, $v\in V_{k}$ and $T\in\Omega_{h}$. Then there exists a constant $C$ depending only on $\alpha$ and the shape regularity such that

[TABLE]

where $h_{T}$ is the diameter of $T$.

Proof.

First, we consider the one dimensional case, i.e., $T=[0,h]$. Let $\hat{T}=[0,1]$ be the reference triangle. If $v$ is a function defined on $T$, then $\hat{v}$ is defined on $\hat{T}$ by $\hat{v}(\hat{x})=v(h\hat{x})~{}~{}\forall\hat{x}\in\hat{T}$.

For $v$ vanishing outside of $T$, there exists

[TABLE]

Actually, taking $x=hy,~{}\xi=hs$

[TABLE]

Similarly,

[TABLE]

Thus, (A.1) is implied by

[TABLE]

On the other hand,

[TABLE]

Using the equivalence of any two norms on the finite dimensional function space $V_{k}$,

[TABLE]

Next, we consider the case that $T$ is a triangle. Let $\hat{T}$ be a reference triangle, and $F$ be an affine map from $\hat{T}$ onto $T$. If $v$ is a function defined on $T$, then $\hat{v}$ is defined on $\hat{T}$ by $\hat{v}(\hat{x})=v(x)~{}~{}\forall F(\hat{x})=x$.

Note that

[TABLE]

in which two terms can be decomposed $\alpha-$derivative in $x/y$ direction and no-derivative in $y/x$ direction.

Then inspired by the one dimensional case above, one obtains

[TABLE]

When $h\leq 1$, we have

[TABLE]

∎