A highly accurate collocation algorithm for 1 + 1 and 2 + 1 fractional percolation equations

Abstract

This paper reports two successive spectral collocation methods, that enable easy and highly accuracy discretization, for 1 + 1 and 2 + 1 fractional percolation equations (FPEs). The first step depends mainly on the shifted Legendre Gauss–Lobatto collocation method for spatial discretization. An expansion in a series of shifted Legendre polynomials for the approximate solution and its spatial derivatives occurring in the FPE is investigated. In addition, the Legendre-Gauss–Lobatto quadrature rule is established to treat the boundary conditions. Thereby, the expansion coefficients are then determined by reducing the FPE with its boundary conditions to a system of ordinary differential equations for these coefficients. The second step is to propose the shifted Chebyshev Gauss–Radau collocation scheme, for temporal discretization, to reduce such a system to a system of algebraic equations, which is far easier to solve. The proposed collocation scheme, both in temporal and spatial discretizations, is successfully extended to the numerical solution of two-dimensional FPEs. An upper bound of the absolute error is obtained of the approximate solution for the two-dimensional case. Convergence properties of the method are discussed through numerical examples. Several numerical examples with comparisons are reported to highlight the high accuracy of the present method over other numerical techniques.

Keywords

Space fractional percolation equations two-dimensional fractional percolation equations collocation method orthogonal polynomials Gauss-type quadrature

1. Introduction

Spectral methods (Heinrichs, 1989; Canuto et al., 2006; Doha et al., 2014b; Eslahchi et al., 2014; Kayedi-Bardeh et al., 2014; Khalil and Khan, 2014; Abdelkawy et al., 2015a) are often efficient and highly accurate schemes when compared with local methods. The speed of convergence is one of the great advantages of spectral approximations. The spectral collocation method is a specific type of spectral method that is more applicable and widely used to solve almost all types of differential equations (Bhrawy, 2014; Doha et al., 2014a; Bhrawy and Zaky, 2015b). In pseudo-spectral techniques for partial differential equations, the boundary conditions play a much more crucial role than for ordinary differential equations, and it becomes important to decide whether to satisfy the related conditions implicitly, in the formulation of the collocation scheme, or to investigate the related conditions as additional constraints. The boundary conditions were treated implicitly in some recent pseudo-spectral approximations (Bhrawy, 2014; Doha et al., 2014a).

The development and analysis of fractional calculus began in recent decades, when the fractional differential equation emerged as a tool for the description of phenomena in nature. Fractional differential equations (Giona and Roman, 1992; Kirchner et al., 2000; Magin, 2006; Li and Deng, 2007; Garrappa and Popolizio, 2011; Alipour et al., 2012; Baleanu et al., 2012; Machado et al., 2013; Rostamy et al., 2013) are used to model many phenomena in several fields (Pooseh et al., 2013; Dehghan et al., 2014; Lazo and Torres, 2014; Pinto and Carvalho, 2015; Raja and Chaudhary, 2015; Xu et al., 2015). Numerical techniques are widely used by scientists and engineers to solve fractional PDEs. A major advantage of numerical techniques is that a numerical solution can be obtained even when a problem has no analytical solution. In some cases, PDEs of fractional order can be solved analytically, where finding their solutions in the general case is hard and limited to the linear one. Therefore, there has been a growing interest in recent decades in developing numerical methods for solving FDEs. Several numerical treatments based on the improvement of finite difference schemes for FDEs are given in Meerschaert and Tadjeran (2006), Ding et al. (2010) and Wang and Du (2014). Also, efficient spectral algorithms (He, 1998; Chou et al., 2006; Chen et al., 2011b; Doha et al., 2011; Doha et al., 2012; Chen et al., 2013; Bhrawy, 2014; Bhrawy et al., 2015a; Chen et al., 2014; Irandoust-Pakchin et al., 2014; Abdelkawy et al., 2015b; Ezz-Eldien et al., 2015) have been designed and developed for solving different kinds of fractional differential equation.

Percolation flow problems (Chou et al., 2006) have been discussed in many fields including groundwater dynamics, seepage hydraulics, and fluid dynamics in porous media. He (1998) introduced analytical solutions for fractional percolation equations (FPEs) by using the variational iteration method. Fractional derivatives are becoming widely used and accepted in models of percolation flow problems. Recently, three advanced implicit finite difference methods have been investigated in Chen et al. (2011b, 2013, 2014) to discretize the numerical solutions of FPEs in one- and two-dimensional spaces. Along the same line of thought, Chen et al. (2011b) proposed a new implicit finite difference scheme for solving the 1 + 1 FPE. A developed numerical algorithm has been achieved and analyzed by Chen et al. (2013), for solving the 2 + 1 FPE. Several authors have also improved and developed efficient numerical methods for approximating the solution of other similar fractional PDEs (Chen et al., 2011a; Al-Khaled and Alquran, 2014; Irandoust-Pakchin et al., 2014; Mohebbi et al., 2014; Bhrawy and Zaky, 2015a; Bhrawy et al., 2015a, 2015b).

The main aim of the present paper is to propose a numerical method that improves the accuracy of the numerical solutions of FPEs in one- and two-dimensional space. The main advantage of the present method is that it proposes a collocation scheme for both temporal and spatial discretizations. Firstly, the shifted Legendre Gauss–Lobatto collocation (SL-GL-C) is proposed, with a suitable modification for treating the boundary conditions, for spatial discretization. This treatment, for the conditions, improves the accuracy of the scheme greatly. Therefore, the FPE with its boundary conditions is reduced to a system of ordinary differential equations (SODEs) subject to a vector of initial values. Secondly, the SC-GR-C is then investigated for temporal discretization, which is more reasonable for solving initial value problems. Thereby, the problem is reduced to a system of algebraic equations which makes it far easier to solve. In addition, this algorithm is developed for the numerical solution of FPEs in two dimensions. An upper bound of the absolute error is obtained for the approximate solution for the two-dimensional case. Finally, several numerical examples with comparisons highlighting the high accuracy and effectiveness of the present method are included. By choosing relatively limited Legendre Gauss–Lobatto and Chebyshev Gauss–Radau collocation nodes, we are able to obtain highly accurate solutions, confirming the applicability and high accuracy of the proposed method over other numerical schemes in the literature.

This article is organized as follows. Section 2 presents some fractional calculus preliminaries, and shifted Legendre and shifted Chebyshev polynomials. Spectral approximation schemes of 1 + 1 and 2 + 1 FPEs based on a combination of SL-GL-C and SC-GR-C methods are presented in Sections 3 and 4, respectively. Numerical results are introduced in Section 5. Finally, we end the paper with some concluding remarks.

2. Preliminaries and notation

This section presents several useful fractional definitions, and shifted Legendre and shifted Chebyshev polynomials.

2.1. Fractional calculus

There are many definitions of the fractional-order derivative, which are not necessarily equivalent to each other (Miller and Ross, 1993; Kayedi-Bardeh et al., 2014; Ezz-Eldien et al., 2015; Jafari and Tajadodi, 2015). Riemann–Liouville and Caputo fractional definitions are the two most popular definitions from all the other definitions of fractional calculus which have been introduced recently.

Definition 2.1

The Riemann–Liouville fractional integral of order ν ≥ 0 is obtained from

J ν f (x) = 1 Γ (ν) \int_{0}^{x} (x - ζ) ν - 1 f (ζ) d ζ, ν > 0, x > 0, J 0 f (x) = f (x)

(1)

where

Γ (ν) = \int_{0}^{\infty} x ν - 1 e - x d x

(2)

is the gamma function.

The fractional operator J^ν satisfies

J ν J μ f (x) = J ν + μ f (x), J ν J μ f (x) = J μ J ν f (x), J ν x β = Γ (β + 1) Γ (β + 1 + ν) x β + ν

(3)

The Riemann–Liouville fractional derivative of order μ satisfies

D μ f (x) = d n x n (J n - μ f (x))

(4)

where n − 1 < μ ≤ n, n ∈ N and n is the smallest integer greater than μ.

Lemma 2.1

If n − 1 < μ ≤ n, n ∈ N, then

D μ J μ f (x) = f (x), J μ D μ f (x) = f (x) - \sum j = 0 n - 1 f (j) (0 +) x j j!, x > 0

(5)

2.2. Shifted Legendre Gauss–Lobatto interpolation

In this subsection, we recall some approximation results for the shifted Legendre Gauss–Lobatto (SL-GL) interpolation, which play important roles in the proposed collocation scheme. The Legendre polynomials P_k(x) (k = 0,1…) satisfy the following Rodrigue's formula

P k (x) = (- 1) k 2 k k! D k ((1 - x 2) k)

(6)

also we recall that P_k(x) is a polynomial of degree k. Accordingly,

P_{k}^{(p)} (x)

(the pth derivative of P_k(x)) is given by

P_{k}^{(p)} (x) = \sum_{i = 0 (i + k = even)}^{k - p} C p (k, i) P i (x)

(7)

where

C p (k, i) = 2 p - 1 (2 i + 1) Γ (p + k - i 2) Γ (p + k + i + 12) Γ (p) Γ (2 - p + k - i 2) Γ (3 - p + k + i 2)

Next, denoting by

‖ u ‖

and (u,v) the norm and inner product of space L²[ − 1,1]. The set of P_k(x) is a complete orthogonal system in L²[ − 1,1], namely

(P j (x), P k (x)) = \int_{- 1}^{1} P j (x) P k (x) d x = h k δ jk

(8)

where

h i = 22 i + 1

and δ_jk is the Dirac function. Thus for any v ∈ L²[ − 1,1]

v (x) = \sum_{i = 0}^{\infty} a i P i (x), a i = 1 h i \int_{- 1}^{1} v (x) P i (x) d x

(9)

For any positive integer N, let S_N[ − 1,1] be the set of all polynomials of degree at most N, due to the L-GL quadrature. Thus, for any φ ∈ S_2N−1[ − 1,1] we obtain

\int_{- 1}^{1} φ (x) d x = \sum_{i = 0}^{N} ϖ N, i φ (x N, i)

(10)

where x_N,i (0 ≤ i ≤ N) and ϖ_N,i (0 ≤ i ≤ N) are the nodes and Christoffel numbers of the L-GL interpolation on the interval [ − 1,1], respectively. Now, it is quite useful in the sequel to define the following norm and discrete inner product

‖ u ‖ N = (u, v) N 12, (u, v) N = \sum_{j = 0}^{N} u (x N, j) v (x N, j) ϖ N, j

(11)

Let us denote by P_L,i(x) the shifted Legendre polynomials which are defined on the interval [0,L]. These polynomials can be generated from the following recurrence relation

(j + 1) P L, j + 1 (x) = (2 j + 1) \times (2 x L - 1) P L, j (x) - jP L, j - 1 (x), j = 1, 2, \dots

(12)

The analytical form of P_L,i(x) can be written as

P L, i (x) = \sum_{k = 0}^{i} (- 1) i + k (i + k)! (i - k)! (k!) 2 L k x k

(13)

The Riemann–Liouville fractional integration of P_L,i(x) may be obtained from

J ν P L, i (x) = \sum_{k = 0}^{i} (- 1) k + i (k + i)! (- k + i)! (k!) 2 t k J ν x k = \sum_{k = 0}^{i} (- 1) k + i (k + i)! k! (- k + i)! (k!) 2 t k Γ (k + ν + 1) \times x k + ν, i = 0, 1, \dots, N

(14)

where P_L,i(0) = ( − 1)ⁱ.

The orthogonality condition is

\int_{0}^{L} P L, j (x) P L, k (x) w L (x) d x = h_{k}^{L} δ jk

(15)

where

w L (x) = 1 and h_{k}^{L} = L 2 k + 1 .

If a function u(t) ∈ L²[0,L], then one can express it by means of P_L,i(t) as

u (t) = \sum_{i = 0}^{\infty} c i P L, i (t)

where c_i is given by

c i = 1 h_{i}^{L} \int_{0}^{L} u (t) P L, i (t) x, i = 0, 1, 2, \dots

(16)

In the approximation, u(x) may be expanded as

u N (t) \approx \sum_{i = 0}^{N} c i P L, i (t)

(17)

2.3. Shifted Chebyshev Gauss–Radau interpolation

The Chebyshev polynomials are defined on the interval [ − 1,1], by

T k (t) = \cos (k \arccos (t)), k \geq 0

(18)

Also

T k (\pm 1) = (\pm 1) k, T k (- t) = (- 1) k T k (t)

(19)

Let

w c (t) = 1 \sqrt{1 - t 2}

, then we introduce the following norm and inner product of the the weighted space

L_{w c}^{2} [- 1, 1]

‖ u ‖ w c = (u, u) 12 w c, (u, v) w c = \int_{- 1}^{1} u (t) v (t) w c (t) d t

(20)

The set of Chebyshev polynomials satisfies

‖ T k ‖ 2 w c = h_{k}^{c} = {ς k 2 π, k = j, 0, k \neq j, ς 0 = 2, ς k = 1, k \geq 1

(21)

Now, we introduce the following norm and discrete inner product

‖ u ‖ w c = (u, u) 12 w c, (u, v) w c = \sum_{j = 0}^{N} u (t N, j) v (t N, j) ϖ_{N, j}^{c}

(22)

Let us denote by T_T,n(t) the shifted Chebyshev polynomials which are defined on the interval [0,T]. The analytical form of T_T,n(t) is obtained from

T T, n (t) = n \sum_{k = 0}^{n} (- 1) n - k (n + k - 1)! 2 2 k (n - k)! (2 k)! T k t k

(23)

where T_T,n(0) = ( − 1)ⁿ and T_T,n(T) = 1.

The orthogonality condition is

\int_{0}^{L} T T, m (t) T T, n (x) w T (x) d x = δ_{mn}^{c} h_{n}^{T}

(24)

where

w T (t) = 1 \sqrt{T t - t 2} and c h_{n}^{T} = c n 2 π, with c 0 = 2, c i = 1, and i \geq 1 .

As in the previous subsection, if $u (t) \in L_{w T (t)}^{2} [0, T]$ , then one can express it by means of T_T,i(t) as

u (t) = \sum_{j = 0}^{\infty} a j T T, j (t)

(25)

where

a j = 1 c h_{j}^{T} \int_{0}^{T} u (t) T T, j (t) w T (t) d t, j = 0, 1, 2, \dots

(26)

3. One-dimensional space fractional percolation equation

This section outlines the key aspects of discretization for the numerical solution of the one-dimensional FPE using a new spectral technique.

Consider the following 1 + 1 FPE

\partial t u (x, t) = \partial ν 1 x (g 1 (x) \partial ν 2 x u (x, t)) + f (x, t), \times (x, t) \in [0, L] \times [0, T]

(27)

subject to

u (x, 0) = g 2 (x), x \in [0, L]

(28)

and

u (0, t) = g 3 (t), u (L, t) = g 4 (t), t \in [0, T]

(29)

where ν₁ and ν₂ (0 < ν₁, ν₂ < 1) are the fractional orders, and f(x,t), g₁(x), g₂(x), g₃(t) and g₄(t) are given functions. We discretize equation (27) for the spatial variable x by the SL-GL-C method. While the SC-GR-C method is applied to discretize the resulting SODEs in time advance.

3.1. SL-GL-C scheme for the space variable

We first propose the SL-GL-C method to transform the 1 + 1 FPE into SODEs. To this end, we approximate the spatial variable using the SL-GL-C method at some nodes. The nodes are the set of points in a specified domain. To acquire high accuracy, the choice of nodes is usually related to some Gaussian integration formula, see Canuto et al. (2006) for more details. The collocation points are taken to be the SL-GL quadrature nodes which we denote by $x_{N, i}^{L}$ .

Equation (27) may be restated as

\partial t u (x, t) = g 5 (x) \partial ν 2 x u (x, t) + g 1 (x) \partial ν 1 x (\partial ν 2 x u (x, t)) + f (x, t), 0 < ν 1, ν 2 < 1, (x, t) \in [0, L] \times [0, T]

(30)

where

g 5 (x) = \partial x ν 1 g 1 (x)

Now, we present the main steps of applying the SL-GL-C scheme to solve the FPE. Let the approximate solution of equation (27) be

u N (x, t) = \sum_{j = 0}^{N} a j (t) P L, j (x)

(31)

the unknown coefficients a_j(t) may be approximated by

a j (t) = 1 h_{j}^{L} \sum_{i = 0}^{N} P L, j (x_{N, i}^{L}) ϖ_{N, i}^{L} u (x_{N, i}^{L}, t)

(32)

Accordingly, the approximate solution (31) may be expanded to

u N (x, t) = \sum_{i = 0}^{N} (\sum_{j = 0}^{N} 1 h_{j}^{L} P L, j (x_{N, i}^{L}) P L, j (x) ϖ_{N, i}^{L}) u (x_{N, i}^{L}, t)

(33)

where

ϖ_{N, j}^{L}

(0 ≤ j ≤ N) are the weights of the SL-GL quadrature on [0,L].

The fractional derivative of the approximate solution u_N(x,t) is then estimated as

\partial ν 2 x u N (x, t) = \sum_{i = 0}^{N} (\sum_{j = 0}^{N} 1 h_{j}^{L} P L, j (x_{N, i}^{L}) \partial ν 2 x (P L, j (x)) ϖ_{N, i}^{L}) u (x_{N, i}^{L}, t),

(34)

For 0 < ν₂ < 1, the Riemann–Liouville fractional derivative is

\partial ν 2 x x k = 1 Γ (1 - ν 2) \partial x (\int_{0}^{x} χ k (x - χ) ν 2 d χ) = x k - ν 2 Γ (1 + k) Γ (1 + k - ν 2)

(35)

thus

\partial ν 2 x P L, i (x) = P_{L, i}^{(ν 2)} (x) = \sum_{k = 0}^{i} (- 1) i + k (k + i)! Γ (k + 1) (i - k)! (k!) 2 Γ (1 + k - ν 2) L k x k - ν 2

(36)

Consequently

\partial ν 2 x u N (x, t) = \sum_{i = 0}^{N} (\sum_{j = 0}^{N} 1 h_{j}^{L} P L, j (x_{N, i}^{L}) P_{L, i}^{(ν 2)} (x) ϖ_{N, i}^{L}) u (x_{N, i}^{L}, t)

(37)

Similarly, we obtain

\partial ν 1 x (\partial ν 2 x u N (x, t)) = \sum_{i = 0}^{N} (\sum_{j = 0}^{N} 1 h_{j}^{L} P L, j (x_{N, i}^{L}) \partial ν 1 x (P_{L, i}^{(ν 2)} (x)) ϖ_{N, i}^{L}) u (x_{N, i}^{L}, t) = \sum_{i = 0}^{N} (\sum_{j = 0}^{N} 1 h_{j}^{L} P L, j (x_{N, i}^{L}) (P_{L, i}^{(ν 1 + ν 2)} (x)) ϖ_{N, i}^{L}) u (x_{N, i}^{L}, t)

(38)

The fractional derivative of the numerical solution u_N(x,t) is then computed at the collocation points as

(\partial x ν 2 u N (x, t)) x = x_{N, n}^{L} = \sum_{i = 0}^{N} γ n, i u i (t), n = 1, 2, \dots, N - 1,

(39)

(\partial x ν 1 (\partial x ν 2 u N (x, t))) x = x_{N, n}^{L} = \sum_{i = 0}^{N} ɛ n, i u i (t), n = 1, 2, \dots, N - 1

(40)

where

u i (t) = u N (x_{N, i}^{L}, t),

γ n, i = \sum_{j = 0}^{N} ϖ_{N, i}^{L} h_{j}^{L} P L, j (x_{N, i}^{L}) (P_{L, j}^{(ν 2)} (x)) x = x_{N, n}^{L},

ɛ n, i = \sum_{j = 0}^{N} ϖ_{N, i}^{L} h_{j}^{L} P L, j (x_{N, i}^{L}) (P_{L, i}^{(ν 1 + ν 2)} (x)) x = x_{N, n}^{L}

Since the approximate solution and its derivatives in equations (33)–(40) do not satisfy the FPE (27) exactly, we enforce the approximate solution and its derivatives to satisfy it exactly at the chosen collocation points, by making the residuals equal to zero. In addition, the boundary conditions (29) will be satisfied at the two nodes 0 and L. Thus, we ensure that the boundary conditions are satisfied without adding additional equations. Therefore, adopting equations (30)–(40), enables one to write equation (27) in the form

. u n (t) = g 5 (x_{N, n}^{L}) (\sum_{i = 1}^{N - 1} γ n, i u i (t) + γ n, 0 u 0 (t) + γ n, N u N (t)) + f (x_{N, n}^{L}, t) + g 1 (x) (\sum_{i = 1}^{N - 1} ɛ n, i u i (t) + ɛ n, 0 u 0 (t) + ɛ n, N u N (t))

(41)

Obviously, the problem (27)–(29) is equivalent to the following SODEs

. u n (t) = g 5 (x_{N, n}^{L}) (\sum_{i = 1}^{N - 1} γ n, i u i (t) + γ n, 0 g 3 (t) + γ n, N g 4 (t)) + f (x_{N, n}^{L}, t) + g 1 (x) (\sum_{i = 1}^{N - 1} ɛ n, i u i (t) + ɛ n, 0 g 3 (t) + ɛ n, N g 4 (t)), n = 1, \dots, N - 1

(42)

subject to the initial values

u n (0) = g 2 (x_{N, n}^{L}), n = 1, \dots, N - 1

(43)

We can reformulate the previous equations in matrix form as

(. u 1 (t) . u 2 (t) \dots \dots \dots . u N - 2 (t) . u N - 1 (t)) = (F 1 (t, u 1 (t), u 2 (t), \dots, u N - 1 (t)) F 2 (t, u 1 (t), u 2 (t), \dots, u N - 1 (t)) \dots \dots \dots F N - 2 (t, u 1 (t), u 2 (t), \dots, u N - 1 (t)) F N - 1 (t, u 1 (t), u 2 (t), \dots, u N - 1 (t)))

(44)

subject to a vector of initial values

(u 1 (0) u 2 (0) \dots \dots \dots u N - 2 (0) u N - 1 (0)) = (g 2 (x_{N, 1}^{L}) g 2 (x_{N, 2}^{L}) \dots \dots \dots g 2 (x_{N, N - 2}^{L}) g 2 (x_{N, N - 1}^{L}))

(45)

where

F n (t, u 1 (t), u 2 (t), \dots, u N - 1 (t)) = g 5 (x_{N, n}^{L}) (\sum_{i = 1}^{N - 1} γ n, i u i (t) + γ n, 0 g 3 (t) + γ n, N g 4 (t)) + f (x_{N, n}^{L}, t) + g 1 (x) (\sum_{i = 1}^{N - 1} ɛ n, i u i (t) + ɛ n, 0 g 3 (t) + ɛ n, N g 4 (t))

(46)

3.2. SC-GR-C scheme for the time variable

The second step of our algorithm is to propose an efficient collocation scheme to approximate the solution of the following SODEs

(. u 1 (t) . u 2 (t) \dots \dots \dots . u N - 2 (t) . u N - 1 (t)) = (F 1 (t, u 1 (t), u 2 (t), \dots, u N - 1 (t)) F 2 (t, u 1 (t), u 2 (t), \dots, u N - 1 (t)) \dots \dots \dots F N - 2 (t, u 1 (t), u 2 (t), \dots, u N - 1 (t)) F N - 1 (t, u 1 (t), u 2 (t), \dots, u N - 1 (t)))

(47)

subject to the initial values

(u 1 (0) u 2 (0) \dots \dots \dots u N - 2 (0) u N - 1 (0)) = (g 2 (x_{N, 1}^{L}) g 2 (x_{N, 2}^{L}) \dots \dots \dots g 2 (x_{N, N - 2}^{L}) g 2 (x_{N, N - 1}^{L}))

(48)

In order to approximate the solution of the previous system, we use the SC-GR-C method to deal with the temporal variable t. We choose the approximate solution of the form

u n, M (t) = \sum_{j = 0}^{M} a n, j T T, j (t), n = 1, \dots, N - 1

(49)

Furthermore, the approximation of the time derivative can be computed as

\partial t u n, M (t) = \sum_{j = 0}^{M} a n, j \partial t (T T, j (t)) = \sum_{j = 0}^{M} a n, j T_{T, j}^{(1)} (t), n = 1, \dots, N - 1

(50)

where

T_{T, k}^{(1)} (t)

represents the first time derivative of the shifted Chebyshev polynomials, which can be easily evaluated at any point

t_{K, s}^{T}

(shifted Chebyshev Gauss–Radau points).

Therefore, adopting equations (49)–(50) enables one to write equations (47)–(48) in the form

(\sum_{j = 0}^{M} a 1, j T_{T, j}^{(1)} (t) \sum_{j = 0}^{M} a 2, j T_{T, j}^{(1)} (t) \dots \dots \dots \sum_{j = 0}^{M} a N - 2, j T_{T, j}^{(1)} (t) \sum_{j = 0}^{M} a N - 1, j T_{T, j}^{(1)} (t)) = (F 1 (t, \sum_{j = 0}^{M} a 1, j T T, j (t), \dots, \sum_{j = 0}^{M} a N - 1, j T T, j (t)) F 2 (t, \sum_{j = 0}^{M} a 1, j T T, j (t), \dots, \sum_{j = 0}^{M} a N - 1, j T T, j (t)) \dots \dots \dots F N - 2 (t, \sum_{j = 0}^{M} a 1, j T T, j (t), \dots, \sum_{j = 0}^{M} a N - 1, j T T, j (t)) F N - 1 (t, \sum_{j = 0}^{M} a 1, j T T, j (t), \dots, \sum_{j = 0}^{M} a N - 1, j T T, j (t))),

(51)

(\sum_{j = 0}^{M} a 1, j T T, j (0) \sum_{j = 0}^{M} a 2, j T T, j (0) \dots \dots \dots \sum_{j = 0}^{M} a N - 2, j T T, j (0) \sum_{j = 0}^{M} a N - 1, j T T, j (0)) = (g 2 (x_{N, 1}^{L}) g 2 (x_{N, 2}^{L}) \dots \dots \dots g 2 (x_{N, N - 2}^{L}) g 2 (x_{N, N - 1}^{L}))

(52)

In the proposed technique, the residual of equation (51) has to be enforced to zero at M × (N − 1) collocation points. In other words, we have to collocate equation (51) at the M × (N − 1) shifted Chebyshev Gauss–Radau collocation nodes, which immediately yields

\sum_{j = 0}^{M} a n, j T_{T, j}^{(1)} (t_{M, s}^{T}) = F n (t_{M, s}^{T}, \sum_{j = 0}^{M} a 1, j T T, j (t_{M, s}^{T}), \dots, \times \sum_{j = 0}^{M} a N - 1, j T T, j (t_{M, s}^{T})), n = 1, \dots, N - 1, s = 1, \dots, M

(53)

In virtue of equation (52), we get

\sum_{j = 0}^{M} a n, j T T, j (0) = τ n, n = 1, \dots, N - 1

(54)

The set of previous equations is equivalent to a system of (N − 1) × (M + 1) algebraic equations in the unknowns a_i,j, i = 1,…,N − 1; j = 0,…,M

(κ 1, 0 \dots κ 1, M κ 2, 0 \dots κ 2, M \dots \dots \dots \dots \dots \dots κ N - 2, 0 \dots κ N - 2, M κ N - 1, 0 \dots κ N - 1, M) = (ξ 1, 0 \dots ξ 1, M ξ 2, 0 \dots ξ 2, M \dots \dots \dots \dots \dots \dots ξ N - 2, 0 \dots ξ N - 2, M ξ N - 1, 0 \dots ξ N - 1, M)

(55)

where

κ n, s = {\sum_{j = 0}^{M} a n, j T T, j (0), s = 0, n = 1, \dots, N - 1, \sum_{j = 0}^{M} a n, j T_{T, j}^{(1)} (t_{M, s}^{T}), n = 1, \dots, N - 1, s = 1, \dots, M

(56)

and

ξ l, m = {τ n, s = 0, n = 1, \dots, N - 1, F n (t_{M, s}^{T}, \sum_{j = 0}^{M} a 1, j T T, j (t_{M, s}^{T}), \dots, \sum_{j = 0}^{M} a N - 1, j T T, j (t_{M, s}^{T})), n = 1, \dots, N - 1, s = 1, \dots, M

(57)

The system of algebraic equations can be solved using Newton's iterative method. After the coefficients a_i,j are determined, the approximate solution u_N,M(x,t) can be computed at any value of (x,t) in the given domain from

u N, M (x, t) = \sum_{k = 0}^{M} \sum_{i = 0}^{N} \sum_{j = 0}^{N} a i, k (P L, j (x_{N, i}^{L}) ϖ_{N, i}^{L} h_{j}^{L}) P L, j (x) T T, k (t)

(58)

4. Two-dimensional space fractional percolation equations

In this section, we develop the algorithm for one-dimensional FPEs to handle the following two-dimensional space FPE

\partial t u (x, y, t) = \partial x ν 1 (g 1 (x, y) \partial x ν 2 u (x, y, t)) + \partial y ν 3 (g 2 (x, y) \partial y ν 4 u (x, y, t)) + f (x, y, t), 0 < ν 1, ν 2, ν 3, ν 4 < 1, (x, y, t) \in [0, L 1] \times [0, L 2] \times [0, T]

(59)

with the initial condition

u (x, y, 0) = g 3 (x, y), (x, y) \in [0, L 1] \times [0, L 2]

(60)

and four boundary conditions

u (0, y, t) = g 4 (y, t), u (L 1, y, t) = g 5 (y, t), (y, t) \in [0, L 2] \times [0, T], u (x, 0, t) = g 6 (x, t), u (x, L 2, t) = g 7 (x, t) (x, t) \in [0, L 1] \times [0, T]

(61)

4.1. SL-GL-C scheme for the space variable

The SL-GL-C method will be extended to reduce the solution of the previous FPE into SODEs. Let us expand the dependent variable in the form

u N, M (x, y, t) = \sum_{j = 0}^{M} \sum_{i = 0}^{N} a i, j (t) P L 1, i (x) P L 2, j (y)

(62)

The coefficients a_i,j(t) can be approximated by

a i, j (t) = 1 h_{i}^{L 1} h_{j}^{L 2} \sum_{l = 0}^{N} \sum_{k = 0}^{M} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) \times u (x_{N, l}^{L 1}, y_{M, k}^{L 2}, t)

(63)

Hence we can rewrite u_N,M(x,y,t) as

u N, M (x, y, t) = \sum_{j = 0}^{M} \sum_{i = 0}^{N} \sum_{l = 0}^{N} \sum_{k = 0}^{M} \times (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times P L 1, i (x) P L 2, j (y) u l, k (t)

(64)

where u_N,M(x_N,n,y_M,m,t) = u_n,m(t).

Similarly to equations (34)–(38), we compute the fractional spatial partial derivatives as

\partial x ν 2 u (x, y, t) = \sum_{l = 0}^{N} \sum_{k = 0}^{M} \sum_{j = 0}^{M} \sum_{i = 0}^{N} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times P_{L 1, i}^{(ν 2)} (x) P L 2, j (y) u l, k (t),

(65)

\partial y ν 4 u (x, y, t) = \sum_{l = 0}^{N} \sum_{k = 0}^{M} \sum_{j = 0}^{M} \sum_{i = 0}^{N} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times P_{L 2, j}^{(ν 4)} (y) P L 1, i (x) u l, k (t),

(66)

\partial x ν 1 (\partial x ν 2 u (x, y, t)) = \sum_{l = 0}^{N} \sum_{k = 0}^{M} \sum_{j = 0}^{M} \sum_{i = 0}^{N} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times P_{L 1, i}^{(ν 1 + ν 2)} (x) P L 2, j (y) u l, k (t),

(67)

\partial y ν 3 (\partial y ν 4 u (x, y, t)) = \sum_{l = 0}^{N} \sum_{k = 0}^{M} \sum_{j = 0}^{M} \sum_{i = 0}^{N} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times P_{L 2, j}^{(ν 3 + ν 4)} (y) P L 1, i (x) u l, k (t)

(68)

The previous fractional derivatives can be evaluated at the SL-GL interpolation nodes for n = 1, 2, …, N−1, m = 1, 2, … , M−1, as follows

(\partial x ν 2 u (x, y, t)) y = y_{M, m}^{L 2} x = x_{N, n}^{L 1}, = \sum_{l = 0}^{N} \sum_{k = 0}^{M} ρ_{l, k}^{n, m} u l, k (t)

(69)

(\partial x ν 1 (\partial x ν 2 u (x, y, t))) y = y_{M, m}^{L 2} x = x_{N, n}^{L 1}, = \sum_{l = 0}^{N} \sum_{k = 0}^{M} σ_{l, k}^{n, m} u l, k (t)

(70)

(\partial y ν 4 u (x, y, t)) y = y_{M, m}^{L 2} x = x_{N, n}^{L 1}, = \sum_{l = 0}^{N} \sum_{k = 0}^{M} λ_{l, k}^{n, m} u l, k (t)

(71)

(\partial y ν 3 (\partial y ν 4 u (x, y, t))) y = y_{M, m}^{L 2} x = x_{N, n}^{L 1}, = \sum_{l = 0}^{N} \sum_{k = 0}^{M} μ_{l, k}^{n, m} u l, k (t)

(72)

where

ρ_{l, k}^{n, m} = \sum_{i = 0}^{N} \sum_{j = 0}^{M} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times (P_{L 1, i}^{(ν 2)} (x)) x = x_{N, n}^{L 1} P L 2, j (y_{M, m}^{L 2})

(73)

σ_{l, k}^{n, m} = \sum_{i = 0}^{N} \sum_{j = 0}^{M} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times (P_{L 1, i}^{(ν 1 + ν 2)} (x)) x = x_{N, n}^{L 1} P L 2, j (y_{M, m}^{L 2}),

(74)

λ_{l, k}^{n, m} = \sum_{i = 0}^{N} \sum_{j = 0}^{M} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times (P_{L 2, j}^{(ν 4)} (y)) y = y_{M, m}^{L 2} P L 1, i (x_{N, n}^{L 1}),

(75)

μ_{l, k}^{n, m} = \sum_{i = 0}^{N} \sum_{j = 0}^{M} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times (P_{L 2, j}^{(ν 3 + ν 4)} (y)) y = y_{M, m}^{L 2} P L 1, i (x_{N, n}^{L 1})

(76)

In the proposed SL-GL-C method, the residual of equation (59) is set to zero at (N − 1) × (M − 1) SL-GL points. Moreover, the values of u_0,k(t), u_N,k(t), u_l,0(t) and u_l,N(t) can be given by

u 0, k (t) = g 7 (y_{M, k}^{L 2}, t), u N, k (t) = g 8 (y_{M, k}^{L 2}, t), k = 0, \dots, M, u l, 0 (t) = g 9 (x_{N, l}^{L 1}, t), u l, N (t) = g 10 (x_{N, l}^{L 1}, t), l = 0, \dots, N

(77)

the combination of equations (64)–(77), with equations (59)–(61) leads us to the (N − 1) × (M − 1) SODEs

. u n, m (t) = g 8 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) \sum_{l = 0}^{N} \sum_{k = 0}^{M} ρ_{l, k}^{n, m} u l, k (t) + g 1 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) \times \sum_{l = 0}^{N} \sum_{k = 0}^{M} σ_{l, k}^{n, m} u l, k (t) + g 9 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) \times \sum_{l = 0}^{N} \sum_{k = 0}^{M} λ_{l, k}^{n, m} u l, k (t) + g 2 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) \times \sum_{l = 0}^{N} \sum_{k = 0}^{M} μ_{l, k}^{n, m} u l, k (t) + f (x_{N, n}^{L 1}, y_{M, m}^{L 2}, t), n = 1, \dots, N - 1, m = 1, \dots, M - 1

(78)

subject to the initial conditions

u n, m (0) = g 3 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) = g_{3}^{n, m}, n = 1, \dots, N - 1, m = 1, \dots, M - 1

(79)

where

g 8 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) = (\partial x ν 1 g 1 (x, y)) y = y_{M, m}^{L 2} x = x_{N, n}^{L 1},, g 9 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) = (\partial y ν 3 g 2 (x, y)) y = y_{M, m}^{L 2} x = x_{N, n}^{L 1},

(80)

Finally, we can rearrange equations (78)–(79) to their matrix formulation

(. u 1, 1 (t) \dots . u 1, M - 1 (t) . u 2, 1 (t) : . u 2, M - 1 (t) \dots \dots \dots \dots \dots \dots . u N - 2, 1 (t) : . u N - 2, M - 1 (t) . u N - 1, 1 (t) \dots . u N - 1, M - 1 (t)) = (H 1, 1 (t, u 1, \dots, u N - 1) \dots H 1, M - 1 (t, u 1, \dots, u N - 1) H 2, 1 (t, u 1, \dots, u N - 1) : H 2, M - 1 (t, u 1, \dots, u N - 1) \dots \dots \dots \dots \dots \dots H N - 2, 1 (t, u 1, \dots, u N - 1) : H N - 2, M - 1 (t, u 1, \dots, u N - 1) H N - 1, 1 (t, u 1, \dots, u N - 1) \dots H N - 1, M - 1 (t, u 1, \dots, u N - 1))

(81)

(u 1, 1 (0) \dots u 1, M - 1 (0) u 2, 1 (0) : u 2, M - 1 (0) \dots \dots \dots \dots \dots \dots u N - 2, 1 (0) : u N - 2, M - 1 (0) u N - 1, 1 (0) \dots u N - 1, M - 1 (0)) = (g_{3}^{1, 1} \dots g_{3}^{1, M - 1} g_{3}^{2, 1} : g_{3}^{2, M - 1} \dots \dots \dots \dots \dots \dots g_{3}^{N - 2, 1} : g_{3}^{N - 2, M - 1} g_{3}^{N - 1, 1} \dots g_{3}^{N - 1, M - 1})

(82)

where

H n, m (t, u 1, 1, \dots, u N - 1, M - 1) = g 8 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) \times \sum_{l = 0}^{N} \sum_{k = 0}^{M} ρ_{l, k}^{n, m} u l, k (t) + g 1 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) \times \sum_{l = 0}^{N} \sum_{k = 0}^{M} σ_{l, k}^{n, m} u l, k (t) + g 9 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) \sum_{l = 0}^{N} \sum_{k = 0}^{M} λ_{l, k}^{n, m} u l, k (t) + g 2 (x_{N, n}^{L 1}, y_{M, m}^{L 2}) \sum_{l = 0}^{N} \sum_{k = 0}^{M} μ_{l, k}^{n, m} u l, k (t) + f (x_{N, n}^{L 1}, y_{M, m}^{L 2}, t), n = 1, \dots, N - 1, m = 1, \dots, M - 1

(83)

4.2. SC-GR-C scheme for the time variable

We are interested in using the SC-GR-C method to transform the SODEs (81) subject to equation (82) into a system of algebraic equations. To this end, we approximate the time variable using the following approximation

u n, m, K (t) = \sum_{k = 0}^{K} a_{k}^{n, m} T T, k (t), n = 1, \dots, N - 1, m = 1, \dots, M - 1

(84)

This allows us to immediately obtain the time derivative of the approximate solution in the form

\partial t u n, m (t) = \sum_{k = 0}^{K} a_{k}^{n, m} \partial t (T T, k (t)) = \sum_{k = 0}^{K} a_{k}^{n, m} T_{T, k}^{(1)} (t), n = 1, \dots, N - 1, m = 1, \dots, M - 1

(85)

Combining equations (81)–(85), enables us to write

(\sum_{k = 0}^{K} a_{k}^{1, 1} T_{T, k}^{(1)} (t) \dots \sum_{k = 0}^{K} a_{k}^{1, M - 1} T_{T, k}^{(1)} (t) \sum_{k = 0}^{K} a_{k}^{2, 1} T_{T, k}^{(1)} (t) : \sum_{k = 0}^{K} a_{k}^{2, M - 1} T_{T, k}^{(1)} (t) \dots \dots \dots \dots \dots \dots \sum_{k = 0}^{K} a_{k}^{N - 2, 1} T_{T, k}^{(1)} (t) : \sum_{k = 0}^{K} a_{k}^{N - 2, M - 1} T_{T, k}^{(1)} (t) \sum_{k = 0}^{K} a_{k}^{N - 1, 1} T_{T, k}^{(1)} (t) \dots \sum_{k = 0}^{K} a_{k}^{N - 1, M - 1} T_{T, k}^{(1)} (t)) = (η 1, 1 (t) \dots η 1, M - 1 (t) η 2, 1 (t) \dots η 2, M - 1 (t) \dots \dots \dots \dots \dots \dots η N - 2, 1 (t) \dots η N - 2, M - 1 (t) η N - 1, 1 (t) \dots η N - 1, M - 1 (t)),

(86)

(\sum_{k = 0}^{K} a_{k}^{1, 1} T T, k (0) \dots \sum_{k = 0}^{K} a_{k}^{1, M - 1} T T, k (0) \sum_{k = 0}^{K} a_{k}^{2, 1} T T, k (0) \dots \sum_{k = 0}^{K} a_{k}^{2, M - 1} T T, k (0) \dots \dots \dots \dots \dots \dots \sum_{k = 0}^{K} a_{k}^{N - 2, 1} T T, k (0) \dots \sum_{k = 0}^{K} a_{k}^{N - 2, M - 1} T T, k (0) \sum_{k = 0}^{K} a_{k}^{N - 1, 1} T T, k (0) \dots \sum_{k = 0}^{K} a_{k}^{N - 1, M - 1} T T, k (0)) = (g_{3}^{1, 1} \dots g_{3}^{1, M - 1} g_{3}^{2, 1} : g_{3}^{2, M - 1} \dots \dots \dots \dots \dots \dots g_{3}^{N - 2, 1} : g_{3}^{N - 2, M - 1} g_{3}^{N - 1, 1} \dots g_{3}^{N - 1, M - 1})

(87)

where

η n, m (t) = H n, m (t, \sum_{k = 0}^{K} a_{k}^{1, 1} T T, k (t), \dots, \sum_{k = 0}^{K} a_{k}^{N - 1, M - 1} T T, k (t)) .

In the proposed method, the residual of equation (86) is set to be zero at K × (N − 1) × (M − 1) collocation points.

\sum_{k = 0}^{K} a_{k}^{n, m} T_{T, k}^{(1)} (t_{K, s}^{T}) = η n, m (t_{K, s}^{T}), n = 1, \dots, N - 1, m = 1, \dots, M - 1, s = 1, \dots, K

(88)

according to initial conditions, we have another (N − 1) × (M − 1) algebraic equations

\sum_{k = 0}^{K} a_{k}^{n, m} T T, k (0) = g_{3}^{n, m}, n = 1, \dots, N - 1, m = 1, \dots, M - 1

(89)

Finally, we can merge the previous equations (88) and (89) in the following matrix system form

Ω = Ψ

(90)

where

ℜ n, m s = {\sum_{k = 0}^{K} a_{k}^{n, m} T T, k (0), s = 0, n = 1, \dots, N - 1, m = 1, \dots, M - 1, \sum_{k = 0}^{K} a_{k}^{n, m} T_{T, k}^{(1)} (t_{K, s}^{T}), n = 1, \dots, N - 1, m = 1, \dots, M - 1, s = 1, \dots, K

(93)

and

ξ l, m = {g_{3}^{n, m}, s = 0, n = 1, \dots, N - 1, m = 1, \dots, M - 1; η n, m (t_{K, s}^{T}), n = 1, \dots, N - 1, m = 1, \dots, M - 1, s = 1, \dots, K

(94)

Finally, from equations (91)–(94), we get a system of (N − 1) × (M − 1) × (K + 1) algebraic equations. This system can be solved by any iteration technique. Then, the approximate solution u_N,M,K(x,y,t) can be evaluated as

u N, M, K (x, y, t) = \sum_{r = 0}^{K} \sum_{j = 0}^{M} \sum_{i = 0}^{N} \sum_{l = 0}^{N} \sum_{k = 0}^{M} a_{r}^{l, k} (P L 2, j (y_{M, k}^{L 2}) ϖ_{M, k}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) h_{i}^{L 1} h_{j}^{L 2} \times P L 1, i (x) P L 2, j (y) T T, r (t)

(95)

4.3. Error bound

In this subsection, we present an analytic expression for the error norm of the best approximation for a smooth function u(x,y,t) ∈ Ω ≡ [0,L₁] × [0,L₂] × [0,T] by its expansion

u N, M, K (x, y, t) = \sum_{k = 0}^{K} \sum_{i = 0}^{N} \sum_{j = 0}^{M} π i, j, k P L 1, i (x) P L 2, j (y) T T, k (t)

(96)

where

π i, j, k = \sum_{l = 0}^{N} \sum_{r = 0}^{M} a_{k}^{l, r} (P L 2, j (y_{M, r}^{L 2}) ϖ_{M, r}^{L 2} P L 1, i (x_{N, l}^{L 1}) ϖ_{N, l}^{L 1}) ℎ i ℎ j

(97)

at any value of (x,y,t) in the given domain. Here we provide an upper bound on the error expected in our approximations in the two-dimensional case and for the one-dimensional case we refer the reader to Bhrawy and Zaky (2015a).

Let us first consider the space

Π N, M, K = span {P L 1, i (x) P L 2, j (y) T T, k (t), i = 0, 1, \dots, N, j = 0, 1, \dots, M, k = 0, 1, \dots, K}

In the following analysis it will always be assumed that u_N,M,K(x,y,t) ∈ Π^N,M,K is the best approximation of u(x,y,t), then by the definition of the best approximation, we have

\forall v N, M, K (x, y, t) \in Π N, M, K, ‖ u (x, y, t) - u N, M, K (x, y, t) ‖ \infty \leq ‖ u (x, y, t) - v N, M, K (x, y, t) ‖ \infty

(98)

It turns out that the previous inequality is also true if v_N,M,K(x,y,t) denotes the interpolating polynomial for u(x,y,t) at points

(x_{N, i}^{L 1}, y_{M, j}^{L 2}, t_{K, k}^{T})

, where

t_{K, k}^{T}, (0 \leq k \leq K)

are the shifted Chebyshev Gauss–Radau points and

x_{N, i}^{L 1}, (0 \leq i \leq N)

and

y_{M, j}^{L 2}, (0 \leq j \leq M)

are the shifted Legendre Gauss–Lobatto collocation points. Then by similar procedures to those in Bhrawy and Zaky (2015a), we can write

u (x, y, t) - v N, M, K (x, y, t) = \partial N + 1 u (ε 1, y, t) \partial x N + 1 (N + 1)! Π_{i = 0}^{N} (x - x_{N, i}^{L 1}) + \partial M + 1 u (x, ε 2, t) \partial y M + 1 (M + 1)! Π_{j = 0}^{M} (y - y_{M, j}^{L 2}) + \partial K + 1 u (x, y, ε 3) \partial t K + 1 (K + 1)! Π_{k = 0}^{K} (t - t_{K, k}^{T}) - \partial N + M + K + 3 u (ε 1', ε 2', ε 3') \partial x N + 1 \partial y M + 1 \partial t K + 1 (M + 1)! (N + 1)! (K + 1)! \times Π_{i = 0}^{N} (x - x_{N, i}^{L 1}) Π_{j = 0}^{M} (y - y_{M, j}^{L 2}) Π_{k = 0}^{K} (t - t_{K, k}^{T})

(99)

where ɛ₁, ɛ′₁ ∈ [0,L₁], ɛ₂, ɛ′₂ ∈ [0,L₂] and ɛ₃, ɛ′₃ ∈ [0,T], and thereby we obtain that

‖ u (x, y, t) - v N, M, K (x, y, t) ‖ \infty \leq max (x, y, t) \in Ω | \partial N + 1 u (ε 1, y, t) \partial x N + 1 | ‖ Π_{i = 0}^{N} (x - x_{N, i}^{L 1}) ‖ \infty (N + 1)! + max (x, y, t) \in Ω | \partial M + 1 u (x, ε 2, t) \partial y M + 1 | ‖ Π_{j = 0}^{M} (y - y_{M, j}^{L 2}) ‖ \infty M + 1)! + max (x, y, t) \in Ω | \partial K + 1 u (x, y, ε 3) \partial t K + 1 | ‖ Π_{k = 0}^{K} (t - t_{K, k}^{T}) ‖ \infty (K + 1)! + max (x, y, t) \in Ω | \partial N + M + K + 3 u (ɛ 1', ɛ 2', ɛ 3') \partial x N + 1 \partial y M + 1 \partial t K + 1 | \times ‖ Π_{i = 0}^{N} (x - x_{N, i}^{L 1}) ‖ \infty ‖ Π_{j = 0}^{M} (y - y_{M, j}^{L 2}) ‖ \infty ‖ Π_{k = 0}^{K} (t - t_{K, k}^{T}) ‖ \infty (K + 1)! (N + 1)! (M + 1)!

(100)

Since u(x,y,t) is a smooth function on Ω, then there exists the constants C₁, C₂, C₃ and C₄, such that

max (x, y, t) \in Ω | \partial N + 1 u (x, y, t) \partial x N + 1 | \leq C 1, max (x, y, t) \in Ω | \partial M + 1 u (x, y, t) \partial y M + 1 | \leq C 2, max (x, y, t) \in Ω | \partial K + 1 u (x, y, t) \partial t K + 1 | \leq C 3, max (x, y, t) \in Ω | \partial N + M + K + 3 u (x, t) \partial x N + 1 \partial y M + 1 \partial t K + 1 | \leq C 4

(101)

Here, we try to reduce the error in the polynomial approximation by minimizing the factor

‖ Π_{i = 0}^{N} (x - x_{N, i}^{L 1}) ‖ \infty

. Let us use the one-to-one mapping

x = L 12 (z + 1)

between the intervals [ − 1,1] and [0,L₁] to deduce that

min x_{N, i}^{L 1} \in [0, L 1] max 0 \leq x \leq L 1 | Π_{i = 0}^{N} (x - x_{N, i}^{L 1}) | = min z i \in [- 1, 1] max - 1 \leq z \leq 1 | Π_{i = 0}^{N} L 12 (z - z i) | = (L 12) N + 1 min z i \in [- 1, 1] max - 1 \leq z \leq 1 | Π_{i = 0}^{N} (z - z i) | = (L 12) N + 1 max - 1 \leq z \leq 1 | P N + 1 (z) κ N |

(102)

where

κ N = (2 N)! 2 N (N!) 2

is the leading coefficient of P_N+1(z) and z_i are the Legendre Gauss–Lobatto collocation points. It is a well-known fact (Canuto et al., 2006), that the Legendre and Chebyshev polynomials satisfy

max - 1 \leq z \leq 1 | P N + 1 (z) | \leq 1 and max - 1 \leq z \leq 1 | T K + 1 (z) | \leq 1

We can obtain

‖ Π_{i = 0}^{N} (x - x_{N, i}^{L 1}) ‖ \infty \leq (L 1 / 2) N + 1 κ N

(103)

Also we can obtain

‖ Π_{k = 0}^{K} (t - t_{K, k}^{T}) ‖ \infty \leq (T / 2) K + 1 κ' K

(104)

where

κ' K = 2 K

is the leading coefficient of T_K+1(z).

Combining equations (100), (101), (103) and (104), yields

‖ u (x, y, t) - u N, M, K (x, y, t) ‖ \infty \leq C' 1 (L 1 / 2) N + 1 κ N (N + 1)! + C' 2 (L 2 / 2) M + 1 κ M (M + 1)! + C' 3 (T / 2) K + 1 κ' K (K + 1)! + C' 4 (L 1 / 2) N + 1 (L 2 / 2) M + 1 (T / 2) K + 1 κ N κ M κ' K (N + 1)! (M + 1)! (K + 1)!

(105)

Hence, an upper bound of the absolute error is obtained for the approximate solutions. The convergence of the proposed method depends basically on the above error bound.

5. Numerical results and comparisons

After the construction of the spectral collocation methods, we now carry out some numerical examples of 1 + 1 and 2 + 1 FPEs to study the performance of the methods, which are discussed and developed in the current paper, and to compare our results with those proposed in Chen et al. (2011b, 2013, 2014). We divide this section into two main parts. In the first one, we introduce two examples of 1 + 1 FPEs. Then, we deal with three numerical examples of 2 + 1 FPEs in the second part.

5.1. Examples for one-dimensional space fractional percolation equations

Example 1

Let us start with the following FPE (Chen et al., 2011b)

\partial t u (x, t) = \partial x ((30 - x 2) \partial 0.5 x u (x, t)) - e - t (45 \sqrt{x} - 3.5 x 2.5) Γ (3) Γ (2.5) - e - t x 2, (x, t) \in [0, 1] \times [0, 1]

(106)

subject to

The exact solution is u(x,t) = e^−t x².

Recently, Chen et al. (2011b) introduced approximate solutions to the previous problem with various choices of N and M using a novel implicit finite difference (NIFD) method. Now, we list in Table 1 a comparison based on the maximum absolute error (M_E) between the numerical method introduced in Chen et al. (2011b), and our results with choices of N and M, which are less than those in Chen et al. (2011b). We observe that the error, using the present method, decays exponentially with an increase in the collocation nodes N,M. In Table 2, we see the high accuracy of the present method based on absolute error at N = M = 12.

Table 1.

Comparing M_E of the proposed method with the NIFD (Chen et al., 2011b) method for problem (106).

N = M	SL-GL-C Method			NIFD (Chen et al., 2011b)
N = M	4	8	12	40	80
M _E	1.24 × 10⁻⁵	1.94 × 10⁻¹¹	6.66 × 10⁻¹⁶	1.52 × 10⁻⁷	7.68 × 10⁻⁸

Table 2.

Absolute error at N = M = 12 and various choices of (x,t) for problem (106).

x	t	E(x,t)	x	t	E(x,t)	x	t	E(x,t)
0.1	0.1	3.33 × 10⁻¹⁶	0.1	0.5	3.05 × 10⁻¹⁶	0.1	1	1.11 × 10⁻¹⁶
0.2		2.22 × 10⁻¹⁶	0.2		3.33 × 10⁻¹⁶	0.2		1.11 × 10⁻¹⁶
0.3		3.05 × 10⁻¹⁶	0.3		4.44 × 10⁻¹⁶	0.3		1.38 × 10⁻¹⁶
0.4		4.16 × 10⁻¹⁶	0.4		5.13 × 10⁻¹⁶	0.4		1.38 × 10⁻¹⁶
0.5		3.33 × 10⁻¹⁶	0.5		5.07 × 10⁻¹⁶	0.5		5.89 × 10⁻¹⁷
0.6		3.54 × 10⁻¹⁶	0.6		4.61 × 10⁻¹⁶	0.6		8.50 × 10⁻¹⁷
0.7		3.05 × 10⁻¹⁶	0.7		3.89 × 10⁻¹⁶	0.7		1.04 × 10⁻¹⁶
0.8		3.33 × 10⁻¹⁶	0.8		5.27 × 10⁻¹⁶	0.8		1.39 × 10⁻¹⁶
0.9		1.67 × 10⁻¹⁶	0.9		3.33 × 10⁻¹⁶	0.9		0
1		1.11 × 10⁻¹⁶	1		2.22 × 10⁻¹⁶	1		1.11 × 10⁻¹⁶

A space–time graph of the error function between the exact and approximate solutions at N = 12, is sketched in Figure 1. In addition, to confirm the high accuracy and convergence of the present scheme, in Figure 2, we plot the logarithmic graph of M_E(log₁₀ M_E) at various values of N(N = M), which shows that the proposed method provides an accurate approximation and yields exponential convergence rates.

Figure 1.

Absolute error of problem (106).

Figure 2.

Convergence of problem (106).

Example 2

Consider the following FPE (Chen et al., 2011b)

\partial t u (x, t) = \partial 0.9 x (\partial 0.8 x u (x, t)) - e - t (x 2 + Γ (0.3) Γ (1.3) x 0.3), (x, t) \in [0, 1] \times [0, 1]

(108)

subject to

u (x, 0) = x 2, x \in [0, 1], u (0, t) = 0, u (1, t) = e - t, t \in [0, 1]

(109)

The exact solution is u(x,t) = e^−t x².

In Table 3, the numerical results based on M_E, obtained using the proposed algorithm, are compared with those presented by the NIFD method (Chen et al., 2011b). From this table, we observe that the proposed method is more accurate than the NIFD method (Chen et al., 2011b). A comparison of the x-directional curves of the exact solution and the numerical solution at t = 0.1, 0.5, 0.9 is shown in Figure 3. Meanwhile, the absolute error in the x-directional and at time value t = 0 and N = M = 14, is plotted in Figure 4. In Figure 5, we sketch the logarithmic graph of M_E (log₁₀M_E) at various values of N(N = M) to demonstrate the convergence of the present method.

Table 3.

Comparing M_E of the proposed method and the NIFD (Chen et al., 2011b) method for problem (108).

N = M	SL-GL-C Method			NIFD (Chen et al., 2011b)
N = M	6	10	14	40	80
M _E	1.77 × 10⁻⁸	8.55 × 10⁻¹⁵	6.66 × 10⁻¹⁶	7.55 × 10⁻⁴	3.82 × 10⁻⁴

Figure 3.

x-Directional curves of the approximate and exact solutions of problem (108) at N = M = 14.

Figure 4.

Absolute error of problem (108) at t = 0 and N = M = 14.

Figure 5.

Convergence of problem (108).

5.2. Examples for two-dimensional space fractional percolation equations

Example 3

Consider the two-dimensional FPE (Chen et al., 2013)

\partial t u (x, y, t) = \partial β 1 x ((2 - x 2) \partial α 1 x u (x, y, t)) + \partial β 2 y ((2 - y 2) \partial α 2 y u (x, y, t)) - e - t x 2 y 2 - e - t y 2 Γ (3) Γ (3 - α 1) (2 Γ (3 - α 1) Γ (3 - α 1 - β 1) x 2 - α 1 - β 1 - 2 Γ (5 - α 1) Γ (5 - α 1 - β 1) x 4 - α 1 - β 1) - e - t x 2 Γ (3) Γ (3 - α 2) (2 Γ (3 - α 2) Γ (3 - α 2 - β 2) y 2 - α 2 - β 2 - 2 Γ (5 - α 2) Γ (5 - α 2 - β 2) y 4 - α 2 - β 2), \times (x, y, t) \in [0, 1] \times [0, 1] \times [0, 1]

(110)

subject to the initial condition

u (x, y, 0) = x 2 y 2, (x, y) \in [0, 1] \times [0, 1]

(111)

and the boundary conditions

u (0, y, t) = 0, u (1, y, t) = e - t y 2, (y, t) \in [0, 1] \times [0, 1], u (x, 0, t) = 0, u (x, 1, t) = e - t x 2, (x, t) \in [0, 1] \times [0, 1]

(112)

The exact solution is u(x,t) = e^−t x²,y².

In spite of using small values of N, M and K, we obtain highly accurate solutions to this problem using the present method. In Table 4, we present the maximum absolute error (M_E at t = 1) obtained by our method, at different values of N, M and K, along with a comparison with the results given by using the implicit finite difference method (Chen et al., 2013), at the special values α₁ = 0.5, β₁ = 1, α₂ = 0.5 and β₂ = 1. In Table 5, we list the absolute error at N = M = K = 8 for problem (110), where α₁ = 0.5, β₁ = 1, α₂ = 0.5 and β₂ = 1. We also plot the space graph of absolute error in Figure 6, with the parameter values listed in its caption. In order to demonstrate the exponential convergence of the method, in Figure 7, we sketch the logarithmic graphs of M_E at various values of N(N = M = K).

Table 4.

Comparing M_E of the proposed method and the NIFD (Chen et al., 2013) method for problem (110) at α₁ = 0.5, β₁ = 1, α₂ = 0.5 and β₂ = 1.

N = M = K	SL-GL-C Method			NIFD (Chen et al., 2013)
N = M = K	4	6	8	80	160
M _E	2.94 × 10⁻⁶	6.38 × 10⁻⁹	1.44 × 10⁻¹¹	–	–
M_E (at t = 1)	6.00 × 10⁻⁶	1.57 × 10⁻⁸	6.82 × 10⁻¹²	1.17 × 10⁻³	5.92 × 10⁻⁴

Table 5.

Absolute error for problem (110) with α₁ = 0.5, β₁ = 1, α₂ = 0.5, β₂ = 1 at N = M = K = 8.

x	y	t	E(x,y,t)	x	y	t	E(x,y,t)	x	y	t	E(x,y,t)
0	0.5	0.1	1.60 × 10⁻¹³	0	0.5	0.5	3.31 × 10⁻¹³	0	0.5	1	6.05 × 10⁻¹³
0.1			2.84 × 10⁻¹³	0.1			8.05 × 10⁻¹³	0.1			5.40 × 10⁻¹³
0.2			2.62 × 10⁻¹³	0.2			1.18 × 10⁻¹²	0.2			8.12 × 10⁻¹³
0.3			6.44 × 10⁻¹⁴	0.3			1.44 × 10⁻¹²	0.3			1.43 × 10⁻¹²
0.4			1.02 × 10⁻¹³	0.4			1.99 × 10⁻¹²	0.4			1.36 × 10⁻¹²
0.5			4.63 × 10⁻¹³	0.5			2.30 × 10⁻¹²	0.5			1.30 × 10⁻¹²
0.6			1.04 × 10⁻¹²	0.6			2.23 × 10⁻¹²	0.6			1.37 × 10⁻¹²
0.7			1.51 × 10⁻¹²	0.7			2.17 × 10⁻¹²	0.7			2.16 × 10⁻¹³
0.8			1.83 × 10⁻¹²	0.8			1.75 × 10⁻¹²	0.8			2.24 × 10⁻¹²
0.9			2.24 × 10⁻¹²	0.9			4.12 × 10⁻¹³	0.9			2.36 × 10⁻¹²
1			1.60 × 10⁻¹³	1			3.31 × 10⁻¹³	1			6.05 × 10⁻¹³

Figure 6.

The space graph of the absolute error for problem (110) at α₁ = 0.5, β₁ = 1, α₂ = 0.5, β₂ = 1 and N = M = K = 8.

Figure 7.

Convergence of problem (110).

Example 4

Consider the two-dimensional FPE (Chen et al., 2014)

\partial t u (x, y, t) = \partial x (Γ (3 - ν 1) (3 - 2 x) 2 \partial ν 1 x u (x, y, t)) + \partial y (Γ (4 - ν 2) (4 - y) 6 \partial ν 2 y u (x, y, t)) - e - t x 2 y 3 + e - t x 2 y 2 - ν 2 ((ν 2 - 3) (4 - y) + y) + e - t x 1 - ν 1 y 3 ((ν 1 - 2) (3 - 2 x) + 2 x), (x, y, t) \in [0, 1] \times [0, 1] \times [0, 1]

(113)

subject to

u (x, y, 0) = x 2 y 3, (x, y) \in [0, 1] \times [0, 1]

(114)

and

u (0, y, t) = 0, u (1, y, t) = e - t y 3, (y, t) \in [0, 1] \times [0, 1], u (x, 0, t) = 0, u (x, 1, t) = e - t x 2, (x, t) \in [0, 1] \times [0, 1]

(115)

The exact solution of equation (113) is u(x,y,t) = e^−t x²y³.

In order to show that our method is more accurate than the implicit difference method (IDM) (Chen et al., 2014), in Table 6 the numerical results for M_E using the proposed algorithm are compared with those presented by the IDM (Chen et al., 2014). From this table, we observe that the proposed method is more accurate than the IDM (Chen et al., 2014).

Table 6.

Comparing M_E of the proposed method and IDM (Chen et al., 2014) for problem (113).

N = M = k	SL-GL-C Method			IDM (Chen et al., 2014)
N = M = k	4	6	8	80	160
M _E	5.04 × 10⁻⁶	9.82 × 10⁻⁹	1.19 × 10⁻¹¹	–	–
M_E (at t = 1)	5.90 × 10⁻⁶	1.67 × 10⁻⁸	1.48 × 10⁻¹¹	2.07 × 10⁻³	1.06 × 10⁻³

Example 5

Finally, consider the two-dimensional FPE

\partial t u (x, y, t) = \partial x 0.8 (\partial x 0.7 u (x, y, t)) + \partial y 0.9 (\partial y 0.6 u (x, y, t)) - e - t x 2 y 2 - Γ (3) Γ (1.5) e - t x 0.5 y 2 - Γ (3) Γ (1.5) e - t x 2 y 0.5, (x, y, t) \in [0, 1] \times [0, 1] \times [0, 1]

(116)

subject to the initial condition

u (x, y, 0) = x 2 y 2, (x, y) \in [0, 1] \times [0, 1]

(117)

and the boundary conditions

u (0, y, t) = 0, u (1, y, t) = e - t y 2, (y, t) \in [0, 1] \times [0, 1], u (x, 0, t) = 0, u (x, 1, t) = e - t x 2, (x, t) \in [0, 1] \times [0, 1]

(118)

The exact solution of equation (116) is u(x,y,t) = e^−t x²y².

Table 7 lists the M_E of u(x,y,t) for problem (116) with various choices of N, M, and K. The numerical results presented in this table show that the results are very accurate for small values of N, M and K. Figure 8 demonstrates that the absolute error E(x,y,t) is very small for even the small number of grid points taken.

Table 7.

Maximum absolute error for problem (116).

N = M = K	4	6	8
M _E	3.11 × 10⁻⁶	6.57 × 10⁻⁹	7.83 × 10⁻¹²
M_E (at t = 1)	5.78 × 10⁻⁶	1.59 × 10⁻⁸	1.44 × 10⁻¹¹

Figure 8.

t-Directional curve of the absolute error of problem (116) at x = y = 0.5 and N = M = K = 8.

6. Conclusion

As one of the few articles dealing with FPEs, we have proposed and developed a spectral collocation method to obtain accurate numerical solutions for one- and two-dimensional FPEs. The core of the proposed method was to discretize the FPE in the spatial direction by the SL-GL-C method, along with a new treatment for the subjected conditions, to create SODEs of the unknown coefficients of spectral expansion in the time direction. An efficient numerical integration process for SODEs was investigated based on the SC-GR-C method. The proposed method has been successfully applied to numerically solve one- and two-dimensional FPEs. The main advantage of the present approach is, by adding a few terms of the shifted SL-GL and SC-GR collocation nodes, a highly accurate solution of the problem can be obtained. Comparison between our approximate solutions and the approximate solutions achieved by other methods were presented to demonstrate the high accuracy and validity of our method.

Footnotes

Acknowledgements

The authors are very grateful to the reviewers for carefully reading the paper and for their comments and suggestions which have improved the paper.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Abdelkawy

Ahmed

Sanchez

(2015a) newblock A method based on Legendre pseudo-spectral approximations for solving inverse problems of parabolic types equations. Mathematical Sciences Letters 4(1): 81–90.

Abdelkawy

Ezz-Eldien

Amim

AZM

(2015b) A Jacobi spectral collocation scheme for solving Abel's integral equations. Progress in Fractional Differentiation and Applications 1(3): 1–14.

Al-Khaled

Alquran

(2014) An approximate solution for a fractional model of generalized Harry Dym equation. Mathematical Sciences 8(4): 125–130.

Alipour

Rostamy

Baleanu

(2012) Solving multidimensional FOCPs with inequality constraint by BPs operational matrices. Journal of Vibration and Control 19(16): 2523–2540.

Baleanu

Diethelm

Scalas

Trujillo

(2012) Fractional Calculus Models and Numerical Methods, New York: World Scientific Publishing Company.

Bhrawy

(2014) An efficient Jacobi pseudospectral approximation for nonlinear complex generalized Zakharov system. Applied Mathematics and Computation 247: 30–46.

Bhrawy

Doha

Baleanu

Ezz-Eldien

(2015a) A spectral tau algorithm based on Jacobi operational matrix for numerical solution of time fractional diffusion-wave equations. Journal of Computational Physics 293: 142–156.

Bhrawy

Zaky

(2015a) A method based on the Jacobi tau approximation for solving multi-term time–space fractional partial differential equations. Journal of Computational Physics 281: 876–895.

Bhrawy

Zaky

(2015b) Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation. Nonlinear Dynamics 80(1): 101–116.

10.

Bhrawy

Zaky

Machado

(2015b) Efficient Legendre spectral tau algorithm for solving the two-sided space–time Caputo fractional advection–dispersion equation. Journal of Vibration and Control. Epub ahead of print 22 January 2015. DOI: 10.1177/1077546314566835.

11.

Canuto

Hussaini

Quarteroni

Zang

(2006) Spectral Methods: Fundamentals in Single Domains, New York: Springer-Verlag.

12.

Chen

Liu

Anh

Turner

(2011a) Numerical simulation for the variable-order Galilei invariant advection diffusion equation with a nonlinear source term. Applied Mathematics and Computation 217(12): 5729–5742.

13.

Chen

Liu

Anh

(2011b) A novel implicit finite difference method for the one-dimensional fractional percolation equation. Numerical Algorithms 56(4): 517–535.

14.

Chen

Liu

Burrage

(2014) Numerical simulation of a new two-dimensional variable-order fractional percolation equation in non-homogeneous porous media. Computers & Mathematics with Applications 68(12): 2133–2141.

15.

Chen

Liu

Turner

Anh

(2013) An implicit numerical method for the two-dimensional fractional percolation equation. Applied Mathematics and Computation 219(9): 4322–4331.

16.

Chou H, Lee B and Chen C (2006) The transient infiltration process for seepage flow from cracks. Advances in Subsurface Flow and Transport: Eastern and Western Approaches III. Eos Trans. AGU, Western Pacific Meeting, 2006.

17.

Dehghan

Hamedi

Khosravian-Arab

(2014) A numerical scheme for the solution of a class of fractional variational and optimal control problems using the modified Jacobi polynomials. Journal of Vibration and Control. Epub ahead of print 3 September 2014. DOI: 10.1177/1077546314543727.

18.

Ding

Xiao

(2010) Weighted finite difference methods for a class of space fractional partial differential equations with variable coefficients. Journal of Computational and Applied Mathematics 233(8): 1905–1914.

19.

Doha

Bhrawy

Abdelkawy

Van Gorder

(2014a) Jacobi–Gauss–Lobatto collocation method for the numerical solution of 1 + 1 nonlinear Schrödinger equations. Journal of Computational Physics 261: 244–255.

20.

Doha

Bhrawy

Baleanu

Hafez

(2014b) A new Jacobi rational–Gauss collocation method for numerical solution of generalized pantograph equations. Applied Numerical Mathematics 77: 43–54.

21.

Doha

Bhrawy

Ezz-Eldien

(2011) Efficient Chebyshev spectral methods for solving multi-term fractional orders differential equations. Applied Mathematical Modelling 35(12): 5662–5672.

22.

Doha

Bhrawy

Ezz-Eldien

(2012) A new Jacobi operational matrix: An application for solving fractional differential equations. Applied Mathematical Modelling 36(10): 4931–4943.

23.

Eslahchi

Dehghan

Parvizi

(2014) Application of the collocation method for solving nonlinear fractional integro-differential equations. Journal of Computational and Applied Mathematics 257: 105–128.

24.

Ezz-Eldien

Doha

Baleanu

Bhrawy

(2015) A numerical approach based on Legendre orthonormal polynomials for numerical solutions of fractional optimal control problems. Journal of Vibration and Control. Epub ahead of print 3 March 2015. DOI: 10.1177/1077546315573916.

25.

Garrappa

Popolizio

(2011) On the use of matrix functions for fractional partial differential equations. Mathematics and Computers in Simulation 81(5): 1045–1056.

26.

Giona

Roman

(1992) Fractional diffusion equation for transport phenomena in random media. Physica A: Statistical Mechanics and its Applications 185(1): 87–97.

27.

(1998) Approximate analytical solution for seepage flow with fractional derivatives in porous media. Computer Methods in Applied Mechanics and Engineering 167(1): 57–68.

28.

Heinrichs

(1989) Spectral methods with sparse matrices. Numerische Mathematik 56(1): 25–41.

29.

Irandoust-Pakchin

Dehghan

Abdi-Mazraeh

Lakestani

(2014) Numerical solution for a class of fractional convection–diffusion equations using the flatlet oblique multiwavelets. Journal of Vibration and Control 20(6): 913–924.

30.

Jafari

Tajadodi

(2015) Numerical solutions of the fractional advection-dispersion equation. Progress in Fractional Differentiation and Applications 1(1): 37–45.

31.

Kayedi-Bardeh

Eslahchi

Dehghan

(2014) A method for obtaining the operational matrix of fractional Jacobi functions and applications. Journal of Vibration and Control 20(5): 736–748.

32.

Khalil

Khan

(2014) A new method based on Legendre polynomials for solutions of the fractional two-dimensional heat conduction equation. Computers & Mathematics with Applications 67(10): 1938–1953.

33.

Kirchner

Feng

Neal

(2000) Fractal stream chemistry and its implications for contaminant transport in catchments. Nature 403(6769): 524–527.

34.

Lazo

Torres

(2014) The Legendre condition of the fractional calculus of variations. Optimization 63(8): 1157–1165.

35.

Deng

(2007) Remarks on fractional derivatives. Applied Mathematics and Computation 187(2): 777–784.

36.

Machado

Galhano

Trujillo

(2013) Science metrics on fractional calculus development since 1966. Fractional Calculus and Applied Analysis 16(2): 479–500.

37.

Magin

(2006) Fractional Calculus in Bioengineering, West Redding: Begell House Publishers.

38.

Meerschaert

Tadjeran

(2006) Finite difference approximations for two-sided space-fractional partial differential equations. Applied Numerical Mathematics 56(1): 80–90.

39.

Miller

Ross

(1993) An Introduction to the Fractional Calculus and Fractional Differential Equations, New York: John Wiley & Sons.

40.

Mohebbi

Abbaszadeh

Dehghan

(2014) High-order difference scheme for the solution of linear time fractional Klein–Gordon equations. Numerical Methods for Partial Differential Equations 30(4): 1234–1253.

41.

Pinto

CMA

Carvalho

(2015) Fractional modeling of typical stages in HIV epidemics with drug-resistance. Progress in Fractional Differentiation and Applications 1(2): 111–122.

42.

Pooseh

Almeida

Torres

(2013) A discrete time method to the first variation of fractional order variational functionals. Central European Journal of Physics 11(10): 1262–1267.

43.

Raja

MAZ

Chaudhary

(2015) Two-stage fractional least mean square identification algorithm for parameter estimation of CARMA systems. Signal Processing 107: 327–339.

44.

Rostamy

Alipour

Jafari

Baleanu

(2013) Solving multi-term orders fractional differential equations by operational matrices of BPs with convergence analysis. Romanian Reports in Physics 65(2): 334–349.

45.

Wang

(2014) Fast alternating-direction finite difference methods for three-dimensional space-fractional diffusion equations. Journal of Computational Physics 258: 305–318.

46.

Hesthaven

Chen

(2015) A parareal method for time-fractional differential equations. Journal of Computational Physics 293: 173–183.