Numerical Solution of the Fractional Partial Differential Equations by the Two-Dimensional Fractional-Order Legendre Functions

Abstract

A numerical method is presented to obtain the approximate solutions of the fractional partial differential equations (FPDEs). The basic idea of this method is to achieve the approximate solutions in a generalized expansion form of two-dimensional fractional-order Legendre functions (2D-FLFs). The operational matrices of integration and derivative for 2D-FLFs are first derived. Then, by these matrices, a system of algebraic equations is obtained from FPDEs. Hence, by solving this system, the unknown 2D-FLFs coefficients can be computed. Three examples are discussed to demonstrate the validity and applicability of the proposed method.

1. Introduction

Fractional partial differential equations play a significant role in modeling physical and engineering processes. Therefore, there is an urgent need to develop efficient and fast convergent methods for FPDEs. Recently, several different techniques, including Adomian’s decomposition method (ADM) [1, 2], homotopy perturbation method (HPM) [3–5], variational iteration method (VIM) [6–8], spectral methods [9–13], orthogonal polynomials method [14–17], and wavelets method [18–21] have been presented and applied to solve FPDEs.

The method based on the orthogonal functions is a wonderful and powerful tool for solving the FDEs and has enjoyed many successes in this realm. The operational matrix of fractional integration has been determined for some types of orthogonal polynomials, such as Chebyshev polynomials [16], Legendre polynomials [22], Laguerre polynomials [23–25], and Jacobi polynomials [26]. Moreover, the operational matrix of fractional derivative for Chebyshev polynomials [9] and Legendre polynomials [9, 14] also has been derived. However, since these polynomials using integer power series to approximate fractional ones, it cannot accurately represent properties of fractional calculus. Recently, Rida and Yousef [27] presented a fractional extension of the classical Legendre polynomials by replacing the integer order derivative in Rodrigues formula with fractional order derivatives. The defect is that the complexity of these functions made them unsuitable for solving FDEs. Subsequently, Kazem et al. [28] presented the orthogonal fractional order Legendre functions based on shifted Legendre polynomials to find the numerical solution of FDEs and drew a conclusion that their method is accurate, effective, and easy to implement.

Benefiting from their “exponential-convergence” property when smooth solutions are involved, spectral methods have been widely and effectively used for the numerical solution of partial differential equations. The basic idea of spectral methods is to expand a function into sets of smooth global functions, called the trial functions. Because of their special properties, the orthogonal polynomials are usually chosen to be trial functions. Spectral methods can obtain very accurate approximations for a smooth solution while only need a few degrees of freedom. Recently, Chebyshev spectral method [9], Legendre spectral method [10], and adaptive pseudospectral method [11] were proposed for solving fractional boundary value problems. Moreover, generalized Laguerre spectral algorithms and Legendre spectral Galerkin method were developed by Baleanu et al. [12] and Bhrawy and Alghamdi [13] for fractional initial value problems, respectively.

Motivated and inspired by the ongoing research in orthogonal polynomials methods and spectral methods, we construct two-dimensional fractional-order Legendre functions and derive the operational matrices of integration and derivative for the solution of FPDEs. To the best of the authors’ knowledge, such approach has not been employed for solving FPDEs.

The rest of the paper is organized as follows. In Section 2, we introduce some mathematical preliminaries of the fractional calculus theory and fractional-order Legendre functions. In Section 3, a basis of 2D-FLFs is defined and some properties are given. Section 4 is devoted to the operational matrices of fractional derivative and integration for 2D-FLFs. Some numerical examples are presented in Section 5. Finally, we conclude the paper with some remarks.

2. Preliminaries and Notations

2.1. Fractional Calculus Theory

Some necessary definitions and Lemma of the fractional calculus theory [29, 30] are listed here for our subsequent development.

Definition 1. A real function ℎ(𝑡), 𝑡>0, is said to be in the space 𝐶_𝜇, 𝜇∈𝑅, if there exists a real number 𝑝>𝜇, such that ℎ(𝑡)=𝑡^𝑝ℎ₁(𝑡), where ℎ₁(𝑡)∈𝐶(0,∞), and it is said to be in the space 𝐶^𝑛_𝜇 if and only if ℎ^(𝑛)∈𝐶_𝜇, 𝑛∈𝑁.

Definition 2. Riemann-Liouville fractional integral operator (𝐽^𝛼) of order 𝛼≥0, of a function 𝑓∈𝐶_𝜇, 𝜇≥−1 is defined as 𝐽^𝛼𝑓(𝑡)=1Γ(𝛼)^𝑡₀(𝑡−𝜏)^𝛼−1𝑓(𝜏)𝐽𝑑𝜏,𝑡>0,₀𝑓(𝑡)=𝑓(𝑡), (1) where Γ(𝛼) is the well-known Gamma function. Some properties of the operator 𝐽^𝛼 can be found, for example, in [29, 30].

Definition 3. The fractional derivative of 𝑓(𝑥) in the Caputo sense is defined as 𝐷^𝛼𝑓(𝑥)=1Γ()×∫𝑚−𝛼_𝑥0𝑓^(𝑚)(𝜉)(𝑥−𝜉)^{𝛼−𝑚+1}(),𝑑𝑑𝜉,𝛼>0,𝑚−1<𝛼<𝑚_𝑚𝑓(𝑥)𝑑𝑥^𝑚,𝛼=𝑚, (2) where 𝑓∶𝑅→𝑅, 𝑥→𝑓(𝑥) denotes a continuous (but not necessarily differentiable) function.

Lemma 4. Let 𝑛−1<𝛼≤𝑛, 𝑛∈𝑁, 𝑡>0, ℎ∈𝐶^𝑛_𝜇, 𝜇≥−1. Then 𝐽^𝛼𝐷^𝛼ℎ(𝑡)=ℎ(𝑡)−^𝑛−1_𝑘=0ℎ^(𝑘)0⁺𝑡^𝑘.𝑘! (3)

2.2. Fractional-Order Legendre Functions

In this section, we introduce the fractional-order Legendre functions which were first proposed by Kazem et al. [28]. The normalized eigenfunctions problem for FLFs is 𝑥−𝑥^1+𝛼𝐿^′𝛼_𝑖(𝑥)^′+𝛼²𝑖(𝑖+1)𝑥^𝛼−1𝐿^𝛼_𝑖(𝑥)(),=0,𝑥∈0,1 (4) which is a singular Sturm-Liouville problem. The fractional-order Legendre polynomials, denoted by FL^𝛼_𝑖(𝑥), are defined on the interval [0,1] and can be determined with the aid of following recurrence formulae: FL^𝛼₀(𝑥)=1,FL^𝛼₁(𝑥)=2𝑥^𝛼−1,FL^𝛼_𝑖+1(𝑥)=(2𝑖+1)2𝑥^𝛼−1𝑖+1FL^𝛼_𝑖(𝑥)−𝑖𝑖+1FL^𝛼_𝑖−1(𝑥),𝑖=1,2,…, (5) and the analytic form of FL^𝛼_𝑖(𝑥) of degree 𝑖 is given by FL^𝛼_𝑖(𝑥)=^𝑖_𝑠=0𝑏_𝑠,𝑖𝑥^𝑠𝛼,𝑏_𝑠,𝑖=(−1)^𝑖+𝑠(𝑖+𝑠)!(𝑖−𝑠)!(𝑠!)², (6) where FL^𝛼_𝑖(0)=(−1)^𝑖 and FL^𝛼_𝑖(1)=1. The orthogonality condition is ¹₀FL^𝛼_𝑛(𝑥)FL^𝛼_𝑚(𝑥)𝜔(𝑥)1𝑑𝑥=()𝛼𝛿2𝑛+1^𝑛𝑚, (7) where 𝜔(𝑥)=𝑥^𝛼−1 is the weight function and 𝛿 is the Kronecker delta. For more details, please see [28].

3. 2D-FLFs

In this section, the definitions and theorems of 2D-FLFs are given by Liu’s method described in [31].

3.1. Definitions and Properties of the 2D-FLFs

Definition 5. Let {FL^𝛼_𝑛(𝑥)}^∞_𝑛=0 be the fractional Legendre polynomials on [0,1]; we call {FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗(𝑦)}^∞_𝑖,𝑗=0 the two-dimensional fractional Legendre polynomials on [0,1]×[0,1].

Theorem 6. The basis {𝐹𝐿^𝛼_𝑖(𝑥)𝐹𝐿^𝛽_𝑗(𝑦)}^∞_𝑖,𝑗=0 is orthogonal on [0,1]×[0,1] with the weight function 𝜔(𝑥,𝑦)=𝜔(𝑥)𝜔(𝑦)=𝑥^𝛼−1𝑦^𝛽−1.

Proof. Let 𝑖≠𝑚 or 𝑗≠𝑛¹₀¹₀𝜔𝑥,𝑦FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦FL^𝛼_𝑚(𝑥)FL^𝛽_𝑛𝑦=𝑑𝑥𝑑𝑦₁₀𝜔(𝑥)FL^𝛼_𝑖(𝑥)FL^𝛼_𝑚(𝑥)×𝑑𝑥₁₀𝜔𝑦FL^𝛽_𝑗𝑦FL^𝛽_𝑛𝑦𝑑𝑦=0. (8)

Theorem 7. Consider ¹₀¹₀𝜔𝑥,𝑦𝐹𝐿^𝛼_𝑖(𝑥)𝐹𝐿^𝛽_𝑗𝑦²=1𝑑𝑥𝑑𝑦()𝛼12𝑖+1𝛽,2𝑗+1₁₀₁₀𝜔𝑥,𝑦𝐹𝐿^𝛼_𝑖(𝑥)𝐹𝐿^𝛽_𝑗𝑦²=𝑑𝑥𝑑𝑦₁₀𝜔(𝑥)𝐹𝐿^𝛼_𝑖(𝑥)²𝑑𝑥¹₀𝜔𝑦𝐹𝐿^𝛽_𝑗𝑦²=1𝑑𝑦()𝛼12𝑖+1𝛽.2𝑗+1 (9)

3.2. 2D-FLFs Expansion

Definition 8. A function of two independent variables 𝑓(𝑥,𝑦) which is integrable in square [0,1]×[0,1] can be expanded as 𝑓𝑥,𝑦=^∞∞_𝑖=0𝑗=0𝑎_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦, (10) where 𝑎_𝑖𝑗=()×2𝑖+12𝑗+1𝛼𝛽₁₀₁₀𝑓𝑥,𝑦𝜔𝑥,𝑦FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦𝑑𝑥𝑑𝑦. (11)

Theorem 9. If the series ∑^∞_𝑖=0∑^∞_𝑗=0𝑎_𝑖𝑗𝐹𝐿^𝛼_𝑖(𝑥)𝐹𝐿^𝛽_𝑗(𝑦) converges uniformly to 𝑓(𝑥,𝑦) on the square [0,1]×[0,1], then we have 𝑎_𝑖𝑗=()×2𝑖+12𝑗+1𝛼𝛽₁₀₁₀𝑓𝑥,𝑦𝜔𝑥,𝑦𝐹𝐿^𝛼_𝑖(𝑥)𝐹𝐿^𝛽_𝑗𝑦𝑑𝑥𝑑𝑦. (12)

Proof. By multiplying 𝜔(𝑥,𝑦)FL^𝛼_𝑛(𝑥)FL^𝛽_𝑚(𝑦) on both sides of (10), where 𝑛 and 𝑚 are fixed and integrating termwise with regard to 𝑥 and 𝑦 on [0,1]×[0,1], then ¹₀¹₀𝑓𝑥,𝑦𝜔𝑥,𝑦FL^𝛼_𝑛(𝑥)FL^𝛽_𝑚𝑦=𝑑𝑥𝑑𝑦_{∞∞𝑖=0𝑗=0}𝑎_𝑖𝑗¹₀¹₀𝜔𝑥,𝑦FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦×FL^𝛼_𝑛(𝑥)FL^𝛽_𝑚𝑦=𝑑𝑥𝑑𝑦_{∞∞𝑖=0𝑗=0}𝑎_𝑖𝑗¹₀𝜔(𝑥)FL^𝛼_𝑖(𝑥)FL^𝛼_𝑛(𝑥)×𝑑𝑥₁₀𝜔𝑦FL^𝛽_𝑗𝑦FL^𝛽_𝑚𝑦𝑑𝑦=𝑎_𝑛𝑚¹₀𝜔(𝑥)FL^𝛼_𝑛(𝑥)²𝑑𝑥¹₀𝜔𝑦FL^𝛽_𝑚𝑦²𝑑𝑦=𝑎_𝑛𝑚1()𝛼12𝑛+1()𝛽.2𝑚+1 (13) Finally one can get (11).

If the infinite series in (10) is truncated, then it can be written as 𝑓𝑥,𝑦≈^𝑚𝑚_𝑖=0^′_𝑗=0𝑎_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦=𝐶^𝑇Ψ𝑥^𝛼,𝑦^𝛼, (14) where 𝐶 and Ψ(𝑥^𝛼,𝑦^𝛽) are given by 𝐶=𝑐_0,0,𝑐_0,1,…,𝑐_{0,𝑚′−1},𝑐_1,0,𝑐_1,1𝑐,…,_{1,𝑚′−1},…,𝑐_𝑚−1,0,𝑐_𝑚−1,1,…,𝑐_{𝑚−1,𝑚′−1}^𝑇,Ψ𝑥^𝛼,𝑦^𝛽=𝜓_0,0,𝜓_0,1,…,𝜓_{0,𝑚′−1},𝜓_1,0,𝜓_1,1𝜓,…,_{1,𝑚′−1},…,𝜓_𝑚−1,0,𝜓_𝑚−1,1,…,𝜓_{𝑚−1,𝑚′−1}^𝑇, (15)(16) where 𝜓_𝑖𝑗=FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗(𝑦), 𝑖=0,1,…,𝑚, and 𝑗=0,1,…,𝑚^′.

According to the definition of FLFs, one can find that fractional Legendre polynomials are identical to Legendre polynomials shifted to [0,1] when using the transform 𝑥^𝛼→𝑥, 𝑦^𝛽→𝑦. Therefore, in a similar method described in [31], we can easily get the convergence and stability theorems of proposed method.

Lemma 10. If the function 𝑓(𝑥,𝑦) is a continuous function on [0,1]×[0,1] and the series ∑^∞_𝑖=0∑^∞_𝑗=0𝑎_𝑖𝑗𝐹𝐿^𝛼_𝑖(𝑥)𝐹𝐿^𝛽_𝑗(𝑦) converges uniformly to 𝑓(𝑥,𝑦), then ∑^∞_𝑖=0∑^∞_𝑗=0𝑎_𝑖𝑗𝐹𝐿^𝛼_𝑖(𝑥)𝐹𝐿^𝛽_𝑗(𝑦) is the 2D-FLFs expansion of 𝑓(𝑥,𝑦).

Proof (by contradiction). Let 𝑓𝑥,𝑦=^∞∞_𝑖=0𝑗=0𝑏_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦,𝑓𝑥,𝑦∼^∞∞_𝑖=0𝑗=0𝑎_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦. (17) Then there is at least one coefficient such that 𝑎_𝑛𝑚≠𝑏_𝑛𝑚. However, 𝑏_𝑛𝑚=()×2𝑛+1)(2𝑚+1𝛼𝛽₁₀₁₀𝑓𝑥,𝑦𝜔𝑥,𝑦FL^𝛼_𝑛(𝑥)FL^𝛽_𝑚𝑦𝑑𝑥𝑑𝑦=𝑎_𝑛𝑚. (18)

Lemma 11. If two continuous functions defined on [0,1]×[0,1] have the identical 2D-FLFs expansions, then these two functions are identical.

Proof. Suppose that 𝑓(𝑥,𝑦) and 𝑔(𝑥,𝑦) can be expanded by 2D-FLFs as follows: 𝑓𝑥,𝑦∼^∞∞_𝑖=0𝑗=0𝑎_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦,𝑔𝑥,𝑦∼^∞∞_𝑖=0𝑗=0𝑎_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦. (19) By subtracting the above two equations with each other, one has 𝑓𝑥,𝑦−𝑔𝑥,𝑦∼^∞∞_𝑖=0𝑗=0𝑎_𝑖𝑗−𝑎_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦=0=^∞∞_𝑖=0𝑗=00×FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦. (20) Then Lemma 11 can be proved.

Theorem 12. If the 2D-FLFs expansion of a continuous function 𝑓(𝑥,𝑦) converges uniformly, then the 2D-FLFs expansion converges to the function 𝑓(𝑥,𝑦).

Proof. Theorem 12 can be proved by Theorems 7 and 9.

Theorem 13. If the sum of the absolute values of the 2D-FLFs coefficients of a continuous function 𝑓(𝑥,𝑦) forms a convergent series, then the 2D-FLFs expansion is absolutely uniformly convergent, and converges to the function 𝑓(𝑥,𝑦).

Proof. Consider ^∞∞_𝑖=0𝑗=0𝑎_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦≤^∞∞_𝑖=0𝑗=0𝑎_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦≤^∞∞_𝑖=0𝑗=0𝑎_𝑖𝑗. (21) Then ∑^∞_𝑖=0∑^∞_𝑗=0𝑎_𝑖𝑗FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗(𝑦) converges uniformly to the function 𝑓(𝑥,𝑦).

Theorem 14. If a continuous function 𝑓(𝑥,𝑦), defined on [0,1]×[0,1], has bounded mixed partial derivative 𝐷_𝑥^2𝛼𝐷_𝑦^2𝛽𝑓(𝑥,𝑦), then the 2D-FLFs expansion of the function converges uniformly to the function.

Proof. Let 𝑓(𝑥,𝑦) be a function defined on [0,1]×[0,1] such that 𝐷_𝑥^2𝛼𝐷_𝑦^2𝛽𝑓𝑥,𝑦≤𝑀, (22) where 𝑀 is a positive constant and 𝑎_𝑖𝑗=()×2𝑖+12𝑗+1𝛼𝛽₁₀₁₀𝑓𝑥,𝑦𝜔𝑥,𝑦FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦𝑑𝑥𝑑𝑦. (23) By employing the transform 𝑋=2𝑥^𝛼−1 and 𝑌=2𝑦^𝛽−1, one can obtain 𝑎_𝑖𝑗=2𝑖+122𝑗+12¹₋₁¹₋₁𝑓(𝑋,𝑌)𝑝_𝑖(𝑋)𝑝_𝑗(𝑌)𝑑𝑋𝑑𝑌. (24) Consequently, in a similar method described in [31], Theorem 14 can be proved.

4. Operational Matrices of 2D-FLFs

4.1. Integration Operational Matrices of 2D-FLFs

Lemma 15. The Riemann-Liouville fractional integration of order 𝛾>0 of the 2D-FLFs 𝜓_𝑖𝑗 can be obtained in the form of 𝐽^𝛾_𝑥𝜓_𝑖𝑗𝑥^𝛼,𝑦^𝛽=𝐹𝐿^𝛽_𝑗𝑦^𝑖_𝑠=0𝑏_𝑠𝑖Γ(1+𝑠𝛼)Γ𝑥1+𝑠𝛼+𝛾^{𝑠𝛼+𝛾}. (25)

Proof. Consider 𝐽^𝛾_𝑥𝜓_𝑖𝑗𝑥^𝛼,𝑦^𝛽=𝐽^𝛾_𝑥FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦=𝐽^𝛾_𝑥FL^𝛼_𝑖(𝑥)FL^𝛽_𝑗𝑦=𝐽^𝛾_𝑥^𝑖_𝑠=0𝑏_𝑠𝑖𝑥^𝑠𝛼FL^𝛽_𝑗𝑦=FL^𝛽_𝑗𝑦^𝑖_𝑠=0𝑏_𝑠𝑖Γ(1+𝑠𝛼)Γ𝑥1+𝑠𝛼+𝛾^{𝑠𝛼+𝛾}. (26)

Lemma 16. Let 𝛾>0; then one has ¹₀¹₀𝐽^𝛾_𝑥𝜓_𝑖𝑗𝜓_{𝑖′𝑗′}𝜔=𝑥,𝑦𝑑𝑥𝑑𝑦_𝑖∑∑^𝑖_𝑠=0^′_𝑠′=0𝑏_𝑠𝑖𝑏_{𝑠′𝑖′}𝑠+𝑠^′×Γ()+1𝛼+𝛾1+𝑠𝛼Γ11+𝑠𝛼+𝛾2𝑗+1𝛽,𝑗=𝑗^′0,𝑗≠𝑗^′. (27)

Proof. Using previous Lemma 15 and (6), one can have ¹₀¹₀𝐽^𝛾_𝑥𝜓_𝑖𝑗𝑥^𝛼,𝑦^𝛽𝜓_{𝑖′𝑗′}𝑥^𝛼,𝑦^𝛽𝜔=𝑥,𝑦𝑑𝑥𝑑𝑦₁₀₁₀𝜔𝑥,𝑦FL^𝛼_𝑖′𝑥^𝛼FL^𝛽_𝑗′𝑦^𝛽FL^𝛽_𝑗𝑦^𝛽×^𝑖_𝑠=0𝑏_𝑠𝑖Γ(1+𝑠𝛼)Γ𝑥1+𝑠𝛼+𝛾^{𝑠𝛼+𝛾}=𝑑𝑥𝑑𝑦₁₀₁₀𝜔𝑦FL^𝛽_𝑗𝑦^𝛽FL^𝛽_𝑗′𝑦^𝛽×^𝑖𝑖_𝑠=0^′_𝑠′=0𝑏_𝑠𝑖𝑏_{𝑠′𝑖′}Γ(1+𝑠𝛼)Γ1+𝑠𝛼+𝛾×𝑥^{(𝑠+𝑠′+1)𝛼+𝛾−1}=𝑑𝑥𝑑𝑦₁₀𝜔𝑦FL^𝛽_𝑗𝑦^𝛽FL^𝛽_𝑗′𝑦^𝛽×¹₀^𝑖𝑖_𝑠=0^′_𝑠′=0𝑏_𝑠𝑖𝑏_{𝑠′𝑖′}Γ(1+𝑠𝛼)Γ1+𝑠𝛼+𝛾×𝑥^{(𝑠+𝑠′+1)𝛼+𝛾−1}=𝑑𝑥𝑑𝑦_{𝑖𝑖𝑠=0′𝑠′=0}𝑏_𝑠𝑖𝑏_{𝑠′𝑖′}𝑠+𝑠^′Γ()+1𝛼+𝛾1+𝑠𝛼Γ1+𝑠𝛼+𝛾×¹₀𝜔𝑦FL^𝛽_𝑗𝑦^𝛽FL^𝛽_𝑗′𝑦^𝛽=𝑑𝑦_𝑖∑∑^𝑖_𝑠=0^′_𝑠′=0𝑏_𝑠𝑖𝑏_{𝑠′𝑖′}𝑠+𝑠^′×Γ()+1𝛼+𝛾1+𝑠𝛼Γ11+𝑠𝛼+𝛾2𝑗+1𝛽,𝑗=𝑗^′0,𝑗≠𝑗^′. (28)

Theorem 17. Let Ψ(𝑥^𝛼,𝑦^𝛽) be the 2D-FLFs vector defined in (16); then one has 𝐽^𝛾_𝑥Ψ𝑥^𝛼,𝑦^𝛽≃𝐏^𝛾_𝑥Ψ𝑥^𝛼,𝑦^𝛽, (29) where 𝐏^𝛾_𝑥 is the 𝑚𝑚^′×𝑚𝑚^′ operational matrix of Riemann-Liouville fractional integration of order 𝛾>0, and has the form as follows: 𝐏^𝛾_𝑥=𝐸_0,0𝐸_0,1⋯𝐸_0,𝑚−1𝐸_1,0𝐸_1,1⋯𝐸_1,𝑚−1𝐸⋮⋮⋱⋮_𝑚−1,0𝐸_𝑚−1,1⋯𝐸_{𝑚−1,𝑚−1}, (30) in which 𝐸_𝑖,𝑖′ is 𝑚^′×𝑚^′ matrix and the elements are defined as follows: 𝐸_𝑖,𝑖′=𝐼^𝑖𝑖_𝑠=0^′_𝑠′=0𝑏_𝑠𝑖𝑏_{𝑠′𝑖′}2𝑖^′+1𝛼𝑠+𝑠^′Γ()+1𝛼+𝛾1+𝑠𝛼Γ,1+𝑠𝛼+𝛾𝑖,𝑖^′=0,1,…,𝑚−1, (31) and 𝐼 is 𝑚^′×𝑚^′ identity matrix.

Proof. Using (29) and orthogonality property of FLFs, one can get 𝐏^𝛾_𝑥=𝐽^𝛾_𝑥Ψ𝑥^𝛼,𝑦^𝛽,Ψ^𝑇𝑥^𝛼,𝑦^𝛽𝐻⁻¹, (32) where ⟨𝐽^𝛾_𝑥Ψ(𝑥^𝛼,𝑦^𝛽),Ψ^𝑇(𝑥^𝛼,𝑦^𝛽)⟩ and 𝐻⁻¹ are two 𝑚^′×𝑚^′ matrices defined as 𝐽^𝛾_𝑥Ψ𝑥^𝛼,𝑦^𝛽,Ψ^𝑇𝑥^𝛼,𝑦^𝛽=¹₀¹₀𝐽^𝛾_𝑥Ψ_𝑘𝑥^𝛼,𝑦^𝛽×Ψ_𝑘′𝑥^𝛼,𝑦^𝛽𝜔𝑥,𝑦𝑑𝑥𝑑𝑦^𝑚𝑚′_𝑘,𝑘′=^𝑖𝑖_𝑠=0^′_𝑠′=0𝑏_𝑠𝑖𝑏_{𝑠′𝑖′}Γ(1+𝑠𝛼)Γ×11+𝑠𝛼+𝛾𝑠+𝑠^′1+1𝛼+𝛾2𝑗+1𝛽^𝑚;𝑚′_{𝑖,𝑖′;𝑗=𝑗′},𝐻⁻¹=2𝑖^′+12𝑗+1𝛼𝛽^𝑚;𝑚′_{𝑖,𝑖′;𝑗=𝑗′}. (33)
Now by substituting above equations in (32), Theorem 12 can be proved.

In a similar way as previous, one can obtain the operational matrix of Riemann-Liouville fractional integration with respect to variable 𝑦.

Theorem 18. Let Ψ(𝑥^𝛼,𝑦^𝛽) be the 2D-FLFs vector defined in (16); one has 𝐽^𝛾_𝑦Ψ𝑥^𝛼,𝑦^𝛽≃𝐏^𝛾_𝑦Ψ𝑥^𝛼,𝑦^𝛽, (34) where 𝐏^𝛾_𝑦 is the 𝑚𝑚^′×𝑚𝑚^′ operational matrix of Riemann-Liouville fractional integration of order 𝛾>0, and has the form as follows: 𝐏^𝛾_𝑦=𝐸𝑂⋯𝑂𝑂𝐸⋯𝑂⋮⋮⋱⋮𝑂𝑂⋯𝐸, (35) in which 𝐸 is 𝑚^′×𝑚^′ matrix and the elements are defined as follows: 𝐸_𝑗,𝑗′=^𝑗𝑗_𝑟=0^′_𝑟′=0𝑏_𝑟𝑗𝑏_{𝑟′𝑗′}2𝑗^′+1𝛽𝑟+𝑟^′Γ+1𝛽+𝛾1+𝑟𝛽Γ,1+𝑟𝛽+𝛾𝑗,𝑗^′=0,1,…,𝑚^′−1. (36)

4.2. Derivative Operational Matrices of 2D-FLFs

Lemma 19. The FLFs Caputo fractional derivative of 𝛾>0 can be obtained in the form of 𝐷^𝛾_𝑥𝜓_𝑖𝑗𝑥^𝛼,𝑦^𝛽=𝐹𝐿^𝛽_𝑗𝑦^𝛽𝑖_𝑠=0𝑏^′_𝑠𝑖Γ(1+𝑠𝛼)Γ𝑥1+𝑠𝛼−𝛾^{𝑠𝛼−𝛾}, (37) where 𝑏^′_𝑠,𝑖=0 when 𝑠𝛼∈𝑁₀ and 𝑠𝛼<𝛾 in other case 𝑏^′_𝑠,𝑖=𝑏_𝑠,𝑖.

Proof. Consider 𝐷^𝛾_𝑥𝜓_𝑖𝑗𝑥^𝛼,𝑦^𝛽=𝐷^𝛾_𝑥FL^𝛼_𝑖𝑥^𝛼FL^𝛽_𝑗𝑦^𝛽=FL^𝛽_𝑗𝑦^𝛽𝐷^𝛾_𝑥FL^𝛼_𝑖𝑥^𝛼=𝐷^𝛾_𝑥^𝑖_𝑠=0𝑏_𝑠𝑖𝑥^𝑠𝛼FL^𝛽_𝑗𝑦^𝛽=FL^𝛽_𝑗𝑦^𝛽𝑖_𝑠=0𝑏^′_𝑠𝑖Γ(1+𝑠𝛼)Γ𝑥1+𝑠𝛼−𝛾^{𝑠𝛼−𝛾}. (38)

Lemma 20. Let 𝛾>0, 𝛼∉𝑁; then one has ¹₀¹₀𝐷^𝛾_𝑥𝜓_𝑖𝑗𝜓_{𝑖′𝑗′}𝜔=𝑥,𝑦𝑑𝑥𝑑𝑦_𝑖∑∑^𝑖_𝑠=0^′_𝑠′=0𝑏_𝑠𝑖𝑏_{𝑠′𝑖′}𝑠+𝑠^′×Γ()+1𝛼−𝛾1+𝑠𝛼Γ11+𝑠𝛼−𝛾2𝑗+1𝛽,𝑗=𝑗^′0,𝑗≠𝑗^′. (39)

Proof. Using previous Lemma 19 and (6), one can have ¹₀¹₀𝐷^𝛾_𝑥𝜓_𝑖𝑗𝑥^𝛼,𝑦^𝛽𝜓_{𝑖′𝑗′}𝑥^𝛼,𝑦^𝛽𝜔=𝑥,𝑦𝑑𝑥𝑑𝑦₁₀₁₀𝜔𝑥,𝑦FL^𝛼_𝑖′𝑥^𝛼FL^𝛽_𝑗′𝑦^𝛽FL^𝛽_𝑗𝑦^𝛽×^𝑖_𝑠=0𝑏^′_𝑠𝑖Γ(1+𝑠𝛼)Γ𝑥1+𝑠𝛼−𝛾^{𝑠𝛼−𝛾}=𝑑𝑥𝑑𝑦₁₀₁₀𝜔𝑦FL^𝛽_𝑗𝑦^𝛽FL^𝛽_𝑗′𝑦^𝛽×^𝑖𝑖_𝑠=0^′_𝑠′=0𝑏^′_𝑠𝑖𝑏_{𝑠′𝑖′}Γ(1+𝑠𝛼)Γ1+𝑠𝛼−𝛾×𝑥^{(𝑠+𝑠′+1)𝛼−𝛾−1}=𝑑𝑥𝑑𝑦₁₀𝜔𝑦FL^𝛽_𝑗𝑦^𝛽FL^𝛽_𝑗′𝑦^𝛽×¹₀^𝑖𝑖_𝑠=0^′_𝑠′=0𝑏^′_𝑠𝑖𝑏_{𝑠′𝑖′}Γ(1+𝑠𝛼)Γ1+𝑠𝛼−𝛾×𝑥^{(𝑠+𝑠′+1)𝛼−𝛾−1}=𝑑𝑥𝑑𝑦_{𝑖𝑖𝑠=0′𝑠′=0}𝑏^′_𝑠𝑖𝑏_{𝑠′𝑖′}𝑠+𝑠^′Γ()+1𝛼−𝛾1+𝑠𝛼Γ1+𝑠𝛼−𝛾×¹₀𝜔𝑦FL^𝛽_𝑗𝑦^𝛽FL^𝛽_𝑗′𝑦^𝛽=𝑑𝑦_𝑖∑∑^𝑖_𝑠=0^′_𝑠′=0𝑏^′_𝑠𝑖𝑏_{𝑠′𝑖′}𝑠+𝑠^′×Γ()+1𝛼−𝛾1+𝑠𝛼Γ11+𝑠𝛼−𝛾2𝑗+1𝛽,𝑗=𝑗^′0,𝑗≠𝑗^′. (40)

Theorem 21. Let Ψ(𝑥^𝛼,𝑦^𝛽) be the 2D-FLFs vector defined in (16); one has 𝐷^𝛾_𝑥Ψ𝑥^𝛼,𝑦^𝛽≃𝐃^𝛾_𝑥Ψ𝑥^𝛼,𝑦^𝛽, (41) where 𝐃^𝛾_𝑥 is the 𝑚𝑚^′×𝑚𝑚^′ operational matrix of Caputo fractional derivative of order 𝛾>0, and has the form as follows: 𝐃^𝛾_𝑥=𝐹𝑂𝑂⋯𝑂_1,0𝐹𝑂⋯𝑂⋮⋮⋱⋮_𝑚−1,0𝐹_𝑚−1,1⋯𝑂 (42) in which 𝐹_𝑖,𝑖′ is 𝑚^′×𝑚^′ matrix and the elements are defined as follows: 𝐹_𝑖,𝑖′=𝐼^𝑖𝑖_𝑠=0^′_𝑠′=0𝑏^′_𝑠𝑖𝑏_{𝑠′𝑖′}2𝑖^′+1𝛼𝑠+𝑠^′Γ()+1𝛼−𝛾1+𝑠𝛼Γ,1+𝑠𝛼−𝛾𝑖,𝑖^′=0,1,…,𝑚−1, (43) and 𝐼 is a 𝑚^′×𝑚^′ identity matrix.

Proof. Using (41) and the orthogonality property of FLFs, one can have 𝐃^𝛾_𝑥=𝐷^𝛾_𝑥Ψ𝑥^𝛼,𝑦^𝛽,Ψ^𝑇𝑥^𝛼,𝑦^𝛽𝐻⁻¹, (44) where ⟨𝐷^𝛾_𝑥Ψ(𝑥^𝛼,𝑦^𝛽),Ψ^𝑇(𝑥^𝛼,𝑦^𝛽)⟩ and 𝐻⁻¹ are two 𝑚𝑚^′×𝑚𝑚^′ matrices defined as 𝐷^𝛾_𝑥Ψ𝑥^𝛼,𝑦^𝛽,Ψ^𝑇𝑥^𝛼,𝑦^𝛽=¹₀¹₀𝐷^𝛾_𝑥Ψ_𝑘𝑥^𝛼,𝑦^𝛽×Ψ_𝑘′𝑥^𝛼,𝑦^𝛽𝜔𝑥,𝑦𝑑𝑥𝑑𝑦^𝑚𝑚′_𝑘,𝑘′=^𝑖𝑖_𝑠=0^′_𝑠′=0𝑏^′_𝑠𝑖𝑏_{𝑠′𝑖′}𝑠+𝑠^′×Γ()+1𝛼−𝛾1+𝑠𝛼Γ11+𝑠𝛼−𝛾2𝑗+1𝛽^𝑚;𝑚′_{𝑖,𝑖′;𝑗=𝑗′}𝐻⁻¹=2𝑖^′+12𝑗+1𝛼𝛽^𝑚;𝑚′_{𝑖,𝑖′;𝑗=𝑗′}. (45)
Now by substituting above equations in (44), Theorem 21 can be proved.

In a similar way as above, one can get Caputo fractional derivative of order 𝛾>0 with respect to variable 𝑦.

Theorem 22. Let Ψ(𝑥^𝛼,𝑦^𝛽) be the 2D-FLFs vector defined in (16); one can have 𝐷^𝛾_𝑦Ψ𝑥^𝛼,𝑦^𝛽≃𝐃^𝛾_𝑦Ψ𝑥^𝛼,𝑦^𝛽, (46) where 𝐃^𝛾_𝑦 is the 𝑚𝑚^′×𝑚𝑚^′ operational matrix of Caputo fractional derivative of order 𝛾>0, and has the form as follows: 𝐃^𝛾_𝑦=𝐹𝑂⋯𝑂𝑂𝐹⋯𝑂⋮⋮⋱⋮𝑂𝑂⋯𝐹, (47) in which 𝐹 is 𝑚^′×𝑚^′ matrix and the elements are defined as follows: 𝐹_𝑗,𝑗′=^𝑗𝑗_𝑟=0^′_𝑟′=0𝑏^′_𝑟𝑗𝑏_{𝑟′𝑗′}2𝑗^′+1𝛽𝑟+𝑟^′Γ+1𝛽+𝛾1+𝑟𝛽Γ,1+𝑟𝛽+𝛾𝑗,𝑗^′=0,1,…,𝑚^′−1. (48)

5. Applications and Results

Consider the following FPDEs: 𝐷^𝛼_𝑥𝑢(𝑥,𝑡)+𝐷^𝛽_𝑡𝑢()[𝑢([𝑢(()(],𝑥,𝑡+𝑁𝑥,𝑡)]+𝐿𝑥,𝑡)]=𝑔𝑥,𝑡,𝛼,𝛽∈0,1 (49) where 𝐿 and 𝑁 are linear operator and nonlinear operator; respectively. 𝐷^𝛼 and 𝐷^𝛽 are the Caputo fractional derivatives of order 𝛼 and 𝛽, respectively; 𝑔 is a known analytic function.

By employing operator 𝐽^𝛽_𝑡 on both sides of (49) and then using the Lemma 4, one can have 𝑢(𝑥,𝑡)+𝐽^𝛽_𝑡𝐷^𝛼_𝑥𝑢()()(−𝑥,𝑡+𝑁𝑢𝑥,𝑡+𝐿𝑢𝑥,𝑡)_{𝑚−1𝑘=0}𝑢^(𝑘)()𝑥𝑥,0^𝑘𝑘!−𝐽^𝛽_𝑡𝑔(𝑥,𝑡)=0. (50)

We first express unknown function 𝑢(𝑥,𝑡) and derivative term 𝐷^𝛼_𝑥𝑢(𝑥,𝑡) as 𝑢(𝑥,𝑡)=𝐶^𝑇Ψ𝑥^𝛼,𝑡^𝛽,𝐷^𝛼_𝑥𝑢(𝑥,𝑡)=𝐶^𝑇𝐃^𝛼_𝑥Ψ𝑥^𝛼,𝑡^𝛽. (51)

Now for the nonlinear part, by employing the nonlinear term approximation method described in [32] and then by using transform 𝑥→𝑥^𝛼, 𝑡→𝑡^𝛽, one can get the 2D-FLFs expansion of nonlinear term as 𝑁𝑢(𝑥,𝑡)=𝑁^𝑇Ψ𝑥^𝛼,𝑡^𝛽, (52) where 𝑁^𝑇 is coefficient matrix of nonlinear term which must be computed and its order is 𝑚𝑚^′×𝑚𝑚^′.

For the linear part, we have 𝐿𝑢(𝑥,𝑡)=𝐿^𝑇Ψ𝑥^𝛼,𝑡^𝛽, (53) where 𝐿 is a matrix of order 𝑚𝑚^′×𝑚𝑚^′.

After substituting (51)–(53) into (50), one can obtain 𝐶^𝑇+𝐶^𝑇𝐃^𝛼_𝑥+𝑁^𝑇+𝐿^𝑇𝐏^𝛽_𝑦−𝐶^𝑇_guess=0. (54) According to the Wu’s [33] technology for determining the initial iteration value, the initial iteration value is chosen as 𝑢_guess=∑^𝑚−1_𝑘=0𝑢^(𝑘)(𝑥,0)(𝑥^𝑘/𝑘!)+𝐽^𝛽_𝑡{𝑔(𝑥,𝑡)}=𝐶^𝑇_guessΨ(𝑥^𝛼,𝑡^𝛽). The coefficient matrix 𝐶^𝑇 can be computed by using the MATLAB function fsolve( ) or the method described in [34].

Now, the present method is applied to solve the linear and nonlinear FPDEs, and their results are compared with the solution of other methods. The accuracy of our approach is estimated by the following error functions: 𝑒_𝑗=𝑢_exact𝑗−𝑢_approx𝑗,𝑒=𝑢_exact−𝑢_approx,‖𝑒‖_𝐿∞=max_{1≤𝑗≤𝑁}𝑒_𝑗,‖𝑒‖_𝐿2=^𝑁_𝑗=1𝑒_𝑗²,‖𝑒‖_RMS=1𝑁^𝑁_𝑗=1𝑒_𝑗². (55)

Example 23. Consider the one-dimensional linear inhomogeneous fractional Burger’s equation [35]: 𝜕^𝛽𝑢(𝑥,𝑡)𝜕𝑡^𝛽+𝜕𝑢(𝑥,𝑡)−𝜕𝜕𝑥²𝑢(𝑥,𝑡)𝜕𝑥²=2𝑡^2−𝛽Γ3−𝛽+2𝑥−2,0<𝛽≤1, (56) with the initial condition 𝑢(𝑥,0)=𝑥² and the exact solution being 𝑢(𝑥,𝑡)=𝑥²+𝑡².
By employing 2D-FLFs method, one can get 𝐶^𝑇𝐼+𝐃^𝛼_𝑥−𝐃^𝛼_𝑥²𝐏^𝛽_𝑡=𝐶^𝑇_guess, (57) where 𝛼=1. Then we can get 𝐶^𝑇=𝐶^𝑇_guessinv(𝐼+(𝐃^𝛼_𝑥−(𝐃^𝛼_𝑥)²)𝐏^𝛽_𝑡).
Figures 1(a) and 1(b) show the numerical results for 𝛽=0.25 with 𝑚=3, 𝑚^′=9 and 𝛽=0.5 with 𝑚=3, 𝑚^′=5, respectively. It should be found that the accuracy of 2D-FLFs method is very high while only a small number of 2D-FLFs are needed.

(a)

(b)

(c)

(d)

Figure 1 

Numerical results for Example 23 when 𝛽 = 0.25, 0.50.

Example 24. Consider nonlinear fractional Klein-Gordon equation [36, 37]: 𝐷^𝛽_𝑡𝑢(𝑥,𝑡)−𝐷^𝛼_𝑥𝑢(𝑥)+𝑢³(𝑥)(),(],=𝑔𝑥,𝑡𝑥≥0,𝑡>0,𝛼,𝛽∈1,2 (58) subject to the initial conditions 𝑢(𝑥,0)=0,𝑢_𝑡(𝑥,0)=0, (59) and 𝑔(𝑥,𝑡)=Γ(𝛽+1)𝑥^𝛼−Γ(𝛼+1)𝑡^𝛽+𝑥^3𝛼𝑡^3𝛽. The exact solution of (58) is 𝑢(𝑥,𝑡)=𝑥^𝛼𝑡^𝛽.
By employing 2D-FLFs method with 𝑚=3 and 𝑚^′=3, one can have 𝐶^𝑇+−𝐶^𝑇𝐃^𝛼_𝑥+𝑁^𝑇𝐏^𝛽_𝑡−𝐶^𝑇_guess=0. (60)
The numerical results of Example 24 for different values of 𝛼 and 𝛽 are shown in Figure 2. In addition, 𝐿₂ and 𝐿_∞ errors are presented in Table 1. From Table 1, one can conclude that the solutions of 2D-FLFs method are in good agreement with the exact results. Compared with homotopy analysis method (HAM) [36] and homotopy perturbation method (HPM) [37], 2D-FLFs method can get high accuracy solution while only need a few terms of 2D-FLFs.

Table 1 

Errors of Example 24 for different values of 𝛼 and 𝛽 with 𝑀=𝑀^′=4.

(a)

(b)

(c)

(d)

Figure 2 

Numerical results of Example 24 for different values of 𝛼 and 𝛽.

Example 25. Consider the nonlinear time-fractional advection partial differential equation [37–39] 𝐷^𝛽_𝑡𝑢(𝑥,𝑡)+𝑢(𝑥,𝑡)𝑢_𝑥(𝑥,𝑡)=𝑥+𝑥𝑡²,𝑡>0,𝑥∈𝑅,0<𝛽≤1, (61) subject to the initial condition 𝑢(𝑥,0)=0. (62)
Figure 3 gives the approximation solutions of (61) for 𝛽=0.50 with 𝑚=4, 𝑚^′=5 and 𝛽=0.75 with 𝑚=4, 𝑚^′=9. Moreover, Table 2 shows the approximate solutions for (61) obtained for different values of 𝛽 using the fractional variational iteration method (FVIM) [39] and 2D-FLFs method. The values of 𝛽=1 are the only case for which we know the exact solution 𝑢(𝑥,𝑡)=𝑥𝑡. It should be noted that only the fourth-order term of the FVIM was used in evaluating the approximate solutions for Table 2. From Table 2, it clearly appears that 2D-FLFs method is more accurate than FVIM and the obtained results are in good agreement with exact solution.

Table 2 

Numerical values when 𝛽 = 0.50, 0.75, and 1.0 for (61).

(a)

(b)

Figure 3 

Numerical results of Example 25 for different value of 𝛽.

Example 26. We finally consider the linear time-fractional wave equation: 𝜕^2𝛽𝑢𝜕𝑡^2𝛽=12𝑥²𝜕²𝑢𝜕𝑥²,𝑡>0,𝑥∈𝑅,0.5<𝛽≤1, (63) subject to the initial conditions 𝑢()()𝑥,0=𝑥,𝜕𝑢𝑥,0𝜕𝑡=𝑥². (64)
Table 3 gives a comparison of the approximate solutions at different values of 𝛽 using the FVIM [39] and 2D-FLFs method. Figure 4 shows the numerical solutions of 2D-FLFs method for (63) at different values of 𝛽 with 𝑚=3, 𝑚^′=9. The values of 𝛽=1 are the only case for which we know the exact solution 𝑢(𝑥,𝑡)=𝑥+𝑥²sinh(𝑡). As previous, only the fourth-order term of the FVIM was used in evaluating the numerical solutions for Table 3. In the case of 𝛽=1, it can be found that absolute error of 2D-FLFs is not bigger than 1.0𝑒−10 which is very small compared with that obtained by FVIM.

Table 3 

Numerical values when 𝛽 = 0.750, 0.875, and 1.000 for (63).

(a)

(b)

(c)

(d)

Figure 4 

Numerical results of Example 26 for different value of 𝛽.

6. Conclusion

We define a basis of 2D-FLFs and derived its operational matrices of fractional derivative and integration, which are used to approximate the numerical solution of FPDEs. Compared with other numerical methods, 2D-FLFs method can accurately represent properties of fractional calculus. Moreover, only a small number of 2D-FLFs are needed to obtain a satisfactory result. The obtained results demonstrate the validity and applicability of proposed method for solving the FPFEs.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant no. 11272352). The authors are grateful to the anonymous referees for their comments which substantially improved the quality of this paper.