A Least Squares Differential Quadrature Method for a Class of Nonlinear Partial Differential Equations of Fractional Order

: In this paper a new method called the least squares differential quadrature method (LSDQM) is introduced as a straightforward and efﬁcient method to compute analytical approximate polynomial solutions for nonlinear partial differential equations with fractional time derivatives. LSDQM is a combination of the differential quadrature method and the least squares method and in this paper it is employed to ﬁnd approximate solutions for a very general class of nonlinear partial differential equations, wherein the fractional derivatives are described in the Caputo sense. The paper contains a clear, step-by-step presentation of the method and a convergence theorem. In order to emphasize the accuracy of LSDQM we included two test problems previously solved by means of other, well-known methods, and observed that our solutions present not only a smaller error but also a much simpler expression. We also included a problem with no known exact solution and the solutions computed by LSDQM are in good agreement with previous ones.


Introduction
Fractional calculus has its beginnings in the late part of the 17th century, when the concept of fractional derivative was developed in a letter between the mathematicians Leibniz and L'Hospital regarding the derivative of order 1 2 [1]. During the 18th and 19th centuries, many well-known mathematicians, including L.Euler, P.-S. Laplace, J. Fourier, N.H. Abel and J. Liouville contributed to the study of the fractional derivatives [2].
During the 20th century, fractional calculus proved itself to be an essential instrument in the study of many important problems arising in various areas of science and engineering. Particular problems wherein fractional calculus proved to be especially successful include problems in rheology and fluid flow, dynamical processes in self-similar and porous structures, control theory of dynamical systems, frequency-dependent acoustic wave propagation in porous media, modeling of the behavior of viscoelastic materials, micro-grids and decentralized wireless networks, electrochemistry of corrosion, epidemiology (spread of diseases) and many more [1][2][3].
The increased interest observed recently in the study of fractional partial differential equations is motivated by the fact that many engineering and practical phenomena can be modeled better by using fractional equations compared to the corresponding classical differential equations. The study of the fractional equations usually gives better insights into the problem under study. For example, some fractional differential models allow taking into account nonlocal space or time effects, and in the study of properties of materials the models of fractional order are more useful that the ones of integer order since the derivatives of fractional order allow for the modeling of memory and hereditary properties [4][5][6][7].
For several particular problems modeled by fractional order partial differential equations it is possible to find exact solutions. For example, in [3] some implicit analytical solutions for a nonlinear fractional partial differential beam equation are presented, in [8] new optical solutions of the resonant nonlinear Schrödinger equation are constructed in terms of Jacobi elliptic functions, in [9] a generalized bifurcation method is applied in order to find exact solutions of the space-time fractional Drinfel'd-Sokolov-Wilson equation and in [10] the invariant subspace method is generalized and is used to derive exact solutions for a class of nonlinear fractional partial differential equations with mixed partial derivatives.
Unfortunately in the cases of most of the problems modeled by fractional order partial differential equations it may be not possible to find exact solutions for the problems, hence the need to apply numerical or approximate analytical methods. In recent years many classical methods were adapted for the case of fractional derivatives and many more new methods or variants of methods were proposed.
In the following we present a few of the recently employed methods for obtaining numerical and analytical approximate solutions of fractional partial differential equations (in no particular order; some overlapping is present since several methods may be included in more than one category): • Numerical methods include a difference scheme for time fractional heat equations based on the Crank-Nicholson method [11], a high-order predictor-corrector method [12] and an application of Picard's method [13]. • Homotopy-type methods include the well-known homotopy perturbation method [2], the local fractional homotopy perturbation method [14] and the discrete homotopy analysis method [15].

•
Collocation-type methods include the generalized Taylor collocation method [16], a wavelet based collocation method [17], the Legendre collocation method for the nonlinear spacetime fractional partial differential equations [18] and the Sinc collocation method for variable-order fractional integro-partial differential equations [19].

•
Polynomial approximation methods include a Bernstein polynomials method for solving fractional heat and wave-like equations [20] and a Ritz-Galerkin method based on two-dimensional Genocchi polynomials [21]. • Differential transform methods include a method based on the Sumudu transform and the homotopy analysis method [6], the fractional Sumudu decomposition method [22], the fractional reduced differential transform method [23] and the fractional Laplace differential transform method [24]. • Decomposition methods include the Laplace-Adomian decomposition method [25] and the natural transform decomposition method [26].

•
Other methods such as a neural network method based on the sine and the cosine functions [1].
The motivation for this paper was to present a new method obtained from the combination of the differential quadrature method (see [27]) and the least squares method, a method which we named the least squares differential quadrature method (LSDQM). The LSDQM can be used to solve nonlinear partial differential equations of fractional order of the type: where with a, b, c, d real constants, f is a continuous function and D α * denotes the variable order Caputo fractional derivative.
The Equation (1) will be solved together with the initial conditions: where f 0 (x) and f d 0 (x) are continuous functions.

The Least Squares Differential Quadrature Method
For the numerical discretization of the domain D [27]) we will consider a partition consisting of N and M grid points in the x and t directions respectively. In the following we will work under the hypothesis that the functions f , f 0 (x) and f d 0 (x) satisfy the necessary conditions such that the problem (1) and (2) has an unique continuous solution on D.
For the problem (1) and the initial conditions (2) we consider the operator: Ifũ is an approximate solution of Equation (1); the error obtained by replacing the exact solution with the approximationũ is given by the remainder:

Definition 1.
We call an ε-approximate of the problem (1) and (2) an approximate polynomial solutionũ satisfying the following relations in the nodes of the domain's grid: Definition 2. We consider the sequence of polynomials We call the sequence of polynomials P N,M (x, t) convergent to the solution of the problem (1) and (2) if: We consider the sequence of polynomials (7), satisfying the conditions: We will find ε-approximate polynomial solutions of the type: where the constantsc ij are calculated using the following steps: • We attach to the problem (1) and (2) the following real functional: J (c 02 , c 12 , ..., c N2 , c 03 , c 13 , ..., c N3 , ..., c 0M , c 1M , ..., where the constants c 00 , c 10 , ..., c N0 and c 01 , c 11 , ..., c N1 are calculated by solving the linear systems obtained from the initial conditions as follows.

•
From the first initial condition we have: By imposing that the relation is satisfied in the nodes (x 0 , 0), (x 1 , 0), ..., (x N , 0) we obtain a linear system of N + 1 equations in N + 1 unknowns which has the unique solution c 00 , c 10 , ..., c N0 .

•
From the second initial condition we have: Again, by imposing that the relation is satisfied in the corresponding nodes, we obtain a linear system which has the unique solution c 01 , c 11 , ..., c N1 . The following convergence theorem takes place: Proof of Theorem 1. Let u be the exact solution of the problem (1) and (2). From the hypothesis it follows that there exists a sequence of polynomials P N,M (x, t) convergent to u; i.e., We have: It follows that:

Numerical Examples
In order to illustrate the accuracy of our method, we demonstrate herein analytical approximate solutions for two initial value fractional nonlinear problems. These problems were previously solved in [28,29] by using the well-known homotopy perturbation method and the Adomian decomposition method, respectively. We performed the computations by using the Wolfram Mathematica software.

Nonlinear Time-Fractional Advection Partial Differential Equation
This problem was previously solved in [28] by using the homotopy perturbation method and in [29] by using the Adomian decomposition method.

The Case α = 1
For the case α = 1 the exact solution of this problem is u(x, t) = x · t. Since this solution is polynomial, we expect to find the exact solution by using our method, and indeed the exact solution can be found by choosing for the domain D = [0, 1] × [0, 1] a grid step as large as 1. Taking the steps highlighted in the previous section, we seek a polynomial approximate solution of third degree:ũ(x, t) = c 00 + c 10 · x + c 20 · x 2 + c 01 · t + c 02 · t 2 + c 11 · x · t.
Taking into account the initial condition, we find the following values of the constants: Thus the approximate solution has the form: We compute the remainder (4) corresponding to the problem (15) as: We attach to the problem (15) the corresponding functional (11): The minimum of this functional can be easily found in the usual way (find the stationary points, compute the sign of the second differential, etc.). We obtain the corresponding values of the constants as c 01 = 0, c 02 = 0, c 11 = 1 and thus the approximate solution given by LSDQM is actually the exact solution:ũ (x, t) = x · t.
We must mention the fact that homotopy perturbation method and Adomian decomposition method were not able to find the exact solution, only approximate ones with errors ranging from 10 −7 to 10 −1 .

The Case α = 0.75
Since for α < 1 a comparison with an exact solution is not possible, we will compare our solution with the previously computed solutions using other methods in [28,29] by means of the remainders (4) (computed as the expressions obtained by replacing the approximate solutions in the operator (3)) corresponding to each solution.
For the case α = 0.75, using the same steps presented in the previous section, we computed the following approximate solution: u(x, t) = 0.18723135124960114 · t − 0.18069828901844043 · t 2 + 1.0082978876756237 · t · x.  Table 1 present for the case α = 0.75 the absolute values of the remainders (4) corresponding to our approximate solution (green surface), to the solution given by the homotopy perturbation method (HPM) in [28] (yellow surface) and to the solution given by using the Adomian decomposition method (ADM) in [29] (red surface).  For the case α = 0.5 we computed the following approximate solution: u(x, t) = 0.3840538045550427 · t − 0.4015428715206559 · t 2 + 1.0293138444906609 · x · t. Table 2 present for the case α = 0.5 the absolute values of the remainders (4) corresponding to our approximate solution (green surface), to the solution given by the homotopy perturbation method (HPM) in [28] (yellow surface) and to the solution given by using the Adomian decomposition method (ADM) in [29] (red surface).

Nonlinear Time-Fractional Hyperbolic Equation
This problem was also previously solved in [28] by using the homotopy perturbation method and in [29] by using the Adomian decomposition method.
Thus the approximate solution has the form: We compute the remainder (4) corresponding to the problem (16) and attach to the problem the corresponding functional (11) (too long to be included here). Minimizing this functional, we obtain the corresponding values of the constants and thus the following approximate solution: Figure 3 and Table 3 present for the case α = 2 the absolute errors (as the differences in absolute value between the approximate solution and the exact solution) corresponding to our approximate solution (green surface), to the solution given by the homotopy perturbation method (HPM) in [28] (yellow surface) and to the solution given by using the Adomian decomposition method (ADM) in [29] (red surface).   Since for α < 2 a comparison with an exact solution is again not possible, we will compare our solution with the solutions previously computed by the other methods in [28,29] by means of the remainders (4) corresponding to each solution.
For the case α = 1.75, using the same steps presented in the previous sections, we computed the following approximate solution: Table 4 present for the case α = 1.75 the absolute values of the remainders (4) (as expressions obtained by replacing the approximate solutions in the operator (3)) corresponding to our approximate solution (green surface), to the solution given by the homotopy perturbation method (HPM) in [28] (yellow surface) and to the solution given by using the Adomian decomposition method (ADM) in [29] (red surface).  For the case α = 1.5, we computed the following approximate solution:
The exact solution of this problem can not be computed.
In [30] approximate polynomial solutions were computed using the analysis method.
We computed using LSDQM approximate polynomial solutions of for various values of α, solutions which are in good agreement with previous ones.

Conclusions
The present paper introduces the least squares differential quadrature method as a straightforward and efficient method with which to compute analytical approximate polynomial solutions for nonlinear partial differential equations with fractional time derivatives.
The included applications clearly illustrate the accuracy of the method. The comparison with the previously computed approximate solutions by well-known methods, such as the homotopy perturbation method [28] and the Adomian decomposition method [29], revealed not only a smaller error but also a much simpler expression of the solution.