Vieta–Lucas polynomials for solving a fractional-order mathematical physics model

In this article, a fractional-order mathematical physics model, advection–dispersion equation (FADE), will be solved numerically through a new approximative technique. Shifted Vieta–Lucas orthogonal polynomials will be considered as the main base for the desired numerical solution. These polynomials are used for transforming the FADE into an ordinary differential equations system (ODES). The nonstandard finite difference method coincidence with the spectral collocation method will be used for converting the ODES into an equivalence system of algebraic equations that can be solved numerically. The Caputo fractional derivative will be used. Moreover, the error analysis and the upper bound of the derived formula error will be investigated. Lastly, the accuracy and efficiency of the proposed method will be demonstrated through some numerical applications.

lutions depend on several techniques such as finite difference, finite volume, variational iteration, Legendre polynomials, Chebyshev collocation, homotopy perturbation, operational matrix, variational iteration, Adams-Bashforth, nonstandard finite difference, sinccollocation, compact finite difference, tau method, block pulse, decomposition, radial basis, Taylor collocation, and wavelets spectral (see, for example, [5,11,24,25,35,38]). All these methods introduce numerical solutions for many types of fractional-order differential equations. In this work, we strive to solve the fractional-order advection-dispersion equation numerically via a nonstandard finite difference method besides a collocation method that depends on a new class of orthogonal polynomials (Vieta-Lucas).
Here u(x, t) refers to the solute concentration, the fractional derivative term ∂ α u(x,t) ∂x α is a dispersion function with the dispersion coefficient λ, the second fractional derivative ∂ β u(x,t) ∂x β is the advection term with the average fluid velocity coefficient μ, and the last term s(x, t) is the source/sink term. Besides that, the fractional-order terms are considered as the fractional derivatives in the Caputo operator of differentiation.
The classical order advection-dispersion equation can be obtained by using the values α = 2, β = 1 in Eq. (1) to have: Equation (1) is a generalization of Eq. (4), therefore, the researchers have a big chance for modeling many problems in various areas of science such as anomalous diffusion, biology problems, petroleum engineering, heat transfer, physical problems, and others [3,4,33,37]. The fractional-order advection-dispersion equation is solved numerically through different approximation methods (see, for instance, [8,14,15,21,23,29]). Despite the big efforts of the researchers and mathematicians for solving this type of equations, many other researchers still research new techniques and methods that give high accuracy solutions for the same equations. We believe that the proposed shifted Vieta-Lucas orthogonal polynomials are a new method to solve FADE, and no one used this approach to date. Consequently, the main idea behind this work is to introduce a new numerical technique for solving the FADE. This technique depends on the shifted Vieta-Lucas polynomials as a family of the desired numerical solution, besides the nonstandard finite difference method, and the spectral collocation technique. As the final step, the obtained algebraic system of equations will be solved numerically via any iteration method. This paper is structured as follows. In the next section, some necessary mathematical tools required for the construction of this work will be given. In addition, the derivation of some important relations for shifted Vieta-Lucas polynomials will be introduced in the same section. The desired series solution will be presented in Sect. 3. In Sect. 4, the error bound are investigated. Construction methodology for solving the FADE that is given in Eq. (1) via the spectral collocation method, together with the nonstandard method, will be proposed in Sect. 5. For the demonstration of the accuracy, efficiency, and applicability of the suggested technique, some numerical applications will be presented in Sect. 6. Lastly, concluding remarks are reported.

Caputo's differentiation operator
The Caputo operator is linear. Moreover, via Definition 2.1, the following can be claimed: where the ceiling function of r is r .

Nonstandard finite difference scheme notations
The discrete first derivative can be described by: where ψ(h) and φ(h) are functions in the step-size discretization h = t and This formula for the first derivative is called the nonstandard finite difference method presentation. Also, the denominator function satisfies the condition 0 < φ(h) < 1, h → 0. There is no determined base for the best choice of the function φ(h), but we can introduce the most popular functions in the nonstandard finite difference technique as follows: For more details, see some of Mickens's publications [22].

Vieta-Lucas polynomials
In this part of the paper, we study a class of orthogonal polynomials, which, to the best of our knowledge, are presented here for the first time. These polynomials can be created by means of the recurrence relations and analytical formula to build a new family of orthogonal polynomials that will be well-known as Vieta-Lucas polynomials.
Polynomial VL n (x) can be created by means of the following iterative formula: Also, VL n (x) can be obtained through the following explicit power series formula: where n 2 is the ceiling function. Moreover, VL n (x) are orthogonal polynomials with respect to the following integral: where 1 √ 4-x 2 is the weight function corresponding to VL n (x).

Shifted Vieta-Lucas polynomials
In this subsection, the relevant properties and relations of Vieta-Lucas polynomials and their shifts will be concluded. The shifted Vieta-Lucas polynomials (VL * n (x)) can be considered as a new class of orthogonal polynomials defined on the closed interval [0, 1].

Definition 2.3
The shifted Vieta-Lucas polynomials of degree n on [0, 1] can be obtained from VL n (x) as follows: Also, VL * n (x) are created by utilizing the following recurrence formula: with the starting values Moreover, the explicit analytical formula for VL * n (x) can be obtained through the following expression: Polynomials VL * n (x) have the orthogonality property with respect to the following inner product: where ω(x) = 1 √ x-x 2 is the wight function. Let the function u(x) be Lebesgue-square-integrable on the interval [0, 1] and suppose that it can be expressed as a linear combination of the independent power functions in terms of VL * n (x) as follows: where c i are the unknown coefficients. Generally, only the first n + 1 terms of the series in Eq. (14) is appropriate in all cases of the approximation theory. Therefore, we have where the undetermined coefficients c i , i = 0, 2, . . . , n can be obtained via one of the following expressions: or where

Derivation of the main scheme
Theorem 1 Consider the approximate solution of the main problem described in Eq. (1) expressed in the terms of shifted Vieta-Lucas polynomials as in Eq. (15). Then, the fractional-order terms can be transformed into algebraic equations as follows: where Proof Consider the approximate solution for Eq. (1) as given in Eq. (15). Then applying the fractional operator gives Using Caputo's operator definition (7), we have In addition, Using Eqs. (6) and (7), Eq. (23) is reformulated as Collocating Eqs. (21), (22), and (24), we have Here, Eq. (25) can be rewritten in the following form: where The proof is completed.
Remark The second term of fractional-order given in our problem Eq. (1) can be transferred into an algebraic term by using the same theorem (Theorem (1)).

Discussion of error estimate
where M is a constant. Then u(x) can be expressed as an infinite linear combination of shifted Vieta-Lucas polynomials, and u n (x) has only n + 1 terms of this expression. Also, this numerical solution converges uniformly to the function u(x) (u n (x) → u(x) as n → ∞). Moreover, the coefficients given in Eq. (17) are bounded, i.e., where Proof Consider the function u(x) satisfying the stated conditions and having expression as in Eq. (28). Use the approximation theory to take only n + 1 terms of this series as follows: where the coefficients c i , i = 0, 1, . . . , n can be determined via Eq. (17). To compute the integrals, use the substitution 4x -2 = 2 cos(θ ) in Eq. (17), and then the following is gained: Integrating twice by parts in Eq. (31) gives where Applying the absolute value properties on Eq. (32) produces From the properties of the function u(x) and of the trigonometric functions, Eq. (34) is transformed into the inequality By rearranging the integration and performing some analytical manipulations, we have Hence, the proof is completed. • c n is convergent, and R n = ∞ k=n+1 c k . Then Lemma 4.1 is necessary for the following theorem.
Proof The error in L 2 ω [0, 1] is defined by Using Eqs. (14) and (15), in addition to the orthogonality property of VL * i (x) that was introduced in Eq. (13), we have the following result: By Eq. (40) of Theorem 2, we acquire the following inequality: An application of Lemma 4.1 leads to Finally, Eq. (42) describes the error bound as the following formula: The theorem is proved.

Theorem 4
If the function u(x) ∈ [0, 1] is n times continuously differentiable and the most square-suitable approximation for this function is u n (x) defined in Eq. (15), then we have where Proof The function u(x) can be expanded in the following series: Then If u n (x) is as in Theorem (4), then Since Hence, The proof is complete.

Suggested technique for the main problem (FADE)
This section concerns the numerical solution scheme for the fractional-order advectiondispersion equation via our proposed technique in this work. Therefore, assume the problem given in Eqs.
(1)- (3). Firstly, propose the approximate solution of this problem presented in terms of shifted Vieta-Lucas polynomials as follows: Based on Theorem 1, along with substituting Eq. (52) into Eq. (1), we have the following equation: Now, collocate Eq. (53) at the roots of the function VL * n+1α (x), which are named the collocation points x p , as follows: where c i (i = 0, 1, . . . , n) are n+1 unknown coefficients. For obtaining these coefficients n+1 algebraic equations are needed, but Eq. (54) gives us only n + 1α ordinary differential equations that are transformed into algebraic equations. So α equations are required for completing the desired system. For this aim, substitute the numerical approach (52)  Lastly, the nonstandard finite difference method will be used for the transformation of Eq. (54) into n + 1α algebraic equations. Hence, combining these results with Eqs. (55) and (56) will constitute a linear system of n + 1 equations. This system can be solved numerically via any suitable procedure for obtaining the unknowns c i (i = 0, 1, . . . , n) and therefore the desired approximate solution u n (x) of FADE can be calculated.

Numerical experiments
In this section, we introduce some examples that are solved numerically via the suggested method. Also, comparisons of our results with other published numerical results will be made. These examples and numerical benchmarks will be constructed to illustrate the accuracy, applicability, and efficiency of our proposed method for solving the FADE numerically.
Example 6.1 Consider the following fractional-order advection-dispersion equation [15,30,31,34]: with the source function the initial condition and the boundary conditions The analytic solution of Eqs. (58)-(60) is given by Apply our proposed method as follows; Use Eq. (61) together with (54) to obtain where the collocation points x p are the roots for the shifted Vieta-Lucas polynomials of the second degree. Also, Eqs. (55) and (56) are used respectively to have The nonstandard finite difference method given in Sect. 2.2 is applied here to convert the two generated equations of Eq. (62) into two algebraic equations. These equations are combined with Eqs. (63) and (64) to have a system of algebraic equations. Now, we have four equations in four unknown coefficients c i (i = 0, 1, 2, 3). This system can be presented in the following matrix formula: where The numerical results obtained through our suggested method for the first Example 6.1 in different cases are reported in Table 1 in addition to Figs. 1, 2, 3, and 4. In Table 1, Table 1 For Example 6.1, a comparison the absolute error results in [15] and [31], in the case α = 2, β = 1, T = 2 at n = 5, and those in our setting when n = 3, φ(h) = 0.5(exp(2h) -1)  the numerical results were obtained in published articles [15] and [31], and are compared with ours. These numerical results are reported when the parameters are α = 2, β = 1, T = 2 at n = 5 for [15] and [31] Table 1, and the results in Figs. 1, 2, 3, and 4, the suggested method is applicable and more accurate than those in [15] and [31]. Also, it is obvious that our proposed method can be used for large enough times of the problem and gives good accuracy.
Finally, this discussion supported and proved the truth of the theoretical technique and theorems introduced in this article.

Conclusions
In the present article, we have introduced a trustworthy method for solving a mathematical physics model of fractional-order, advection-dispersion equation, numerically. This method was based on a class of orthogonal polynomials which is called the shifted Vieta-Lucas polynomials. The principle strategy of the suggested treatment is transforming the original equation into a system of ODEs by applying the spectral collocation method. Next, we used the nonstandard finite difference method to transform these equations into algebraic equations. After that, the Gaussian elimination method was used (or any appropriate iterative technique). Actually, we have proposed an accurate technique that converges fast to the real solution for solving the FADEs. Also, an error estimate of our method was derived. Some numerical applications were introduced for demonstrating the accuracy, efficiency, and applicability of the current method. Moreover, our numerical results were compared with other published approximation data in the literature, showing reliability and accuracy of the proposed technique. Furthermore, this method gives us flexibility for solving different equations appearing in applications in many fields such as physics, engineering, chemistry, and others, because the given problems are solved at different values of T. All computed results were obtained using MATLAB program. In the future, we can apply this technique for solving many equations from applications in one and two dimensions, especially in physics, such as diffusion equation, wave equation, telegraph equation, and so on.