Numerical Solution of Advection–Diffusion Equation of Fractional Order Using Chebyshev Collocation Method

: This work presents a highly accurate method for the numerical solution of the advection– diffusion equation of fractional order. In our proposed method, we apply the Laplace transform to handle the time-fractional derivative and utilize the Chebyshev spectral collocation method for spatial discretization. The primary motivation for using the Laplace transform is its ability to avoid the classical time-stepping scheme and overcome the adverse effects of time steps on numerical accuracy and stability. Our method comprises three primary steps: (i) reducing the time-dependent equation to a time-independent equation via the Laplace transform, (ii) employing the Chebyshev spectral collocation method to approximate the solution of the transformed equation, and (iii) numerically inverting the Laplace transform. We discuss the convergence and stability of the method and assess its accuracy and efﬁciency by solving various problems in two dimensions.


Introduction
Fractional calculus (FC) is the generalization of classical calculus, which includes integrals and derivatives in non-integer order [1,2].In recent years, the concept of fractional derivatives has been applied with great success to model various real-life phenomena in many scientific fields [3].Fractional order operators include the history of a physical phenomenon from the initial state to the current state.Therefore, fractional order operators are often applied to model systems that describe the influence of memory effects [4].Fractional order operators have broadened the concept of integer order operators [5].In the last ten years, PDEs of fractional order have attracted the remarkable attention of the research community due to their ability to model various problems, such as neural networks [6], optimal voltage controllers for electronic vehicle charging stations [7], robotics [8], optimal controls for impulsive systems [9], hydrology [10], medical sciences [11], and many others [12].
The time-fractional advection-diffusion equations (TFADEs) have been widely used in many fields of science because of their ability to model memory and non-local properties [13].ADEs combine diffusion and advection terms to describe physical events in which particles, energy, or other physical quantities are moved inside a physical system by two processes, i.e., diffusion and advection.The concentration of substances for mass and heat transfer is represented by the advection and diffusion models.This equation has been used to model heat transfer in draining films [14], flow in porous media [15], and air pollution [16], etc.There have been many attempts to solve the TFADEs accurately.
In [17], a second-order implicit difference method was proposed for solving a time-spacefractional ADE.Mardani et al. [18] studied a meshless method for the solution of ADEs with variable coefficients.A fully implicit finite difference scheme based on cubic B-splines was developed in [19] for solving TFADEs.Abbaszadeh and Amjadian [20] studied the spectral element method for the solution of TFADEs.In the current study, we investigate the numerical solution of two-dimensional TFADEs using the Chebyshev spectral collocation method (CSCM) coupled with the Laplace transform (LT).
Spectral collocation methods have attracted the attention of the research community due to their highly accurate solutions to partial differential equations.They converge exponentially fast as compared to the algebraic convergence rates for FEM and FDM.Spectral methods are global in nature.This means that the optimal accuracy can be obtained with few grid points.For the first time, spectral methods were introduced by Gottlieb and Orszag [21].Their work paved the way for more advanced extensions of the method for a larger variety of problems [22,23].The most important tool for spectral methods is the computation of differentiation matrices, discussed briefly by Welfert [24].The spectral collocation methods have been used by many authors.For example, in [25], the authors discussed the Hermite spectral methods for unbounded domains.The solution of the fractional advection-dispersion equation by the Chebyshev spectral method was incorporated by Sweilam et al. [26].In [27], the spectral Legendre-Chebyshev collocation method for variable coefficients was examined.Tian et al. [28] investigated the polynomial spectral collocation method for space-fractional ADEs.The spectral collocation method and the non-standard finite difference technique for the solution of TFADEs were introduced by [29].
In the aforementioned methods, the time discretization is carried out using the finite difference time-stepping method.The main shortcoming of the time-stepping technique is that it may not always result in a stable solution.A finite difference time-stepping scheme is stable if the errors calculated at a one-time step do not cause the errors to be increased during the calculations.In other words, the time-stepping scheme is stable whether the errors remain constant or decay during the computations.Furthermore, for optimal accuracy, we need to decrease the time step, which results in an increased computational time, and, thus, has low efficiency in simulating the time-fractional problems.To overcome this drawback of the finite difference time-stepping scheme, the LT has been used as one of the best alternative mathematical tools [30,31].However, while using the LT, the main difficulty is its inversion.In the literature, many algorithms have been developed for the numerical inversion of the LT [32][33][34].In the current study, we use Stehfest's method [35] and Talbot's method [32,36] for the numerical inversion of the LT.

Numerical Scheme
In our numerical scheme, a two-dimensional TFADE is taken into account.First, the LT is employed to reduce the TFADE to a time-independent problem in the LT domain.Then, the CSCM is employed to obtain the approximate solution to the transformed problem.Finally, the time domain solution is retrieved from the LT domain solution via the inverse LT.

Two-Dimensional TFADE
We consider a two-dimensional time-fractional advection-diffusion equation of the form with boundary conditions and the initial condition where − → u = (u 1 , u 2 ) is a known vector, λ 1 and λ 2 are arbitrary constants, Q( x, τ), U 0 ( x), h( x, τ) are given continuous functions, Ω ⊂ R 2 is a bounded domain with a smooth boundary ∂Ω, and D α τ represents Caputo's derivative, defined as [3]

Time Discretization
In this section, the LT is used to reduce the TFADE to an equivalent time-independent problem in the LT domain.The LT of the function U( x, τ) is denoted by U( x, s) and is defined as and the LT of Caputo's derivative is defined as and Equations ( 4) and ( 5) can be rearranged as and where and I denotes the identity operator.To obtain the approximate solution of the system defined in ( 6) and (7), first, the linear spatial operator L is discretized using the CSCM.Then, the fully discrete system is solved for each node s in the LT domain.Finally, the numerical inversion of the LT is performed to obtain the solution of the problem defined in Equations ( 1)-(3).

Chebyshev Spectral Collocation Method
In the CSCM, we interpolate the data x j , U(x j ) by using the Lagrange interpolation polynomial (LIP) l j (x) of a degree of at most n [23,37] where l j (x) is the LIP at the point x j (j = 0, 1, . . ., n), given as where U j = U(x j ).In the CSCM, the Chebyshev nodes are selected as the interpolation nodes.The domain [−1, 1] is discretized using the Chebyshev points defined as The approximation of the first derivative ∂x on the Chebyshev points can be obtained as where the entries of the matrix D n are The off-diagonal entries of [D n ] i,j are obtained as where α −1 j = ∏ n i =j (x i − x j ), and the diagonal elements are obtained as Next, the entries of the µth-order differentiation matrix D µ n can be obtained analytically as A more accurate and stable evaluation can be found in [24,38].In [24], the author derived a useful recursion relation for the differentiation matrices, given as The LIP corresponding to the above Chebyshev points are given as where l ij ( xij ) = δ ij , i, j = 0, 1, 2, . . ., n.In terms of x and y, the second derivatives of the LIPs (10) are given as where D 2 n denotes the second-order Chebyshev differentiation matrix.Applying the linear operator L on the LIP at the Chebyshev nodes, { xrs } gives Therefore, the discretized form of the linear differential operator L via the CSCM can be expressed as By using Equation ( 12) in Equation ( 6), we have where we incorporate the boundary conditions in Equation ( 7) by taking the matrix L D , based on all collocation points x, and then replace the rows of L D corresponding to collocation at boundary points with unit vectors that have a one in the position corresponding to the diagonal of L D .Thus, the boundary condition U( x, s) = h( x, s) in ( 7) will be explicitly enforced at this point as soon as we set the right-hand side to the corresponding value of h( x, s) [23].After incorporating the boundary conditions and solving Equation ( 13) for each node s, we obtain the approximate solution U( x, s) in the LT domain.In Section 2.4, we describe the numerical inverse LT method used to evaluate the approximate time domain solution U( x, τ) for the original problem in Equations ( 1)-( 3).

Stability and Error Analysis
Given the Chebyshev nodes in Equation ( 9) and the Lagrange interpolation polynomials defined in Equation ( 8), the interpolation operator is defined as follows [37]: The steps proposed by Borm et al. [39] for constructing the error bound are utilized here.Let M n be a constant with the stability estimate given as For the Chebyshev interpolation, we have Thus, the stability constant depends on n and develops extremely slowly.Additionally, the following approximation error bound [39] holds for a function Theorem 1 ([39]).If (15) and (18) hold for U ∈ C (n+1) [−1, 1], and considering k ∈ {0, 1, 2, . . ., n}, then depending on the stability constant Now, use ( 18) and ( 19) for the TFADE (1) in one dimension.For the 1 − D case, the linear differential operator L is of the form L = λ 1 . Thus, we find the following error estimate The time derivative is accurately evaluated.Therefore, the error bound of (D α t (U − I n (U))) ∞ is the same order of (U − I n (U)) ∞ .Finally, we have [37]: where C is constant, resulting from calculation of the coefficients of U (n+1) ∞ .We can use the tensor product interpolation operators for two-dimensional problems.

Numerical Inversion of Laplace Transform
In this section, we implement the inverse Laplace transform method to convert the CSCM solution U( x, s) from the Laplace domain to the time domain as follows: Here the transform U( x, s) needs to be inverted, ρ 0 is the converging abscissa, and ρ > ρ 0 .This means that all of the singularities of U( x, s) lie in the open half-plane Re(s) < ρ.We aim to approximate the integral defined in Equation (20).In most cases, it is quite difficult to evaluate the integral defined in Equation (20) analytically; thus, a numerical method must be used.There are several numerical algorithms in the literature that can be used to evaluate the integral defined in Equation (20).Among them, we can list the Fourier series method [40], the de Hoog method [41], Stehfest's method [35,42], Talbot's method [32,36], and the Weeks method [43,44], etc.Each approach has an identifiable use and is appropriate for a particular function.In this article, we use two popular inversion algorithms: the improved Talbot's and Stehfest's methods will be presented in the following sections.

TheImproved Talbot's Method
Here, we use Talbot's method to approximate the U( x, τ) where C is a suitably chosen line.It is very difficult to solve the above integral due to the highly oscillatory exponential function e sτ and the slow-decaying transform U( x, s).
Talbot [36] suggested that by utilizing contour integration, this problem can be resolved.He further said that the integration might be carried out so that the real parts of the contours begin and stop on the left side of the complex plane, which contains all of the singularities in the transformed function U( x, s).Due to the exponential element, the integrand decays quickly on such contours, making integral Equation ( 21) suitable for approximation with the midpoint rule.Consider a contour of the form [32]: where Res(±π) = −∞, and s(η) is given as where µ, , and γ will be chosen precisely.Using Equation (23) in Equation ( 21) , we obtain The integral defined in ( 24) is approximated via the midpoint rule with spacing h = 2π M T as

Convergence
The approximate value of the defined integral in Equation ( 24) depends on the contour of integration Γ; different rates of convergence are attained for the suggested numerical scheme.Additionally, the quadrature step h determines that the suggested scheme will converge.In order to have optimal results, we need to search for the optimal contour of integration, which can be found by using the optimal values for the parameters involved in (23).In [32], the authors have proposed optimal values for the parameters as δ = 0.6122, , γ = 0.2645, µ = 0.6407, and = 0.5017, For the improved Talbot's method, the error estimate is given as

TheStehfest's Method
The Gaver-Stehfest method is one of the most important techniques for the numerical inversion of the Laplace transform.It was designed in the second half of the 1960s.It has gained prominence in a number of disciplines, such as chemistry, finance, economics, and computational physics, due to its effectiveness and ease.The series of Gaver approximants, as determined by Gaver [42], is the foundation of the Gaver-Stehfest method.Acceleration was required because the Gaver approximants' convergence was basically logarithmic.Stehfest [35] used the Salzer acceleration method to offer a linear acceleration method.The Gaver-Stehfest approach uses a series of functions to approximate U( x, τ) as where the coefficients α j are defined by Solve Equation ( 6) for the corresponding Laplace parameters s = ln2 τ j, j = 1, 2, 3, . . ., M S .The solution to the given problem in Equation ( 1) can be obtained by (26).There are a few appealing features of the Gaver-Stehfest algorithm: (i) U App ( x, τ) are linear in terms of the values of U( x, s); (ii) the values of U( x, s) are needed only for the real value of s; (iii) the process of computing the coefficients is quite simple; and (iv) for constant functions, this approach yields highly accurate approximations.In the literature, this technique has been used by many authors in [45,46], where it has been demonstrated that this method converges very fast to U App ( x, τ) (provided U App ( x, τ) is non-oscillatory).

Convergence
In [45], the author has derived two conditions for U( x, τ), which guarantee the convergence of U App ( x, τ).The conditions are given in the following theorem.
Theorem 2. Let U : (0, ∞) → R be a locally integrable function and the approximate solution U App ( x, τ) be defined by (26), and let its LT U( x, s) be defined for s > 0 : 1.
U App ( x, τ) converges for the values of U App ( x, τ) in the neighborhood of τ.

2.
Let for some real number m and some Then, U → m as M S → +∞.

3.
Allow that U App ( x, τ) has a bounded variation in the neighborhood of τ.Then, Corollary 1.Under the above assumptions, if ∀ ξ and some υ in the neighborhood of τ, then U App ( x, τ) → U( x, τ), as M S → +∞.
Additionally, the authors in [47] performed various experiments for the parameter effect on the accuracy of the numerical scheme, and they concluded in their findings that "If p significant digits are required, then let M S be the positive integer 2.2p .Set the system precision to q = 1.1M S and, for a given M S , calculate α i , 1 ≤ i ≤ M S , using (27).Then, for the given transform U( x, s) and the argument τ, calculate the U App ( x, τ) in (26)".According to these conclusions, if the error of input data is 10 −(q+1) ≤ U−U U ≤ 10 −q , with an even positive integer M S via q = 1.1M S , then the final error is 10 −(p+1) ≤ U−U U ≤ 10 −p , where M S = 2.2p [48].Next, we present Algorithm 1 for the proposed numerical scheme.
Algorithm 1 : Algorithm for the proposed numerical scheme 1: Input: The computational domain, the fractional order derivative, the final time, the optimal parameters for the inverse laplace transform methods, the inhomogeneous function, and other conditions.2: Step i: Apply the Laplace transform to problems (1) and (2), and obtain the timeindependent problems ( 6) and ( 7).

3:
Step ii: Obtain L D by discretizing the linear differential operator L using the CSCM via (12).

4:
Step iii: Solve Equation ( 13) with the boundary conditions given in Equation ( 7) for each point in the LT domain.

Numerical Results and Discussion
The numerical results and discussions are addressed in this section.Three examples are used to evaluate and validate the effectiveness and accuracy of the Laplace-transformed CSCM.We performed our experiments in MATLAB R2019a on a Windows 10 (64-bit) PC equipped with an Intel(R) Core TM i5-4310M CPU @ 2.70 GHz with 4 GB of RAM.Three error norms-the relative L 2 error, the maximum absolute error, and the RMS error-are used to evaluate the accuracy of the proposed method and are defined below: where U( xj , τ) and U App ( xj , τ) represent the analytic and numerical solutions, respectively.
Example 1.Consider the 2D time-fractional advection-diffusion Equation (1) with an exact solution given as where λ 1 = 1, − → u = (1, 1) T , and λ 2 = 0.The forcing term Q(x, y, τ), and the initial-boundary data are extracted from the exact solution.The numerical results obtained via Talbot's and Stehfest's methods are presented in Tables 1 and 2 by using different values of α, M S , M T , and n at τ = 1.
The approximate solution obtained by the proposed method is shown in Figure 1a.The plots for the comparison of errors er 2 , er ∞ , and er rms using Stehfest's method for several values of α, M S , and τ at n = 25 are shown in Figure 1b-d, respectively.Figure 2a-c show the plots of errors er 2 , er ∞ , and er rms obtained using Talbot's method with different values of α, M T , and τ at n = 25.Figure 2d presents the contour plot of the absolute error using Talbot's method for different values of α with M T = 26, τ = 1, and n = 25.The accuracy of Stehfest's method enhances for M S = 8, 10, 12, 14 and then gradually diminishes for values exceeding M S = 14.From all of the obtained results presented in the figures and tables, we conclude that the proposed method is stable, accurate, and efficient.Clearly, the numerical results show that the accuracy of Talbot's method is better than that of Stehfest's method.
Example 2. Consider the 2D time-fractional advection-diffusion Equation (1) with an exact solution given as where λ 1 = 1, − → u = (1, 1) T , and λ 2 = 0.The forcing term Q(x, y, τ), and the initial-boundary data are extracted from the exact solution.The numerical results obtained via Talbot's and Stehfest's methods are presented in Tables 3 and 4 by using different values for α, M S , M T , and n at τ = 1.
The approximate solution obtained by the proposed method is shown in Figure 3a.The plots of the comparison of errors er 2 , er ∞ , and er rms using Stehfest's method for several values of α, M S , and τ at n = 25 are shown in Figure 3b-d, respectively.Figure 4a-c show the plots of errors er 2 , er ∞ , and er rms by using Talbot's method with different values of α, M T , and τ at n = 25.Figure 3e presents the contour plot of the absolute error obtained by using Stehfest's method for different values of α with M S = 12, τ = 1, and n = 25.The accuracy of Stehfest's method enhances for M S = 8, 10, 12, 14 and then gradually diminishes for values exceeding M S = 14.From all of the obtained results presented in the figures and tables, we conclude that the proposed method is stable, accurate, and efficient.Clearly, the results show that the accuracy of Talbot's method is better than that of Stehfest's method.Example 3. Consider the 2D time-fractional advection-diffusion Equation (1) with an exact solution given as where λ 1 = 1, − → u = (1, 1) T , and λ 2 = 0.The forcing term Q(x, y, τ), and the initial-boundary data are extracted from the exact solution.The numerical results obtained via Talbot's and Stehfest's methods are presented in Tables 5 and 6 by using different values of α, M S , andM T , with n at τ = 1.The approximate solution obtained by the proposed method is shown in Figure 5a.The plots of the comparison of the errors er 2 , er ∞ , and er rms using Stehfest's method for several values of α, M S , and τ at n = 25 are shown in Figure 5b-d, respectively.Figure 6a-c show the plots of errors er 2 , er ∞ , and er rms via Talbot's method with different values of α, M T , and τ at n = 25.Figure 6d presents the contour plot of the absolute error by using Talbot's method for different values of α with M T = 26, τ = 1, and n = 25.The accuracy of Stehfest's method enhances for M S = 8, 10, 12, 14 and then gradually diminishes for values exceeding than M S = 14.From all the obtained results presented in the figures and tables, we conclude that the proposed method is stable, accurate, and efficient.Clearly, the results show that the accuracy of Talbot's method is better than that of Stehfest's method.

Conclusions
In this article, a numerical method based on the LT coupled with the CSCM for the numerical solution of two-dimensional TFADEs in the Caputo sense has been developed.It first uses the LT to convert a TFADE into a time-independent inhomogeneous equation in the LT space.Then, it uses the CSCM to discretize the spatial derivatives of this Laplacetransformed inhomogeneous equation.Finally, it chooses Stehfest's method and Talbot's method for the numerical inversion of the LT to retrieve the numerical solutions of the TFADE from the corresponding CSCM solutions in the LT domain.Compared with the finite difference time-stepping scheme, the proposed method introduces the LT technique to avoid the calculation of a costly convolution integral in the approximation of time-fractional derivation and to avoid the effect of the time step on the stability and accuracy.Hence, it can successfully solve the long time history TFADEs.The proposed advection-diffusion equation of fractional order presents significant potential for extension into uncertainty propagation and sensitivity analysis, enabling a more comprehensive and nuanced understanding of system responses to parameter variations and initial conditions [49][50][51].

Figure 1 .
Figure 1.(a) Numerical solution of Example 1.(b) The plot shows a comparison of er 2 , er ∞ , and er rms using Stehfest's method for different α, with M S = 10 and n = 25 at τ = 1 for Example 1. (c) The plot shows a comparison of er 2 , er ∞ , and er rms using Stehfest's method for different M S with α = 0.3 and n = 25 at τ = 1 for Example 1.(d) The plot shows a comparison of er 2 , er ∞ , and er rms using Stehfest's method for different τ with M S = 10, α = 0.3, and n = 25 for Example 1.

Figure 2 .
Figure 2. (a) The plot shows a comparison of er 2 , er ∞ , and er rms using Talbot's method for different α with M T = 22 and n = 25 at τ = 1 for Example 1.(b) The plot shows a comparison of er 2 , er ∞ , and er rms using Talbot's method for different M T with α = 0.3 and n = 25 at τ = 1 for Example 1. (c) The plot shows a comparison of er 2 , er ∞ , and er rms using Talbot's method for different τ with M T = 22, al pha = 0.3, and n = 25 for Example 1.(d) Contour plot of the absolute error obtained with M T = 26, τ = 1 al pha = 0.3, and n = 25 using Talbot's method for Example 1.

Figure 4 .
Figure 4. (a) The plot shows a comparison of er 2 , er ∞ , and er rms using Talbot's method for different α, with M T = 22 and n = 25 at τ = 1 for Example 2. (b) The plot shows a comparison of er 2 , er ∞ , and er rms using Talbot's method for different M T with α = 0.3 and n = 25 at τ = 1 for Example 2. (c) The plot shows a comparison of er 2 , er ∞ , and er rms using Talbot's method for different τ with M T = 22, al pha = 0.3, and n = 25 for Example 2.

Figure 6 .
Figure 6.(a) The plot shows a comparison of er 2 , er ∞ , and er rms using Talbot's method for different α, with M T = 26 and n = 25 at τ = 1 for Example 3. (b) The plot shows a comparison of er 2 , er ∞ , and er rms using Talbot's method for different M T with α = 0.3 and n = 25 at τ = 1 for Example 3. (c) The plot shows a comparison of er 2 , er ∞ , and er rms using Talbot's method for different τ with M T = 26, al pha = 0.3, and n = 25 for Example 3. (d) Contour plot of the absolute error obtained with M T = 26, τ = 1 al pha = 0.5, and n = 25 using Talbot's method for Example 3.

Table 1 .
The er 2 , er ∞ , and er rms obtained using Stehfest's method for different values of α, n, and M S at τ = 1 for Example 1.

Table 2 .
The er 2 , er ∞ , and er rms obtained via Talbot's method for different values of α, n, and M T at τ = 1 for Example 1.

Table 3 .
The er 2 , er ∞ , ander rms obtained via Stehfest's method for different values of α, n, and M S at τ = 1 of Example 2.

Table 4 .
The er 2 , er ∞ , ander rms obtained using Talbot's method for different values of α, n, and M T at τ = 1 of Example 2.

Table 5 .
The er 2 , er ∞ , ander rms obtained via Stehfest's method for different values of α, n, and M S at τ = 1 of Example 3.

Table 6 .
The er 2 , er ∞ , ander rms obtained using Talbot's method for different values of α, n, and M T at τ = 1 for Example 3.