Optimal order finite element approximations for variable-order time-fractional diffusion equations

TL;DR

This paper develops and analyzes a finite element method for variable-order time-fractional diffusion equations, achieving optimal convergence rates without requiring full regularity of solutions, validated by numerical experiments.

Contribution

It introduces a fully discrete finite element scheme for variable-order fractional diffusion equations with proven optimal convergence rates.

Findings

First-order accuracy in time and second-order in space achieved.

Optimal convergence estimates proved without full regularity assumptions.

Numerical experiments confirm theoretical results.

Abstract

We study a fully discrete finite element method for variable-order time-fractional diffusion equations with a time-dependent variable order. Optimal convergence estimates are proved with the first-order accuracy in time (and second order accuracy in space) under the uniform or graded temporal mesh without full regularity assumptions of the solutions. Numerical experiments are presented to substantiate the analysis.

Tables

Table 1. Table 1: Convergence rates under α ​ ( 1 ) = 0.4 𝛼 1 0.4 \alpha(1)=0.4 and α ​ ( 0 ) = 1 𝛼 0 1 \alpha(0)=1 or 0.6 0.6 0.6 .

	Uniform		Graded		Uniform			Uniform
	$\alpha \left(0\right)=0.6$		$\alpha \left(0\right)=0.6$		$\alpha \left(0\right)=1$			$\alpha \left(0\right)=1$
$N$	$h=1/32$	$\kappa$	$h=1/32$	$\kappa$	$h=1/32$	$\kappa$	$h$	$N=1/h^2$	$\gamma$
1/8	3.26E-02		3.54E-02		1.00E-02		1/8	1.60E-03
1/16	2.09E-02	0.65	1.87E-02	0.92	4.90E-03	1.03	1/16	3.97E-04	2.01
1/32	1.34E-02	0.63	9.60E-03	0.97	2.41E-03	1.03	1/24	1.76E-04	2.01
1/64	8.74E-03	0.62	4.85E-03	0.98	1.20E-03	1.01	1/32	9.89E-05	2.00

Table 2. Table 2: Convergence rates under α ​ ( 1 ) = 0.6 𝛼 1 0.6 \alpha(1)=0.6 and α ​ ( 0 ) = 1 𝛼 0 1 \alpha(0)=1 or 0.8 0.8 0.8 .

	Uniform		Graded		Uniform			Uniform
	$\alpha \left(0\right)=0.8$		$\alpha \left(0\right)=0.8$		$\alpha \left(0\right)=1$			$\alpha \left(0\right)=1$
$N$	$h=1/32$	$\kappa$	$h=1/32$	$\kappa$	$h=1/32$	$\kappa$	$h$	$N=1/h^2$	$\gamma$
1/8	1.37E-02		1.72E-02		6.65E-03		1/8	1.11E-03
1/16	7.98E-03	0.78	9.03E-03	0.93	3.34E-03	0.99	1/16	2.76E-04	2.00
1/32	4.65E-03	0.78	4.63E-03	0.96	1.68E-03	0.99	1/24	1.22E-04	2.00
1/64	2.72E-03	0.77	2.36E-03	0.97	8.60E-04	0.97	1/32	6.89E-05	2.00

Equations

$$W^{m}_{p}({\cal I};{\cal X}):=\big{\{}f:f_{t}^{(l)}(\cdot,t)\in{\cal X},~{}t\in{\cal I},~{}\bigl{\|}f_{t}^{(l)}(\cdot,t)\bigr{\|}_{\cal X}\in L_{p}({\cal I}),~{}0\leq l\leq m\big{\}}.$$

$$\|f\|_{C^{m}(\cal I;\cal X)}:=\max_{0\leq l\leq m}\sup_{t\in\cal I}\,\bigl{\|}f_{t}^{(l)}(\cdot,t)\bigr{\|}_{\cal X}.$$

$$\begin{array}[]{c}u_{t}+k(t)~{}{}_{0}^{R}D_{t}^{1-\alpha(t)}u+\mathcal{L}u=f(\bm{x},t),~{}~{}(\bm{x},t)\in\Omega\times(0,T];\\[3.61371pt] \displaystyle u(\bm{x},0)=u_{0}(\bm{x}),~{}\bm{x}\in\Omega;\quad u(\bm{x},t)=0,~{}(\bm{x},t)\in\partial\Omega\times[0,T].\end{array}$$

$$\alpha ,k\in C[0,T],k_m:=\underset{t\in [0,T]}{\min }k(t)>0;$$

$$0<{\alpha }_m:=\underset{t\in [0,T]}{\inf }\alpha (t)\le \alpha (t)\le 1,t\in [0,T],$$

$$\displaystyle\displaystyle\lim_{t\rightarrow 0^{+}}\big{(}\alpha(t)-\alpha(0)\big{)}\ln t=0;$$

$$0<K_m\le {\mathbf{\xi }}^T\mathbf{K}\mathbf{\xi }\le K_M<\infty ,\forall \mathbf{\xi }\in {\mathbb{R}}^d,|\mathbf{\xi }|=1,$$

$$k_{ij}\in C^1(\Omega ),1\le i,j\le d.$$

$$\displaystyle{}_{a}^{R}D_{t}^{1-\alpha(t)}g(t):=\bigg{[}\frac{1}{\Gamma\big{(}\alpha(t)\big{)}}\frac{d}{d\xi}\int_{a}^{\xi}\frac{g(s)}{(\xi-s)^{1-\alpha(t)}}ds\bigg{]}\bigg{|}_{\xi=t}.$$

$${}_a{I}_t^{\alpha (t)}g(t):=\frac{1}{\Gamma (\alpha (t))}{\int }_a^t\frac{g(s)}{(t-s{)}^{1-\alpha (t)}}ds,{}_a^C{D}_t^{1-\alpha (t)}g(t)={}_a{I}_t^{\alpha (t)}{g}^{\prime }(t).$$

$${}_0^R{D}_t^{1-\alpha (t)}g(t)={}_0^C{D}_t^{1-\alpha (t)}g(t)+\frac{g(0)t^{\alpha (t)-1}}{\Gamma (\alpha (t))},t\in (0,T].$$

$$\frac{\partial u}{\partial t}+k(t){}_0^C{D}_t^{1-\alpha (t)}u+\mathcal{L}u=-\frac{k(t)u_0(\mathbf{x})t^{\alpha (t)-1}}{\Gamma (\alpha (t))}+f(\mathbf{x},t)$$

$$rT\frac{(n-1)^{r-1}}{N^{r}}\leq\tau_{n}=T\Big{(}\frac{n}{N}\Big{)}^{r}-T\Big{(}\frac{n-1}{N}\Big{)}^{r}\leq rT\frac{n^{r-1}}{N^{r}},\quad r\geq 1,\quad 1\leq n\leq N.$$

$$\displaystyle\big{(}\bm{K}(\cdot)\nabla(g-\Pi g),\nabla\chi\big{)}=0,~{}~{}\forall\chi\in S_{h},~{}\mbox{for}~{}g\in H_{0}^{1}(\Omega),$$

$$\Vert g-\Pi g{\Vert }_{L_2(\Omega )}\le Q{h}^2\Vert g{\Vert }_{H^2(\Omega )},\forall g\in H^2(\Omega )\cap H_0^1(\Omega ).$$

$$\begin{array}[]{rl}\displaystyle u_{t}(\bm{x},t_{n})&\displaystyle=\delta_{\tau_{n}}u(\bm{x},t_{n})+r_{1,n}\\[7.22743pt] &\displaystyle:=\frac{u(\bm{x},t_{n})-u(\bm{x},t_{n-1})}{\tau_{n}}+\frac{1}{\tau_{n}}\int_{t_{n-1}}^{t_{n}}u_{tt}(\bm{x},t)(t-t_{n-1})dt,\\[10.84006pt] \displaystyle{}_{0}^{C}D_{t}^{1-\alpha(t_{n})}&\displaystyle u(\bm{x},t_{n})\\ &\displaystyle=\frac{1}{\Gamma(\alpha(t_{n}))}\sum_{k=1}^{n}\Big{[}\int_{t_{k-1}}^{t_{k}}\frac{\delta_{\tau_{k}}u(\bm{x},t_{k})}{(t_{n}-t)^{1-\alpha(t_{n})}}dt+\int_{t_{k-1}}^{t_{k}}\frac{u_{t}-\delta_{\tau_{k}}u(\bm{x},t_{k})}{(t_{n}-t)^{1-\alpha(t_{n})}}dt\Big{]}\\[12.64746pt] &\displaystyle=:\delta_{\tau}^{1-\alpha(t_{n})}u(\bm{x},t_{n})+r_{2,n},\\ \end{array}$$

$$\begin{array}[]{l}\displaystyle\delta_{\tau}^{1-\alpha(t_{n})}u(\bm{x},t_{n})\\[7.22743pt] \displaystyle\quad\quad:=\frac{1}{\Gamma(1+\alpha(t_{n}))}\sum_{k=1}^{n}\Big{[}\big{(}t_{n}-t_{k-1}\big{)}^{\alpha(t_{n})}-\big{(}t_{n}-t_{k}\big{)}^{\alpha(t_{n})}\Big{]}\delta_{\tau_{k}}u(\bm{x},t_{k})\\[10.84006pt] \displaystyle\quad\quad=\frac{1}{\Gamma(1+\alpha(t_{n}))}\sum_{k=1}^{n}b^{n}_{k}\big{(}u(\bm{x},t_{k})-u(\bm{x},t_{k-1})\big{)},\\[10.84006pt] r_{2,n}\displaystyle:=\frac{1}{\Gamma(\alpha(t_{n}))}\sum_{k=1}^{n}\int_{t_{k-1}}^{t_{k}}\frac{u_{t}-\delta_{\tau_{k}}u(\bm{x},t_{k})}{(t_{n}-t)^{1-\alpha(t_{n})}}dt\\[10.84006pt] \displaystyle\qquad=\frac{1}{\Gamma(\alpha(t_{n}))}\sum_{k=1}^{n}\int_{t_{k-1}}^{t_{k}}\frac{1}{\tau_{k}(t_{n}-t)^{1-\alpha(t_{n})}}\Big{[}\int_{t_{k-1}}^{t_{k}}\int_{s}^{t}u_{tt}(\bm{x},\theta)d\theta ds\Big{]}dt,\end{array}$$

$$b_k^n:=\frac{(t_n-t_{k-1}{)}^{\alpha (t_n)}-(t_n-t_k{)}^{\alpha (t_n)}}{{\tau }_k},1\le k\le n\le N$$

$$\left\{\begin{array}[]{l}\tau_{n}^{\alpha(t_{n})-1}=b^{n}_{n}>b^{n}_{n-1}>\cdots>b^{n}_{k}>\ldots b^{n}_{1}>0,\\[3.61371pt] \displaystyle\alpha(t_{n})(t_{n}-t_{k-1})^{\alpha(t_{n})-1}\leq b^{n}_{k}\leq\alpha(t_{n})(t_{n}-t_{k})^{\alpha(t_{n})-1}.\end{array}\right.$$

$$\begin{array}[]{l}\displaystyle(\delta_{\tau_{n}}u^{n},\chi)+\big{(}\bm{K}(\cdot)\nabla u^{n},\nabla\chi\big{)}=-k(t_{n})(\delta_{\tau}^{1-\alpha(t_{n})}u^{n},\chi)-\frac{k(t_{n})t_{n}^{\alpha(t_{n})-1}}{\Gamma(\alpha(t_{n}))}(u_{0},\chi)\\[10.84006pt] \quad\quad\displaystyle+(f(\cdot,t_{n}),\chi)-\big{(}k(t_{n})r_{2,n}+r_{1,n},\chi\big{)},~{}~{}\chi\in H_{0}^{1}(\Omega),~{}~{}n=1,\cdots,N.\end{array}$$

$$\begin{array}[]{l}\displaystyle(\delta_{\tau_{n}}u_{h}^{n},\chi)+\big{(}\bm{K}(\cdot)\nabla u_{h}^{n},\nabla\chi\big{)}=-k(t_{n})(\delta_{\tau}^{1-\alpha(t_{n})}u_{h}^{n},\chi)-\frac{k(t_{n})t_{n}^{\alpha(t_{n})-1}}{\Gamma(\alpha(t_{n}))}(u_{0},\chi)\\[10.84006pt] \quad\quad\displaystyle+(f(\cdot,t_{n}),\chi),~{}~{}\chi\in S_{h},~{}~{}n=1,\cdots,N.\end{array}$$

$$\begin{array}[]{l}\displaystyle\|r_{1}\|_{\hat{L}_{\infty}(0,T;L_{2})}:=\max_{1\leq n\leq N}\|r_{1,n}\|_{L_{2}}\leq QQ_{0}N^{-1},~{}~{}\|r_{2}\|_{\hat{L}_{\infty}(0,T;L_{2})}\leq QQ_{0}N^{-1},\\[7.22743pt] \displaystyle\quad\quad Q_{0}:=\big{(}\|u_{0}\|_{\check{H}^{4}}+\|f\|_{H^{1}(0,T;\check{H}^{2})}+\|f\|_{H^{2}(0,T;L_{2})}\big{)}.\end{array}$$

$$\begin{array}[]{l}\displaystyle\|r_{1,n}\|_{\hat{L}_{\infty}(0,T;L_{2})}\leq QQ_{0}n^{-1-r(1-\alpha(0))}N^{r(1-\alpha(0))},\\[7.22743pt] \displaystyle\|r_{2,n}\|_{\hat{L}_{\infty}(0,T;L_{2})}\leq QQ_{0}\Big{(}\hat{\delta}_{n,1}N^{r(1-\alpha(0)-\alpha(1))}\\[7.22743pt] \displaystyle\qquad~{}~{}\qquad\qquad\qquad+(1-\hat{\delta}_{n,1})N^{r(1-\alpha(0)-\alpha(t_{n}))}n^{-1-r(1-\alpha(0)-\alpha(t_{n}))}\Big{)}.\end{array}$$

$$\displaystyle r_{3,n}:=\delta_{\tau_{n}}\big{(}u(\bm{x},t_{n})-\Pi u(\bm{x},t_{n})\big{)},~{}~{}r_{4,n}:=\delta_{\tau}^{1-\alpha(t_{n})}\big{(}u(\bm{x},t_{n})-\Pi u(\bm{x},t_{n})\big{)},$$

$$\Vert r_{3,n}{\Vert }_{{\hat{L}}_{\infty }(0,T;L_2)}+\Vert r_{4,n}{\Vert }_{{\hat{L}}_{\infty }(0,T;L_2)}\le Q{Q}_1{h}^2,Q_1:=\Vert u_0{\Vert }_{{\check{H}}^4}+\Vert f{\Vert }_{H^1(0,T;{\check{H}}^2)}$$

$$\begin{array}[]{l}\displaystyle\|r_{3,n}\|_{\hat{L}_{\infty}(0,T;L_{2})}\leq QQ_{1}h^{2}\big{(}\delta_{n,1}+(1-\delta_{n,1})n^{r(\alpha(0)-1)}\big{)}N^{-r(\alpha(0)-1)},\\[10.84006pt] \displaystyle\|r_{4,n}\|_{\hat{L}_{\infty}(0,T;L_{2})}\\[10.84006pt] \displaystyle\qquad\leq QQ_{1}h^{2}\big{(}\delta_{n,1}+(1-\delta_{n,1})n^{r(\alpha(0)+\alpha(t_{n})-1)}\big{)}N^{-r(\alpha(0)+\alpha(t_{n})-1)}.\end{array}$$

$$\begin{array}[]{rl}\displaystyle\|r_{3,n}\|&\displaystyle=\tau_{n}^{-1}\|(I-\Pi)(u^{n}-u^{n-1})\|\\[7.22743pt] &\displaystyle\leq Qh^{2}\tau_{n}^{-1}\|u^{n}-u^{n-1}\|_{H^{2}}\leq Qh^{2}\|u\|_{W^{1}_{\infty}(0,T;H^{2})}\end{array}$$

$$\begin{array}[]{l}\displaystyle\|r_{4,n}\|=\frac{1}{\Gamma(1+\alpha(t_{n}))}\bigg{\|}\sum_{k=1}^{n}b^{n}_{k}(I-\Pi)(u^{k}-u^{k-1})\bigg{\|}\\[10.84006pt] \displaystyle~{}~{}~{}~{}~{}\,~{}~{}\leq Qh^{2}\sum_{k=1}^{n}b^{n}_{k}\|u^{k}-u^{k-1}\|_{H^{2}}\leq Qh^{2}\|u\|_{W^{1}_{\infty}(0,T;H^{2})}\sum_{k=1}^{n}b^{n}_{k}\tau_{k}\\[10.84006pt] \displaystyle~{}~{}~{}\,~{}~{}~{}~{}\leq Qh^{2}\|u\|_{W^{1}_{\infty}(0,T;H^{2})},\end{array}$$

$$\begin{array}[]{rl}\displaystyle\|r_{3,n}\|&\displaystyle=\frac{1}{\tau_{n}}\|(I-\Pi)(u^{n}-u^{n-1})\|\\[10.84006pt] &\displaystyle\leq\frac{Qh^{2}}{\tau_{n}}\|u^{n}-u^{n-1}\|_{H^{2}}=\frac{Qh^{2}}{\tau_{n}}\bigg{\|}\int_{t_{n-1}}^{t_{n}}u_{t}dt\bigg{\|}_{H^{2}}\leq\frac{QQ_{1}h^{2}}{\tau_{n}}\int_{t_{n-1}}^{t_{n}}t^{\alpha(0)-1}dt\\[10.84006pt] &\displaystyle\leq\left\{\begin{array}[]{l}\displaystyle QQ_{1}h^{2}t_{1}^{\alpha(0)-1}=\frac{QQ_{1}h^{2}}{N^{r(\alpha(0)-1)}},~{}n=1,\\[7.22743pt] \displaystyle QQ_{1}h^{2}t_{n-1}^{\alpha(0)-1}\leq\frac{QQ_{1}h^{2}(n-1)^{r(\alpha(0)-1)}}{N^{r(\alpha(0)-1)}},~{}n>1.\end{array}\right.\end{array}$$

$$\begin{array}[]{rl}\displaystyle\|r_{4,n}\|&\displaystyle=\frac{1}{\Gamma(1+\alpha(t_{n}))}\bigg{\|}\sum_{k=1}^{n}b^{n}_{k}(I-\Pi)(u^{k}-u^{k-1})\bigg{\|}\\[10.84006pt] &\displaystyle\leq Qh^{2}\sum_{k=1}^{n}b^{n}_{k}\|u^{k}-u^{k-1}\|_{H^{2}}\leq Qh^{2}\sum_{k=1}^{n}b^{n}_{k}\bigg{\|}\int_{t_{k-1}}^{t_{k}}u_{t}dt\bigg{\|}_{H^{2}}\\[14.45377pt] &\displaystyle\leq QQ_{1}h^{2}\sum_{k=1}^{n}b^{n}_{k}\int_{t_{k-1}}^{t_{k}}t^{\alpha(0)-1}dt:=QQ_{1}h^{2}\sum_{k=1}^{n}J_{n,k}.\end{array}$$

$$J_{1,1}=\frac{b_1^1{t}_1^{\alpha (0)}}{\alpha (0)}=\frac{t_1^{\alpha (1)+\alpha (0)-1}}{\alpha (0)}\le Q{N}^{-r(\alpha (1)+\alpha (0)-1)}.$$

$$\begin{array}[]{rl}\displaystyle|J_{n,1}|&\displaystyle=\frac{b^{n}_{1}t_{1}^{\alpha(0)}}{\alpha(0)}=\frac{(t_{n}^{\alpha(t_{n})}-(t_{n}-t_{1})^{\alpha(t_{n})})t_{1}^{\alpha(0)}}{\tau_{1}\alpha(0)}\\[10.84006pt] &\displaystyle\leq Q(t_{n}-t_{1})^{\alpha(t_{n})-1}t_{1}^{\alpha(0)}=Q\bigg{(}\frac{n^{r}-1}{N^{r}}\bigg{)}^{\alpha(t_{n})-1}\frac{1}{N^{r\alpha(0)}}\\[10.84006pt] &\displaystyle\leq Q\frac{n^{r(\alpha(t_{n})-1)}}{N^{r(\alpha(0)+\alpha(t_{n})-1)}},\end{array}$$

Full text

1 Introduction

Fractional partial differential equations (FPDEs) have shown to provide adequate descriptions for the challenging phenomena such as the anomalous diffusive transport and the memory effect MeeSik ; MetKla00 ; Pod . For instance, in the diffusive transport in heterogeneous porous media, a large amount of particles may get absorbed to the surface of the rock. Thus, the travel time of the adsorbed particles may deviate from that of the particles in the bulk phase ZhoStr , which leads to a subdiffusive transport that can be modeled by a time-fractional diffusion PDE (tFPDE) MeeSik ; MetKla00 .

Extensive mathematical and numerical analysis of FPDEs has been conducted Die ; DieFor ; ErvRoo05 ; KilSri ; LeMclMus ; Luc ; Pod ; SakYam ; StyOriGra , and it is gradually getting clear that the FPDEs introduce mathematical issues that are not common in the context of integer-order PDEs. For instance, the smoothness of the coefficients and right-hand side of a linear elliptic or parabolic fractional PDE in one space dimension cannot ensure the smoothness of its solutions JinLazPas ; StyOriGra ; WanYanZhu ; WanZha . Hence, many error estimates in the literature that were proved under full regularity assumptions of the true solutions are inappropriate.

Variable-order tFPDEs, in which the order of the fractional derivatives varies in time as $t\rightarrow 0$ to accommodate the impact of the local initial condition at time $t=0$, should be a natural candidate to eliminate the nonphysical singularity of the solutions to (constant-order) tFPDEs and open up opportunities for modeling multiphysics phenomena from nonlocal to local dynamics and vice versa LiWan ; PanPer ; SunCheChe ; ZenZhaKar ; ZhuLiu .

Due to the difficulties of solving variable-order tFPDEs analytically, several numerical methods have been developed (see e.g. ZhuLiu ; ZenZhaKar ) under certain smoothness assumptions of the solutions. It was shown in StyOriGra that the first order time derivatives of solutions to the $\alpha$-order time-fractional diffusion equations (tFDEs) exhibit the singularity of $\mathcal{O}(t^{\alpha-1})$ at the initial time $t=0$, which leads to a sub-optimal convergence of the fully discrete finite difference method. It was also proved that by using the graded temporal mesh with a proper chosen mesh grading parameter according to the singularity of the solutions, the optimal convergence rate of the proposed finite difference method can be recovered.

Recently, the wellposedness of a variable-order tFDE model and the regularity of its solutions were studied in WanZheJMAA . In particular, the solutions have full regularity like those to the integer-order tFDEs if the variable order has an integer limit at $t=0$ or exhibit singularity at $t=0$ like in the case of the constant-order tFDEs if the variable order has a non-integer value at time $t=0$. Based on these theoretical results, we present a first order time-discretized finite element method for this variable-order tFDE model. When the variable order smoothly transit the fractional order model to the integer order ones near the initial time, the solutions have full regularity and the optimal convergence is proved under the uniform temporal mesh. Otherwise, a graded mesh with a properly chosen mesh grading parameter in terms of the singularity of the solutions at the initial time is applied to recover the optimal convergence rate.

The rest of the papers are organized as follows: In §2 a variable-order tFDE model and the auxiliary results to be used subsequently were presented. In §3 a first order time-discretized finite element method was developed for the proposed model and we proved the corresponding optimal error estimates under the uniform or graded temporal mesh in terms of the regularity of the solutions in §4. Several numerical experiments were presented in §5 to demonstrate the theoretical analysis.

2 Model problem and preliminaries

Let $m\in\mathbb{N}$, $1\leq p\leq\infty$, $\Omega\subset\mathbb{R}^{d}$ ($d=1,2,3$) be a simply-connected bounded domain with smooth boundary $\partial\Omega$ and ${\cal I}\subset[0,\infty)$ be a bounded interval. Let $L_{p}(\Omega)$ be the spaces of the $p$-th power Lebesgue integrable functions on $\Omega$ and $W^{m}_{p}(\Omega)$ be the Sobolev spaces of functions with derivatives of order up to $m$ in $L_{p}(\Omega)$. Let $H^{m}(\Omega):=W^{m}_{2}(\Omega)$ and $H^{m}_{0}(\Omega)$ be the completion of $C^{\infty}_{0}(\Omega)$, the space of infinitely many time differentiable functions with compact support in $\Omega$, in $H^{m}(\Omega)$ AdaFou . For the case of non-integer order $s$, the fractional Sobolev spaces $H^{s}(\Omega)$ are defined by interpolation, see AdaFou . Furthermore, for the Banach space ${\cal X}$, we introduce the Sobolev spaces involving time AdaFou ; Eva

[TABLE]

In particular, $W^{0}_{p}({\cal I};{\cal X})=L_{p}({\cal I};{\cal X})$ for $1\leq p\leq\infty$. We also let $C^{m}({\cal I};{\cal X})$ be the spaces of functions with continuous derivatives up to order $m$ on ${\cal I}$ equipped with the norm

[TABLE]

In this paper we study the initial-boundary value problem of a variable-order linear tFDE

[TABLE]

Here $u_{t}$ refers to the first-order partial derivative in time, ${\bm{x}}:=(x_{1},\cdots,x_{d})$, \mathcal{L}:=-\nabla\cdot\big{(}\bm{K}(\bm{x})\nabla\big{)} with \nabla:=(\partial/\partial x_{1},\cdots,\partial/\partial x_{d}\big{)}^{T} and ${\bm{K}}({\bm{x}}):=(k_{ij}(\bm{x}))_{i,j=1}^{d}$ the diffusion tensor. We make the following assumptions throughout the paper.

Assumption A.

[TABLE]

The variable-order Riemann-Liouville fractional derivative is defined by ZenZhaKar ; ZhuLiu

[TABLE]

Moreover, we also use the variable-order fractional integral operator ${}_{a}I_{t}^{\alpha(t)}$ and the Caputo fractional differential operator ${}_{a}^{C}D_{t}^{\alpha(t)}$ ZenZhaKar ; ZhuLiu

[TABLE]

Remark 1

The constant-order analogue of the proposed model is known as the mobile-immobile time-fractional diffusion equations, see e.g., LiuLiu .

The relation between the Riemann-Liouville and Caputo fractional derivatives Die ; ErvRoo05 ; Pod was extended to the variable-order analogues ZhuLiu .

Lemma 1

Let $g\in W^{1}_{1}(0,T)$. Then

[TABLE]

In this paper we use $Q$ to denote generic positive constants that may assume different values at different occurrences. For convenience, we may drop the subscript $L_{2}$ in $(\cdot,\cdot)_{L_{2}}$ and $\|\cdot\|_{L_{2}}$ as well as the notation $\Omega$ in the Sobolev spaces and norms, and abbreviate $W^{m}_{p}(0,T;{\cal X})$ and $W^{m}_{p}({\cal X})$, when no confusion occurs.

3 Fully discrete finite element method for variable-order tFDEs

By Lemma 1 and Theorem 6.1, the Riemann-Liouville variable-order tFDE (1) and the following Caputo variable-order tFDE

[TABLE]

coincide. So we will develop and analyze the corresponding finite element schemes for (3).

Let $t_{n}:=T(n/N)^{r}$, $0\leq n\leq N$ be a partition of $[0,T]$, which forms a graded mesh when $r>1$ and reduces to a uniform partition for $r=1$. Applying the mean-value theorem we bound $\tau_{n}:=t_{n}-t_{n-1}$ by

[TABLE]

Define $\Omega_{e}$ a quasi-uniform partition of $\Omega$ with parameter $h$ and $S_{h}$ the space of piece-wise linear functions on $\Omega$ with compact support. The Ritz projection $\Pi:H^{1}_{0}(\Omega)\rightarrow S_{h}$ defined by

[TABLE]

has the following approximation property Tho

[TABLE]

We discretize $u_{t}$ and ${}_{0}^{C}D_{t}^{1-\alpha(t)}u$ at $t=t_{n}$, $1\leq n\leq N$ by

[TABLE]

with

[TABLE]

where

[TABLE]

has the following properties StyOriGra

[TABLE]

Let $u^{n}:=u(\bm{x},t_{n})$. We plug (7) into (3), multiply $\chi\in H_{0}^{1}(\Omega)$ on both sides and integrate the resulting equation on $\Omega$ to get the weak formulation of (3)

[TABLE]

We drop the truncation error terms to obtain a first order time-discretized finite element scheme for (3): find $u_{h}^{n}\in S_{h}$ such that

[TABLE]

4 Convergence estimates of the finite element approximations

We prove the optimal error estimates of the finite element approximations under the uniform or graded mesh in terms of the regularity of the solutions to the proposed model.

4.1 Analysis of truncation errors

The estimates of the local truncation errors $r^{n}_{1}$ and $r^{n}_{2}$ in (7) and (8) are given in the following theorem.

Theorem 4.1

Suppose that $u_{0}\in\check{H}^{4}$ and $f\in H^{1}(0,T;\check{H}^{2})\cap H^{2}(0,T;L_{2})$. For the case of $\alpha(0)=1$, $\alpha^{\prime}(0)=0$ and $\lim_{t\rightarrow 0^{+}}\alpha^{\prime}(t)\ln t$ is finite, the following estimate holds under the uniform temporal partition

[TABLE]

Otherwise, the following estimates hold under the graded mesh with $r\alpha(0)>1$

[TABLE]

Here $\hat{\delta}_{m,n}$ is the Kronecker delta function.

Proof

The proof is exactly the same as that of Theorem 7 in WanZhe with $\alpha_{*}$ in that theorem replaced by $\alpha(0)$ according to Theorem 6.2.

We bound another two truncation terms for $1\leq n\leq N$

[TABLE]

for the convenience of the convergence estimates.

Theorem 4.2

Suppose $u_{0}\in\check{H}^{4}$, $f\in H^{1}(0,T;\check{H}^{2})$ and the Assumption A holds. Then the following estimates hold under the uniform temporal partition for the case of $\alpha(0)=1$, $\alpha^{\prime}(0)=0$ and $\lim_{t\rightarrow 0^{+}}\alpha^{\prime}(t)\ln t$ is finite

[TABLE]

and under the graded mesh otherwise

[TABLE]

Proof

When $\alpha(0)=1$, $\alpha^{\prime}(0)=0$ and $\lim_{t\rightarrow 0^{+}}\alpha^{\prime}(t)\ln t$ is finite, u\in C^{1}\big{(}[0,T];H^{2}(\Omega)\big{)} by Theorem 6.1 so a uniform partition of $[0,T]$ suffices. We apply (6) to obtain

[TABLE]

and

[TABLE]

where $I$ refers to the identity operator. Then an application of Theorem 6.1 leads to (15).

For other cases, we only need to consider the case that $\alpha(0)<1$. The graded mesh with mesh grading $r$ will be used to capture the singularity of the solutions at the initial time. By (25) and the mean-value theorem we bound $r_{3,n}$ by

[TABLE]

We remain to bound $r_{4,n}$, which requires a careful argument. From (18) we have

[TABLE]

When $n=1$, $J_{1,1}$ can be bounded by

[TABLE]

For $n>1$, we first bound $J_{n,1}$ and $J_{n,n}$ by (4) and the mean-value theorem

[TABLE]

[TABLE]

We remain to consider the case $n\geq 3$ since the estimates (19) and (20) have covered the case $n=2$. By (4) and the mean-value theorem we obtain

[TABLE]

and

[TABLE]

We summarize the above estimates to finish the proof.

4.2 Convergence estimates of the finite element approximations

We prove the optimal error estimate of the fully discrete finite element method (11) by the following theorem.

Theorem 4.3

Suppose that $u_{0}\in\check{H}^{4}$ and $f\in H^{1}(0,T;\check{H}^{2})\cap H^{2}(0,T;L_{2})$. We set $r=1$ for the case of $\alpha(0)=1$, $\alpha^{\prime}(0)=0$ and $\lim_{t\rightarrow 0^{+}}\alpha^{\prime}(t)\ln t$ is finite and $r=2/\alpha(0)$ otherwise. Then the following optimal order error estimate holds

[TABLE]

Here $Q=Q(T,\alpha_{m},\|k\|_{C[0,T]})$.

Proof

We split the error by $e^{n}:=u_{h}^{n}-u^{n}=\xi^{n}+\eta^{n}$ where $\xi^{n}:=u_{h}^{n}-\Pi u^{n}$ and $\eta^{n}:=\Pi u^{n}-u^{n}$. The estimate of $\eta^{n}$ is given by (6) so we remain to bound $\xi^{n}$. We subtract (11) from (10) with $\chi=\xi^{n}$ and apply (5) and (14) into the resulting equation to obtain the following error equation in terms of $\xi^{n}$

[TABLE]

We rearrange $\delta_{\tau}^{1-\alpha(t_{n})}$ by

[TABLE]

and apply $u_{h}^{0}:=\Pi u^{0}$ to reformulate (21) as

[TABLE]

from which we use (9) to obtain

[TABLE]

where $A_{n}:=k(t_{n})/\Gamma(1+\alpha(t_{n}))$.

We turn to evaluate the truncation error terms on the right-hand side of (22). In the case $\alpha(0)=1$, $\alpha^{\prime}(0)=0$ and $\lim_{t\rightarrow 0^{+}}\alpha^{\prime}(t)\ln t$ is finite, the uniform partition is applied and by (12) and (15) we directly obtain

[TABLE]

with $Q_{0}$ defined in (12). Otherwise, the graded mesh with $r=2/\alpha(0)$ is chosen. We use (13), (16) and the fact that the bound of $r_{3,n}$ dominates that of $r_{4,n}$ (see Theorem 4.2) to obtain

[TABLE]

and

[TABLE]

Therefore, in any case of $\alpha(t)$, we obtain the following estimates of the truncation errors under the appropriate temporal partition

[TABLE]

for some fixed constant $Q^{\prime}$. We then prove the convergence estimates by mathematical induction. Applying (23) to (22) with $n=1$ yields

[TABLE]

Assume

[TABLE]

Plugging (23) and (24) with $2\leq m\leq n-1$ into (22) leads to

[TABLE]

in which we divide $1+A_{n}\tau_{n}^{\alpha(t_{n})}$ on both sides to obtain (24) for $m=n$ and thus for any $m\geq 2$ by mathematical induction. Then the proof is finished by applying $m\tau_{N}\leq rTm/N\leq rT$ into (24).

5 Numerical experiments

We substantiate the analysis numerically by investigating the impact of $\alpha(t)$ on the convergence rate of the fully discrete finite element scheme (11).

Let $(\bm{x},t)\in(0,1)^{3}\times[0,1]$, $k(t)=1$, $\bm{K}=\text{diag}(0.001,0.001,0.001)$ and

[TABLE]

with

[TABLE]

which satisfies $\alpha^{\prime}(0)=0$ and $\lim_{t\rightarrow 0^{+}}\alpha^{\prime}(t)\ln t=0$ and $f$ evaluated accordingly. We measure the convergence rates $\kappa$ and $\gamma$ such that

[TABLE]

We select the uniform partition on space domain and the uniform temporal mesh is used for the case of $\alpha(0)=1$ (i.e., $u_{tt}$ is continuous on $[0,T]$). For the case of $\alpha(0)<1$ (i.e., the solutions exhibit singularity at $t=0$), both uniform and graded meshes with $r=2/\alpha(0)$ for time are applied. We present results in Table 1 and 2, which reveal that the scheme (11) with a uniform mesh has an optimal-order convergence rate for the case of smooth solutions (i.e., $\alpha(0)=1$), but only a sub-optimal order for the case of $\alpha(0)<1$. Instead, the scheme (11) with the temporal graded mesh of $r=2/\alpha(0)$ achieves an optimal-order convergence rate. These results coincide with Theorem 4.3.

6 Appendix: Wellposedness of the variable-order tFDE and regularity of its solutions

It is known CouHil ; Eva that the eigenfunctions $\{\phi_{i}\}_{i=1}^{\infty}$ of the Sturm-Liouville problem

[TABLE]

form an orthonormal basis in $L_{2}(\Omega)$. The eigenvalues $\{\lambda_{i}\}_{i=1}^{\infty}$ are positive and form a nondecreasing sequence that tend to $\infty$ with $i$. We use the theory of sectorial operators to define the fractional Sobolev spaces SakYam ; Tho

[TABLE]

with the norm being defined by \|v\|_{\check{H}^{\gamma}}:=\big{(}\|v\|^{2}+|v|_{\check{H}^{\gamma}}^{2}\big{)}^{1/2}. Note that $\check{H}^{\gamma}(\Omega)$ is a subspace of the fractional Sobolev space $H^{\gamma}(\Omega)$ characterized by AdaFou ; SakYam ; Tho

[TABLE]

and the seminorms $|v|_{\check{H}^{\gamma}}$ and $|v|_{{H}^{\gamma}}$ are equivalent in $\check{H}^{\gamma}$.

Then using the approaches of variable seperation, the wellposedness of the model (1) and the regularity estimates of its solutions $u(x)$ are proved by the following theorems WanZheJMAA .

Theorem 6.1

Suppose that $u_{0}\in\check{H}^{\gamma+2}$, $f\in H^{1}(0,T;\check{H}^{\gamma})$ for some $d/2<\gamma\in\mathbb{R}^{+}$ and Assumption A holds. If $\alpha(0)=1$ then problem (1) has a unique solution u\in C^{1}\big{(}[0,T];\check{H}^{\gamma}\big{)} with the following stability estimates for any $0\leq s\leq\gamma$

[TABLE]

If $\alpha(0)<1$ problem (1) has a unique solution u\in C\big{(}[0,T];\check{H}^{\gamma}\big{)}\cap C^{1}\big{(}(0,T];\check{H}^{\gamma}\big{)} with the stability estimate

[TABLE]

Here Q=Q\big{(}\alpha_{m},\|k\|_{C[0,T]},T\big{)}.

Theorem 6.2

Suppose that $u_{0}\in\check{H}^{s+4}$, $f\in H^{1}(0,T;\check{H}^{s+2})\cap H^{2}(0,T;\check{H}^{s})$ for $s\geq 0$, and that $\alpha,k\in C^{1}[0,T]$ and (2) holds. Then the following conclusions hold:

Case 1.

If $\alpha(0)=1$, $\alpha^{\prime}(0)=0$ and $\lim_{t\rightarrow 0^{+}}\alpha^{\prime}(t)\ln t$ is finite, then $u\in C^{2}([0,T];$ $\check{H}^{s})$ and

[TABLE]

Case 2.

If $\alpha(0)=1$ but $\alpha^{\prime}(0)\neq 0$ or $\lim_{t\rightarrow 0^{+}}\alpha^{\prime}(t)\ln t$ is not finite, then $u\in C^{1}([0,T];\check{H}^{s})\cap C^{2}((0,T];\check{H}^{s})$ and for any $0<\varepsilon\ll 1$

[TABLE]

Case 3.

If $\alpha(0)<1$, then $u\in C^{2}((0,T];\check{H}^{s}(\Omega))$ and for any $0<\varepsilon\ll 1$

[TABLE]

Here Q=Q\big{(}\alpha_{m},\|k\|_{C^{1}[0,T]},T\big{)}.

Acknowledgements

This work was funded by the OSD/ARO MURI Grant W911NF-15-1-0562 and the National Science Foundation under Grant DMS-1620194.