Application of a Legendre collocation method to the space–time variable fractional-order advection–dispersion equation

ABSTRACT

In this paper, an efficient spectral collocation method is presented for solving the space–time variable-order fractional advection–dispersion equations (ST-VO-FADE). The proposed method is based on the Legendre collocation spectral procedure together with the Legendre operational matrices for fractional derivatives, described in the sense of Riemann–Liouville and Caputo. The main characteristic behind this approach is to reduce such problems to those of solving systems of algebraic equations in the unknown expansion coefficients of the sought-for spectral approximations, which greatly simplifies the solution process. The validity of the method is demonstrated by solving two numerical examples. Finally, comparisons between the algorithm derived in this paper and the existing algorithms are given, which show that our numerical schemes exhibit better performances than the existing ones.

KEYWORDS:

• Variable-order derivative
• Space–time fractional advection–dispersion equation
• Collocation method
• Legendre polynomials

1. Introduction

Fractional derivative operators are found to be more real in modelling a variety of engineering processes, physical behaviours, biological models and financial applications, such as viscoelastic materials, anomalous diffusion and non-exponential relaxation patterns, among others (see, e.g. [1–9] and the references therein). Regarding their importance, there are difficulties in getting the exact solutions of fractional differential equations (FDEs) due to the property of non-locality of the fractional derivative operators. Hence, the numerical methods are important tools to understand the physical behaviour of these equations. A great deal of work is done on various numerical methods to solve space or time constant-order fractional partial differential equations during the past decade [10–22]. Introduction by Samko et al. [23] in 1993 of the concept of variable-order operator was started. A generalization to the classical fractional calculus is given in their work by introducing the study of fractional integration and differentiation when the order is a function instead of a constant of arbitrary order [24,25]. As a result, a new generation of mathematicians and physicist is concerned with studying physical problems involving the variable order derivatives due to the property of memory incorporation for changes with time or spatial location (see, for example, [26–28]). Lorenzo and Hartley [29] gave the idea where the variable order operator is a varying function of the independent variables of differentiation or other unrelated variables which lead to the introduction of distributed order fractional operators. Most of the variable-order fractional differential equations in general do not have exact solutions, hence, the numerical methods to obtain an approximate solutions appeared as the better approach for such equations [30–36]. Recently, Abdelkawy et al. [37] suggested a novel spectral scheme to get a high precision solution for time variable fractional order mobile–immobile advection–dispersion model. Bhrawy and Zaky [34] used a numerical method for solving the variable-order nonlinear cable equation based on shifted Jacobi collocation in combination with the shifted Jacobi operational matrix for variable-order fractional derivatives. In another paper, they also suggested an accurate and effective approach to approximate the solution of functional Dirichlet boundary value problem based on shifted Chebyshev collocation procedure in combination with the shifted Chebyshev operational matrix for variable-order fractionalderivatives [38].

In recent decades, FDEs have been the focus of attention as a probable representation for description of anomalous diffusion and relaxation phenomena which are seen in a wide range of science and engineering fields [39–44], with applications in transport of fluid in porous media, diffusion of plasma, diffusion at liquid surfaces, surface growth and two-dimensional rotating flow. However, many recent researches [45,46] showed that fractional diffusion equations cannot totally represent some more complicated diffusion processes, whose diffusion behaviours depend on the spatial variation or time evolution. To handle these issues, variable-order fractional diffusion equations [35,45–47] have been recommended, for which the variable-order time fractional operator can be spatial and/or time dependent.

Taking as a guidance [34] and [35], we shall present the space–time variable-order fractional advection–dispersion equation: (1) $\begin{array}{rl}{}_0^{\mathrm{C}}{D}_t^{\alpha \left(x\right)}w(x,t)& =\mu {}_0^{\mathrm{RL}}{D}_x^{\beta (x,t)}w(x,t)\\ & -\lambda {}_0^{\mathrm{RL}}{D}_x^{\gamma (x,t)}w(x,t)+g(x,t),\\ (x,t)\in \mathrm{\Omega }& =[0,L]\times [0,T],\end{array}$(1) subject to the initial boundary conditions (2) $\begin{array}{rl}w(x,0)& =\phi \left(x\right),x\in [0,L],\\ w(0,t)& =f_1\left(t\right),w(L,t)=f_2\left(t\right),t\in [0,T],\end{array}$(2) where $\mu ,\lambda >0,0<\underset{\_}{\alpha }\le \alpha \left(x\right)\le \overline{\alpha }\le 1,1<\underset{\_}{\beta }\le \beta (x,t)\le \overline{\beta }\le 2,0<\underset{\_}{\gamma }\le \gamma (x,t)\le \overline{\gamma }\le 1.$ Here   ${}_0^{\mathrm{RL}}{D}_x^{\beta (x,t)}w(x,t)$ is the variable-order Riemann–Liouville fractional derivative defined as (3) $\begin{array}{rl}{}_0^{\mathrm{RL}}{D}_x^{\beta (x,t)}h\left(x\right)& =\left[\frac{1}{\mathrm{\Gamma }(n-\beta (x,t))}\frac{d^n}{d{\xi }^n}\right.\\ & {\left.{\int }_0^{\xi }(\xi -\eta {)}^{n-\beta (x,t)-1}h\left(\eta \right)d\eta \right]}_{\xi =x},\end{array}$(3) and ${}_0^{\mathrm{C}}{D}_t^{\alpha \left(x\right)}w(x,t)$ is the variable-order Caputo derivative defined as (4) ${}_0^{\mathrm{C}}{D}_x^{\alpha \left(x\right)}h\left(x\right)=\frac{1}{\mathrm{\Gamma }(n-\alpha \left(x\right))}{\int }_0^x\frac{h^{\left(n\right)}\left(\eta \right)d\eta }{(x-\eta {)}^{\alpha \left(x\right)-n+1}},$(4) where $n-1<\alpha \le n\in \mathbb{N}$. This study aims to create a numerical algorithm to present a better accuracy of the numerical solutions of the ST-VO-FADE (1(1) $\begin{array}{rl}{}_0^{\mathrm{C}}{D}_t^{\alpha \left(x\right)}w(x,t)& =\mu {}_0^{\mathrm{RL}}{D}_x^{\beta (x,t)}w(x,t)\\ & -\lambda {}_0^{\mathrm{RL}}{D}_x^{\gamma (x,t)}w(x,t)+g(x,t),\\ (x,t)\in \mathrm{\Omega }& =[0,L]\times [0,T],\end{array}$(1) ). The applied algorithm converts the ST-VO-FADE into a system of algebraic equations by using the so-called operational matrix of variable-order differentiation and the shifted Legendre–Gauss collocation approach. As a result, the effort performed in calculations is reduced.

The rest of this paper is organized in this way. Section 2 includes some preliminary definitions of fractional calculus and properties of the shifted Legendre polynomials. Section 3 gives a representation of the operational matrices for the variable-order fractional derivatives of the shifted Legendre polynomials. Section 4 puts the development of a collocation scheme to solve the ST-VO-FADE. In Section 5, two trial examples are applied using the proposed method. At the end, a summary of concluding remarks is given in Section 6.

2. Shifted Legendre polynomials

In this section, we introduce the main properties of the shifted Legendre polynomials which help us in what follows [16,48–51]. The classical Legendre polynomials are defined in $[-1,1]$ and may be generated from the three term recurrence relation $\begin{array}{rl}{\mathcal{L}}_0\left(z\right)& =1,{\mathcal{L}}_1\left(x\right)\\ & =x,{\mathcal{L}}_{j+1}\left(z\right)=\frac{2j+1}{j+1}z{\mathcal{L}}_j\left(z\right)\\ & -\frac{j}{j+1}{\mathcal{L}}_{j-1}\left(z\right),j\ge 2.\end{array}$

taken into the consideration the shifted Legendre polynomials ${\mathcal{L}}_j(2x/h-1)$ be represented by ${\mathcal{L}}_j^h\left(x\right).$ Then ${\mathcal{L}}_j^h\left(x\right)$ can be formed by using the following recurrence equation given by $\begin{array}{rl}{\mathcal{L}}_0^h\left(x\right)& =1,{\mathcal{L}}_1^h\left(x\right)=\frac{2x}{h}-1,\\ {\mathcal{L}}_{j+1}^h\left(x\right)& =\frac{(2j+1)(2x-h)}{(j+1)h}{\mathcal{L}}_j^h\left(x\right)\\ & -\frac{j}{j+1}{\mathcal{L}}_{j-1}^h\left(x\right),j\ge 2.\end{array}$ The orthogonality property of the shifted Legendre polynomials is (5) ${\int }_0^h{\mathcal{L}}_i^h\left(x\right){\mathcal{L}}_j^h\left(x\right)dx={\rho }_j,$(5) where ${\rho }_j={\displaystyle \frac{{\delta }_{ij}h}{2j+1}}$.

The analytic representation of ${\mathcal{L}}_j^h\left(x\right)$ of degree j is given explicitly by [16] (6) ${\mathcal{L}}_j^h\left(x\right)=\sum _{k=0}^j{\epsilon }_{j,k}^h{x}^k,$(6) where (7) ${\epsilon }_{j,k}^h=\frac{{(-1)}^{j+k}(j+k)!}{(j-k)!(k!{)}^2{h}^k},$(7) which can be rewritten in the matrix form (8) ${\mathrm{\Theta }}_{h,M}\left(x\right)=E_h{X}_M\left(x\right),$(8) where ${\epsilon }_{j,k}^h$ for $j,k=0,1,\dots ,M$ are the entries of matrix $E_h,$ (9) $\begin{array}{rl}{\mathrm{\Theta }}_{h,M}\left(x\right)& =[{\mathcal{L}}_0^h\left(x\right),{\mathcal{L}}_1^h\left(x\right)),\dots ,{\mathcal{L}}_M^h\left(x\right){]}^T,\\ X_M\left(x\right)& =[1,x,x^2,\dots ,x^M{]}^T.\end{array}$(9) Regarding the orthogonality property of the shifted Legendre polynomials (5(5) ${\int }_0^h{\mathcal{L}}_i^h\left(x\right){\mathcal{L}}_j^h\left(x\right)dx={\rho }_j,$(5) ), it is found that the matrix $E_h$ is invertible and the vector $X_M\left(x\right)$ can be expressed in terms of ${\mathrm{\Theta }}_{h,M}\left(x\right)$ in the form (10) $X_M\left(x\right)=E_h^{-1}{\mathrm{\Theta }}_{h,M}\left(x\right).$(10) The values of the shifted Legendre polynomial at the endpoint are given by (11) ${\mathcal{L}}_j^h\left(0\right)={(-1)}^j,{\mathcal{L}}_j^h\left(h\right)=1,$(11) which have significant importance which will be shown later.

Assume $w\left(x\right)$ is a square integrable function in $[0,L]$, then it can be expressed in terms of shifted Legendre polynomials as (12) $w\left(x\right)=\sum _{j=0}^{\mathrm{\infty }}{c}_j{\mathcal{L}}_j^h\left(x\right),$(12) where the coefficients $c_j$ are given by (13) $c_j=\frac{1}{{\rho }_j}{\int }_0^{h}w\left(x\right){\mathcal{L}}_j^h\mathrm{d}x,j=0,1,\dots .$(13) Considering the approximation of $w\left(x\right)$ by the first $(M+1)$-terms, then we can write (14) $w_M\left(x\right)=\sum _{j=0}^M{c}_j{\mathcal{L}}_j^h\left(x\right)={\mathbf{C}}^T{\mathrm{\Theta }}_{h,M}\left(x\right),$(14) where the shifted Legendre coefficient vector $\mathbf{C}$ is given by ${\mathbf{C}}^T=[c_0,c_1,\dots ,c_M].$

3. Differentiation matrices

In this section, we present the fractional derivative of variable order for the shifted Legendre vector in the Caputo definition sense. The expression of the first-order derivative of the shifted Legendre vector ${\mathrm{\Theta }}_{h,M}\left(x\right)$ has the form (15) $\frac{\mathrm{d}}{\mathrm{d}x}{\mathrm{\Theta }}_{h,M}\left(x\right)={\mathbf{D}}_h^{\left(1\right)}{\mathrm{\Theta }}_{h,M}\left(x\right),$(15) where ${\mathbf{D}}_h^{\left(1\right)}$ is the operational matrix of the first derivative of ${\mathrm{\Theta }}_{h,M}\left(x\right)$ with dimension of $(M+1)\times (M+1)$. (16) $\frac{\mathrm{d}}{\mathrm{d}x}{\mathrm{\Theta }}_{h,M}\left(x\right)=E_h\frac{\mathrm{d}}{\mathrm{d}x}{X}_M\left(x\right)=E_h{\mathrm{\Xi }}_M{X}_M\left(x\right),$(16) where ${\mathrm{\Xi }}_M$ is the operational matrix of the first derivative of $X_M\left(x\right)$ with dimension of $(M+1)\times (M+1)$ which is a result from (17) ${\mathrm{\Xi }}_M=\left({\xi }_{ij}\right)=\left\{\begin{array}{ll}j+1,& \mathrm{for}i=j+1,j=0,1,\dots ,M,\\ 0,& \mathrm{otherwise}.\end{array}\right.$(17) Now, by (16(16) $\frac{\mathrm{d}}{\mathrm{d}x}{\mathrm{\Theta }}_{h,M}\left(x\right)=E_h\frac{\mathrm{d}}{\mathrm{d}x}{X}_M\left(x\right)=E_h{\mathrm{\Xi }}_M{X}_M\left(x\right),$(16) ) and (10(10) $X_M\left(x\right)=E_h^{-1}{\mathrm{\Theta }}_{h,M}\left(x\right).$(10) ), then it is easy to write (18) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}x}{\mathrm{\Theta }}_{h,M}\left(x\right)& =E_h{\mathrm{\Xi }}_M{E}_h^{-1}{\mathrm{\Theta }}_{h,M}\left(x\right)\\ & ={\mathbf{D}}_h^{\left(1\right)}{\mathrm{\Theta }}_{h,M}\left(x\right).\end{array}$(18) Accordingly, it can be deduce that (19) ${\mathbf{D}}_h^{\left(1\right)}=E_h{\mathrm{\Xi }}_M{E}_h^{-1}.$(19) Using (18(18) $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}x}{\mathrm{\Theta }}_{h,M}\left(x\right)& =E_h{\mathrm{\Xi }}_M{E}_h^{-1}{\mathrm{\Theta }}_{h,M}\left(x\right)\\ & ={\mathbf{D}}_h^{\left(1\right)}{\mathrm{\Theta }}_{h,M}\left(x\right).\end{array}$(18) ) repeatedly, gives the relation (20) $\begin{array}{rl}\frac{{\mathrm{d}}^p}{\mathrm{d}{x}^p}{\mathrm{\Theta }}_{h,M}\left(x\right)& =\left({\mathbf{D}}_h^{\left(1\right)}{)}^p{\mathrm{\Theta }}_{h,M}\right(x)={\mathbf{D}}_h^{\left(p\right)}{\mathrm{\Theta }}_{h,M}(x),\\ p& =1,2,\dots ,\end{array}$(20) where $p\in \mathbb{N}$.

In Theorem 3.1 [34], the generalization of the operational matrix of the derivative of the shifted Legendre polynomials can be extended for variable-order fractional derivatives as (21) ${}_0^{\mathrm{C}}{D}_x^{\alpha (x,t)}{\mathrm{\Theta }}_{h,M}\left(x\right)={\mathbf{D}}_{h,\alpha (x,t)}{\mathrm{\Theta }}_{h,M}\left(x\right),$(21) where $n-1<{\alpha }_{min}<\alpha (x,t)<{\alpha }_{max}<n$ and ${\mathbf{D}}_{h,\alpha (x,t)}$ is an $(M+1)\times (M+1)$ matrix of the following form: ${\mathbf{D}}_{h,\alpha (x,t)}=x^{-\alpha (x,t)}{E}_h\mathbf{B}{E}_h^{-1},$ where $E_h$ is defined in (8(8) ${\mathrm{\Theta }}_{h,M}\left(x\right)=E_h{X}_M\left(x\right),$(8) ) and $\mathbf{B}$ is a $(M+1)\times (M+1)$ matrix and its elements, $b_{ij};0\le i,j\le M$ are given as follows: $b_{ij}=\left\{\begin{array}{ll}\frac{\mathrm{\Gamma }(i+1)}{\mathrm{\Gamma }(i+1-\alpha (x,t\left)\right)},& \mathrm{for}i=j,j=n,n+1,\dots ,M,\\ 0,& \mathrm{otherwise}.\end{array}\right.$

4. Legendre spectral collocation method

In this section, the shifted Legendre collocation method is applied to solve ST-VO-FADE (1(1) $\begin{array}{rl}{}_0^{\mathrm{C}}{D}_t^{\alpha \left(x\right)}w(x,t)& =\mu {}_0^{\mathrm{RL}}{D}_x^{\beta (x,t)}w(x,t)\\ & -\lambda {}_0^{\mathrm{RL}}{D}_x^{\gamma (x,t)}w(x,t)+g(x,t),\\ (x,t)\in \mathrm{\Omega }& =[0,L]\times [0,T],\end{array}$(1) )–(2(2) $\begin{array}{rl}w(x,0)& =\phi \left(x\right),x\in [0,L],\\ w(0,t)& =f_1\left(t\right),w(L,t)=f_2\left(t\right),t\in [0,T],\end{array}$(2) ).An approximate solution of $w(x,t)$ can be expressed by the series form of the double shifted Legendre polynomials as (22) $\begin{array}{rl}w(x,t)\simeq w_{N,M}(x,t)& =\sum _{i=0}^N\sum _{j=0}^M{a}_{i,j}{\mathcal{L}}_i^{\tau }\left(t\right){\mathcal{L}}_j^L\left(x\right)\\ & ={\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right),\end{array}$(22) where $\mathbf{A}$ is a matrix of unknown elements of order $(N+1)\times (M+1)$.

Now, using Equations (20(20) $\begin{array}{rl}\frac{{\mathrm{d}}^p}{\mathrm{d}{x}^p}{\mathrm{\Theta }}_{h,M}\left(x\right)& =\left({\mathbf{D}}_h^{\left(1\right)}{)}^p{\mathrm{\Theta }}_{h,M}\right(x)={\mathbf{D}}_h^{\left(p\right)}{\mathrm{\Theta }}_{h,M}(x),\\ p& =1,2,\dots ,\end{array}$(20) ), (21(21) ${}_0^{\mathrm{C}}{D}_x^{\alpha (x,t)}{\mathrm{\Theta }}_{h,M}\left(x\right)={\mathbf{D}}_{h,\alpha (x,t)}{\mathrm{\Theta }}_{h,M}\left(x\right),$(21) ) and (22(22) $\begin{array}{rl}w(x,t)\simeq w_{N,M}(x,t)& =\sum _{i=0}^N\sum _{j=0}^M{a}_{i,j}{\mathcal{L}}_i^{\tau }\left(t\right){\mathcal{L}}_j^L\left(x\right)\\ & ={\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right),\end{array}$(22) ), it can be written as (23) $\begin{array}{rl}& {}_0^{\mathrm{RL}}{D}_x^{\gamma (x,t)}w(x,t)\\ & ={}_0^{\mathrm{C}}{D}_x^{\gamma (x,t)}w(x,t)+\frac{f_1\left(t\right)x^{-\gamma (x,t)}}{\mathrm{\Gamma }(1-\gamma (x,t\left)\right)}\\ & \simeq {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_{L,\gamma (x,t)}{\mathrm{\Theta }}_{L,M}\left(x\right)\\ & +\frac{f_1\left(t\right)x^{-\gamma (x,t)}}{\mathrm{\Gamma }(1-\gamma (x,t\left)\right)}.\end{array}$(23) (24) $\begin{array}{rl}& {}_0^{\mathrm{RL}}{D}_x^{\beta (x,t)}w(x,t)\\ & ={}_0^{\mathrm{C}}{D}_x^{\beta (x,t)}w(x,t)+\frac{f_2\left(t\right)x^{-\beta (x,t)}}{\mathrm{\Gamma }(1-\beta (x,t\left)\right)}\\ & +\frac{{\mathrm{\partial }}_{x}w(0,t)x^{1-\beta (x,t)}}{\mathrm{\Gamma }(2-\beta (x,t\left)\right)}\\ & \simeq {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_{L,\beta (x,t)}{\mathrm{\Theta }}_{L,M}\left(x\right)\\ & +\frac{f_2\left(t\right)x^{-\beta (x,t)}}{\mathrm{\Gamma }(1-\beta (x,t\left)\right)}\\ & +{\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_L^{\left(1\right)}{\mathrm{\Theta }}_{L,M}\left(0\right)\frac{x^{1-\beta (x,t)}}{\mathrm{\Gamma }(2-\beta (x,t\left)\right)},\end{array}$(24) (25) $\begin{array}{rl}{}_0^{\mathrm{C}}{D}_x^{\alpha \left(x\right)}w(x,t)& ={\mathrm{\Theta }}_{\tau ,N}^T\left(t\right){\mathbf{D}}_{\tau ,\alpha \left(x\right)}^T\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right),\end{array}$(25) (26) $\begin{array}{rl}& \begin{array}{rl}{w}_{N,M}(x,0)& ={\mathrm{\Theta }}_{\tau ,N}^T\left(0\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right),\\ w_{N,M}(0,t)& ={\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(0\right),\\ w_{N,M}(L,t)& ={\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(L\right).\end{array}\end{array}$(26) Employing Equations (22(22) $\begin{array}{rl}w(x,t)\simeq w_{N,M}(x,t)& =\sum _{i=0}^N\sum _{j=0}^M{a}_{i,j}{\mathcal{L}}_i^{\tau }\left(t\right){\mathcal{L}}_j^L\left(x\right)\\ & ={\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right),\end{array}$(22) )–(26(26) $\begin{array}{rl}& \begin{array}{rl}{w}_{N,M}(x,0)& ={\mathrm{\Theta }}_{\tau ,N}^T\left(0\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right),\\ w_{N,M}(0,t)& ={\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(0\right),\\ w_{N,M}(L,t)& ={\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(L\right).\end{array}\end{array}$(26) ) in Equations (1(1) $\begin{array}{rl}{}_0^{\mathrm{C}}{D}_t^{\alpha \left(x\right)}w(x,t)& =\mu {}_0^{\mathrm{RL}}{D}_x^{\beta (x,t)}w(x,t)\\ & -\lambda {}_0^{\mathrm{RL}}{D}_x^{\gamma (x,t)}w(x,t)+g(x,t),\\ (x,t)\in \mathrm{\Omega }& =[0,L]\times [0,T],\end{array}$(1) )–(2(2) $\begin{array}{rl}w(x,0)& =\phi \left(x\right),x\in [0,L],\\ w(0,t)& =f_1\left(t\right),w(L,t)=f_2\left(t\right),t\in [0,T],\end{array}$(2) ) yields (27) $\begin{array}{rl}& {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right){\mathbf{D}}_{\tau ,\alpha \left(x\right)}^T\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right)\\ & =\mu {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_{L,\gamma (x,t)}{\mathrm{\Theta }}_{L,M}\left(x\right)+\frac{f_1\left(t\right)x^{-\gamma (x,t)}}{\mathrm{\Gamma }(1-\gamma (x,t\left)\right)}\\ & -\lambda {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_{L,\beta (x,t)}{\mathrm{\Theta }}_{L,M}\left(x\right)+\frac{f_2\left(t\right)x^{-\beta (x,t)}}{\mathrm{\Gamma }(1-\beta (x,t\left)\right)}\\ & +{\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_L^{\left(1\right)}{\mathrm{\Theta }}_{L,M}\left(0\right)\frac{x^{1-\beta (x,t)}}{\mathrm{\Gamma }(2-\beta (x,t\left)\right)}+g(x,t),\end{array}$(27) (28) $\begin{array}{rl}& \begin{array}{rl}& {\mathrm{\Theta }}_{\tau ,N}^T\left(0\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right)=\phi \left(x\right),\\ & {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(0\right)=g_1\left(t\right),\\ & {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(L\right)=g_2\left(t\right).\end{array}\end{array}$(28) Now, we apply directly the collocation method to solve (27(27) $\begin{array}{rl}& {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right){\mathbf{D}}_{\tau ,\alpha \left(x\right)}^T\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right)\\ & =\mu {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_{L,\gamma (x,t)}{\mathrm{\Theta }}_{L,M}\left(x\right)+\frac{f_1\left(t\right)x^{-\gamma (x,t)}}{\mathrm{\Gamma }(1-\gamma (x,t\left)\right)}\\ & -\lambda {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_{L,\beta (x,t)}{\mathrm{\Theta }}_{L,M}\left(x\right)+\frac{f_2\left(t\right)x^{-\beta (x,t)}}{\mathrm{\Gamma }(1-\beta (x,t\left)\right)}\\ & +{\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_L^{\left(1\right)}{\mathrm{\Theta }}_{L,M}\left(0\right)\frac{x^{1-\beta (x,t)}}{\mathrm{\Gamma }(2-\beta (x,t\left)\right)}+g(x,t),\end{array}$(27) )–(28(28) $\begin{array}{rl}& \begin{array}{rl}& {\mathrm{\Theta }}_{\tau ,N}^T\left(0\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right)=\phi \left(x\right),\\ & {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(0\right)=g_1\left(t\right),\\ & {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(L\right)=g_2\left(t\right).\end{array}\end{array}$(28) ). Using the nodes $x_i(0\le i\le M)$ which are the shifted Legendre–Gauss–Lobatto roots of ${\mathcal{L}}_M^L\left(x\right)$ and $t_j(0\le j\le N-1)$ is the shifted Legendre roots of ${\mathcal{L}}_N^{\tau }\left(x\right)$. We substitute these nodes in (27(27) $\begin{array}{rl}& {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right){\mathbf{D}}_{\tau ,\alpha \left(x\right)}^T\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right)\\ & =\mu {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_{L,\gamma (x,t)}{\mathrm{\Theta }}_{L,M}\left(x\right)+\frac{f_1\left(t\right)x^{-\gamma (x,t)}}{\mathrm{\Gamma }(1-\gamma (x,t\left)\right)}\\ & -\lambda {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_{L,\beta (x,t)}{\mathrm{\Theta }}_{L,M}\left(x\right)+\frac{f_2\left(t\right)x^{-\beta (x,t)}}{\mathrm{\Gamma }(1-\beta (x,t\left)\right)}\\ & +{\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathbf{D}}_L^{\left(1\right)}{\mathrm{\Theta }}_{L,M}\left(0\right)\frac{x^{1-\beta (x,t)}}{\mathrm{\Gamma }(2-\beta (x,t\left)\right)}+g(x,t),\end{array}$(27) )–(28(28) $\begin{array}{rl}& \begin{array}{rl}& {\mathrm{\Theta }}_{\tau ,N}^T\left(0\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right)=\phi \left(x\right),\\ & {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(0\right)=g_1\left(t\right),\\ & {\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(L\right)=g_2\left(t\right).\end{array}\end{array}$(28) ); therefore, the collocation scheme can be written as (29) $\begin{array}{rl}& {\mathrm{\Theta }}_{\tau ,N}^T\left(t_i\right){\mathbf{D}}_{\tau ,\alpha \left(x_j\right)}^T\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x_j\right)\\ & =\mu {\mathrm{\Theta }}_{\tau ,N}^T\left(t_i\right)\mathbf{A}{\mathbf{D}}_{L,\gamma (x_j,t_i)}{\mathrm{\Theta }}_{L,M}\left(x_j\right)+\frac{f_1\left(t_i\right)x_j^{-\gamma (x_j,t_i)}}{\mathrm{\Gamma }(1-\gamma (x_j,t_i\left)\right)}\\ & -\lambda {\mathrm{\Theta }}_{\tau ,N}^T\left(t_i\right)\mathbf{A}{\mathbf{D}}_{L,\beta (x_j,t_i)}{\mathrm{\Theta }}_{L,M}\left(x_j\right)+\frac{f_2\left(t_i\right)x_j^{-\beta (x_j,t_i)}}{\mathrm{\Gamma }(1-\beta (x_j,t_i\left)\right)}\\ & +{\mathrm{\Theta }}_{\tau ,N}^T\left(t_i\right)\mathbf{A}{\mathbf{D}}_L^{\left(1\right)}{\mathrm{\Theta }}_{L,M}\left(0\right)\frac{x_j^{1-\beta (x_j,t_i)}}{\mathrm{\Gamma }(2-\beta (x_j,t_i\left)\right)}+g(x_j,t_i),\\ & (1\le j\le M-1),(0\le i\le N-1),\end{array}$(29) (30) $\begin{array}{rl}& \begin{array}{rl}& {\mathrm{\Theta }}_{\tau ,N}^T\left(0\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x_j\right)=\phi \left(x_j\right),0\le j\le M,\\ & {\mathrm{\Theta }}_{\tau ,N}^T\left(t_i\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(0\right)=f_1\left(t_i\right),0\le i\le N-1,\\ & {\mathrm{\Theta }}_{\tau ,N}^T\left(t_i\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(L\right)=f_2\left(t_i\right),0\le i\le N-1.\end{array}\end{array}$(30) This generates a system of $(N+1)\times (M+1)$ nonlinear algebraic equations in the required double shifted Legendre coefficients $a_{i,j},(i=0,1,\dots ,M;j=0,1,\dots ,N)$, which is solved by using any standard iteration technique, like Newton's iteration method. As a result, the approximate solution (22(22) $\begin{array}{rl}w(x,t)\simeq w_{N,M}(x,t)& =\sum _{i=0}^N\sum _{j=0}^M{a}_{i,j}{\mathcal{L}}_i^{\tau }\left(t\right){\mathcal{L}}_j^L\left(x\right)\\ & ={\mathrm{\Theta }}_{\tau ,N}^T\left(t\right)\mathbf{A}{\mathrm{\Theta }}_{L,M}\left(x\right),\end{array}$(22) ) can be obtained.

5. Numerical results

To demonstrate the effectiveness of the proposed method, two test examples are carried out in this section.

Example 1

Consider the following initial boundary value problem of ST-VO-FADE [35] (31) $\begin{array}{rl}& {}_0^{\mathrm{C}}{D}_t^{\alpha \left(x\right)}w(x,t)\\ & ={}_0^{\mathrm{C}}{D}_x^{\beta (x,t)}w(x,t)-{}_0^{\mathrm{C}}{D}_x^{\gamma (x,t)}u(x,t)+g(x,t),\\ (x,t)& \in (0,1)\times (0,1],\end{array}$(31) subject to the initial boundary conditions $\begin{array}{rl}w(x,0)& =10x^2(1-x),x\in (0,1),\\ w(0,t)& =u(1,t)=0,t\in (0,T],\end{array}$

where $\begin{array}{rl}g(x,t)& =\frac{10x^2(1-x)t^{1-\alpha \left(x\right)}}{\mathrm{\Gamma }(2-\alpha (x\left)\right)}\\ & -10(t+1)[\frac{2x^{2-\beta (x,t)}}{\mathrm{\Gamma }(3-\beta (x,t\left)\right)}\\ & -\frac{6x^{3-\beta (x,t)}}{\mathrm{\Gamma }(4-\beta (x,t\left)\right)}]\\ & +10(t+1)[\frac{2x^{2-\gamma (x,t)}}{\mathrm{\Gamma }(3-\gamma (x,t\left)\right)}-\frac{6x^{3-\gamma (x,t)}}{\mathrm{\Gamma }(4-\gamma (x,t\left)\right)}].\end{array}$

The exact solution is $w(x,t)=10(t+1)x^2(1-x).$ We solve the equation with $\begin{array}{rl}\alpha \left(x\right)& =1-0.5e^{-x},\\ \beta (x,t)& =1.7+0.5e^{-x^2/1000-t/50-1},\\ \gamma (x,t)& =0.7+0.5e^{-x^2/1000-t/50-1}.\end{array}$ The space–time graphs of the approximate solution and the absolute errors (AE) at N=M=4 are shown in Figure 1(left) and (right), respectively.

Figure 1. The space–time graph of approximate solution (left) and AE (right) at N = M = 4 for Example 1.

A comparison of the presented method at N=M=5 with the numerical method proposed in [35] is listed in Table 1.

Table 1. Comparing maximum absolute errors of the proposed method and Implicit Euler approximation (IEA)[35] at t=1 for Example 1.

$x_i$	IEA [35]	Our results $(N=M=4)$
0.1	$2.391\times {10}^{-3}$	$2.220\times {10}^{-15}$
0.2	$5.601\times {10}^{-3}$	$4.440\times {10}^{-15}$
0.3	$9.237\times {10}^{-3}$	$1.231\times {10}^{-15}$
0.4	$1.282\times {10}^{-2}$	$2.775\times {10}^{-15}$
0.5	$1.582\times {10}^{-2}$	$2.220\times {10}^{-15}$
0.6	$1.766\times {10}^{-2}$	$8.881\times {10}^{-15}$
0.7	$1.772\times {10}^{-2}$	$1.665\times {10}^{-15}$
0.8	$1.532\times {10}^{-2}$	$7.216\times {10}^{-15}$
0.9	$9.687\times {10}^{-3}$	$2.220\times {10}^{-15}$

Example 2

Consider the following ST-VO-FADE [35] (32) $\begin{array}{rl}{}_0^{\mathrm{C}}{D}_t^{\alpha \left(x\right)}w(x,t)& =2{}_0^{\mathrm{C}}{D}_x^{\beta (x,t)}w(x,t)\\ & -{}_0^{\mathrm{C}}{D}_x^{\gamma (x,t)}w(x,t)+g(x,t),\\ (x,t)& \in (0,1)\times (0,1],\end{array}$(32) with the initial and boundary conditions: $\begin{array}{rl}w(x,0)& =5x(1-x),x\in (0,1),\\ w(0,t)& =w(1,t)=0,t\in (0,T],\end{array}$ where $\begin{array}{rl}g(x,t)& =\frac{10x^2(1-x)t^{2-\alpha \left(x\right)}}{\mathrm{\Gamma }(3-\alpha (x\left)\right)}\\ & -10(t^2+1)\left[\frac{2x^{2-\beta (x,t)}}{\mathrm{\Gamma }(3-\beta (x,t\left)\right)}-\frac{6x^{3-\beta (x,t)}}{\mathrm{\Gamma }(4-\beta (x,t\left)\right)}\right]\\ & +5(t^2+1)\left[\frac{2x^{2-\gamma (x,t)}}{\mathrm{\Gamma }(3-\gamma (x,t\left)\right)}-\frac{6x^{3-\gamma (x,t)}}{\mathrm{\Gamma }(4-\gamma (x,t\left)\right)}\right].\end{array}$ The exact solution is $w(x,t)=5(t^2+1)x^2(1-x).$

We solve the equation with $\begin{array}{rl}\alpha \left(x\right)& =0.8-0.01\ln \left(5x\right),\\ \beta (x,t)& =1.8+0.01x^2{t}^2,\\ \gamma (x,t)& =0.8+0.01x^2\sin t.\end{array}$

Figure 2 shows the space–time graph of approximate solution (left) and AE (right) of (32(32) $\begin{array}{rl}{}_0^{\mathrm{C}}{D}_t^{\alpha \left(x\right)}w(x,t)& =2{}_0^{\mathrm{C}}{D}_x^{\beta (x,t)}w(x,t)\\ & -{}_0^{\mathrm{C}}{D}_x^{\gamma (x,t)}w(x,t)+g(x,t),\\ (x,t)& \in (0,1)\times (0,1],\end{array}$(32) ) at $N=M=4.$

Figure 2. The space–time graph of approximate solution (left) and AE (right) at N=M=4 for Example 2.

Figure 3 shows the solution behaviour of (32(32) $\begin{array}{rl}{}_0^{\mathrm{C}}{D}_t^{\alpha \left(x\right)}w(x,t)& =2{}_0^{\mathrm{C}}{D}_x^{\beta (x,t)}w(x,t)\\ & -{}_0^{\mathrm{C}}{D}_x^{\gamma (x,t)}w(x,t)+g(x,t),\\ (x,t)& \in (0,1)\times (0,1],\end{array}$(32) ) at $T=0.25,T=0.5$ and T=1 with N=M=6.

Figure 3. The solution behaviour of (32) at T=0:25; T=0:5 and T=1 with N=M=6.

6. Conclusion

In this paper, we have proposed fast and precise algorithm based on shifted Legende collocation technique combined with the associated operational matrices of variable-order fractional derivatives. This algo-rithm was employed for solving the space–time variable-order fractional advection–dispersion model with Caputo time variable fractional derivative and Riemann–Liouville space variable fractional derivatives. This algorithm has the advantage of transforming the problem into the solution of a system of algebraic equations which greatly simplifying the problem. Finally, two numerical examples have been presented to demonstrate the efficiency of the proposed algorithm.

Authors' contributions

The authors have equal contributions to each part of this paper. All the authors read and approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

F. Mallawi  http://orcid.org/0000-0002-3434-443X

J. F. Alzaidy  http://orcid.org/0000-0003-2318-8176

R. M. Hafez  http://orcid.org/0000-0001-9533-3171

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.