Numerical Solution of Fractional Differential Equations

In this article, two numerical techniques, namely, the homotopy perturbation and the matrix approach methods have been proposed and implemented to obtain an approximate solution of the linear fractional differential equation. To test the effectiveness of these methods, two numerical examples with known exact solution are illustrated. Numerical experiments show that the accuracy of these methods is in a good agreement with the exact solution. However, a comparison between these methods shows that the matrix approach method provides more accurate results.


Introduction
Fractional differential equations appear frequently in various fields involving science and engineering, namely, in signal processing, control theory, diffusion, thermodynamics, biophysics, blood flow phenomena, rheology, electrodynamics, electrochemistry, electromagnetism, continuum and statistical mechanics and dynamical systems. For more details on the applications of fractional differential Equations (see [1] [2] [3] [4] [5]). The concept of fractional calculus is now considered as a partial technique in many branches of science including physics (Oldham and Spanier [6]). Srivastava et al. [7] gave the model of under actuated mechanical system with fractional order derivative and Sharma [8] studied advanced generalized fractional kinetic equation in Astrophysics. Caputo [9] reformulated the more "classic" definition of the Riemann-Liouville fractional derivatives in order to use integer order initial conditions to solve his fractional order differential equation. Kowankar and Gangel [10] reformulated the Rie-mann-Liouville fractional derivative in order to differentiate nowhere differentiable fractal functions.
In general, most of fractional differential equations do not have exact solutions. Instead, analytical and numerical methods become increasingly important for finding solution of fractional differential equations. In recent years, many efficient methods for solving FDEs have been developed. Among these methods are the monotone iterative technique [11] [12], topological degree theory [13], and fixed point theorems [14] [15] [16]. Moreover, numerical solutions are obtained by the following methods: the Adomian decomposition method and the variational iteration method [17], homotopy perturbation method [18] [19], Haar wavelet operational method [20], Neural networks [21], and so forth. Very recently, Hamdan et al. [22] applied Haar Wavelet and the product integration methods to solve the fractional Volterra integral equation of the second kind. In addition, Saadeh [23] in her master thesis has employed several numerical methods for solving fractional differential equations. These methods are the Adomian decomposition method, Homotopy perturbation method, Variational iteration method and Matrix approach method. A comparison between these methods has been carried out. In spirit, our numerical methods, namely the homotopy perturbation and the matrix approach methods are an improvement to those methods presented in the master thesis by one of the authors of this article.
A comparison between these methods is carried out by solving some test examples using MAPLE software.
The paper is organized as follows: In Section (2) we recall some basic definitions and notions concerning fractional calculus. In Section (3), we introduce the homotopy perturbation method. The matrix approach method is addressed in Section (4). The proposed methods are implemented using numerical examples with known analytical solution by applying MAPLE software in Section (5).
Conclusions are given in Section (6).

Mathematical Preliminaries and Notions
In this section, we review some necessary definitions and mathematical preliminaries concerning fractional calculus that will be used in this work. Definition 1. [24] The Riemann-Liouville fractional integral operator of order The Riemann-Liouville derivative of fractional order p is defined as: Proof. See [26] for more details.

Homotopy Perturbation Method (HPM)
The fractional initial value problem in the operator form is: where i c is the initial conditions, L is the linear operator which might include other fractional derivative operators ( ) D β β α < , while the function g, the source function is assumed to be in 1 c − if α is an integer, and in 1 is not an integer. The solution ( ) f t is to be determined in 1 n c − . In virtue of [16], we can write Equation (8) in the homotopy form is an embedding parameter. If 0 p = , Equations (10) and (11) become when 1 p = , both Equations (10) and (11) yields the original FDE Equation (8). The solution of Equation (8) is: Substituting 1 p = in Equation (13) then we get the solution of Equation (8) in the form: Substituting Equation (13) into Equation (11) and collecting all the terms with the same powers of p, we get: and so on.
Following [6], we can write the first three terms of the homotopy perturbation method solution as: The general form of the HPM solution is: The homotopy perturbation solution takes the general form:

Left-Sided Fractional Derivative
Consider a function ( ) Let us take equidistant nodes with the step size ( ) Using the backward fractional difference approximation for the derivative at the points , 0,1, , Equation (21) is equivalent to the following matrix form [11]: In Equation (22), the column vector of functions the result is the column vector of approximated values of the fractional derivatives ( ), 0,1, , The generating function for the matrix is Since N A α and N A β are lower triangular matrices then we have

Right-Sided Fractional Derivative
Then its right sided fractional derivative of real order α ( ) Thus we get the discrete analogue of the right sided fractional differentiation with the step size The generating function for the matrix N G α is the same for N A α ( ) ( ) The transposition of the matrix N A α gives the matrix N G α and the opposite holds:

Numerical Examples and Results
In this section, in order to examine the accuracy of the proposed methods, we solve two numerical examples of fractional differential equations. Moreover, the numerical results will be compared with exact solution.
x t x t px t p x t p x t = + + + + : Applying α Ω and the inverse operator of D α , on both sides of Equation (33) and using the definition of Riemann-Liouville fractional integral operator ( α Ω ) of order 0 α ≥ we obtain: x Hence the solution of Equation (30) is: Now, we implement Algorithm 1 to solve Equation (30) using the matrix approach method. Table 1 displays the exact and numerical results using the matrix approach method with 1.9 α = and 51 N = . The maximum error with 51 N = is 0.039016195358901. Figure 1(a) compares both the exact and numerical solutions for the fractional differential Equation (30). Moreover, Figure 1(b) shows the absolute error between exact and numerical solutions.
The solution of Equation (38) has the form: x t x t px t p x t p x t = + + + + : 0 : : : Applying α Ω and the inverse operator of D α , on both sides of Equation (40), then we get: and using the definition of Riemann-Liouville fractional integral operator ( α Ω ) of order 0 α ≥ we obtain:   Figure 4(a) compares both the exact and numerical solutions for the (37). Moreover, Figure 4(b) shows the absolute error between exact and numerical solutions. For 60 N = and 70 N = , Figure 5(a) and Figure 6(a) compare both the exact and numerical solutions for the fractional differential Equation (30). Moreover, Figure 5(b) and Figure 6(b) show the absolute error between exact and numerical solutions.

Conclusion
In this article, two numerical techniques namely, the homotopy perturbation method and the matrix approach method have been proposed and implemented to solve fractional differential equations. The accuracy and the validity of these techniques are tested with some numerical examples. The results show clearly that both techniques are in a good agreement with the analytical solution. According to numerical results mentioned in tables and figures, we conclude that the matrix approach method provides more accurate results than its counterpart and therefore is more advantageous. In addition, we strongly believe that the matrix approach method is regarded to be one of the most effective methods among the other methods mentioned in the literature. It is known for its fast converges and accuracy.