A New Parallelized Computation Method of HASC-N Difference Scheme for Inhomogeneous Time Fractional Fisher Equation

Abstract

The fractional Fisher equation has a wide range of applications in many engineering fields. The rapid numerical methods for fractional Fisher equation have momentous scientific meaning and engineering applied value. A parallelized computation method for inhomogeneous time-fractional Fisher equation (TFFE) is proposed. The main idea is to construct the hybrid alternating segment Crank-Nicolson (HASC-N) difference scheme based on alternating segment difference technology, using the classical explicit scheme and classical implicit scheme combined with Crank-Nicolson (C-N) scheme. The unique existence, unconditional stability and convergence are proved theoretically. Numerical tests show that the HASC-N difference scheme is unconditionally stable. The HASC-N difference scheme converges to $O({\tau }^{2-\alpha }+h^2)$ under strong regularity and $O({\tau }^{\alpha }+h^2)$ under weak regularity of fractional derivative discontinuity. The HASC-N difference scheme has high precision and distinct parallel computing characteristics, which is efficient for solving inhomogeneous TFFE.

1. Introduction

The time fractional Fisher equation (TFFE) is a nonlinear physical model with linear diffusion and nonlinear growth. Derived from population dynamics, chemical dynamics and other fields, it describes phenomena such as mutant gene reproduction, nonlinear evolution of population and autocatalytic chemical reaction [1,2]. Exact solutions of TFFE are difficult to be given explicitly and most of them contain special functions, such as the multivariable Mittag-Leffler function [3,4,5]. In the past two decades, with the deepening of the application of TFFE, the rapid numerical solution for TFFE has become an urgent research work [6,7].

At present, the finite difference method is still the more widely used and mature numerical method for solving TFFE. The finite difference method can achieve the precision and stability of simulation requirements well [8]. Zhang et al. (2014) [9] constructed a fully discrete scheme of TFFE by combining the finite difference method and locally discontinuous Galerkin finite element method, and discussed the stability and error estimation of the method. Alquran et al. (2015) [10] numerically solved the TFFE based on the self-collocation method and finite difference method, and analyzed the analytical and numerical solutions of the equation. Mejía and Piedrahita (2019) [11] proposed an implicit finite difference scheme for approximating TFFE with variable coefficients, and the numerical results verified the correctness of the theoretical analysis. There are also many research results on other numerical solutions of TFFE [12,13,14,15], but in most of the above numerical methods, computational efficiency has not been paid enough attention.

Due to the improvement of cluster technology and the increasing number of CPU cores, the parallelized numerical method has become one of the main methods for fast computing [16,17]. For the past few years, parallel computing has been widely used in the field of rapid numerical solutions for fractional partial differential equations (FPDEs). At present, there are three kinds of parallel algorithms using fractional differential equations: algebraic parallel algorithm, Parareal algorithm and parallel difference scheme.

Based on the algebraic parallel algorithm, Gong et al. (2013) [18] came up with a parallelized calculation method of explicit difference schemes for fractional reaction-diffusion equations, which was mainly used for parallel calculation of matrix and vector in algebraic matrix equation. Sweilam et al. (2014) [19] proposed an algebraic parallel algorithm for the time-fractional parabolic equation. This method solved the algebraic equation matrix after discrete-time in parallel. Biala and Khaliq (2018) [20] developed a C-N scheme similar to integer-order parabolic equations for nonlinear spatio-temporal fractional parabolic equations, and used the precursor-correction method respectively in MPI, OpenMP and Hybrid Version.

Using Parareal algorithm, Fu and Wang (2019) [21] constructed a Parareal algorithm to solve the space-time FPDE that models an anomalous diffusion process in a one-dimensional tube. The numerical advantages of the traditional Parareal algorithm were well preserved in this method. Yue et al. (2019) [22] proposed a multi-grid time reduction (MGRIT) algorithm based on time-varying time-grid propagators for two-dimensional fractional diffusion equations, and presented the two-level convergence theory of the algorithm. Liu et al. (2020) [23] proposed the finite volume method for time-varying fractional parabolic equations, and parallelized it with the parallel-In-time method to improve the computational efficiency of the finite volume method. Based on the Parareal method, Lorin (2020) [24] constructed the Parareal-Gorenflo algorithm for space-time FPDEs, and the spatial parallelization of this method relied on the parallelization of Riesz derivative and fast Fourier transform.

For the study of parallel difference schemes, Wang et al. (2016) [25] parallelized the implicit difference scheme of the Caputo fractional reaction-diffusion equation, and changed the serial algorithm to parallel as far as possible without changing the original serial difference scheme, to reasonably allocate computing tasks. Yang and Wu (2020) [26] proposed a parallel nature difference method for a multi-term time-fractional diffusion equation and proved that the method was unconditionally stable and convergent through theoretical analysis. Numerical experiments showed that the scheme proposed by Yang and Wu is an efficient scheme for the multi-term time-fractional diffusion equation.

To solve the problem of large computation of fractional Fisher parabolic equation, we explore the parallelization of the difference scheme for the inhomogeneous TFFE. A new parallelized computation method is proposed by using an alternate technique appropriately, which ensures the unconditional stability and spatial convergence order $O\left(h^2\right)$ of the new algorithm, and is easy to be used in many types of parallel machines.

2. HASC-N Difference Scheme for Inhomogeneous TFFE

2.1. Inhomogeneous Time Fractional Fisher Equation

Consider the inhomogeneous TFFE as follows [27,28,29]:

(1)$\left\{\begin{array}{c}\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u(x,t)}{\partial x^2}+u(x,t)(1-u(x,t))+g(x,t),(x,t)\in \left(0,L\right)\times \left(0,T\right],\hfill \\ u(x,0)=\varphi \left(x\right),x\in \left[0,L\right],\hfill \\ u(0,t)={\phi }_1\left(t\right),u(L,t)={\phi }_2\left(t\right),t\in \left(0,T\right],\hfill \end{array}\right.$

where $\varphi \left(x\right),{\phi }_1\left(t\right),{\phi }_2\left(t\right)$ are the given functions with suitable smoothness. The nonlinear source term $u(1-u)$ is a nonlinear function, $0<\alpha \le 1$.

For brevity, let $f\left(u\left(x,t\right),x,t\right)=u\left(1-u\right)$ be Lipschitz continuous with respect to u, and there exists a Lipschitz coefficient l such that $\left|f\left(u_1\right)-f\left(u_2\right)\right|\le l\left|u_1-u_2\right|$. The inhomogeneous term $g\left(x,t\right)$ is a known function independent of u. $\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}$ is the fractional derivative in the sense of Caputo:

(2)$\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{\partial u(x,s)}{\partial \tau }\frac{ds}{{(t-s)}^{\alpha }},0<\alpha <1,$

where $\Gamma (•)$ is the Gamma function. When $\alpha =1$, equation is

(3)$\frac{\partial u\left(x,t\right)}{\partial t}=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+u\left(x,t\right)(1-u\left(x,t\right)).$

The above Equation (3) is called the classical Fisher equation in general. As $\alpha$ tends to 1, according to the conclusion in reference [30,31], solution $u\left(x,t\right)$ of TFFE tends to $\tilde{u}\left(x,t\right)$ ($\tilde{u}\left(x,t\right)$ is the solution of the classical Fisher equation).

2.2. Construction of HASC-N Difference Scheme for Inhomogeneous TFFE

To construct the HASC-N difference scheme of inhomogeneous TFFE (1), the solution region $\Omega =\left\{(x,t)|0\le x\le L,0\le t\le T\right\}$ is meshed. Take the space step $h=\frac{L}{M}$ and time step $\tau =\frac{T}{N}$, where M and N are positive integers. Thus, $x_j=jh(j=1,2,\dots ,M),$ $Mh=L, t_k=k\tau (k=1,2,\dots ,N), N\tau =T$ and the grid node is $(x_j,t_k)$. Define $u_j^k=u(x_j,t_k)$, $f_j^k=f(u(x_j,t_k),x_j,t_k)$, $g_j^k=g(x_j,t_k)$.

Lemma 1

([7,32]). Suppose $0<\alpha <1$, let $y\in C^2[0,t_n]$. Then

(4)$\begin{array}{c}\frac{1}{\Gamma (1-\alpha )}{\int }_0^{t_n}\frac{y^{\prime }\left(\xi \right)d\xi }{{(t-\xi )}^{\alpha }}-\frac{{\tau }^{-\alpha }}{\Gamma (2-\alpha )}[y^n-{\displaystyle \sum _{k=1}^{n-1}}(l_{n-k-1}-l_{n-k})y^k-l_{n-1}{y}^0]\hfill \\ \le \frac{1}{\Gamma (2-\alpha )}[\frac{1-\alpha }{12}+\frac{2^{2-\alpha }}{2-\alpha }-(1+2^{\alpha })]\underset{0\le t\le t_n}{max}\left|y^{″}\left(t\right)\right|{\tau }^{2-\alpha },\hfill \end{array}$

where $l_i={(i+1)}^{(1-\alpha )}-i^{(1-\alpha )}$, $i=0,1,2,\dots ,N$.

The discrete formula is defined by Lemma 1:

(5)$D_t^{\alpha }u(x_j,t_{k+1})=\frac{{\tau }^{-\alpha }}{\Gamma (2-\alpha )}\left(l_{0}u(x_j,t_{k+1})-\sum _{i=1}^k(l_{i-1}-l_i)u(x_j,t_{k-i+1})-l_{k}u(x_j,t_0)\right).$

The method of processing nonlinear source term $f\left(u\right)$ is derived from references [33,34]:

(6)$f_j^k=2f_j^{k-1}-f_j^{k-2}+O\left({\tau }^2\right).$

Define the space second derivative discrete formula:

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

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

(9)$D{u}_j^k:=\frac{1}{2h^2}\left(u_{j-1}^{k+1}-2u_j^{k+1}+u_{j+1}^{k+1}+u_{j-1}^k-2u_j^k+u_{j+1}^k\right).$

Three difference schemes are obtained:

Classical explicit scheme,

(10)$D_t^{\alpha }u(x_j,t_{k+1})={\delta }_x^2{u}_j^k+f_j^k+g_j^k.$

Classical implicit scheme,

(11)$D_t^{\alpha }u(x_j,t_{k+1})={\delta }_x^2{u}_j^{k+1}+f_j^{k+1}+g_j^{k+1}.$

Classical C-N scheme,

(12)$D_t^{\alpha }u(x_j,t_{k+1})=D{u}_j^k+f_j^{k+\frac{1}{2}}+g_j^{k+\frac{1}{2}}.$

Further collate the above three difference schemes, we get

(13)$u_j^{k+1}=a{u}_{j-1}^k+(b_1-2a)u_j^k+a{u}_{j+1}^k+\sum _{i=2}^k{b}_i{u}_j^{k-i+1}+l_k{u}_j^0+c{f}_j^k+c{g}_j^k,$

(14)$-a{u}_{j-1}^{k+1}+(1+2a)u_j^{k+1}-a{u}_{j+1}^{k+1}=b_1{u}_j^k+\sum _{i=2}^k{b}_i{u}_j^{k-i+1}+l_k{u}_j^0+c{f}_j^{k+1}+c{g}_j^{k+1},$

(15)$-\frac{a}{2}{u}_{j-1}^{k+1}+(1+a)u_j^{k+1}-\frac{a}{2}{u}_{j+1}^{k+1}=\frac{a}{2}{u}_{j-1}^k+(b_1-a)u_j^k+\frac{a}{2}{u}_{j+1}^k+\sum _{i=2}^k{b}_i{u}_j^{k-i+1}+l_k{u}_j^0+c{f}_j^{k+\frac{1}{2}}+c{g}_j^{k+\frac{1}{2}},$

where $b_j=l_{j-1}-l_j,c={\tau }^{\alpha }\Gamma (2-\alpha ),a=\frac{c}{h^2}$.

According to the thought in references [35,36], the HASC-N difference scheme for inhomogeneous TFFE (1) is constructed:

$M+1$ points are taken at each time layer, except for the first point and the $M+1$ point on the boundary, the remaining $M-1$ points to be calculated at the same layer are divided into B segments (B is odd without losing generality). If there are n points in each segment, $nB=M-1$ (n and B are positive integers and $n\ge 3$, $B\ge 3$). The classical explicit scheme and classical implicit scheme are used alternately at the boundary points of two adjacent time layers. At the inner boundary points of two adjacent time layers, the classical explicit scheme and the classical implicit scheme are used alternately. The C-N scheme is used at the remaining points of two adjacent time layers. ● is the classical explicit scheme, ◯ is the classical implicit scheme, and ■ is the classical C-N scheme. HASC-N difference scheme construction principle is shown in Figure 1:

The HASC-N difference scheme for inhomogeneous TFFE (1) can be as follows:

(16)$\left\{\begin{array}{c}\left(I+A_{1}G\right)U^{k+1}=\left(b_{1}I-A_{2}G\right)U^k+c^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{U}^{k-i+1}+l_k{U}^0+c{A}_1{F}^{k+1}+c{A}_2{F}^k,\hfill \\ \left(I+A_{2}G\right)U^{k+2}=\left(b_{1}I-A_{1}G\right)U^{k+1}+c^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{U}^{k-i+2}+l_{k+1}{U}^0+c{A}_2{F}^{k+2}+c{A}_1{F}^{k+1},\hfill \end{array}\right.k=0,2,4\dots$

where

$G={\left[\begin{array}{cccccc}2a& -a& & & & \\ -a& 2a& -a& & & \\ & -a& 2a& -a& & \\ & & \ddots & \ddots & \ddots & \\ & & & -a& 2a& -a\\ & & & & -a& 2a\end{array}\right]}_{\left(M-1\right)\times \left(M-1\right)}$

$A_1={\left[\begin{array}{cccccc}{\theta }_1& & & & & \\ & {\theta }_2& & & & \\ & & {\theta }_3& & & \\ & & & \ddots & & \\ & & & & {\theta }_{M-2}& \\ & & & & & {\theta }_{M-1}\end{array}\right]}_{\left(M-1\right)\times \left(M-1\right)}$

${\theta }_j=\left\{\begin{array}{c}0,j=n,2n,\dots ,(B-1)n,\hfill \\ 1,j=n+1,2n+1,\dots ,(B-1)n+1,\hfill \\ \frac{1}{2},elsewhere.\hfill \end{array}\right.$, $U^k={\left(u_1^k,u_2^k,\dots ,u_{M-1}^k\right)}^T$, $c^k={\left(a{u}_0^k,0,\dots ,0,a{u}_M^k\right)}^{\mathrm{T}}$, $f^k={(f_1^k,f_2^k,\dots ,f_{M-1}^k)}^T$, $g^k={(g_1^k,g_2^k,\dots ,g_{M-1}^k)}^T$, $F^k=f^k+g^k$, $A_2=I-A_1$, I is identity matrix.

3. Existence and Uniqueness of Solution to HASC-N Difference Scheme for Inhomogeneous TFFE

Lemma 2.

The matrices $A_{1}G$ and $A_{2}G$ in HASC-N difference scheme (16) are non-negative definite matrices.

Proof. 

According to

(17)$A_{1}G={\left[\begin{array}{cccccc}2{\theta }_{1}a& -{\theta }_{1}a& & & & \\ -{\theta }_{2}a& 2{\theta }_{2}a& -{\theta }_{2}a& & & \\ & -{\theta }_{3}a& 2{\theta }_{3}a& -{\theta }_{3}a& & \\ & & \ddots & \ddots & \ddots & \\ & & & -{\theta }_{M-2}a& 2{\theta }_{M-2}a& -{\theta }_{M-2}a\\ & & & & -{\theta }_{M-1}a& 2{\theta }_{M-1}a\end{array}\right]}_{(M-1)\times (M-1)},$

we get that $A_{1}G$ is the diagonally dominant tridiagonal matrix, and the main diagonal elements are non-negative real numbers, So $A_{1}G$ is a non-negative definite matrix. Similarly, $A_{2}G$ is a non-negative definite matrix. Lemma 2 is proved. □

Theorem 1.

The solution of HASC-N difference scheme (16) for inhomogeneous TFFE (1) is existing and unique.

Proof. 

According to Lemma 2, the inverse matrices ${\left(I+A_{1}G\right)}^{-1}$ and ${\left(I+A_{2}G\right)}^{-1}$ of $I+A_{1}G$ and $I+A_{2}G$ exist, the HASC-N difference scheme (16) has a unique solution. Therefore, the Theorem 1 is proved. □

4. Stability of HASC-N Difference Scheme for Inhomogeneous TFFE

Theorem 2.

The HASC-N difference scheme (16) for inhomogeneous TFFE (1) is unconditionally stable.

Proof. 

Assume that $u_j^k$ is the HASC-N difference scheme solution for inhomogeneous TFFE (1), ${\overline{U}}_j^k$ is the approximate solution of HASC-N difference scheme for inhomogeneous TFFE (1). Error ${\epsilon }_j^k$ is defined as ${\epsilon }_j^k={\overline{U}}_j^k-u_j^k$, let ${\epsilon }_0^k={\epsilon }_M^k=0$, $E^k=\left({\epsilon }_1^k,{\epsilon }_2^k,\dots ,{\epsilon }_{M-1}^k\right),$ $k=0,1,2,\dots ,N$. Substitute the approximate solution ${\overline{U}}_j^k$ of HASC-N difference scheme and the HASC-N difference scheme solution $u_j^k$ into scheme (16), respectively, to get two equations, and make the difference between the two equations, we get

(18)$\left\{\begin{array}{c}\left(I+A_{1}G\right)E^{k+1}=\left(b_{1}I-A_{2}G\right)E^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{E}^{k-i+1}+l_k{E}^0+c{A}_1({\overline{f}}^{k+1}-f^{k+1})+c{A}_2({\overline{f}}^k-f^k),\hfill \\ \left(I+A_{2}G\right)E^{k+2}=\left(b_{1}I-A_{1}G\right)E^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{E}^{k-i+2}+l_{k+1}{E}^0+c{A}_2({\overline{f}}^{k+2}-f^{k+2})+c{A}_1({\overline{f}}^{k+1}-f^{k+1}),\hfill \end{array}\right.k=0,2,4\dots$

where $f^k=(f_1^k,f_2^k,\dots ,f_M^k)$, ${\overline{f}}^k=({\overline{f}}_1^k,{\overline{f}}_2^k,\dots ,{\overline{f}}_M^k)$, ${\overline{f}}_j^k=f({\overline{u}}_j^k,x_j,t_k)$.

Since $f(u,x,t)$ satisfifies the Lipschitz condition l, there is the Lipschitz conditional constant, we have

(19)$\left|{\overline{f}}^k-f^k\right|\le l\left|{\overline{U}}^k-u^k\right|\le l\left|E^k\right|.$

Substitute Equation (19) into Equation (18) and we get

(20)$\left\{\begin{array}{c}\left(\left(I-lc{A}_1\right)+A_{1}G\right)E^{k+1}\le \left(\left(b_{1}I+lc{A}_2\right)-A_{2}G\right)E^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{E}^{k-i+1}+l_k{E}^0,\hfill \\ \left(\left(I-lc{A}_2\right)+A_{2}G\right)E^{k+2}\le \left(\left(b_{1}I+lc{A}_1\right)-A_{1}G\right)E^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{E}^{k-i+2}+l_{k+1}{E}^0.\hfill \end{array}\right.k=0,2,4\dots$

So for simplicity, let ${\alpha }_{{\theta }_j}=lc{\theta }_j$, ${\beta }_{{\theta }_j}=lc\left(1-{\theta }_j\right)$, we have

(21)$\left\{\begin{array}{c}\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)E^{k+1}\le \left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{E}^{k-i+1}+l_k{E}^0,\hfill \\ \left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)E^{k+2}\le \left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)E^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{E}^{k-i+2}+l_{k+1}{E}^0.\hfill \end{array}\right.k=0,2,4\dots$

Define norm $∥U^k∥={∥U^k∥}_{\infty }=\underset{1\le j\le (M-1)}{max}\left\{\left|u_j^k\right|\right\}$. Known by the definition of matrices $A_1,A_2,G$, $A_{1}G$ and $A_{2}G$5 are non-negative definite matrices, and they have different non-negative characteristic values. Let the characteristic value of $A_{1}G$ be ${\lambda }_j$ and the characteristic value of $A_{2}G$ be ${\gamma }_j$, $\left|{\lambda }_j\right|\le H_1$, $\left|{\gamma }_j\right|\le H_2$, $H_1$ and $H_2$ are constants, ${\gamma }_j={\lambda }_j+K_j$, $K_j$ is constant, $j=1,2,\dots ,M-1$.

According to reference [14], there is an unequal relationship between the time process T of the TFFE and Lipschitz coefficient l. Assume that in the unequal relation between time process T and Lipschitz coefficient l, the following inequality holds: $2l_1-1\le {\alpha }_{{\theta }_j}\le min\left\{1,{\lambda }_j\right\}$, $\left(2l_1-1\right)+K_j\le {\beta }_{{\theta }_j}\le min\left\{1+K_j,{\gamma }_j\right\}$, where ${\gamma }_j={\lambda }_j+K_j$. This is bound to affect the length of time process T, however, in order to ensure the stability of HASC-N scheme, the following proofs are carried out under the premise that the above assumption is true. The results of numerical tests also confirm the feasibility of this assumption.

We will prove $∥E^k∥\le ∥E^0∥$ by mathematical induction.

When $k=0$, namely $\left\{\begin{array}{c}\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)E^1\le \left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^0,\hfill \\ \left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)E^2\le \left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)E^1+l_1{E}^0.\hfill \end{array}\right.$

Firstly, We discuss $\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)E^1\le \left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^0$. Solve for $E^1$ and take the norm of both sides, we get

$∥E^1∥\le ∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^0∥\le max\left\{\left|\frac{(b_1+{\beta }_{{\theta }_j})-{\gamma }_j}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\right|\right\}∥E^0∥.$

Case 1, $b_1+{\beta }_{{\theta }_j}>{\gamma }_j$,

(22)$max\left\{\left|\frac{(b_1+{\beta }_{{\theta }_j})-{\gamma }_j}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\right|\right\}\le \frac{(b_1+{\beta }_{{\theta }_j})-{\gamma }_j}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\le \frac{1-({\gamma }_j-{\beta }_{{\theta }_j})}{1+({\lambda }_j-{\alpha }_{{\theta }_j})}\le 1.$

Case 2, $b_1+{\beta }_{{\theta }_j}\le {\gamma }_j$,

(23)$max\left\{\left|\frac{(b_1+{\beta }_{{\theta }_j})-{\gamma }_j}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\right|\right\}\le \frac{{\gamma }_j-(b_1+{\beta }_{{\theta }_j})}{{\lambda }_j+(1-{\alpha }_{{\theta }_j})}\le 1.$

According to (22) and (23), we have $max\left\{\left|\frac{(b_1+{\beta }_{{\theta }_j})-{\lambda }_j}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\right|\right\}\le 1$, $∥E^1∥\le ∥E^0∥.$

Secondly, we discuss $\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)E^2\le \left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)E^1+l_1{E}^0$. Solve for $E^2$ and take the norm of both sides, we get

$\begin{array}{c}∥E^2∥\le ∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}\left[\right((b_1+{\alpha }_{{\theta }_j})I-A_{1}G)E^1+l_1{E}^0]∥\hfill \\ \begin{array}{ccc}& & \end{array}\le max\left\{\left|\frac{\left|(b_1+{\alpha }_{{\theta }_j})-{\lambda }_j\right|+l_1}{(1-{\beta }_{{\theta }_j})+{\gamma }_j}\right|\right\}∥E^0∥.\hfill \end{array}$

Case 1, $b_1+{\alpha }_{{\theta }_j}>{\lambda }_j$,

(24)$max\left\{\left|\frac{\left|(b_1+{\alpha }_{{\theta }_j})-{\lambda }_j\right|+l_1}{(1-{\beta }_{{\theta }_j})+{\gamma }_j}\right|\right\}\le \frac{(b_1+{\alpha }_{{\theta }_j})-{\lambda }_j+l_1}{(1-{\beta }_{{\theta }_j})+{\gamma }_j}=\frac{1-({\lambda }_j-{\alpha }_{{\theta }_j})}{1+({\gamma }_j-{\beta }_{{\theta }_j})}\le 1.$

Case 2, $b_1+{\alpha }_{{\theta }_j}\le {\lambda }_j$,

(25)$max\left\{\left|\frac{\left|(b_1+{\alpha }_{{\theta }_j})-{\lambda }_j\right|+l_1}{(1-{\beta }_{{\theta }_j})+{\gamma }_j}\right|\right\}\le \frac{{\lambda }_j-(b_1+{\alpha }_{{\theta }_j})+l_1}{(1-{\beta }_{{\theta }_j})+{\gamma }_j}=\frac{{\lambda }_j-\left((1-2l_1)+{\alpha }_{{\theta }_j}\right)}{{\lambda }_j+(1+K_j-{\beta }_{{\theta }_j})}\le 1.$

According to (24) and (25), we have $max\left\{\left|\frac{\left|(b_1+{\alpha }_{{\theta }_j})-{\lambda }_j\right|+l_1}{(1-{\beta }_{{\theta }_j})+{\gamma }_j}\right|\right\}\le 1$, $∥E^2∥\le ∥E^0∥.$

Finally, assuming that the previous layers are all true, namely $∥E^k∥\le ∥E^0∥.$

When the time layer is layer $k+1$ and layer $k+2$, $\left\{\begin{array}{c}\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)E^{k+1}\le \left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{E}^{k-i+1}+l_k{E}^0,\hfill \\ \left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)E^{k+2}\le \left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)E^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{E}^{k-i+2}+l_{k+1}{E}^0.\hfill \end{array}\right.$

Solve for $E^{k+1}$ and $E^{k+2}$ and take the norm of both sides, we get

$\left\{\begin{array}{c}∥E^{k+1}∥\le ∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{E}^{k-i+1}+l_k{E}^0]∥,\hfill \\ ∥E^{k+2}∥\le ∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)E^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{E}^{k-i+2}+l_{k+1}{E}^0]∥.\hfill \end{array}\right.$

According to $∥E^k∥\le ∥E^0∥$ and $b_j=l_{j-1}-l_j$, we get

$\begin{array}{c}\begin{array}{c}∥E^{k+1}∥\end{array}=∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{E}^{k-i+1}+l_k{E}^0]∥\hfill \\ =∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^k+b_2{E}^{k-1}+\dots +b_k{E}^1+l_k{E}^0]∥\hfill \\ \le ∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^0+b_2{E}^0+\dots +b_k{E}^0+l_k{E}^0]∥\hfill \\ =∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)E^0+l_1{E}^0]∥\hfill \\ \le ∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)+l_1]∥∥E^0∥\hfill \\ \le max\left\{\left|\frac{\left|(b_1+{\beta }_{{\theta }_j})-{\gamma }_j\right|+l_1}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\right|\right\}∥E^0∥.\hfill \end{array}$

Case 1, $b_1+{\alpha }_{{\theta }_j}>{\lambda }_j$,

(26)$max\left\{\left|\frac{\left|(b_1+{\beta }_{{\theta }_j})-{\gamma }_j\right|+l_1}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\right|\right\}\le \frac{(b_1+{\beta }_{{\theta }_j})-{\gamma }_j+l_1}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}=\frac{1-({\gamma }_j-{\beta }_{{\theta }_j})}{1+({\lambda }_j-{\alpha }_{{\theta }_j})}\le 1.$

Case 2, $b_1+{\alpha }_{{\theta }_j}\le {\lambda }_j$,

(27)$max\left\{\left|\frac{\left|(b_1+{\beta }_{{\theta }_j})-{\gamma }_j\right|+l_1}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\right|\right\}\le \frac{{\gamma }_j-(b_1+{\beta }_{{\theta }_j})+l_1}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}=\frac{{\lambda }_j-\left((1-2l_1)+{\beta }_{{\theta }_j}-K_j\right)}{{\lambda }_j+(1-{\alpha }_{{\theta }_j})}\le 1.$

According to (26) and (27), we have $max\left\{\left|\frac{\left|(b_1+{\beta }_{{\theta }_j})-{\gamma }_j\right|+l_1}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\right|\right\}\le 1$, $∥E^{k+1}∥\le ∥E^0∥.$

Similarly, according to $∥E^{k+1}∥\le ∥E^0∥$ and $b_j=l_{j-1}-l_j$, we get

$\begin{array}{c}\begin{array}{c}∥E^{k+2}∥\end{array}=∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)E^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{E}^{k-i+2}+l_{k+1}{E}^0]∥\hfill \\ =∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)E^{k+1}+b_2{E}^k+\dots +b_{k+1}{E}^1+l_{k+1}{E}^0]∥\hfill \\ \le ∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)E^0+b_2{E}^0+\dots +b_{k+1}{E}^0+l_{k+1}{E}^0]∥\hfill \\ =∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)E^0+l_1{E}^0]∥\hfill \\ \le ∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)+l_1]∥∥E^0∥\hfill \\ \le max\left\{\left|\frac{\left|(b_1+{\alpha }_{{\theta }_j})-{\lambda }_j\right|+l_1}{(1-{\beta }_{{\theta }_j})+{\gamma }_j}\right|\right\}∥E^0∥\le ∥E^0∥.\hfill \end{array}$

In summary, stability is proved. □

5. Convergence of HASC-N Difference Scheme for Inhomogeneous TFFE

Lemma 3

([6]). Suppose $0<\alpha <1$, let $y\in C^2[0,t_{n+1}]$. Then we have

$\begin{array}{c}\frac{{\partial }^{\alpha +1}y\left(t_{n+1}\right)}{\partial t^{\alpha +1}}=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{{\partial }^{2}y\left(\xi \right)}{\partial {\xi }^2}\frac{d\xi }{{\left(t_{n+1}-\xi \right)}^{\alpha }}=\frac{1}{\Gamma (1-\alpha )}{\displaystyle \sum _{j=0}^n}{\int }_{j\tau }^{\left(j+1\right)\tau }\frac{{\partial }^{2}y\left(\xi \right)}{\partial {\xi }^2}\frac{d\xi }{{\left(t_{n+1}-\xi \right)}^{\alpha }}\hfill \\ \le \frac{1}{\Gamma (1-\alpha )}\underset{0\le t\le t_{n+1}}{max}\left\{\left|\frac{{\partial }^{2}y\left(t\right)}{\partial t^2}\right|\right\}{\displaystyle \sum _{j=0}^n}{\int }_{j\tau }^{\left(j+1\right)\tau }\frac{d\xi }{{\left(t_{n+1}-\xi \right)}^{\alpha }}=\frac{C_y}{\Gamma (1-\alpha )}{\displaystyle \sum _{j=0}^n}{\int }_{j\tau }^{\left(j+1\right)\tau }\frac{d\xi }{{\left(t_{n+1}-\xi \right)}^{\alpha }}.\hfill \\ =\frac{C_y{\tau }^{1-\alpha }}{\Gamma (2-\alpha )}{\displaystyle \sum _{j=0}^n}{\int }_{j\tau }^{\left(j+1\right)\tau }\left[{\left(n+1-j\right)}^{1-\alpha }-{\left(n-j\right)}^{1-\alpha }\right]\le \frac{{\left(n+1\right)}^{1-\alpha }{C}_y}{\Gamma (2-\alpha )}{\tau }^{1-\alpha },\hfill \end{array}$ where $C_y=\underset{0\le t\le t_{n+1}}{max}\left\{\left|\frac{{\partial }^{2}y\left(t\right)}{\partial t^2}\right|\right\}.$

The solution of the inhomogeneous TFFE (1) satisfies the strong regularity condition as follows,

(28)$\frac{{\partial }^{\gamma }u}{\partial t^{\gamma }}\in C\left(\left[0,L\right]\times \left[0,T\right]\right),\frac{{\partial }^{\delta }u}{\partial x^{\delta }}\in C\left(\left[0,L\right]\times \left[0,T\right]\right),$

where $\gamma \in \left\{1,2\right\}$ and $\delta \in \left\{0,1,2,3,4\right\}.$

Consider the explicit scheme on the time layer $k+1$,

(29)$D_t^{\alpha }u(x_j,t_{k+1})=\frac{1}{h^2}(u_{j-1}^k-2u_j^k+u_{j+1}^k)+F_j^k,$

and the implicit scheme on the time layer $k+2$,

(30)$D_t^{\alpha }u(x_j,t_{k+2})=\frac{1}{h^2}(u_{j-1}^{k+2}-2u_j^{k+2}+u_{j+1}^{k+2})+F_j^{k+2}.$

Taylor expansion is performed at $u_j^{k+1}$ for scheme (29) and scheme (30) to obtain truncation error,

(31)$R_1(\tau ,h)=\frac{{\partial }^{\alpha }u(x_j,t_{k+1})}{\partial t^{\alpha }}-u_{xx}+\tau u_{xxt}-\frac{1}{12}{h}^2{u}_{xxxx}-\tau \frac{\partial F\left(u\right)}{\partial t}+O({\tau }^{2-\alpha }+h^2),$

(32)$R_2(\tau ,h)=\frac{{\partial }^{\alpha }u(x_j,t_{k+1})}{\partial t^{\alpha }}+\tau \frac{{\partial }^{\alpha +1}u(x_j,t_{k+1})}{\partial t^{\alpha +1}}-u_{xx}-\tau u_{xxt}-\frac{1}{12}{h}^2{u}_{xxxx}+\tau \frac{\partial F\left(u\right)}{\partial t}+O({\tau }^{2-\alpha }+h^2).$

Consider the C-N scheme on the time layer $k+1$,

(33)$D_t^{\alpha }u(x_j,t_{k+1})=\frac{1}{2h^2}(u_{j-1}^{k+1}-2u_j^{k+1}+u_{j+1}^{k+1}+u_{j-1}^k-2u_j^k+u_{j+1}^k)+\frac{1}{2}\left(F_j^{k+1}+F_j^k\right),$

and the C-N scheme on the time layer $k+2$,

(34)$D_t^{\alpha }u(x_j,t_{k+2})=\frac{1}{2h^2}(u_{j-1}^{k+2}-2u_j^{k+2}+u_{j+1}^{k+2}+u_{j-1}^{k+1}-2u_j^{k+1}+u_{j+1}^{k+1})+\frac{1}{2}\left(F_j^{k+2}+F_j^{k+1}\right).$

Taylor expansion is performed at $u_j^{k+1}$ for scheme (33) and scheme (34) to obtain truncation error,

(35)$R_3(\tau ,h)=\frac{{\partial }^{\alpha }u(x_j,t_{k+1})}{\partial t^{\alpha }}-u_{xx}+\frac{\tau }{2}{u}_{xxt}-\frac{1}{12}{h}^2{u}_{xxxx}-\frac{\tau }{2}\frac{\partial F\left(u\right)}{\partial t}+O({\tau }^{2-\alpha }+h^2),$

(36)$R_4(\tau ,h)=\frac{{\partial }^{\alpha }u(x_j,t_{k+1})}{\partial t^{\alpha }}+\tau \frac{{\partial }^{\alpha +1}u(x_j,t_{k+1})}{\partial t^{\alpha +1}}-u_{xx}-\frac{\tau }{2}{u}_{xxt}-\frac{1}{12}{h}^2{u}_{xxxx}+\frac{\tau }{2}\frac{\partial F\left(u\right)}{\partial t}+O({\tau }^{2-\alpha }+h^2).$

According to Lemma 1, the calculation precision of $\frac{{\partial }^{\alpha }u(x_j,t_{k+1})}{\partial t^{\alpha }}$ is $O\left({\tau }^{2-\alpha }\right)$, the calculation precision of $\tau \frac{{\partial }^{\alpha +1}u(x_j,t_{k+1})}{\partial t^{\alpha +1}}$ is also $O\left({\tau }^{2-\alpha }\right)$ according to Lemma 3. By using explicit and implicit schemes alternately at the inner boundary points of adjacent time layers, two basic error components with opposite signs are generated, and the two partially cancel each other, so as to obtain ideal calculation precision.

Add (31) and (32) to get

$R_1(\tau ,h)+R_2(\tau ,h)=2\frac{{\partial }^{\alpha }u(x_j,t_{k+1})}{\partial t^{\alpha }}+\tau \frac{{\partial }^{\alpha +1}u(x_j,t_{k+1})}{\partial t^{\alpha +1}}-2u_{xx}-\frac{1}{6}{h}^2{u}_{xxxx}+O({\tau }^{2-\alpha }+h^2).$

At the inner boundary points, the calculation precision is $O({\tau }^{2-\alpha }+h^2)$. Similarly, C-N scheme is alternately used at interior points of adjacent time layers, (35) and (36) are added to obtain $R_3(\tau ,h)+R_4(\tau ,h)=2\frac{{\partial }^{\alpha }u(x_j,t_{k+1})}{\partial t^{\alpha }}+\tau \frac{{\partial }^{\alpha +1}u(x_j,t_{k+1})}{\partial t^{\alpha +1}}-2u_{xx}-\frac{1}{6}{h}^2{u}_{xxxx}+O({\tau }^{2-\alpha }+h^2).$

So the precision at the interior points is also $O({\tau }^{2-\alpha }+h^2)$.

Theorem 3.

Assuming that the solution of Equation (1) satisfies the strong regularity condition (28), the HASC-N difference scheme (16) for inhomogeneous TFFE (1) is convergent, and $∥e^n∥\le C({\tau }^{2-\alpha }+h^2),n=1,2,\dots ,N,C>0.$

Proof. 

Let $U_j^k=U(x_j,t_k)$ be the exact solution of inhomogeneous TFFE (1) at $t=t_k,$ $x=x_j$ under strong regularity. Define $e_j^k=U_j^k-u_j^k, 1\le j\le M-1, e_0^k=e_M^k=0,$ $e^k=(e_1^k,\dots ,e_{M-1}^k), e^0=0.$

Substitute the exact solution $U_j^k$ and the HASC-N difference scheme solution $u_j^k$ into scheme (16), respectively, to get two equations, and make the difference between the two equations, we get

(37)$\left\{\begin{array}{c}\left(I+A_{1}G\right)e^{k+1}=\left(b_{1}I-A_{2}G\right)e^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{e}^{k-i+1}+c{A}_1({\overline{f}}^{k+1}-f^{k+1})+c{A}_2({\overline{f}}^k-f^k)+{\tilde{R}}^{k+1},\hfill \\ \left(I+A_{2}G\right)e^{k+2}=\left(b_{1}I-A_{1}G\right)e^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{e}^{k-i+2}+c{A}_2({\overline{f}}^{k+2}-f^{k+2})+c{A}_1({\overline{f}}^{k+1}-f^{k+1})+{\tilde{R}}^{k+2},\hfill \end{array}\right.k=0,2,4\dots$

where ${\tilde{R}}^{k+1}={\tau }^{\alpha }O({\tau }^{2-\alpha }+h^2)$, $∥{\tilde{R}}^{k+1}∥\le C_1{\tau }^{\alpha }({\tau }^{2-\alpha }+h^2)=C_1({\tau }^2+{\tau }^{\alpha }{h}^2)$, $C_1$ is a real constant.

Similar to stability analysis, convergence is studied:

When $k=0,\left\{\begin{array}{c}\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)e^1\le \left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)e^0+{\tilde{R}}^1,\hfill \\ \left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)e^2\le \left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)e^1+{\tilde{R}}^2.\hfill \end{array}\right.$

Solve for $e^1$ and $e^2$ and take the norm, we get

$\left\{\begin{array}{c}∥e^1∥\le ∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)e^0+{\tilde{R}}^1]∥,\hfill \\ ∥e^2∥\le ∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)e^1+{\tilde{R}}^2]∥.\hfill \end{array}\right.$

Firstly, $∥e^1∥\le ∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)e^0+{\tilde{R}}^1]∥,$ and, $e^0=0$, we get $∥e^1∥\le ∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}{\tilde{R}}^1∥\le max\left\{\frac{1}{1+\left|{\lambda }_j-{\alpha }_{{\theta }_j}\right|}\right\}∥{\tilde{R}}^1∥\le l_0^{-1}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2).$

Secondly, $\begin{array}{c}∥e^2∥\le ∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)e^1+{\tilde{R}}^2]∥\hfill \\ \le l_1^{-1}∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)+l_1]∥||{\tilde{R}}^2||\hfill \\ \le l_1^{-1}max\left\{\left|\frac{\left|\left(b_1+{\alpha }_{{\theta }_j}\right)-{\lambda }_j\right|+l_1}{\left(1-{\beta }_{{\theta }_j}\right)+{\gamma }_j}\right|\right\}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2)\le l_1^{-1}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2).\hfill \end{array}$

Assuming that $∥e^k∥\le l_{k-1}^{-1}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2)$ is true for all the previouslayers. When the time layer are $k+1$ and $k+2$,

$\left\{\begin{array}{c}\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)e^{k+1}\le \left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)e^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{e}^{k-i+1}+{\tilde{R}}^{k+1},\hfill \\ \left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)e^{k+2}\le \left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)e^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{e}^{k-i+2}+{\tilde{R}}^{k+2}.\hfill \end{array}\right.$

Solve for $e^{k+1}$ and $e^{k+2}$ and take the norm, we get

$\left\{\begin{array}{c}∥e^{k+1}∥\le ∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)e^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{e}^{k-i+1}+{\tilde{R}}^{k+1}]∥,\hfill \\ ∥e^{k+2}∥\le ∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)e^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{e}^{k-i+2}+{\tilde{R}}^{k+2}]∥.\hfill \end{array}\right.$

Then $\begin{array}{c}∥e^{k+1}∥\le ∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)e^k+{\displaystyle \sum _{i=2}^k}{b}_{\mathrm{i}}{e}^{k-i+1}+{\tilde{R}}^{k+1}]∥\hfill \\ \le l_k^{-1}∥{\left((1-{\alpha }_{{\theta }_j})I+A_{1}G\right)}^{-1}[\left((b_1+{\beta }_{{\theta }_j})I-A_{2}G\right)+b_2+\dots +b_k+l_k]∥||{\tilde{R}}^{k+1}||\hfill \\ \le l_k^{-1}max\left\{\left|\frac{\left|(b_1+{\beta }_{{\theta }_j})-{\gamma }_j\right|+l_1}{(1-{\alpha }_{{\theta }_j})+{\lambda }_j}\right|\right\}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2)\le l_k^{-1}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2),\hfill \end{array}$

$\begin{array}{c}∥e^{k+2}∥\le ∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)e^{k+1}+{\displaystyle \sum _{i=2}^{k+1}}{b}_{\mathrm{i}}{e}^{k-i+2}+{\tilde{R}}^{k+2}]∥\hfill \\ \le l_{k+1}^{-1}∥{\left((1-{\beta }_{{\theta }_j})I+A_{2}G\right)}^{-1}[\left((b_1+{\alpha }_{{\theta }_j})I-A_{1}G\right)+b_2+\dots +b_{k+1}+l_{k+1}]∥||{\tilde{R}}^{k+2}||\hfill \\ \le l_{k+1}^{-1}max\left\{\frac{\left|(b_1+{\alpha }_{{\theta }_j})-{\lambda }_j\right|+l_1}{(1-{\beta }_{{\theta }_j})+{\gamma }_j}\right\}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2)\le l_{k+1}^{-1}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2).\hfill \end{array}$

In conclusion, we prove that $∥e^n∥\le l_{n-1}^{-1}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2), n=1,2,\dots ,N.$ From $\underset{n\to \infty }{lim}\frac{l_n^{-1}}{n^{\alpha }}=\underset{n\to \infty }{lim}\frac{n^{-\alpha }}{n^{(1-\alpha )}-{(n-1)}^{(1-\alpha )}}=\underset{n\to \infty }{lim}\frac{n^{-1}}{1-{(1-\frac{1}{n})}^{(1-\alpha )}}=\frac{1}{1-\alpha },$ there exists $C=\frac{{\left(n\tau \right)}^{\alpha }}{1-\alpha }{C}_1$ such that $∥e^n∥\le l_{n-1}^{-1}{C}_1({\tau }^2+{\tau }^{\alpha }{h}^2)\le C({\tau }^{2-\alpha }+h^2), n=1,2,\dots ,N.$ Therefore, proof is completed, and the convergence order is $O({\tau }^{2-\alpha }+h^2).$ □

Remark 1.

The exact solution $u(x,t)$ of the inhomogeneous TFFE (1) satisfies the strong regularity condition (28). In this case, the convergence order of HASC-N difference scheme (16) is $O({\tau }^{2-\alpha }+h^2)$. In general, the condition of strong regularity is too harsh. The exact solution of inhomogeneous TFFE (1) cannot meet this requirement under some conditions, such as solving the solution of inhomogeneous TFFE with initial singularities, and the corresponding theoretical analysis has obvious limitations. Nevertheless, the conclusion is significant because Theorem 3 at least rigorously proves the theoretical correctness of HASC-N difference scheme (16) in a certain range.

Remark 2.

The time fractional derivative of the exact solution $u(x,t)$ for the inhomogeneous TFFE (1) is a discontinuous function at initial time, namely $\frac{{\partial }^{\gamma }u}{\partial t^{\gamma }}$ does not exist in some regions of $[0,L]\times [0,T]$, where $\gamma \in \{0,1,2\}$. In this case, the strong regularity condition (28) cannot be satisfied, resulting in the initial singularity of inhomogeneous TFFE (1). Consider the following two cases:

(1) The partial derivative of the solution $u(x,t)$ in the spatial direction satisfies $\frac{{\partial }^{\delta }u}{\partial x^{\delta }}\in C([0,L]\times$$[0,T]), \delta \in \{0,1,2,3,4\}$. The HASC-N difference scheme (16) converges to $O\left({\tau }^{\alpha }\right)$ in the temporal direction (consistent with the conclusions of references [37,38,39]), and it converges to $O\left(h^2\right)$ in the spatial direction.

(2) The partial derivative of the solution $u(x,t)$ in the spatial direction is a discontinuous function, that is, $\frac{{\partial }^{\delta }u}{\partial x^{\delta }}$ does not exist in some regions of $[0,L]\times [0,T]$, where $\delta \in \{0,1,2,3,4\}$. In this case, the local truncation error of HASC-N difference scheme (16) lacks clear overall control. Even if the loose discrete $L_2$ norm is used as a measure, the order of local truncation error is not clear [7,40]. Therefore, the analytic path of spatial and temporal convergence order based on strong regularity condition is no longer effective.

6. Numerical Tests

The numerical tests are based on Intel Core I5-5200 CPU @2.20 GHz, dual-core processor, and carried out in MatlabR2018b environment. Numerical tests verify the correctness of the above theoretical analysis.

Example 1

([41]). Consider the inhomogeneous TFFE with a smooth solution:

(38)$\left\{\begin{array}{c}\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u(x,t)}{\partial x^2}+u(x,t)(1-u(x,t))+g_1\left(x,t\right),(x,t)\in \left(0,1\right)\times \left(0,1\right],\hfill \\ u(x,0)=0,x\in \left[0,1\right],\hfill \\ u(0,t)=u(1,t)=0,t\in \left(0,1\right].\hfill \end{array}\right.$

where $g_1\left(x,t\right)=24t^{(4-\alpha )}sin\left(2\pi x\right)/\Gamma (5-\alpha )+4{\pi }^2{t}^{4}sin\left(2\pi x\right)-t^{4}sin\left(2\pi x\right)(1-t^{4}sin\left(2\pi x\right)), 0<\alpha \le 1.$ Exact solution of the inhomogeneous TFFE (38): $u(x,t)=t^{4}sin\left(2\pi x\right).$

When $\alpha =0.7, N=100, M=71$, the exact solution surface, C-N scheme solution surface and HASC-N scheme solution surface are as follows:

According to Figure 2, Figure 3 and Figure 4, the surfaces of the two schemes are consistent with those of the exact solution and the surface of the HASC-N difference scheme is smooth. It is shown below that when α is of different values, the HASC-N scheme solution is compared with the exact solution at $t=0.5$. The HASC-N scheme solution approximates the exact solution well, and the calculation results are shown in Table 1:

Let $U_j^k$ be the exact solution, and ${\tilde{U}}_j^k$ be the solution of C-N scheme and HASC-N difference scheme. Define the Sum of Relative Error for every Time layer, $SRET\left(k\right)={\displaystyle \sum _{j=1}^M}\frac{\left|{\tilde{U}}_j^k-U_j^k\right|}{U_j^k}.$ The purpose is to test the stability of HASC-N difference scheme.

The SRET values of the two schemes are shown in Figure 5. When $\alpha =0.7, N=100, M=71$, the SRET values of the two schemes tend to 0 with the increase of the time grid numbers. Therefore, the C-N scheme and HASC-N difference scheme of inhomogeneous TFFE (38) are stable, and the results in Figure 5 verify the correctness of Theorem 2.

The spatial convergence order and the temporal convergence order of HASC-N difference scheme are compared. The error $E_{\infty }\left(m,\tau \right)$, the error $E_{\infty }\left(h,n\right)$, the spatial convergence order $Order1$ and temporal convergence order $Order2$ are defined as follows [42,43]:

$E_{\infty }\left(m,\tau \right)=\underset{0\le k\le N}{max}\left|{\tilde{U}}_m^k-U_m^k\right|,\left(0\le m\le M\right),E_{\infty }\left(h,n\right)=\underset{0\le j\le M}{max}\left|{\tilde{U}}_j^n-U_j^n\right|,\left(0\le n\le N\right).$

$Order1=\frac{ln\left(E_{\infty }\left(h_1,n\right)/\phantom{E_{\infty }\left(h_1,n\right)E_{\infty }\left(h_2,n\right)}{E}_{\infty }\left(h_2,n\right)\right)}{ln\left(h_1/\phantom{h_1{h}_2}{h}_2\right)}$, $Order2=\frac{ln\left(E_{\infty }(m,{\tau }_1)/\phantom{E_{\infty }(m,{\tau }_1)E_{\infty }(m,{\tau }_2)}{E}_{\infty }(m,{\tau }_2)\right)}{ln\left({\tau }_1/\phantom{{\tau }_1{\tau }_2}{\tau }_2\right)}$.

To verify the spatial convergence order of HASC-N difference scheme, take $M=21,41,81,161$ and $\tau =h^2/4.$ Table 2 shows that the spatial convergence order of HASC-N difference scheme is $O\left(h^2\right)$, and its error decreases gradually with the increase of space step. The theoretical analysis is validated by numerical test data.

Calculate the temporal convergence order of HASC-N difference scheme. Fixed space step h = 1/101, namely, M = 101 and let N = 16, 32, 64, 128. As can be seen from Table 3, the temporal convergence order of HASC-N scheme reaches $O\left({\tau }^{2-\alpha }\right)$, and the error of HASC-N difference scheme decreases gradually with the increase of time step.

The numerical experimental data in Table 2 and Table 3 can correspond to the conclusion of spatial convergence order $O\left(h^2\right)$ and temporal convergence order $O\left({\tau }^{2-\alpha }\right)$ in Theorem 3.

Speed-up ratio $Sp=T/Tp$ (T is the CPU time in C-N scheme, $Tp$ is the CPU time in HASC-N difference scheme) and efficiency $Ep=Sp/P$ (p is the number of processor cores) [17]. Take $\alpha =0.7, N=100$, space grid points $M=201$, 401, 601, 801, 1001, 1201. Table 4 shows the CPU time of C-N scheme solution and HASC-N scheme solution, speed-up ratio ($Sp$) and efficiency ($Ep$) of HASC-N scheme solution.

According to the comparative analysis in Table 4, the computational efficiency of the HASC-N difference scheme of inhomogeneous TFFE (38) is obviously better than that of the C-N scheme. With the encryption of the spatial grid, the computational time advantage of the HASC-N difference scheme is more and more prominent than that of the C-N scheme. The speed-up ratio of the HASC-N difference scheme and C-N scheme is above 3, and the efficiency is about $1.5$. The results show that the HASC-N difference scheme has obvious parallel computing characteristics.

Example 2.

Consider the inhomogeneous TFFE for the discontinuity of the time fractional derivative at the initial time:

(39) $\left\{\begin{array}{c}\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u(x,t)}{\partial x^2}+u(x,t)(1-u(x,t))+g_2\left(x,t\right),(x,t)\in \left(0,1\right)\times \left(0,1\right],\hfill \\ u(x,0)=0,x\in \left[0,1\right],\hfill \\ u(0,t)=u(1,t)=0,t\in \left(0,1\right].\hfill \end{array}\right.$

where $g_2\left(x,t\right)=\Gamma (1+\alpha )x\left(x-1\right)-2t^{\alpha }-x\left(x-1\right)t^{\alpha }(1-x\left(x-1\right)t^{\alpha }), 0<\alpha \le 1.$ Exact solution of the inhomogeneous TFFE (39): $u(x,t)=x\left(x-1\right)t^{\alpha }.$

When $\alpha =0.5, N=100, M=71$, the exact solution surface and HASC-N scheme solution surface are as follows:

When $\alpha =0.5,$ the solution of inhomogeneous TFFE (39) has initial singularity near $t=0$, and the solution is smooth away from $t=0$, as shown in Figure 6 and Figure 7.

In Table 5, the number of space grids $M=21,41,81,161$ and let $\tau =h^2/4.$ In Table 6, the number of time grids $N=8,16,32,64,$ and the fixed spatial step $h=1/101$, namely, $M=101.$ When α is set to different values, it can be seen from Table 5 that the spatial convergence order of HASC-N difference scheme is $O\left(h^2\right)$, and the temporal convergence order given in Table 6 can reach $O\left({\tau }^{\alpha }\right)$. Therefore, the fractional derivative $\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}$ of inhomogeneous TFFE (39) is a discontinuous function (satisfying the weak regularity conditions), the partial derivative of spatial direction satisfies $\frac{{\partial }^{\delta }u}{\partial x^{\delta }}\in C\left(\left[0,L\right]\times \left[0,T\right]\right)(\delta \in 0,1,2,3,4)$, and the solution of HASC-N difference scheme converges to $O({\tau }^{\alpha }+h^2)$, which verifies the first statement in Remark 2.

Example 3.

Consider the inhomogeneous TFFE for the discontinuity of the temporal fractional derivative at the initial time and the discontinuity of the spatial derivative:

(40) $\left\{\begin{array}{c}\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u(x,t)}{\partial x^2}+u(x,t)(1-u(x,t))+g_3\left(x,t\right),(x,t)\in \left(0,1\right)\times \left(0,1\right],\hfill \\ u(x,0)=0,x\in \left[0,1\right],\hfill \\ u(0,t)=u(1,t)=0,t\in \left(0,1\right].\hfill \end{array}\right.$

where $0<\alpha \le 1$,

$\begin{array}{c}{g}_3\left(x,t\right)=\frac{\Gamma \left(2\right)}{\Gamma (2-\alpha )}sin\left(\pi x\right)+\frac{\Gamma \left(3\right)}{\Gamma (3-\alpha )}{\left[x\left(1-x\right)\right]}^{0.5}{t}^{2-\alpha }+{\pi }^2{t}^{\alpha }sin\left(\pi x\right)\hfill \\ \begin{array}{ccc}& & \end{array}-\frac{t^2}{2}[-\frac{1}{2}{\left[x\left(1-x\right)\right]}^{-1.5}{(1-2x)}^2-2{\left[x\left(1-x\right)\right]}^{-0.5}]\hfill \\ \begin{array}{ccc}& & \end{array}-(t^{\alpha }sin\left(\pi x\right)+t^2{\left[x\left(1-x\right)\right]}^{0.5})(1-(t^{\alpha }sin\left(\pi x\right)+t^2{\left[x\left(1-x\right)\right]}^{0.5})).\hfill \end{array}$

Exact solution of the inhomogeneous TFFE (40):

(41) $u(x,t)=t^{\alpha }sin\left(\pi x\right)+t^2{\left[x\left(1-x\right)\right]}^{0.5}.$

When α is set to different values, $N=100, M=71$, the HASC-N scheme solution is compared with the exact solution at $t=0.5$, and the calculation results are shown in Table 7:

According to the inhomogeneous TFFE (40) and its exact solution $u(x,t)$ (41), the equation has the initial singularity. In addition, it has singularity near the boundary of $u(0,t)$ and $u(1,t)$. The Equation (40) meets the weak regularity condition and produces a certain disturbance to the HASC-N difference scheme. According to the analysis of Table 7, although there is some error between the HASC-N scheme solution and the exact solution, the approximation effect is still satisfactory.

The solution (41) of the inhomogeneous TFFE (40) has an initial singularity (satisfying the weak regularity condition), and the partial derivative of the spatial direction $\frac{{\partial }^{\delta }u}{\partial x^{\delta }}(\delta \in 0,1,2,3,4)$ is a discontinuous function on $\left[0,L\right]\times \left[0,T\right]$. In order to explore whether the truncation error of HASC-N scheme solution for inhomogeneous TFFE (40) has a clear global control, The loose $L_2$ norm is used as a measure, and the $L_2$ norm is defined as follows [44]:

$E_2\left(m,\tau \right)={\left\{{\displaystyle \sum _{k=1}^N}{\left({\tilde{U}}_m^k-U_m^k\right)}^2\tau \right\}}^{\frac{1}{2}},\left(0\le m\le M\right),$

$E_2\left(h,n\right)={\left\{{\displaystyle \sum _{j=1}^M}{\left({\tilde{U}}_j^n-U_j^n\right)}^{2}h\right\}}^{\frac{1}{2}},\left(0\le n\le N\right).$

So the spatial convergence order $Order3$ and temporal convergence order $Order4$ are defined as:

$Order3=\frac{ln\left(E_2\left(h_1,n\right)/\phantom{E_2\left(h_1,n\right)E_2\left(h_2,n\right)}{E}_2\left(h_2,n\right)\right)}{ln\left(h_1/\phantom{h_1{h}_2}{h}_2\right)}$, $Order4=\frac{ln\left(E_2(m,{\tau }_1)/\phantom{E_2(m,{\tau }_1)E_2(m,{\tau }_2)}{E}_2(m,{\tau }_2)\right)}{ln\left({\tau }_1/\phantom{{\tau }_1{\tau }_2}{\tau }_2\right)}$.

As shown in Table 8 and Table 9, the local truncation error of the HASC-N difference scheme lacks a clear overall control. Even if the loose discrete $L_2$ norm is used as the measure, the local truncation error has no definite order. Therefore, the analysis path of convergence order based on strong regularity condition (28) is no longer effective. This verifies the second statement in Remark 2.

7. Conclusions

Most schemes with parallelism are not unconditionally stable for a long time, or the stability meets the requirements but the space has only precision $O\left(h\right)$ [45,46]. The HASC-N difference scheme for inhomogeneous TFFE is constructed in this paper, which is unconditionally stable. The convergence order of HASC-N difference scheme is $O({\tau }^{2-\alpha }+h^2)$ under the strong regularity condition, and $O({\tau }^{\alpha }+h^2)$ under the weak regularity condition that the time-fractional derivative is discontinuous at the initial time and the space derivative is continuous. Under the weak regularity condition that the time-fractional derivative is discontinuous at the initial time and the spatial derivative is discontinuous, the error of the HASC-N difference scheme lacks a clear global control and does not specify the convergence order. Therefore, the analysis path of convergence order based on strong regularity conditions is no longer effective.

The HASC-N difference scheme has obvious parallel computing properties. The localization characteristics of the HASC-N difference scheme in computing and communication will become more and more remarkable with the continuous encryption of space grid points, which is suitable for parallelized computing systems with distributed storage. The numerical tests verify the theoretical analysis and show that the HASC-N difference scheme in this paper is high-efficient in solving inhomogeneous TFFE.

Author Contributions

Conceptualization, R.L., X.Y. and P.L.; methodology, R.L., X.Y. and P.L.; software, R.L., X.Y. and P.L.; validation, R.L., X.Y. and P.L.; formal analysis, R.L., X.Y. and P.L.; investigation, R.L., X.Y. and P.L.; resources, R.L., X.Y. and P.L.; data curation, R.L., X.Y. and P.L.; writing—original draft preparation, R.L., X.Y. and P.L.; writing—review and editing, R.L., X.Y. and P.L.; visualization, R.L., X.Y. and P.L.; supervision, X.Y.; project administration, X.Y.; funding acquisition, X.Y. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this research. Data sharing does not apply to this research.

Acknowledgments

We would like to thank Lifei Wu of North China Electric Power University for many helpful discussions.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Construction principle of HASC-N difference scheme.

Figure 2. Exact solution surface for Example 1.

Figure 3. C-N scheme solution surface for Example 1. ([Forumla omitted. See PDF.]).

Figure 4. HASC-N scheme solution surface for Example 1. ([Forumla omitted. See PDF.]).

Figure 5. SRET values in two schemes for Example 1. ([Forumla omitted. See PDF.]).

Figure 6. Exact solution surface for Example 2. ([Forumla omitted. See PDF.]).

Figure 7. HASC-N scheme solution surface for Example 2. ([Forumla omitted. See PDF.]).

Table 1

Comparison between exact solution and HASC-N scheme solution for Example 1 ($t=0.5$).

$\mathit{\alpha }$		x
		0.25	0.5	0.75
0.3	Exact solution	5.752108 × 10${}^{-2}$	2.549964 × 10${}^{-3}$	−5.763390 × 10${}^{-2}$
	HASC-N scheme solution	5.517192 × 10${}^{-2}$	2.529212 × 10${}^{-3}$	−5.600727 × 10${}^{-2}$
0.5	Exact solution	5.752108 × 10${}^{-2}$	2.549964 × 10${}^{-3}$	−5.763390 × 10${}^{-2}$
	HASC-N scheme solution	5.565180 × 10${}^{-2}$	2.490347 × 10${}^{-3}$	−5.608690 × 10${}^{-2}$
0.7	Exact solution	5.752108 × 10${}^{-2}$	2.549964 × 10${}^{-3}$	−5.763390 × 10${}^{-2}$
	HASC-N scheme solution	5.605673 × 10${}^{-2}$	2.486836 × 10${}^{-3}$	−5.629021 × 10${}^{-2}$

Table 2

Numerical error and spatial convergence order of HASC-N difference scheme for Example 1.

$\mathit{\alpha }$	M	N	${\mathit{E}}_{\mathit{\infty }}\left(\mathit{h},\mathit{n}\right)$	$\mathit{Order}1$
0.5	21	100	5.805049 × 10${}^{-3}$	$--$
	41	400	1.535960 × 10${}^{-3}$	1.987252
	81	1600	3.930984 × 10${}^{-4}$	2.001614
	161	6400	9.953251 × 10${}^{-5}$	1.999513
0.7	21	100	5.272847 × 10${}^{-3}$	$--$
	41	400	1.326208 × 10${}^{-3}$	2.062997
	81	1600	3.261035 × 10${}^{-4}$	2.060378
	161	6400	7.997470 × 10${}^{-5}$	2.045991
0.9	21	100	4.857505 × 10${}^{-3}$	$--$
	41	400	1.220010 × 10${}^{-3}$	2.065118
	81	1600	3.029576 × 10${}^{-4}$	2.045922
	161	6400	7.518864 × 10${}^{-5}$	2.028652

Table 3

Numerical error and temporal convergence order of HASC-N difference scheme for Example 1.

$\mathit{\alpha }$	M	N	${\mathit{E}}_{\mathit{\infty }}\left(\mathit{m},\mathit{\tau }\right)$	$\mathit{Order}2$
0.5	101	16	1.035923 × 10${}^{-2}$	$--$
		32	3.687553 × 10${}^{-3}$	1.490181
		64	1.303842 × 10${}^{-3}$	1.499895
		128	4.594546 × 10${}^{-4}$	1.504774
0.7	101	16	9.512878 × 10${}^{-3}$	$--$
		32	3.879871 × 10${}^{-3}$	1.293873
		64	1.578013 × 10${}^{-3}$	1.297899
		128	6.392140 × 10${}^{-4}$	1.303738
0.9	101	16	8.555828 × 10${}^{-3}$	$--$
		32	3.996333 × 10${}^{-3}$	1.098231
		64	1.865814 × 10${}^{-3}$	1.098872
		128	8.681042 × 10${}^{-4}$	1.103865

Table 4

CPU time, speed-up ratio and efficiency of the schemes for Example 1. ($\alpha =0.7, N=100$).

M	201	401	601	801	1001	1201
$T\left(s\right)$	0.123251	0.325116	0.466755	0.879479	1.357376	1.891728
$Tp\left(s\right)$	0.041364	0.106737	0.151507	0.282110	0.421785	0.582479
$Sp$	2.979668	3.045954	3.080749	3.117504	3.218170	3.247719
$Ep$	1.489834	1.522977	1.540374	1.558752	1.609085	1.623860

Table 5

Numerical error and spatial convergence of HASC-N difference scheme for Example 2.

$\mathit{\alpha }$	M	N	${\mathit{E}}_{\mathit{\infty }}\left(\mathit{h},\mathit{n}\right)$	$\mathit{Order}1$
0.5	21	100	1.820529 × 10${}^{-3}$	$--$
	41	400	4.481606 × 10${}^{-4}$	2.095107
	81	1600	1.114969 × 10${}^{-4}$	2.043181
	161	6400	2.783025 × 10${}^{-5}$	2.020325
0.7	21	100	2.197871 × 10${}^{-3}$	$--$
	41	400	5.454330 × 10${}^{-4}$	2.083051
	81	1600	1.359905 × 10${}^{-4}$	2.040008
	161	6400	3.395738 × 10${}^{-5}$	2.019752
0.9	21	100	2.440139 × 10${}^{-3}$	$--$
	41	400	6.093245 × 10${}^{-4}$	2.073776
	81	1600	1.522336 × 10${}^{-4}$	2.036982
	161	6400	3.804131 × 10${}^{-5}$	2.018681

Table 6

Numerical error and temporal convergence of HASC-N difference scheme for Example 2.

$\mathit{\alpha }$	M	N	${\mathit{E}}_{\mathit{\infty }}\left(\mathit{m},\mathit{\tau }\right)$	$\mathit{Order}2$
0.5	101	16	1.939613 × 10${}^{-2}$	$--$
		32	1.355892 × 10${}^{-2}$	0.516527
		64	9.510714 × 10${}^{-3}$	0.511617
		128	6.674831 × 10${}^{-3}$	0.510822
0.7	101	16	1.137173 × 10${}^{-2}$	$--$
		32	6.929028 × 10${}^{-3}$	0.714727
		64	4.230193 × 10${}^{-3}$	0.711930
		128	2.590579 × 10${}^{-3}$	0.707449
0.9	101	16	1.424087 × 10${}^{-2}$	$--$
		32	7.566100 × 10${}^{-3}$	0.912416
		64	4.024197 × 10${}^{-3}$	0.910849
		128	2.154974 × 10${}^{-3}$	0.901031

Table 7

Comparison between exact solution and HASC-N scheme solution for Example 3 ($t=0.5$).

$\mathit{\alpha }$		x
		0.25	0.5	0.75
0.3	Exact solution	6.540829 × 10${}^{-1}$	9.271849 × 10${}^{-1}$	6.816085 × 10${}^{-1}$
	HASC-N scheme solution	6.538316 × 10${}^{-1}$	9.323716 × 10${}^{-1}$	6.823412 × 10${}^{-1}$
0.5	Exact solution	5.807393 × 10${}^{-1}$	8.198668 × 10${}^{-1}$	6.048696 × 10${}^{-1}$
	HASC-N scheme solution	5.844902 × 10${}^{-1}$	8.303472 × 10${}^{-1}$	6.094448 × 10${}^{-1}$
0.7	Exact solution	5.171474 × 10${}^{-1}$	7.268177 × 10${}^{-1}$	5.383339 × 10${}^{-1}$
	HASC-N scheme solution	5.199174 × 10${}^{-1}$	7.356311 × 10${}^{-1}$	5.417396 × 10${}^{-1}$

Table 8

Numerical error and spatial convergence order of HASC-N difference scheme for Example 3.

$\mathit{\alpha }$	M	N	${\mathit{E}}_2\left(\mathit{h},\mathit{n}\right)$	$\mathit{Order}3$
0.1	21	100	4.456840 × 10${}^{-2}$	$--$
	41	400	3.011746 × 10${}^{-2}$	0.565421
	81	1600	2.088088 × 10${}^{-2}$	0.528418
	161	6400	1.448888 × 10${}^{-2}$	0.527236
0.2	21	100	2.061921 × 10${}^{-2}$	$--$
	41	400	8.043848 × 10${}^{-2}$	1.358031
	81	1600	3.662202 × 10${}^{-2}$	1.135174
	161	6400	1.709010 × 10${}^{-3}$	1.099551
0.3	21	100	1.064295 × 10${}^{-2}$	$--$
	41	400	5.183752 × 10${}^{-3}$	1.037830
	81	1600	2.538664 × 10${}^{-3}$	1.029927
	161	6400	1.330512 × 10${}^{-3}$	0.932088
0.4	21	100	1.206877 × 10${}^{-2}$	$--$
	41	400	5.867005 × 10${}^{-3}$	1.040582
	81	1600	2.594804 × 10${}^{-3}$	1.176999
	161	6400	1.241889 × 10${}^{-3}$	1.063090
0.5	21	100	1.505094 × 10${}^{-2}$	$--$
	41	400	6.260401 × 10${}^{-3}$	1.265527
	81	1600	2.542224 × 10${}^{-3}$	1.300164
	161	6400	1.017199 × 10${}^{-3}$	1.321489
0.6	21	100	1.108872 × 10${}^{-2}$	$--$
	41	400	6.254087 × 10${}^{-3}$	0.826221
	81	1600	3.460837 × 10${}^{-3}$	0.853678
	161	6400	2.050544 × 10${}^{-3}$	0.755114
0.7	21	100	1.210763 × 10${}^{-2}$	$--$
	41	400	5.822212 × 10${}^{-3}$	1.056277
	81	1600	2.310618 × 10${}^{-3}$	1.333288
	161	6400	1.009494 × 10${}^{-3}$	1.194646
0.8	21	100	1.428816 × 10${}^{-2}$	$--$
	41	400	5.795162 × 10${}^{-3}$	1.301900
	81	1600	2.306303 × 10${}^{-3}$	1.329267
	161	6400	1.063531 × 10${}^{-3}$	1.116719
0.9	21	100	1.789562 × 10${}^{-2}$	$--$
	41	400	7.830313 × 10${}^{-3}$	1.192465
	81	1600	2.261498 × 10${}^{-3}$	1.791791
	161	6400	9.787323 × 10${}^{-4}$	1.208293

Table 9

Numerical error and temporal convergence order of HASC-N difference scheme for Example 3.

$\mathit{\alpha }$	M	N	${\mathit{E}}_2\left(\mathit{m},\mathit{\tau }\right)$	$\mathit{Order}4$
0.1	101	16	3.235754 × 10${}^{-1}$	$--$
		32	1.322441 × 10${}^{-1}$	1.290898
		64	3.805932 × 10${}^{-2}$	1.796882
		128	1.137961 × 10${}^{-2}$	1.741799
0.2	101	16	3.400465 × 10${}^{-1}$	$--$
		32	1.126743 × 10${}^{-2}$	1.593573
		64	3.902410 × 10${}^{-2}$	1.529722
		128	1.266599 × 10${}^{-2}$	1.623405
0.3	101	16	1.441621 × 10${}^{-1}$	$--$
		32	4.722499 × 10${}^{-2}$	1.610070
		64	1.306512 × 10${}^{-2}$	1.853830
		128	4.114175 × 10${}^{-3}$	1.667045
0.4	101	16	4.485998 × 10${}^{-2}$	$--$
		32	2.509118 × 10${}^{-2}$	0.838249
		64	1.128387 × 10${}^{-2}$	1.152918
		128	3.700395 × 10${}^{-3}$	1.608511
0.5	101	16	8.566342 × 10${}^{-2}$	$--$
		32	4.305125 × 10${}^{-2}$	0.992624
		64	1.101626 × 10${}^{-2}$	1.966421
		128	3.663051 × 10${}^{-3}$	1.588516
0.6	101	16	1.175347 × 10${}^{-1}$	$--$
		32	4.492863 × 10${}^{-2}$	1.387380
		64	1.211080 × 10${}^{-2}$	1.891341
		128	3.668064 × 10${}^{-3}$	1.723204
0.7	101	16	1.384476 × 10${}^{-1}$	$--$
		32	4.957889 × 10${}^{-2}$	1.481543
		64	1.458671 × 10${}^{-2}$	1.765071
		128	4.035364 × 10${}^{-3}$	1.853884
0.8	101	16	1.572810 × 10${}^{-1}$	$--$
		32	5.660295 × 10${}^{-2}$	1.474396
		64	1.851979 × 10${}^{-2}$	1.611810
		128	5.576710 × 10${}^{-3}$	1.731581
0.9	101	16	1.769765 × 10${}^{-1}$	$--$
		32	6.605010 × 10${}^{-2}$	1.421925
		64	2.405686 × 10${}^{-2}$	1.457112
		128	8.758845 × 10${}^{-3}$	1.457636

References

1. Uchaikin, V.V. Fractional Derivatives for Physicists and Engineers, Volume II: Applications; Springer: Berlin, Germany, 2013.

2. Rashid, S.; Hammouch, Z.; Aydi, H.; Ahmad, A.G.; Alsharif, A.M. Novel computations of the time-fractional Fisher’s model via generalized fractional integral operators by means of the Elzaki transform. Fractal Fract.; 2021; 5, 94. [DOI: https://dx.doi.org/10.3390/fractalfract5030094]

3. Li, C.P.; Zeng, F.H. Numerical Methods for Fractional Calculus; CRC Press: New York, NY, USA, 2015.

4. Sabatier, J.; Agrawal, O.P. Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering; Machado, J.A.T. Beijing World Publishing Corporation: Beijing, China, 2014.

5. Diethelm, K. The Analysis of Fraction Differential Equations; Springer: Berlin, Germany, 2010.

6. Liu, F.W.; Zhuang, P.H.; Liu, Q.X. Numerical Methods and its Application of Fractional Partial Differential Equation; Science Press: Beijing, China, 2015; (In Chinese)

7. Sun, Z.Z.; Gao, G.H. Finite Difference Methods for Fractional Differential Equations; 2nd ed. Science Press: Beijing, China, 2021; (In Chinese)

8. Guo, B.L.; Pu, X.K.; Huang, F.H. Fractional Partial Differential Equations and Their Numerical Solutions; Science Press: Beijing, China, 2015.

9. Zhang, X.D.; He, Y.N.; Wei, L.L.; Tang, B.; Wang, S.L. A fully discrete local discontinuous Galerkin method for one-dimensional time-fractional Fisher’s equation. Int. J. Comput. Math.; 2014; 91, pp. 2021-2038. [DOI: https://dx.doi.org/10.1080/00207160.2013.866233]

10. Alquran, M.; Al-Khaled, K.; Sardar, T.; Chattopadhyay, J. Revisited Fisher’s equation in a new outlook: A fractional derivative approach. Phys. A; 2015; 438, pp. 81-93. [DOI: https://dx.doi.org/10.1016/j.physa.2015.06.036]

11. Mejía, C.E.; Piedrahita, A. A numerical method for a time-fractional advection–dispersion equation with a nonlinear source term. J. Appl. Math. Comput.; 2019; 61, pp. 593-609. [DOI: https://dx.doi.org/10.1007/s12190-019-01266-x]

12. Al Qurashi, M.M.; Korpinar, Z.; Baleanu, D. Mustafa Inc. A new iterative algorithm on the time-fractional Fisher equation: Residual power series method. Adv. Mech. Eng.; 2017; 9, pp. 1-8. [DOI: https://dx.doi.org/10.1177/1687814017716009]

13. Khader, M.M.; Saad, K.M. A numerical approach for solving the fractional Fisher equation using Chebyshev spectral collocation method. Chaos Solitons Fractals; 2018; 110, pp. 169-177. [DOI: https://dx.doi.org/10.1016/j.chaos.2018.03.018]

14. Li, D.F.; Liao, H.L.; Sun, W.W.; Wang, J.L.; Zhang, J.W. Analysis of L1-Galerkin FEMs for time-fractional nonlinear parabolic problems. Commun. Comput. Phys.; 2018; 24, pp. 86-103. [DOI: https://dx.doi.org/10.4208/cicp.OA-2017-0080]

15. Veeresha, P.; Prakasha, D.G.; Baskonus, H.M. Novel simulations to the time-fractional Fisher’s equation. Math. Sci.; 2019; 13, pp. 33-42. [DOI: https://dx.doi.org/10.1007/s40096-019-0276-6]

16. Petter, B.; Mitchell, L. Parallel Solution of Partial Differential Equations; Springer: New York, NY, USA, 2000.

17. Pacheco, P. An Introduction to Parallel Programming; Morgan Kaufmann: Burlington, VT, USA, 2011.

18. Gong, C.Y.; Bao, W.M.; Tang, G.J. A parallel algorithm for the Riesz fractional reaction-diffusion equation with explicit finite difference method. Fract. Calc. Appl. Anal.; 2013; 16, pp. 654-669. [DOI: https://dx.doi.org/10.2478/s13540-013-0041-8]

19. Sweilam, N.H.; Moharram, H.; Moniem, N.K.A.; Ahmed, S. A parallel Crank–Nicolson finite difference method for time-fractional parabolic equation. J. Numer. Math.; 2014; 22, pp. 363-382. [DOI: https://dx.doi.org/10.1515/jnma-2014-0016]

20. Biala, T.A.; Khaliq, A.Q.M. Parallel algorithms for nonlinear time-space fractional parabolic PDEs. J. Comput. Phys.; 2018; 375, pp. 135-154. [DOI: https://dx.doi.org/10.1016/j.jcp.2018.08.034]

21. Fu, H.F.; Wang, H. A preconditioned fast parareal finite difference method for space-time fractional partial differential equation. J. Sci. Comput.; 2019; 78, pp. 1724-1743. [DOI: https://dx.doi.org/10.1007/s10915-018-0835-2]

22. Yue, X.Q.; Shu, S.; Xu, X.W.; Bu, W.P.; Pan, K.J. Parallel-in-time multigrid for space-time finite element approximations of two-dimensional space-fractional diffusion equations. Comput. Math. Appl.; 2019; 78, pp. 3471-3484. [DOI: https://dx.doi.org/10.1016/j.camwa.2019.05.017]

23. Liu, H.; Cheng, A.J.; Wang, H. A parareal finite volume method for variable-order time-fractional diffusion equations. J. Sci. Comput.; 2020; 85, 19. [DOI: https://dx.doi.org/10.1007/s10915-020-01321-x]

24. Lorin, E. A parallel algorithm for space-time-fractional partial differential equations. Adv. Differ. Equ.; 2020; 2020, 283. [DOI: https://dx.doi.org/10.1186/s13662-020-02744-4]

25. Wang, Q.L.; Liu, J.; Gong, C.Y.; Tang, X.T.; Fu, G.T.; Xing, Z.C. An efficient parallel algorithm for Caputo fractional reaction-diffusion equation with implicit finite-difference method. Adv. Differ. Equ.; 2016; 2016, 207. [DOI: https://dx.doi.org/10.1186/s13662-016-0929-9]

26. Yang, X.; Wu, L. A New Kind of Parallel Natural Difference Method for Multi-Term Time Fractional Diffusion Model. Mathematics; 2020; 8, 596. [DOI: https://dx.doi.org/10.3390/math8040596]

27. Henry, B.I.; Langlands, T.A.M.; Wearne, S.L. Anomalous diffusion with linear reaction dynamics: From continuous time random walks to fractional reaction-diffusion equations. Phys. Rev. E; 2006; 74, 031116. [DOI: https://dx.doi.org/10.1103/PhysRevE.74.031116]

28. Angstmann, C.N.; Henry, B.I. Time Fractional Fisher–KPP and Fitzhugh–Nagumo Equations. Entropy; 2020; 22, 1035. [DOI: https://dx.doi.org/10.3390/e22091035]

29. Sandev, T.; Tomovski, Z. Fractional Equations and Models: Theory and Applications; Springer: Berlin, Germany, 2019.

30. Ngoc, T.B.; Tri, V.V.; Hammouch, Z.; Can, N.H. Stability of a class of problems for time-space fractional pseudo-parabolic equation with datum measured at terminal time. Appl. Numer. Math.; 2021; 167, pp. 308-329. [DOI: https://dx.doi.org/10.1016/j.apnum.2021.05.009]

31. Dang, D.T.; Nane, E.; Nguyen, D.M.; Tuan, N.H. Continuity of Solutions of a Class of Fractional Equations. Potential Anal; 2018; 49, pp. 423-478. [DOI: https://dx.doi.org/10.1007/s11118-017-9663-5]

32. Zhang, Y.N.; Sun, Z.Z. Alternating direction implicit schemes for the two-dimensional fractional sub-diffusion equation. J. Comput. Phys.; 2011; 230, pp. 8713-8728. [DOI: https://dx.doi.org/10.1016/j.jcp.2011.08.020]

33. Liu, N.; Liu, Y.; Li, H.; Wang, J. Time second-order finite difference/finite element algorithm for nonlinear time-fractional diffusion problem with fourth-order derivative term. Comput. Math. Appl.; 2018; 75, pp. 3521-3536. [DOI: https://dx.doi.org/10.1016/j.camwa.2018.02.014]

34. Wang, Y.; Liu, Y.; Li, H.; Wang, J. Finite element method combined with second-order time discrete scheme for nonlinear fractional cable equation. Eur. Phys. J. Plus; 2016; 131, 61. [DOI: https://dx.doi.org/10.1140/epjp/i2016-16061-3]

35. Zhang, B.L.; Gu, T.X.; Mo, Z.Y. Principles and Methods of Numerical Parallel Computation; National Defence Industry Press: Beijing, China, 1999; (In Chinese)

36. Zhou, Y.L. Difference schemes with intrinsic parallelism for quasi-linear parabolic systems. Sci. China Math.; 1997; 40, pp. 270-278. [DOI: https://dx.doi.org/10.1007/BF02874519]

37. Stynes, M.; O’Riordan, E.; Gracia, J.L. Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. SIAM J. Numer. Anal.; 2017; 55, pp. 1057-1079. [DOI: https://dx.doi.org/10.1137/16M1082329]

38. Shen, J.Y.; Sun, Z.Z.; Du, R. Fast finite difference schemes for time-fractional diffusion equations with a weak singularity at initial time. East Asian J. Appl. Math.; 2018; 8, pp. 834-858. [DOI: https://dx.doi.org/10.4208/eajam.010418.020718]

39. Yang, X.Z.; Liu, X.L. Numerical analysis of two new finite difference methods for time-fractional telegraph equation. Discrete Contin. Dyn. Syst. Ser. B; 2021; 26, pp. 3921-3942. [DOI: https://dx.doi.org/10.3934/dcdsb.2020269]

40. Evans, G.; Blackledge, J.; Yardley, P. Numerical Methods for Partial Differential Equations; Springer: Berlin, Germany, 2012.

41. Kumar, D.; Chaudhary, S.; Srinivas Kumar, V.V.K. Fractional Crank–Nicolson–Galerkin finite element scheme for the time-fractional nonlinear diffusion equation. Numer. Methods Partial Differ. Equ.; 2019; 35, pp. 2056-2075. [DOI: https://dx.doi.org/10.1002/num.22399]

42. Vong, S.; Lyu, P.; Wang, Z.B. A compact difference scheme for fractional sub-diffusion equations with the spatially variable coefficient under neumann boundary conditions. J. Sci. Comput.; 2016; 66, pp. 725-739. [DOI: https://dx.doi.org/10.1007/s10915-015-0040-5]

43. Tverdyi, D.; Parovik, R. Investigation of Finite-Difference Schemes for the Numerical Solution of a Fractional Nonlinear Equation. Fractal Fract.; 2022; 6, 23. [DOI: https://dx.doi.org/10.3390/fractalfract6010023]

44. Long, L.D.; Luc, N.H.; Tatar, S.; Baleanu, D.; Can, N.H. An inverse source problem for pseudo-parabolic equation with Caputo derivative. J. Appl. Math. Comput.; 2022; 68, pp. 739-765. [DOI: https://dx.doi.org/10.1007/s12190-021-01546-5]

45. Yuan, G.W.; Sheng, Z.Q.; Hang, X. The unconditional stability of parallel difference schemes with second order convergence for nonlinear parabolic system. J. Partial Differ. Equ.; 2007; 20, pp. 45-64.

46. Deng, W.H.; Zhang, Z.J. High Accuracy Algorithms for the Differential Equation Governing Anomalous Diffusion, Algorithm and Models for Anomalous Diffusion; World Scientific: Singapore, 2019.