Abstract

A new spectral shifted Legendre Gauss-Lobatto collocation (SL-GL-C) method is developed and analyzed to solve a class of two-dimensional initial-boundary fractional diffusion equations with variable coefficients. The method depends basically on the fact that an expansion in a series of shifted Legendre polynomials , for the function and its space-fractional derivatives occurring in the partial fractional differential equation (PFDE), is assumed; the expansion coefficients are then determined by reducing the PFDE with its boundary and initial conditions to a system of ordinary differential equations (SODEs) for these coefficients. This system may be solved numerically by using the fourth-order implicit Runge-Kutta (IRK) method. This method, in contrast to common finite-difference and finite-element methods, has the exponential rate of convergence for the two spatial discretizations. Numerical examples are presented in the form of tables and graphs to make comparisons with the results obtained by other methods and with the exact solutions more easier.

1. Introduction

In recent years there has been a high level of interest of employing spectral methods for numerically solving many types of integral and differential equations, due to their ease of applying them for finite and infinite domains [1–8]. The speed of convergence is one of the great advantages of spectral method. Besides, spectral methods have exponential rates of convergence; they also have high level of accuracy. From the overview of spectral approximation to differential equations, the spectral methods have been divided to four types, namely, collocation [9–11], tau [12, 13], Galerkin [14, 15], and Petrov Galerkin [16, 17] methods. The main idea of all versions of spectral methods is to express the spectral solution of the problem as a finite sum of certain basis functions (combination of orthogonal polynomials or functions) and then to choose the coefficients in order to minimize the difference between the exact and numerical solutions as possible. The spectral collocation method is a specific type of spectral methods, which is more applicable and widely used to solve almost types of differential equations [18–21].

Fractional differential equations [22–31] model many phenomena in several fields such as biology, viscoelasticity, finance, fluid mechanics, physics, chemistry, and engineering. Most fractional differential equations do not have exact solutions, so approximation and numerical techniques must be used. Here, we focus on the application of SL-GL-C scheme for the numerical solution of the time dependent fractional diffusion equation in two dimensions. This problem models phenomena exhibiting [32–34] anomalous diffusion that cannot be modeled accurately by the classical two-dimensional diffusion equations with variable coefficients. For the previous numerical methods used for time dependent fractional diffusion equation, see, for example, the papers in [35–40].

The solution of this equation is approximated as which can be expressed as a finite expansion of shifted Legendre polynomials for the space variables and then evaluate the spatial fractional partial derivatives of the finite expansion of at the shifted Legendre-Gauss-Lobatto (SL-GL) quadrature points. Substituting these approximations in the underlined PFDE provides a SODEs in time. This system may be solved by IRK scheme. This scheme is one of the suitable schemes for solving SODEs of first order. Since the IRK scheme has excellent stability properties, then it allows us to reduce computational costs. Indeed, this is the first work for employing SL-GL-C scheme to solve two-dimensional fractional diffusion equations.

This paper is organized as follows. We present few relevant properties of fractional integration in the Riemann-Liouville and shifted Legendre polynomials in the coming section. Section 3 presents the numerical solution of the two-dimensional fractional diffusion equations with initial-boundary conditions. Section 4 is devoted to solve two test examples. Finally, some concluding remarks are given in the last section.

2. Preliminaries and Notation

2.1. The Fractional Integration in the Riemann-Liouville Sense

There are several definitions of the fractional integration of order and not necessarily equivalent to each other; see [41, 42]. The most used definition is due to Riemann-Liouville, which is defined as One of the basic properties of the operator is The Riemann-Liouville fractional derivative of order will be denoted by . The next relation defines the Riemann-Liouville fractional derivative of order as where , and is the smallest integer greater than .

Lemma 1. If , then

2.2. Properties of Shifted Legendre Polynomials

The well-known Legendre polynomials are defined on the interval . Some properties about the standard Legendre polynomials will be recalled in this section. The Legendre polynomials satisfy the following Rodrigue's formula: We recall also that is a polynomial of degree , and therefore (the th derivative of ) is where The Legendre polynomials satisfy the following relations: and the orthogonality relation where , , and is the Kronecker delta function. The Legendre-Gauss-Lobatto quadrature is used to evaluate the previous integral accurately. For any , we get We introduce the following discrete inner product: where () and () are used as usual nodes and the corresponding Christoffel numbers in the interval , respectively. Let the shifted Legendre polynomials be denoted by ; . Then can be obtained with the aid of the following recurrence formula: The analytic form of the shifted Legendre polynomials of degree is given by The Riemann-Liouville fractional integration of shifted Legendre polynomials of degree is given by where . The orthogonality condition is where   and  .

The function , square integrable in , may be expressed in terms of shifted Legendre polynomials as where the coefficients are given by In practice, only the first -terms shifted Legendre polynomials are considered. Hence can be expressed in the form

3. Shifted Legendre Spectral Collocation Method

In this section, we derive the SL-GL-C method for solving two-dimensional fractional diffusion equation of the form which is subject to the initial condition and boundary conditions

Denote by , and , the nodes and Christoffel numbers of the standard (shifted) Legendre-Gauss-Lobatto quadrature on the intervals , respectively. Then one can clearly deduce that and if denotes the set of all polynomials of degree at most , then it follows for any that

Here, the SL-GL-C method will be employed to transform the previous two-dimensional fractional diffusion equation into SODEs. Let us expand the dependent variables in the form We observe from (15) that the coefficients can be approximated by The expression (24) can be rewritten as where . In what follows, the approximation of the fractional spatial partial derivative with respect to for can be computed as

Also, the fractional spatial partial derivative of fractional order with respect to of the approximate solution is The previous fractional derivatives can be computed at the SL-GL interpolation nodes as where In the proposed SL-GL-C method the residual of (19) is set to zero at of SL-GL points. Moreover, the values of , , , and are given by Therefore, adapting (26)–(31) enables one to write (19)–(21) as a SODEs: which is subject to the initial conditions: where Finally, we can arrange (32) and (33) to their matrix formulation, which can be solved by IRK scheme: where

The Runge-Kutta method can be expressed as one of powerful numerical integrations tools used for initial value SODEs of first order. The IRK method presents a subclass of the well-known family of Runge-Kutta methods and has many applications in the numerical solution of systems of ordinary differential equations.

4. Numerical Results

In this section, we present two numerical examples to show the accuracy, robustness, and applicability of the proposed method. We compare the results obtained by our method with those obtained in [43].

The difference between the measured value of the approximate solution and the exact solution (absolute error) is given by where , , , , and are the space vectors, time, exact, and numerical solutions, respectively. Moreover, the maximum absolute error (MAEs) is given by Also, we define the norm infinity as