Numerical Solution for Initial and Boundary Value Problems of Fractional Order

A. M. Kawala

Department of Mathematics, Faculty of Science, Helwan University, Cairo, Egypt.

DOI: 10.4236/apm.2018.812051                

Abstract

Fractional calculus has been used in many fields, such as engineering, population, medicine, fluid mechanics and different fields of chemistry and physics. These fields were found to be best described using fractional differential equations (FDEs) to model their processes and equations. One of the well-known methods for solving fractional differential equations is the Shifted Legendre operational matrix (LOM) method. In this article, I proposed a numerical method based on Shifted Legendre polynomials for solving a class of fractional differential equations. A fractional order operational matrix of Legendre polynomials is also derived where the fractional derivatives are described by the Caputo derivative sense. By using the operational matrix, the initial and boundary equations are transformed into the products of several matrixes and by scattering the coefficients and the products of matrixes. I got a system of linear equations. Results obtained by using the proposed method (LOM) presented here show that the numerical method is very effective and appropriate for solving initial and boundary value problems of fractional ordinary differential equations. Moreover, some numerical examples are provided and the comparison is presented between the obtained results and those analytical results achieved that have proved the method’s validity.

Keywords

Fractional Differential Equations, Shifted Legendre Polynomials, Operational Matrix, Caputo Derivative

1. Introduction

Fractional calculus, the theory of differentiation and integration to non-integer order, is very useful for the description of various physical phenomena, such as damping laws, diffusion process, etc. Fractional differential equations extend prominent tools that are perfectly representing many engineering and physical problems. However, solving fractional differential equations are a challenging and stimulating area of research in mathematics and engineering since these fractional equations do not have exact and analytic solutions and that is the main reason made accurate numerical techniques preferable for solving fractional differential equations.

Many applications of shifted Legendre polynomials have been exemplified in research [1] [2] [3] [4]. In this article, an application of Legendre polynomials to solve fractional differential equations is provided. The main aim is to generalize Legendre operational matrix to fractional calculus and based on using the operational matrix of an orthogonal function for solving differential equations. In the last decade or so, extensive research has been done in the numerical development methods that are numerically stable for both linear and nonlinear fractional differential equations [5] - [13]. Orthogonal functions and polynomial series have received appreciable attention in dealing with various problems of dynamic systems. Legendre polynomials are well known family of orthogonal polynomials on the interval [−1, 1] that have many applications [14]. They are widely used because of their good properties in the approximation of functions.

The proposed method is to obtain the numerical solution for FDE of the form [15] presented in Equation (I)-(III) based on the shifted Legendre method:

$A_1{D}^{2}y\left(x\right)+A_2{D}^{1}y\left(x\right)+A_3{D}^{\alpha }y\left(x\right)+A_4{D}^{\beta }y\left(x\right)+A_{5}y\left(x\right)=f\left(x\right)$ (I)

subject to the conditions

$y\left(0\right)=\delta ,y{\left(0\right)}^{\prime }={\delta }_1$ (II)

or the boundary conditions

$y\left(0\right)={\gamma }_0,y\left(R\right)={\gamma }_1$ (III)

where $A_1,A_2,A_3,A_4,A_5,{\delta }_0,{\delta }_1,{\gamma }_0,{\gamma }_1$ are constants.

$m-1<\alpha ,\beta <m$ and f(x) are the source terms.

The fractional derivatives are defined in the Caputo sense. The main idea in the current work is to apply the shifted Legendre polynomials and the operational matrix of fractional derivative together to discretize Equation (1) to get a satisfactory result. Figures 1-3 and results in next sections show the effectiveness of the proposed method in comparison with the analytical results.

The remainder of the article is organized as follows. In the next section, some mathematical preliminaries of the fractional calculus theory are introduced in addition to some relevant properties of the Legendre polynomials. Section 3 summarizes the application of the shifted Legendre method to the solution of problems (I)-(III). As a result, a system of algebraic equations is obtained and the solution of the considered problem is given. Section 4 illustrates applying the Legendre operational matrix of fractional derivative for solving multi-order fractional differential equation. In Section 5, the proposed method is applied to several examples. Finally, conclusion is given in Section 6. The numerical results are all computed using Mathematica Development Environment.

Figure 1. Comparison of y(x) Analytical and Approximate for m = 3.

Figure 2. Comparison of y(x) Analytical and Approximate for m = 5.

Figure 3. Comparison of y(x) Analytical and Approximate for m = 3.

2. LOM Preliminaries and Notations

In this section, some basic definitions and properties of fractional calculus theory are given that are further used in this article.

2.1. The Fractional Derivative in the Caputo Sense

The fractional calculus is a name for the theory of integrals and derivatives of arbitrary order, which unifies and generalizes the notions of integer-order differentiation and n-fold integration [16] [17]. There are various definitions of fractional integration and differentiation, such as Grunwald_Letnikov’s definition and Riemann_Liouville’s definition. In the present work, the fractional derivatives are considered in the Caputo sense. The reason for adopting the Caputo definition, as pointed by [18] , is as follows: to solve differential equations (both classical and fractional), I need to specify additional conditions in order to produce a unique solution. For the case of the Caputo fractional differential equations, these additional conditions are just the traditional conditions, which are akin to those of classical differential equations, and are therefore familiar to us. In contrast, for the Riemann_Liouville fractional differential equations, these additional conditions constitute certain fractional derivatives (and/or integrals) of the unknown solution at the initial point x = 0, which are functions of x. For more details see [19].

Definition 2.1. The Caputo definition of the fractional-order derivative is defined as

$D^{\alpha }f\left(x\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_0^x\frac{f^{\left(n\right)}\left(t\right)}{{\left(x-t\right)}^{\alpha +1-n}}\text{d}t},n-1<\alpha \le n,n\in N,$ (1)

where, $\alpha >0$ is the order of the derivative and n is the smallest integer greater than $\alpha$ For the Caputo derivative I have [20]

$D^{\alpha }C=0$ (C is a constant), (2)

$D^{\alpha }{x}^{\beta }=\{\begin{array}{l}0,\text{for}\beta \in N_0\text{and}\beta <\alpha \\ \frac{\Gamma \left(\beta +1\right)}{\Gamma \left(\beta +1-\alpha \right)}{x}^{\beta -\alpha },\text{for}\beta \in N_0\text{and}\beta \ge \alpha \text{or}\beta \notin N\text{and}\beta >\alpha \end{array}$ (3)

I use the ceiling function, $\lceil \alpha \rceil$ to denote the smallest integer greater than or equal to $\alpha$, and the floor function $\lfloor \alpha \rfloor$ to denote the largest integer less than or equal to $\alpha$ . Also $N=\left\{1,2,\cdots \right\}$ and $N_0=\left\{0,1,2,\cdots \right\}$ . Recall that for $\alpha \in N$, the Caputo differential operator coincides with the usual differential operator of an integer order.

2.2. Properties of Shifted Legendre Polynomials

The well-known Legendre polynomials are defined on the interval $\left[-1,1\right]$ and can be determined with the aid of the following recurrence formulae:

$L_{i+1}\left(z\right)=\frac{2i+1}{i+1}z{L}_i\left(z\right)-\frac{i}{i+1}{L}_{i-1}\left(z\right),i=1,2,\cdots$

where $L_0\left(z\right)=1$ and $L_1\left(z\right)=z$ . In order to use these polynomials on the interval $x\in \left[0,1\right]$ I define the so-called shifted Legendre polynomials by introducing the change of variable.

$z=2x-1$ . Let the shifted Legendre polynomials $L_i\left(2x-1\right)$ be denoted by $P_i\left(x\right)$ . Then $P_i\left(x\right)$ can be obtained as follows:

$P_{i+1}\left(x\right)=\frac{\left(2i+1\right)\left(2x-1\right)}{i+1}{P}_i\left(x\right)-\frac{i}{i+1}{P}_{i-1}\left(x\right),i=1,2,\cdots ,$ (5)

where $P_0\left(x\right)=1$ and $P_1\left(x\right)=2x-1$ . The analytic form of the shifted Legendre polynomial $P_i\left(x\right)$ of degree i given by

$P_i\left(x\right)=\underset{k=0}{\overset{i}{{\displaystyle \sum }}}\frac{{\left(-1\right)}^{i+k}\left(i+k\right)!}{\left(i-k\right)!{\left(k!\right)}^2}{x}^k$ (6)

Note that $P_i\left(0\right)={\left(-1\right)}^i$ and $P_i\left(1\right)=1$ The orthogonality condition is

${\displaystyle {\int }_0^1{p}_i\left(x\right)p_j\left(x\right)\text{d}x}=\{\begin{array}{l}\frac{1}{2i+1}\text{for}i=j\\ 0\text{for}i\ne j\end{array}$ (7)

A function $y\left(x\right)$, square integrable in $\left[0,1\right]$, may be expressed in terms of shifted Legendre polynomials as

$y\left(x\right)={\displaystyle {\sum }_{j=0}^{\infty }{c}_j{P}_j\left(x\right)}$,

where $c_j$ the coefficients are given by $P_j(x)$

$c_j=\left(2j+1\right){\displaystyle {\int }_0^{1}y\left(x\right)p_j\left(x\right)\text{d}x},j=1,2,\cdots$

In practice, only the first $\left(m+1\right)$ - terms shifted Legendre polynomials are considered. Then I have

$y\left(x\right)=\underset{j=0}{\overset{m}{{\displaystyle \sum }}}{c}_j{P}_j\left(x\right)=C^{\text{T}}\varphi \left(x\right),$

where the shifted Legendre coefficient vector C and the shifted Legendre vector $\varphi \left(x\right)$ are given by

$C^{\text{T}}=\left[c_0,\cdots ,c_m\right]$,

$\varphi \left(x\right)={\left[P_0\left(x\right),P_1\left(x\right),\cdots ,P_m\left(x\right)\right]}^{\text{T}}$ (8)

The derivative of the vector $\varphi \left(x\right)$ can be expressed by

$\frac{\text{d}\varphi \left(x\right)}{\text{d}x}=D^{\left(1\right)}\varphi \left(x\right)$, (9)

where $D^{\left(1\right)}$ is the $\left(m+1\right)\times \left(m+1\right)$ operational matrix of derivative given by

$D^{(1)}=\left(d_{ij}\right)=\{\begin{array}{l}2\left(2j+1\right),\text{for}j=i-k,\{\begin{array}{l}k=1,3,\cdots ,m,\text{if}m\text{odd}\\ k=1,3,\cdots ,m-1,\text{if}m\text{even}\end{array}\\ 0,\text{otherwise},\end{array}$

for example for even m I had

$D^{\left(1\right)}=2\left(\begin{array}{cccccccc}0& 0& 0& 0& \cdots & 0& 0& 0\\ 1& 0& 0& 0& \cdots & 0& 0& 0\\ 0& 3& 0& 0& \cdots & 0& 0& 0\\ 1& 0& 5& 0& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 1& 0& 5& 0& \cdots & 2m-3& 0& 0\\ 0& 3& 0& 7& \cdots & 0& 2m-1& 0\end{array}\right)$

3. Generalized Legendre Operational Matrix to Fractional Calculus

By using Equation (9), it is clear that

$\frac{{\text{d}}^n\varphi \left(x\right)}{\text{d}{x}^n}={\left(D^{\left(1\right)}\right)}^n\varphi \left(x\right)$, (10)

where $n\in N$ and the superscript, in $D^{\left(1\right)}$, denotes matrix powers. Thus

$D^{\left(n\right)}={\left(D^{\left(1\right)}\right)}^n,n=1,2,\cdots$ (11)

Lemma 1. Let $P_i\left(x\right)$ be a shifted Legendre polynomial then

$D^{\left(\alpha \right)}{P}_i\left(x\right)=0,i=0,1,\cdots ,\lceil \alpha \rceil -1,\alpha >0$ (12)

Proof. Using Equations (2), (3) in Equation (6) the lemma can be proved.

In the following theorem I generalize the operational matrix of derivative of shifted Legendre polynomials given in (9) for fractional derivative.

Theorem 1. Let $\varphi \left(x\right)$ be shifted Legendre vector defined in (8) and also suppose, $\alpha >0$ then

$D^{\alpha }\varphi \left(x\right)\simeq D^{\left(\alpha \right)}\varphi \left(x\right)$ (13)

where $D^{\left(\alpha \right)}$ is the $\left(m+1\right)\times \left(m+1\right)$ operational matrix of fractional derivative of order $\alpha$ in the Caputo sense and is defined as follows:

$D^{\left(\alpha \right)}=\left(\begin{array}{cccc}0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0\\ {\displaystyle {\sum }_{k=\lceil \alpha \rceil }^{\lceil \alpha \rceil }{\theta }_{\lceil \alpha \rceil ,0,k}}& {\displaystyle {\sum }_{k=\lceil \alpha \rceil }^{\lceil \alpha \rceil }{\theta }_{\lceil \alpha \rceil ,1,k}}& \cdots & {\displaystyle {\sum }_{k=\lceil \alpha \rceil }^{\lceil \alpha \rceil }{\theta }_{\lceil \alpha \rceil ,m,k}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle {\sum }_{k=\lceil \alpha \rceil }^i{\theta }_{i,0,k}}& {\displaystyle {\sum }_{k=\lceil \alpha \rceil }^i{\theta }_{i,1,k}}& \cdots & {\displaystyle {\sum }_{k=\lceil \alpha \rceil }^i{\theta }_{i,m,k}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle {\sum }_{k=\lceil \alpha \rceil }^m{\theta }_{m,0,k}}& {\displaystyle {\sum }_{k=\lceil \alpha \rceil }^m{\theta }_{m,1,k}}& \cdots & {\displaystyle {\sum }_{k=\lceil \alpha \rceil }^m{\theta }_{m,m,k}}\end{array}\right)$ (14)

where ${\theta }_{i,j,k}$ is given by

${\theta }_{i,j,k}=\left(2j+1\right){\displaystyle {\sum }_{i=0}^j\frac{{\left(-1\right)}^{i+j+k+l}\left(i+k\right)!\left(l+j\right)!}{\left(i-k\right)!\left(k!\right)\text{Γ}\left(k-\alpha +1\right)\left(j-l\right)!{\left(l!\right)}^2\left(k+l-\alpha +1\right)}}$ (15)

Note that in $D^{\left(\alpha \right)}$ the first $\lceil \alpha \rceil$ rows, are all zero.

Proof. Using Equations (3), (4) and (6) I have

$\begin{array}{c}{D}^{\alpha }{p}_i\left(x\right)={\displaystyle {\sum }_{k=0}^i\frac{{\left(-1\right)}^{i+k}\left(i+k\right)!}{\left(i-k\right)!{\left(k!\right)}^2}{D}^{\alpha }\left(x^k\right)}\\ ={\displaystyle {\sum }_{k=\lceil \alpha \rceil }^i\frac{{\left(-1\right)}^{i+k}\left(i+k\right)!}{\left(i-k\right)!\left(k!\right)\text{Γ}\left(k-\alpha +1\right)}{x}^{k-\alpha }},i=\lceil \alpha \rceil ,\cdots ,m\end{array}$ (16)

Now, approximate $x^{k-\alpha }$ by $\left(m+1\right)$ terms of shifted Legendre series, I have

$x^{k-\alpha }\simeq {\displaystyle {\sum }_{j=0}^m{b}_{k,j}{P}_j\left(x\right)}$, (17)

where

$\begin{array}{c}{b}_{k,j}=\left(2j+1\right){\displaystyle \underset{0}{\overset{1}{\int }}{x}^{k-\alpha }{p}_j\left(x\right)\text{d}x}\\ =\left(2j+1\right)\underset{l=0}{\overset{j}{{\displaystyle \sum }}}\frac{{\left(-1\right)}^{j+l}\left(j+l\right)!}{\left(j-l\right)!{\left(l!\right)}^2}{\displaystyle \underset{0}{\overset{1}{\int }}{x}^{k+l-\alpha }\text{d}x}\\ =\left(2j+1\right)\underset{l=0}{\overset{j}{{\displaystyle \sum }}}\frac{{\left(-1\right)}^{j+l}\left(j+l\right)!}{\left(j-l\right)!{\left(l!\right)}^2\left(k+l-\alpha +1\right)}\end{array}$ (18)

Employing Equations (16)-(18) I get

$\begin{array}{c}{D}^{\alpha }{P}_i\left(x\right)\approx \underset{k=\lceil \alpha \rceil }{\overset{i}{{\displaystyle \sum }}}\underset{j=0}{\overset{m}{{\displaystyle \sum }}}\frac{{\left(-1\right)}^{i+k}\left(i+k\right)!}{\left(i-k\right)!\left(k!\right)\Gamma \left(k-\alpha +1\right)}{b}_{k,j}{P}_j\left(x\right)\\ ={\displaystyle {\sum }_{j=0}^m\left({\displaystyle {\sum }_{k=\lceil \alpha \rceil }^i{\theta }_{i,j,k}}\right)P_j\left(x\right)},i=\lceil \alpha \rceil ,\cdots ,m\end{array}$ (19)

where ${\theta }_{i,j,k}$ is given in Equation (15). Rewrite Equation (19) as a vector form I have

$D^{\left(\alpha \right)}{P}_i\left(x\right)\simeq \left[{\displaystyle {\sum }_{k=\lceil \alpha \rceil }^i{\theta }_{i,0,k}},{\displaystyle {\sum }_{k=\lceil \alpha \rceil }^i{\theta }_{i,1,k}},\cdots ,{\displaystyle {\sum }_{k=\lceil \alpha \rceil }^i{\theta }_{i,m,k}}\right]\varphi \left(x\right),i=\lceil \alpha \rceil ,\cdots ,m$ (20)

Also according to Lemma 1, I can write

$D^{\left(\alpha \right)}{P}_i\left(x\right)=\left[0,0,\cdots ,0\right]\varphi \left(x\right),i=0,1,\cdots ,\lceil \alpha \rceil -1$ (21)

A combination of Eqs. (20) and (21) leads to the desired result.

Remark. If $\alpha =n\in N$ Then Theorem 1 gives the same result as Equation (11).

4. Applications of the Operational Matrix of Fractional Derivative

Operational matrix of fractional derivative is applied to solve multi-order fractional differential equation. The existence and uniqueness and continuous dependence of the solution to this problem are discussed in [21].

4.1. Linear Multi-Order Fractional Differential Equation [22]

Consider the linear multi-order fractional differential equation

$D^{\alpha }y\left(x\right)=a_1{D}^{{\beta }_1}y\left(x\right)+\cdots +a_k{D}^{{\beta }_k}y\left(x\right)+a_{k+1}y\left(x\right)-a_{k+2}g\left(x\right)$ (22)

with initial conditions

$y^{\left(i\right)}\left(0\right)=d_i,i=0,\cdots ,n$ (23)

where $j=1,\cdots ,k+2$ are real constant coefficients and also $n<\alpha \le n+1$, $0<{\beta }_1<{\beta }_2<\cdots <{\beta }_k<\alpha$ and $D^{\alpha }$ denotes the Caputo fractional derivative of order $\alpha$

To solve problem (22) and (23) I approximate $y\left(x\right)$ and $g\left(x\right)$ by the shifted Legendre polynomials as

$y\left(x\right)\simeq {\displaystyle {\sum }_{i=0}^m{c}_i{P}_i\left(x\right)}=C^{\text{T}}\varphi \left(x\right)$ (24)

$g\left(x\right)\simeq {\displaystyle {\sum }_{i=0}^m{g}_i{P}_i\left(x\right)}=G^{\text{T}}\varphi \left(x\right)$ (25)

where vector $G={\left[g_0,\cdots ,g_m\right]}^{\text{T}}$ is known but $C={\left[c_0,\cdots ,c_m\right]}^{\text{T}}$ is an unknown vector. By using Equations (13) and (24) I have

$D^{\alpha }y\left(x\right)\simeq c^{\text{T}}{D}^{\alpha }\varphi \left(x\right)\simeq c^{\text{T}}{D}^{\left(\alpha \right)}\varphi \left(x\right)$ (26)

$D^{\beta j}y\left(x\right)\simeq c^{\text{T}}{D}^{\beta j}\varphi \left(x\right)\simeq c^{\text{T}}{D}^{\left(\beta j\right)}\varphi \left(x\right),j=1,\cdots ,k$ (27)

Employing Equations (24)-(27) the residual $R_m\left(x\right)$ for Equation (22) can be written as

$R_m\left(x\right)\simeq \left(c^{\text{T}}{D}^{\left(\alpha \right)}-c^{\text{T}}{\displaystyle {\sum }_{i=1}^k{a}_j{D}^{\left(\beta j\right)}}-c^{\text{T}}{a}_{k+1}-G^{\text{T}}{a}_{k+2}\right)\varphi \left(x\right)$ (28)

As in a typical tau method [23] I generate m − n linear equations by applying

$\langle R_m\left(x\right),P_j\left(x\right)\rangle ={\displaystyle {\int }_0^1{R}_m\left(x\right)P_j\left(x\right)\text{d}x}=0,j=0,1,\cdots ,m-n-1$ (29)

Also, by substituting Equations (11) and (24) in Equation (23) I get

$\begin{array}{l}y\left(0\right)=c^{\text{T}}\varphi \left(0\right)=d_0\\ y^{\left(1\right)}\left(0\right)=c^{\text{T}}{D}^{\left(1\right)}\varphi \left(0\right)=d_1\\ ⋮\\ y^{\left(n\right)}\left(0\right)=c^{\text{T}}{D}^{\left(n\right)}\varphi \left(0\right)=d_n\end{array}$ (30)

Equations (29) and (30) generate m − n and n +1 set of linear equations, respectively. These linear equations can be solved for unknown coefficients of the vector C. Consequently, y(x) given in Equation (24) can be calculated.

4.2. Treatment of Nonhomogeneous Boundary Conditions

To solve Equation (22) with respect to the following boundary conditions (for n is even),

$u^{\left(i\right)}\left(0\right)=a_i,u^{\left(i\right)}\left(L\right)=b_i,i=0,1,\cdots ,\frac{n}{2}-1$ (31)

The same technique described in Section 4.1 was applied, but the (n) set of linear equations resulting from (31) is changed to be obtained from

$u^{\left(i\right)}\left(0\right)=c^{\text{T}}{D}^{\left(i\right)}\varphi \left(0\right)=a_i,u^{\left(i\right)}\left(L\right)=c^{\text{T}}{D}^{\left(i\right)}\varphi \left(L\right)=b_i,i=0,1,\cdots ,\frac{n}{2}-1$ (32)

Equations (29) and (31) generate (N + 1) system of linear equations. This system can be solved to determine the unknown coefficients of the vector C.

4.3. Nonlinear Multi-Order Fractional Differential Equation

4.3.1. Consider the Nonlinear Multi-Order Fractional Differential Equation

$D^{\alpha }y\left(x\right)=F\left(x,y\left(x\right),D^{{\beta }_1}y\left(x\right),\cdots ,D^{{\beta }_k}y\left(x\right)\right)$, (33)

with initial conditions

$y^{\left(i\right)}\left(0\right)=di,i=0,\cdots ,n$ (34)

where $n<\alpha \le n+1,0<{\beta }_1<{\beta }_2<\cdots <{\beta }_k<\alpha$ and $D^{\alpha }$ denotes the Caputo fractional derivative of order $\alpha$ . It should be noted that F can be nonlinear in general.

In order to use shifted Legendre polynomials for this problem, I first approximate $y\left(x\right),D^{\alpha }y\left(x\right)$ and $D^{\beta j}y\left(x\right)$ for $j=0,\cdots ,k$ as Equations (24), (26) and (27) respectively. By substituting these equations in Equation (33) I get

$C^{\text{T}}{D}^{\left(\alpha \right)}\phi \left(x\right)\approx F\left(x,C^{\text{T}}\phi \left(x\right),C^{\text{T}}{D}^{\left({\beta }_1\right)}\phi \left(x\right),\cdots ,C^{\text{T}}{D}^{\left({\beta }_k\right)}\phi \left(x\right)\right)$ (35)

Also, by substituting Equations (11) and (24) in Equation (34) I obtain

$y\left(0\right)=c^{\text{T}}\varphi \left(0\right)=d_0,$

$y^{\left(i\right)}\left(0\right)=c^{\text{T}}{D}^{\left(i\right)}\varphi \left(0\right)=d_i,i=1,2,\cdots ,n$ (36)

To find the solution y(x), I first collocate Equation (35) at m − n points. For suitable collocation points I use the first (m-n) shifted Legendre roots of $P_{m+1}\left(x\right)$ . These equations together with Equation (36) generate a system of (.m + 1) nonlinear equations which can be solved using Newton’s iterative method. Consequently y(x) given in Equation (24) can be calculated.

4.3.2. Boundary Value Problem

Consider the nonlinear FDE (33) with boundary conditions (31). I apply the same technique described in Section 4.3.1, but Equation (36) shall be changed to be (32)., I have a system of (N + 1) nonlinear algebraic equations, which can be solved using Newton’s iterative method.

5. Numerical Results

The presented method in the previous two sections has been applied to solve some examples. In this section, the results for the examples are shown along with their figures and a comparison with the analytical results was also presented in order to show the effectiveness of the technique.

Example 1. [24] Consider the linear fractional-order IVP:

$y^{\left(2\right)}\left(x\right)+3y^{\left(1\right)}\left(x\right)+2y^{\left(0.1379\right)}\left(x\right)+y^{\left(0.0159\right)}\left(x\right)+5y\left(x\right)=f\left(x\right),x\in \left(0,1\right),$ (37)

Subject to the initial conditions

$y\left(0\right)=1,y^{\prime }\left(0\right)=0$ (38)

where f(x) is chosen such that the exact solution of (37) is $y\left(x\right)=1+\frac{x^2}{2}$

By applying the technique described in Section 4.1 with m = 3, solution approximated as

$y\left(x\right)=c_0{p}_0\left(x\right)+c_1{p}_1\left(x\right)+c_2{p}_2\left(x\right)+c_3{p}_3\left(x\right)=c^{\text{T}}\phi (x)$

Here, I have

$D^{\left(1\right)}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 2& 0& 0& 0\\ 0& 6& 0& 0\\ 2& 0& 10& 0\end{array}\right)$ $D^{\left(2\right)}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 12& 0& 0& 0\\ 0& 60& 0& 0\end{array}\right)$

$D^{\left(1.379\right)}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 1.1313681331& 1.0223463203& -0.0608397666& 0.0199340611\\ -1.0223463203& 0.36356929917& 1.196322451& -0.093190361\\ 1.07052836676& -0.17397613077& 0.3074747759& 1.2963937228\end{array}\right)$

$D^{\left(0.0159\right)}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 1.01472967373& 1.0039162279& -0.00667748802491& 0.00190548426\\ -1.0039162279& 0.03644683& 1.0231320558& -0.00944784045\\ 1.0080521857& -0.01921582796& 0.02742824636& 1.03280174727\end{array}\right)$

$G=\left(\begin{array}{c}9.2199297\\ 4.10949395\\ 0.918660451\\ 0.03295167\end{array}\right)$

Therefore using Equation (29) I obtain

$-9.2199297+5c_0+9.2774659c_1+8.95139113c_2+9.1491089c_3=0$ (39)

$-4.1094939+8.04860886c_1+18.763585c_2+59.6328319c_3=0$ (40)

Now, by applying Equation (30) I have

$-1+c_0-c_1+c_2-c_3=0$ (41)

$2c_1-6c_2+12c_3=0$ (42)

Finally by solving Equations (39)-(42) I got

$c_0=1.17436533$ ; $c_1=0.2662085$ ; $c_2=0.09495018$ ; $c_3=0.003107$

Thus I can write

$\begin{array}{c}y\left(x\right)=\left(\begin{array}{cccc}1.17436533& 0.2662085& 0.09495018& 0.003107\end{array}\right)\\ \times \left(\begin{array}{c}1\\ -1+2x\\ 1-6x+6x^2\\ -1+12x-30x^2+20x^3\end{array}\right)\\ =1+\frac{x^2}{2}\end{array}$

which is the exact solution.

Example 2. Consider the boundary value problem

$D^{\frac{3}{2}}y\left(x\right)+y\left(x\right)=x^5-x^4+\frac{128}{7\sqrt{\text{π}}}{x}^{3.5}-\frac{64}{5\sqrt{\text{π}}}{x}^{2.5}$

subject to boundary conditions (43)

$y\left(0\right)=0,y\left(1\right)=0$

where the exact solution of this problem is $y\left(x\right)=x^4\left(x-1\right)$ .

This fractional boundary value problem is solved by applying the method described in Section 4.2 by using shifted Legendre expansion and its operational matrices of derivatives with m = 5. Using (43) four linear equations obtained, and by applying boundary condition I had two linear equations. By solving this linear system I got the unknown vector C. By substituting this vector in Equation (43), I had the exact solution.

$y\left(x\right)=\underset{i=0}{\overset{5}{{\displaystyle \sum }}}{c}_i{p}_i(x)$

Then, the 6 unknown coefficients will be in the form

$c_0=-0.033333$ ; $c_1=-0.042857$ ; $c_2=0.0119047$ ; $c_3=0.0388888$ ; $c_4=0.02142857$ ; $c_5=0.00396825$

Therefore, I can write

$\begin{array}{c}y\left(x\right)=\left(\begin{array}{cccccc}-0.033333& -0.042857& 0.0119047& 0.0388888& 0.02142857& 0.00396825\end{array}\right)\\ \times \left(\begin{array}{c}1\\ -1+2x\\ 1-6x+6x^2\\ -1+12x-30x^2+20x^3\\ 1-20x+90x^2-140x^3+70x^4\\ -1+30x-210x^2+560x^3-630x^4+252x^5\end{array}\right)\\ =x^4\left(x-1\right)\end{array}$

Numerical results will not be presented since the exact solution is obtained.

Example 3. [25] Consider the equation

$D^{2}u\left(x\right)+D^{\frac{3}{4}}u\left(x\right)+u\left(x\right)=x^3+6x+\frac{128}{15\Gamma \left(\frac{1}{4}\right)}{x}^{\frac{9}{4}}$ (44)

Subject to initial conditions $u\left(0\right)=0,u^{\prime }\left(0\right)=0$

where the exact solution of this problem is $u\left(x\right)=x^3$

By applying the technique described in Section 4.1 with m = 3, the approximate solution and the right hand side may be written in the form

$u\left(x\right)\simeq {\displaystyle {\sum }_{i=1}^3{c}_i{p}_i\left(x\right)}=c^{\text{T}}\varphi (x)$

$g\left(x\right)\simeq {\displaystyle {\sum }_{i=1}^3{g}_i{p}_i\left(x\right)}=G^{\text{T}}\varphi (x)$

Here, I have

$D^{\left(2\right)}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 12& 0& 0& 0\\ 0& 60& 0& 0\end{array}\right)$ $G=\left(\begin{array}{c}\frac{13}{4}+\frac{512}{195\text{Gamma}\left[\frac{1}{4}\right]}\\ 3\left(\frac{23}{20}+\frac{1536}{1105\text{Gamma}\left[\frac{1}{4}\right]}\right)\\ 5\left(\frac{1}{20}+\frac{512}{1547\text{Gamma}\left[\frac{1}{4}\right]}\right)\\ 7\left(\frac{1}{140}+\frac{512}{38675\text{Gamma}\left[\frac{1}{4}\right]}\right)\end{array}\right)$

$D^{\left(3/4\right)}=\left(\begin{array}{cccc}0& 0& 0& 0\\ \frac{8}{5}& \frac{8}{15}& -\frac{8}{39}& \frac{392}{3315}\\ -\frac{8}{15}& \frac{216}{65}& \frac{56}{51}& -\frac{88}{195}\\ \frac{272}{195}& -\frac{48}{85}& \frac{1264}{273}& \frac{5488}{3315}\end{array}\right)$

Therefore, using Equation (4.8) I obtain

$\frac{8c_1}{5}+\frac{172c_2}{15}+\frac{272c_3}{195}-\frac{13}{4}-\frac{512}{195\text{Gamma}\left[\frac{1}{4}\right]}+c_0=0$ (45)

$\left(\frac{8c_1}{15}+\frac{216c_2}{65}+\frac{5052c_3}{85}\right)-\frac{153}{20}-\frac{22016}{3315\text{Gamma}\left[\frac{1}{4}\right]}+c_1=0$ (46)

Now, by applying Equation (4.9), I have

$c^{\text{T}}\phi \left(0\right)=c_0-c_1+c_2-c_3=0$ (47)

$c^{\text{T}}{D}^{\left(1\right)}\phi \left(0\right)=2c_1-6c_2+12c_3=0$ (48)

Finally, by solving linear system of four Equations (45)-(48), I obtained

$u\left(x\right)=\left(\begin{array}{cccc}0.25& 0.45& 0.25& 0.05\end{array}\right)\left(\begin{array}{c}1\\ -1+2x\\ 1-6x+6x^2\\ -1+12x-30x^2+20x^3\end{array}\right)=x^3$

which is the exact solution.

It is clear that in Examples 1 - 3 the present method can be considered as an efficient method.

6. Conclusion

Shifted Legendre approximation method for solving higher order fractional differential equations has been presented. These equations are transformed to a system of algebraic equations to provide a matrix representation. The solution is expressed as a truncated Legendre series, and so it can be easily evaluated for arbitrary values by using computer program. From illustrative examples, it can be seen that this matrix approach can obtain very accurate and satisfactory results. The solution obtained is in very excellent agreement with the already existing ones and shows that this approach can solve the problems effectively. Comparisons between approximate solutions and analytical solutions illustrate the validity and the great potential of the technique.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Saadatmandi, A., Razzaghi, M. and Dehghan, M. (2005) Hartley Series Approximations for the Parabolic Equations. International Journal of Computer Mathematics, 82, 1149-1156.
[2]	Saadatmandi, A. and Dehghan, M. (2006) A Tau Method for the One-Dimensional Parabolic Inverse Problem Subject to Temperature over Specification. Computers & Mathematics with Applications, 52, 933-940.
[3]	Saadatmandi, A. and Dehghan, M. (2008) Numerical Solution of a Mathematical Model for Capillary Formation in Tumor Angiogenesis via the Tau Method. Communications in Numerical Methods in Engineering, 24, 1467-1474.
[4]	Saadatmandi, A. and Dehghan, M. (2007) Numerical Solution of the One-Dimensional Wave Equation with an Integral Condition. Numerical Methods for Partial Differential Equations, 23, 282-292. https://doi.org/10.1002/num.20177
[5]	Momani, S. and Shawagfeh, N.T. (2006) Decomposition Method for Solving Fractional Riccati Differential Equations. Applied Mathematics and Computation, 182, 1083-1092. https://doi.org/10.1016/j.amc.2006.05.008
[6]	Gejji, V.D. and Jafari, H. (2007) Solving a Multi-Order Fractional Differential Equation. Applied Mathematics and Computation, 189, 541-548. https://doi.org/10.1016/j.amc.2006.11.129
[7]	Wang, Q. (2006) Numerical Solutions for Fractional KdV-Burgers Equation by Adomian Decomposition Method. Applied Mathematics and Computation, 182, 1048-1055. https://doi.org/10.1016/j.amc.2006.05.004
[8]	Inc, M. (2008) The Approximate and Exact Solutions of the Space and Time-Fractional Burgers Equations with Initial Conditions by Variational Iteration Method. Journal of Mathematical Analysis and Applications, 45, 476-484. https://doi.org/10.1016/j.jmaa.2008.04.007
[9]	Momani, S. and Odibat, Z. (2006) Analytical Approach to Linear Fractional Partial Differential Equations Arising in Fluid Mechanics. Physics Letters A, 355, 271-279. https://doi.org/10.1016/j.physleta.2006.02.048
[10]	Odibat, Z. and Momani, S. (2006) Application of Variational Iteration Method to Nonlinear Differential Equations of Fractional Order. International Journal of Nonlinear Sciences and Numerical Simulation, 7, 271-279. https://doi.org/10.1515/IJNSNS.2006.7.1.27
[11]	Hashim, I., Abdulaziz, O. and Momani, S. (2009) Homotopy Analysis Method for Fractional IVPs. Communications in Nonlinear Science and Numerical Simulation, 14, 674-684. https://doi.org/10.1016/j.cnsns.2007.09.014
[12]	Liua, F., Anh, V. and Turner, I. (2004) Numerical Solution of the Space Fractional Fokker-Planck Equation. Journal of Computational and Applied Mathematics, 166, 209-219. https://doi.org/10.1016/j.cam.2003.09.028
[13]	Esmaeili, S. and Shamsi, M. (2011) A Pseudo-Spectral Scheme for the Approximate Solution of a Family of Fractional Differential Equations. Communications in Nonlinear Science and Numerical Simulation, 16, 3646-3654. https://doi.org/10.1016/j.cnsns.2010.12.008
[14]	Khader, M.M., Sweilam, N.H. and Mahdy, A.M.S. (2011) An Efficient Numerical Method for Solving the Fractional Diffusion Equation. Journal of Applied Mathematics & Bioinformatics, 1, 1-12.
[15]	Mohammadi, F. (2014) Numerical Solution of Bagley-Torvik Equation Using Chebyshev Wavelet Operational Matrix of Fractional Derivative. International Journal of Applied Mathematics and Mechanics, 2, 83-91.
[16]	Kilbas, A.A., Srivastava, H.M. and Trujillo, J.J. (2006) Theory and Applications of Fractional Differential Equations. Elsevier, San Diego.
[17]	Das, S. (2008) Functional Fractional Calculus for System Identification and Controls. Springer, New York.
[18]	Momani, S. and Noor, M.A. (2006) Numerical Methods for Fourth-Order Fractional Integro-Differential Equations. Applied Mathematics and Computation, 182, 754-760. https://doi.org/10.1016/j.amc.2006.04.041
[19]	Podlubny, I. (2002) Geometric and Physical Interpretation of Fractional Integration and Fractional Differentiation. Fractional Calculus and Applied Analysis, 5, 367-386.
[20]	Diethelm, K., Ford, N.J., Freed, A.D. and Luchko, Yu. (2005) Algorithms for the Fractional Calculus: A Selection of Numerical Methods. Computer Methods in Applied Mechanics and Engineering, 194, 743-773. https://doi.org/10.1016/j.cma.2004.06.006
[21]	Diethelm, K. and Ford, N.J. (2004) Multi-Order Fractional Differential Equations and Their Numerical Solution. Applied Mathematics and Computation, 154, 621-640. https://doi.org/10.1016/S0096-3003(03)00739-2
[22]	Dohaa, E.H., Bhrawyb, A.H. and Ezz-Eldien, S.S. (2011) A Chebyshev Spectral Method Based on Operational Matrix for Initial and Boundary Value Problems of Fractional Order. Computers and Mathematics with Applications, 62, 2364-2373. https://doi.org/10.1016/j.camwa.2011.07.024
[23]	Canuto, C., Hussaini, M.Y., Quarteroni, A. and Zang, T.A. (1988) Spectral Methods in Fluid Dynamic. Prentice-Hall, Englewood Cliffs. https://doi.org/10.1007/978-3-642-84108-8
[24]	Keshavarz, E., Ordokhani, Y. and Razzaghi, M. (2014) Bernoulli Wavelet Operational Matrix of Fractional Order Integration and Its Applications in Solving the Fractional Order Differential Equations. Applied Mathematical Modelling, 38, 6038-6051. https://doi.org/10.1016/j.apm.2014.04.064
[25]	Ford, N.J. and Connolly, J.A. (2009) Systems-Based Decomposition Schemes for the Approximate Solution of Multi-Term Fractional Differential Equations. Journal of Computational and Applied Mathematics, 229, 382-391. https://doi.org/10.1016/j.cam.2008.04.003