A New Numerical Algorithm for Solving a Class of Fractional Advection-Dispersion Equation with Variable Coefficients Using Jacobi Polynomials

Abstract

We propose Jacobi-Gauss-Lobatto collocation approximation for the numerical solution of a class of fractional-in-space advection-dispersion equation with variable coefficients based on Caputo derivative. This approach has the advantage of transforming the problem into the solution of a system of ordinary differential equations in time this system is approximated through an implicit iterative method. In addition, some of the known spectral collocation approximations can be derived as special cases from our algorithm if we suitably choose the corresponding special cases of Jacobi parameters 𝛼 and 𝛽. Finally, numerical results are provided to demonstrate the effectiveness of the proposed spectral algorithms.

1. Introduction

Spectral methods have emerged as powerful techniques used in applied mathematics and scientific computing to numerically solve differential equations [1]. Also, they have became increasingly popular for solving fractional differential equations (see, for instance, [2–6]). The main idea of spectral methods is to put the solution of the problem as a sum of certain basis functions and then to choose the coefficients in the sum in order to minimize the difference between the exact solution and approximate one as well as possible. Spectral collocation method has an exponential convergence rate, which is valuable in providing highly accurate solutions to nonlinear differential equations even using a small number of grids. Moreover, the choice of collocation points is very useful for the convergence and efficiency of the collocation approximation [7, 8].

In recent years, considerable interest in fractional partial differential equations has been motivated because of their growing applications in electromagnetics, acoustics, viscoelasticity, electrochemistry, and material science [9, 10]. Several analytical algorithms have been investigated for treating these equations analytically to obtain closed-form solutions such as variational iteration method, Fourier transform method, homotopy analysis method, the method of separation of variables, Adomian decomposition method, and Laplace transform method [9, 11–14]. However, there are only few types of these equations in which the analytical solutions are available. Therefore, numerical means have to be used in general.

In numerous physical models, an equation commonly used to describe transport diffusive problems is the classical advection-diffusion (or -dispersion) equation which may be generalized to the fractional ones to cover other very interesting physical models. The advection-dispersion equation, which is based on Fick's law, is commonly used to simulate contaminant transport in porous media [15]. The space, time and time-space fractional advection-dispersion equations are presented as a reliable model to simulate the transport of passive tracers carried by fluid flow in a porous media and are used in groundwater hydrology [13, 16]. Moreover, they have been introduced to describe other important physical phenomena (see [17–21]).

In the last few years, theory and numerical analysis of fractional partial differential equations have received an increasing attention. In this direction, Rihan [22] proposed the 𝜃-method for approximating time-fractional parabolic partial differential equations in the Caputo sense. An explicit Euler method, an implicit Euler method, and the fractional Cranck-Nicholson method for solving fractional differential equations are discussed in [23–26]. An explicit difference approach for solving space fractional diffusion equation has been proposed in [27]. Ding et al. [28] investigated a class of weighted finite difference method for tackling a class of time-dependent fractional differential equations based on shifted Grünwald formula. K. Wang and H. Wang [29] and Huang et al. [16] proposed a fast numerical scheme for fractional time-dependent advection-diffusion and advection-dispersion equations based on finite difference method, respectively. Recently, the Sinc-Legendre spectral method has been developed in [30] for the fractional convection-diffusion. Jiang and Lin [31] proposed a new method for a class of fractional advection-dispersion in the reproducing kernel space. Furthermore, Liu et al. [32] proposed an efficient implicit numerical method for a class of fractional advection-dispersion models in which they discussed five fractional models. In the area of numerical methods of fractional partial differential equations, little work has been done by spectral methods compared to finite difference and finite element methods. This partially motivates our interest in such methods.

The main purpose of the this paper is to construct the solution of a class of space fractional advection-dispersion equation with variable coefficients using Jacobi-Gauss-Lobatto collocation (J-GL-C) approximation, based on Jacobi-Gauss-Lobatto quadrature knots, combined with an implicit iterative method for treating the time discretization. More precisely, implementing the J-GL-C approximation to the spatial variable of the fractional advection-dispersion equation and the corresponding boundary conditions reduces the problem to the time integration of a system of ordinary differential equations in respect to the time variable. To the best of our knowledge, such algorithm has not been implemented for solving space fractional initial-boundary problems.

The plan of the paper is as follows. In the next section, we introduce basic properties of Jacobi polynomials. In Section 3, the way of constructing the Gauss-Lobatto collocation technique for space fractional advection-dispersion equation with variable coefficients is described using the Jacobi polynomials, and in Section 4 the proposed method is applied to two problems. Finally, some concluding remarks are given in Section 5.

2. Preliminaries

In this section, we give some definitions and properties of the fractional calculus (see, e.g., [9, 18, 33, 34]) and Jacobi polynomials (see, e.g., [35–37]).

For 𝑚 to be the smallest integer that exceeds 𝜈, Caputo's fractional derivative operator of order 𝜈>0 is defined as 𝜕^𝜈_𝑥𝑓(𝑥)=𝐽^𝑚−𝜈𝐷^𝑚𝑓(𝑥)𝐷,if𝑚−1<𝜈<𝑚,_𝑚𝑓(𝑥),if𝜈=𝑚,𝑚∈𝑁, (1) where 𝐽^𝜈 is the Riemann-Liouville fractional integral operator of order 𝜈(𝜈≥0) and is defined as 𝐽^𝜈𝑓(𝑥)=1Γ(𝜈)^𝑥₀(𝑥−𝑡)^𝜈−1𝑓(𝑡)𝑑𝑡,𝜈>0,𝑥>0. (2)

Form (1), the Caputo fractional derivative of 𝑥^𝛽 is given by 𝜕^𝜈_𝑥𝑥^𝛽=0,for𝛽∈𝑁₀⌈𝜈⌉,Γ,𝛽<𝛽+1Γ𝑥𝛽+1−𝜈^𝛽−𝜈,for𝛽∈𝑁₀⌈𝜈⌉,⌊𝜈⌋,,𝛽≥or𝛽∉𝑁,𝛽> (3) where ⌈𝜈⌉ and ⌊𝜈⌋ are ceiling and floor functions. Also, 𝑁={1,2,…} and 𝑁₀={0,1,2,…}. Caputo's fractional differentiation is a linear operation; that is, 𝜕^𝜈_𝑥𝜆𝑓(𝑥)+𝜇𝑔(𝑥)=𝜆𝜕^𝜈_𝑥𝑓(𝑥)+𝜇𝜕^𝜈_𝑥𝑔(𝑥), (4) where 𝜆 and 𝜇 are constants.

Let 𝛼>−1, 𝛽>−1, and 𝑃_𝑘^(𝛼,𝛽)(𝑥) be the standard Jacobi polynomial of degree 𝑘. We have that 𝑃_𝑘^(𝛼,𝛽)(−𝑥)=(−1)^𝑘𝑃_𝑘^(𝛼,𝛽)(𝑥),𝑃_𝑘^(𝛼,𝛽)()=()−1−1^𝑘Γ𝑘+𝛽+1.𝑘!Γ𝛽+1 (5)

For integer 𝑚, the 𝑚th-order derivative of Jacobi polynomials is 𝐷^𝑚𝑃_𝑘^(𝛼,𝛽)(𝑥)=2^−𝑚Γ𝑚+𝑘+𝛼+𝛽+1Γ𝑃𝑘+𝛼+𝛽+1^{(𝛼+𝑚,𝛽+𝑚)𝑘−𝑚}(𝑥). (6) Let 𝑤^(𝛼,𝛽)(𝑥)=(1−𝑥)^𝛼(1+𝑥)^𝛽; then we define the weighted space 𝐿²_{𝑤(𝛼,𝛽)}(−1,1) as usual, equipped with the following inner product and norm: (𝑢,𝑣)_{𝑤(𝛼,𝛽)}=¹₋₁𝑢(𝑥)𝑣(𝑥)𝑤^(𝛼,𝛽)(𝑥)‖𝑣‖𝑑𝑥,_{𝑤(𝛼,𝛽)}=(𝑣,𝑣)_𝑤^1/2_(𝛼,𝛽). (7) The set of Jacobi polynomials forms a complete 𝐿²_{𝑤𝛼,𝛽}(−1,1)-orthogonal system, and 𝑃_𝑘^(𝛼,𝛽)2_{𝑤(𝛼,𝛽)}=ℎ_𝑘^(𝛼,𝛽)=2^𝛼+𝛽+1Γ(𝑘+𝛼+1)Γ𝑘+𝛽+1Γ()Γ.2𝑘+𝛼+𝛽+1𝑘+1𝑘+𝛼+𝛽+1 (8)

Let 𝐿>0; then the shifted Jacobi polynomial of degree 𝑘 on the interval (0,𝐿) is defined by 𝑃^(𝛼,𝛽)_𝐿,𝑘(𝑥)=𝑃_𝑘^(𝛼,𝛽)((2𝑥/𝐿)−1).

With the aid of (5), we demonstrate that 𝑃^(𝛼,𝛽)_𝐿,𝑗(0)=(−1)^𝑗Γ𝑗+𝛽+1Γ𝛽+1𝑗!. (9)

Next, let 𝑤_𝐿^(𝛼,𝛽)(𝑥)=(𝐿−𝑥)^𝛼𝑥^𝛽; then we define the weighted space 𝐿²_{𝑤𝐿(𝛼,𝛽)}(0,𝐿) in the usual way, with the following inner product and norm: (𝑢,𝑣)_𝑤𝐿^(𝛼,𝛽)=^𝐿₀𝑢(𝑥)𝑣(𝑥)𝑤_𝐿^(𝛼,𝛽)(𝑥)‖𝑣‖𝑑𝑥,_{𝑤𝐿(𝛼,𝛽)}=(𝑣,𝑣)_𝑤^1/2_{𝐿(𝛼,𝛽)}. (10)

The set of shifted Jacobi polynomials is a complete 𝐿²_{𝑤𝐿(𝛼,𝛽)}(0,𝐿)-orthogonal system. Moreover, due to (8), we have 𝑃^(𝛼,𝛽)_𝐿,𝑘²_{𝑤𝐿(𝛼,𝛽)}=𝐿2^𝛼+𝛽+1ℎ_𝑘^(𝛼,𝛽)=ℎ^(𝛼,𝛽)_𝐿,𝑘. (11)

3. Jacobi Spectral Collocation Method

Since the Jacobi spectral collocation method approximates the initial-boundary problems in physical space and it is a global method, it is very easy to implement and adapt to various problems, including variable coefficient and nonlinear problems (see, for instance, [7, 38]). In this section, a new algorithm for solving time-dependent space fractional advection-dispersion equation is proposed based on Jacobi-Gauss-Lobatto spectral collocation approximation and an implicit iterative method in finite space-time domain.

In this section, we consider the space fractional advection-dispersion equations with space and time variable coefficients [19, 23, 28, 39, 40]: 𝜕𝑢(𝑥,𝑡)()()𝜕𝑡+𝑎𝑥,𝑡𝜕𝑢𝑥,𝑡()𝜕𝜕𝑥−𝑏𝑥,𝑡^𝜈𝑢(𝑥,𝑡)𝜕𝑥^𝜈=𝑞(𝑥,𝑡),0≤𝑥≤𝐿,0<𝑡≤𝑇,1<𝜈≤2, (12) subject to the boundary conditions: 𝑢(0,𝑡)=𝑔₀(𝑥)𝑢(),0<𝑡≤𝑇,𝐿,𝑡=𝑔₁(𝑥),0<𝑡≤𝑇, (13) and the initial value: 𝑢(𝑥,0)=𝑓(𝑥),0≤𝑥≤𝐿, (14) where 𝑎 is the drift of the process, that is, the mean advective velocity, 𝑏≥0 is the coefficient of dispersion, 𝜈 is the fractional order in the Caputo sense, 1<𝜈≤2, and 𝑞 is a source function. In particular, if 𝜈=2, (12) is the classical convection-diffusion equation with a source term which has commonly been used to describe the Brownian motion of particles [41]. Moreover, if 𝑎=0, it reduces to the space fractional diffusion equation (cf. [42–46]).

Now we introduce the Jacobi-Gauss-Lobatto quadratures in two different intervals (−1,1), and (0,𝐿). Denote by 𝑥^(𝛼,𝛽)_𝑁,𝑗(𝑥^(𝛼,𝛽)_{𝐿,𝑁,𝑗}), 0⩽𝑗⩽𝑁, and 𝜛^(𝛼,𝛽)_𝑁,𝑗(𝜛^(𝛼,𝛽)_{𝐿,𝑁,𝑗}), (0≤𝑖≤𝑁), the nodes and Christoffel numbers of the standard (shifted) Jacobi-Gauss-Lobatto quadratures on the intervals (−1,1) and (0,𝐿), respectively. Then one can clearly deduce that 𝑥^(𝛼,𝛽)_{𝐿,𝑁,𝑗}=𝐿2𝑥^(𝛼,𝛽)_𝑁,𝑗𝜛+1,0≤𝑗≤𝑁,_{(𝛼,𝛽)𝐿,𝑁,𝑗}=𝐿2^𝛼+𝛽+1𝜛^(𝛼,𝛽)_𝑁,𝑗,0≤𝑗≤𝑁, (15) and if 𝑆_𝑁(0,𝐿) denotes the set of all polynomials of degree at most 𝑁, then it follows that, for any 𝜙∈𝑆_2𝑁+1(0,𝐿), we have ^𝐿₀𝑤_𝐿^(𝛼,𝛽)(𝑥)𝜙(𝑥)=𝐿𝑑𝑥2^𝛼+𝛽+11₋₁(1−𝑥)^𝛼(1+𝑥)^𝛽𝐿×𝜙2()=𝐿𝑥+1𝑑𝑥2^{𝑁𝛼+𝛽+1}_𝑗=0𝜛^(𝛼,𝛽)_𝑁,𝑗𝜙𝐿2𝑥^(𝛼,𝛽)_𝑁,𝑗=+1_𝑁𝑗=0𝜛^(𝛼,𝛽)_{𝐿,𝑁,𝑗}𝜙𝑥^(𝛼,𝛽)_{𝐿,𝑁,𝑗}. (16)

We define the discrete inner product and norm as follows: (𝑢,𝑣)_{𝑤𝛼,𝛽,𝑁}=^𝑁_𝑘=0𝑢𝑥^(𝛼,𝛽)_𝑁,𝑘𝑣𝑥^(𝛼,𝛽)_𝑁,𝑘𝜛^(𝛼,𝛽)_𝑁,𝑘,‖𝑢‖_{𝑤𝛼,𝛽,𝑁}=(𝑢,𝑢)_{𝑤𝛼,𝛽,𝑁}. (17) Obviously, (𝑢,𝑣)_{𝑤𝛼,𝛽,𝑁}=(𝑢,𝑣)_{𝑤𝛼,𝛽}∀𝑢𝑣∈𝑆_2𝑁−1. (18) Thus, for any 𝑢∈𝑆_𝑁, the norms ‖𝑢‖_{𝑤𝛼,𝛽,𝑁} and ‖𝑢‖_{𝑤𝛼,𝛽} coincide.

Associating with this quadrature rule, we denote by 𝐼^{𝑃(𝛼,𝛽)}_𝑁 the Jacobi-Gauss-Lobatto interpolation (cf. [47]): 𝐼^{𝑃(𝛼,𝛽)}_𝑁𝑢𝑥^(𝛼,𝛽)_𝑁,𝑘=𝑢𝑥^(𝛼,𝛽)_𝑁,𝑘,0≤𝑗≤𝑁. (19)

We now derive an efficient algorithm for solving space-fractional advection diffusion equation (12)–(14). We expand the numerical approximation in terms of Jacobi polynomials: 𝑢_𝑁(𝑥,𝑡)=^𝑁_𝑗=0𝑎_𝑗(𝑡)𝑃^(𝛼,𝛽)_𝐿,𝑗(𝑥). (20) If we make use of the orthogonality property of Jacobi polynomials with respect to the weight functions 𝑤^𝛼,𝛽 and the discrete inner product (17), then we get 𝑎_𝑗(𝑡)=1ℎ_𝑗^𝑁_𝑖=0𝑃^(𝛼,𝛽)_𝐿,𝑗𝑥^(𝛼,𝛽)_𝑁,𝑖𝜛^(𝛼,𝛽)_𝑁,𝑖𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖,𝑡, (21) and accordingly, (20) takes the following form: 𝑢_𝑁(𝑥,𝑡)=^𝑁_𝑗=01ℎ_𝑗^𝑁_𝑖=0𝑃^(𝛼,𝛽)_𝐿,𝑗𝑥^(𝛼,𝛽)_𝑁,𝑖×𝜛^(𝛼,𝛽)_𝑁,𝑖𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖,𝑡𝑃^(𝛼,𝛽)_𝐿,𝑗(𝑥), (22) or equivalently 𝑢_𝑁(𝑥,𝑡)=^𝑁_𝑖=0^𝑁_𝑗=01ℎ_𝑗𝑃^(𝛼,𝛽)_𝐿,𝑗𝑥^(𝛼,𝛽)_𝑁,𝑖×𝑃^(𝛼,𝛽)_𝐿,𝑗(𝑥)𝜛^(𝛼,𝛽)_𝑁,𝑖𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖,𝑡. (23)

The first order spatial derivative of the spectral solution can be approximated by the J-GL-C points 𝑢_𝑥𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑡=^𝑁_𝑖=0^𝑁_𝑗=01ℎ_𝑗𝑃_𝑗^(𝛼,𝛽)𝑥^(𝛼,𝛽)_𝑁,𝑖×𝑃_𝑗^(𝛼,𝛽)𝑥^(𝛼,𝛽)_𝑁,𝑛^′𝜛^(𝛼,𝛽)_𝑁,𝑖×𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖,𝑡,𝑛=0,1,…,𝑁. (24)

According to 𝐷𝑃_𝑗^(𝛼,𝛽)(𝑥)=Γ𝑗+𝛼+𝛽+2𝑃2Γ𝑗+𝛼+𝛽+1^{(𝛼+1,𝛽+1)𝑗−1}(𝑥), (25) equation (24) can be written in the following form: 𝑢_𝑥𝑥^(𝛼,𝛽)_𝑁,𝑛=,𝑡_{𝑁𝑖=0𝑁𝑗=0}1ℎ_𝑗𝑃_𝑗^(𝛼,𝛽)𝑥^(𝛼,𝛽)_𝑁,𝑖Γ𝑗+𝛼+𝛽+22Γ𝑗+𝛼+𝛽+1×𝑃^{(𝛼+1,𝛽+1)}_𝑗−1𝑥^(𝛼,𝛽)_𝑁,𝑛𝜛^(𝛼,𝛽)_𝑁,𝑖×𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖=,𝑡_𝑁𝑖=0𝐴_𝑛𝑖𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖,𝑡,𝑛=0,1,…,𝑁, (26) where 𝐴_𝑛𝑖=^𝑁_𝑗=0Γ𝑗+𝛼+𝛽+22Γ𝑗+𝛼+𝛽+1ℎ_𝑗×𝑃_𝑗^(𝛼,𝛽)𝑥^(𝛼,𝛽)_𝑁,𝑖𝑃^{(𝛼+1,𝛽+1)}_𝑗−1𝑥^(𝛼,𝛽)_𝑁,𝑛𝜛^(𝛼,𝛽)_𝑁,𝑖. (27)

The fractional derivative of order 𝜈 in the Caputo sense for the Jacobi polynomials is given by 𝜕^𝜈𝜕𝑥^𝜈𝑃^(𝛼,𝛽)_𝐿,𝑖(𝑥)=^∞_𝑗=0𝑆_𝜈𝑖,𝑗,𝛼,𝛽𝑃^(𝛼,𝛽)_𝐿,𝑗(𝑥),𝑖=⌈𝜈⌉,⌈𝜈⌉+1,…, (28) where 𝑆_𝜈𝑖,𝑗,𝛼,𝛽=^𝑖_{⌈𝜈⌉𝑘=}(−1)^𝑖−𝑘𝐿^{𝛼+𝛽−𝜈+1}ΓΓ×ℎ𝑗+𝛽+1×Γ𝑖+𝛽+1𝑖+𝑘+𝛼+𝛽+1_𝑗ΓΓ()!𝑗+𝛼+𝛽+1𝑘+𝛽+1×Γ𝑖+𝛼+𝛽+1×Γ𝑘−𝜈+1)(𝑖−𝑘⁻¹×^𝑗_𝑙=0(−1)^𝑗−𝑙Γ()Γ×Γ×𝑗+𝑙+𝛼+𝛽+1×Γ𝛼+1𝑙+𝑘+𝛽−𝜈+1𝑙+𝛽+1×Γ𝑙+𝑘+𝛼+𝛽−𝜈+2𝑗−𝑙!𝑙!⁻¹. (29)

The spatial partial fractional derivatives of order 𝜈 for the spectral solution (20) can be evaluated at the J-GL-C points {𝑥^(𝛼,𝛽)_𝑁,𝑛;𝑛=0,1,…,𝑁}. Hence, we have 𝜕^𝜈𝜕𝑥^𝜈𝑢_𝑁𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑡=^𝑁_𝑖=0^𝑁_𝑗=01ℎ_𝑗𝑃^(𝛼,𝛽)_𝐿,𝑗𝑥^(𝛼,𝛽)_𝑁,𝑖×𝜕^𝜈𝜕𝑥^𝜈𝑃^(𝛼,𝛽)_𝐿,𝑗𝑥^(𝛼,𝛽)_𝑁,𝑛𝜛^(𝛼,𝛽)_𝑁,𝑖×𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖=,𝑡_𝑁𝑖=0𝐵_𝑛𝑖𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖,𝑡,𝑛=0,1,…,𝑁, (30) where 𝐵_𝑛𝑖=^𝑁_𝑗=01ℎ_𝑗𝑃^(𝛼,𝛽)_𝐿,𝑗𝑥^(𝛼,𝛽)_𝑁,𝑖𝜕^𝜈𝜕𝑥^𝜈𝑃^(𝛼,𝛽)_𝐿,𝑗𝑥^(𝛼,𝛽)_𝑁,𝑛𝜛^(𝛼,𝛽)_𝑁,𝑖=^𝑁_𝑗=01ℎ_𝑗𝑃^(𝛼,𝛽)_𝐿,𝑗𝑥^(𝛼,𝛽)_𝑁,𝑖^∞_𝑙=0𝑆_𝜈𝑗,𝑙,𝛼,𝛽𝑃^(𝛼,𝛽)_𝐿,𝑙𝑥^(𝛼,𝛽)_𝑁,𝑛×𝜛^(𝛼,𝛽)_𝑁,𝑖≃^𝑁𝑁_𝑗=0𝑙=01ℎ_𝑗𝑃^(𝛼,𝛽)_𝐿,𝑗𝑥^(𝛼,𝛽)_𝑁,𝑖𝑆_𝜈𝑗,𝑙,𝛼,𝛽𝑃^(𝛼,𝛽)_𝐿,𝑙𝑥^(𝛼,𝛽)_𝑁,𝑛×𝜛^(𝛼,𝛽)_𝑁,𝑖, (31) for 1<𝜈≤2, and 𝑆_𝜈(𝑗,𝑙,𝛼,𝛽) is defined in (29).

If we apply the Jacobi-Gauss-Lobatto collocation method of (12) without the two assigned abscissas 0 and 𝐿; {𝑥^(𝛼,𝛽)_𝑁,0=0,𝑥^(𝛼,𝛽)_𝑁,𝑁=𝐿}, which will be necessary used as two points from the collocation nodes for enforcing the boundary conditions (13), and using (30), then (12) may be written as 𝑢_𝑡𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑡=−𝑎𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑡^𝑁_𝑖=0𝐴_𝑛𝑖𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖𝑥,𝑡−𝑏^(𝛼,𝛽)_𝑁,𝑛,𝑡^𝑁_𝑖=0𝐵_𝑛𝑖𝑢𝑥^(𝛼,𝛽)_𝑁,𝑖𝑥,𝑡+𝑞^(𝛼,𝛽)_𝑁,𝑛,𝑡,𝑛=1,…,𝑁−1. (32)

Let us denote that 𝑢_𝑛(𝑡)=𝑢_𝑁𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑎,𝑡_𝑛(𝑡)=𝑎𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑏,𝑡_𝑛(𝑡)=𝑏𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑞,𝑡_𝑛(𝑡)=𝑞𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑢̇,𝑡_𝑛(𝑡)=𝑢_𝑡𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑡. (33) Thus, (32) can be rewritten in the following simple form: 𝑢̇_𝑛(𝑡)=−𝑎_𝑛(𝑡)^𝑁_𝑖=0𝐴_𝑛𝑖𝑢_𝑖(𝑡)−𝑏_𝑛(𝑡)^𝑁_𝑖=0𝐵_𝑛𝑖𝑢_𝑖(𝑡)+𝑞_𝑛(𝑡),𝑛=1,…,𝑁−1. (34) Let us assume that 𝑑_𝑛(𝑡)=𝑎_𝑛(𝑡)𝐴_𝑛0𝑔₁(𝑡)+𝐴_𝑛𝑁𝑔₂(𝑡),𝑒_𝑛(𝑡)=𝑏_𝑛(𝑡)𝐵_𝑛0𝑔₁(𝑡)+𝐵_𝑛𝑁𝑔₂(𝑡), (35) Then (34) and using the two-point boundary conditions (13) generate a system of (𝑁−1) ordinary differential equations in time. 𝑢̇_𝑛(𝑡)=−𝑎_𝑛(𝑡)^𝑁−1_𝑖=1𝐴_𝑛𝑖𝑢_𝑖(𝑡)−𝑏_𝑛(𝑡)^𝑁−1_𝑖=1𝐵_𝑛𝑖𝑢_𝑖(𝑡)−𝑑_𝑛(𝑡)−𝑒_𝑛(𝑡)+𝑞_𝑛(𝑡),𝑛=1,…,𝑁−1. (36) with the initial values 𝑢_𝑛(0)=𝑓𝑥^(𝛼,𝛽)_𝑁,𝑛,𝑛=1,…,𝑁−1, (37) which may be written in the following matrix form: 𝐮̇(𝑡)((𝑡,𝐮(0)=𝐅𝑡,𝑢))=𝐟, (38) where 𝐮̇(𝑡)=𝑢̇₁(𝑡),𝑢̇₂(𝑡),…,𝑢̇_𝑁−1(𝑡)^𝑇,𝐟=𝑓𝑥_𝑁,1,𝑓𝑥_𝑁,2,…,𝑓𝑥_{𝑁,𝑁−1},𝐅(𝑡,𝑢(𝑡))=𝐹₁(𝑡,𝑢(𝑡)),𝐹₂(𝑡,𝑢(𝑡)),…,𝐹_𝑁−1(𝑡,𝑢(𝑡))^𝑇,𝐹_𝑛(𝑡,𝑢(𝑡))=−𝑎_𝑛(𝑡)^𝑁−1_𝑖=1𝐴_𝑛𝑖𝑢_𝑖(𝑡)−𝑏_𝑛(𝑡)^𝑁−1_𝑖=1𝐵_𝑛𝑖𝑢_𝑖(𝑡)−𝑑_𝑛(𝑡)−𝑒_𝑛(𝑡)+𝑞_𝑛(𝑡),𝑛=1,…,𝑁−1. (39)

The system of ordinary differential equations (38) in time may be solved using any standard technique to find 𝑢_𝑛(𝑡) and then 𝑢_𝑁(𝑥,𝑡) from (22).

4. Numerical Results

In order to check the accuracy and reliability of the proposed algorithm, we present two numerical examples using the proposed algorithm. In the first example, we compute the space fractional diffusion equation to check the accuracy, and space fractional advection-dispersion equation with variable coefficients is solved in the second example which confirms the good accuracy of our method. Comparing the results obtained by various choices of Jacobi parameters 𝛼 and 𝛽 and results presented elsewhere reveals that the present method is very effective and convenient for all choices of 𝛼 and 𝛽.

Example 1. Consider the space fractional diffusion equation (see, [42, 48, 49]): 𝜕𝑢(𝑥,𝑡)−Γ()𝑥𝜕𝑡2.2^2.86𝜕^1.8𝑢(𝑥,𝑡)𝜕𝑥^1.8=−(1+𝑥)𝑒^−𝑡𝑥³,𝑥∈(0,1),𝑡∈(0,𝑇], (40) with the initial condition: 𝑢(𝑥,0)=𝑥³for𝑥∈(0,1), (41) and the boundary conditions: 𝑢(0,𝑡)=0,𝑢(1,𝑡)=𝑒^−𝑡,𝑡∈(0,𝑇]. (42) The exact solution to this problem is 𝑢(𝑥,𝑡)=𝑒^−𝑡𝑥³. (43)

In Table 1, we list the maximum absolute errors using J-GL-C method with three choices of the Jacobi parameters 𝛼 and 𝛽 and various choices of the 𝑁.

Table 1 

Comparing maximum absolute errors of the proposed method and [42, 49].

We contrast our results with the corresponding results for the backward Euler finite difference scheme (BEFD [49]), the fractional Crank-Nicholson approach (CN [50]), and the extrapolated fractional Crank-Nicholson approach (Extra CN [50]) which we have presented in the fifth, sixth, and seventh columns of Table 1. We should note that for all values of 𝛼 and 𝛽, the proposed method is always more accurate than the results of CN [50], Extra CN [50], and BEFD [49], which shows the spectral accuracy of our method.

In Figures 1 and 2, the analytical solutions and the numerical solutions for 𝑥∈(0,1), 𝑡=0.1,0.5,0.7,0.9, and 𝑥=0.2,0.6,0.8,1.0, 𝑡∈(0,1), are shown, respectively. Consequently, we see that all numerical solutions are in complete agreement with the analytical solutions.

Figure 1 

The comparison of the curves of analytical solutions and approximate solutions at 𝑁=24 and 𝑡=0.1,0.5,0.7,0.9.

Figure 2 

The comparison of the curves of analytical solutions and approximate solutions at 𝑁=24 and 𝑥=0.2,0.6,0.8,1.0.

Example 2. Consider the space fractional advection-dispersion equation with variable coefficients: 𝜕𝑢(𝑥,𝑡)𝜕𝑡=−𝑡𝑥𝜃𝜕𝑢(𝑥,𝑡)+𝑡𝜕𝑥²𝑥^𝜃Γ()𝜕1+𝜃^𝜃𝑢(𝑥,𝑡)𝜕𝑥^𝜃(),[](],+𝑓𝑥,𝑡𝑥∈0,2,𝑡∈0,𝑇 (44) where 𝑓(𝑥,𝑡)=𝑥^𝜃4𝜋cos(4𝜋𝑡)+𝑡sin(4𝜋𝑡)−𝑡²sin(4𝜋𝑡), (45) with the initial condition: 𝑢(𝑥,0)=0for𝑥∈[0,2], (46) and the boundary conditions: 𝑢(0,𝑡)=0,𝑢(2,𝑡)=2^𝜃sin(4𝜋𝑡),𝑡∈(0,𝑇]. (47) The exact solution to this problem is 𝑢(𝑥,𝑡)=sin(4𝜋𝑡)𝑥^𝜃. (48)

For the sake of comparison of some different values of Jacobi parameters 𝛼 and 𝛽, we introduce in Table 2 the maximum absolute errors between the exact and numerical solutions using Jacobi Gauss-Lobatto collocation method with 𝜃=1.9, 𝑡∈(0,2], and various choices of 𝑁. Consequently, we conclude that all numerical solutions are in good agreement with the analytical solutions in all choices of 𝛼 and 𝛽.

Table 2 

Comparing maximum absolute errors for different choices of 𝑎𝑙 and 𝛽 and 𝑁=3,6,12,24.

In case of Chebyshev polynomials of the first kind 𝛼=𝛽=−1/2, the space-time graphs of approximate solutions at 𝑁=28 for the two choices 𝜃=1.01 and 𝜃=1.81 are shown in Figures 3 and 4, respectively. From these figures, it can be seen that the numerical solutions are in excellent agreement with the exact solutions. Numerical simulation is given in Figure 5 to compare the curves of exact solution and approximate solution (in case of Legendre polynomials 𝛼=𝛽=0 and 𝑁=16) for 𝜃=1.5 with 𝑥∈(0,2) and 𝑡=0,0.4,0.6,1.2,1.8. Moreover, the curves of exact solution and approximate solution (in case of Chebyshev polynomials of the second kind 𝛼=𝛽=1/2 and 𝑁=16) for 𝜃=1.5 with 𝑡∈(0,2) and 𝑥=0.2,0.6,1.2,1.6,2 are sketched in Figure 6. Consequently, we see that the curves of the exact and approximate solutions almost coincide for all chosen values of 𝑡 and 𝑥.

Figure 3 

The space-time graph of approximate solutions at 𝜃=1.01.

Figure 4 

The space-time graph of approximate solutions at 𝜃=1.81.

Figure 5 

The comparison of the curves of analytical solutions and approximate solutions at 𝑡=0,0.4,0.6,1.2,1.8.

Figure 6 

The comparison of the curves of analytical solutions and approximate solutions at 𝑥=0.2,0.6,1.2,1.6,2.

The obtained results of this example show that the Jacobi Gauss-Lobatto collocation method is simple and very accurate for all values of 𝛼 and 𝛽. Also by selecting limited Gauss-Lobatto collocation points, excellent numerical results are obtained.

5. Conclusion and Future Work

In this paper, we have proposed the Jacobi Gauss-Lobatto collocation spectral approximation for tackling fractional-in-space advection-dispersion equation subject to initial-boundary conditions. Applying the collocation method has reduced the problem to system of ordinary differential equations in time. This system may be solved by an implicit iterative technique. One of the main advantages of the proposed method is the Legendre Gauss-Lobatto collocation approximation, and the four kinds of Chebyshev Gauss-Lobatto collocation approximations may be obtained as special cases of the proposed Jacobi Gauss-Lobatto collocation approximation by taking the corresponding special cases of the Jacobi parameters 𝛼 and 𝛽. The numerical results given in Section 4 demonstrate the good accuracy of proposed algorithm.

The implementation of Jacobi Gauss-Lobatto collocation spectral approximation for time-space fractional advection-dispersion equations may also constitute another line of our future lines of research. We also conclude that this algorithm can be useful in dealing with coupled nonlinear partial differential equations.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, Saudi Arabia under Grant no. (130-141-D1433). The author, therefore, acknowledges the DSR technical and financial support.