A Hahn computational operational method for variable order fractional mobile–immobile advection–dispersion equation

In this paper, we consider the discrete Hahn polynomials f H n g and investigate their application for numerical solutions of the time fractional variable order mobile–immobile advection–dispersion model which is advantageous for modeling dynamical systems. This paper presented the operational matrix of derivative of discrete Hahn polynomials. The main advantage of approximating a continuous function by Hahn polynomials is that they have a spectral accuracy in the interval [0, N ]. Furthermore, for computing the coefﬁcients of the expansion u ð x Þ ¼ P 1 n ¼ 0 c n H n ð x Þ , we have to only compute a summation and the calculation of coefﬁcients is exact. Also an upper bound for the error of the presented method, with equidistant nodes, is investigated. Illustrative examples are provided to show the accuracy and efﬁciency of the presented method. Using a small number of Hahn polynomials, signiﬁcant results are achieved which are compared to other methods.


Introduction
Fractional calculus is the generalization of the ordinary calculus which studies integrals and derivatives of noninteger (real or complex) order. Some numerical methods have been used to solve functional equations containing fractional derivatives, such as [18,22,27,31]. In a wide class of physical, dynamical and biological models and complex problems, the constant order fractional equation cannot describe the characterization of problems. In 1993, intellectual curiosity of Samko and Ross led to an extension of the classical fractional calculus [32]. They presented continuously varying order q(x) for differential and integral operators and defined these operators in two ways. The first way was a direct definition and the second way was based on Fourier transforms. Variable-order (VO) fractional calculus is good at representing the memory property which changes with time or location. Applications of the VO idea have been developed by several researchers such as Lorenzo and Hartley [25], Ingman and Suzdalnitsky [17], Coimbra [9] and others. VO fractional calculus as a powerful tool has recently been applied in the fluid dynamics and control of nonlinear viscoelasticity [29], mechanics [9], medical imaging [38], processing of geographical data [10], etc. The analytic results on the existence and uniqueness of the solutions for a generalized fractional differential equations with VO operators have been discussed in [30].
Solving VO fractional problems is quite difficult. Therefore, efficient numerical techniques are necessary to be developed. A finite difference technique is used for solving VO fractional integro-differential equations [37]. Recently, spectral methods using continuous orthogonal polynomials, such as Jacobi, Chebyshev and Legendre polynomials, have been developed for solving different kinds of VO fractional differential equations. Chen et al. proposed Legendre wavelets functions to solve a class of nonlinear VO fractional differential equations [8]. Bernstein polynomials are used to numerically solve the VO fractional partial differential equations by Wang et al. [35]. In [7], numerical solutions of the VO linear cable equation with adopted Bernstein polynomials basis on the interval [0, R] are presented. Although in most cases continuous orthogonal polynomials are used as basis functions for the approximate solution of equations, recently discrete orthogonal polynomials have been noticed for solving stochastic differential equations [36] and in numerical fluid dynamics problems [13] because of the behavior and property of these polynomials [15].
Discrete orthogonal polynomials are orthogonal with respect to a weighted discrete inner product. Classical orthogonal polynomial (continuous/discrete) are a class of the orthogonal polynomials associated with the Wiener-Askey polynomials including Hermite, Laguerre, Jacobi as continuous, and Charlier, Meixner, Krawtchouk and Hahn as discrete polynomials [24]. That in this paper we focus on Hahn polynomials.
Advection-dispersion equation and its extensions (such as the mobile-immobile equation) are a combination of the diffusion and advection equations. These equations are used to the modeling of transformation of pollutants, energy, subsurface water flows, deeper river flows, streams, and groundwater [4,11,12,19,23]. Recently, VO fractional diffusion equation is expanded to describe time-dependent anomalous diffusion and diffusion process in inhomogeneous porous media and for the description of complex dynamical problems [40].
In this work, we investigate the VO fractional problem arises in a mobile-immobile advection-dispersion equation (VOFMIE) [39]. This model is obtained from the standard advection-dispersion equation by adding the variable timefractional derivative in the Caputo sense of order 0\qðx; tÞ 1. In [26], this model is used to simulate solute transport in watershed catchments and rivers. Some numerical methods have been developed for the solution of VOFMIE. Authors of [20] are used reproducing kernel theory and collocation method for solving the VOFMIE. An implicit Euler approximation for the VOFMIE has been explained in [39] Also, Chebyshev wavelets method [16] is employed to solve this VO equation. In [34], the VOFMIE in 2-D arbitrary domains is introduced and a meshless method based on the MLS approach is proposed to solve it. In [1], a method, based on shifted Jacobi Gauss-Lobatto and shifted Jacobi Gauss-Radau spectral methods, is presented for solving the mobile-immobile advection-dispersion model with a Coimbra time variable fractional derivative. with the following initial and boundary conditions: uð0; tÞ ¼ w 1 ðtÞ; uðX; tÞ ¼ w 2 ðtÞ; ð3Þ

Definition and properties of Hahn polynomials
First, we state some definitions.

Definition 1
The Pochhammer symbol is defined by the following relations: From the explicit formula (6) and the definition of Pochhammer symbol, it can be seen that H n ðx; a; b; NÞ is a polynomial of variable x with degree exactly n. The Hahn polynomials H n ðx; a; b; NÞ are classical discrete orthogonal polynomials of degree n on the interval I :¼ ½0; N. They are orthogonal on I with respect to the discrete inner product: f ðxÞgðxÞwðxÞ: where wðxÞ : ½a; b ! R is a non-negative weight-function given by They are normalized by hH m ðx; a; b; NÞ; H n ðx; a; b; NÞi The orthogonality property of these basis polynomials is given in [21]. Also, the Hahn polynomials and Jacobi polynomials P k ¼ P k ða; bÞ are related by the following:

Expansion of Hahn polynomials in terms of Taylor basis
To obtain the expansion of Hahn polynomials as a Taylor series, we use the relation: where S ðmÞ k are the Stirling numbers of the first kind (see page 824 of [2]). Closed forms for Stirling numbers of the first are as the following [2]: where s ðiÞ k are the Stirling numbers of the second kind, namely Using Eq. (8), we have Remark 1 Matrix A is an upper triangular matrix and also from definition (), for i ¼ 0; 1; . . .; N, and a; b ! À 1, we have It can be easily seen from (9) that S ð0Þ 0 ¼ 1. Therefore, detðAÞ 6 ¼ 0 and A is an invertible matrix.

Function approximations by Hahn polynomials
In Theorem 1.1 from [14], Goertz and Ö ffner describe the expansion of a function by Hahn polynomials and they presented the following result: Let u 2 C½a; b be a function of bounded variation. Then, the series expansion P N i¼0 c i H i ðxÞ of a function u, by Hahn polynomials, converges pointwise, if the series expansion P N i¼0 d i P i ðxÞ of the function u, by Jacobi polynomials, converges pointwise and n 4 =N ! 0 for n; N ! 1. From Theorem 3.1 in [15], spectral accuracy follows directly for Hahn polynomials. Therefore, a continuous function u(x), of bounded variation, can be expanded by truncated Hahn series as follows: where In the approximation of a function u(x) on [a, b], in terms of continuous orthogonal polynomials / i ðxÞ, we use the expansion uðxÞ ¼ P N i¼0 c i / i ðxÞ, where the coefficients u i are computed by The integral of (15) is computed by numerical quadrature where usually is not exact. But using discrete orthogonal polynomials, we have only to compute summation (14) and the calculation of coefficients is exact. This is our motivation for using discrete Hahn polynomials in the presented method.
Similarly, a function of two independent variables uðx; tÞ 2 X :¼ ½0; X Â ½0; T may be expanded in terms of the double Hahn polynomials as: where coefficient matrix U is given by where u ij ¼ hHðxÞ; hU; HðtÞii:

Operational matrix of derivatives for Hahn polynomials
Here, we present the ordinary and VO differentiation matrix of the Hahn polynomials in the Caputo sense. We need the following formula of the VO fractional derivative of function t n , in Caputo's sense of order 0\qðx; tÞ 1 [28]: : The last vector can be written as: Substituting Eqs. (25), (26), (27), and (24) into Eq. (1) leads to the following relation: To discretize the Eq. (28) Equation AUB þ CUD ¼ E is a linear system in terms of unknown matrix U. To solve this system, we use the theorem 2.12 and corollary 2.4 in [33].

Error bound
In this section, an upper bound for the error of the presented approximation scheme, with equidistant nodes, is obtained with a similar procedure as in [5].
Since u(x, t) is a smooth function on X, then there exist the constants K 1 , K 2 , and K 3 such that On the other hand, for equidistant points, we have and Now, from Eqs. (32)-(37), Eq. (31) becomes kuðx; tÞ À u I;J ðx; tÞk 1 In this manner, an upper bound of the absolute errors is obtained for the approximate solutions with equidistant point.

Numerical examples
In this section, numerical examples are provided to illustrate the efficiency of the mentioned approach. We compare our results to the exact solutions at T, via the following norm: where ke N k 1 denotes the error corresponding to polynomial of degree N. Also, to show the efficiency, we report the convergence order (CO) of our method for two last examples. CO is defined by [6] CO ¼ with the following initial and boundary conditions: where X ¼ ½0; 1 Â ½0; T, qðx; tÞ ¼ 1 À 1 2 e Àxt , and f ¼ f 1 þ f 2 , with: f 1 ðx; tÞ ¼ 10x 2 ð1 À xÞ 2 þ 10x 2 ð1 À xÞ 2 t 1Àqðx;tÞ Cð2 À qðx; tÞÞ ; The exact solution is u E ðx; tÞ ¼ 10ðt þ 1Þx 2 ð1 À xÞ 2 .   The absolute errors of the numerical solutions, at T ¼ 1; 2; 4, and for a ¼ 0:01, b ¼ 0:01, and N ¼ 5 are shown in Table 1. A comparison of the numerical solutions is made by the results reported in [39] and [20] at T ¼ 1, by Finite Difference Method (FD) and Reproducing Kernel Method (RKM), respectively. Also, the results of the presented method at T ¼ 2 and T ¼ 4 show that the method is accurate even for larger domains. Fig. 1 shows the results of the approximate and exact solutions at T ¼ 1; 2; 4. The absolute error for this approach is shown in Fig. 2 where X ¼ ½0; 1 Â ½0; T, with the following initial and boundary conditions:  f 1 ðx; tÞ ¼ 5xð1 À xÞ þ 5xð1 À xÞt 1Àqðx;tÞ Cð2 À qðx; tÞÞ ; The exact solution is u E ðx; tÞ ¼ 5ðt þ 1Þxð1 À xÞ.
The errors of the numerical solutions for a ¼ 1, b ¼ 1, and N ¼ 10, at T ¼ 1; 2; 4, are shown in Table 2. Also, a comparison of the numerical solutions, at T ¼ 1, is made by the results of RKM [20], in Table 2. In Fig. 3, the results of the approximate and the exact solutions at T ¼ 1; 2; 4 are compared. The absolute errors for this approach are shown in Fig. 4.
where X ¼ ½0; 1 Â ½0; 1, with the following initial and boundary conditions: The exact solution is u E ðx; tÞ ¼ t 3 e x . The errors of the numerical solutions for a ¼ b ¼ 1 and different N are shown in Table 3. We investigate the convergence order of our method, and a comparison of the numerical solutions is made by the results of Chebyshev wavelets method (CWs) [16], in Table 4. The absolute errors for this approach, at T ¼ 1, are shown in Fig. 6. The exact solution is u E ðx; tÞ ¼ t 2 sinðpx=LÞ. The errors of the numerical solutions for a ¼ b ¼ 1 and different N are shown in Table 5. The convergence order of our method is investigated in Table 6. Also, the absolute errors for this approach, at T ¼ 0:5, are shown in Fig. 7.

Discussion and conclusion
In this paper, we present an operational matrix method based on Hahn polynomials for solving the time VO fractional mobile-immobile advection-dispersion model. This method converts the VO fractional equation into an algebraic system which is solved by a technique of linear algebra. Using the discrete Hahn polynomials in a projection approach will not result in any numerical calculation errors.
As the main result, a new operational matrix of VO fractional derivative, in the Caputo sense, is derived for the Hahn polynomials. Then, it is used to obtain an approximate solution to the problem under study. Also, an upper bound for the error of the presented method, with equidistant nodes, is investigated. Numerical examples show the accuracy and the efficiency of the presented method. An advantage of using the Hahn polynomials is that the numerical results achieved only using a small number of bases are accurate in a larger intervals and significant results are achieved. A comparison of the numerical solutions by FDM [39] and RKM [20] shows that this technique is accurate enough to be known as a powerful device. This method can be applied to solve other types of VO fractional functional equations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.