An efficient high order numerical scheme for the time-fractional diffusion equation with uniform accuracy

Nomenclature

$u$: Change in temperature (K); $f$: The heat of source (K/s); $x$: Change in space (m); $t$: Change in time (s); $\Lambda$: Region; $T$: Maximum time; $\beta$: The order of fractional derivative of time; ${}_0D_t^{\beta }$: Fractional derivative of time in the sense of Caputo; $\partial_x^2$: The second derivative of space; $\Delta t$: The step of time (s); $h$: The space division step size (m); $K$: Maximum number of time divisions; $N$: Maximum number of space divisions

1.   Introduction

Fractional calculus is widely used in natural phenomena because many problems in engineering and science can be well described by fractional differential equation model. Recently, the research of solutions and numerical solutions of fractional differential equations is a very important research topic. For instance, finite difference method ^[1], spectral method ^[2,3], wavelets method ^[4], etc.

In this paper, we consider the TFDE as following form:

with the initial and boundary conditions as follows:

where $T > 0, 0 < \beta < 1, \Lambda = (a, b)$ and $I = (0, T]$. The Caputo fractional derivatives ${}_0D_t^{\beta } u(x, t)$ of $\beta$ order in (1.1) is defined by

where $\Gamma(\cdot)$ being Euler's gamma function. It is easy to see that the Caputo fractional derivative can be understood as the weighted integral average of the classical derivatives in the past time. This means that the change rate of function at the current time is affected by the past time. This property is useful to describe some materials mechanical behaviors. Therefore, the study of numerical algorithms for TFDE has attracted many researchers. In ^[5], an alternating direction implicit scheme was used to solve the TFDE of 2D with Dirichlet boundary condition with initial weak singularity solution by L1 scheme on uniform mesh. In ^[6], it proposed a novel stability and convergence numerical scheme for TFDE with $\alpha (1 < \alpha < 2)$ Caputo fractional derivative. The stability and convergence L2 type numerical scheme of the Caputo fractional derivative was constructed in ^[7]. A high-order convergence compact finite difference method for solving the TFDE of 2D was constructed in ^[8] by the L1 approximation with operator-splitting technique. In ^[9], they used L1 and L2 type approximation to construct a high order stability and convergence numerical method for the TFDE with rigorous analysis. In ^[10], it used Legendre polynomial to solve nonlinear fractional diffusion equation with advection and reaction terms. In ^[11], they constructed a space-time Petrov-Galerkin spectral method for TFDEs based on eigen-decomposition. In ^[12], it used the spectral collocation methods to solve distributed order TFDEs. In ^[13], they used the Legendre spectral tau method to solve the multi-term TFDEs and gave the error estimate and rigorous convergence analysis. For efficient numerical scheme constructing, an popular numerical method is based on sum-of-exponentials technique ^[14,15,16]. Many researchers use non-uniform grid or graded grid numerical schemes to solve problems when the solution of TFDE is initial value singularity ^{[17,18,19,20,21]}. In ^[22], they proposed an approximate spectral method for the nonlinear time-fractional partial integro-differential equation with a weakly singular kernel by using new basis functions based on shifted first-kind Chebyshev polynomials. From partial integro-differential equations, an efficient alternating direction implicit scheme were constructed for the nonlinear TFDEs in ^[23]. Based on the block-by-block method, a new general efficient technique to construct efficient numerical schemes for the fractional derivative was constructed in ^[24,25]. In ^[26], it presented finite difference and finite element scheme for space-time fractional diffusion equations. In ^[27], the time-stepping spectral method for the TFDEs was constructed by using the temporal second-order difference scheme and spatial spectral method discrete. In ^[28], it gave an efficient numerical scheme for the variable coefficient multiterm TFDE by the Crank-Nicolson method in temporal part and the exponential B-splines in spatial part with the unconditional stability and convergence rates analysis. Readers can refer to more references such as ^[29,30,31]. In this paper, we construct an efficient numerical scheme for TFDE with uniform accuracy by block-by-block approach to overcome the lack of the theoretical convergence order at the initial time of the above algorithms.

The organizational structure of the rest article is as follows: We first describe the time discretization for the time-fractional derivative and derive a sharp estimate for the truncation error in Section 2. In Section 3, the full discrete numerical scheme for TFDE is given by combing the spatial center second order scheme and the efficient time scheme of Section 2. In Section 4, the efficient scheme's stability and error estimates are strictly theoretical analysis to obtain that the scheme is unconditionally stable and $3-\beta$ convergence order in time with uniform accuracy. Some typical numerical examples are given to show the effectiveness of numerical algorithm in Section 5. We provide some concluding remarks in the final section.

2.   Discretization in time: a finite difference scheme

In order to simplify the symbols without losing generality, we set $f(x, t)\equiv0$ during the construction and analysis of numerical scheme. Now, we construct an efficient high order scheme for the fractional derivative in time, and divide the interval $(0, T]$ into $K$ equal subintervals with $\Delta t = \frac{T}{K}, t_k = k\Delta t, k = 0, 1, \cdots, K$.

Next, we discuss the numerical scheme for time-fractional derivative ${}_0D_t^{\beta } u(x, t)$ in (1.1). Firstly, we determine the values of $u(\cdot, t)$ on $t_1$ and $t_2$. The approximate formula of $u(\cdot, t)$ on the interval $[t_0, t_2]$ is given by Lagrange interpolation ^[32]:

where $\bar{\omega }_{i.0}(t)$ is defined by

When $k = 1, 2$, based on (2.1) we can approximate ${}_0D_t^\beta u(x, t_1), {}_0D_t^\beta u(x, t_2)$ by replacing $u(\cdot, t)$ with $J_{\left[t_0, t_2\right]} u(\cdot, t)$ as following,

where

Divide the integration interval into several subintervals, one can obtain ${}_0D_t^\beta u(x, t_k)$ for $k \geq 3$ as follows:

For every subintervals $[t_j, t_{j+1}], j = 1, 2, \cdots, k-1$, the approximation of $u(x, t)$ is as follows:

where $w_{i, j}, i = 0, 1, 2;j = 1, 2, \cdots k-1$ are defined by

Similar to (2.2), based on (2.3) for $k\geq 3$, $_0D_t^\beta u(x, t_k)$ can be approximated by replacing $u(\cdot, t)$ with $J_{\left[t_j, t_{j+1}\right]} u(\cdot, t)$ in every $[t_j, t_{j+1}]$ as follows:

where

Combining (2.2), (2.4), let $\beta_0 = \Gamma(3-\beta) \Delta t^\beta$, one can immediately obtain the efficient scheme for $_0D_t^\beta u(x, t_k)$ as follows:

where

In the following Theorem 2.1, we give the error estimate to proposed numerical scheme (2.5) of fractional derivative.

Theorem 2.1. Suppose $u(\cdot, t) \in C^{3}[0, T]$ and denote

then

where $C_{u}$ is a constant and only depending on the function $u$.

Proof. Based on the Taylor theorem, one can obtain the residue of (2.1) as follows

where $\xi(t) \in [t_{0}, t_{2}]$.

Firstly, we will use (2.7) to estimate the case $k = 1$ of (2.6) as follows:

Secondly, the estimate of the case $k = 2$ of (2.6) is similar to $k = 1$, so the proof process is omitted without losing generality.

Thirdly, based on the direct calculation it is easy to split the estimate of the case of (2.6) for $k\geq 3$ into two parts in the following

Repeating the similar proof process of $|r_1(\Delta t)|$, one can immediately obtain the following estimate

For the estimate of the second part $N$ of $|r_k(\Delta t)|$, using the Taylor theorem one can obtain

In the above proof, we use $|(\tau-t_{j})(\tau-t_{j+1})(\tau-t_{j-1})|\leq\frac{2\sqrt{3}}{9}\Delta t^{3}.$ Combining (2.8) and (2.9), one can get the following conclusions

To sum up, the proof of Theorem 2.1 is completed. □

3.   Full discretization

For convenience and generality, we divide $\Lambda = (a, b)$ into $N$ equal parts and denote $h = \frac{b-a}{N}$, and $x_j = a+jh, j = 0, 1, \cdots, N$. The numerical solution (1.1) at $(x_j, t_k)$ is denoted as $u_j^k$. We use the following differential notations

Similar to ^[16], for the sake stability and convergence we introduce the discrete inner product and norm as follows,

Though direct calculation, it is easy to prove the discrete Green formula:

Let $u = y, v = z_{\bar{x}}$, the Eq (3.3) immediately becomes

Therefore, using (3.1), the second-order central difference scheme in space as follows:

Combining (2.5) with (3.5), one can immediately obtain the full discrete scheme for TFDEs as follows:

In order to analyze the numerical scheme (3.6), we firstly rewrite (3.6) for $k\geq4$ into the following equivalent form

where

Similarly, we secondly rewrite (3.6) for $k = 3$ into the following equivalent form

where

According to (3.7) and (3.9), the full discrete scheme of (3.6) can be rewritten as the following equivalent form

Before carrying out the full discrete scheme (4.1)'s stability and convergence, the properties of the coefficients $d_{k-i}^k$ will be analyze firstly in the Lemma 3.1.

Lemma 3.1. For all $0 < \beta < 1$, as $k \geq 4$, the coefficients in the scheme (4.1d) satisfy

(I) $C_1 = \frac{4-\beta}{2} \in(\frac{3}{2}, 2)$,

(II) ${ }\sum_{i = 1}^k d_{k -i}^k\equiv1$,

(III) $d_{k-i}^k > 0, i = 3, \cdots, k$,

(IV) $0 < d_{k-1}^k < \frac{4}{3}$,

(V) $d_{k-2}^k$ there are positive and negative,

(VI) $d_{k-2}^k+\frac{1}{4}(d_{k-1}^k)^2 > 0.$

Proof. (I) Based on the direct calculation, it is easy to prove directly.

(II) Based on the fact that the fractional derivative of the constant is zero, one can prove it immediately by combining (2.5) and (3.8).

(III) For $k\geq4$, we can easy to observe that

Denote $i-2 = \bar{i}$ for $i\geq 6$, and we use Taylor formula yields

where

It is easy to check that $\sum_{k = 0}^{+\infty}a_{k}$ is an alternating series with a positive first term by carefully calculate, and we have $0 < \sum_{k = 0}^{+\infty}a_{k} < 4(\beta+1)$. Therefore, we have

Similarly, based on the directly calculating it is easy to prove as follows:

Combining (3.11) and (3.12), one can easily obtain that

For $d_{1}^{k} = -(A_{k-2}+B_{k-1}+\bar{B}_{k})\beta_{0}^{-1}$, we let $W_{1} = -(A_{k-2}+B_{k-1}+\bar{B}_{k})$, then

Let $k-2 = \bar{k}$, we use a Taylor expansion and get following as

where

Therefore, ${ } \sum\limits_{n = 0}^{+\infty}\prod\limits_{i = 0}^{n}(-i-\beta)(\frac{1}{\bar{k}})^{n}\cdot a_{n}$ is an alternating series with a positive first term and satisfy

For $d_{2}^{k} = -(\bar{C}_{k}+A_{k-3}+B_{k-2}+C_{k-1})\beta_{0}^{-1}$, taking $W_{2} = -(A_{k-3}+B_{k-2}+C_{k-1}+\bar{C}_{k})$, we have

Let $k-2 = \hat{k}$, we still use Taylor formula and get

where

As $n\geq 2$, $b_{n} < 0$ for odd number $n$ and $b_{n} = (n-1)[4\cdot2^{n}-\frac{3}{2}]-3 > 0$ for even number $n$. It is easy to verify that $b_{0}$ and $b_{1}$ are all positive. Therefore, it is an alternating series for $n\geq 2$, i.e.,

Therefore, we have

(IV) According to (3.8), we have

Therefore,

Let $f(\beta) = 5(4-\beta)+3(\beta-6)2^{1-\beta}$, by carefully calculate, we have

Therefore, $f^{\prime}(\beta)$ is a monotone decreasing function and $f^{\prime}(1) = 15\ln2-2 > 0, 0 < f^{\prime}(1) < f^{\prime}(\beta) < f^{\prime}(0)$. So $f(\beta)$ is a monotone increasing function and $f(0) < f(\beta) < f(1) = 0$. To sum up, we obtain $\frac{4}{3} > d_{k-1}^{k} > 0$.

We take $g(\beta) = 3(4-\beta)+(\beta-6)2^{1-\beta}$, due to

We can deduce $g(\beta)$ is a monotone increasing function and $0 < g(0) < g(\beta) < g(1)$. Therefore, $0 < d_{k-1}^{k} < \frac{4}{3}.$

(V) As $k\geq4$, we have

Next, we discuss the sign of the $f(\beta)$. $\forall {\beta} \in (0, 1)$, $4-\beta > 0$ always holds. So the sign of the $d_{k-2}^{k}$ is determined by $f(\beta)$. After calculation, $f^{\prime}(0) > 0, f^{\prime}(1) < 0$. Therefore, $f(\beta)$ increases first and then decreases and $f(0) = 0, f(1) = -1 < 0$, we know $f(\beta)$ has only one zero $\beta_{0}^{*}\in(0, 1)$. Furthermore, when $\beta\in(0, \beta_{0}^{*})$, $f(\beta) > 0$. When $\beta\in(\beta_{0}^{*}, 1)$, $f(\beta) < 0$. That is, $d_{k-2}^{k}$ has positive and negative on $\beta \in(0, 1)$.

(VI) The following equation can be obtained through direct calculation

where

According to Lagrange mean value theorem, we have $4^{1-\beta} > \frac{4}{3+\beta}\cdot 3^{1-\beta}$, and

where

So $\sum_{k = 0}^{+\infty}b_{k}$ is an alternating and convergent series, i.e., $0 < \sum_{k = 0}^{+\infty}b_{k} < b_{0}$, we have

where

It can be obtained by careful calculation, we obtain

Next, we estimate $\bar{f}_{3}(\beta)$. $\forall\beta\in(0, 1)$, $\beta^{3} < \beta^{2}, \beta^{6} > 0$, we have

where

After careful calculation, $f^{\prime\prime}_{3}(\beta) < 0$. Therefore, $f^{\prime}_{3}(\beta)$ is a monotone decreasing function, so $f^{\prime}_{3}(0) > f^{\prime}_{3}(\beta) > f^{\prime}_{3}(1)$. Similarly,

we find the sign of $f^{\prime}_{3}(\beta)$ become positive to negative, i.e., $f_{3}(\beta)$ increases first and then decreases. $f_{3}(0) = 0$, $f_{3}(1) = 191 > 0$. Therefore we prove

Combining (3.14), (3.15) with (3.16), (3.17), we have $f(\beta) > f_{1}(\beta) > f_{3}(\beta) > 0$. From (3.13), we find

To sum up, we already proved Lemma 3.1. □

Remark 3.1. From Lemma 3.1, the coefficient $d_{k-2}^{k}$ has positive and negative on $\beta \in(0, 1)$: When $\beta\in(0, \beta_{0}^{*})$, $d_{k-2}^{k} > 0$; when $\beta\in(\beta_{0}^{*}, 1)$, $d_{k-2}^{k} < 0$. This brings great difficulties to the convergence and stability analysis of the proposed scheme, and the classical theoretical analysis method are invalid. Therefore, here a novel technique for the stability analysis $\forall\beta \in(0, 1)$ will be given.

Next, we will introduce a technique for numerical scheme (4.1):

By recombining of the Eq (4.1d), we have

Next, we denote

Thus (4.1d) can be equivalent as follows:

In the following Lemma 3.2, we give some good properties of the coefficients in numerical scheme (4.1).

Lemma 3.2. When $k \geq 4$, $\forall0 < \beta < 1$, the coefficients in numerical scheme (3.21) satisfy

(I) $0 < \rho < \frac{2}{3}$,

(II) $\bar{d}_{k- i}^k > 0, i = 2, 3, \cdots, k$,

(III) ${ }\rho+\sum_{i = 2}^{k-1} \bar{d}_{k-i}^k+\bar{d}_0^k \leq 1.$

Proof. (I) According to (3.18) and (IV) in Lemma 3.1, we immediately give the estimate for $\rho$.

(II) As $i = 2$, we have

According to (VI) in Lemma 3.1, we get $\bar d_{k-2}^{k} > 0$. Furthermore,

Because of $\rho > 0$ and (III) in Lemma 3.1, we obtain

(III) Let $P_{k} = { }\rho + \sum_{i = 2}^{k-1}\bar d_{k-i}^{k}+\bar d_{0}^{k}$, we can obtain from (3.19) that

That is,

Then, we use (II), (III), (VI) in Lemma 3.1 gives

Therefore, we already proved ${ }\rho+\sum\limits_{i = 2}^{k-1} \bar{d}_{k-i}^k+\bar{d}_0^k \leq 1$. □

As $k = 3$, the equivalent of (4.1c) as follows:

where

We find $\rho\neq d_2^3-\rho$. So we will give some properties of the coefficients in numerical scheme (3.22) in the following Lemma 3.3 as $k = 3$.

Lemma 3.3. The coefficients of (3.22) have the following properties when $0 < \beta < 1$ and $k = 3$,

(I) $\bar{d}_{3-i}^3 > 0, i = 1, 2, 3$,

(II) $\bar{d}_{2}^{3}-\rho < 0$,

(III) $\bar{d}_{0}^{3}+\bar{d}_{1}^{3}+\bar{d}_{2}^{3}\leq 1$.

Proof. (I) Let's firstly prove that,

Based on calculating carefully, one can immediately obtain that

Because of $d_{1}^{3} = \frac{2}{4-\beta}[(2\beta-2)3^{1-\beta}-6+\frac{3}{2}\beta]$, so

where

Next, using a Taylor expansion yields

Next, we will estimate $a_{1}\cdot2^{1-\beta}+a_{2}\cdot2^{1-2\beta} > 0$.

where $f(\beta) = (5\beta-4)(4-\beta)+(4+\beta)(6-\beta)$. Because of $f'(\beta) = 26-12\beta > 0, f(\beta)$ monotonic increase, $f(\beta)\geq f(0) = 8 > 0$, so $f(\beta) > 0$.

Because of $\frac{1-\beta}{1!}\cdot\frac{1}{2}+\frac{(1-\beta)(-\beta)}{2!}\cdot(\frac{1}{2})^{2} +\frac{(1-\beta)(-\beta)(-\beta-1)}{3!}(\frac{1}{2})^{3}+\cdots\doteq\sum_{k = 0}^{+\infty}a_{k}$ is an alternating series with positive first term, and $\sum_{k = 0}^{+\infty}a_{k} = a_{0}+a_{1}+\sum_{k = 2}^{+\infty}a_{k}$, where $\sum_{k = 2}^{+\infty}a_{k}$ is an alternating series which the first term is positive number, so $0 < \sum_{k = 2}^{+\infty}a_{k} < a_{2}$, and

where

Next, we will prove $\bar{d}_{1}^{3} > 0$, that is $-\frac{3}{4}+\frac{2^{-\beta}}{8(4-\beta)^{2}}[f_{1}(\beta)+f_{2}(\beta)] > 0\Longleftrightarrow f_{1}(\beta)+f_{2}(\beta) > \frac{3}{4}\cdot8(4-\beta)^{2}2^{\beta} = 6(4-\beta)^{2}2^{\beta}$, that is $f_{1}(\beta)+f_{2}(\beta) > 6(4-\beta)^{2}2^{\beta}\doteq6f_{3}(\beta)$. So to prove $\bar{D}_{1}^{3} > 0$, just prove $f_{1}(\beta)+f_{2}(\beta)-6f_{3}(\beta) > 0$. Let's remember $\bar{f}(\beta) = f_{1}(\beta)+f_{2}(\beta)-6f_{3}(\beta)$. Since $\bar{f}(\beta)$ is an increasing firstly and then decreasing function, and $\bar{f}(0) = 0$, $\bar{f}(1) = 16 > 0$, so $\bar{f}(\beta) > 0$, therefore $\bar{d}_{1}^{3} > 0$.

Because $d_{0}^{3} = [2-(3^{2-\beta}+1)\frac{\beta}{2}] > 0$, so $\bar{d}_{0}^{3} = \bar{d}_{1}^{3}\rho+d_{0}^{3} > 0.$ To sum up, (3.24) is completed proved.

(II) Because of

That is $\bar{d}_{2}^{3}-\rho < 0$.

(III) According to (3.23), we have

Therefore, we have

where $d_{0}^{3} > 0$, $\bar{d}_{1}^{3} > 0$. By carefully calculate, we have $d_{0}^{3}+d_{1}^{3}+d_{2}^{3} = 1$. According to (3.25), we obtain $(1-\rho)P_{3}\leq 1-\rho$, i.e.,

The proof is then completed. □

4.   Stability and error estimates

From (3.21) and (3.22), thus the equivalent of (4.1) as follows:

We will give the estimation of $\|\bar{u}^{1}\|_0^2+\beta_0C_1^{-1}\|u_{\bar{x}}^1]|^2, \|\bar{u}^{2}\|_0^2+\beta_0C_1^{-1}\|u_{\bar{x}}^2]|^2$ as following Lemma 4.1,

Lemma 4.1. Let $\widehat{\beta} = \min\{-\widehat{D_1}\widetilde{D_1}, \widetilde{D_2}\widehat{D_2}, -\widetilde{D_1}, \widehat{D_2}\}$ and $\widehat{\alpha} = \max\{\widehat{D_0}\widetilde{D_1}, |-\widetilde{D_0}\widehat{D_2}|\}$, we have

where $M_0$ satisfies

Proof. Multiplying $-\widetilde{D_1}u_j^1h$ on both sides of (4.1a) for $k = 1$ and taking the sum over $j$, one immediately gets

Multiplying $\widehat{D_2}u_j^2h$ on both sides of (4.1b) by for $k = 2$ and taking the sum over $j$ with (3.2), one immediately gets

Add Eqs (4.5) and (4.6) correspondingly, one can obtain by combining similar terms,

Because of $-\widehat{D_1}\widetilde{D_1}, \widetilde{D_2}\widehat{D_2}, -\widetilde{D_1}, \widehat{D_2}$ and $\widehat{D_0}\widetilde{D_1}$ are all positive number depend on $\beta$, therefore $\widehat{\beta} > 0, \widehat{\alpha} > 0$, we have

That is,

Because of $C_1 = \frac{4-\beta}{2}$, therefore $C_1^{-1} = \frac{2}{4-\beta}\in (\frac{1}{2}, \frac{2}{3})$, so $C_1^{-1} < 1, \beta_0C_1^{-1} < \beta_0$. From (4.7), we can get

According to $\overline{u}^{1} = u^{1}-\rho u^0$ and trigonometric inequality, one can get

Using (4.8) and (4.10), it is easy to obtain that

Similarly, we will estimate $\|\bar{u}^{2}\|_0^2+\beta_0C_1^{-1}\|u_{\bar{x}}^2]|^2$. Using (4.8) and (4.9), one can get

In summary, by using (4.4), combining (4.11) with (4.12), then

Lemma 4.1 already completed proved. □

Next, for $k\geq3$, we give the estimate for $\|\bar{u}^k\|_0^2+\beta_0 C_1^{-1} \| u_{\bar x}^k]|^2$ in the following Lemma 4.2, which is an important conclusion for the stability analysis of the proposed scheme.

Lemma 4.2. $\|\bar{u}^k\|_0^2+\beta_0 C_1^{-1} \| u_{\bar x}^k]|^2\leq M_0\|u^0\|_0^2, \ \ 3\leq k\leq K,$ where $M_0$ defined in (4.4).

Proof. First, we deduce the general formula. Using (1.3) and (3.4), we have

That is,

Multiplying $2\bar{u}_j^3h$ on both sides of (4.1c) for $k = 3$ and taking the sum over $j$ and using (4.13), we get

According to (I) in Lemma 3.2 and (II) in lemma 3.3, we have

By directly computing, it can be easy to deduce that $\bar{d}_{0}^{3}+\bar{d}_{1}^{3}+\rho\leq 1$. By using (4.2) and (4.3), (4.14) is becoming as follows

For $k\geq4$, multiplying both sides of the (4.1d) by $2\bar{u}_j^kh$ and taking the sum over $j$ and using (4.13), one gets

According to (III) in Lemma 3.3, the above inequality (4.16) becomes

By (I) in Lemma 3.2: $0 < \rho < 1$ and (4.17) yields

That is,

Using the mathematics induction, it is easy to prove the following inequality

As $k = 4$, by (4.2), (4.3) and (4.15), from (4.18) we can obtain

According to (4.20), this is the case of (4.19) when $k = 4$. Assuming that (19) establish for $k = 5, 6, \cdots, K-1$, and from (4.18) one immediately obtain that,

The proof of Lemma 4.2 is completed. □

We will give numerical scheme (4.1) is unconditionally stable in the following Theorem 4.1.

Theorem 4.1. The numerical scheme (4.1) is unconditionally stable, and its solution satisfies the following estimate for all $\Delta t > 0, h > 0$:

Proof. According to Lemma 4.1 and Lemma 4.2, we immediately give the following estimate

From (4.21), we have

According to (3.20), (I) in Lemma 3.2, we have

According to the above estimation, we have

Theorem 4.1 is then completed. □

In the next Theorem 4.2, we give the convergence analysis of the full discrete scheme.

Theorem 4.2. Assume that $u(x, t)$ which is the solution of Eq (1.1), $u_j^k$ be the solution of the problem (4.1), Suppose ${ } \max\limits_{t\in (0, T]}|\partial_t^3u| \leq M$. Then for $k = 1, 2, \cdots, K$, we have

where $0 < \beta < 1,$ and $C$ is a positive constant that does not depend on $\Delta t$, $h$.

Proof. According to Theorem 2.1 and (3.5), similar to the Theorem 4.1, we can obtain the proof of Theorem 4.2. Here we omit it. □

Remark 4.1 The convergence of the proposed scheme can be analysis with the above analysis technique and the idea of ^[21] on the graded mesh. For the stability analysis of the proposed scheme is very difficult at present by using the analysis method in this paper. The mainly difficulty is whether the coefficients $d^k_{k-i}$ in Lemma 3.1 is nonnegativity for the graded temporal mesh, which is still an open problem up to now.

5.   Numerical validation

5.1.   Algorithm implementation

For the sake of simplicity, we only implement algorithm (4.1) as following:

$(1)$ Denote the $I_{(N-1)\times(N-1)}$ as a identity matrix of $(N-1)\times(N-1)$, and $D_{(N-1)\times(N-1)}$ as a difference matrix of second derivative in space, i.e.,

$(2)$ For $k = 1$ and $k = 2$, based on (4.1a) and (4.1b), we have

where $(u^1)^{\mathit{T}}$ represent the transpose of $(u^1)$. Solving the above Eq (5.1), we obtain $u^1$ and $u^2$.

$(3)$ For $k = 3$, based on (4.1c), we have

By solving the above Eq (5.2), we can obtain $u^3$. {1.} $(4)$ For $k\geq4$, based on (4.1d), we have

We solve the Eq (5.3) and get $u^k, k = 4, 5, \cdots, k$.

5.2.   Numerical results

In this part, we carry out in this section a series of numerical experiments and present some results to confirm our theoretical statements. The main purpose is to check the convergence behavior of the numerical solution with respect to the time step $\Delta t$ and space size $h$ used in the calculation.The first two numerical examples are proposed to show the efficiency of the $3-\beta$ order in time and second order in space of one and two dimension in space, respectively. The last numerical example, we choose the graded grid numerical schemes to solve problems when the solution of TFDEs is initial value singularity.

Example 5.1. Consider $f(x, t)$ and $u_0(x)$ in (1.2) as the following form

it is easy to check that the exact solution is given by $u(x, t) = t^4\sin(2\pi x)$. In this example, we choose $a = 0, b = 1, T = 1, \beta = 0.2, 0.5, 0.8$ and the step size $\Delta t = \frac{1}{K}, h = \frac{1}{N}$. Here we take two different sets of step size parameters. In the first case, we verify that the spatial convergence order and choose $K = 2^{10}$, and $N = 2^k, k = 2, 3, 4, 5, 6, 7$. In the second case, we verify that the convergence order in time and choose $K = 2^k, k = 2, 3, 4, 5, 6, 7$; and $N = [K^{\frac{3-\beta}{2}}]$ where $[x]$ denote the maxinum integer part of $x$. We denote the max error as $e^{\Delta t, h} = \mathop {\max }\limits_{j, k} | {u_j^k - u({x_j}, {t_k})}|$, here $u_j^k$ and $u({x_j}, {t_k})$ are the numerical solution and exact solution for (1.1) at point $(x_j, t_k)$, respectively.

Firstly, we plot the error distribution. In Figure 1, the error distribution of $K = 2^6, N = [K^{(1.5-0.5*\beta)}]$ and $\beta = 0.5$ is shown, where $[\cdot]$ indicates rounding up. From Figure 1, we find that the errors can be as small as $10^{-4}$.

Secondly, we plot log-log graph of error. In Figure 2, A logarithmic scale has been used for both $\Delta t$-axis and error-axis in this figure. For $\beta = 0.2, 0.8$, we find the temporal approximation order close to $3-\beta$, i.e. the slopes of the error curves in these log-log plots are 2.8, 2.2 respectively for $\beta = 0.2, 0.8$.

It is easy to see that the convergence order in Table 1 is very close to $2$. This shows that the convergence order of space is 2 and not related to the selection of $\beta$, which is consistent with the theoretical result of the Theorem 4.2.

From Table 2, it is easy to see that convergence order is almost $2.8, 2.5, 2.2$ with respect to $\beta = 0.2, 0.5, 0.8$, respectively. This indicates that the time convergence order is $3-\beta$ which is consistent with the theoretical result of the Theorem 4.2.

Example 5.2. Consider $f(x, y, t)$ and $u_0(x, y)$ in (1.1) as the following form

it is easy to check that the exact solution of $(1.1)$ is following

In this numerical example, we choose $T = 1$ and space domain is $[0, 1]^2$. Denote the step size $\Delta t = \frac{1}{K}$, and spatial size $\Delta x = \Delta y = h = \frac{1}{N}$, and the error is

Here we take two different sets of step size parameters. In the first case, we verify that the spatial convergence order and choose $K = 5\cdot 2^8$ and $N = 5\cdot 2^k, k = 2, 3, 4, 5, 6$. In the second case, we verify that the temporal convergence order and choose $K = 5\cdot 2^k, k = 2, 3, 4, 5, 6$, and $N = [K^{\frac{3-\beta}{2}}]$.

From Table 3, it is easy to show that the convergence order in space is close to $2$ which is consistent with the convergence order analysis in space of the Theorem 4.2.

From Table 4, one can be easy to see that convergence order is almost $2.8, 2.5, 2.2$ for $\beta = 0.2, 0.5, 0.8$, respectively. This indicates that the temporal convergence order is $3-\beta$ which is consistent with the theoretical result of the Theorem 4.2.

In the above two numerical examples, we assume that the solution is sufficiently smooth, which is also the basic requirement of all high order numerical schemes. In order to illustrate the effectiveness of the numerical scheme of this paper, the third numerical example is the numerical result of nonsmooth solution by using the graded mesh.

Example 5.3. In this example, we consider the problem (1.1) with exact solution is

which the third derivative of time is singular at the initial value. It is easy to check that the right hand function is

and $u_0(x) = 0.$

It is easy to prove that the exact solution of this example is not satisfied the condition of Theorem 2.1. In order to ensure the time convergence order unreduced, we choose graded mesh is $t_j = (j/K)^{\gamma}, j = 0, 1, \cdots, K$, with a graded parameter $\gamma$ base on the idea of ^[21] and $\beta = 0.3, 0.5, 0.7$. Here we take two different sets of step size parameters. In the first case, we verify that the spatial convergence order and choose $K = 2^9$ and $N = 2^k, k = 2, 3, 4, 5, 6, 7$. In the second case, we verify that the temporal convergence order and choose $K = 4, 8, 16, 32, 64, 128$ with $\gamma = \frac{3-\beta}{\beta}$, $N = [K^{\frac{3-\beta}{2}}]$.

From Table 5, one can see that the convergence order of spatial is almost $2$. And under the condition of $\gamma = \frac{3-\beta}{\beta}$, it is easy to obtain that the convergence order in time is $3-\beta$ with respect to the theoretical analysis Theorem 4.1 in ^[21]. From the Table 6, we see that the convergence order in time is almost $3-\beta$.

6.   Conclusions

On the idea of ^[24,25], we propose a uniform accuracy $3-\beta$ order for fractional derivatives based on piecewise quadratic interpolation, and apply it to solve TFDE. In particular, the proposed numerical scheme overcome the problem of order reduction at the initial value of the high-order numerical scheme, and also provides a general construction method for the numerical scheme with uniform convergence order. The stability and error estimates of high order uniform accuracy are strictly theoretical analysis. In terms of numerical implementation, we firstly use the full discrete scheme to solve one and two dimensional TFDE with the sufficiently smooth solution. We secondly use the numerical scheme to solve the TFDE with nonsmooth solution on the graded mesh. These two kinds of numerical examples fully illustrate that the full discrete scheme of this paper is very effective. The method of analyzing the convergence and stability of schemes in this paper can provide an effective analysis tool for analyzing the convergence and stability of numerical schemes of fractional integro-differential equations. In the future, we will use the proposed effective scheme to solve TFDE optimal control problem, and expands proposed effective scheme for solving the nonlinear multi-dimensional fractional integro-differential equations. Based on the ideas of references ^[33], we also consider to solve the nonlinear Lane-Emden equation with fractal-fractional derivative.

Acknowledgments

This work was supported by National Natural Science Foundation of China (Grant Nos. 11961009, 11901135), Foundation of Guizhou Science and Technology Department (Grant No. [2020]1Y015), Natural Science Research Project of Department of Education of Guizhou Province (Grant Nos. QJJ2022015, QJJ2022047).

Conflict of interest

The authors declare that they have no competing interests.