THE ALTERNATIVE LEGENDRE TAU METHOD FOR SOLVING NONLINEAR MULTI-ORDER FRACTIONAL DIFFERENTIAL EQUATIONS

In this paper, the alternative Legendre polynomials (ALPs) are used to approximate the solution of a class of nonlinear multi-order fractional differential equations (FDEs). First, the operational matrix of fractional integration of an arbitrary order and the product operational matrix are derived for ALPs. These matrices together with the spectral Tau method are then utilized to reduce the solution of the mentioned equations into the one of solving a system of nonlinear algebraic equations with unknown ALP coefficients of the exact solution. The fractional derivatives are considered in the Caputo sense and the fractional integration is described in the Riemann-Liouville sense. Numerical examples illustrate that the present method is very effective for linear and nonlinear multi-order FDEs and high accuracy solutions can be obtained only using a small number of ALPs.


Introduction
Many problems in physics, chemistry and engineering give rise to the linear and nonlinear fractional differential equations [1,11,16,18,21,34,35]. Due to the lack of analytical solutions to fractional differential equations in most cases, various numerical methods developed by many authors for obtaining approximate solutions of these equations. These methods include spectral methods [2,5,9,13,22,25,29,37], Laplace transform method [19], differential transform method [12,36], Adomian decomposition method [23], variational iteration method [24,26], wavelet method [17,[28][29][30][31]33,38] and etc. Among them, the spectral methods based on the application of different sets of basis functions have special attention of many researchers. The main feature of these methods lies in their accuracy for a given number of unknowns. It is proved that, for smooth problems, they offer exponential rates of convergence [5,39].
In this paper, we intend to use the operational Tau method based on ALPs, a new family of orthogonal polynomials were recently introduced by Vladimir Chelyshkov [6], to numerically solve the following nonlinear multi-order FDE: where γ i (i = 1, . . . , p) and δ i (i = 1, . . . , q) are real constant coefficients, m = ⌈ν⌉, and 0 < β 1 < β 2 < · · · < β p < ν. In Eq. (1.1), D ν u(x) denotes the Caputo fractional derivative of order ν for u(x) and the values of d i (i = 0, . . . , m − 1) describe the initial state of u(x) and f (x) is a given source function.
In the present method, by applying multiple fractional integration in the Riemann-Liouville sense, the FDE (1.1) is first converted to a fully integrated form. Then, the various signals involved in the integrated form equation are approximated by representing them as linear combinations of ALPs. Finally, the integrated form equation is converted to an algebraic equation by introducing the operational matrix of the product and the operational matrix of fractional integration for ALPs.
The reminder of this paper is organized as follows. First, in Section 2, we give necessary knowledge about the fractional calculus theory. In Section 3, we briefly describe some characteristics of the ALPs and derive the ALP operational matrix of the product and the ALP operational matrix of fractional integration in the Riemann-Liouville sense. Section 4 is devoted to the application of alternative Legendre Tau method using ALPs for the solution of problem (1.1)-(1.2). Some numerical examples are given in Section 5 to demonstrate the efficiency and accuracy of the method. Conclusions of the work are given in Section 6.

Preliminaries of the fractional calculus
In this section, we give some definitions and basic properties of fractional calculus theory that will be used in this paper [8,27].
is called the Caputo fractional differential operator of order ν.
The Riemann-Liouville fractional integral operator J ν and the Caputo fractional differential operator D ν satisfy the following properties

Definition and properties
Let n be a fixed integer number. The set P n = {P nk (x)} n k=0 of alternative Legendre polynomials is defined over interval [0, 1] as [6] These polynomials are orthogonal over the interval [0, 1] with respect to the weight function w(x) = 1, and satisfy the orthogonality relationships where δ kl denotes the Kronecker delta. In contrast to common sets of orthogonal polynomials, every member in P n has degree n. Relation (3.1) yields Rodrigues's type representation It follows from (3.4) that By using Eq. (3.1), the ALP vector (3.5) can be written in the form and Q is the following upper triangular matrix of order (n + 1) × (n + 1) Using the orthogonal property (3.3), the dual matrix of Φ(x) is obtained as

Suppose that
Since H n is a finite dimensional subspace of H, then it is closed [20, and for every given f ∈ H there exist a unique best approximationf n ∈ H n [20, Moreover, we have [7,Theorem 4.14] Therefore, any function f (x) ∈ H may be approximated by its projection in the subspace H n as and In the following theorem, we obtain an estimation of the error norm of the best approximation of smooth functions by some polynomials in H n .
where C is a constant such that Proof. Assume that p n (x) is the interpolating polynomial to f at points x l , where x l , l = 0, 1, . . . , n, are the roots of the (n + 1)-degree shifted Chebyshev polynomial in Ω. Then, for any x ∈ [0, 1], we have [15] f By taking into account the estimates for Chebyshev interpolation nodes [15], we obtain Sincef is the unique best approximation of f in H n , then we have Taking the square root of both sides of (3.13) gives the result.

Operational matrices of fractional integration and product for ALPs
Many papers have employed the operational matrix for spectral methods [2,4,5,14, 17, 22, 25, 28-31, 33, 38]. In this section, we will derive the Riemann-Liouville operational matrix of fractional integration and also the operational matrix of product for ALPs. These operational matrices will be applied in the next section for solving nonlinear FDEs by the alternative Legendre spectral Tau method. We need the following two lemmas-the proofs are quite easy and so we omit them.
Lemma 3.1. Let P nk (x) and P nj (x) are kth and jth ALPs, respectively and v 1 = min{j, k} and v 2 = max{j, k}. Then the product of P nk (x) and P nj (x) is a polynomial of degree at most 2n that can be written as

Lemma 3.2. Let r be a nonnegative integer, then
By combining Lemmas 3.1 and 3.2, we can easily obtain the following result.

Lemma 3.3.
Let P ni (x), P nj (x) and P nk (x) are ith, jth and kth ALPs, respectively.
If, according to Lemma 3.1, we define P nk (x)P nj (x) = 2n ∑ r=k+j λ (k,j) r x r , then In the next theorem, the ALP operational matrix of fractional integration in the Riemann-Liouville sense is first derived.

kr ] n k,r=0 is the ALP operational matrix of integration of fractional order ν in the Riemann-Liouville sense in which
In (3.14), the sign ≃ means that any element in the right-hand side vector is an approximation to the corresponding one in the left-hand side vector.

Proof. Since, by relation (2.4), the Riemann-Liouville's fractional integration is a linear operation, then using (2.1) and (3.1) we have
Using Eq. (3.10), one can approximate x ν+j in terms of ALPs as where, by Eq. (3.9) and Lemma 3.2, we have ε jr = (2r + 1) Therefore, Varying k from 0 to n gives the result.
The following theorem is of great importance, when we want to reduce the integral equation (1.1) to a set of algebraic equations in the next section.

k=0 is the ALP operational matrix of the product in whicĥ
The values τ ijk in (3.16) are computed by means of Lemma 3.3. The sign ≃ is previously defind in Theorem 3.2.

Alternative Legendre spectral Tau method
In this section, we use the spectral Tau method based on ALPs to solve nonlinear FDE (1.1) with initial conditions (1.2). To this end, we first apply the Riemann-Liouville fractional integral of order ν on (1.1). Then, making use of (2.1)-(2.4), we get the integrated form of equation (1.1) as where in which m i = ⌈β i ⌉, i = 1, . . . , p.
In order to use the alternative Legendre Tau method for solving the fully integrated problem (4.1) with initial conditions (4.2), we approximate u(x) and g(x) by the ALPs as where the vector G = [g 0 , . . . , g n ] T is known but U = [u 0 , . . . , u n ] T is an unknown vector. We can also use (3.15) to approximate the nonlinear term u i (x), for any integer i ≥ 1, in terms of ALPs as whereÛ is the (n + 1) × (n + 1) ALP operational matrix of product for U . Now, using (2.4), (4.5) and Theorems 3.2 and 3.3, the Riemann-Liouville fractional integrals of orders ν − β i and ν can be written, respectively, as and where P (ν) is the ALP operational matrix of fractional integration of order ν in the Riemann-Liouville sense which is derived in Theorem 3.2. Applying Eqs. (4.3)-(4.7), the residual function R n (x) for Eq. (4.1) can be written as As in a typical Tau method, see [5], we generate n − m + 1 nonlinear algebraic equations by employing Also from Eqs. (1.2) and (3.6), we get where Q is the upper triangular matrix defined in (3.7) and e i+1 is the i-th standard unit vector of order n + 1.
Eqs. (4.9) and (4.10) form a system of n + 1 nonlinear algebraic equations. These nonlinear equations can be solved for unknown coefficients of the vector U . Consequently,ū n (x) given in Eq.

Illustrative examples
In this section, the presented method is applied to five examples to show its accuracy and applicability. Experiments were performed on a personal computer using a 2.50 GHz processor and the codes were written in Mathematica 9. In the case of nonlinear equations, we have solved system (4.9)-(4.10) using the Mathematica function FindRoot, which uses Newton's method as the default method. In all the cases, this function has succeeded to obtain an accurate approximate solution of the system, even starting with a zero initial approximation.
Example 5.1. Consider the following nonlinear FDE The exact solution of this problem is u(x) = x. If we use the technique described in Section 4 with n = 3, we get ] .
Thus, by relation (4.4), we have which is the exact solution.
Example 5.2. [10,19,32] Consider the following nonlinear initial value problem, The exact solution of this problem is u(x) = x 2 . If we use the technique described in Section 4 with n = 4, we get 9 7 ] .
Thus, by relation (4.4), we have which is the exact solution. The exact solution of this problem is u(x) = x+1. If we use the technique described in Section 4 with n = 2, we get 35 12 ] .
Thus, by relation (4.4), we have which is the exact solution. Also, for n = 3 we obtain ] .
Therefore, by relation (4.4), we get which is again the exact solution. The exact solution of this problem is u(x) = x 3 . If we use the technique described in Section 4 with n = 4, we get U T = [ 0, 0, 1 8 , 9 8 ] .