gL1 Scheme for Solving a Class of Generalized Time-Fractional Diffusion Equations

Abstract

In this paper, a numerical scheme based on a general temporal mesh is constructed for a generalized time-fractional diffusion problem of order $\alpha$. The main idea involves the generalized linear interpolation and so we term the numerical scheme the gL1 scheme. The stability and convergence of the numerical scheme are analyzed using the energy method. It is proven that the temporal convergence order is $(2-\alpha )$ for a general temporal mesh. Simulation is carried out to verify the efficiency of the proposed numerical scheme.

1. Introduction

The fractional model has been shown to be a powerful tool in modeling various memory process in applications [1,2], such as diffusion process in a porous media. The essence of this tool lies in fractional derivative (or integral) which is an extension of integer derivative (or integral). Many well known mathematicians such as Euler, Lagrange, Liouville, Riemann, Grüwald, Letnikov and Caputo have devoted their efforts to fractional calculus and have made great contributions to this topic. The extension from integer derivative (or integral) to fractional derivative (or integral) may not be unique due to the different techniques applied. To unify the different approaches, Agrawal [3] has proposed some generalized fractional operators which may unify some well known fractional operators such as Caputo, Riemann-Louville and Hadamard. The generalized fractional operators incorporate a scale function $z\left(t\right)$ and a weight function $w\left(t\right)$. With these two functions, many equations can be written in a general form and thus can be solved in an elegant way as shown in [3]. Moreover, by choosing different $z\left(t\right)$ and $w\left(t\right)$, one can readily obtain the well known fractional operators.

The generalized fractional operators naturally lead to generalized fractional problems. As in the fractional situation, the analytical solutions may not be easily derived and hence the corresponding numerical solutions are both necessary and useful in applications. In fact, some pioneer works have been done on the numerical treatment of certain generalized fractional problems [4,5,6,7,8,9,10]. In the earlier papers [8,9,10], the convergence of the numerical scheme is established using the Lax–Richtmyer theorem but the order of convergence is not explicitly given. Inspired by [11,12], some generalized weighted shifted Grünwald-Letnikov (gWSGL) type approximations [4,5,6] and generalized Alikhanov’s approximation [7] were recently proposed, which improve the accuracy of previous work. In fact, the convergence order of these methods are shown to be $O\left({\tau }_z^2\right)$ (or higher) based on a particular choice of the temporal mesh that closely depends on the scale function $z\left(t\right)$. It is noteworthy that the higher convergence order comes at the expense of more computations and restrictions in reality.

In this paper, we shall consider a generalized time-fractional diffusion problem that is more general than that in [9]. As one of the pioneer numerical treatments for generalized time-fractional diffusion problems, the authors in [9] first present analytical solution of a generalized time-fractional problem involving second order spatial derivative. Then, based on a uniform temporal mesh, they construct a numerical scheme by finite difference method. The stability of the numerical scheme is investigated via an estimate of the inverse of the coefficient matrix, and the convergence then follows from the Lax–Richtmyer theorem without giving explicit convergence order. Different from the work of [9], we aim to derive a numerical scheme based on a more general temporal mesh than [9] and give the convergence order of the proposed scheme explicitly using energy method. Our major contributions in this paper are as follows:

• We consider a problem involving an operator $\mathcal{L}$ and propose a numerical scheme based on a general temporal mesh using generalized linear interpolation.
• We establish the stability and convergence of the proposed scheme, which is based on a general temporal mesh, via energy method. The analysis in the context of general temporal mesh is not trivial.

We consider the following generalized time-fractional diffusion equation with weight function $w\left(t\right)\equiv 1$

$$\left\{\begin{array}{c}{\displaystyle {}_0^C{D}_{t;\left[z\right(t),1]}^{\alpha }u(x,t)=\mathcal{L}u(x,t)+f(x,t),(x,t)\in (0,1)\times (0,1)}\hfill \\ \hfill \\ {\displaystyle u(x,0)=\psi \left(x\right),x\in [0,1]}\hfill \\ \hfill \\ {\displaystyle u(0,t)={\varphi }_1\left(t\right),u(1,t)={\varphi }_2\left(t\right),t\in (0,1]}\hfill \end{array}\right.$$

(1)

where $\mathcal{L}$ is a linear operator defined by

$$\mathcal{L}u(x,t)={\left(p(x,t)u_x(x,t)\right)}_x-q(x,t)u(x,t),$$

with $p(x,t)\ge p_0>0$ and $q(x,t)\ge 0$, ${}_0^C{D}_{t;\left[z\right(t),1]}^{\alpha }u(x,t)(0<\alpha <1)$ is the generalized Caputo fractional derivative given by

$${}_0^C{D}_{t;\left[z\right(t),1]}^{\alpha }u(x,t)=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{1}{{[z\left(t\right)-z\left(s\right)]}^{\alpha }}\frac{\partial u(x,s)}{\partial s}ds,$$

(2)

and ${\varphi }_1\left(t\right),{\varphi }_2\left(t\right),\psi \left(x\right),f(x,t)$ are given functions that are sufficiently smooth. We remark that (i) a generalized fractional problem of type (1) with any weight function $w\left(t\right)$ (not necessarily 1) can be converted to (1) by a simple formula $u(x,t)=w\left(t\right)v(x,t)$ (see [4,5,6,7]); and (ii) the generalized fractional equation considered in [9] is a special case of (1) when $p(x,t)\equiv 1$ and $q(x,t)\equiv 0.$

The plan of the paper is as follows. In Section 2, we shall develop a numerical scheme for the problem (1) based on a more general temporal mesh than [9]. Then, the stability as well as the convergence of the proposed numerical scheme will be established rigorously using energy method in Section 3. In Section 4, we carry out experiments to verify as well as to demonstrate the efficiency of the proposed numerical scheme. Finally, a brief conclusion is given in Section 5.

2. Numerical Scheme

In this section, we shall derive a numerical scheme for the problem (1) using the key idea of generalized linear interpolation. To begin, let

$$\Delta :0=t_0<t_1<\cdots <t_{N-1}<t_N=1 \mathrm{and} {\Delta }^{\prime }:0=x_0<x_1<\cdots <x_{M-1}<x_M=1$$

be any mesh in the temporal dimension and a uniform mesh in the spatial dimension with step size $h=\frac{1}{M}$, respectively. Throughout, assume that the scale function $z\left(t\right)$ is a strictly increasing function. Let $Z:t\to z\left(t\right)$ be a map from $[0,1]$ to $\left[z\right(0),z(1\left)\right]$. Denote $z\left(t_n\right)=z_n.$ Moreover, denote by $U_j^n$ the exact solution of (1) at $(x_j,t_n)$, and by $u_j^n$ an approximation of $U_j^n$.

$$L_{j,k-1}\left(z\left(t\right)\right)=\frac{z\left(t\right)-z_{k-1}}{z_k-z_{k-1}}{u}_j^k+\frac{z\left(t\right)-z_k}{z_{k-1}-z_k}{u}_j^{k-1}.$$

(3)

Note that $L_{j,k-1}\left(z\left(t\right)\right)$ is a generalized linear polynomial of $z\left(t\right)$ and it will be used to approximate $u(x_j,t)$ over the interval $[t_{k-1},t_k]$ in (2). Indeed, we derive the following approximation scheme at $(x_j,t_n)$

$$\begin{array}{cc}{\displaystyle {}_0^C{D}_{t;\left[z\right(t),1]}^{\alpha }u(x_j,t_n)}\hfill & {\displaystyle =\frac{1}{\Gamma (2-\alpha )}\sum _{k=1}^n{\int }_{t_{k-1}}^{t_k}\frac{1}{{[z\left(t_n\right)-z\left(s\right)]}^{\alpha }}\frac{\partial u(x,s)}{\partial s}ds}\hfill \\ \hfill \\ \hfill & {\displaystyle \approx \frac{1}{\Gamma (2-\alpha )}\sum _{k=1}^n{\int }_{t_{k-1}}^{t_k}\frac{1}{{[z\left(t_n\right)-z\left(s\right)]}^{\alpha }}{\left[L_{j,k-1}\left(z\left(s\right)\right)\right]}^{\prime }ds}\hfill \\ \hfill \\ \hfill & {\displaystyle =\frac{1}{\Gamma (2-\alpha )}\sum _{k=1}^n{\int }_{t_{k-1}}^{t_k}\frac{1}{{[z\left(t_n\right)-z\left(s\right)]}^{\alpha }}\frac{u_j^k-u_j^{k-1}}{z_k-z_{k-1}}{z}^{\prime }\left(s\right)ds}\hfill \\ \hfill \\ \hfill & {\displaystyle =\mu \left[{\omega }_z^0{u}_j^n-\sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)u_j^k-{\omega }_z^{n-1}{u}_j^0\right]:=gL1\left[u(x_j,t_n)\right],}\hfill \end{array}$$

(4)

where $\mu =\frac{1}{\Gamma (2-\alpha )}$ and the coefficients

$${\omega }_z^{n-k-1}=\frac{{(z_n-z_k)}^{1-\alpha }-{(z_n-z_{k+1})}^{1-\alpha }}{z_{k+1}-z_k},0\le k\le n-1.$$

(5)

We shall call (4) the gL1 approximation of the generalized Caputo fractional derivative.

Remark 1.

(a) 

Note that (4) is derived for any temporal mesh Δ and this is an extension of [9] which considers uniform temporal mesh. Moreover, the technique used to obtain (4) involves generalized linear interpolation whereas in [9] finite difference is used to approximate the derivatives in the integrand of the fractional derivative.

(b) 

The Formula (4) is also different from the nonuniform $L1$ formula in [13] since a scale function $z\left(t\right)$ such that $z_j\ne z_i,i\ne j$ is considered.

To proceed further, we shall investigate the accuracy of the $gL1$ approximation (4) and the properties of the coefficients ${\omega }_z^k$ which are both vital in subsequent analysis. For the former one, we introduce the following definition.

Definition 1

([14]). Given the mesh $\Delta :0=t_0<t_1<\cdots <t_{N-1}<t_N=1,$ denote $z\left(t_k\right)=z_k.$ The mesh

$${\Delta }_z:z\left(0\right)=z_0<z_1<\cdots <z_{N-1}<z_N=z\left(1\right)$$

is said to be quasi uniform if

$$\frac{{\tau }_{z,max}}{{\tau }_{z,min}}\le \rho ,$$

where ${\tau }_{z,max}={max}_{1\le k\le N}|z_k-z_{k-1}|$, ${\tau }_{z,min}={min}_{1\le k\le N}|z_k-z_{k-1}|$ and $\rho >0$ is a constant.

Theorem 1

(gL1 approximation). Assume that for any fixed $x=x_j$, $u(x_j,Z^{-1}\left(z\right))=g\left(z\right)\in C^2[z\left(0\right),z\left(1\right)]$. Suppose that the mesh ${\Delta }_z:z\left(0\right)=z_0<z_1<\cdots <z_{N-1}<z_N=z\left(1\right)$ is quasi uniform. Then, we have for any fixed $\alpha \in (0,1)$,

$${}_0^C{D}_{t;\left[z\right(t),1]}^{\alpha }u(x_j,t_n)=gL1\left[u(x_j,t_n)\right]+O\left({\tau }_{z,max}^{2-\alpha }\right).$$

(6)

Proof. 

From (4) we see that the error term $R_j^n$ satisfies

$${}_0^C{D}_{t;\left[z\right(t),1]}^{\alpha }u(x_j,t_n)=gL1\left[u(x_j,t_n)\right]+R_j^n,$$

and

$$R_j^n=\frac{1}{\Gamma (1-\alpha )}\sum _{k=1}^n{\int }_{t_{k-1}}^{t_k}{[z\left(t_n\right)-z\left(s\right)]}^{-\alpha }\left\{\frac{\partial u(x_j,s)}{\partial s}-{\left[L_{j,k-1}\left(z\left(s\right)\right)\right]}^{\prime }\right\}ds.$$

Noting that $L_{j,k-1}\left(z\left(t_k\right)\right)=u_j^k=u(x_j,t_k)$ and the relation $Z^{-1}\left(z\right)=t$, after applying integration by parts, we get

$$R_j^n=-\frac{\alpha }{\Gamma (1-\alpha )}\sum _{k=1}^n{\int }_{z_{k-1}}^{z_k}{(z_n-z)}^{-\alpha -1}\left[u(x_j,Z^{-1}\left(z\right))-L_{j,k-1}\left(z\right)\right]dz.$$

(7)

Since $u(x_j,Z^{-1}\left(z\right))=g\left(z\right)\in C^2[z\left(0\right),z\left(1\right)]$, it is well known that

$$u(x_j,Z^{-1}\left(z\right))-L_{j,k-1}\left(z\right)=(z-z_{k-1})(z-z_k)\frac{g_{zz}\left({\xi }_k\right)}{2},z_{k-1}<{\xi }_k<z_k.$$

(8)

Denote $M_g={max}_{z\left(0\right)\le z\le z\left(1\right)}\left|g_{zz}\left(z\right)\right|$. Upon substituting (8) into (7), we find

$$\begin{array}{cc}|R_j^n|\hfill & {\displaystyle \le \frac{\alpha M_g}{2\Gamma (1-\alpha )}\sum _{k=1}^n{\int }_{z_{k-1}}^{z_k}{(z_n-z)}^{-\alpha -1}(z-z_{k-1})(z_k-z)dz}\hfill \\ \hfill \\ & {\displaystyle =\frac{M_g}{2\Gamma (2-\alpha )}\left[\sum _{k=1}^{n-1}{\int }_{z_{k-1}}^{z_k}{(z_n-z)}^{-\alpha -1}(z-z_{k-1})(z_k-z)dz+{\int }_{z_{n-1}}^{z_n}\frac{z-z_{n-1}}{{(z_n-z)}^{\alpha }}dz\right]}\hfill \\ \hfill \\ & {\displaystyle \le \frac{M_g}{2\Gamma (2-\alpha )}\left[\sum _{k=1}^{n-1}\frac{{(z_k-z_{k-1})}^2}{4}{\int }_{z_{k-1}}^{z_k}{(z_n-z)}^{-\alpha -1}dz+\frac{{(z_n-z_{n-1})}^{2-\alpha }}{1-\alpha }\right].}\hfill \end{array}$$

This further gives

$$\begin{array}{cc}{\displaystyle |R_j^n|}\hfill & {\displaystyle \le \frac{M_g}{2\Gamma (2-\alpha )}\left[\frac{{\tau }_{z,max}^2}{4}{\int }_{z_0}^{z_{n-1}}{(z_n-z)}^{-\alpha -1}dz+\frac{{\tau }_{z,max}^{2-\alpha }}{1-\alpha }\right]}\hfill \\ \hfill \\ & {\displaystyle =\frac{M_g}{2\Gamma (2-\alpha )}\left[\frac{{\tau }_{z,max}^2}{4\alpha }\left({(z_n-z_{n-1})}^{-\alpha }-{(z_n-z_0)}^{-\alpha }\right)+\frac{{\tau }_{z,max}^{2-\alpha }}{1-\alpha }\right]}\hfill \\ \hfill \\ & {\displaystyle \le \frac{M_g}{2\Gamma (2-\alpha )}\left[\frac{{\tau }_{z,max}^2}{4\alpha }{(z_n-z_{n-1})}^{-\alpha }+\frac{{\tau }_{z,max}^{2-\alpha }}{1-\alpha }\right]}\hfill \\ \hfill \\ & {\displaystyle \le \frac{M_g}{2\Gamma (2-\alpha )}\left[\frac{{\rho }^{\alpha }}{4\alpha }+\frac{1}{1-\alpha }\right]{\tau }_{z,max}^{2-\alpha }.}\hfill \end{array}$$

Hence, for any fixed $\alpha \in (0,1)$, we get

$$|R_j^n|=O\left({\tau }_{z,max}^{2-\alpha }\right).$$

This completes the proof. □

Remark 2.

There is some relation between the commonly used uniform mesh of t and the quasi uniform mesh of $z\left(t\right)$.

(a)

For finite intervals of t, say $t\in [0,1],$ in practice the partition of $[0,1]$ always results in finite number of subintervals, so we are able to find a constant ρ such that $\frac{{\tau }_{z,max}}{{\tau }_{z,min}}\le \rho .$ Hence, it is clear that any general mesh of t yields a quasi uniform mesh of $z\left(t\right)$. In particular, we can say that for finite intervals, the commonly used uniform mesh of t is a special case of the quasi uniform mesh of $z\left(t\right)$.

(b)

For infinite intervals of t, in general the uniform mesh of t may not yield the quasi uniform mesh of $z\left(t\right)$. However, if $0<c<|z^{\prime }\left(t\right)|<C$ for all t in the infinite interval, then the uniform mesh of t will give the quasi uniform mesh of $z\left(t\right)$, and so we can conclude here that the uniform mesh of t is a special case of the quasi uniform mesh of $z\left(t\right)$.

Our next result gives the properties of the coefficients ${\omega }_z^k$ in (5) that is vital in subsequent analysis.

Lemma 1.

For fixed n, we have

$$\begin{array}{ccc}& & {\omega }_z^{n-k-1}\ge (1-\alpha ){(z_n-z_k)}^{-\alpha }>0,0\le k\le n-1\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& & {\omega }_z^{n-k-1}\ge {\omega }_z^{n-k},0\le k\le n-1.\hfill \end{array}$$

(10)

Proof. 

Let $F\left(z\right)={(z_n-z)}^{1-\alpha }$. Applying mean value theorem and noting the Definition (5), we get

$$F^{\prime }\left(\xi \right)=-(1-\alpha ){(z_n-\xi )}^{-\alpha }=\frac{F\left(z_{k+1}\right)-F\left(z_k\right)}{z_{k+1}-z_k}=-{\omega }_z^{n-k-1},z_k\le \xi \le z_{k+1}$$

i.e.,

$${\omega }_z^{n-k-1}=(1-\alpha ){(z_n-\xi )}^{-\alpha },z_k\le \xi \le z_{k+1}.$$

(11)

Since ${(z_n-·)}^{-\alpha }$ is an increasing function, (11) immediately leads to (9) and (10). □

We are now ready to construct a numerical scheme for the generalized time-fractional diffusion Equation (1). Discretizing (1) at $(x_j,t_n)$ and using the approximation (4) together with finite difference method in the spatial dimension yields

$$\left\{\begin{array}{c}{\displaystyle \mu \left[{\omega }_z^0{u}_j^n-\sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)u_j^k-{\omega }_z^{n-1}{u}_j^0\right]=\mathrm{\Lambda }{u}_j^n+f_j^n,1\le j\le M-1}\hfill \\ \hfill \\ {\displaystyle u_j^0=\psi \left(x_j\right),1\le j\le M-1}\hfill \\ \hfill \\ u_0^n={\varphi }_1\left(t_n\right),u_M^n={\varphi }_1\left(t_n\right),0\le n\le N\hfill \end{array}\right.$$

(12)

where $\mu =\frac{1}{\Gamma (2-\alpha )}$ and $\mathrm{\Lambda }$ is an operator given by

$$\mathrm{\Lambda }{u}_j^n={\delta }_x{\left(p{\delta }_{x}u\right)}_j^n-q_j^n{u}_j^n,$$

with

$${\delta }_x{\left(p{\delta }_{x}u\right)}_j^n=\frac{1}{h}\left(p_{j+\frac{1}{2}}^n{\delta }_x{u}_{j+\frac{1}{2}}^n-p_{j-\frac{1}{2}}^n{\delta }_x{u}_{j-\frac{1}{2}}^n\right),$$

$p_{j+\frac{1}{2}}^n=p(x_{j+\frac{1}{2}},t_n),q_j^n=q(x_j,t_n)$ and ${\delta }_x{u}_{j+\frac{1}{2}}^n=\left(u_{j+1}^n-u_j^n\right)/h$. We shall call (12) the gL1 scheme of the generalized time-fractional diffusion Equation (1).

Remark 3.

It is easy to verify that the coefficient matrix of the system (12) is strictly diagonally dominated. Therefore, it is uniquely solvable.

3. Theoretical Results

In this section, we shall analyze the stability as well as the convergence of the numerical scheme (12). To begin, let $U_h=\{u=(u_0,u_1,\cdots ,u_M)|u_0=u_M=0\}$ and for any $u,v\in U_h$, define

$${\delta }_x{u}_{j+\frac{1}{2}}=\frac{1}{h}(u_{j+1}-u_j),{\delta }_x{u}_j=\frac{1}{h}\left(u_{j+\frac{1}{2}}-u_{j-\frac{1}{2}}\right),{\delta }_x^2{u}_j=\frac{1}{h}\left({\delta }_x{u}_{j+\frac{1}{2}}-{\delta }_x{u}_{j-\frac{1}{2}}\right)$$

and

$$(u,v)=h\sum _{j=1}^{M-1}{u}_j{v}_j{,\left|u\right|}_1=\sqrt{h\sum _{j=0}^{M-1}{\left({\delta }_x{u}_{j+\frac{1}{2}}\right)}^2}{,\parallel u\parallel }_{\infty }=\underset{1\le j\le M-1}{max}\left|u_j\right|.$$

(13)

Obviously, $\parallel ·\parallel =\sqrt{(·,·)}$ is a norm defined over the space $U_h$.

Next, we shall present three lemmas that will be used later to establish the stability and convergence results. The first one gives a relation between ${\left|u\right|}_1$ and ${\parallel u\parallel }_{\infty }$.

Lemma 2

([15]). For any $u\in U_h$, we have the inequality ${\parallel u\parallel }_{\infty }\le \frac{1}{2}{\left|u\right|}_1.$

The next lemma is from [16] which reveals a relation between ${\left|u\right|}_1$ and $\parallel u\parallel$.

Lemma 3

(Discrete Poincare inequality [16]). Suppose that $u\in U_h$, then

$$\frac{2}{h}sin\left(\frac{\pi h}{2}\right)\parallel u\parallel \le {\left|u\right|}_1.$$

(14)

Remark 4.

If $h<1$ which is the case in our method, from the fact $sin\left(\frac{\pi }{2}h\right)\ge h$, the inequality (14) yields ${\left|u\right|}_1\ge 2\parallel u\parallel$ that will be used in subsequent analysis.

The last lemma gives an explicit expression of $-(\mathrm{\Lambda }u,u)$.

Lemma 4.

For any $u\in U_h$, we have

$$-(\mathrm{\Lambda }u,u)=h\sum _{j=0}^{M-1}{p}_{j+\frac{1}{2}}{\left({\delta }_x{u}_{j+\frac{1}{2}}\right)}^2+h\sum _{j=1}^{M-1}{q}_j{u}_j^2.$$

Proof. 

Using the definition in (13), it is found that

$$\begin{array}{cc}-(\mathrm{\Lambda }u,u)\hfill & {\displaystyle =-h\sum _{j=1}^{M-1}\left[{\delta }_x{\left(p{\delta }_{x}u\right)}_j-q_j{u}_j\right]u_j=-h\sum _{j=1}^{M-1}\left[{\delta }_x{\left(p{\delta }_{x}u\right)}_j\right]u_j+h\sum _{j=1}^{M-1}{q}_j{\left(u_j\right)}^2}\hfill \\ \hfill \\ & {\displaystyle =-\sum _{j=1}^{M-1}\left(p_{j+\frac{1}{2}}{\delta }_x{u}_{j+\frac{1}{2}}-p_{j-\frac{1}{2}}{\delta }_x{u}_{j-\frac{1}{2}}\right)u_j+h\sum _{j=1}^{M-1}{q}_j{\left(u_j\right)}^2}\hfill \\ \hfill \\ & {\displaystyle =-\sum _{j=1}^{M-1}{p}_{j+\frac{1}{2}}{u}_j{\delta }_x{u}_{j+\frac{1}{2}}+\sum _{j=1}^{M-1}{p}_{j-\frac{1}{2}}{u}_j{\delta }_x{u}_{j-\frac{1}{2}}+h\sum _{j=1}^{M-1}{q}_j{\left(u_j\right)}^2}\hfill \\ \hfill \\ & {\displaystyle =\sum _{j=0}^{M-1}{p}_{j+\frac{1}{2}}(u_{j+1}-u_j){\delta }_x{u}_{j+\frac{1}{2}}+h\sum _{j=1}^{M-1}{q}_j{\left(u_j\right)}^2,}\hfill \end{array}$$

where we have used $u_0=u_M=0$ in the last equality. The result is immediate from the above equation. □

Remark 5.

Since $p(x,t)\ge p_0>0$ and $q(x,t)\ge 0$, it is clear from Lemma 4 that $-(\mathrm{\Lambda }u,u)\ge p_0{\left|u\right|}_1^2$ for any $u\in U_h.$

Now, let us present the stability and convergence of the numerical scheme (12). To this aim, we shall first consider (12) with zero boundary conditions.

Theorem 2.

Let $\{u_j^n,1\le j\le M-1,1\le n\le N\}$ be the solution of the system (12) with zero boundary conditions. Then, we have

$${\parallel u^n\parallel }^2+\frac{p_0}{{\omega }_z^0\mu }{\left|u^n\right|}_1^2\le E,1\le n\le N$$

(15)

where $\mu =\frac{1}{\Gamma (2-\alpha )},{\omega }_z^0={(z_n-z_{n-1})}^{-\alpha }, p_0>0$ and

$$E=\parallel u^0{\parallel }^2+\frac{\Gamma (1-\alpha ){(z\left(1\right)-z\left(0\right))}^{\alpha }}{4p_0}\underset{1\le n\le N}{max}{\parallel f^n\parallel }^2.$$

Proof. 

Multiplying both sides of the first equation of (12) by $-u_j^n$ and summing j from 1 to $(M-1)$ gives

$$-\mu \left[{\omega }_z^0{u}_j^n-\sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)u_j^k-{\omega }_z^{n-1}{u}_j^0\right]=-(\mathrm{\Lambda }{u}^n,u^n)-(f^n,u^n),$$

which is rearranged to

$$\mu {\omega }_z^0(u^n,u^n)-(\mathrm{\Lambda }{u}^n,u^n)=\mu \sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)(u^k,u^n)+\mu {\omega }_z^{n-1}(u^0,u^n)+(f^n,u^n).$$

(16)

Since we consider zero boundary conditions here, it is obvious that $u^n\in U_h$. Therefore, using Remark 5, we get a lower bound for the left side of (16) below

$$\mu {\omega }_z^0(u^n,u^n)-(\mathrm{\Lambda }{u}^n,u^n)\ge \mu {\omega }_z^0\parallel u^n{\parallel }^2+p_0{\left|u^n\right|}_1^2.$$

Noting (10) in Lemma 1 and using $xy\le \frac{1}{2}(x^2+y^2)$, an upper bound for the first two terms on the right side of (16) is found as follows

$$\begin{array}{cc}& {\displaystyle \mu \sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)(u^k,u^n)+\mu {\omega }_z^{n-1}(u^0,u^n)}\hfill \\ \hfill \\ \le \hfill & {\displaystyle \frac{\mu }{2}\left[\sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)(u^k,u^k)+{\omega }_z^{n-1}(u^0,u^0)\right]+\frac{\mu }{2}{\omega }_z^0(u^n,u^n)}\hfill \\ \hfill \\ =\hfill & {\displaystyle \frac{\mu }{2}\left[\sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)\parallel u^k{\parallel }^2+{\omega }_z^{n-1}{\parallel u^0\parallel }^2\right]+\frac{\mu }{2}{\omega }_z^0{\parallel u^n\parallel }^2.}\hfill \end{array}$$

For the third term on the right side of (16), the Young’s inequality gives the following upper bound

$$(f^n,u^n)\le \frac{\parallel f^n{\parallel }^2}{4ϵ}+ϵ{\parallel u^n\parallel }^2,\forall ϵ>0.$$

Upon substituting the above upper and lower bounds into (16), we immediately get

$$\frac{\mu }{2}{\omega }_z^0\parallel u^n{\parallel }^2+p_0{\left|u^n\right|}_1^2-ϵ\parallel u^n{\parallel }^2\le \frac{\mu }{2}\left[\sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)\parallel u^k{\parallel }^2+{\omega }_z^{n-1}{\parallel u^0\parallel }^2\right]+\frac{1}{4ϵ}{\parallel f^n\parallel }^2.$$

(17)

Next, using Lemma 3 and noting Remark 4, we have $|u^n{|}_1\ge 2\parallel u^n\parallel$. Then, with $ϵ=2p_0$, the left side of (17) gives the following lower bound

$$\frac{\mu }{2}{\omega }_z^0\parallel u^n{\parallel }^2+p_0{\left|u^n\right|}_1^2-ϵ\parallel u^n{\parallel }^2\ge \frac{\mu }{2}{\omega }_z^0\parallel u^n{\parallel }^2+\left(p_0-\frac{ϵ}{4}\right)|u^n{|}_1^2=\frac{\mu }{2}{\omega }_z^0\parallel u^n{\parallel }^2+\frac{p_0}{2}{\left|u^n\right|}_1^2.$$

Substituting the above into (17) leads to

$${\omega }_z^0\parallel u^n{\parallel }^2+\frac{p_0}{\mu }{\left|u^n\right|}_1^2\le \sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right){\parallel u^k\parallel }^2+{\omega }_z^{n-1}\left(\parallel u^0{\parallel }^2+\frac{1}{4p_0\mu {\omega }_z^{n-1}}{\parallel f^n\parallel }^2\right).$$

(18)

Further, from (9) we see that ${\omega }_z^{n-1}\ge (1-\alpha ){(z_n-z_0)}^{-\alpha }>0$ and hence

$$\frac{1}{{\omega }_z^{n-1}}\le \frac{{(z_n-z_0)}^{\alpha }}{1-\alpha }\le \frac{{(z\left(1\right)-z\left(0\right))}^{\alpha }}{1-\alpha }.$$

Using the above inequality and $\mu =\frac{1}{\Gamma (2-\alpha )}$ in (18), we find

$${\omega }_z^0\left(\parallel u^n{\parallel }^2+\frac{p_0}{{\omega }_z^0\mu }{\left|u^n\right|}_1^2\right)\le \sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right){\parallel u^k\parallel }^2+{\omega }_z^{n-1}E,$$

where E is defined in the theorem. Noting (10), the above inequality readily leads to

$${\omega }_z^0\left(\parallel u^n{\parallel }^2+\frac{p_0}{{\omega }_z^0\mu }{\left|u^n\right|}_1^2\right)\le \sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)\left[\parallel u^k{\parallel }^2+\frac{p_0{(z_k-z_{k-1})}^{\alpha }}{\mu }{\left|u^k\right|}_1^2\right]+{\omega }_z^{n-1}E,$$

or equivalently

$$\begin{array}{cc}\hfill & {\displaystyle {\omega }_z^0\left[\parallel u^n{\parallel }^2+\frac{p_0{(z_n-z_{n-1})}^{\alpha }}{\mu }{\left|u^n\right|}_1^2\right]}\hfill \\ \hfill \\ \le \hfill & {\displaystyle \sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)\left[\parallel u^k{\parallel }^2+\frac{p_0{(z_k-z_{k-1})}^{\alpha }}{\mu }{\left|u^k\right|}_1^2\right]+{\omega }_z^{n-1}E.}\hfill \end{array}$$

(19)

Now, we shall show by mathematical induction that

$${\parallel u^n\parallel }^2+\frac{p_0{(z_n-z_{n-1})}^{\alpha }}{\mu }{\left|u^n\right|}_1^2\le E,1\le n\le N.$$

(20)

In fact, let $n=1$ in (19) and we get

$${\omega }_z^0\left[\parallel u^1{\parallel }^2+\frac{p_0{(z_1-z_0)}^{\alpha }}{\mu }{\left|u^1\right|}_1^2\right]\le {\omega }_z^{0}E,$$

which implies that (20) holds for $n=1$. Suppose that (20) is true up to $(n-1)$. Then, from (19), we have

$$\begin{array}{cc}& {\displaystyle {\omega }_z^0\left[\parallel u^n{\parallel }^2+\frac{p_0{(z_n-z_{n-1})}^{\alpha }}{\mu }{\left|u^n\right|}_1^2\right]}\hfill \\ \hfill \\ \le \hfill & {\displaystyle \sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)\left[\parallel u^k{\parallel }^2+\frac{p_0{(z_k-z_{k-1})}^{\alpha }}{\mu }{\left|u^k\right|}_1^2\right]+{\omega }_z^{n-1}E}\hfill \\ \hfill \\ \le \hfill & {\displaystyle \sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)E+{\omega }_z^{n-1}E={\omega }_z^{0}E,}\hfill \end{array}$$

which immediately gives (20), or equivalently (15). This completes the proof. □

Remark 6

(Stability). Using a similar argument as in [4,5,6,7], it can readily be deduced from Theorem 2 that the numerical scheme (12) is robust (or stable) with respect to the initial data $\psi \left(x\right)$ and the non-homogeneous data $f(x,t).$

We are now ready to establish the convergence of the proposed scheme (12).

Theorem 3

(Convergence). Assume that $u(x,Z^{-1}\left(z\right))=\overline{u}(x,z)\in C^{4,2}([0,1]\times [z\left(0\right),z\left(1\right)])$. Suppose that the mesh ${\Delta }_z:z\left(0\right)=z_0<z_1<\cdots <z_{N-1}<z_N=z\left(1\right)$ is quasi uniform. Let $\{U_j^n=u(x_j,t_n)\}$ be the exact solution of the problem (1), $\left\{u_j^n\right\}$ be the numerical solution obtained from the scheme (12) and $e_j^n=U_j^n-u_j^n$ be the error at $(x_j,t_n)$. Then, we have

$$\parallel e^n{\parallel }^2+c(\alpha ,z){\left|e^n\right|}_1^2\le {\left[O({\tau }_{z,max}^{2-\alpha }+h^2)\right]}^2,1\le n\le N$$

(21)

where $c(\alpha ,z)=\frac{p_0}{{\omega }_z^0\mu }$.

Proof. 

Since $\left\{U_j^n\right\}$ is the exact solution of (1), it is clear that

$$\left\{\begin{array}{c}{\displaystyle \mu \left[{\omega }_z^0{U}_j^n-\sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)U_j^k-{\omega }_z^{n-1}{U}_j^0\right]=\mathrm{\Lambda }{U}_j^n+f_j^n+T_j^n,1\le j\le M-1}\hfill \\ \hfill \\ {\displaystyle U_j^0=\psi \left(x_j\right),1\le j\le M-1}\hfill \\ \hfill \\ U_0^n={\varphi }_1\left(t_n\right),U_M^n={\varphi }_1\left(t_n\right),0\le n\le N\hfill \end{array}\right.$$

(22)

where $T_j^n$ is the local truncation error of the j-th equation.

Noting (6) in Theorem 1, and using Taylor expansion at $x=x_j$ in (22), we find that

$$T_j^n=O({\tau }_{z,max}^{2-\alpha }+h^2),1\le j\le M-1,1\le n\le N.$$

(23)

Next, from (12) and (22) it is obvious that $\left\{e_j^n\right\}$ is the solution of the system

$$\left\{\begin{array}{c}{\displaystyle \mu \left[{\omega }_z^0{e}_j^n-\sum _{k=1}^{n-1}\left({\omega }_z^{n-k-1}-{\omega }_z^{n-k}\right)e_j^k-{\omega }_z^{n-1}{e}_j^0\right]=\mathrm{\Lambda }{e}_j^n+T_j^n,1\le j\le M-1}\hfill \\ \hfill \\ {\displaystyle e_j^0=0,1\le j\le M-1}\hfill \\ \hfill \\ e_0^n=e_M^n=0,0\le n\le N.\hfill \end{array}\right.$$

(24)

Hence, $e^n\in U_h$. Finally, applying Theorem 2 to (24) and noting $e_j^0=0$, we obtain for $1\le n\le N,$

$${\parallel e^n\parallel }^2+c(\alpha ,z){\left|e^n\right|}_1^2\le \frac{\Gamma (1-\alpha ){(z\left(1\right)-z\left(0\right))}^{\alpha }}{4p_0}\underset{1\le n\le N}{max}{\parallel T^n\parallel }^2={\left[O({\tau }_{z,max}^{2-\alpha }+h^2)\right]}^2$$

(25)

which completes the proof. □

Remark 7.

From Theorem 3, it is easily seen that

$$\parallel e^n\parallel =O({\tau }_{z,max}^{2-\alpha }+h^2).$$

Hence, the temporal convergence order in the norm $\parallel ·\parallel$ is $(2-\alpha )$, which is optimal. On the other hand, the convergence order in ${\parallel ·\parallel }_{\infty }$ is not optimal. In fact, from Lemma 2, (21) and the properties of quasi uniform mesh, we get

$$\begin{array}{cc}{\displaystyle \parallel e^n{\parallel }_{\infty }\le \frac{1}{2}{\left|e^n\right|}_1}\hfill & {\displaystyle \le \sqrt{\frac{{(z_n-z_{n-1})}^{-\alpha }}{4p_0\Gamma (2-\alpha )}}O({\tau }_{z,max}^{2-\alpha }+h^2)}\hfill \\ \hfill \\ & {\displaystyle \le \sqrt{\frac{{\tau }_{min,z}^{-\alpha }}{4p_0\Gamma (2-\alpha )}}O({\tau }_{z,max}^{2-\alpha }+h^2)=O({\tau }_{z,max}^{2-\frac{3}{2}\alpha }+h^2).}\hfill \end{array}$$

Remark 8.

Theorem 3 is an extension of the work [9] in the sense that

• the problem (1) involves an operator $\mathcal{L}$ and is more general than the problem considered in [9];
• we consider a general temporal mesh which is quasi uniform in terms of the scale function $z\left(t\right)$, in view of Remark 2 this is more general than the uniform temporal mesh considered in [9];
• (21) and Remark 7 give the explicit convergence order of the proposed scheme (12), while convergence is proven without giving the explicit convergence order in [9].

4. Numerical Simulation

In this section, we shall use two examples to demonstrate the efficacy of the proposed numerical scheme (12) and to verify the theoretical result in Remark 7. To be specific, we shall

• compute errors at $t=1$ using

  $$e(N,h)=\parallel e^N\parallel (or\parallel e^N{\parallel }_{\infty })$$

  as well as the corresponding temporal and spatial convergence orders by

  $$\mathrm{TCO}={log}_2\frac{e(N,h)}{e(2N,h)}\mathrm{and}\mathrm{SCO}={log}_2\frac{e(N,h)}{e(N,h/2)},$$

  respectively;
• demonstrate that the proposed numerical scheme (12) works well for three types of temporal meshes, namely

  –

  Uniform: uniform mesh with respect to t;

  –

  Graded: graded mesh [17,18] with $t_j=Z^{-1}\left((z\left(1\right)-z\left(0\right)){\left(\frac{j}{N}\right)}^r+z\left(0\right)\right),$$0\le j\le N$ (let $r=\frac{2-\alpha }{\alpha }$ in the experiment to get optimal accuracy, refer to [17,18] for details);

  –

  Uniform_z: uniform mesh with respect to $z\left(t\right)$.

  In view of Remark 2, we note that all the above types of temporal meshes are particular cases of quasi uniform mesh of $z\left(t\right)$.

Clearly, the exact solution is required to compute $e(N,h)$. When the exact solution is not available (which is commonly encountered in applications), we shall use ‘approximate’ exact solution, which is obtained by the numerical scheme (12) with sufficiently small mesh sizes (e.g., $M=N=2000$ in our experiments), as ‘exact’ solution to compute errors. This is reasonable as the numerical scheme (12) is convergent.

Example 1

([6,7,9]). Consider the generalized time-fractional diffusion equation

$$\left\{\begin{array}{c}{\displaystyle {}_0^C{D}_{t;\left[z\right(t),1]}^{\alpha }u(x,t)=u_{xx}(x,t)+f(x,t),(x,t)\in (0,1)\times (0,1)}\hfill \\ \hfill \\ {\displaystyle u(x,0)=sin\left(\pi x\right),x\in [0,1]}\hfill \\ \hfill \\ {\displaystyle u(0,t)=0,u(1,t)=0,t\in (0,1]}\hfill \end{array}\right.$$

(26)

where $0<\alpha <1$, $z\left(t\right)$ is a strictly increasing scale function and

$$f(x,t)=\frac{2}{\Gamma \left(2.15\right)}(x^2-x)t^{1.15}+{\pi }^{2}sin\left(\pi x\right)-2t^2.$$

Note that when $z\left(t\right)=t$ and $\alpha =0.85$, the exact solution of Equation (26) is

$$u(x,t)=sin\left(\pi x\right)+x(x-1)t^2.$$

In this example, $p(x,t)\equiv 1$ and $q(x,t)\equiv 0$. Let us start with $\alpha =0.85$ and $z\left(t\right)=t,t^{0.5},t^2$. Applying the numerical scheme (12) with fixed $h=\frac{1}{512}$ and varied N, we compute the errors and temporal convergence orders for three types of temporal meshes—Uniform, Graded, Uniform_z. The results are displayed in Table 1. From the table, it is easily seen that the experimental temporal convergence orders (≈1.15) are consistent with the theoretical ones (=$2-\alpha$) for various scale functions $z\left(t\right)$ and different types of temporal meshes. It is a pleasant surprise that the numerical performance in terms of maximum norm is better than the theoretical result in Remark 7.

Table 1. (Example 1) Temporal convergence order when $\alpha =0.85,h=\frac{1}{512}$.

$\mathit{z}\left(\mathit{t}\right)$	N	Uniform Mesh				Graded Mesh				Uniformz Mesh
		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$
t	16	8.55 $\times {10}^{-4}$	-	6.07 $\times {10}^{-4}$	-	1.15 $\times {10}^{-3}$	-	8.17 $\times {10}^{-4}$	-	-	-	-	-
	32	3.84 $\times {10}^{-4}$	1.15	2.73 $\times {10}^{-4}$	1.15	5.22 $\times {10}^{-4}$	1.14	3.71 $\times {10}^{-4}$	1.14	-	-	-	-
	64	1.71 $\times {10}^{-4}$	1.16	1.22 $\times {10}^{-4}$	1.16	2.34 $\times {10}^{-4}$	1.16	1.66 $\times {10}^{-4}$	1.16	-	-	-	-
	128	7.55 $\times {10}^{-5}$	1.18	5.36 $\times {10}^{-5}$	1.18	1.03 $\times {10}^{-4}$	1.18	7.33 $\times {10}^{-5}$	1.18	-	-	-	-
$t^{0.5}$	16	1.86 $\times {10}^{-3}$	-	1.32 $\times {10}^{-3}$	-	4.87 $\times {10}^{-3}$	-	3.46 $\times {10}^{-3}$	-	3.61 $\times {10}^{-3}$	-	2.56 $\times {10}^{-3}$	-
	32	8.38 $\times {10}^{-4}$	1.15	5.95 $\times {10}^{-4}$	1.15	2.27 $\times {10}^{-3}$	1.10	1.61 $\times {10}^{-3}$	1.10	1.66 $\times {10}^{-3}$	1.12	1.18 $\times {10}^{-3}$	1.12
	64	3.74 $\times {10}^{-4}$	1.16	2.66 $\times {10}^{-4}$	1.16	1.03 $\times {10}^{-3}$	1.13	7.35 $\times {10}^{-4}$	1.13	7.51 $\times {10}^{-4}$	1.14	5.34 $\times {10}^{-4}$	1.14
	128	1.65 $\times {10}^{-4}$	1.18	1.17 $\times {10}^{-4}$	1.18	4.60 $\times {10}^{-4}$	1.17	3.27 $\times {10}^{-4}$	1.17	3.33 $\times {10}^{-4}$	1.17	2.37 $\times {10}^{-4}$	1.17
$t^2$	16	2.66 $\times {10}^{-5}$	-	1.88 $\times {10}^{-5}$	-	1.50 $\times {10}^{-5}$	-	1.06 $\times {10}^{-5}$	-	5.68 $\times {10}^{-6}$	-	4.02 $\times {10}^{-6}$	-
	32	1.22 $\times {10}^{-5}$	1.12	8.66 $\times {10}^{-6}$	1.12	6.86 $\times {10}^{-6}$	1.12	4.86 $\times {10}^{-6}$	1.12	2.49 $\times {10}^{-6}$	1.19	1.76 $\times {10}^{-6}$	1.19
	64	5.52 $\times {10}^{-6}$	1.15	3.90 $\times {10}^{-6}$	1.15	3.08 $\times {10}^{-6}$	1.15	2.18 $\times {10}^{-6}$	1.15	1.07 $\times {10}^{-6}$	1.23	7.55 $\times {10}^{-7}$	1.23
	128	2.44 $\times {10}^{-6}$	1.18	1.73 $\times {10}^{-6}$	1.18	1.36 $\times {10}^{-6}$	1.18	9.63 $\times {10}^{-7}$	1.18	4.50 $\times {10}^{-7}$	1.24	3.19 $\times {10}^{-7}$	1.24

Table 1. (Example 1) Temporal convergence order when $\alpha =0.85,h=\frac{1}{512}$.

$\mathit{z}\left(\mathit{t}\right)$	N	Uniform Mesh				Graded Mesh				Uniformz Mesh
		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$
t	16	8.55 $\times {10}^{-4}$	-	6.07 $\times {10}^{-4}$	-	1.15 $\times {10}^{-3}$	-	8.17 $\times {10}^{-4}$	-	-	-	-	-
	32	3.84 $\times {10}^{-4}$	1.15	2.73 $\times {10}^{-4}$	1.15	5.22 $\times {10}^{-4}$	1.14	3.71 $\times {10}^{-4}$	1.14	-	-	-	-
	64	1.71 $\times {10}^{-4}$	1.16	1.22 $\times {10}^{-4}$	1.16	2.34 $\times {10}^{-4}$	1.16	1.66 $\times {10}^{-4}$	1.16	-	-	-	-
	128	7.55 $\times {10}^{-5}$	1.18	5.36 $\times {10}^{-5}$	1.18	1.03 $\times {10}^{-4}$	1.18	7.33 $\times {10}^{-5}$	1.18	-	-	-	-
$t^{0.5}$	16	1.86 $\times {10}^{-3}$	-	1.32 $\times {10}^{-3}$	-	4.87 $\times {10}^{-3}$	-	3.46 $\times {10}^{-3}$	-	3.61 $\times {10}^{-3}$	-	2.56 $\times {10}^{-3}$	-
	32	8.38 $\times {10}^{-4}$	1.15	5.95 $\times {10}^{-4}$	1.15	2.27 $\times {10}^{-3}$	1.10	1.61 $\times {10}^{-3}$	1.10	1.66 $\times {10}^{-3}$	1.12	1.18 $\times {10}^{-3}$	1.12
	64	3.74 $\times {10}^{-4}$	1.16	2.66 $\times {10}^{-4}$	1.16	1.03 $\times {10}^{-3}$	1.13	7.35 $\times {10}^{-4}$	1.13	7.51 $\times {10}^{-4}$	1.14	5.34 $\times {10}^{-4}$	1.14
	128	1.65 $\times {10}^{-4}$	1.18	1.17 $\times {10}^{-4}$	1.18	4.60 $\times {10}^{-4}$	1.17	3.27 $\times {10}^{-4}$	1.17	3.33 $\times {10}^{-4}$	1.17	2.37 $\times {10}^{-4}$	1.17
$t^2$	16	2.66 $\times {10}^{-5}$	-	1.88 $\times {10}^{-5}$	-	1.50 $\times {10}^{-5}$	-	1.06 $\times {10}^{-5}$	-	5.68 $\times {10}^{-6}$	-	4.02 $\times {10}^{-6}$	-
	32	1.22 $\times {10}^{-5}$	1.12	8.66 $\times {10}^{-6}$	1.12	6.86 $\times {10}^{-6}$	1.12	4.86 $\times {10}^{-6}$	1.12	2.49 $\times {10}^{-6}$	1.19	1.76 $\times {10}^{-6}$	1.19
	64	5.52 $\times {10}^{-6}$	1.15	3.90 $\times {10}^{-6}$	1.15	3.08 $\times {10}^{-6}$	1.15	2.18 $\times {10}^{-6}$	1.15	1.07 $\times {10}^{-6}$	1.23	7.55 $\times {10}^{-7}$	1.23
	128	2.44 $\times {10}^{-6}$	1.18	1.73 $\times {10}^{-6}$	1.18	1.36 $\times {10}^{-6}$	1.18	9.63 $\times {10}^{-7}$	1.18	4.50 $\times {10}^{-7}$	1.24	3.19 $\times {10}^{-7}$	1.24

Next, we investigate the performance of the proposed numerical scheme for different values of $\alpha$ and scale functions $z\left(t\right)$. Here, we use the uniform mesh of t (Uniform) and compute $\parallel e^N\parallel$.

• In Table 2, we fix $h=\frac{1}{512}$ and let N vary. It is observed that the temporal convergence order is $(2-\alpha )$ which is consistent with the theoretical one.
• In Table 3, we fix $N=2000$ and let h vary to investigate the spatial convergence. Obviously, the scheme (12) achieves second order spatial convergence order which is the same as that stated in Remark 7.

Table 2. (Example 1) Temporal convergence order for various $\alpha$ when $h=\frac{1}{512}$.

$\mathit{z}\left(\mathit{t}\right)$	N	$\mathit{\alpha }=0.2$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.8$
$t^{0.5}$	8	1.2119 $\times {10}^{-4}$	-	6.4821 $\times {10}^{-4}$	-	2.3862 $\times {10}^{-3}$	-
	16	3.6868 $\times {10}^{-5}$	1.72	2.3301 $\times {10}^{-4}$	1.48	1.0446 $\times {10}^{-3}$	1.19
	32	1.1124 $\times {10}^{-5}$	1.73	8.3262 $\times {10}^{-5}$	1.48	4.5475 $\times {10}^{-4}$	1.20
	64	3.3301 $\times {10}^{-6}$	1.74	2.9572 $\times {10}^{-5}$	1.49	1.9656 $\times {10}^{-4}$	1.21
$t^2$	8	3.6792 $\times {10}^{-6}$	-	7.9954 $\times {10}^{-6}$	-	2.7299 $\times {10}^{-5}$	-
	16	1.1025 $\times {10}^{-6}$	1.74	2.8681 $\times {10}^{-6}$	1.48	1.2814 $\times {10}^{-5}$	1.09
	32	3.2908 $\times {10}^{-7}$	1.74	1.0251 $\times {10}^{-6}$	1.48	5.7126 $\times {10}^{-6}$	1.17
	64	9.7656 $\times {10}^{-8}$	1.75	3.6452 $\times {10}^{-7}$	1.49	2.4954 $\times {10}^{-6}$	1.19

Table 2. (Example 1) Temporal convergence order for various $\alpha$ when $h=\frac{1}{512}$.

$\mathit{z}\left(\mathit{t}\right)$	N	$\mathit{\alpha }=0.2$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.8$
$t^{0.5}$	8	1.2119 $\times {10}^{-4}$	-	6.4821 $\times {10}^{-4}$	-	2.3862 $\times {10}^{-3}$	-
	16	3.6868 $\times {10}^{-5}$	1.72	2.3301 $\times {10}^{-4}$	1.48	1.0446 $\times {10}^{-3}$	1.19
	32	1.1124 $\times {10}^{-5}$	1.73	8.3262 $\times {10}^{-5}$	1.48	4.5475 $\times {10}^{-4}$	1.20
	64	3.3301 $\times {10}^{-6}$	1.74	2.9572 $\times {10}^{-5}$	1.49	1.9656 $\times {10}^{-4}$	1.21
$t^2$	8	3.6792 $\times {10}^{-6}$	-	7.9954 $\times {10}^{-6}$	-	2.7299 $\times {10}^{-5}$	-
	16	1.1025 $\times {10}^{-6}$	1.74	2.8681 $\times {10}^{-6}$	1.48	1.2814 $\times {10}^{-5}$	1.09
	32	3.2908 $\times {10}^{-7}$	1.74	1.0251 $\times {10}^{-6}$	1.48	5.7126 $\times {10}^{-6}$	1.17
	64	9.7656 $\times {10}^{-8}$	1.75	3.6452 $\times {10}^{-7}$	1.49	2.4954 $\times {10}^{-6}$	1.19

Table 3. (Example 1) Spatial convergence order for various $\alpha$ when $N=2000$.

$\mathit{z}\left(\mathit{t}\right)$	h	$\mathit{\alpha }=0.2$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.8$
$t^{0.5}$	$\frac{1}{10}$	5.3034 $\times {10}^{-3}$	-	5.5169 $\times {10}^{-3}$	-	5.7900 $\times {10}^{-3}$	-
	$\frac{1}{20}$	1.3214 $\times {10}^{-3}$	2.00	1.3746 $\times {10}^{-3}$	2.00	1.4425 $\times {10}^{-3}$	2.01
	$\frac{1}{40}$	3.2999 $\times {10}^{-4}$	2.00	3.4325 $\times {10}^{-4}$	2.00	3.6020 $\times {10}^{-4}$	2.00
	$\frac{1}{80}$	8.2382 $\times {10}^{-5}$	2.00	8.5692 $\times {10}^{-5}$	2.00	8.9924 $\times {10}^{-5}$	2.00
$t^2$	$\frac{1}{10}$	5.2701 $\times {10}^{-3}$	-	5.4177 $\times {10}^{-3}$	-	5.6052 $\times {10}^{-3}$	-
	$\frac{1}{20}$	1.3131 $\times {10}^{-3}$	2.00	1.3497 $\times {10}^{-3}$	2.01	1.3961 $\times {10}^{-3}$	2.01
	$\frac{1}{40}$	3.2791 $\times {10}^{-4}$	2.00	3.3704 $\times {10}^{-4}$	2.00	3.4862 $\times {10}^{-4}$	2.00
	$\frac{1}{80}$	8.1862 $\times {10}^{-5}$	2.00	8.4140 $\times {10}^{-5}$	2.00	8.7030 $\times {10}^{-5}$	2.00

Table 3. (Example 1) Spatial convergence order for various $\alpha$ when $N=2000$.

$\mathit{z}\left(\mathit{t}\right)$	h	$\mathit{\alpha }=0.2$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.8$
$t^{0.5}$	$\frac{1}{10}$	5.3034 $\times {10}^{-3}$	-	5.5169 $\times {10}^{-3}$	-	5.7900 $\times {10}^{-3}$	-
	$\frac{1}{20}$	1.3214 $\times {10}^{-3}$	2.00	1.3746 $\times {10}^{-3}$	2.00	1.4425 $\times {10}^{-3}$	2.01
	$\frac{1}{40}$	3.2999 $\times {10}^{-4}$	2.00	3.4325 $\times {10}^{-4}$	2.00	3.6020 $\times {10}^{-4}$	2.00
	$\frac{1}{80}$	8.2382 $\times {10}^{-5}$	2.00	8.5692 $\times {10}^{-5}$	2.00	8.9924 $\times {10}^{-5}$	2.00
$t^2$	$\frac{1}{10}$	5.2701 $\times {10}^{-3}$	-	5.4177 $\times {10}^{-3}$	-	5.6052 $\times {10}^{-3}$	-
	$\frac{1}{20}$	1.3131 $\times {10}^{-3}$	2.00	1.3497 $\times {10}^{-3}$	2.01	1.3961 $\times {10}^{-3}$	2.01
	$\frac{1}{40}$	3.2791 $\times {10}^{-4}$	2.00	3.3704 $\times {10}^{-4}$	2.00	3.4862 $\times {10}^{-4}$	2.00
	$\frac{1}{80}$	8.1862 $\times {10}^{-5}$	2.00	8.4140 $\times {10}^{-5}$	2.00	8.7030 $\times {10}^{-5}$	2.00

Our next example is modified from an example in [19] and it involves a general operator $\mathcal{L}$.

Example 2

([19]). Consider the generalized time-fractional diffusion equation

$$\left\{\begin{array}{c}{\displaystyle {}_0^C{D}_{t;\left[z\right(t),1]}^{\alpha }u(x,t)=\mathcal{L}u+f(x,t),(x,t)\in (0,1)\times (0,1)}\hfill \\ \hfill \\ {\displaystyle u(x,0)=sin\left(\pi x\right),x\in [0,1]}\hfill \\ \hfill \\ {\displaystyle u(0,t)=0,u(1,t)=0,t\in (0,1]}\hfill \end{array}\right.$$

(27)

where $\mathcal{L}u={\partial }_x\left(p(x,t){\partial }_{x}u\right)-q(x,t)u$,

$$p(x,t)=2-cos\left(xt\right),q(x,t)=1-sin\left(xt\right)$$

and

$$\begin{array}{cc}f(x,t)=\hfill & {\displaystyle \left[\frac{\Gamma (4+\alpha )}{6}{t}^3+\frac{2}{\Gamma (3-\alpha )}{t}^{2-\alpha }+\left(t^{3+\alpha }+t^2+1\right)\left({\pi }^{2}p(x,t)+q(x,t)\right)\right]sin\left(\pi x\right)}\hfill \\ \hfill \\ & {\displaystyle -\pi \left(t^{4+\alpha }+t^3+t\right)sin\left(xt\right)cos\left(\pi x\right).}\hfill \end{array}$$

When $z\left(t\right)=t$, the exact solution of (27) is

$$u(x,t)=sin\left(\pi x\right)(t^{3+\alpha }+t^2+1).$$

First, consider (27) when $\alpha =0.5$. With fixed $h=\frac{1}{2000}$, we shall apply the numerical scheme (12) and compute the errors $\parallel e^N\parallel$, $\parallel e^N{\parallel }_{\infty }$ and temporal convergence orders. The results are presented in Table 4, and it is clear that the numerical scheme (12) works well for different types of temporal meshes as well as for a wide range of problems (i.e., different $z\left(t\right)$). The experimental temporal convergence orders in $\parallel ·\parallel$ (≈1.5) agree with the theoretical ones (=$2-\alpha$), while once again it is a pleasant surprise that the numerical performance in terms of maximum norm is better than the theoretical result in Remark 7.

Table 4. (Example 2) Temporal convergence order when $\alpha =0.5,h=\frac{1}{2000}$.

$\mathit{z}\left(\mathit{t}\right)$	N	Uniform Mesh				Graded Mesh				Uniformz Mesh
		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$
t	40	7.02 $\times {10}^{-4}$	-	4.93 $\times {10}^{-4}$	-	3.13 $\times {10}^{-3}$	-	2.20 $\times {10}^{-3}$	-	-	-	-	-
	80	2.53 $\times {10}^{-4}$	1.47	1.78 $\times {10}^{-4}$	1.47	1.18 $\times {10}^{-3}$	1.41	8.26 $\times {10}^{-4}$	1.41	-	-	-	-
	160	8.94 $\times {10}^{-5}$	1.50	6.28 $\times {10}^{-5}$	1.50	4.27 $\times {10}^{-4}$	1.46	3.00 $\times {10}^{-4}$	1.46	-	-	-	-
	320	3.06 $\times {10}^{-5}$	1.55	2.15 $\times {10}^{-5}$	1.55	1.49 $\times {10}^{-4}$	1.52	1.04 $\times {10}^{-4}$	1.52	-	-	-	-
$t^{0.5}$	40	1.15 $\times {10}^{-3}$	-	8.05 $\times {10}^{-4}$	-	1.25 $\times {10}^{-2}$	-	8.80 $\times {10}^{-3}$	-	2.92 $\times {10}^{-3}$	-	2.05 $\times {10}^{-3}$	-
	80	4.12 $\times {10}^{-4}$	1.48	2.89 $\times {10}^{-4}$	1.48	4.91 $\times {10}^{-3}$	1.35	3.45 $\times {10}^{-3}$	1.35	1.07 $\times {10}^{-3}$	1.44	7.55 $\times {10}^{-4}$	1.44
	160	1.46 $\times {10}^{-4}$	1.50	1.02 $\times {10}^{-4}$	1.50	1.83 $\times {10}^{-3}$	1.42	1.29 $\times {10}^{-3}$	1.42	3.86 $\times {10}^{-4}$	1.48	2.71 $\times {10}^{-4}$	1.48
	320	4.97 $\times {10}^{-5}$	1.55	3.49 $\times {10}^{-5}$	1.55	6.47 $\times {10}^{-4}$	1.50	4.55 $\times {10}^{-4}$	1.50	1.33 $\times {10}^{-4}$	1.53	9.36 $\times {10}^{-5}$	1.53
$t^2$	40	2.56 $\times {10}^{-4}$	-	1.80 $\times {10}^{-4}$	-	4.48 $\times {10}^{-4}$	-	3.15 $\times {10}^{-4}$	-	1.03 $\times {10}^{-4}$	-	7.22 $\times {10}^{-5}$	-
	80	9.26 $\times {10}^{-5}$	1.47	6.50 $\times {10}^{-5}$	1.47	1.64 $\times {10}^{-4}$	1.45	1.15 $\times {10}^{-4}$	1.45	3.64 $\times {10}^{-5}$	1.50	2.55 $\times {10}^{-5}$	1.50
	160	3.28 $\times {10}^{-5}$	1.50	2.31 $\times {10}^{-5}$	1.50	5.86 $\times {10}^{-5}$	1.48	4.12 $\times {10}^{-5}$	1.48	1.27 $\times {10}^{-5}$	1.51	8.94 $\times {10}^{-6}$	1.51
	320	1.12 $\times {10}^{-5}$	1.55	7.90 $\times {10}^{-6}$	1.55	2.02 $\times {10}^{-5}$	1.54	1.42 $\times {10}^{-5}$	1.54	4.32 $\times {10}^{-6}$	1.56	3.04 $\times {10}^{-6}$	1.56

Table 4. (Example 2) Temporal convergence order when $\alpha =0.5,h=\frac{1}{2000}$.

$\mathit{z}\left(\mathit{t}\right)$	N	Uniform Mesh				Graded Mesh				Uniformz Mesh
		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$		$\parallel {\mathbf{e}}^{\mathbf{N}}{\parallel }_{\infty }$		$\parallel {\mathbf{e}}^{\mathbf{N}}\parallel$
t	40	7.02 $\times {10}^{-4}$	-	4.93 $\times {10}^{-4}$	-	3.13 $\times {10}^{-3}$	-	2.20 $\times {10}^{-3}$	-	-	-	-	-
	80	2.53 $\times {10}^{-4}$	1.47	1.78 $\times {10}^{-4}$	1.47	1.18 $\times {10}^{-3}$	1.41	8.26 $\times {10}^{-4}$	1.41	-	-	-	-
	160	8.94 $\times {10}^{-5}$	1.50	6.28 $\times {10}^{-5}$	1.50	4.27 $\times {10}^{-4}$	1.46	3.00 $\times {10}^{-4}$	1.46	-	-	-	-
	320	3.06 $\times {10}^{-5}$	1.55	2.15 $\times {10}^{-5}$	1.55	1.49 $\times {10}^{-4}$	1.52	1.04 $\times {10}^{-4}$	1.52	-	-	-	-
$t^{0.5}$	40	1.15 $\times {10}^{-3}$	-	8.05 $\times {10}^{-4}$	-	1.25 $\times {10}^{-2}$	-	8.80 $\times {10}^{-3}$	-	2.92 $\times {10}^{-3}$	-	2.05 $\times {10}^{-3}$	-
	80	4.12 $\times {10}^{-4}$	1.48	2.89 $\times {10}^{-4}$	1.48	4.91 $\times {10}^{-3}$	1.35	3.45 $\times {10}^{-3}$	1.35	1.07 $\times {10}^{-3}$	1.44	7.55 $\times {10}^{-4}$	1.44
	160	1.46 $\times {10}^{-4}$	1.50	1.02 $\times {10}^{-4}$	1.50	1.83 $\times {10}^{-3}$	1.42	1.29 $\times {10}^{-3}$	1.42	3.86 $\times {10}^{-4}$	1.48	2.71 $\times {10}^{-4}$	1.48
	320	4.97 $\times {10}^{-5}$	1.55	3.49 $\times {10}^{-5}$	1.55	6.47 $\times {10}^{-4}$	1.50	4.55 $\times {10}^{-4}$	1.50	1.33 $\times {10}^{-4}$	1.53	9.36 $\times {10}^{-5}$	1.53
$t^2$	40	2.56 $\times {10}^{-4}$	-	1.80 $\times {10}^{-4}$	-	4.48 $\times {10}^{-4}$	-	3.15 $\times {10}^{-4}$	-	1.03 $\times {10}^{-4}$	-	7.22 $\times {10}^{-5}$	-
	80	9.26 $\times {10}^{-5}$	1.47	6.50 $\times {10}^{-5}$	1.47	1.64 $\times {10}^{-4}$	1.45	1.15 $\times {10}^{-4}$	1.45	3.64 $\times {10}^{-5}$	1.50	2.55 $\times {10}^{-5}$	1.50
	160	3.28 $\times {10}^{-5}$	1.50	2.31 $\times {10}^{-5}$	1.50	5.86 $\times {10}^{-5}$	1.48	4.12 $\times {10}^{-5}$	1.48	1.27 $\times {10}^{-5}$	1.51	8.94 $\times {10}^{-6}$	1.51
	320	1.12 $\times {10}^{-5}$	1.55	7.90 $\times {10}^{-6}$	1.55	2.02 $\times {10}^{-5}$	1.54	1.42 $\times {10}^{-5}$	1.54	4.32 $\times {10}^{-6}$	1.56	3.04 $\times {10}^{-6}$	1.56

Next, we investigate the temporal convergence and spatial convergence of the numerical scheme (12) for different values of $\alpha$ and scale functions $z\left(t\right).$ Here, we use the uniform mesh of t and compute $\parallel e^N\parallel .$ The results are presented in

Table 5. (Example 2) Temporal convergence order for various $\alpha$ when $h=\frac{1}{2000}$.

$\mathit{z}\left(\mathit{t}\right)$	N	$\mathit{\alpha }=0.2$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.8$
$t^{0.5}$	10	8.7309 $\times {10}^{-4}$	-	5.9471 $\times {10}^{-3}$	-	2.5396 $\times {10}^{-2}$	-
	20	2.7332 $\times {10}^{-4}$	1.68	2.2093 $\times {10}^{-3}$	1.43	1.1397 $\times {10}^{-2}$	1.16
	40	8.3985 $\times {10}^{-5}$	1.70	8.0538 $\times {10}^{-4}$	1.46	5.0256 $\times {10}^{-3}$	1.18
	80	2.5426 $\times {10}^{-5}$	1.72	2.8938 $\times {10}^{-4}$	1.48	2.1831 $\times {10}^{-3}$	1.20
$t^2$	10	2.0452 $\times {10}^{-4}$	-	1.3063 $\times {10}^{-3}$	-	5.1492 $\times {10}^{-3}$	-
	20	6.5603 $\times {10}^{-5}$	1.64	4.9048 $\times {10}^{-4}$	1.41	2.3170 $\times {10}^{-3}$	1.15
	40	2.0515 $\times {10}^{-5}$	1.68	1.8006 $\times {10}^{-4}$	1.45	1.0227 $\times {10}^{-3}$	1.18
	80	6.2940 $\times {10}^{-6}$	1.70	6.5018 $\times {10}^{-5}$	1.47	4.4443 $\times {10}^{-4}$	1.20

and

Table 6. (Example 2) Spatial convergence order for various $\alpha$ when $N=2000$.

$\mathit{z}\left(\mathit{t}\right)$	h	$\mathit{\alpha }=0.2$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.8$
$t^{0.5}$	$\frac{1}{10}$	1.3899 $\times {10}^{-2}$	-	1.3094 $\times {10}^{-2}$	-	1.1890 $\times {10}^{-2}$	-
	$\frac{1}{20}$	3.4684 $\times {10}^{-3}$	2.00	3.2680 $\times {10}^{-3}$	2.00	2.9682 $\times {10}^{-3}$	2.00
	$\frac{1}{40}$	8.6645 $\times {10}^{-4}$	2.00	8.1642 $\times {10}^{-4}$	2.00	7.4155 $\times {10}^{-4}$	2.00
	$\frac{1}{80}$	2.1633 $\times {10}^{-4}$	2.00	2.0384 $\times {10}^{-4}$	2.00	1.8515 $\times {10}^{-4}$	2.00
$t^2$	$\frac{1}{10}$	1.4330 $\times {10}^{-2}$	-	1.4555 $\times {10}^{-2}$	-	1.4952 $\times {10}^{-2}$	-
	$\frac{1}{20}$	3.5756 $\times {10}^{-3}$	2.00	3.6316 $\times {10}^{-3}$	2.00	3.7304 $\times {10}^{-3}$	2.00
	$\frac{1}{40}$	8.9322 $\times {10}^{-4}$	2.00	9.0719 $\times {10}^{-4}$	2.00	9.3187 $\times {10}^{-4}$	2.00
	$\frac{1}{80}$	2.2301 $\times {10}^{-4}$	2.00	2.2650 $\times {10}^{-4}$	2.00	2.3266 $\times {10}^{-4}$	2.00

. We observe that the experimental temporal/spatial convergence orders are consistent with the theoretical result.

5. Conclusions

In this paper, we derive a numerical scheme based on a general temporal mesh for the generalized time-fractional diffusion problem. The main idea involves the generalized linear interpolation. The stability and convergence of the proposed numerical scheme is established rigorously using energy method. More importantly, it is shown that the global convergence order is $O({\tau }_{z,max}^{2-\alpha }+h^2)$ that extends the previous work [9]. For future work, we plan to investigate (i) the validity of the proposed scheme for nonsmooth data; (ii) high order methods based on general temporal mesh and spatial mesh for generalized fractional problems with smooth as well as nonsmooth data. We believe this will make the numerical scheme more applicable in reality.