High order approximation to Caputo derivative on graded mesh and time-fractional diffusion equation for non-smooth solutions

Shweta Kumari Abhishek Kumar Singh Vaibhav Mehandiratta Mani Mehra

Abstract

In this paper, a high-order approximation to Caputo-type time-fractional diffusion equations involving an initial-time singularity of the solution is proposed. At first, we employ a numerical algorithm based on the Lagrange polynomial interpolation to approximate the Caputo derivative on the non-uniform mesh. Then truncation error rate and the optimal grading constant of the approximation on a graded mesh are obtained as $\min\{4-\alpha,r\alpha\}$ and $\frac{4-\alpha}{\alpha}$ , respectively, where $\alpha\in(0,1)$ is the order of fractional derivative and $r\geq 1$ is the mesh grading parameter. Using this new approximation, a difference scheme for the Caputo-type time-fractional diffusion equation on graded temporal mesh is formulated. The scheme proves to be uniquely solvable for general $r$ . Then we derive the unconditional stability of the scheme on uniform mesh. The convergence of the scheme, in particular for $r=1$ , is analyzed for non-smooth solutions and concluded for smooth solutions. Finally, the accuracy of the scheme is verified by analyzing the error through a few numerical examples.

Keywords. Time-fractional diffusion equation, Caputo derivative, non-smooth solution, graded mesh, stability.

1 Introduction

The classical integer-order differential equations are incapable of framing several daily life occurrences due to their complex structure and dependency of the physical system on previous impressions. However, the fractional operators are known mainly for their nonlocal behavior and past memory storage. Hence, the involvement of differential equations with fractional order, namely the fractional differential equations (FDEs), is crucial in molding such phenomena. They can be particularised by generalizing the order of the integer-order derivative to some arbitrary value. Although both calculus and fractional calculus trace their roots back to the same time period, i.e. the late 17th century, the latter did not receive ample attention as no practical usage of it was apparent at the time. But in the past few years, the study of FDEs upsurged due to their wide applications in the fields of science and engineering [7, 20]. FDEs are generally used to describe the dynamic systems that represent many natural processes in physics [4], fluid mechanics [15], biology [5, 26], system control [18, 17], wavelet [8, 24], network like structures [16], finance [19, 23] and various other engineering fields [9, 25, 22]. Among fractional partial differential equations (FPDEs), fractional diffusion equations are attracting strong interests because of their ability to well specify anomalous diffusion processes [14] in complex media since classical diffusion is not the best description for them.

This paper is concerned with developing a numerical scheme for the Caputo-type time-fractional diffusion equation (TFDE) for the following initial-boundary value problem (IBVP):

\textsubscript{C}D_{0,t}^{\alpha}u(x,t)-\varrho\partial_{x}^{2}u(x,t)=f(x,t),~{}~{}0<t\leq T,~{}0<x<L,\\

(1.1)

u(0,t)=u(L,t)=0,~{}~{}~{}0<t\leq T,\\

(1.2)

u(x,0)=\varphi(x),~{}~{}~{}0\leq x\leq L,

(1.3)

where $\textsubscript{C}D_{0,t}^{\alpha}$ denotes the left Caputo fractional derivative of order $\alpha\in(0,1)$ , with respect to time variable $t$ , $\partial_{x}^{2}$ is the spatial Laplacian operator, (1.2) denotes the Dirichlet-type boundary conditions, and (1.3) denotes the initial condition. The notation $\varrho$ $(>0)$ is a constant diffusion coefficient. Here, $f\in C(\bar{\Omega})$ and $\varphi\in C[0,L]$ , where $\bar{\Omega}=[0,L]\times[0,T]$ .

Recently, intense focus has been placed on the approximation to Caputo derivative as it models phenomena with productive memory functions that take account of past interactions and problems with nonlocal properties [10]. The classic $L1$ format to approximate the Caputo derivative on a uniform grid is derived by Lin and Xu in [13], where the integrand of the Caputo derivative is approximated by a linear interpolation polynomial in each time subinterval. For the subdiffusion equation with non-smooth data, the $L1$ scheme was analysed recently in [6]. Gao et al. [3] proposed the $L1-2$ format for Caputo fractional derivative on a uniform grid, i.e. a linear and quadratic interpolation polynomial is used to discretize the integrand in first and other time subintervals, respectively, with a convergence order of $(3-\alpha)$ . Other works for higher order numerical algorithms having $(3-\alpha)$ convergence rate are presented in [12, 1]. In [12], an approximation of order $(3-\alpha)$ was derived for the Caputo-type advection-diffusion equation. Alikhanov [1] gave the $L2-1_{\sigma}$ format to approximate the fractional derivatives on uniform grid, where $\sigma=1-\frac{\alpha}{2}$ is the offset parameter. A scheme of order $(4-\alpha)$ is developed by Cao et al. [2] for time fractional advection-diffusion equations having smooth solutions throughout the time interval on a uniform mesh.

Many researchers have obtained the numerical solution of the TFDE (1.1)-(1.3) with the assumption that the solution and its consecutive derivatives satisfy certain smooth conditions. Nevertheless, Stynes et al. [27] were the first to observe that the solution of the problem (1.1)-(1.3) usually has an initial-time singular behavior. In that case, however, the numerical schemes proposed over the uniform grid converge, but at the cost of a significantly reduced error convergence rate. So, they presented their work [27] for the time-fractional reaction-diffusion problem to deal with this initial-time singularity on a non-uniform grid with the ideal convergence order of $\min\{2-\alpha,r\alpha\}$ , where $\alpha\in(0,1)$ is the fractional order and $r\geq 1$ is the mesh grading parameter. They used the $L1$ format to discretize the Caputo fractional derivative on graded temporal meshes, and the diffusion term is discretized by the central difference approximation on a uniform spatial grid. Following a similar path, Qiao and Cheng [21] successfully got a truncation error of order $\min\{3-\alpha,r\alpha\}$ using the $L1-2$ type format to discretize the Caputo derivative on the graded mesh. But the examination of the stability and convergence of the scheme was missing.

Taking the aforementioned works as motivation, we attempt to develop an even higher order approximation of the Caputo derivative on graded mesh and use it to construct a difference scheme to TFDE (1.1) having non-smooth solutions. To achieve this, we apply a numerical algorithm to discretize the Caputo derivative on graded temporal mesh in which the linear, quadratic, and cubic interpolation polynomials are used to approximate the integrand of the Caputo derivative in the first, second, and third onwards time subintervals, respectively. The diffusion term is discretized by second order central difference scheme on a uniform spatial grid. This algorithm requires the solution to be at least fourth-order continuously differentiable in time. However, solutions usually do not satisfy such high regularity assumptions, which results in a low error convergence rate of the scheme. Hence, the grading of mesh is considered so that the initial-time singularity does not affect the overall convergence of the scheme. Surprisingly, there was a lack of literature on the convergence analysis of the scheme for TFDEs with smooth solutions as well. So, we also investigate the convergence order for TFDEs with smooth solutions on a uniform mesh.

The rest organization of this paper goes as follows: Section 2 contains the preliminary content required for the reader to understand the concepts of this paper. The higher order approximation of Caputo derivative on graded mesh has been constructed and thoroughly analyzed in Section 3. It also comprises the estimation of the local truncation error of the approximation. Section 4 builds a numerical scheme for the TFDE problem (1.1)-(1.3) using the approximation developed in Section 3. Further, we check the uniqueness of the solution of the scheme. We dedicate a subsection to study the unconditional stability and convergence of the scheme developed for $r=1$ , i.e. on uniform mesh and for both non-smooth as well as smooth solutions. In Section 5, a few numerical examples are displayed to verify the accuracy and efficiency of the scheme by studying the errors. Finally, Section 6 concludes the paper.

2 Preliminaries

In this part, we present some fundamental concepts that are prerequisites for a better understanding of our work in this paper.

Definition 1.

The left Caputo derivative with fractional order $\alpha>0$ of the function $f(t),\\ f\in C^{m}(a,b)$ , is defined as [11]

\textsubscript{C}D_{a,t}^{\alpha}f(t)=\frac{1}{\Gamma(m-\alpha)}\int_{a}^{t}(t-s)^{m-\alpha-1}f^{(m)}(s)ds,

(2.1)

where $m\in\mathbb{N}$ satisfying $m-1<\alpha\leq m$ and $\Gamma(.)$ is the Euler’s gamma function .

Definition 2.

Euler’s gamma function is denoted and defined as

\Gamma(z)=\int_{0}^{\infty}x^{z-1}e^{-x}dx,~{}~{}z\in\mathbb{C},~{}\Re(z)>0.

Here $\Gamma(z+1)=z\Gamma(z)$ , and in particular, $\Gamma(z)=(z-1)!$ for $z\in\mathbb{Z}^{+}$ .

Definition 3.

The linear interpolation of a function $f(t)$ on some interval $[t_{k-1},t_{k}],~{}k\geq 1,$ is defined as

I_{1k}f(t)=\frac{t-t_{k}}{t_{k-1}-t_{k}}f(t_{k-1})+\frac{t-t_{k-1}}{t_{k}-t_{k-1}}f(t_{k}).

(2.2)

The quadratic interpolation of a function $f(t)$ on some interval $[t_{k-1},t_{k}],~{}k\geq 2,$ is defined as

	$\displaystyle I_{2k}f(t)$	$\displaystyle=\frac{(t-t_{k-1})(t-t_{k})}{(t_{k-2}-t_{k-1})(t_{k-2}-t_{k})}f(t_{k-2})+\frac{(t-t_{k-2})(t-t_{k})}{(t_{k-1}-t_{k-2})(t_{k-1}-t_{k})}f(t_{k-1})$		(2.3)
		$\displaystyle\hskip 14.22636pt+\frac{(t-t_{k-2})(t-t_{k-1})}{(t_{k}-t_{k-2})(t_{k}-t_{k-1})}f(t_{k}),$		(2.3)

and the cubic interpolation of a function $f(t)$ on some interval $[t_{k-1},t_{k}],~{}k\geq 3,$ is constructed in the following form

I_{3k}f(t)=\sum_{l=0}^{3}{f(t_{k-l})}\prod_{i=0,i\neq l}^{3}\frac{t-t_{k-i}}{t_{k-l}-t_{k-i}}.

(2.4)

Thus, we have the interpolation error from [21] as

f(t)-I_{3k}f(t)=\frac{f^{(iv)}(\xi_{k})}{12}(t-t_{k-3})(t-t_{k-2})(t-t_{k-1})(t-t_{k}),

(2.5)

where $\xi_{k}\in(t_{k-3},t_{k}),~{}k\geq 3$ .

Definition 4.

The Supremum norm of a function $f(x)$ , $x\in X$ , is defined as

\|f\|_{\infty}=\sup_{x\in X}|f(x)|.

Then for any mesh function ${u_{i}^{n}}=u(x_{i},t_{n})$ , $(x_{i},t_{n})\in\bar{\Omega}_{h}\times\bar{\Omega}_{\tau}$ ,

\|u^{n}\|_{\infty}=\max_{0\leq i\leq K_{x}}|u_{i}^{n}|\ \text{ and }\ \|u\|_{\infty}=\max_{0\leq n\leq K_{t}}\|u^{n}\|_{\infty}.

3 Derivation and analysis of the higher order approximation to Caputo derivative on graded mesh

To approximate the solution of TFDE (1.1)-(1.3), the Caputo derivative is discretized on graded mesh, and the diffusion term is discretized on a uniform spatial grid. For graded temporal discretization, we discretize the time domain $[0,T]$ into $K_{t}$ subintervals with $t_{0}=0<t_{1}<\cdots<t_{K_{t}}=T$ . From here onwards, we denote $\mathbb{Z}_{[0,M]}:=\{z\in\mathbb{Z}$ , $0\leq z\leq M\}$ . We define $\tau_{n}={t_{n}}-t_{n-1}$ as the grid size in time direction with $t_{n}=T{\left(\frac{n}{K_{t}}\right)}^{r}$ , where $n\in\mathbb{Z}_{[0,K_{t}]}$ and $r\geq 1$ is the mesh grading parameter. Meanwhile, the spatial interval $[0,L]$ is divided into $K_{x}$ equal subintervals of uniform length $h=\frac{L}{K_{x}}$ , where $K_{x}\in\mathbb{Z^{+}}$ . Thus, the temporal and spatial meshes are given by

\bar{\omega}_{\tau}=\{t_{n}:t_{n}=T{\Big{(}\frac{n}{K_{t}}\Big{)}}^{r},n\in\mathbb{Z}_{[0,K_{t}]}\},~{}\text{and}~{}\bar{\omega}_{h}=\{x_{i}:x_{i}=ih,i\in\mathbb{Z}_{[0,K_{x}]}\},

respectively. Then the domain $\bar{\Omega}=[0,L]\times[0,T]$ is covered by $\bar{\omega}_{h\tau}:=\bar{\omega}_{h}\times\bar{\omega}_{\tau}$ . We define the following grid functions as

u(x_{i},t_{n})=U_{i}^{n}\text{ and }f(x_{i},t_{n})=f_{i}^{n}.

Moreover, we use the second-order central finite difference scheme to approximate the diffusion term in problem (1.1) as

	$\displaystyle{\partial_{x}^{2}}u(x_{i},t_{n})$	$\displaystyle=\frac{U_{i+1}^{n}-2U_{i}^{n}+U_{i-1}^{n}}{h^{2}}+\mathcal{O}(h^{2}),$		(3.1)
		$\displaystyle=\delta_{x}^{2}U_{i}^{n}+\mathcal{O}(h^{2}).$		(3.1)

From (1), the Caputo time-fractional derivative of a function $u(x,t)$ at mesh point $(x_{i},t_{n})$ and order $\alpha\in(0,1)$ is denoted and defined as

	$\displaystyle\textsubscript{C}D_{0,t}^{\alpha}u(x_{i},t_{n})$	$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t_{n}}(t_{n}-s)^{-\alpha}{\frac{\partial u(x_{i},s)}{\partial{s}}}ds,$
		$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\sum_{k=1}^{n}\int_{t_{k-1}}^{t_{k}}(t_{n}-s)^{-\alpha}{\frac{\partial u(x_{i},s)}{\partial{s}}}ds.$		(3.2)

3.1 Approximation to Caputo fractional derivative

For the foremost subinterval $[t_{0},t_{1}]$ , we apply linear interpolation (2.2) to approximate (3.2), thus we get

$\displaystyle\frac{1}{\Gamma(1-\alpha)}{\int_{t_{0}}^{t_{1}}}({{t_{n}}-s})^{-\alpha}{\frac{\partial{u}(x_{i},s)}{\partial{s}}}ds$	$\displaystyle\approx\frac{1}{\Gamma(1-\alpha)}{\int_{t_{0}}^{t_{1}}}({{t_{n}}-s})^{-\alpha}{\frac{\partial{I_{11}u}(x_{i},s)}{\partial{s}}}ds,$
	$\displaystyle=\frac{U_{i}^{1}-U_{i}^{0}}{\tau_{1}\Gamma(1-\alpha)}{\int_{t_{0}}^{t_{1}}}{({t_{n}}-s)}^{-\alpha}ds,$
	$\displaystyle=\frac{a_{n}}{\Gamma(2-\alpha)}\big{[}U_{i}^{1}-U_{i}^{0}\big{]},$	(3.3)

where $a_{n}={\frac{1}{\tau_{1}}}\left[({t_{n}}-{t_{0}})^{1-\alpha}-({t_{n}}-{t_{1}})^{1-\alpha}\right].$

Now for the second subinterval $[t_{1},t_{2}]$ , using quadratic interpolation $\eqref{quad_int}$ , we approximate (3.2) as follows

	$\displaystyle\frac{1}{\Gamma(1-\alpha)}{\int_{t_{1}}^{t_{2}}}({{t_{n}}-s})^{-\alpha}{\frac{\partial u(x_{i},s)}{\partial{s}}}ds$	$\displaystyle\approx\frac{1}{\Gamma(1-\alpha)}{\int_{t_{1}}^{t_{2}}}({{t_{n}}-s})^{-\alpha}{\frac{\partial{I_{22}u}(x_{i},s)}{\partial{s}}}ds,$		(3.4)
		$\displaystyle=\frac{1}{\Gamma(2-\alpha)}\big{[}b_{n}{U_{i}^{0}}+c_{n}{U_{i}^{1}}+d_{n}{U_{i}^{2}}\big{]},$		(3.4)

where

	$\displaystyle b_{n}=$	$\displaystyle{\frac{1}{\tau_{1}({\tau_{1}}+{\tau_{2}})}}\bigg{[}{\frac{2}{2-\alpha}}\big{\{}({t_{n}}-{t_{1}})^{2-\alpha}-({t_{n}}-{t_{2}})^{2-\alpha}\big{\}}-\tau_{2}\big{\{}({t_{n}}-{t_{1}})^{1-\alpha}+({t_{n}}-{t_{2}})^{1-\alpha}\big{\}}\bigg{]},$
	$\displaystyle c_{n}=$	$\displaystyle-{\frac{1}{{\tau_{1}}{\tau_{2}}}}\bigg{[}{\frac{2}{2-\alpha}}\big{\{}({t_{n}}-{t_{1}})^{2-\alpha}-({t_{n}}-{t_{2}})^{2-\alpha}\big{\}}+({\tau_{1}}-{\tau_{2}})({t_{n}}-{t_{1}})^{1-\alpha}-({\tau_{1}}+{\tau_{2}})({t_{n}}-{t_{2}})^{1-\alpha}\bigg{]},$
	$\displaystyle d_{n}=$	$\displaystyle{\frac{1}{\tau_{2}({\tau_{1}}+{\tau_{2}})}}\bigg{[}{\frac{2}{2-\alpha}}\big{\{}({t_{n}}-{t_{1}})^{2-\alpha}-({t_{n}}-{t_{2}})^{2-\alpha}\big{\}}+\tau_{1}({t_{n}}-{t_{1}})^{1-\alpha}-(2\tau_{2}+\tau_{1})({t_{n}}-{t_{2}})^{1-\alpha}\bigg{]}.$

Finally for other subintervals $[t_{j-1},t_{j}]$ , $3\leq j\leq n$ , using a cubic interpolation function (2.4) to approximate (3.2) provides

$\displaystyle{\frac{1}{\Gamma(1-\alpha)}}\sum_{j=3}^{n}{\int_{t_{j-1}}^{t_{j}}}({{t_{n}}-s})^{-\alpha}{\frac{\partial u(x_{i},s)}{\partial{s}}}ds$	$\displaystyle\approx{\frac{1}{\Gamma(1-\alpha)}}\sum_{j=3}^{n}{\int_{t_{j-1}}^{t_{j}}}({{t_{n}}-s})^{-\alpha}{\frac{\partial{I_{3j}u}(x_{i},s)}{\partial{s}}}ds,$
	$\displaystyle=\frac{1}{\Gamma(2-\alpha)}\sum_{k=3}^{n}\Big{[}w_{1,k}^{n}U_{i}^{k}+w_{2,k}^{n}U_{i}^{k-1}+w_{3,k}^{n}U_{i}^{k-2}$
	$\displaystyle\hskip 73.97733pt+w_{4,k}^{n}U_{i}^{k-3}\Big{]},$	(3.5)

where

	$\displaystyle w_{1,k}^{n}$	$\displaystyle=A_{1}\Big{[}\tau_{k-1}(\tau_{k-1}+\tau_{k-2}){(t_{n}-t_{k-1})}^{(1-\alpha)}-\big{\{}(\tau_{k}+\tau_{k-1})(2\tau_{k}+\tau_{k-1}+\tau_{k-2})$
		$\displaystyle\hskip 11.38092pt+\tau_{k}(\tau_{k}+\tau_{k-1}+\tau_{k-2})\big{\}}{(t_{n}-t_{k})}^{1-\alpha}+\frac{2}{(2-\alpha)}\big{\{}{(2\tau_{k-1}+\tau_{k-2})}{(t_{n}-t_{k-1})^{2-\alpha}}$
		$\displaystyle\hskip 11.38092pt-(3\tau_{k}+2\tau_{k-1}+\tau_{k-2})(t_{n}-t_{k})^{2-\alpha}\big{\}}+\frac{6}{(2-\alpha)(3-\alpha)}\big{\{}{(t_{n}-t_{k-1})}^{3-\alpha}-{(t_{n}-t_{k})}^{3-\alpha}\big{\}}\Big{]},$

	$\displaystyle w_{2,k}^{n}$	$\displaystyle=A_{2}\Big{[}\big{\{}\tau_{k-1}(\tau_{k}-\tau_{k-1}-\tau_{k-2})+\tau_{k}(\tau_{k-1}+\tau_{k-2})\big{\}}{(t_{n}-t_{k-1})}^{(1-\alpha)}$
		$\displaystyle\hskip 11.38092pt+(\tau_{k}+\tau_{k-1}+\tau_{k-2})(\tau_{k}+\tau_{k-1}){(t_{n}-t_{k})}^{1-\alpha}+\frac{2}{(2-\alpha)}\big{\{}(\tau_{k}-2\tau_{k-1}-\tau_{k-2})(t_{n}-t_{k-1})^{2-\alpha}$
		$\displaystyle\hskip 11.38092pt+{(2\tau_{k}+2\tau_{k-1}+\tau_{k-2})}{(t_{n}-t_{k})^{2-\alpha}}\big{\}}-\frac{6}{(2-\alpha)(3-\alpha)}\big{\{}{(t_{n}-t_{k-1})}^{3-\alpha}-{(t_{n}-t_{k})}^{3-\alpha}\big{\}}\Big{]},$
	$\displaystyle w_{3,k}^{n}$	$\displaystyle=A_{3}\Big{[}-\tau_{k}(\tau_{k-1}+\tau_{k-2}){(t_{n}-t_{k-1})}^{1-\alpha}-\tau_{k}(\tau_{k}+\tau_{k-1}+\tau_{k-2}){(t_{n}-t_{k})}^{(1-\alpha)}$
		$\displaystyle\hskip 11.38092pt+\frac{2}{(2-\alpha)}\big{\{}(\tau_{k-1}+\tau_{k-2}-\tau_{k})(t_{n}-t_{k-1})^{2-\alpha}-(2\tau_{k}+\tau_{k-1}+\tau_{k-2}){(t_{n}-t_{k})^{2-\alpha}}\big{\}}$
		$\displaystyle\hskip 11.38092pt+\frac{6}{(2-\alpha)(3-\alpha)}\big{\{}{(t_{n}-t_{k-1})}^{3-\alpha}-{(t_{n}-t_{k})}^{3-\alpha}\big{\}}\Big{]},$
	$\displaystyle w_{4,k}^{n}$	$\displaystyle=A_{4}\Big{[}\tau_{k-1}\tau_{k}{(t_{n}-t_{k-1})}^{(1-\alpha)}+\tau_{k}(\tau_{k}+\tau_{k-1}){(t_{n}-t_{k})}^{1-\alpha}$
		$\displaystyle\hskip 11.38092pt+\frac{2}{(2-\alpha)}\big{\{}{(\tau_{k}-\tau_{k-1})}{(t_{n}-t_{k-1})^{2-\alpha}}+(2\tau_{k}+\tau_{k-1})(t_{n}-t_{k})^{2-\alpha}\big{\}}$
		$\displaystyle\hskip 11.38092pt-\frac{6}{(2-\alpha)(3-\alpha)}\big{\{}{(t_{n}-t_{k-1})}^{3-\alpha}-{(t_{n}-t_{k})}^{3-\alpha}\big{\}}\Big{]},$

and

\displaystyle\begin{cases}\vspace{.2cm}A_{1}&=~{}\frac{1}{\tau_{k}({\tau_{k}}+{\tau_{k-1}})({\tau_{k}}+{\tau_{k-1}}+{\tau_{k-2}})},\\ \vspace{.2cm}A_{2}&=~{}\frac{1}{\tau_{k}\tau_{k-1}(\tau_{k-1}+\tau_{k-2})},\\ \vspace{.2cm}A_{3}&=~{}\frac{1}{\tau_{k-1}\tau_{k-2}({\tau_{k}}+{\tau_{k-1}})},\\ A_{4}&=~{}\frac{1}{\tau_{k-2}({\tau_{k-1}}+{\tau_{k-2}})({\tau_{k}}+{\tau_{k-1}}+{\tau_{k-2}})}.\\ \end{cases}

Now adjoining (3.3), (3.4), and (3.5), we obtain an approximation of the Caputo derivative of order $\alpha\in(0,1)$ as below

$\displaystyle\textsubscript{C}D_{0,t}^{\alpha}u(x_{i},t_{n})$	$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t_{n}}(t_{n}-s)^{-\alpha}{\frac{\partial u(x_{i},s)}{\partial{s}}}ds,$
	$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\bigg{[}{\int_{t_{0}}^{t_{1}}}({{t_{n}}-s})^{-\alpha}{\frac{\partial u(x_{i},s)}{\partial{s}}}ds+{\int_{t_{1}}^{t_{2}}}({{t_{n}}-s})^{-\alpha}{\frac{\partial u(x_{i},s)}{\partial{s}}}ds$
	$\displaystyle\hskip 56.9055pt+\sum_{j=3}^{n}{\int_{t_{j-1}}^{t_{j}}}({{t_{n}}-s})^{-\alpha}{\frac{\partial u(x_{i},s)}{\partial{s}}}ds\bigg{]},$
	$\displaystyle\approx\frac{1}{\Gamma(2-\alpha)}\sum_{j=0}^{n}{p_{j}}{U_{i}^{n-j}},$	(3.6)
	$\displaystyle=\textsubscript{C}D_{0,t_{K_{t}}}^{\alpha}u(x_{i},t_{n}).$

Here the coefficients $p_{j}$ have different values for different $n$ ,

		$\displaystyle\text{for }n=1,\ p_{0}=a_{1},\ \ \ p_{1}=-a_{1};$
		$\displaystyle\text{for }n=2,\ p_{0}=d_{2},\ \ \ p_{1}=a_{2}+c_{2},\ \ \ \ \ \ p_{2}=b_{2}-a_{2};$
		$\displaystyle\text{for }n=3,\ p_{0}=w_{1,3}^{3},\ p_{1}=w_{2,3}^{3}+d_{3},\ \ \ p_{2}=w_{3,3}^{3}+a_{3}+c_{3},\ \ \ \hskip 8.5359ptp_{3}=w_{4,3}^{3}+b_{3}-a_{3};$
		$\displaystyle\text{for }n=4,\ p_{0}=w_{1,4}^{4},\ p_{1}=w_{1,3}^{4}+w_{2,4}^{4},\ \hskip 0.28436ptp_{2}=w_{3,4}^{4}+w_{2,3}^{4}+d_{4},\ \ \ p_{3}=w_{3,3}^{4}+w_{4,4}^{4}+a_{4}+c_{4},$
		$\displaystyle\hskip 48.36958ptp_{4}=w_{4,3}^{4}+b_{4}-a_{4};$
		$\displaystyle\text{for }n=5,\ p_{0}=w_{1,5}^{5},\ p_{1}=w_{1,4}^{5}+w_{2,5}^{5},\ p_{2}=w_{1,3}^{5}+w_{2,4}^{5}+w_{3,5}^{5},\ p_{3}=w_{2,3}^{5}+w_{3,4}^{5}+w_{4,5}^{5}+d_{5},$
		$\displaystyle\hskip 48.36958ptp_{4}=w_{3,3}^{5}+w_{4,4}^{5}+a_{5}+c_{5},\ \ \ p_{5}=w_{4,3}^{5}+b_{5}-a_{5};$

for $6\leq n\leq K_{t}$ , the relation between the coefficients takes form as follows:

		$\displaystyle\hskip 8.5359ptp_{0}=w_{1,n}^{n},\ \ p_{1}=w_{1,n-1}^{n}+w_{2,n}^{n},\ \ p_{2}=w_{1,n-2}^{n}+w_{2,n-1}^{n}+w_{3,n}^{n},\hskip 56.9055pt$		(3.7)
		$\displaystyle\hskip 8.5359ptp_{k}=w_{1,n-k}^{n}+w_{2,n-k+1}^{n}+w_{3,n-k+2}^{n}+w_{4,n-k+3}^{n},\ (3\leq k\leq n-3),$
		$\displaystyle\hskip 8.5359ptp_{n-2}=w_{2,3}^{n}+w_{3,4}^{n}+w_{4,5}^{n}+d_{n},$
		$\displaystyle\hskip 8.5359ptp_{n-1}=w_{3,3}^{n}+w_{4,4}^{n}+a_{n}+c_{n},\ \ \ \hskip 5.69046ptp_{n}=w_{4,3}^{n}+b_{n}-a_{n}.$

3.2 Truncation error estimate

This portion covers the estimation of truncation error for the approximation to the Caputo fractional derivative.

Lemma 1.

Let $u(\cdot,t)\in C[0,T]\cap C^{4}(0,T]$ . Suppose that $\left|\frac{\partial^{l}u(\cdot,t)}{\partial t^{l}}\right|\leq C(1+t^{\alpha-l})$ for $l=0,1,\ldots,4,$ there exists a constant C such that $\forall$ $(x_{i},t_{n})\in\Omega_{h}\times\bar{\Omega}_{\tau},$ we have

\displaystyle\left|\textsubscript{C}D_{0,t_{K_{t}}}^{\alpha}u(x_{i},t_{n})-\textsubscript{C}D_{0,t}^{\alpha}u(x_{i},t_{n})\right|\leq Cn^{-\min\{4-\alpha,r\alpha\}},

(3.8)

where C is a generic positive constant independent of $\tau_{n}$ .

Proof.

According to Taylor’s expansion formula,

	$\displaystyle\tau_{k+1}=t_{k+1}-t_{k}$	$\displaystyle=T\Big{(}\frac{k+1}{K_{t}}\Big{)}^{r}-T\Big{(}\frac{k}{K_{t}}\Big{)}^{r},$
		$\displaystyle\leq CT{K_{t}}^{-r}k^{r-1},\ k=0,1,\ldots,K_{t}-1.$		(3.9)

For $n=1,2,\ldots,K_{t},\text{ and }i=1,2,\ldots,K_{x}-1,$ we have

		$\displaystyle\textsubscript{C}D_{0,t_{K_{t}}}^{\alpha}u_{i}^{n}-\textsubscript{C}D_{0,t}^{\alpha}u(x_{i},t_{n})=\sum_{k=1}^{n}T_{nk},$
		$\displaystyle\left\|\textsubscript{C}D_{0,t_{K_{t}}}^{\alpha}u_{i}^{n}-\textsubscript{C}D_{0,t}^{\alpha}u(x_{i},t_{n})\right\|\leq\sum_{k=1}^{n}\left\|T_{nk}\right\|.$		(3.10)

According to our algorithm, we need to form at least three time subintervals to discretize the Caputo derivative, i.e., $n\geq 3$ . Here, for $k=1$ in (3.10), we get

\displaystyle T_{n1}=\frac{1}{\Gamma(1-\alpha)}{\int_{t_{0}}^{t_{1}}}({{t_{n}}-s})^{-\alpha}\bigg{[}{\frac{\partial I_{11}u}{\partial s}(x_{i},s)-\frac{\partial u}{\partial s}(x_{i},s)}\bigg{]}ds.

Moreover, from [27], we have

\displaystyle\left|T_{n1}\right|\leq Cn^{-r\alpha},\text{ where }n=3,4,\ldots,K_{t}.

(3.11)

Now for $k=2$ ,

\displaystyle T_{n2}=\frac{1}{\Gamma(1-\alpha)}{\int_{t_{1}}^{t_{2}}}({{t_{n}}-s})^{-\alpha}\bigg{[}{\frac{\partial I_{22}u}{\partial s}(x_{i},s)-\frac{\partial u}{\partial s}(x_{i},s)}\bigg{]}ds.

From [21], we have

\displaystyle\left|T_{n2}\right|\leq Cn^{-r(1+\alpha)},\text{ where }n=3,4,\ldots,K_{t}.

(3.12)

For $k=3,4,\ldots,K_{t},$

\displaystyle T_{nk}=\frac{1}{\Gamma(1-\alpha)}{\int_{t_{k-1}}^{t_{k}}}({{t_{n}}-s})^{-\alpha}\bigg{[}{\frac{\partial I_{3k}u}{\partial s}(x_{i},s)-\frac{\partial u}{\partial s}(x_{i},s)}\bigg{]}ds.

For $n=3,4,\ldots,K_{t},\text{ and }k=3,4,\ldots,\lceil{\frac{n-1}{2}}\rceil$ , using (2.5) here we get,

	$\displaystyle T_{nk}$	$\displaystyle=\frac{1}{\Gamma(1-\alpha)}{\int_{t_{k-1}}^{t_{k}}}({{t_{n}}-s})^{-\alpha}\bigg{[}{\frac{\partial I_{3k}u}{\partial s}(x_{i},s)-\frac{\partial u}{\partial s}(x_{i},s)}\bigg{]}ds,$
		$\displaystyle=\frac{-\alpha}{\Gamma(1-\alpha)}{\int_{t_{k-1}}^{t_{k}}}({{t_{n}}-s})^{-\alpha-1}\frac{u_{tttt}(x_{i},\phi_{k})}{12}(s-t_{k-3})(s-t_{k-2})(s-t_{k-1})(s-t_{k})ds.$

Then, using the assumption $\left|\frac{\partial^{l}u(\cdot,t)}{\partial t^{l}}\right|\leq C(1+t^{\alpha-l})$ for $l=0,1,\ldots,4$ , we get

\left|T_{nk}\right|\leq\ Ct_{k-1}^{\alpha-4}\tau_{k}^{4}\int_{t_{k-1}}^{t_{k}}(t_{n}-s)^{-\alpha-1}ds\ \leq Ct_{k-1}^{\alpha-4}\tau_{k}^{5}(t_{n}-t_{k})^{-\alpha-1},

(3.13)

where

$\displaystyle(t_{n}-t_{k})^{-\alpha-1}={\bigg{(}T\frac{n^{r}-k^{r}}{{K_{t}}^{r}}\bigg{)}}^{-\alpha-1}$	$\displaystyle=T^{-\alpha-1}\bigg{(}{\frac{{K_{t}}^{r}}{n^{r}-k^{r}}}\bigg{)}^{\alpha+1},$
	$\displaystyle\leq T^{-\alpha-1}\bigg{(}{\frac{{K_{t}}^{r}}{n^{r}-{\lceil\frac{n-1}{2}\rceil}^{r}}}\bigg{)}^{\alpha+1},$
	$\displaystyle\leq CT^{-\alpha-1}{\bigg{(}\frac{K_{t}}{n}\bigg{)}}^{r(1+\alpha)}.$	(3.14)

Now, using (3.9) and (3.14) in (3.13), we get

	$\displaystyle\left\|T_{nk}\right\|$	$\displaystyle\leq C\bigg{[}T\bigg{(}\frac{k-1}{K_{t}}\bigg{)}^{r}\bigg{]}^{\alpha-4}T^{5}{K_{t}}^{-5r}(k-1)^{5(r-1)}T^{-\alpha-1}\left(\frac{K_{t}}{n}\right)^{r(1+\alpha)},$
		$\displaystyle\leq C{K_{t}}^{r(4-\alpha)-5r+r(1+\alpha)}n^{-r(1+\alpha)}(k-1)^{r(\alpha-4)+5(r-1)},$
		$\displaystyle\leq Cn^{-r(1+\alpha)}(k-1)^{r(1+\alpha)-5}.$

Hence,

		$\displaystyle\sum_{k=3}^{\lceil\frac{n-1}{2}\rceil}T_{nk}\leq Cn^{-r(1+\alpha)}\sum_{k=3}^{\lceil\frac{n-1}{2}\rceil}(k-1)^{r(1+\alpha)-5},$
		$\displaystyle\hskip 44.10185pt\leq\begin{cases}Cn^{-r(1+\alpha)},\ \ \ r(1+\alpha)<4,\\ Cn^{-4}\ln(n),\ \ r(1+\alpha)=4,\\ Cn^{-4},\ \ \ \ \ \ \ \ \ r(1+\alpha)>4,\end{cases}$		(3.15)

and for $k=\lceil{\frac{n-1}{2}}\rceil+1,\ldots,n$ , we have

$\displaystyle\sum_{k=\lceil\frac{n-1}{2}\rceil+1}^{n-1}T_{nk}$	$\displaystyle=\sum_{k=\lceil\frac{n-1}{2}\rceil+1}^{n-1}t_{k-1}^{\alpha-4}\tau_{k}^{4}\int_{t_{k-1}}^{t_{k}}(t_{n}-s)^{-\alpha-1}ds,$
	$\displaystyle\leq Ct_{k-1}^{\alpha-4}\tau_{k}^{4}\sum_{k=\lceil\frac{n-1}{2}\rceil+1}^{n-1}\frac{1}{\alpha}\bigg{[}(t_{n}-t_{k})^{-\alpha}-(t_{n}-t_{k-1})^{-\alpha}\bigg{]},$
	$\displaystyle\leq Ct_{k-1}^{\alpha-4}\tau_{k}^{4}\bigg{[}(t_{n}-t_{n-1})^{-\alpha}-(t_{n}-t_{\lceil\frac{n-1}{2}\rceil})^{-\alpha}\bigg{]},$
	$\displaystyle\leq Ct_{k-1}^{\alpha-4}\tau_{k}^{4}\tau_{n}^{-\alpha}\leq Ct_{k-1}^{\alpha-4}\bigg{[}TK_{t}^{-r}(k-1)^{r-1}\bigg{]}^{4}\tau_{n}^{-\alpha},$
	$\displaystyle\leq Ct_{n}^{\alpha-4}T^{4}K_{t}^{-4r}n^{4(r-1)}\bigg{[}TK_{t}^{-r}n^{r-1}\bigg{]}^{-\alpha},$
	$\displaystyle\leq CT^{4}K_{t}^{-4r}n^{4(r-1)}T^{\alpha-4}\bigg{(}\frac{n}{K_{t}}\bigg{)}^{r(\alpha-4)}T^{-\alpha}K_{t}^{r\alpha}n^{-\alpha(r-1)},$
	$\displaystyle\leq Cn^{-(4-\alpha)}.$	(3.16)

Finally, we consider $T_{nn}$ for $n>2$ as

$\displaystyle T_{nn}$	$\displaystyle=\frac{1}{\Gamma(1-\alpha)}{\int_{t_{n-1}}^{t_{n}}}({{t_{n}}-s})^{-\alpha}\bigg{[}{\frac{\partial I_{3k}u}{\partial s}(x_{i},s)-\frac{\partial u}{\partial s}(x_{i},s)}\bigg{]}ds,$
	$\displaystyle\leq Ct_{n-1}^{\alpha-4}\tau_{n-1}^{3}\int_{t_{n-1}}^{t_{n}}({{t_{n}}-s})^{-\alpha}ds,$
	$\displaystyle\leq Ct_{n-1}^{\alpha-4}\tau_{n-1}^{3}({t_{n}}-t_{n-1})^{1-\alpha},$
	$\displaystyle\leq C\big{[}TK_{t}^{-r}(n-1)^{r}\big{]}^{\alpha-4}\big{[}TK_{t}^{-r}(n-2)^{r-1}\big{]}^{3}\big{[}TK_{t}^{-r}(n-1)^{r-1}\big{]}^{1-\alpha},$
	$\displaystyle\leq Cn^{-(4-\alpha)}.$	(3.17)

Combining formulas (3.11), (3.12), and (3.15)-(3.17), we obtain

\displaystyle\left|\textsubscript{C}D_{0,t_{K_{t}}}^{\alpha}u(x_{i},t_{n})-\textsubscript{C}D_{0,t}^{\alpha}u(x_{i},t_{n})\right|\leq Cn^{-\min\{4-\alpha,r\alpha\}}.

∎

Remark 1.

We get the truncation error of the approximation to Caputo derivative of a non-smooth function on the graded mesh as $\mathcal{O}(\tau^{\min\{4-\alpha,r\alpha\}})$ . However, in [2], the truncation error of Caputo derivative approximation for a smooth function on the uniform mesh is obtained as $\mathcal{O}(\tau^{4-\alpha})$ . Hence it can be easily observed that for a non-smooth function, grading improves the convergence of the approximation and for $r\geq\frac{4-\alpha}{\alpha}$ , the convergence will always be optimum and equal to that of the smooth solution, i.e., $(4-\alpha)$ .

In this next section, with the help of the above approximation of the Caputo derivative, we construct a high-order numerical scheme for TFDE (1.1)-(1.3) and perform some important analysis on it.

4 Numerical scheme for TFDE

For the function $u(x,t)$ , its analytical and approximate solutions at the mesh point $(x_{i},t_{n})$ are denoted by $U_{i}^{n}$ and $u_{i}^{n}$ , respectively. In view of the obtained results (3.1) and (3.8), we consider the following approximation of the TFDE (1.1):

\displaystyle L_{K_{x},K_{t}}u_{i}^{n}=\textsubscript{C}D_{0,t_{K_{t}}}^{\alpha}u_{i}^{n}-\varrho\delta_{x}^{2}u_{i}^{n}=f_{i}^{n},\ 1\leq i\leq K_{x}-1,1\leq n\leq K_{t},

(4.1)

Therefore, we derive the following finite difference scheme from (4.1):

		$\displaystyle\frac{1}{\Gamma(2-\alpha)}\Big{[}{p_{0}}u_{i}^{n}+{p_{1}}u_{i}^{n-1}+{p_{2}}u_{i}^{n-2}+\sum_{k=3}^{n-3}{p_{k}}u_{i}^{n-k}+{p_{n-2}}u_{i}^{2}+{p_{n-1}}u_{i}^{1}+{p_{n}}u_{i}^{0}\Big{]}$		(4.2)
		$\displaystyle\hskip 42.67912pt=\varrho\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{h^{2}}+f_{i}^{n},$		(4.2)

where $1\leq n\leq K_{t}$ , and $1\leq i\leq K_{x}-1$ . The initial and boundary conditions (1.2)-(1.3) can be discretized as

	$\displaystyle u_{0}^{n}$	$\displaystyle=0=u_{K_{x}}^{n},\ 0<n\leq K_{t},$		(4.3)
	$\displaystyle u_{i}^{0}$	$\displaystyle=\varphi(x_{i}),\ 0\leq i\leq K_{x}.$		(4.4)

Multiplying (4.2) by $h^{2}$ and introducing the parameter $\mu=\frac{h^{2}}{\Gamma(2-\alpha)}$ to simplify the expression, we get

$\displaystyle\text{ for }n=1,\ -\varrho u_{i-1}^{1}+(\mu p_{0}+2\varrho)u_{i}^{1}-\varrho u_{i+1}^{1}$	$\displaystyle=-\mu{p_{1}}u_{i}^{0}+h^{2}{f_{i}^{1}};$	(4.5)
$\displaystyle\text{ for }n=2,\ -\varrho u_{i-1}^{2}+(\mu p_{0}+2\varrho)u_{i}^{2}-\varrho u_{i+1}^{2}$	$\displaystyle=-\mu{p_{1}}u_{i}^{1}-\mu{p_{2}}u_{i}^{0}+h^{2}{f_{i}^{2}};$
$\displaystyle\text{ for }n\geq 3,\ -\varrho u_{i-1}^{n}+(\mu p_{0}+2\varrho)u_{i}^{n}-\varrho u_{i+1}^{n}$	$\displaystyle=-\mu\sum_{k=1}^{n}{p_{k}}u_{i}^{n-k}+h^{2}{f_{i}^{n}}.$

The system of linear equations (4.5) can also be rewritten in the following matrix form:

	$\displaystyle\text{ for }n=1,\ {A}u^{1}$	$\displaystyle=-\mu p_{1}{u^{0}}+h^{2}{F^{1}};$
	$\displaystyle\text{ for }n=2,\ {A}u^{2}$	$\displaystyle=-\mu p_{1}{u^{1}}-\mu p_{2}{u^{0}}+h^{2}{F^{2}};$
	$\displaystyle\text{ for }n\geq 3,\ {A}u^{n}$	$\displaystyle=-\mu{\sum_{k=1}^{n}}p_{k}{u^{n-k}}+h^{2}{F^{n}};$

where the matrix and vectors are defined by

$\displaystyle A$	$\displaystyle=\ tri\big{[}-\varrho,\mu{p_{0}}+2\varrho,-\varrho\big{]},$	(4.6)
$\displaystyle u^{n}$	$\displaystyle=\ (u_{1}^{n},u_{2}^{n},\ldots,u_{K_{x}-1}^{n})^{T},$	(4.7)
$\displaystyle F^{n}$	$\displaystyle=\ (f_{1}^{n},f_{2}^{n},\ldots,f_{K_{x}-1}^{n})^{T},\text{ for }1\leq n\leq K_{t}.$	(4.8)

Theorem 1.

The solution of the difference scheme (4.2)-(4.4) is unique.

Proof.

The coefficient matrix $A$ in (4.6) is strictly diagonally dominant irrespective of any possible values of $h,\tau_{n}$ , and $\alpha$ , so it is non-singular. Hence, $A$ is invertible. Consequently, the difference scheme (4.2)-(4.4) has a unique solution. ∎

4.1 Stability and convergence analysis of the scheme on uniform mesh

In this segment, we discuss the stability and convergence of the proposed difference scheme (4.2)-(4.4) on a uniform mesh. To accomplish that, we first rewrite the difference scheme with $r=1$ in a desired manner.
For the case $r=1$ , the coefficients become $p_{j}=\tau^{-\alpha}g_{j}$ . Hence, the scheme (4.2) is now

	$\displaystyle\frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}\left[\sum_{k=0}^{n}{g_{k}}u_{i}^{n-k}\right]=\varrho\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{h^{2}}+f_{i}^{n},$
	$\displaystyle-\eta\epsilon u_{i+1}^{n}+(g_{0}+2\eta\epsilon)u_{i}^{n}-\eta\epsilon u_{i-1}^{n}=-\sum_{k=1}^{n}{g_{k}}u_{i}^{n-k}+\epsilon f_{i}^{n},$

where $\eta=\frac{\varrho}{h^{2}}$ and $\epsilon=\tau^{\alpha}\Gamma(2-\alpha)$ . So the matrix equation of the above scheme is

Tu^{n}=-{\sum_{k=1}^{n}}g_{k}{u^{n-k}}+\epsilon{F^{n}},

(4.9)

where the matrix $T$ is given by

T=\ tri\big{[}-\eta\epsilon,g_{0}+2\eta\epsilon,-\eta\epsilon\big{]},

and other vectors $u^{n}$ and $F^{n}$ are from (4.7) and (4.8), respectively.

Following is a result from [2], which will be used to establish the stability and convergence of the proposed scheme on uniform mesh, i.e., $r=1$ .

Lemma 2.

For the case $r=1$ , the following properties of the coefficients $g_{j}$ hold:

$(1)$

When $n\geq 3$ , for any order $\alpha\in(0,1)$ ,
1. (a)
  
  $g_{0}=\frac{1}{3}+\frac{1}{2-\alpha}+{\frac{1}{(2-\alpha)(3-\alpha)}}\in(1,\frac{11}{6})$ ,
2. (b)
  
  $g_{2}>0$ ,
3. (c)
  
  $g_{j}<0,\ j\geq 1,j\neq 2$ .
$(2)$

$\sum_{j=0}^{n}{g_{j}}=0$ , for any order $\alpha\in(0,1)$ and $\text{for all }n$ .

Proof.

For proof, we refer [2, Lemma 2.1 and 2.2] ∎

Next, we derive a lemma that will be further utilized to obtain the stability estimate of the scheme.

Lemma 3.

The norm $\|T^{-1}\|_{\infty}\leq 1$ where $\|T\|_{\infty}=\max_{1\leq i\leq K_{x}-1}\Big{\{}\sum_{j=1}^{K_{x}-1}|t_{i,j}|\Big{\}}$ .

Proof.

The above-given statement can be proved with the help of the Gerschgorin theorem. It states that all the eigenvalues of the matrix $T$ should lie in the union of $K_{x}-1$ circles centered at $t_{i,i}$ with radius $r_{i}=\sum^{K_{x}-1}_{k=1,k\neq i}|t_{i,k}|$ . We need to show that every eigenvalue of the matrix $T$ has a magnitude greater than 1. Hence, using the Gerschgorin theorem and the properties of matrix $T$ , we have

		$\displaystyle t_{i,i}=g_{0}+2\eta\epsilon,\ \ r_{i}=\sum^{K_{x}-1}_{k=1,k\neq i}\|t_{i,k}\|\leq 2\eta\epsilon,$
	and so	$\displaystyle\|\lambda-(g_{0}+2\eta\epsilon)\|\leq 2\eta\epsilon.$

Since for $r=1$ from Lemma 2, we have $g_{0}\in(1,11/6)$ , and therefore we conclude that the eigenvalues of the matrix $T$ have a magnitude greater than $1$ . Also, the matrix $T$ is diagonally dominant. Thus $T$ is invertible and the eigenvalues of $T^{-1}$ are less than or equal to $1$ in magnitude. Therefore,

\|T^{-1}\|_{\infty}\leq 1.

∎

Now the stability of the proposed difference scheme derived above is proved in the following theorem with the help of this lemma.

Theorem 2.

The proposed finite difference scheme defined by (4.2)-(4.4) to the TFDE (1.1)-(1.3) for $r=1$ is unconditionally stable.

Proof.

Let $u_{i}^{n}$ and $\tilde{u}_{i}^{n}$ be the approximate solutions of (4.2) corresponding to different initial conditions. Therefore, $v_{i}^{n}=u_{i}^{n}-\tilde{u}_{i}^{n}$ be another solution of (4.2) where $v_{i}^{0}\neq 0$ and $v^{n}=(v_{1}^{n},v_{2}^{n},\ldots,v_{K_{x}-1}^{n})^{T}$ . The difference solution $v^{n}$ will also satisfy (4.9) with $F^{n}=\textbf{0}$ .

Tv^{n}=-{\sum_{k=1}^{n}}g_{k}{v^{n-k}}.

From the technique of mathematical induction, we prove $\|v^{n}\|\leq C\|v^{0}\|,\ n=0,1,2,\ldots,K_{t}$ , where C is some positive constant. Hence for $n=1$ in the above equation, we have

	$\displaystyle Tv^{1}$	$\displaystyle=-g_{1}v^{0},$
	$\displaystyle\\|v^{1}\\|_{\infty}$	$\displaystyle\leq\left\\|T^{-1}\right\\|_{\infty}\ \left\\|\left(-g_{1}v^{0}\right)\right\\|_{\infty},$
		$\displaystyle\leq\|-g_{1}\|\ \left\\|v^{0}\right\\|_{\infty}.$

For $n=1$ , we have $g_{1}=-g_{0}$ and $g_{0}=1$ . So

	$\displaystyle\left\\|v^{1}\right\\|_{\infty}$	$\displaystyle\leq\|g_{0}\|\ \left\\|v^{0}\right\\|_{\infty},$
	$\displaystyle\left\\|v^{1}\right\\|_{\infty}$	$\displaystyle\leq\left\\|v^{0}\right\\|_{\infty}.$

Now assume that $\|v^{s}\|_{\infty}\leq C\|v^{0}\|_{\infty},\ \forall s\leq n$ . We will prove it is also true for $s=n+1$ . We have,

	$\displaystyle Tv^{n+1}$	$\displaystyle=-{\sum_{k=1}^{n+1}}g_{k}{v^{n+1-k}},$
	$\displaystyle\left\\|v^{n+1}\right\\|_{\infty}$	$\displaystyle=\left\\|T^{-1}\left(-\sum_{k=1}^{n+1}{g_{k}}v^{n+1-k}\right)\right\\|_{\infty},$
		$\displaystyle\leq\left\\|T^{-1}\right\\|_{\infty}\ \left\\|\left(-\sum_{k=1}^{n+1}{g_{k}}v^{n+1-k}\right)\right\\|_{\infty},$
		$\displaystyle\leq\sum_{k=1}^{n+1}\|{g_{k}}\|\ \left\\|v^{n+1-k}\right\\|_{\infty},$
		$\displaystyle\leq C\left\\|v^{0}\right\\|_{\infty}\ \left(\sum_{k=1}^{n+1}\|{g_{k}}\|\right),$

Using the identities $(1)$ and $(2)$ of Lemma 2, we get

	$\displaystyle\left\\|v^{n+1}\right\\|_{\infty}$	$\displaystyle\leq C(g_{0}+2g_{2})\left\\|v^{0}\right\\|_{\infty},$
		$\displaystyle\leq C_{1}\left\\|v^{0}\right\\|_{\infty}.$

Hence the scheme is unconditionally stable. ∎

Theorem 3.

The solutions $u(x_{i},t_{n})$ of the TFDE (1.1)-(1.3) and $u_{i}^{n}$ of the difference scheme (4.2)-(4.4) at grid point $(x_{i},t_{n})$ , for $r=1$ and sufficiently small $h$ and $\tau$ , satisfy

\|e^{n}\|_{\infty}=\mathcal{O}\big{(}{h^{2}}+\tau^{\alpha}\big{)},\ n=1,2,\ldots,K_{t},

(4.10)

where $e_{i}^{n}=u(x_{i},t_{n})-u_{i}^{n}$ , $i=1,2,\ldots,K_{x}-1$ .

Proof.

Let us fix $(x_{i},t_{n})\in\Omega$ . Using (4.1), (4.3)-(4.4), we obtain the following difference equation for the error

		$\displaystyle L_{K_{x},K_{t}}e_{i}^{n}=r_{i}^{n},~{}1\leq i\leq K_{x}-1,\ 1\leq n\leq K_{t},$
		$\displaystyle e_{0}^{n}=0=e_{K_{x}}^{n},~{}0<n\leq K_{t},$
		$\displaystyle e_{i}^{0}=0,~{}0\leq i\leq K_{x}-1.$

We have from our error equation, (3.1) and Lemma 1,

	$\displaystyle\|r_{i}^{n}\|=\|L_{K_{x},K_{t}}e_{i}^{n}\|$	$\displaystyle=\big{\|}\textsubscript{C}D_{0,K_{t}}^{\alpha}e_{i}^{n}-\delta_{x}^{2}e_{i}^{n}\big{\|},$
		$\displaystyle=\Big{\|}\textsubscript{C}D_{0,K_{t}}^{\alpha}\left[u(x_{i},t_{n})-u_{i}^{n}\right]-\delta_{x}^{2}\left[u(x_{i},t_{n})-u_{i}^{n}\right]\Big{\|},$
		$\displaystyle\hskip 2.84544pt\leq C\big{(}{h^{2}}+n^{-\min\{4-\alpha,r\alpha\}}\big{)},$

For $r=1$ , it will convert into

	$\displaystyle\|r_{i}^{n}\|$	$\displaystyle\leq C\big{(}{h^{2}}+n^{-\min\{4-\alpha,\alpha\}}\big{)},$
		$\displaystyle=C\big{(}{h^{2}}+n^{-\alpha}\big{)}.$

Now for $n=1$ , we have from error equation

\displaystyle-\eta\epsilon e_{i+1}^{1}+(g_{0}+2\eta\epsilon)e_{i}^{1}-\eta\epsilon e_{i-1}^{1}=r_{i}^{1},

And for $n>1$

-\eta\epsilon e_{i+1}^{n}+(g_{0}+2\eta\epsilon)e_{i}^{n}-\eta\epsilon e_{i-1}^{n}+\sum_{k=1}^{n}{g_{k}}e_{i}^{n-k}=r_{i}^{n},

We use the mathematical induction method to prove the convergence of the scheme. If $n=1$ , we choose $|e_{l}^{1}|=\max_{1\leq i\leq K_{x}-1}|e_{i}^{1}|$ and hence,

	$\displaystyle\|e_{l}^{1}\|$	$\displaystyle\leq(g_{0}+2\eta\epsilon)\|e_{l}^{1}\|-\eta\epsilon\|e_{l+1}^{1}\|-\eta\epsilon\|e_{l-1}^{1}\|,$
		$\displaystyle\leq\|(g_{0}+2\eta\epsilon)e_{l}^{1}-\eta\epsilon e_{l+1}^{1}-\eta\epsilon e_{l-1}^{1}\|,$
		$\displaystyle=\|r_{l}^{1}\|\leq C({h^{2}}+\tau^{\alpha}),$
	$\displaystyle\\|e^{1}\\|_{\infty}$	$\displaystyle\leq C({h^{2}}+\tau^{\alpha}).$

Suppose that if $s\leq n-1$ , $||e^{s}||_{\infty}\leq M({h^{2}}+\tau^{\alpha})$ holds, then for $s=n$ , let $|e_{l}^{n}|=\max_{1\leq i\leq K_{x}-1}|e_{i}^{n}|$ , we have

	$\displaystyle\|e_{l}^{n}\|$	$\displaystyle\leq(g_{0}+2\eta\epsilon)\|e_{l}^{n}\|-\eta\epsilon\|e_{l+1}^{n}\|-\eta\epsilon\|e_{l-1}^{n}\|,$
		$\displaystyle\leq\left\|(g_{0}+2\eta\epsilon)e_{l}^{n}-\eta\epsilon e_{l+1}^{n}-\eta\epsilon e_{l-1}^{n}\right\|,$
		$\displaystyle\leq\left\|-\sum_{k=1}^{n}{g_{k}}e_{l}^{n-k}+r_{l}^{n}\right\|,$
		$\displaystyle\leq\sum_{k=1}^{n}\|{g_{k}}\|\|e_{l}^{n-k}\|+\|r_{l}^{n}\|,$
		$\displaystyle\leq M_{1}({h^{2}}+\tau^{\alpha})\left(\sum_{k=1}^{n}\|{g_{k}}\|+1\right).$

Similar as above, using the identities $(1)$ and $(2)$ of Lemma 2, we get

\displaystyle|e_{l}^{n}|

\displaystyle\leq(1+g_{0}+2g_{2})M_{1}({h^{2}}+\tau^{\alpha}).

Thus $\|e^{n}\|_{\infty}\leq M_{0}\big{(}{h^{2}}+\tau^{\alpha}\big{)}$ . ∎

Remark 2.

As discussed in Remark 1, we know from [2, Theorem 4.1] that for the TFDE problem (1.1)-(1.3) having a smooth solution, i.e., $u(\cdot,t)\in C^{4}[0,T]$ on uniform mesh, the local truncation error for the approximation of Caputo fractional derivative is $\mathcal{O}(\tau^{4-\alpha})$ . Therefore, keeping in view the proof of the above theorem, one can deduce the rate of convergence for the proposed scheme (4.2)-(4.4) in case of smooth solution on a uniform mesh as $\mathcal{O}(\tau^{4-\alpha}+h^{2})$ . Also, the stability result in [2] is proved via the Fourier series method for $\alpha\in(0,\bar{\bar{\bar{\alpha}}}\approx 0.368)$ , whereas we have shown the stability through the matrix method for any $\alpha\in(0,1)$ .

Remark 3.

The above stability and convergence results for the difference scheme (4.2)-(4.4) have been proved for $r=1$ i.e., on a uniform mesh. The rationale behind it is that the relation between the coefficients $p_{j}$ for $r=1$ , i.e. $g_{j}$ , is known, but for $r>1$ , the behavior of the coefficients becomes very complex. From the table below we can see that the signum behavior of the coefficients fluctuates for different $r$ . For the same number of iterations, the number of positive coefficients varies for different values of $r$ which indicates some of the coefficients are changing their signs. This behavioral change also depends on $n$ as the change in the number of iterations for a fix $r$ changes the number of coefficients of a particular sign. Thus, the behavior of coefficients alters for different $n$ and $r$ for a fixed value of $\alpha$ (here $\alpha=0.5$ ) and hence we are unable to conclude about the stability of the scheme for $r>1$ .

	$r=4$		$r=6$		$r=10$
$n$	+ve sign	-ve sign	+ve sign	-ve sign	+ve sign	-ve sign
$10$	$4$	$7$	$5$	$6$	$4$	$7$
$20$	$4$	$17$	$5$	$16$	$6$	$15$
$30$	$3$	$28$	$6$	$25$	$9$	$22$
$40$	$4$	$37$	$7$	$34$	$12$	$29$
$50$	$6$	$45$	$9$	$42$	$16$	$35$

Table 1: Number of +ve and -ve coefficients

(p_{j})

for

N=50

and

\alpha=0.5

for different values of

n

5 Numerical experiments

In this section, we verify the theoretical analysis provided in the previous section by producing numerical results for two test examples of TFDE problem (1.1)-(1.3) with smooth and non-smooth solutions, respectively, using the proposed difference scheme (4.2)-(4.4). We examine the accuracy and efficiency of the scheme for different values of $K_{x}$ , $K_{t}$ , $\alpha$ , and $r$ . The exact solution to each case is known in advance so that the results can be compared and the convergence rate is verified. The following error norm is used to check the accuracy of the proposed scheme:

E_{\infty}(K_{t},K_{x})=\max_{1\leq i\leq K_{x}-1}\left|u(x_{i},t_{K_{t}})-u_{i}^{K_{t}}\right|.

(5.1)

Using the error norm defined above, we denote the numerical convergence orders in time and space by $T_{rate}$ and $S_{rate}$ , respectively, and are defined as

T_{rate}=\log_{2}\bigg{(}\frac{E_{\infty}(K_{t}/2,K_{x})}{E_{\infty}(K_{t},K_{x})}\bigg{)},\ \ S_{rate}=\log_{2}\bigg{(}\frac{E_{\infty}(K_{t},K_{x}/2)}{E_{\infty}(K_{t},K_{x})}\bigg{)}.

(5.2)

Here, $K_{x}$ and $K_{t}$ are correspondingly the total numbers of subintervals in space and time directions.

We begin by analysing the error of the considered problem (1.1)-(1.3) through an example having a smooth solution, i.e., $u(\cdot,t)\in C^{4}[0,T]$ .

Example 1.

Let us consider TFDE problem (1.1) with the following data:

	$\displaystyle u$	$\displaystyle(0,t)=u(\pi,t)=0,~{}~{}~{}0<t\leq 1,$
		$\displaystyle u(x,0)=0,~{}~{}~{}0\leq x\leq\pi.$

Here $f(x,t)=t^{5}\sin x\big{(}1+\frac{120}{\Gamma(6-\alpha)}t^{-\alpha}\big{)}$ and thus, the exact solution is $u(x,t)=t^{5}\sin x$ .

Clearly, the above example has a smooth solution throughout the time interval i.e., $u(\cdot,t)\in C^{4}[0,T]$ . Therefore we have calculated the errors on a uniform mesh, i.e., $r=1$ , using the relation (5.1) at final time $T=1$ . Numerical results are given in Table 2 and 3 for different numbers of temporal splits with two different values of diffusion coefficient $\varrho=1,30$ ; fixing $K_{x}=10000$ . The results indicate that the maximum nodal errors are minimal and the convergence order of scheme (4.1)-(4.2) for a smooth solution is $(4-\alpha)$ , which is the expected rate as discussed in Remark 2. It is also observed that as the fractional order $\alpha$ increases, simultaneously, the magnitudes of the maximum errors increase and the convergence rate decreases. One can note in Table 3 that with a high diffusion coefficient, the errors are relatively smaller, which specifies that a higher diffusion coefficient fastens the physical process.

	$\alpha=0.3$		$\alpha=0.6$		$\alpha=0.8$
$K_{t}$	$E_{\infty}$	Order	$E_{\infty}$	Order	$E_{\infty}$	Order
$10$	$4.3693e-04$		$1.9000e-03$		$4.3000e-03$
$20$	$3.8295e-05$	$3.5122$	$2.0232e-04$	$3.2313$	$5.2108e-04$	$3.0448$
$40$	$3.1936e-06$	$3.5839$	$2.0257e-05$	$3.3201$	$5.9698e-05$	$3.1258$
$80$	$2.6190e-07$	$3.6081$	$1.9790e-06$	$3.3556$	$6.6679e-06$	$3.1624$
$160$	$2.0187e-08$	$3.6975$	$1.9197e-07$	$3.3658$	$7.3596e-07$	$3.1795$

Table 2: Maximum temporal errors and order of convergence

(T_{rate})

for

r=1

and

\varrho=1

	$\alpha=0.3$		$\alpha=0.6$		$\alpha=0.8$
$K_{t}$	$E_{\infty}$	Order	$E_{\infty}$	Order	$E_{\infty}$	Order
$10$	$3.0835e-05$		$1.4990e-04$		$3.6707e-04$
$20$	$2.6811e-06$	$3.5237$	$1.5458e-05$	$3.2776$	$4.3006e-05$	$3.0934$
$40$	$2.3099e-07$	$3.5369$	$1.5369e-06$	$3.3303$	$4.8561e-06$	$3.1467$
$80$	$1.9091e-08$	$3.5969$	$1.5544e-07$	$3.3056$	$5.4559e-07$	$3.1539$
$160$	$1.4818e-09$	$3.6875$	$1.5646e-08$	$3.3125$	$6.6199e-08$	$3.0421$

Table 3: Maximum temporal errors and order of convergence

(T_{rate})

for

r=1

and

\varrho=30

The graphical illustration of comparison between the exact solution and numerical solution of Example 1 on uniform mesh with diffusion coefficient $\varrho=1$ and $\varrho=30$ is given in Figure 1 and Figure 1, respectively.

Refer to caption — Figure 1: Example 1. (a) Exact solution and numerical solution for $\alpha=0.3$ and $\varrho=1$ ; (b) Exact solution and numerical solution for $\alpha=0.3$ and $\varrho=30$ .

Table 4 displays that $\mathcal{O}(h^{2})$ is the rate of convergence in spatial direction by fixing $K_{t}=1000$ and changing the numbers of splits of spatial interval.

	$\alpha=0.3$		$\alpha=0.6$		$\alpha=0.8$
$K_{x}$	$E_{\infty}$	Order	$E_{\infty}$	Order	$E_{\infty}$	Order
$10$	$3.1000e-03$		$2.1000e-03$		$1.6000e-03$
$20$	$7.6951e-04$	$2.0103$	$5.3696e-04$	$1.9675$	$4.0646e-04$	$1.9769$
$40$	$1.9239e-04$	$1.9999$	$1.3427e-04$	$1.9997$	$1.0165e-04$	$1.9995$
$80$	$4.8097e-05$	$2.0000$	$3.3570e-05$	$1.9999$	$2.5416e-05$	$1.9998$
$160$	$1.2024e-05$	$2.0000$	$8.3930e-06$	$1.9999$	$6.3558e-06$	$1.9996$

Table 4: Maximum spatial errors and order of convergence

(S_{rate})

for

r=1

and

\varrho=1

Next, we perform another numerical examination for TFDE problem (1.1)-(1.3) having a non-smooth behavior of the solution at initial time moment, i.e., $u(\cdot,t)\notin C^{4}[0,T]$ . We observe the maximum nodal errors for different numbers of time intervals and confirm our theoretical convergence analysis. We also analyse the spatial convergence so that $\mathcal{O}(h^{2}+\tau^{\min\{4-\alpha,r\alpha\}})$ is proved numerically as the rate of convergence for the proposed scheme (4.2)-(4.4).

Example 2.

Let us consider the TFDE problem (1.1) with following data:

	$\displaystyle u$	$\displaystyle(0,t)=u(1,t)=0,~{}~{}~{}0<t\leq 1,$
		$\displaystyle u(x,0)=0,~{}~{}~{}0\leq x\leq 1.$

Here $f(x,t)=t^{2}\sin\pi x\big{(}\frac{\Gamma(3+\alpha)}{2}+\pi^{2}t^{\alpha}\big{)}$ , and $u(x,t)=t^{2+\alpha}\sin\pi x$ is the exact solution. We have considered the value of diffusion coefficient $\varrho$ to be 1 in this example. It is evident here that the exact solution has an initial-time singular behaviour, i.e., $u(\cdot,t)\notin C^{4}[0,T]$ .

We numerically solve the problem by keeping $K_{x}=10000$ fixed and taking different splits of temporal intervals on a uniform $(r=1)$ and graded mesh $(r=\frac{4-\alpha}{\alpha})$ using our newly developed scheme (4.2)-(4.4). The maximum nodal errors (5.1) and convergence orders (5.2) using the proposed difference scheme are displayed in tables below at final time $t_{K_{t}}=T=1$ . Results in Table 5 are calculated on the uniform mesh, which shows that the errors are minimal, but the convergence rate is not desired. It is clear from the results that Example 2 with a non-smooth solution on uniform mesh has error convergence of order $(\alpha)$ , which is much less than the order of convergence $(4-\alpha)$ of the previous example having a smooth solution.

In Table 6, the numerical results agree precisely with the theoretical rate of convergence proved in Lemma 1 using our proposed numerical scheme (4.2)-(4.4). We can see that the maximum errors are nominal and the optimal convergence order of $(4-\alpha)$ is achieved for the choice of grading constant $r=(4-\alpha)/\alpha$ . In fact, for a choice of $r\geq(4-\alpha)/\alpha$ , we can always get the optimal convergence order of $(4-\alpha)$ .

	$\alpha=0.3$		$\alpha=0.6$		$\alpha=0.8$
$K_{t}$	$E_{\infty}$	Order	$E_{\infty}$	Order	$E_{\infty}$	Order
$10$	$2.6610e-06$		$2.5619e-06$		$5.4254e-07$
$20$	$3.0269e-07$	$3.1361$	$3.0002e-07$	$3.0932$	$1.1337e-07$	$2.2587$
$40$	$3.9856e-08$	$2.9250$	$4.0843e-08$	$2.8769$	$2.6454e-08$	$2.0995$
$80$	$1.2020e-08$	$1.7294$	$1.0281e-08$	$1.9901$	$8.8073e-09$	$1.5867$
$160$	$9.3683e-09$	$0.3596$	$6.6545e-09$	$0.6276$	$5.0582e-09$	$0.8001$

Table 5: Maximum temporal errors and order of convergence

(T_{rate})

for

r=1

	$\alpha=0.4$		$\alpha=0.6$		$\alpha=0.8$
$K_{t}$	$E_{\infty}$	Order	$E_{\infty}$	Order	$E_{\infty}$	Order
$10$	$9.7738e-04$		$7.2723e-04$		$4.3321e-04$
$20$	$1.4429e-04$	$2.7599$	$9.7544e-05$	$2.8983$	$5.8820e-05$	$2.8807$
$40$	$1.5987e-05$	$3.1740$	$1.0963e-05$	$3.1534$	$7.1073e-06$	$3.0489$
$80$	$1.5462e-06$	$3.3701$	$1.1312e-06$	$3.2767$	$8.0959e-07$	$3.1340$
$160$	$1.2864e-07$	$3.5873$	$1.0580e-07$	$3.4184$	$8.4873e-08$	$3.2538$

Table 6: Maximum temporal errors and order of convergence

(T_{rate})

for

r=\frac{4-\alpha}{\alpha}

The graphical illustration of comparison between the exact solution and numerical solution of Example 2 with diffusion coefficient $\varrho=1$ on uniform mesh $r=1$ and graded mesh $r=\frac{4-\alpha}{\alpha}$ is given in Figure 2 and Figure 2, respectively.

Table 7 is the display of convergence in the spatial direction which shows that $\mathcal{O}(h^{2})$ convergence rate is attained. This evaluation is done keeping $K_{t}=1000$ fixed and taking different splits of space interval.

	$\alpha=0.3$		$\alpha=0.6$		$\alpha=0.8$
$K_{x}$	$E_{\infty}$	Order	$E_{\infty}$	Order	$E_{\infty}$	Order
$10$	$7.2000e-03$		$6.9000e-03$		$6.6000e-03$
$20$	$1.8000e-03$	$2.0000$	$1.7000e-03$	$2.0211$	$1.6000e-03$	$2.0444$
$40$	$4.4627e-04$	$2.0120$	$4.3120e-04$	$1.9791$	$4.1171e-04$	$1.9584$
$80$	$1.1158e-04$	$1.9998$	$1.0778e-04$	$2.0003$	$1.0291e-04$	$2.0002$
$160$	$2.7924e-05$	$1.9985$	$2.6944e-05$	$2.0001$	$2.5727e-05$	$2.0000$

Table 7: Maximum spatial errors and order of convergence

(S_{rate})

for

r=1

6 Conclusion

In this paper, we derive a new high order approximation to the Caputo fractional derivative of order $\alpha$ , $0<\alpha<1$ , over non-uniform mesh. The algorithm applied here is based on the approximation of the integrand of the Caputo derivative by the Lagrange interpolating polynomials. This algorithm relies on the solution being continuously differentiable up to fourth order. However, the considered Caputo-type TFDE problem (1.1)-(1.3) typically has a solution that exhibits initial-time singularity, i.e. $u(\cdot,t)\notin C^{4}[0,T]$ . Therefore the discretization on graded mesh produces a local truncation error of order $\min\{4-\alpha,r\alpha\}$ , $r$ being the mesh grading parameter. The quantity $r=\frac{(4-\alpha)}{\alpha}$ is observed as the optimal grading parameter. Next, a high-order difference scheme for TFDE problem (1.1)-(1.3) is developed. The proposed scheme is proved to be unconditionally stable and convergent on the uniform mesh $(r=1)$ . It is also found that the proposed scheme is uniquely solvable. At present, we are unable to provide proof for the stability of the scheme for $r>1$ due to the complexity of coefficients, but we are committed to making efforts in that direction. Numerical experiments are carried out to verify the theoretical results for smooth and non-smooth solutions. The results of the second example demonstrate clearly that the proposed scheme for a non-smooth solution on a uniform mesh converges with an order ( $\alpha$ ), which is much lower than the order of convergence in the first example, i.e. ( $4-\alpha$ ), where the solution was smooth. Also, our scheme for a non-smooth solution reflects better convergence on the graded mesh rather than the uniform mesh. This scheme outperforms all existing schemes for graded mesh in terms of accuracy and efficiency. The proposed scheme is developed for a general non-uniform mesh and can be applied and analyzed for any kind of grading over meshes.

Declarations

Funding and acknowledgements

The first author acknowledges the support provided by the Council of Scientific and Industrial Research, India, under grant number 09/086(1483)/2020-EMR-I. The fourth author acknowledges the support provided by the SERB, a statutory body of DST, India, under the award SERB–POWER fellowship (grant number SPF/2021/000103).

Conflict of interest

The authors declare that there are no financial and non-financial competing interests that are relevant to the content of this article.

References

[1] A. A. Alikhanov. A new difference scheme for the time fractional diffusion equation. Journal of Computational Physics, 280:424–438, 2015.
[2] J. Cao, C. Li, and Y. Q. Chen. High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations (II). Fractional Calculus and Applied Analysis, 18(3):735–761, 2015.
[3] G. Gao, Z. Sun, and H. Zhang. A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications. Journal of Computational Physics, 259:33–50, 2014.
[4] R. Hilfer. Fractional calculus and regular variation in thermodynamics. In Applications of fractional calculus in physics, pages 429–463. World Scientific, 2000.
[5] A. Iomin, S. Dorfman, and L. Dorfman. On tumor development: fractional transport approach. arXiv preprint q-bio/0406001, 2004.
[6] B. Jin, R. Lazarov, and Z. Zhou. An analysis of the ${L1}$ scheme for the subdiffusion equation with nonsmooth data. IMA Journal of Numerical Analysis, 36(1):197–221, 2016.
[7] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo. Theory and applications of fractional differential equations. Elsevier Science Limited, 2006.
[8] N. Kumar and M. Mehra. Legendre wavelet collocation method for fractional optimal control problems with fractional Bolza cost. Numerical Methods for Partial Differential Equations, 37(2):1693–1724.
[9] N. Kumar and M. Mehra. Collocation method for solving nonlinear fractional optimal control problems by using Hermite scaling function with error estimates. Optimal Control Applications and Methods, 42(2):417–444, 2021.
[10] C. Li, A. Chen, and J. Ye. Numerical approaches to fractional calculus and fractional ordinary differential equation. Journal of Computational Physics, 230(9):3352–3368, 2011.
[11] C. Li and F. Zeng. Numerical methods for fractional calculus. Chapman and Hall/CRC, 2015.
[12] C. P. Li, R. F. Wu, and H. F. Ding. High-order approximation to Caputo derivative and Caputo-type advection-diffusion equation (1), communications in applied and industrial mathematics, in press. 2014.
[13] Y. Lin and C. Xu. Finite difference/spectral approximations for the time-fractional diffusion equation. Journal of Computational Physics, 225(2):1533–1552, 2007.
[14] Y. Luchko. Anomalous diffusion: models, their analysis, and interpretation. In Advances in Applied Analysis, pages 115–145. Springer, 2012.
[15] F. Mainardi and P. Paradisi. A model of diffusive waves in viscoelasticity based on fractional calculus. Decision and Control, Proceedings of the 36th IEEE Conference, 5:4961–4966, 1997.
[16] V. Mehandiratta and M. Mehra. A difference scheme for the time-fractional diffusion equation on a metric star graph. Applied Numerical Mathematics, 158:152–163, 2020.
[17] V. Mehandiratta, M. Mehra, and G. Leugering. Fractional optimal control problems on a star graph: optimality system and numerical solution. Mathematical Control & Related Fields, 11(1):189, 2021.
[18] V. Mehandiratta, M. Mehra, and G. Leugering. Optimal control problems driven by time-fractional diffusion equations on metric graphs: Optimality system and finite difference approximation. SIAM Journal on Control and Optimization, 59(6):4216–4242, 2021.
[19] K. S. Patel and M. Mehra. Fourth order compact scheme for space fractional advection–diffusion reaction equations with variable coefficients. Journal of Computational and Applied Mathematics, 380:112963, 2020.
[20] I. Podlubny. Fractional differential equations, volume 198. Academic Press, 1998.
[21] Haili Qiao and Aijie Cheng. A fast high order method for time fractional diffusion equation with non-smooth data. Discrete and Continuous Dynamical Systems-B, 27(2):903–920, 2022.
[22] A. Shukla, M. Mehra, and G. Leugering. A fast adaptive spectral graph wavelet method for the viscous Burgers’ equation on a star-shaped connected graph. Mathematical Methods in the Applied Sciences, 43(13):7595–7614, 2020.
[23] A. K. Singh and M. Mehra. Uncertainty quantification in fractional stochastic integro-differential equations using Legendre wavelet collocation method. In International Conference on Computational Science, pages 58–71. Springer, 2020.
[24] A. K. Singh and M. Mehra. Wavelet collocation method based on Legendre polynomials and its application in solving the stochastic fractional integro-differential equations. Journal of Computational Science, 51:101342, 2021.
[25] Abhishek Kumar Singh, Mani Mehra, and Samarth Gulyani. Learning parameters of a system of variable order fractional differential equations. Numerical Methods for Partial Differential Equations, 39(3):1962–1976, 2023.
[26] Abhishek Kumar Singh, Mani Mehra, and Samarth Gulyani. A modified variable-order fractional sir model to predict the spread of covid-19 in india. Mathematical Methods in the Applied Sciences, 46(7):8208–8222, 2023.
[27] M. Stynes, E. O’Riordan, and J. L. Gracia. Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. SIAM J. Numer. Anal., 55(2):1057–1079, 2017.

	$\displaystyle Tv^{1}$	$\displaystyle=-g_{1}v^{0},$
	$\displaystyle\\|v^{1}\\|_{\infty}$	$\displaystyle\leq\left\\|T^{-1}\right\\|_{\infty}\ \left\\|\left(-g_{1}v^{0}\right)\right\\|_{\infty},$
		$\displaystyle\leq\|-g_{1}\|\ \left\\|v^{0}\right\\|_{\infty}.$

	$\displaystyle\left\\|v^{1}\right\\|_{\infty}$	$\displaystyle\leq\|g_{0}\|\ \left\\|v^{0}\right\\|_{\infty},$
	$\displaystyle\left\\|v^{1}\right\\|_{\infty}$	$\displaystyle\leq\left\\|v^{0}\right\\|_{\infty}.$

	$\displaystyle Tv^{n+1}$	$\displaystyle=-{\sum_{k=1}^{n+1}}g_{k}{v^{n+1-k}},$
	$\displaystyle\left\\|v^{n+1}\right\\|_{\infty}$	$\displaystyle=\left\\|T^{-1}\left(-\sum_{k=1}^{n+1}{g_{k}}v^{n+1-k}\right)\right\\|_{\infty},$
		$\displaystyle\leq\left\\|T^{-1}\right\\|_{\infty}\ \left\\|\left(-\sum_{k=1}^{n+1}{g_{k}}v^{n+1-k}\right)\right\\|_{\infty},$
		$\displaystyle\leq\sum_{k=1}^{n+1}\|{g_{k}}\|\ \left\\|v^{n+1-k}\right\\|_{\infty},$
		$\displaystyle\leq C\left\\|v^{0}\right\\|_{\infty}\ \left(\sum_{k=1}^{n+1}\|{g_{k}}\|\right),$

	$\displaystyle\|r_{i}^{n}\|=\|L_{K_{x},K_{t}}e_{i}^{n}\|$	$\displaystyle=\big{\|}\textsubscript{C}D_{0,K_{t}}^{\alpha}e_{i}^{n}-\delta_{x}^{2}e_{i}^{n}\big{\|},$
		$\displaystyle=\Big{\|}\textsubscript{C}D_{0,K_{t}}^{\alpha}\left[u(x_{i},t_{n})-u_{i}^{n}\right]-\delta_{x}^{2}\left[u(x_{i},t_{n})-u_{i}^{n}\right]\Big{\|},$
		$\displaystyle\hskip 2.84544pt\leq C\big{(}{h^{2}}+n^{-\min\{4-\alpha,r\alpha\}}\big{)},$

	$\displaystyle\|e_{l}^{1}\|$	$\displaystyle\leq(g_{0}+2\eta\epsilon)\|e_{l}^{1}\|-\eta\epsilon\|e_{l+1}^{1}\|-\eta\epsilon\|e_{l-1}^{1}\|,$
		$\displaystyle\leq\|(g_{0}+2\eta\epsilon)e_{l}^{1}-\eta\epsilon e_{l+1}^{1}-\eta\epsilon e_{l-1}^{1}\|,$
		$\displaystyle=\|r_{l}^{1}\|\leq C({h^{2}}+\tau^{\alpha}),$
	$\displaystyle\\|e^{1}\\|_{\infty}$	$\displaystyle\leq C({h^{2}}+\tau^{\alpha}).$