An Efficient Numerical Scheme for Variable-Order Fractional Sub-Diffusion Equation

Abstract

The variable-order (VO) fractional calculus can be seen as a natural extension of the constant-order, which can be utilized in physical and biological applications. In this study, we derive a new numerical approximation for the VO fractional Riemann–Liouville integral formula and developed an implicit difference scheme (IDS) for the variable-order fractional sub-diffusion equation (VO-FSDE). The derived approximation used in the VO time fractional derivative with the central difference approximation for the space derivative. Investigated the unconditional stability by the van Neumann method, consistency, and convergence analysis of the proposed scheme. Finally, a numerical example is presented to verify the theoretical analysis and effectiveness of the proposed scheme.

1. Introduction

Fractional differential equations (FDEs) got the researchers attention and this fact states the phenomena in various areas such as biochemistry, wave propagation, medicine, anomalous diffusion, and electrical engineering, etc. [1,2,3,4]. The fractional diffusion equation is an important FDE, which has been hugely used in fractional diffusive systems, diffusion in transport processes, unification of diffusion phenomena and fractional random walk [5,6]. The fractional sub-diffusion equation (FSDE) is a subclass of the fractional diffusion equation and obtained by replacing the time derivative by fractional-orderαlying between 0 and 1. The one-dimensional FSDE can be written as:

∂u(x,t)∂t=^∂1−α∂^t1−α∂2u(x,t)∂^x2+f(x,t).

where^∂1−α∂^t1−αrepresent the Riemann-Liouville fractional derivative (RL-FD) of order1−α.

Most of the FDEs are not easy even impossible to get an analytical solution, which contains complicated functions. Therefore, numerical methods have become the main way for the solution of such types of fractional order differential equations. Many researchers have constructed the various numerical methods to solve FDEs such as Yuste [7] considered the weighted average implicit difference scheme (IDS) for the solution of FSDE. They discussed the theoretical analysis by the well known Fourier series method. Ali et al. [8] constructed a new numerical scheme for the modified FSDE. They investigated the theoretical analysis and found that the scheme is unconditionally stable and highly accurate. In another survey Ali et al. [9] derived a new approximation for RL-FD and developed compact IDS for modified FSDE. They successfully discussed the theoretical analysis and accuracy. Mohebbi et al. [10] studied the compact order IDS and solved the modified FSDE. They approximated the fractional derivative by Grünwald–Letnikov and the space derivative by fourth order approximation. The theoretical analysis are discussed by Fourier method and the tested examples shown high accuracy. Zhuang et al. [11] formulated the physical-mathematical method for FSDE and investigated a new numerical scheme. They analyzed the theoretical analysis based on energy method and the numerical values fully support the theoretical analysis. Ali et al. [12] proposed the Crank–Nicolson method for 2D-FSDE. They replaced the fractional derivative with Grünwald–Letnikov approximation and the space derivatives with finite difference approximation. They discussed the theoretical analysis with unconditional stability and convergence. Khan and Rasheed [13] studied the mixed convection in an incompressible fractional Maxwell nanofluid with thermal transmission. They tackled the arising mathematical problem numerically via finite difference procedure in their reported study they have mentioned that escalating values of Prandtl number reduces the thermal transport. Hamid et al. [14] investigated the modeling of unsteady radiating fractional nanofluid problem in a channel. Fractional derivative in Caputo sense is engaged in the model problem. They computed the solution of the boundary layer equations with the help of Crank–Nicolson numerical scheme. They found that the magnetic parameter reduces the velocity, whereas, Grashof number boosts fluid velocity. Other related literature can be seen in [15,16,17,18]. The above cited studies are reported the fractional derivative of constant-order.

The constant-order fractional derivatives are sometimes not more suitable to explain the complex diffusion processes as in porous medium and medium structure because it changes with time [19,20]. In VO operator the order may vary as function of independent variables such as time or space or both time-space. Lorenzo and Hartly [21] discussed some physical phenomena which show that the fractional-order behavior may change with time and space. To solve such type of VO fractional model and developed some numerical techniques. Here, we include only the VO-FSDE, which can be written as:

∂u(x,t)∂t=^{∂1−α(x,t)}∂^{t1−α(x,t)∂2}u(x,t)∂^x2+f(x,t).

There are many researchers who have worked on a VO-FSDE and find the solution by different numerical techniques such as Chen et al [22] developed a numerical scheme for VO-FSDE with first-order temporal and fourth-order spatial accuracy. they investigated the stability, convergence, and solvability analysis by Fourier series method. Lin et al [23] consider the nonlinear VO-FSDE and construct a new numerical finite difference scheme, also discussed stability and convergence analysis. Sun et al. [24] discussed different complex diffusion processes that can be described by the VO model. They differentiate this model into four types such as time-dependent, space-dependent, concentration-dependent, and system parameter dependent models. Sweilam et al. [25] solved the VO Riesz space fractional wave equation by explicit difference method. They discussed the Riesz space fractional derivatives based on shifted Grünwald–Letnikov formula and also established the stability and convergence analysis. Chen et al. [26] formulated the numerical solution for the nonlinear VO-FSDE. They proposed the piecewise function for VOFD and then by operational matrices transformed into an algebraic equation. Bhrawy and Zaky [27] investigated a numerical method based on the collocation method combined with Jacobi operational matrix and successfully solved the one and 2D-VO fractional Cable equations. They analyzed the convergence analysis and the numerical results demonstrated that the algorithm is powerful with high accuracy. Chen [28] solved the 2D-VO modified FSDE by numerical method. They discussed the stability, convergence, and solvability by the Fourier method. The two new second order numerical approximations are derived for VO time-fractional operator by Zhao et al. [29]. They tested on VO-FSDE and super-diffusion problems, the results are different than standard diffusion and constant-order diffusion for the VOα(x,t) . Wang et al. [30] considered the VO-FDE with variable coefficient. They solved the problem by finite difference method and homotopy regularization method based on Legendre polynomials as the basis functions of approximation space. Ma et al. [31] discussed the numerical Adams-Bashforth-Moulton method for the VO fractional financial system of equations. The tested examples reported that the proposed method for VO-FDE is simple and effective. Shekari et al. [32] considered the moving least squares method to find the solution of 2D VO fractional diffusion-wave equation. They utilized the mesh free and collocation method based on moving least square approximation. Xu et al. [33] solved the multi-term VO space-time fractional diffusion equation on a finite domain by numerical method. Investigated the theoretical analysis via mathematical induction and numerical examples confirmed the efficiency and accuracy. Shen et al. [34] considered the Coimbra VO time derivative which is a Caputo type definition and used in the numerical scheme for variable -order fractional diffusion equation. The theoretical analyses are discussed by the well know Fourier series. Further, other related studies of fractional and VO-FDE can see in [35,36,37,38,39,40].

In this article, we developed a new implicit difference method to solve the VO-FSDE by applying the new discretized approximation for VO-RL fractional integral. The first order time derivative is replaced by backward difference approximation. Further, investigated the stability analysis by van Neumann method, consistency and convergence of the proposed method. The proposed scheme is very easy to implement, it decreases the computational complexity and increase efficiency. Another positive point of this scheme is, the easiness in finding the stability and convergence analysis.

This investigation is arranged as follows: Section 1, contains the comprehensive literature survey relevant to the reported material. In Section 2, briefly explain the mathematical preliminaries, we discussed the implicit scheme in Section 3. In Section 3.1 and Section 3.2, we use the van Neumann method to prove the stability and consistency respectively. In Section 4, provided the numerical experiments to confirm the theoretical analysis and in the last section concluded the conclusion.

2. Mathematical Preliminaries

In this section, we discuss the VO fractional calculus [41] and deriving new approximation for VO-RL fractional integral, which is used in this paper:

Definition 1.

The VO RL-FD formula can be written as:

_{0RL Dt1−α(x,t)}u(x,t)=1Γ(α(x,t))∂∂t_∫0tu(x,η)^{(t−η)1−α(x,t)}dη.

Definition 2.

The VO-RL fractional integral formula can be written as:

_I0α(x,t)u(x,t)=1Γα(x,t)_∫0tu(x,η)^{(t−η)1−α(x,t)}dη.

Here, we formulate a new numerical approximation for VO-RL fractional integral of orderα(x,t)and(0<α(x,t)<1)at the grid point(_xi,_tk), as following:

_I0α(xi,tk)u(_xi,_tk)=1Γα(_xi,_tk)_∫0tk u(_xi,ξ)^{(_tk−ξ)1−α(_xi,_tk)}dξ,=1Γα(_xi,_tk)_∫0tk ^{(_tk−ξ)α(_xi,_tk)−1}u(_xi,ξ)dξ,=∑j=0k−1_∫tjtj+1 ^{(_tk−ξ)α(_xi,_tk)−1}u(_xi,ξ)Γα(_xi,_tk)dξ,=∑j=0k−1u(_xi,_tk−j)Γα(_xi,_tk)_∫tjtj+1 ^{(_tk−ξ)α(_xi,_tk)−1}dξ,=∑j=0k−1u(_xi,_tk−j)Γα(_xi,_tk)−^{(_tk−ξ)α(_xi,_tk)}α(_xi,_tk)_|tjtj+1 ,=∑j=0k−1^{τα(_xi,_tk)}u(_xi,_tk−j)Γ(1+α(_xi,_tk))^{(j+1)α(_xi,_tk)}−^{(j)α(_xi,_tk)},=∑j=0k−1^{τα(_xi,_tk)}Γ(1+α(_xi,_tk))(^{(j+1)α(_xi,_tk)}−^{(j)α(_xi,_tk)})u(_xi,_tk−j).

Lemma 1.

The VOα(x,t)(0<α(x,t)<1)fractional RL integral of the functionu(x,t)on[0,t]can be defined in discretized form at the grid point(_xi,_tk)as:

_I0α(xi,tk)u(_xi,_tk)=_Cik∑j=0k−1_{bj(αik) uik−j},

where_Cik=^τ_αik Γ(1+_αik),_bj=(^(j+1)_αik −^(j)_αik ).

Lemma 2.

The coefficients_bj(αik)(j=0,1,2,…) satisfy the following properties [39]:

(i)

_b0(αik)=1,_bj(αik)>0,j=0,1,2,….

(ii)

_bj−1(αik)>_bj(αik),j=1,2,….

(iii)

There exists a positive constantC>0,such thatτ≤C_bj(αik) ^τ_αik ,j=1,2,….

(iv)

_{∑j=0k bj(αik)} ^τ_αik =^(k+1)_αik ≤^T_αik .

3. The Implicit Scheme

To construct an IDS for VO-FSDE (2), first we use Equation (4), then applying lemma 1 for the discretization of RL fractional integral and for space derivative utilizing the central difference approximation. Here, the x variable steps as_xi=iΔx, in the x-direction withi=1,…,M−1,Δx=LMand for t variable the step is_tk=kτ,k=1,…,Nwhereτ=TN. Let_uikbe the numerical approximation tou(_xi,_tk) , put Equation (4) into Equation (2), we have

∂u(_xi,_tk)∂t=∂∂t_I0α(xi,tk)^∂2u(_xi,_tk)∂^x2+f(_xi,_tk),

by applying lemma 1 and backward difference approximation to Equation (7), we get

_uik−_uik−1=^{τ_αik (Δx)2}Γ(1+_αik)∑j=0k−1_{bj(αik) δx2}(_uik−j−_uik−j−1)+τ_fik.

Simplifying Equation (8), we get a new IDS for the one-dimensional VO-FSDE (2) with the conditions, as follows

−_{Cik ui+1k}+(1+2_Cik)_uik−_{Cik ui−1k}=_uik−1−_{Cik bk−1(αik)}(_ui+10−2_ui0+_ui−10)+_Cik∑j=1k−1(_bj(αik)−_bj−1(αik))(_ui+1k−j−2_uik−j+_ui−1k−j)+τ_fik.

Here,

_Cj=^{τ_αik (Δx)2}Γ(1+_αik),_{δx2 uik}=_ui+1k−2_uik+_ui−1k,

andi=1,2,…,_Mx−1, andk=1,2,…,N−1,with

_ui0=φ(_xi),

_u0k=_φ1(_tk),_uM,k=_φ2(_tk),0≤x≤L,0≤t≤T.

3.1. Stability

To analyze the stability of the proposed IDS using van Neumann method. Let_Uik represent the approximated solution for (9), we obtained

−_{Cik Ui+1k}+(1+2_Cik)_Uik−_{Cik Ui−1k}=_Uik−1−_{Cik bk−1(αik)}(_Ui+10−2_Ui0+_Ui−10)+_Cik∑j=1k−1(_bj(αik)−_bj−1(αik))(_Ui+1k−j−2_Uik−j+_Ui−1k−j)+τ_fik.

The error is defined as

_eik=_uik−_Uik,

where_eik satisfies Equation (13), as

−_{Cik ei+1k}+(1+2_Cik)_eik−_{Cik ei−1k}=_eik−1−_{Cik bk−1(αik)}(_ei+10−2_ei0+_ei−10)+_Cik∑j=1k−1(_bj(αik)−_bj−1(αik))(_ei+1k−j−2_eik−j+_ei−1k−j).

The error initial and boundary conditions are given by

_e0k=_eMk=_ei0=0.

By defining the following grid functions fork=1,2,…,N.

^ek(x)=_eik,when_xi−Δx2<x≤_xi+Δx2,0,when0≤x≤Δx2orL−Δx2≤x≤L,

then^ek(x)can be expanded in Fourier series such as:

^ek(x)=∑_l1=−∞∞^λk(_l1)^{e2−1π(_l1x/L)},

where

^λk(_l1)=1L_{∫0L ∫0L} ^ek(x)^{e−2−1π(_l1x/L)}dx.

From the definition of^l2norm and Parseval equality, we have

‖^ek _‖∞2=∑i=1_Mx−1Δx|_eik ^|2=∑_l1∞^|λk(_l1)|2.

Supposing that

_eik=^{λk e−1(_σ1iΔx)},

here,_σ1=2π_l1/L and substituting (21) in (15), we obtain

^λk=1(1+_μik)^λk−1+_{μik bk−1(αik)} ^λ0−_μik∑j=1k−1(_bj−1(αik)−_bj(αik))^λk−j,

where

_μik=4_Ciksiⁿ²(_σ1Δx2).

Proposition 1.

If^λk(k=1,2,…,N) satisfy (22), then|^λk|≤|^λ0|.

Proof.

By using mathematical induction, we takek=1 in (22)

^λ1=(1+_{μi1 b0(αi1)})^λ0(1+_μi1),

and as_μi1≥0and_b0(αi1)=1, then

|^λ1|≤|^λ0|.

Now, assuming that

|^λm|≤|^λ0|;m=1,2,…,k−1.

As0<_αik<1 , from (22) and lemma 1, we have

|^λk|≤|^λk−1|+_{μik bk−1(αik)}|^λ0|+_{μik ∑j=1k−1}(_bj−1(αik)−_bj(αik))|^λk−j|(1+_μik),≤1+_{μik bk−1(αik)}+_{μik ∑j=1k−1}(_bj−1(αik)−_bj(αik))(1+_μik)|^λ0|,=1+_{μik bk−1(αik)}+_μik(_b0(αi1)−_bk−1(αik))(1+_μik)|^λ0|,=1+_μik1+_μik|^λ0|,|^λk|≤|^λ0|.

From the above proof, it is clear that Proposition 1 and Equation (20), reported that the solution of Equation (9) satisfies

‖^λk _‖2≤_‖^λ0‖2,

this proved that the IDS in (9) is unconditionally stable. ☐

3.2. Consistency

To analyse the consistency of the IDS, let U be the exact and u be the approximate solution andF(U)=0 represent the approximated difference equation of the FDE (2) at point(_xi,_tk), thenF(u)=_Tikrepresented the truncation error at grid point(_xi,_tk).

Theorem 1.

The truncation errorT(x,t)of the proposed difference scheme is:_Tik=O(Δt)+O(Δ^x2).

Proof.

_Tik=_uik−_uik−1−_Cik∑j=0k−1_bj(αik)(_ui+1k−j−2_uik−j+_ui−1k−j)−(_{ui+1k−j−1}−2_uik−j−1+_{ui−1k−j−1}).

By Taylor expansion we can write

_Tik=_uik−_uik−(Δt)∂u∂t_|ik+^(Δt)22^∂2u∂^t2_|ik+…−_Cik∑j=0k−1_bj(αik)[(_uik+(Δx)∂u∂x_|ik−j+^(Δx)22^∂2u∂^x2_|ik−j+^(Δx)36∂^u3∂^x3_|ik−j+…−2_uik−j+_uik−(Δx)∂u∂x_|ik−j+^(Δx)22^∂2u∂^x2_|ik−j−^(Δx)36^∂3u∂^x3_|ik−j+…−(_uik−j−1+(Δx)∂u∂x_|ik−j−1+^(Δx)22^∂2u∂^x2_|ik−j−1+^(Δx)36^∂3u∂^x3_|ik−j−1+…−2_uik−j−1+_uik−j−1−(Δx)∂u∂x_|ik−j−1+^(Δx)22^∂2u∂^x2_|ik−j−1−^(Δx)36^∂3u∂^x3_|ik−j−1+…)],

and

_Tik=−(Δt)∂u∂t_|ik+^(Δt)22^∂2u∂^t2_|ik+…−_Cik∑j=0k−1_bj(αik)[^(Δx)2∂2u∂^x2_|ik−j−^∂2u∂^x2_|ik−j−1−^(Δx)412^∂4u∂^x2_|ik−j+…−^∂4u∂^x4_|ik−j−1+…].

After simplification, we get the following truncation error

_Tik=O(Δt)+O^(Δx)2.

This theorem shows that the proposed scheme is consistent ifΔx→ 0 andΔt→ 0, then the truncation error tends to zero. ☐

Theorem 2.

According to Lax equivalence theorem, consistency and stability are both necessary and sufficient for convergence [39], the proposed VO-IDS is convergent.

4. Numerical Experiments

The numerical example of VO-FSDE are presented in this section to confirm the effectiveness of the new IDS. Here, we compute the_E∞and_E2norm error, is defined as follows:

_E∞=max0≤i≤_Mx−1,0≤k≤N|u(_xi,_tk)−_uik|,_E2=^{(Δx)∑i=1_Mx−1|u(_xi,_tk)−_uik|12}.

Example 1.

Consider the VO-FSDE is as following [19]:

∂u(x,t)∂t=^{∂1−α(x,t)}∂^{t1−α(x,t)∂2}u(x,t)∂^x2+2^ext−^t1+α(x,t)Γ(2+α(x,t)),0≤t≤T,0<x<1,

the initial and the boundary conditions are

u(x,0)=0,0≤x≤1,u(0,t)=^t2,u(1,t)=e^t2,0≤x≤L,0≤t≤T.

The exact solution is

u(x,t)=^{ex t2}.

Table 1 and Table 2 shows the numerical results of the new VO IDS for Equation (26). In Figure 1, the numerical values are compared with the previous study in Cao et al. [19] for various values of time steps and fixed space step. The proposed VO scheme reported better level of accuracy and more efficiency. Table 2 show that the error is reduced as the number of space and time steps increased. Figure 1 and Figure 2 are plotted for numerical scheme and compared it with exact solution for different values ofα(x,t)=^2−xt−sin(xt)50,2+sin(xt)400, andy=0.25,0.083andT=1.0respectively. The above discussion proved that new IDS for VO-FDEs is accurate and numerically more efficient.

5. Conclusions

In this article, first we derived a new approximation for VO-RL fractional integral and for space derivative used the implicit central difference approximation. Successfully utilized the new numerical approximation and the IDS for one-dimensional VO-FSDE. Investigated the theoretical analysis by the well-known von Neumann method. The scheme is unconditionally stable and convergent with orderO(τ+^(Δx)2). The tested example result has compared with the previous study and with the exact solution for various values ofα(x,t). We concluded that the scheme is a more effective, high level of accuracy and the solution is well-matched. We believe that this finding is another contribution to the literature. This numerical approximation can also be applied to other types of variable-order fractional differential equations.

Figure 1. The comparison of IDS for (26) with (28) atα(x,t)=2xt-sin(xt)50,T=1,y=0.25andN=4.

Figure 2. The comparison of IDS for (26) with (28) atα(x,t)=2+sin(xt)400,T=1,y=0.083andN=256.

	Error E∞		Error E2
	τ	[19]	Method	[19]	Method
α(x,t)=e−xt400	1/4	6.0309×10−3	2.1024×10−4	4.4044×10−2	1.5379×10−4
	1/8	1.5166×10−3	6.2311×10−5	1.1078×10−2	4.3215×10−5
	1/16	3.7704×10−4	5.2566×10−5	2.7545×10−3	3.3345×10−5
	1/32	9.3918×10−5	1.1145×10−5	6.8622×10−4	7.7436×10−5
α(x,t)=2+sin(xt)400	1/4	6.0281×10−3	7.9047×10−4	4.4023×10−2	5.7621×10−4
	1/8	1.5171×10−3	5.7892×10−4	1.1081×10−2	4.2120×10−4
	1/16	3.7741×10−4	1.3667×10−4	2.7572×10−3	3.1659×10−4
	1/32	9.4057×10−5	3.4864×10−5	6.8721×10−4	2.5194×10−4
α(x,t)=2xt+1300	1/4	6.0318×10−3	8.1499×10−4	4.4051×10−2	5.9070×10−4
	1/8	1.5182×10−3	7.4512×10−4	1.1089×10−2	5.3748×10−4
	1/16	3.7704×10−4	5.2566×10−5	2.7545×10−3	3.3345×10−5
	1/32	9.4079×10−5	4.5761×10−5	6.8736×10−4	2.7098×10−5

α(x,t)	τ=Δx=10	τ=Δx=20	τ=Δx=40	τ=Δx=80
2cos(xt)+(xt)56	2.6785×10−3	1.5233×10−3	8.6046×10−4	6.1524×10−4
1−(xt)3+cos2(xt)10	2.2091×10−3	1.4258×10−3	8.4040×10−4	4.3821×10−4
30+(xt)4750	2.0303×10−3	1.1525×10−3	6.5034×10−4	3.6267×10−4
15+(sin(xt))8400	1.9153×10−3	1.0868×10−3	6.1368×10−4	2.8031×10−4
10−(xt)4300	1.7106×10−3	9.5885×10−4	5.3223×10−4	2.9679×10−4
2+sin(xt)400	5.3892×10−4	4.3962×10−4	1.4306×10−4	7.7335×10−5
2(xt)−sin(xt)50	8.9348×10−4	3.8637×10−4	1.0794×10−4	5.7894×10−5

Author Contributions

Conceptualization, U.A.; methodology, U.A., and M.S.; software, U.A.; formal analysis, U.A., M.S., and F.A.A.; investigation, U.A.; writing-original draft preparation, U.A., and M.S.; writing-review and editing, U.A., M.S., and F.A.A.; Funding F.A.A. All authors have read and agreed to the published version of the manuscript.

Funding

The research is supported by Publication (Article In Journal) fee fund by USM's Research Creativity and Management Office (RCMO) and School of Mathematical Sciences KPI grant.

Acknowledgments

The authors are thankful to the reviewers for the useful suggestions to improve the quality of this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

List of Abbreviation/Nomenclature

IDS	Implicit difference scheme
FSDE	Fractional sub-diffusion equation
FDEs	Frcational differential rquations
VO	Variable-order
2D	Two-dimensional
RL-FD	Riemann-Liouville fractional derivative

References

1. Li, X.; Wong, P.J. A new implicit numerical scheme for fractional sub-diffusion equation. In Proceedings of the 2016 14th International Conference on Control, Automation, Robotics and Vision (ICARCV), Phuket, Thailand, 13-15 November 2016; pp. 1-6.

2. Zhang, P.; Pu, H. A second-order compact difference scheme for the fourth-order fractional sub-diffusion equation. Numer. Algorithms 2017, 76, 573-598.

3. Oldham, K.B.; Spanier, J. Fractional Calculus; Academic Press: New York, NY, USA; London, UK, 1974.

4. Podlubny, I. Fractional Differential Equations; Academic Press: New York, NY, USA, 1999.

5. Zhang, J.; Ye, C. High order numerical method and its analysis of the anomalous subdiffusion equation. Procedia Eng. 2012, 31, 781-790.

6. Zhang, Y.N.; Sun, Z.Z. Alternating direction implicit schemes for the two-dimensional fractional sub-diffusion equation. J. Comput. Phys. 2011, 230, 8713-8728.

7. Yuste, S.B. Weighted average finite difference methods for fractional diffusion equations. J. Comput. Phys. 2006, 216, 264-274.

8. Ali, U.; Abdullah, F.A.; Mohyud-Din, S.T. Modified implicit fractional difference scheme for 2D modified anomalous fractional sub-diffusion equation. Adv. Diff. Equ. 2017, 185, 1-14.

9. Ali, U.; Sohail, M.; Usman, M.; Abdullah, F.A.; Khan, I.; Nisar, K.S. Fourth-Order Difference Approximation for Time-Fractional Modified Sub-Diffusion Equation. Syemmetry 2010, 12, 691.

10. Mohebbi, A.; Abbaszadeh, M.; Dehghan, M. A high-order and unconditionally stable scheme for the modified anomalous fractional sub-diffusion equation with a nonlinear source term. J. Comp. Phys. 2013, 240, 36-48.

11. Zhuang, P.; Liu, F.; Anh, V.; Turne, I. New solution and analytical techniques of the implicit numerical methods for the anomalous sub-diffusion equation. SIAM J. Numer. Anal. 2008, 46, 1079-1095.

12. Ali, U.; Abdullah, F.A.; Ismail, A.I. Crank-Nicolson finite difference method for two-dimensional fractional sub-diffusion equation. J. Interpolat. Approx. Sci. Comput. 2017, 2017, 18-29.

13. Khan, A.Q.; Rasheed, A. Mixed convection magnetohydrodynamics flow of a nanofluid with heat transfer: A numerical study. Math. Probl. Eng. 2019, 2019, 8129564.

14. Hamid, M.; Zubair, T.; Usman, M.; Haq, R.U. Numerical investigation of fractional-order unsteady natural convective radiating flow of nanofluid in a vertical channel. AIMS Math. 2019, 4, 1416.

15. Ali, U.; Abdullah, F.A. Explicit Saul'yev finite difference approximation for two- dimensional fractional sub-diffusion equation. In Proceedings of the AIP Conference Proceedings, Pahang, Malaysia, 27-29 August 2018; Volume 1, p. 020111.

16. Liangliang, M.A.; Dongbing, L. An implicit difference approximation for fractional cable equation in high-dimensional case. J. Liao. Tech. University (Nat. Sci.) 2014, 4, 024.

17. Zhai, S.; Feng, X.; He, Y. An unconditionally stable compact ADI method for 3D time-fractional convection-diffusion equation. J. Comp. Phys. 2014, 269, 138-155.

18. Zhuang, P.; Liu, F.; Turner, L.; Anh, V. Galerkin finite element method and error analysis for the fractional cable equation. Numer. Algor. 2016, 72, 447-466.

19. Cao, J.; Qiu, Y.; Song, G. A compact finite difference scheme for variable order subdiffusion equation. Commun. Nonlinear Sci. Numer. Simul. 2017, 48, 140-149.

20. Sun, H.; Chen, W.; Li, C.; Chen, Y. Finite difference schemes for variable-order time fractional diffusion equation. Int. J. Bifurc. Chaos 2012, 22, 1250085.

21. Lorenzo, C.F.; Hartley, T.T. Variable order and distributed order fractional operators. Nonlinear Dyn. 2002, 29, 57-98.

22. Chen, C.M.; Liu, F.; Anh, V.; Turner, I. Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation. SIAM J. Sci. Comput. 2010, 32, 1740-1760.

23. Lin, R.; Liu, F.; Anh, V.; Turner, I. Stability and convergence of a new explicit finite-difference approximation for the variable-order nonlinear fractional diffusion equation. Appl. Math. Comput. 2009, 212, 435-445.

24. Sun, H.; Chen, W.; Chen, Y. Variable-order fractional differential operators in anomalous diffusion modeling. Phy. A Stat. Mech. Its Appl. 2009, 388, 4586-4592.

25. Sweilam, N.; Khader, M.; Almarwm, H. Numerical studies for the variable-order nonlinear fractional wave equation. Fract. Calc. Appl. Anal. 2012, 15, 669-683.

26. Chen, C.M.; Liu, F.; Turner, I.; Anh, V. Numerical methods with fourth-order spatial accuracy for variable-order nonlinear Stokes' first problem for a heated generalized second grade fluid. Comput. Math. Appl. 2011, 62, 971-986.

27. Bhrawy, A.H.; Zaky, M.A. Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation. Nonlinear Dyn. 2015, 80, 101-116.

28. Chen, C.M. Numerical methods for solving a two-dimensional variable-order modified diffusion equation. Appl. Math. Comput. 2013, 225, 62-78.

29. Zhao, X.; Sun, Z.Z.; Karniadakis, G.E. Second-order approximations for variable order fractional derivatives: Algorithms and applications. J. Comput. Phys. 2015, 293, 184-200.

30. Wang, S.; Wang, Z.; Li, G.; Wang, Y. A simultaneous inversion problem for the variable-order time fractional differential equation with variable coefficient. Math. Probl. Eng. 2019, 2019, 2562580.

31. Ma, S.; Xu, Y.; Yue, W. Numerical solutions of a variable-order fractional nancial system. J. Appl. Math. 2012, 2012, 417942.

32. Shekari, Y.; Tayebi, A.; Heydari, M.H. A meshfree approach for solving 2D variable-order fractional nonlinear diffusion-wave equation. Comput. Methods Appl. Mech. Eng. 2019, 350, 154-168.

33. Xu, T.; Lü, S.; Chen, W.; Chen, H. Finite difference scheme for multi-term variable-order fractional diffusion equation. Adv. Diff. Equ. 2018, 2018, 103.

34. Shen, S.; Liu, F.; Chen, J.; Turner, I.; Anh, V. Numerical techniques for the variable order time fractional diffusion equation. Appl. Math. Comput. 2012, 218, 10861-10870.

35. Ali, U.; Abdullah, F.A. Modified implicit difference method for one-dimensional fractional wave equation. In Proceedings of the AIP Conference Proceedings, Penang, Malaysia, 10-12 December 2019; Volume 2184, p. 0060021.

36. Bhrawy, A.H.; Zaky, M. Numerical algorithm for the variable-order Caputo fractional functional differential equation. Nonlinear Dyn. 2016, 85, 1815-1823.

37. Chen, Y.M.; Wei, Y.Q.; Liu, D.Y.; Yu, H. Numerical solution for a class of nonlinear variable order fractional differential equations with Legendre wavelets. Appl. Math. Lett. 2015, 46, 83-88.

38. Yaghoobi, S.; Moghaddam, B.P.; Ivaz, K. An efficient cubic spline approximation for variable-order fractional differential equations with time delay. Nonlinear Dyn. 2017, 87, 815-826.

39. Ali, U. Numerical Solutions for Two Dimensional Time-Fractional Differential Sub-Diffusion Equation. Ph.D. Thesis, University Sains Malaysia, Penang, Malaysia, 2019; pp. 1-200.

40. Samko, S.G. Fractional integration and differentiation of variable order. Anal. Math. 1995, 21, 213-236.

41. Samko, S.G.; Kilbas, A.A.; Marichev, O.I. Fractional Integrals and Derivaives, Theory and Applications; Gordon and Breach Science Publishers: New York, NY, USA, 1993.