Theory of Fractional calculus is almost of same age as that of the classical calculus. The concept of a fractional derivative initiates the theory of fractional calculus when this concept first appeared in a letter of L’Hospital to Leibnitz in 1695  . But in 1819, Lacroix first mentioned the derivative of arbitrary order in a text  . Later from 1832 Liouville  starts to provide the initial foundations of the fractional calculus. A brief history and fundamental theory of fractional calculus can be found in  . Though many important physical processes in almost all branches of sciences can be described using various tools of classical calculus, there are also many complex phenomena that cannot be modeled using classical calculus. Being one of the most important generalizations of classical calculus, fractional calculus has been used to describe many such complex processes recently. In last few decades, scientists and engineers have found fractional order derivative very powerful to describe problems in anomalous diffusion, signal processing, control processing, fractional stochastic system and viscoelasticity. In space fractional diffusion equation, spatial derivative of order replaces the classical second order spatial derivative in classical diffusion equation and results into super diffusive flow model. Unlike classical order derivative, fractional derivative has several kinds of definitions such as Caputo, Riemann-Liouville, Grünwald-Letnikov and Riesz. Detailed discussion about various definitions of fractional derivative can be found in   .
Spectral method is very effective and efficient and hence popular among scientists and engineers for numerical approximations of both classical and fractional differential and integral equations. In spectral method unknown solution is first approximated by the trial solution as a finite sum of a set of known basis functions. Then expansion coefficients in trial solution can be determined by several techniques like collocation, Galerkin, Petrov Galerkin, tau etc. In all of these techniques expansion coefficients are determined by minimizing the difference between exact and trial solutions. For spectral approximations of fractional differential equations collocation technique is the most efficient and widely used technique. In collocation technique expansion coefficients are determined by forcing the trial solution to satisfy the differential equation at some suitably chosen points from the domain known as collocation points. For a p parameter solution p collocation points are required within the domain and which results into p residual equations. In recent times, researchers showed immense interests to try out various types of polynomials as basis sets in spectral approximations and different kinds of collocation points. For numerical approximations of 1D space fractional diffusion equations Khader  used Chebyshev polynomials of space variable as the basis set in spectral approximations and roots of shifted Chebyshev polynomials as the collocation points to determine the expansion coefficients. Chebyshev polynomials are also used for both time and space domain as basis set in spectral approximations of 1D fractional diffusion equation by Azizi and Loghmani  and Xie et al.  . But Azizi and Loghmani applied collocation technique with Gauss-Lobatto nodes whereas Xie et al. used tau method to determine expansion coefficients. For numerical approximations of 1D fractional diffusion equation Bahsi and Yalcinbas  chosen Fibonacci polynomials to express the trial solution in both space and time domain and then used evenly spaced collocation points. Pirim and Ayaz  also chosen evenly spaced collocation points but expressed the trial solution in terms of Hermite polynomials for approximations of fractional order system of differential equations. Huang and Zheng  used Jacobi polynomials and Lagrange’s basis polynomials for spectral approximations of 1D space fractional diffusion equations and used Gauss-Lobatto nodes for collocation. Doha et al.  presented spectral tau method for 1D space fractional diffusion equation where they expand the solution by Jacobi polynomials. For spectral expansion of trial solution for time fractional diffusion equations Lagrange interpolation polynomials are used in both 1D space and time domain by Huang  . For collocation purpose he chose Jacobi-Gauss nodes for time domain and for space domain he chosen Gauss-Lobatto nodes. Zayernouri and Karniadakis  introduced fractional Lagrange interpolants and developed a new fractional spectral collocation method for 1D fractional differential equations. Krishnaveni et al.  proposed a hybrid method for 1D space fractional diffusion equation where trial solution is approximated by spectral expansion of fractional shifted Legendre polynomials. To approximate the solution of 2D space fractional diffusion equation by spectral collocation method Bhrawy  used shifted Legendre polynomials and Gauss-Lobatto nodes for 2D space domain. Other numerical methods such as ADI and Finite-Difference methods      are also very efficient for 2D problems.
In this research, we consider two dimensional space fractional diffusion equation with an initial condition and non-homogeneous Dirichlet’s boundary conditions where spatial fractional derivatives are described in Riemann-Liouville sense. Variety of spectral collocation methods are available for 1D fractional diffusion equations but spectral collocation method haven’t been tried with lots of variations in 2D problems. To the best of our knowledge Lagrange’s basis polynomials haven’t been used in spectral collocation method for 2D space fractional diffusion equations. Hence we use Lagrange’s basis polynomials for spectral expansion of trial solution in 2D space domain. To determine the expansion coefficients collocation technique is used which transforms the diffusion equation into a first order system of ODE of time variable. System of ODE is then solved by fourth order Runge-Kutta  method. Lagrange’s polynomials depends on the selected nodes. Also in collocation method choice of collocation points affect the resultant accuracy. So we consider four different types of nodes in Lagrange’s spectral collocation method to generate Lagrange’s basis polynomials and as collocation points to investigate their performances in the proposed method. Four different sets of nodes are derived from shifted roots of Chebyshev polynomials of 1st and 2nd kinds, Legendre polynomials and Jacobi polynomials. The performances in terms of resultant accuracy by four sets of nodes in proposed Lagrange’s spectral collocation method are illustrated and compared through numerical examples. Similar investigation on 1D space fractional diffusion equation is carried out by Nova et al.  recently.
Remaining parts of this paper are organized as follows: Section 2 contains the preliminaries of different polynomials and fractional derivative. Then in Section 3, detailed formulations of Lagrange’s spectral method for 2D space fractional diffusion equation are given. Next in Section 4, numerical examples are considered to demonstrate the accuracy of the proposed method and to compare the performance of four types of nodes. Finally, conclusion is drawn about this research in Section 5.
In this section, we present a brief introduction of Legendre, Jacobi and Chebyshev polynomials of 1st and 2nd kinds. Also short description of Lagrange’s basis polynomials and Caputo’s fractional order derivative are given.
Legendre Polynomials: The nth degree Legendre polynomial over the domain are explicitly defined as
Chebyshev Polynomials of 1st Kind: The nth degree Chebyshev polynomials of 1st kind over the domain are explicitly defined as
Chebyshev Polynomials of 2nd Kind: The nth degree Chebyshev polynomials of 2nd kind over the domain are explicitly defined as
Jacobi Polynomials: The nth degree Jacobi polynomials on the interval are defined as
Lagrange Basis Polynomials: For points Lagrange basis polynomials each of degree p are defined as follows:
with the property , where is the Kronecker delta function. Here is the derivative of .
Riemann-Liouville Fractional Derivative: Riemann-Liouville fractional derivative operator of order α is denoted by and defined by:
Then for and
Like classical integer order derivative, Riemann-Liouville fractional order derivative is also a linear operator.
3. Lagrange’s Spectral Collocation Method
Lagrange’s Spectral Collocation method for numerical approximations of 2D space fractional diffusion equation will be presented in this section. In this research, we consider following form of diffusion equation
subject to the initial condition
and non-homogeneous Dirichlet’s boundary conditions
Here and are the diffusion coefficients, is the force function and is the unknown function to be solved.
For numerical approximations of Equation (8) by Lagrange’s spectral collocation method we first divide the interval of space domain x into p subintervals and interval of y domain into q subinterval by means of following nodes:
These nodes can be chosen with no specific pattern and from anywhere in the domain. Then using these nodes from x-domain we form Lagrange’s basis polynomials each of degree p and using nodes
from y-domain we form Lagrange’s basis polynomials
each of degree q. Now at these nods Lagrange’s basis polynomials satisfies and . These nodes are also used as collocation points in the later part of this section.
Now as spectral approximation of , the solution of Equation (8) we consider trial solution in the form
Here Lagrange’s basis polynomials and are used as the basis
functions in spectral approximations and are the expansion
coefficients needed to be determined. Now we write the residual for Equation (8) as
Here is the order Riemann-Liouville fractional derivative of with respect to x and is the order Riemann-Liouville fractional derivative of with respect to y.
To determine the expansion coefficients we chose collocation technique and for this we force the residual equation at Equation (16) to be zero at
interior points of the spatial domain and
force the trial solution to satisfy the boundary conditions Equation (10) & (11) at the boundary lines.
First by setting the residual at Equation (16) to be zero at the interior points of the spatial domain, we have the equation
for and with
Thus Equation (20) gives us the following system of ODE for time variable t
Now forcing the trial solution to satisfy the initial condition Equation (9) at the interior points gives us the initial condition for system of ODE in Equation (26) as
Now forcing the trial solution to agree with the boundary conditions for boundary lines and , that is, to satisfy the boundary conditions at Equation (10) gives
Then at the points Equation (30) and (31) becomes
and for (32)
Again forcing the trial solution to agree with the boundary conditions, now for boundary lines and , that is, to satisfy the boundary conditions at Equation (11) gives
Now at the points Equation (33) and (34) becomes
and for (35)
System of ODE in Equation (26) with its initial condition Equation (29) can
be solved for the value of expansion coefficients at any .
Rest of the values of expansion coefficients can be obtained from Equation (32) and (35). Replacing the values of expansion coefficients for any particular into trial solution in Equation (14) will gives us the approximate solution of Equation (8) in terms of Lagrange’s basis polynomials and .
4. Numerical Results
In this section, we consider two examples to apply the proposed method with four different types of nodes to investigate the accuracy of the method and nodes. System of ordinary differential equations in Equation (26) with initial condition
Equation (29) is solved for the value of expansion coefficients at
by 4th order Runge-Kutta method with time step size 0.01.
In both examples spatial domain was and we consider . Four different types of nodes are thus derived as follows:
Since all the roots of Legendre, Jacobi and both kinds of Chebyshev polynomials are in the interval , to generate nodes over any arbitrary interval , first we choose and . Then, the remaining nodes are generated by shifting the roots of one of the above degree polynomials into the interval .
Since in both examples 2D space domain is square shaped and we consider , following four sets of nodes are thus used for both x and y domain.
To compare the approximations given by the proposed method with the exact solution we consider absolute local error for any given given by
We also consider maximum absolute error given by
For numerical simulation we develop codes for the proposed method using MATLAB.
Example 1: Consider the diffusion equation  
with boundary conditions
and with initial condition
The exact solution of this problem is given by .
Applying the proposed Lagrange’s spectral collocation technique for into Example 1 with four different types of nodes resulting surface graph of absolute local error are given into Figures 1-4. Then Figure 5 contains the absolute local error curve and Figure 6 contains the absolute local error curve for example 1.
Example 2: Consider the diffusion equation 
with boundary conditions
and with initial condition
Figure 1. Error graph for Chebyshev 1st kind.
Figure 2. Error graph for Chebyshev 2nd kind.
Figure 3. Error graph for Legendre.
Figure 4. Error graph for Jacobi.
Figure 5. Absolute local error curve for .
Figure 6. Absolute local error curve for .
The exact solution of this problem is given by .
Now by proposed Lagrange’s spectral collocation technique with
, surface graph of absolute local error
by four different types of nodes for Example 2
Figure 7. Error graph for Chebyshev 1st kind.
Figure 8. Error graph for Chebyshev 2nd kind.
Figure 9. Error graph for Legendre.
Figure 10. Error graph for Jacobi.
Figure 11. Absolute local error curve for .
local error curve and Figure 12 contains the absolute local error curve for Example 2.
Table 1 contains maximum absolute error by four different sets of nodes for both examples.
Observing the surface graphs and error curves for absolute local error in both examples it is evident that proposed Lagrange’s collocation technique gives very high accuracy approximations of 2D space fractional diffusion equations for all four types of collocation points. The maximum absolute errors in both examples for Jacobi nodes, Chebyshev 2nd, Legendre and Chebyshev 1st nodes are found in increasing order. The same performance of these nodes is clearly visible in Figure 5, Figure 11 and Figure 12. From Figure 5 we see that in the last quarter of the x-domain absolute local error for Chebyshev 1st nodes decreased rapidly and for Legendre nodes error decreased slowly whereas errors by both Jacobi and Chebyshev 2nd nodes keeps increasing. But in first three quarter of the
Figure 12. Absolute local error curve for .
Table 1. Maximum absolute error.
x-domain performance of Jacobi nodes, Chebyshev 2nd, Legendre and Chebyshev 1st nodes are found in decreased order.
Spectral collocation method with Lagrange’s basis polynomials for approximations of 2D space fractional diffusion equations is introduced in this research. We used four different types of nodes generated from roots of Legendre, Jacobi, 1st and 2nd kinds of Chebyshev polynomials. Nodes are used to form Lagrange’s basis polynomials for spectral approximations and then as collocation points to determine the expansion coefficients. Performances of these four types of nodes in Lagrange’s spectral collocation method are investigated through two examples of 2D space fractional diffusion equation. Four types of nodes in proposed method give very high accuracy approximations. We found that among four types of nodes Chebyshev 1st kind gives lowest accuracy and Jacobi nodes give highest accuracy. We compared approximate solutions with exact solution and considered absolute local error and maximum absolute error. Four different types of nodes showed consistent performance in both examples.
 Liouville, J. (1832) Mémoire sur quelques Quéstions de Géometrie et de Mécanique, et sur un nouveau genre de Calcul pour résoudre ces Quéstions. Journal de l'école Polytechnique, tome XIII, XXIe cahier, 1-69.
 Ross, B. (1975) A Brief History and Exposition of the Fundamental Theory of fractional Calculus. In: Ross, B., Ed., Fractional Calculus and Its Applications, Lecture Notes in Mathematics, Springer, Berlin, Heidelberg, 57, 1-36.
 Khader, M.M. (2011) On the Numerical Solutions for the Fractional Diffusion Equation. Communications in Nonlinear Science and Numerical Simulation, 16, 2535-2542.
 Azizi, H. and Loghmani, G.B. (2013) Numerical Approximation for Space Fractional Diffusion Equation Chebyshev Finite Difference Method. Journal of Fractional Calculus and Applications, 4, 303-311.
 Xie, J., Huang, Q. and Yang, X. (2016) Numerical Solution of the One-Dimensional Fractional Convection Diffusion Equations Based on Chebyshev Operational Matrix. SpringerPlus, 5, 1149.
 Bahsi, A.K. and Yalcinbas, S. (2016) Numerical Solution and Error Estimations for the Space Fractional Diffusion Equation with Variable Coefficients via Fibonacci Collocation Method. SpringerPlus, 5, 1375.
 Doha, E.H., et al. (2014) The Operational Matrix Formulation of the Jacobi Tau Approximation for Space Fractional Diffusion Equation. Advances in Difference Equations, 2014, 231.
 Huang, F. (2012) A Time-Space Collocation Spectral Approximation for a Class of Time Fractional Differential Equations. International Journal of Differential Equations, 2012, Article ID: 495202.
 Bhrawy, A.H. (2014) A New Legendre Collocation Method for Solving a Two-Dimensional Fractional Diffusion Equation. Abstract and Applied Analysis, 2014, Article ID: 636191.
 Tadjeran, C. and Meerschaert, M.M. (2007) A Second-Order Accurate Numerical Method for the Two-Dimensional Fractional Diffusion Equation. Journal of Computational Physics, 220, 813-823.
 Liu, F., et al. (2013) Numerical Simulation for Two-Dimensional Riesz Space Fractional Diffusion Equations with Nonlinear Reaction Term. Central European Journal of Physics, 11, 1221-1232.
 Nasrollahzadeh, F. and Hosseini, S.M. (2016) An Implicit Difference ADI Method for Two-Dimensional Space-Time Fractional Diffusion Equation. Iranian Journal of Mathematical Sciences and Informatics, 11, 71-86.
 Sweilam, N.H., Nagy, A.M. and Almajbri, T.F. (2014) Error Analysis of an Explicit Finite Difference Approximation for the Two-Dimensional Space Fractional Diffusion Equation. Journal of Fractional Calculus and Applications, 5, 1-10.
 Nova, M.H., Molla, H.U. and Banu, S. (2017) Comparison of Numerical Approximations of One-Dimensional Space Fractional Diffusion Equation Using Different Types of Collocation Points in Spectral Method Based on Lagrange’s Basis Polynomials. American Journal of Computational Mathematics, 7, 469-480.