On a computational method for non-integer order partial differential equations in two dimensions

This manuscript is concerning to investigate numerical solutions for different classes including parabolic, elliptic and hyperbolic partial differential equations of arbitrary order (PDEs). The proposed technique depends on some operational matrices of fractional order differentiation and integration. To compute the mentioned operational matrices, we apply shifted Jacobi polynomials in two dimension. Thank to these matrices, we convert the PDE under consideration to an algebraic equation which is can be easily solved for unknown coefficient matrix required for the numerical solution. The proposed method is very efficient and need no discretization of the data for the proposed PDE. The approximate solution obtain via this method is highly accurate and the computation is easy. The proposed method is supported by solving various examples from well known articles. 2010 Mathematics Subject Classifications: 35R11, 65M06, 65N30


Introduction
Many phenomenons and problems in biology, physics, chemistry, dynamics, fluid mechanics, optical technology,etc and in engineering disciplines are modeled by partial differential equations. The mentioned equations arise in numbers of scientific models like propagation of water waves, long waves and chemical reaction diffusion equations, for detail see [7,29,36,40,41,41]. Recently, considerable attention has been given to study partial differential equations with non-integer order rather than classical order. Because in many situations, fractional order partial differential equations (FPDEs) modeled real worlds phenomenons more accurately as compared to classical order partial differential equations. It has been proved that many interesting process of electromagnetic, acoustics, viscoelasticity, electrochemistry and material sciences are well described by using fractional order partial differential equations, see [6,14,17,35]. The present work is related to study the following linear FPDE given as + C ∂ α w(x, y) ∂y α + D ∂ β w(x, y) ∂x β + E ∂ β w(x, y) ∂y β + F w(x, y) = g(x, y), α ∈ (1, 2], β ∈ (0, 1], corresponding to the initial conditions w(x, y)| x=0 = φ(y), w x (x, y)| x=0 = ψ(y). (1) Where A, B, C, D, E, F are real constants and w(x, y) is the solution of the system to be computed and also w(x, y), g(x, y), ∈ C([0, δ] × [0, δ]). Under the classical order α = 2, β = 1, the class (1) of FPDEs represents different classes of integer order PDE such as parabolic, ecliptic and hyperbolic partial differential equations. The concerned PDE is parabolic corresponding to integer order if B 2 − 4AC = 0, ecliptic if B 2 − 4AC < 0 and hyperbolic, if B 2 − 4AC > 0, where each category describes different phenomena in various disciplines of science and engineering. The aforesaid classes of partial differential equations have many applications in various field of science and technology. As Laplace and Poisson equations are ecliptic partial differential equations which have many applications in electromagnetic theory, and other branches of physics, mechanics, etc. Similarly, Heat equation is parabolic and wave equation is hyperbolic partial differential equations. For the importance and uses of these equations, we refe [1,15,30,38]. The area which has greatly attracted the attention of researches is devoted to the existence theory and numerical solution of fractional order partial differential equations. Since the differential equations involved fractional order are not easy to solve for theirs exact solutions in many situations, therefore a strong motivation has been found to find the approximate solutions.
In general there does not exist a method which produce exact solutions of FPDEs. In many situations for linear partial differential equations it is quite difficult task to find exact solutions. Therefore, numerous techniques were developed to find approximate solutions to the aforementioned equations. In [8,18], authors applied Homotopy analysis method(HAM) and He's variation iteration method(VIM) to calculate approximate solutions of nonlinear FPDEs. In same line, authors, used Adomain decomposition method(ADM) [16], homotopy perturbation method(HPM) [21], Fourier transform method(FTM) [49], Laplace and Natural transform methods [42,43] to solve fractional order partial differential equation. Also, some authors used fractional differential transform (FDTM) [2,3,37], to find numerical solutions of FPDEs. But all these method have their own advantages and disadvantages in application point of view, as homotopy methods depend on small parameters which restricted these methods. Similarly the methods that are involving integral transform are also limited in applications.
In last few decades some interesting numerical schemes based on radial basis functions(RBFs) and meshless techniques depend on collocation and (RBFs) for solving (PDEs), see [31,47].
Recently, numerical schemes based on operational matrices have attracted the attention of many researches. The mentioned techniques provide highly accurate numerical solutions to both linear and nonlinear ordinary as well as partial differential equations of classical and fractional order. In the concerned schemes some operational matrices of fractional order integration and differentiation are constructed which play central rules in approximation of solutions. In [32,50], authors used Haar wavelets and Legendre wavelets to solve linear FPDEs. With the same line some operational matrices based on Chebyshev wavelets were also used to find numerical solutions to FPDEs, see [28]. In [4,5,9,26,33,48,51], authors constructed operational matrices based on Lagendre, Bernstein and Jacobi, Chebyshev polynomials for numerical solutions of fractional differential equations. The concerned method used Tau-collocation method to form the mentioned matrices. As in these methods discretization of date is required which need extra memory but provided best approximate solutions to (FPDEs). To overcome this difficulty, in [22,23] authors constructed the operational matrices from the afore mentioned polynomial directly without discretizing the polynomials and obtained excellent solutions to (FPDEs). The spectral methods based on operational matrices are the best tools to find approximate solutions for both ordinary and partial fractional differential equations, for detail see [10,12,20,25,27,45,46]. Motivated by the aforesaid work, in this paper we have considered a general multi-terms fractional order partial differential equations which represents different classes of partial differential equations by assigning various values to the constants involved in (1) to satisfied the specific conditions. The operational matrices involved in this paper are based on shifted Jacobi polynomials. With the the help of these matrices of fractional order integration and differentiation, the considered equation is converted to an algebraic equation which is easily soluble for unknown coefficient matrix needed in the approximate solution.
All the computations are done by using Matlab. For the demonstration, we solve different examples of physical interest and also the comparison of our method with other method is provided.

Preliminaries
Some auxiliary results and notions needed throughout in this paper are provided as [14]: Definition 2.1. The Riemann-Liouville integral of arbitrary order α ∈ of a function w(x) is defined by such that the integral on right hand side converges pointwise on [0, ∞). Also the aforesaid integral satisfies the following relations (2) I α I β w(x) = I α+β w(x); where the integral on the right is converges on [0, ∞). Further, the operator D satisfies the following properties in particular for any constant c, D α c = 0 and The following relations are necessary throughout this paper.
For any function w(x), the following results hold.

The shifted Jacobi polynomials
The famous Jacobi polynomials with two parameters say p, q are defined on [0, δ] as (2) The orthogonality condition is where W is the weight function, Thus any any function w(x) which is square integrable in [0, δ] may be approximated interm of shifted Jacobi polynomials as as n → ∞ the approximation is converging to the exact function. Inview of (3) (4),(5), we can easily find the constant D l . Writing (6) in matrix form as such that N = n + 1, K N is the coefficient matrix andΦ N (x) is matrix of vectors function. Extending the notation from uni-dimension to two dimension space in order to define two dimensional shifted Jacobi polynomials of order N over the region The orthogonality condition for where D lm can be computed as while W Writing the notation D n = D ab , then (10) implies where
Theorem 3.2. The fractional order differentiation of Φ N 2 (x, y) as defined in (14) corresponding to 'y as D α y (Φ N 2 (x, y)) H (α,y) where H (α,y) N 2 ×N 2 is the operational matrix given as H (α,y) and v = N i + j + 1, r = N a + b + 1, ℵ v,r,k = Ω i,j,a,b,k for i, j, a, b = 0, 1, 2, . . . , n, In the same line, the fractional order differentiation of Φ N 2 (x, y) corresponding to 'x as provided as where H (α,x) Theorem 3.3. The fractional order derivative of Φ N 2 (x, y) is constructed in (14) corresponding to x and y as where H (α,x,y) N 2 ×N 2 is the operational matrix of given as H (α,x,y) and r = N i + j + 1, v = N a + b + 1, ℵ v,r,k = Ω i,j,a,b,k for i, j, a, b = 0, 1, 2, . . . , n, Proof. The operational matrix H (α,x,y) N 2 ×N 2 can be easily obtain by using the product of two operational matrices, H

Numerical procedure based on the operational matrices for FPDE (1)
With the help of operational matrices constructed in section 3, the considered class of FPDEs (1) is converted to simple algebraic equations which is easily soluble. Assume that the following hold ∂ α w(x, y) Then inview of Definition 2.1, we have which on using initial conditions of (1) yields that Now inview of the relation obtained in (32), we have Therefore, with the help of (30), (31), (32) and (33), our considered class of FPDEs (1) becomes which on simplification yields On writing we, have from (35) AK It is obvious that (36) is an algebraic equation of matrices. Solving the (36), for K N 2 , plugging its value in (32), we can get the numerical solution of the problem under consideration.

Numerical Examples
This section is concerning to test the above techniques by some well known problems of physical interest given bellow. α L ∞ at N = 6 L ∞ at N = 8 L ∞ at N = 10 L ∞ at N = 12 1 At α = 2 the exact solution of (37) is w(x, y) = cos(2x) sin(2y). Here A = 1, B = 0, C = 1, D = E = 0, F = 8, clearly B 2 − 4AC < 0, the given FPDE is ecliptic. We plot the approximate solution corresponding to the exact solution at α = 2. We see from the Figure  1, that the proposed scheme provides close agreement to that of exact solution at very small scale level. Further, we compared the absolute error L ∞ = w−w N ∞ at various fractional order corresponding to different scale level N in the following Table 1. From the Table 1, we see that as the fractional order approaches to its integer values the absolute error decreasing and vice versa and on the other hand enlarging the scale level, the accuracy also increases even for fractional order.
We plot the approximate solution corresponding to the exact solution at α = 2. We see from the Figure 2, that the approximate solution obtained via adapting procedure gives close agreement to that of exact solution at very small scale level. We also compute the absolute error L ∞ = w − w N ∞ for different fractional order using various values of N in the given Table 2. We observe from the Table 2 that as α approaches to its integer value, the absolute error decreasing. Similarly, on the other hand increasing the scale level, the accuracy also increases for using different fractional order α. Example 3. Consider the following fractional partial differential equation  At α = 2, β = 1 the exact solution of (39) is w(x, y) = (xy) 4 − (xy) 3 . Here A = 5, B = 6, C = −4, D = 7, E = 8, F = 9. As B 2 − 4AC > 0, the given FPDE is hyperbolic. In the Figure ??, we have plotted the approximate solution corresponding to the exact solution at α = 2. We see that the approximate solution obtained via proposed method has the close agreement to that of exact solution at sufficiently small scale level. Further, we compute the absolute error L ∞ = w − w N ∞ for different fractional order and using various values of N in the given Table 3. We observe from the Table 3 that as α approaches to its integer value, the absolute error decreasing. Similarly, by increasing the scale level, the accuracy also increases for taking different fractional order.  Remark 1. The numerical solutions of FPDEs by using shifted Jacobi polynomials are also depending on the values of the parameters of the considered polynomials. Therefore the absolute error also depends upon on the values of the Jacobi polynomials parameters. We give the following  6. Comparison of the proposed method with some already existing Haar wavelet method The Shifted Legendre polynomials are the special case of Shifted Jacobi polynomials and they can be obtained by p = q = 0, δ = 1 in the relation (2). We compared our technique with the existing methods available in [39], where the authors constructed the procedure by discretizing the data and using the Haar wavelet method, which need extra memory to obtained the operational matrices of fractional derivative and integration because data need discretization and collocation. While in our prosed method, we have obtained the operational matrices without discretization of data and using any collocation method. In the following examples we provide the comparison of the proposed method with the existing some methods.
Comparison between our proposed method and that given in [   Similarly, we can also bring out comparison for other problems of FDEs and FPDEs. From the Table 5, it is clear that the proposed method provides powerful tools to obtain approximate solutions which is closely related with the exact solution of the problem.

Conclusion
In this study, we have successfully applied the operational matrices methods of Shifted Jacobi polynomials to various classes of FPDEs in two dimensions. Thank to these matrices, the concerned FPDE is converted to a simple algebraic equations which is further solved for coefficient matrix needed in the approximate solutions of the considered equation. The method provides an excellent agreement to the exact solutions a very small scale level. The important point in the method proposed in this study need no extra parameter like perturbation method neither required discretization of the data which need extra memory and waste of time and suffers from computational complexity some times. All these mentioned things are avoided in this method and a direct computational procedure has been adopted. In future this method, we can easily extended to nonlinear FPDEs.