ABSTRACT

Due to the ability to model various complex phenomena where classical calculus failed, fractional calculus is getting enormous attention recently. There are several approaches available for numerical approximations of various types of fractional differential equations. For fractional diffusion equations spectral collocation is one of the efficient and most popular approximation techniques. In this research, we introduce spectral collocation method based on Lagrange’s basis polynomials for numerical approximations of two-dimensional (2D) space fractional diffusion equations where spatial fractional derivative is described in Riemann-Liouville sense. We consider four different types of nodes to generate Lagrange’s basis polynomials and as collocation points in the proposed spectral collocation technique. Spectral collocation method converts the diffusion equation into a system of ordinary differential equations (ODE) for time variable and we use 4^th order Runge-Kutta method to solve the resulting system of ODE. Two examples are considered to verify the efficiency of different types of nodes in the proposed method. We compare approximated solution with exact solution and find that Lagrange’s spectral collocation method gives very high accuracy approximation. Among the four types of nodes, nodes from Jacobi polynomial give highest accuracy and nodes from Chebyshev polynomials of 1^st kind give lowest accuracy in the proposed method.

Keywords:

Lagrange’s Spectral Method, 2D Fractional Diffusion Equation, Collocation Method

1. Introduction

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 [1] . But in 1819, Lacroix first mentioned the derivative of arbitrary order in a text [2] . Later from 1832 Liouville [3] starts to provide the initial foundations of the fractional calculus. A brief history and fundamental theory of fractional calculus can be found in [4] . 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 $1<\alpha <2$ 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 [5] [6] .

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 [7] 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 [8] and Xie et al. [9] . 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 [10] chosen Fibonacci polynomials to express the trial solution in both space and time domain and then used evenly spaced collocation points. Pirim and Ayaz [11] 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 [12] 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. [13] 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 [14] . For collocation purpose he chose Jacobi-Gauss nodes for time domain and for space domain he chosen Gauss-Lobatto nodes. Zayernouri and Karniadakis [15] introduced fractional Lagrange interpolants and developed a new fractional spectral collocation method for 1D fractional differential equations. Krishnaveni et al. [16] 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 [17] used shifted Legendre polynomials and Gauss-Lobatto nodes for 2D space domain. Other numerical methods such as ADI and Finite-Difference methods [18] [19] [20] [21] [22] 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 [23] 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 1^st and 2^nd 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. [24] 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.

2. Preliminaries

In this section, we present a brief introduction of Legendre, Jacobi and Chebyshev polynomials of 1^st and 2^nd kinds. Also short description of Lagrange’s basis polynomials and Caputo’s fractional order derivative are given.

Legendre Polynomials: The nth degree Legendre polynomial $P_n\left(w\right)$ over the domain $\left[-1,1\right]$ are explicitly defined as

$P_n\left(w\right)=\sum _{k=0}^n\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}-n-1\\ k\end{array}\right){\left(\frac{1-w}{2}\right)}^k$ (1)

Chebyshev Polynomials of 1^st Kind: The nth degree Chebyshev polynomials of 1^st kind $T_n\left(w\right)$ over the domain $\left[-1,1\right]$ are explicitly defined as

$T_n\left(w\right)=w^n\sum _{k=0}^{\lfloor \frac{n}{2}\rfloor }\left(\begin{array}{c}n\\ 2k\end{array}\right){\left(1-w^{-2}\right)}^k$ (2)

Chebyshev Polynomials of 2^nd Kind: The nth degree Chebyshev polynomials of 2^nd kind $U_n\left(w\right)$ over the domain $\left[-1,1\right]$ are explicitly defined as

$U_n\left(w\right)=w^n\sum _{k=0}^{\lfloor \frac{n}{2}\rfloor }\left(\begin{array}{c}n+1\\ 2k+1\end{array}\right){\left(1-w^{-2}\right)}^k$ (3)

Jacobi Polynomials: The nth degree Jacobi polynomials $P_n^{i,j}\left(w\right)$ on the interval $\left[-1,1\right]$ are defined as

$P_n^{i,j}\left(w\right)=\frac{\Gamma \left(i+n+1\right)}{i!\Gamma \left(i+j+n+1\right)}\sum _{k=0}^r\left(\begin{array}{c}r\\ k\end{array}\right)\frac{\Gamma \left(i+j+r+k+1\right)}{\Gamma \left(i+k+1\right)}{\left(\frac{w-1}{2}\right)}^k$ (4)

Lagrange Basis Polynomials: For $\left(p+1\right)$ points $x_1,x_2,\cdots ,x_p,x_{p+1}$ Lagrange basis polynomials $L_n\left(x\right);n=1,2,\cdots ,p+1$ each of degree p are defined as follows:

$h\left(x\right)=\prod _{k=1}^{p+1}\left(x-x_k\right)$

$L_n\left(x\right)=\frac{h\left(x\right)}{h^{\prime }\left(x\right)\left(x-x_n\right)};n=1,2,\cdots ,p+1$ (5)

with the property $L_n\left(x_m\right)={\delta }_{mn}$ , where ${\delta }_{mn}$ is the Kronecker delta function. Here $h^{\prime }\left(x\right)$ is the derivative of $h\left(x\right)$ .

Riemann-Liouville Fractional Derivative: Riemann-Liouville fractional derivative operator of order α is denoted by $D^{\alpha }$ and defined by:

$D^{\alpha }f\left(x\right)=\frac{1}{\text{Γ}\left(m-\alpha \right)}\frac{{\text{d}}^m}{\text{d}{x}^m}\underset{0}{\overset{x}{\int }}{\left(x-t\right)}^{m-\alpha -1}f\left(t\right)\text{d}t,\alpha >0$ (6)

with $m-1<\alpha <m,m\in \mathbb{N}$ .

Then for $\alpha >0,m-1<\alpha <m,m\in \mathbb{N}$ and $p\ge 0$

$D^{\alpha }{x}^p=\frac{\text{Γ}\left(p+1\right)}{\text{Γ}\left(p-\alpha +1\right)}{x}^{p-\alpha }$ (7)

where $\mathbb{N}=\left\{1,2,3,\cdots \right\}$ .

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

$\begin{array}{l}\frac{\partial u\left(x,y,t\right)}{\partial t}=g_1\left(x,y\right)\frac{{\partial }^{\alpha }u\left(x,y,t\right)}{\partial x^{\alpha }}+g_2\left(x,y\right)\frac{{\partial }^{\beta }u\left(x,y,t\right)}{\partial y^{\beta }}\\ +f\left(x,y,t\right);0\le x\le a,0\le y\le b,t>0\end{array}$ (8)

subject to the initial condition

$u\left(x,y,0\right)=g_3\left(x,y\right);0\le x\le a,0\le y\le b$ (9)

and non-homogeneous Dirichlet’s boundary conditions

$u\left(0,y,t\right)=g_4\left(y,t\right),u\left(a,y,t\right)=g_5\left(y,t\right);0\le y\le b,t\ge 0$ (10)

$u\left(x,0,t\right)=g_6\left(x,t\right),u\left(x,b,t\right)=g_7\left(x,t\right);0\le x\le a,t\ge 0$ (11)

Here $g_1\left(x,y\right)$ and $g_2\left(x,y\right)$ are the diffusion coefficients, $f\left(x,y,t\right)$ is the force function and $u\left(x,y,t\right)$ is the unknown function to be solved.

For numerical approximations of Equation (8) by Lagrange’s spectral collocation method we first divide the interval $\left[0,a\right]$ of space domain x into p subintervals and interval $\left[0,b\right]$ of y domain into q subinterval by means of following nodes:

$0=x_1,x_2,x_3,\cdots ,x_p,x_{p+1}=a$ (12)

$0=y_1,y_2,y_3,\cdots ,y_p,y_{p+1}=b$ (13)

These nodes can be chosen with no specific pattern and from anywhere in the domain. Then using these nodes ${\left\{x_i\right\}}_{i=1}^{i=p+1}$ from x-domain we form Lagrange’s basis polynomials $L_n\left(x\right);n=1,2,\cdots ,p,\left(p+1\right)$ each of degree p and using nodes

${\left\{y_j\right\}}_{j=1}^{j=q+1}$ from y-domain we form Lagrange’s basis polynomials

${L^{\prime }}_m\left(y\right);m=1,2,\cdots ,q,\left(q+1\right)$ each of degree q. Now at these nods Lagrange’s basis polynomials satisfies $L_n\left(x_i\right)={\delta }_{ni}$ and ${L^{\prime }}_m\left(y_j\right)={\delta }_{mj}$ . These nodes are also used as collocation points in the later part of this section.

Now as spectral approximation of $u\left(x,y,t\right)$ , the solution of Equation (8) we consider trial solution $\tilde{u}\left(x,y,t\right)$ in the form

$\tilde{u}\left(x,y,t\right)=\sum _{n=1}^{p+1}\sum _{m=1}^{q+1}{}_m{}^{n}v\left(t\right)L_n\left(x\right){L^{\prime }}_m\left(y\right)$ (14)

$\begin{array}{l}={\displaystyle \sum _{m=1}^{q+1}\left[{}_m{}^{1}v\left(t\right)L_1\left(x\right)+{}_m{}^{p+1}v\left(t\right)L_{p+1}\left(x\right)\right]{L^{\prime }}_m\left(y\right)}+{\displaystyle \sum _{n=2}^p{\displaystyle \sum _{m=2}^q{}_m{}^{n}v\left(t\right)L_n\left(x\right){L^{\prime }}_m\left(y\right)}}\\ +{\displaystyle \sum _{n=2}^p\left[{}_1{}^{n}v\left(t\right){L^{\prime }}_1\left(y\right)+{}_{q+1}{}^{n}v\left(t\right){L^{\prime }}_{q+1}\left(y\right)\right]L_n\left(x\right)}\end{array}$ (15)

Here Lagrange’s basis polynomials $L_n\left(x\right)$ and ${L^{\prime }}_m\left(y\right)$ are used as the basis

functions in spectral approximations and ${\left\{{}_m{}^{n}v\left(t\right)\right\}}_{m=1,2,\cdots ,q,q+1}^{n=1,2,\cdots ,p,p+1}$ are the expansion

coefficients needed to be determined. Now we write the residual $R\left(x,y,t\right)$ for Equation (8) as

$\begin{array}{c}R\left(x,y,t\right)=\frac{\partial \tilde{u}\left(x,y,t\right)}{\partial t}-g_1\left(x,y\right)\frac{{\partial }^{\alpha }\tilde{u}\left(x,y,t\right)}{\partial x^{\alpha }}\\ -g_2\left(x,y\right)\frac{{\partial }^{\beta }\tilde{u}\left(x,y,t\right)}{\partial y^{\beta }}-f\left(x,y,t\right)\end{array}$ (16)

where

$\begin{array}{c}\frac{{\partial }^{\alpha }\tilde{u}}{\partial x^{\alpha }}=\sum _{m=1}^{q+1}\left[{}_m{}^{1}v\left(t\right)L_1^{\alpha }\left(x\right)+{}_m{}^{p+1}v\left(t\right)L_{p+1}^{\alpha }\left(x\right)\right]{L^{\prime }}_m\left(y\right)\\ +\sum _{n=2}^p\sum _{m=2}^q{}_m{}^{n}v\left(t\right)L_n^{\alpha }\left(x\right){L^{\prime }}_m\left(y\right)\\ +\sum _{n=2}^p\left[{}_1{}^{n}v\left(t\right){L^{\prime }}_1\left(y\right)+{}_{q+1}{}^{n}v\left(t\right){L^{\prime }}_{q+1}\left(y\right)\right]L_n^{\alpha }\left(x\right)\end{array}$ (17)

$\begin{array}{c}\frac{{\partial }^{\beta }\tilde{u}}{\partial y^{\beta }}=\sum _{m=1}^{q+1}\left[{}_m{}^{1}v\left(t\right)L_1\left(x\right)+{}_m{}^{p+1}v\left(t\right)L_{p+1}\left(x\right)\right]{L^{\prime }}_m^{\beta }\left(y\right)\\ +\sum _{n=2}^p\sum _{m=2}^q{}_m{}^{n}v\left(t\right)L_n\left(x\right){L^{\prime }}_m^{\beta }\left(y\right)\\ +\sum _{n=2}^p\left[{}_1{}^{n}v\left(t\right){L^{\prime }}_1^{\beta }\left(y\right)+{}_{q+1}{}^{n}v\left(t\right){L^{\prime }}_{q+1}^{\beta }\left(y\right)\right]L_n\left(x\right)\end{array}$ (18)

$\frac{\partial \tilde{u}}{\partial t}=\sum _{n=1}^{p+1}\sum _{m=1}^{q+1}{L}_n\left(x\right){L^{\prime }}_m\left(y\right)\frac{\text{d}{}_m{}^{n}v\left(t\right)}{\text{d}t}$ (19)

Here $L_n^{\alpha }\left(x\right)$ is the $\alpha$ order Riemann-Liouville fractional derivative of $L_n\left(x\right)$ with respect to x and ${L^{\prime }}_m^{\beta }\left(y\right)$ is the $\beta$ order Riemann-Liouville fractional derivative of ${L^{\prime }}_m\left(x\right)$ 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

$\left(p-1\right)\times \left(q-1\right)$ interior points ${\left\{\left(x_i,y_j\right)\right\}}_{j=2,3,\cdots ,q-1,q}^{i=2,3,\cdots ,p-1,p}$ 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

$\begin{array}{l}R\left(x_i,y_j,t\right)=\frac{\partial \tilde{u}\left(x_i,y_j,t\right)}{\partial t}-g_1\left(x_i,y_j\right)\frac{{\partial }^{\alpha }\tilde{u}\left(x_i,y_j,t\right)}{\partial x^{\alpha }}\\ -g_2\left(x_i,y_j\right)\frac{{\partial }^{\beta }\tilde{u}\left(x_i,y_j,t\right)}{\partial y^{\beta }}-f\left(x_i,y_j,t\right)=0\end{array}$ (20)

for $i=2,3,\cdots ,p-1,p$ and $j=2,3,\cdots ,q-1,q$ with

$\frac{{\partial }^{\alpha }\tilde{u}\left(x_i,y_j,t\right)}{\partial x^{\alpha }}=\sum _{m=1}^{q+1}{k}_1\left(m,x_i,t\right){L^{\prime }}_m\left(y_j\right)+\sum _{n=2}^p\sum _{m=2}^q{}_m{}^{n}v\left(t\right)L_n^{\alpha }\left(x_i\right){L^{\prime }}_m\left(y_j\right)$ (21)

$\frac{{\partial }^{\beta }\tilde{u}\left(x_i,y_j,t\right)}{\partial y^{\beta }}=\sum _{n=2}^p\sum _{m=2}^q{}_m{}^{n}v\left(t\right)L_n\left(x_i\right){L^{\prime }}_m^{\beta }\left(y_j\right)+\sum _{n=2}^p{k}_2\left(n,y_j,t\right)L_n\left(x_i\right)$ (22)

$\frac{\partial \tilde{u}\left(x_i,y_j,t\right)}{\partial t}=\sum _{n=2}^p\sum _{m=2}^q{L}_n\left(x_i\right){L^{\prime }}_m\left(y_j\right)\frac{\text{d}{}_m{}^{n}v\left(t\right)}{\text{d}t}$ (23)

where

$k_1\left(m,x_i,t\right)=\left[{}_m{}^{1}v\left(t\right)L_1^{\alpha }\left(x_i\right)+{}_m{}^{p+1}v\left(t\right)L_{p+1}^{\alpha }\left(x_i\right)\right]$ (24)

$k_2\left(n,y_j,t\right)=\left[{}_1{}^{n}v\left(t\right){L^{\prime }}_1^{\beta }\left(y_j\right)+{}_{q+1}{}^{n}v\left(t\right){L^{\prime }}_{q+1}^{\beta }\left(y_j\right)\right]$ (25)

Thus Equation (20) gives us the following system of ODE for time variable t

$\frac{\text{d}{}_j{}^{i}v\left(t\right)}{\text{d}t}=a_{i,j}\left(t\right)+g_1\left(x_i,y_j\right)\sum _{n=2}^p{}_j{}^{n}v\left(t\right)L_n^{\alpha }\left(x_i\right)+g_2\left(x_i,y_j\right)\sum _{m=2}^q{}_m{}^{i}v\left(t\right){L^{\prime }}_m^{\beta }\left(y_j\right)$ (26)

with

$a_{i,j}\left(t\right)=g_1\left(x_i,y_j\right)\cdot k_1\left(j,x_i,t\right)+g_2\left(x_i,y_j\right)\cdot k_2\left(i,y_j,t\right)+f\left(x_i,y_j,t\right)$ (27)

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

$\sum _{n=2}^p\sum _{m=2}^q{}_m{}^{n}v\left(0\right)L_n\left(x_i\right){L^{\prime }}_m\left(y_j\right)=g_3\left(x_i,y_j\right)$ (28)

$\Rightarrow {}_j{}^{i}v\left(0\right)=g_3\left(x_i,y_j\right)$ (29)

Now forcing the trial solution to agree with the boundary conditions for boundary lines $x=x_1$ and $x=x_{p+1}$ , that is, to satisfy the boundary conditions at Equation (10) gives

$\sum _{m=1}^{q+1}{}_m{}^{1}v\left(t\right){L^{\prime }}_m\left(y\right)=g_4\left(y,t\right)$ (30)

$\sum _{m=1}^{q+1}{}_m{}^{p+1}v\left(t\right){L^{\prime }}_m\left(y\right)=g_5\left(y,t\right)$ (31)

Then at the points $y_j;j=1,2,\cdots ,q,q+1$ Equation (30) and (31) becomes

${}_m{}^{1}v\left(t\right)=g_4\left(y_m,t\right)$ and ${}_m{}^{p+1}v\left(t\right)=g_5\left(y_m,t\right)$ for $m=1,2,\cdots ,q,q+1$(32)

Again forcing the trial solution to agree with the boundary conditions, now for boundary lines $y=y_1$ and $y=y_{q+1}$ , that is, to satisfy the boundary conditions at Equation (11) gives

$\sum _{n=2}^p{}_1{}^{n}v\left(t\right)L_n\left(x\right)=g_6\left(x,t\right)$ (33)

$\sum _{n=2}^p{}_{q+1}{}^{n}v\left(t\right)L_n\left(x\right)=g_7\left(x,t\right)$ (34)

Now at the points $x_i;i=2,3,\cdots ,p$ Equation (33) and (34) becomes

${}_1{}^{n}v\left(t\right)=g_6\left(x_n,t\right)$ and ${}_{q+1}{}^{n}v\left(t\right)=g_7\left(x_n,t\right)$ for $n=2,3,\cdots ,\left(p-1\right),p$ (35)

System of ODE in Equation (26) with its initial condition Equation (29) can

be solved for the value of expansion coefficients ${\left\{{}_m{}^{n}v\left(t\right)\right\}}_{m=2,3,\cdots ,q}^{n=2,3,\cdots ,p}$ at any $t=T_F$ .

Rest of the values of expansion coefficients can be obtained from Equation (32) and (35). Replacing the values of expansion coefficients ${\left\{{}_m{}^{n}v\left(t\right)\right\}}_{m=2,3,\cdots ,q,q+1}^{n=1,2,3,\cdots ,p,p+1}$ for any particular $t=T_F$ into trial solution in Equation (14) will gives us the approximate solution $\tilde{u}\left(x,y,T_F\right)$ of Equation (8) in terms of Lagrange’s basis polynomials $L_n\left(x\right)$ and ${L^{\prime }}_m\left(y\right)$ .

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 ${\left\{{}_m{}^{n}v\left(t\right)\right\}}_{m=2,3,\cdots ,q}^{n=2,3,\cdots ,p}$ at

$T_F=1$ by 4^th order Runge-Kutta method with time step size 0.01.

In both examples spatial domain was $\left[0,1\right]\times \left[0,1\right]$ and we consider $p=q=5$ . 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 $\left(-1,1\right)$ , to generate nodes ${\left[x_k\right]}_{k=1}^{k=p+1}$ over any arbitrary interval $\left[c,d\right]$ , first we choose $x_1=c$ and $x_{p+1}=d$ . Then, the remaining $\left(p-1\right)$ nodes are generated by shifting the roots of one of the above $\left(p-1\right)\text{th}$ degree polynomials into the interval $\left[c,d\right]$ .

Since in both examples 2D space domain is square shaped and we consider $p=q=5$ , following four sets of nodes are thus used for both x and y domain.

$\text{chebyshev}{\text{1}}^{\text{st}}\text{kind}=\left\{0,0.0380602,0.308658,0.691342,0.96194,1\right\}$

$\text{chebyshev}{\text{2}}^{\text{nd}}\text{kind}=\left\{0,0.0954915,0.345492,0.654508,0.904508,1\right\}$

$\text{legendre}=\left\{0,0.0694318,0.330009,0.669991,0.930568,1\right\}$

$\text{jacobi}=\left\{0,0.1174723,0.357384,0.642616,0.882528,1\right\}$

To compare the approximations given by the proposed method with the exact solution we consider absolute local error for any given $\left(x,y,t\right)$ given by

$Er\left(x,y,t\right)=\left|u\left(x,y,t\right)-\tilde{u}\left(x,y,t\right)\right|$

We also consider maximum absolute error given by

$M_{AE}\left(t\right)=\max \left\{Er\left(x,y,t\right):0\le x\le a,0\le y\le b\right\}$

For numerical simulation we develop codes for the proposed method using MATLAB.

Example 1: Consider the diffusion equation [17] [18]

$\begin{array}{l}\frac{\partial u\left(x,y,t\right)}{\partial t}=g_1\left(x,y\right)\frac{{\partial }^{1.8}u\left(x,y,t\right)}{\partial x^{1.8}}+g_2\left(x,y\right)\frac{{\partial }^{1.6}u\left(x,y,t\right)}{\partial y^{1.6}}\\ +f\left(x,y,t\right);0\le x\le 1,0\le y\le 1,t>0\end{array}$

with boundary conditions

$u\left(0,y,t\right)=0,u\left(1,y,t\right)={\text{e}}^{-t}{y}^{3.6};0\le y\le 1,t\ge 0$

$u\left(x,0,t\right)=0,u\left(x,1,t\right)={\text{e}}^{-t}{x}^3;0\le x\le 1,t\ge 0$

and with initial condition

$u\left(x,y,0\right)=x^3{y}^{3.6};0\le x\le 1,0\le y\le 1$

where

$g_1\left(x,y\right)=\frac{\Gamma \left(2.2\right)}{6}{x}^{2.8}y,g_2\left(x,y\right)=\frac{2}{\text{Γ}\left(4.6\right)}x{y}^{2.6},f\left(x,y,t\right)=-\left(1+2xy\right){\text{e}}^{-t}{x}^3{y}^{3.6}$

The exact solution of this problem is given by $u\left(x,y,t\right)={\text{e}}^{-t}{x}^3{y}^{3.6}$ .

Applying the proposed Lagrange’s spectral collocation technique for $p=q=5$ into Example 1 with four different types of nodes resulting surface graph of absolute local error $Er\left(x,y,1\right)$ are given into Figures 1-4. Then Figure 5 contains the absolute local error curve $Er\left(x,0.7,1\right)$ and Figure 6 contains the absolute local error curve $Er\left(0.8,y,1\right)$ for example 1.

Example 2: Consider the diffusion equation [17]

$\begin{array}{l}\frac{\partial u\left(x,y,t\right)}{\partial t}=g_1\left(x,y\right)\frac{{\partial }^{1.5}u\left(x,y,t\right)}{\partial x^{1.5}}+g_2\left(x,y\right)\frac{{\partial }^{1.5}u\left(x,y,t\right)}{\partial y^{1.5}}\\ +f\left(x,y,t\right);0\le x\le 1,0\le y\le 1,t>0\end{array}$

with boundary conditions

$u\left(0,y,t\right)=0,u\left(1,y,t\right)={\text{e}}^{-t}{y}^3;0\le y\le 1,t\ge 0$

$u\left(x,0,t\right)=0,u\left(x,1,t\right)={\text{e}}^{-t}{x}^2;0\le x\le 1,t\ge 0$

and with initial condition

$u\left(x,y,0\right)=x^2{y}^3;0\le x\le 1,0\le y\le 1$

where

$g_1\left(x,y\right)=\frac{\Gamma \left(1.5\right)}{2}\left(3-2x\right),g_2\left(x,y\right)=\frac{\Gamma \left(2.5\right)}{6}\left(4-y\right),$

$f\left(x,y,t\right)={\text{e}}^{-t}\left(x^2\left(-y^{\frac{3}{2}}+y-4\right)y^{\frac{3}{2}}+\sqrt{x}\left(2x-3\right)y^3\right)$

The exact solution of this problem is given by $u\left(x,y,t\right)={\text{e}}^{-t}{x}^2{y}^3$ .

Now by proposed Lagrange’s spectral collocation technique with $p=q=5$ , surface graph of absolute local error $Er\left(x,y,1\right)$ by four different types of nodes for Example 2 are given into Figures 7-10. Then Figure 11 contains the absolute

local error curve $Er\left(x,0.7,1\right)$ and Figure 12 contains the absolute local error curve $Er\left(0.8,y,1\right)$ for Example 2.

Table 1 contains maximum absolute error $M_{AE}\left(1\right)$ 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 2^nd, Legendre and Chebyshev 1^st 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 1^st nodes decreased rapidly and for Legendre nodes error decreased slowly whereas errors by both Jacobi and Chebyshev 2^nd nodes keeps increasing. But in first three quarter of the

x-domain performance of Jacobi nodes, Chebyshev 2^nd, Legendre and Chebyshev 1^st nodes are found in decreased order.

5. Conclusion

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, 1^st and 2^nd 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 1^st 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.

Cite this paper

Molla, H.U. and Nova, M.H. (2018) Lagrange’s Spectral Collocation Method for Numerical Approximations of Two-Dimensional Space Fractional Diffusion Equation. American Journal of Computational Mathematics, 8, 121-136. https://doi.org/10.4236/ajcm.2018.82010

References