Chebyshev Pseudo-Spectral Method for Solving Fractional Advection-Dispersion Equation

Abstract

Fractional differential equations have recently been applied in various areas of engineering, science, finance, applied mathematics, bio-engineering and others. However, many researchers remain unaware of this field. In this paper, an efficient numerical method for solving the fractional Advection-dispersion equation (ADE) is considered. The fractional derivative is described in the Caputo sense. The method is based on Chebyshev approximations. The properties of Chebyshev polynomials are used to reduce ADE to a system of ordinary differential equations, which are solved using the finite difference method (FDM). Moreover, the convergence analysis and an upper bound of the error for the derived formula are given. Numerical solutions of ADE are presented and the results are compared with the exact solution.

Share and Cite:

Sweilam, N. , Khader, M. and Adel, M. (2014) Chebyshev Pseudo-Spectral Method for Solving Fractional Advection-Dispersion Equation. Applied Mathematics, 5, 3240-3248. doi: 10.4236/am.2014.519301.

1. Introduction

Ordinary and partial fractional differential equations (FDEs) have been the focus of many studies due to their frequent appearance in various applications in fluid mechanics, viscoelasticity, biology, physics and engineering [1] [2] . Consequently, considerable attention has been given to the solutions of FDEs of physical interest. Most FDEs do not have exact solutions, so approximate and numerical techniques [3] [4] , must be used. Recently, several numerical methods to solve FDEs have been given such as variational iteration method [5] , homotopy perturbation method [3] [6] , Adomian decomposition method [7] [8] , homotopy analysis method [9] , collocation method [10] [11] and finite difference method [12] - [17] .

We introduce some necessary definitions and mathematical preliminaries of the fractional calculus theory that will be required in the present paper.

Definition 1

The Caputo fractional derivative operator of order is defined in the following form

where, , is the gamma function.

Similar to integer-order differentiation, Caputo fractional derivative operator is a linear operation

where and are constants. For the Caputo’s derivative we have [2]

(1)

(2)

We use the ceiling function to denote the smallest integer greater than or equal to and. Recall that for, the Caputo differential operator coincides with the usual differential operator of integer order.

For more details on fractional derivatives definitions and its properties see [2] .

Anomalous, or non-Fickian, dispersion has been an active area of research in the physics community since the introduction of continuous time random walks (CTRW) by Montroll and Weiss [1965]. These random walks extended the predictive capability of models built on the stochastic process of Brownian motion, which is the basis for the classical advectiondispersion equation ( ADE ).

A fractional ADE . is a generalization of the classical ADE in which the second-order derivative is replaced with a fractional-order derivative. In contrast to the classical ADE , the fractional ADE has solutions that resemble the highly skewed and heavy-tailed breakthrough curves observed in field and laboratory studies.

When a fractional Fick’s law replaces the classical Fick’s law in an Eulerian evaluation of solute transport in a porous medium, the result is a fractional ADE.

It describes the spread of solute mass over large distances via a convolutional fractional derivative.

We consider the initial-boundary value problem of the fractional Advection-dispersion equation which is usually written in the following form

(3)

where, , is the source term, is a constant and is the Caputo fractional derivative with respect to and of order, where.

Under the zero boundary conditions

(4)

and the following initial condition

(5)

In the last few years appeared many papers to study this model (3)-(5) [1] [18] - [22] , the most of these papers study the ordinary case of such problem but in this paper we study the fractional case.

Our idea is to apply the Chebyshev collocation method to discretize (3) to get a linear system of ODEs thus greatly simplifying the problem, and use FDM [12] to solve the resulting system.

The organization of this paper is as follows. In the next section, we obtain the approximation of fractional derivative. In Section 3, we prove the error analysis of the proposed formula. In Section 4, we implement chebyshev collocation method to the solution of (3). As a result a system of ordinary differential equations is formed and the solution of the considered problem is introduced. In Section 5, we give some numerical results to clarify the proposed method. Also a conclusion is given in Section 6.

2. Derivation of the Approximate Formula

The well known Chebyshev polynomials are defined on the interval and can be determined with the aid of the following recurrence formula [23]

.

The analytic form of the Chebyshev polynomials of degree is given by

, (6)

where denotes the integer part of. The orthogonality condition is

In order to use these polynomials on the interval we define the so called shifted Chebyshev polynomials by introducing the change of variable. So, the shifted Chebyshev polynomials are defined as.

The analytic form of the shifted Chebyshev polynomials of degree is given by

. (7)

The function, which belongs to the space of square integrable in, may be expressed in terms of shifted Chebyshev polynomials as

, (8)

where the coefficients are given by

. (9)

In practice, only the first -terms of shifted Chebyshev polynomials are considered. Then we have

(10)

Khader [24] introduced a new approximate formula of the fractional derivative and used it to solve numerically the fractional diffusion equation.

The main approximate formula of the fractional derivative of is given in the following theorem.

Theorem 1 [24]

Let be approximated by Chebyshev polynomials as in (10) and also suppose then

(11)

where is given by

(12)

Also, in this section, special attention is given to study the convergence analysis and evaluate an upper bound of the error of the proposed approximate formula.

Theorem 2 (Chebyshev truncation theorem) [23]

The error in approximating by the sum of its first terms is bounded by the sum of the absolute values of all the neglected coefficients. If

, (13)

then

, (14)

for all, all, and all.

Theorem 3 [25]

The Caputo fractional derivative of order for the shifted Chebyshev polynomials can be expressed in terms of the shifted Chebyshev polynomials themselves in the following form

, (15)

where

.

Theorem 4 [11]

The error in approximating by is bounded by

. (16)

3. Procedure Solution of the Fractional Advection-Dispersion Equation

Consider the fractional Advection-dispersion equation of type given in Equation (3). In order to use Chebyshev collocation method, we first approximate as

. (17)

From Equations (3), (17) and Theorem 1 we have

. (18)

We now collocate Equation (18) at points, as

. (19)

For suitable collocation points we use roots of shifted chebyshev polynomial.

Also, by substituting Equations (17) in the boundary conditions (4) we can obtain equations as follows

. (20)

Equation (19), together with equations of the boundary conditions (20), give of ordinary differential equations which can be solved, for the unknowns, , using the finite difference method, as described in the following section.

4. Numerical Results

In this section, we present a numerical example to illustrate the efficiency and the validation of the proposed numerical method when applied to solve numerically the fractional Advection-dispersion equation. Consider the ADE (3) with and the following source term

, (21)

and the boundary conditions, with the initial condition.

The exact solution of Equation (3) in this case is

. (22)

We apply the proposed method with, and approximate the solution as follows

. (23)

Using Equation (19) we have

, (24)

where are roots of shifted Chebyshev polynomial, i.e.

By using Equations (20) and (24) we can obtain the following system of ODEs

, (25)

, (26)

, (27)

, (28)

where

Now, to use FDM for solving the system (25)-(28), we will use the following notations: to be the integration time for. Define,. Then the system (25)-(28), is discretized in time and takes the following form

, (29)

, (30)

, (31)

. (32)

We can write the above system (29)-(32) in the following matrix form as follows

. (33)

We use the notation for the above system

, (34)

where and

The obtained numerical results by means of the proposed method are shown in Table 1 and Figures 1-4. In Table 1, the absolute error between the exact solution and the approximate solution at and with the final time are given. Also, in Figure 1 and Figure 2, comparison between the exact solution and the approximate solution at with time step, and, are presented, respectively. Also, in Figure 3 and Figure 4, the behavior of the approximate solution at and with different values of and are presented, respectively. From, these figures, we can see that the behavior of the approximate solution depends on the order of the fractional derivative.

5. Conclusion and Remarks

The properties of the Chebyshev polynomials are used to reduce the fractional Advection-dispersion equation to the solution of system of ODEs which solved by using FDM. The fractional derivative is considered in the

Table 1. The absolute error between the exact solution and the approximate solution at, and.

Figure 1. Comparison between the exact solution and the approximate solution at T = 0.5 with = 0.0025; m = 3.

Figure 2. Comparison between the exact solution and the approximate solution at T = 0.5 with = 0.0025; m = 5.

Figure 3. The behavior of the approximate solution at different values of α at β = 0.8.

Figure 4. The behavior of the approximate solution at different values of β at α = 1.8.

Caputo sense. In this article, special attention is given to studying the convergence analysis and estimating an upper bound of the error for the proposed approximate formula of the fractional derivative. The solution obtained using the suggested method is in excellent agreement with the already existing ones and shows that this approach can be solved the problem effectively. From the resulted numerical solution, we can conclude that the used techniques in this work can be applied to many other problems. It is evident that the overall errors can be made smaller by adding new terms from the series (23). Comparisons are made between the approximate solution and the exact solution to illustrate the validity and the great potential of the technique. All computations in this paper are done using Matlab 8.

Conflicts of Interest

The authors declare no conflicts of interest.