Numerical Solutions of Fractional Differential Equations by Using Fractional Explicit Adams Method

Differential equations of fractional order are believed to be more challenging to compute compared to the integer-order differential equations due to its arbitrary properties. This study proposes a multistep method to solve fractional differential equations. The method is derived based on the concept of a third-order Adam–Bashforth numerical scheme by implementing Lagrange interpolation for fractional case, where the fractional derivatives are defined in the Caputo sense. Furthermore, the study includes a discussion on stability and convergence analysis of the method. Several numerical examples are also provided in order to validate the reliability and efficiency of the proposed method. The examples in this study cover solving linear and nonlinear fractional differential equations for the case of both single order as α ∈ (0, 1) and higher order, α ∈ [1, 2), where α denotes the order of fractional derivatives of Dαy(t). The comparison in terms of accuracy between the proposed method and other existing methods demonstrate that the proposed method gives competitive performance as the existing methods.


Introduction
Fractional calculus, particularly fractional differential equation (FDE), has significant applications, which thus plays a crucial role in various fields of science and engineering such as signal processing [1], control theory [2], modeling of materials [3] and diffusion processes [4]. The work by Woon [5] mentions the important implications of mathematical applications in other sciences such as physics. In addition, a study by Bagley and Torvik [6] remarks that fractional modelings are also extensively used in the field of viscoelasticity due to its ability to describe the high-frequency behavior of many viscoelastic materials.
First and foremost, we study the fractional initial value problem (FIVP) in the form [7], C D α t 0 y(t) = f (t, y(t)), y(t 0 ) = y 0 (1) where 0 < α < 1 is the fractional order and C D α t 0 denotes the fractional Caputo's α−derivative operator C D α t 0 = RL D α t 0 (y(t) − y(t 0 )) with RL D α t 0 (y(t)) is the Riemann-Liouville differential operator defined as [8]: The definition in Equation (2) has been applied especially in seeking for an analytical solution. However, when it comes to the real application, it might be very challenging. This is because, as pointed out by [9], specific additional conditions are needed to solve a differential equation in order to obtain a unique solution. These additional conditions for the Riemann-Liouville fractional derivative constitute a certain fractional derivative of unknown solution at the initial points which might result in an unclear physical meaning. Due to this reason, in the present work, we consider the fractional Caputo's α-derivative, C D α t 0 (in sequel, we shall simply denote as D α ), which is defined as [10] D α y(t) where α > 0 and m = α ; see also [11].
In literature, several analytical and numerical methods have been developed over the past few decades, namely the variational iteration method [12,13], homotopy perturbation method [14,15] and the Adams-Bashforth-Moulton method [10] to solve the fractional differential and integral equations. However, in some cases, solving differential equations numerically has proven to be more efficient and convenient compared to analytical solutions especially when dealing with huge and complex problems. Therefore, many researchers have developed numerous numerical methods to solve various kinds of differential equations. For example, research by Podlunby [16] used the first order finite difference numerical method and managed to solve FDE problems with O(h) accuracy. Following that, Gorenflo [17] proposed a second order difference method to solve FDE and managed to achieve the desired accuracy. Another well-known numerical method of a predictor-corrector type has been developed by Diethelm et al. [18] In the study, they developed an algorithm of P(EC) m E (Predict-Evaluate-Correct-Evaluate) where m is the number of iterations to solve linear and nonlinear FDE. In addition, Galeone and Garrappa [19] present a study on multistep methods for differential equations of fractional order that concerns numerical treatment of FDE on both implicit and explicit types. They managed to prove that explicit methods in the treatment of FDE give a good numerical solution with good stability analysis. Later, Blaszczyk and Leszczynski [20] proposed a study on FDE of higher order with a mixture of integer and Caputo derivatives using Euler's method. They modified the discrete form of the Caputo derivative being dependent on a range of the parameter, α and found that, when the range increases, the number of discrete equations occurring in the algorithm also increases.
The focus of this study is to derive an explicit multistep method based on the concept of the third-order Adams-Bashforth method, where the derivation of the proposed method is given in Section 2. Furthermore, the analysis of the stability and convergence is demonstrated in Section 3. Following from there, the implementation of the proposed method is shown in Section 4. The numerical results for solving six examples of FDEs are presented in Section 5, where it also includes a discussion on the numerical results obtained to illustrate the efficiency and effectiveness of the proposed method.

Fractional Explicit Adams Method of Order 3
This section will introduce a derivation of the fractional explicit Adams method of order 3 (FEAM3). In the first step, consider the FIVP is in the form [21]: It is well-known that FIVP in Equation (4) can be rewritten in the form equivalent to the Volterra integral equation [21] as: Note that, according to Diethelm [22], it is common to consruct methods for FDE by taking methods for classical (typically first-order) equations and generalizing the concept in an appropriate way. Therefore, we have Equation (4) as: Next, simplify Equation (5) as [23]: Now, we propose an approximate solution that involves approximation at the points t = t n and t = t n+1 . Thus, i. If t = t n ; ii. If t = t n+1 ; Subtracting Equation (9) from Equation (8) will yield: The proposed method is of order 3; therefore, taking the Lagrange interpolation with three interpolating functions of F n , F n−1 and F n−2 are required to evaluate the approximate solutions as follows: Next, let: Substituting Equations (11) and (12) into Equation (10), we obtain: i. The first fractional integral is evaluated as: Now, consider making substitution y = t n+1 − t, dy = −dt; then: Thus, the computation for the first fractional integral is given by: ii. Similar to the above steps, the second fractional integral can be computed as in the following form: Next, making the changes y = t n − t and dy = −dt gives: Therefore, we obtain the computation for the second fractional integral Now, if we substitute Equations (16) and (19) into Equation (10), then we obtain a numerical scheme for FEAM3 as follows: Therefore, Equation (20) is the proposed numerical scheme known as FEAM3. The method analysis and performance evaluation of FEAM3 are discussed in the latter section.

Definition 1 ([24]
). Linear multistep method is said to be of order p if, C 0 = C 1 = . . . = C p = 0 and C p+1 = 0. The formula used to calculate the constant C p is given by: where k is the order of the proposed method, α and β is the coefficient obtained from the proposed method. Note that C p+1 is the error constant of the method.
In order to check the order of the method Equation (20), we consider the general formulation of fractional linear multistep method for the solution of Equation (4) [19] as: where α j and β j are real parameters and α denotes the fractional order. By comparing Equations (20) and (22), we obtain α j and β j as follows: Substitute Equation (23) into Equation (21) and it gives: Therefore, the method is proven to be of order 3 with error constant C 4 = 3 8 .

Stability Analysis
The following test equation [7] was considered in order to analyze the stability properties of FEAM3: where the exact solution can be expressed in terms of the Mittag-Leffler function: Substitute test Equation (25) into the numerical method in Equation (20), which resulted in the stability polynomial of: whereh = λh α . Next, taking Equation (26) into consideration, the stability region of FEAM3 is drawn as shown in Figures 1-3 by using Maple software for different values of α. The planes are separated between the imaginary and the real field (horizontal axis, labelled as Re, represents Real while vertical axis, labelled as Im, represents Imaginary).

Single Order
A fractional equation is called single order FDE when 0 < α < 1. For comparison purposes, this paper separates the analysis for both single and higher order FDE.
While investigating the stability analysis, we are interested to investigate the values ofh = λh α . To determine the stability region, we take into consideration whether every root of λ i is real and must satisfy |λ i | ≤ 1 as stated in [25]. Thus, as the result, Equation (26) proved to be stable inside the shaded region as shown in Figure 1. In addition, we can also see from Figures 1 and 2 that, as the value of α increases, the stability region approaches to the left and is symmetric to the real axis. As for the area of the stability regions, we can observe from Figure 2 that the region became larger as the α increased.

Higher Order
In order to be considered as higher order FDE, the value of α should be when 1 < α < 2. A similar strategy was applied to determine the stability region for higher order FDE where we were interested in investigating the values ofh = λh α of Equation (26).
As mentioned before, we take into consideration whether every root of λ i is real and must satisfy |λ i | ≤ 1. Therefore, Equation (26) proved to be stable inside the shaded region as shown in Figure 3 when 1 < α < 2. According to Figure 3, it shows that the region for higher order FDE also became larger as the α increased and approached 2. The regions are shown to become larger in the left imaginary axis and symmetric to the real axis.
Next, applying Lipschitz condition as in Theorem 1 and the assumption in Equation (31), we have: Rewriting Equation (32) based on Theorem 2, we obtain: From the above analysis, it can be seen that, as h → 0, it is proven that |d n+1 | ≤ |d n |; thus, |y n+1 | = |y * n+1 | and |y n | = |y * n |. In conclusion, Theorem 2 is satisfied. Hence, the proposed method, FEAM3, is proved to converge. Theorem 1 ([22]). Let f (t, y) be Lipschitz continuous at all points (t, y) in a region D, given by: such that, for every t, y, y * , the coordinates t, y, y * and (t, y * ) are both in D where where K is a constant that depends only on α and p as p ∈ (0, 1), 0 < α < 1 and lim h→0 y n = y(t n ).

Algorithm of the Method
This section includes the algorithm of the proposed method. The inputs of the programming are the values of endpoints, a and b, number of intervals, N, the value of α as well as the initial value, y 0 as 0 < α < 1 and y 0 , y 0 as 1 < α < 2. The developed algorithm for the method is illustrated in Algorithm 1.

9:
Note that, error = |y n − Y n |, where y n is the approximation solution and Y n is the exact 10: solution.

Numerical Results
In order to validate the efficiency of the proposed method, six tested FDE problems which consist of single order as 0 < α < 1 and higher order, 1 ≤ α < 2 were considered. The computation was done using C programming (CodeBlock). Below are the notations used in the tables:
where the exact solution is given by y(t) = t 2α − t 2 .

Example 5.
A problem of higher order nonlinear fractional differential equation [29]. where The results are tabulated when α = 1.50 where the exact solution is y(t) = t 5 − 3t 4 + 2t 3 .
Tables 1 and 2 demonstrate the absolute error for solving a simple linear FDE of variable coefficient where the exact solution is the Mittag-Leffler function. Table 1 shows the absolute error of single order FDE as α = 0.30, 0.50, 0.70, while Table 2 presents the absolute error of higher order FDE as α = 1.30, 1.70, 1.90 at different step size, h = 10 −3 , 10 −4 . Based on these tables, it can be seen that, for each step size, h, the absolute error decreases as the order of FDE, and α increases. In addition, Table 3 shows the comparison in terms of absolute error at point, t = 1.0 when solving Example 1 between FAM and FEAM3 as α = 0.10, 0.90, 1.25, 1.85 at N = 10, 20, 40, 80, 160, 320. From the table, it can be seen that FEAM3 managed to produce comparable results as FAM. The performance graph for solving single and higher order of Example 1 are shown in Figures 4 and 5 respectively. The graphs illustrate that, for both cases, the approximate solution clearly approaches the exact solution when N increases.       Table 4 displays the absolute error at each point, t for solving an initial value problem of FDE in single order as α = 0.30, 0.50, 0.70 at different step size, h = 10 −3 , 10 −4 . The table shows that FEAM3 managed to perform well, whereby, as h decreases and α approaches 1.0, better accuracy was obtained. In order to observe the efficiency of FEAM3, Table 5, which displays the absolute error at point, t = 1 when solving Example 2 between FEAM3 and FLMM-3 when α = 0.7, 0.8 at h = 0.1, 0.01, 0.001 is also included. Based on the table, it shows that a comparable result is obtained between FEAM3 and FLMM-3. For better analysis, the performance graph for Example 2 is included in Figure 6 for α = 0.50 at N = 10, 100, 1000. Based on the graph, it can be seen that, as α increases, the approximate solutions approach the exact solution.  In addition, the numerical results for solving initial value problem of higher order FDE are tabulated in Table 6 in the form of absolute error when α = 1.30, 1.70, 1.90 as h = 10 −3 , 10 −4 . The table proved that FEAM3 is also able to perform well in solving nonlinear FDE, where the absolute error decreases when α increases and h decreases. On the other hand, Table 7 shows the absolute error of various α = 1.25 at different h between FEAM3 and ATPC for solving nonlinear FDE of Example 3, where FEAM3 is able to obtain a comparable result as the existing method, ATPC. Additionally, the graph of approximate solution at each point, t, for solving Example 3 when α = 1.50 for different N is shown in Figure 7. The graph highlights that the approximate solution does indeed approach the exact solution as α increases. Table 4. Absolute error at each point, t when solving Example 2 using FEAM3 for different step size, h and 0 < α < 1.       Next, Table 8 shows the numerical result for solving Example 4 when α = 1.0 , while Table 9 shows the results when α = 1.50 for solving Example 5. Both tables present the approximate solution and absolute errors at each point, t for different N. Based on these tables, it can be seen that better accuracy is obtained as N increases. For comparison purposes, Table 10 presents the comparison of absolute error between FEAM3 and SFMoPF at each point, t when α = 1.0 and N = 10 in solving Example 4, while Table 11 demonstrates the comparison of absolute error at each point, t between FEAM3 and 3-HOFLMSM in solving nonlinear higher order FDE of Example 5 when α = 1.50 and N = 100. From these tables, the results given by FEAM3 are seen to be comparable to the respective existing methods. The performance graph for both examples are shown in Figures 8 and 9, respectively, where, as α increases, the approximate solutions approach the exact solution.      This paper also includes solving fractional Riccati differential equation (FRDE). According to [30], it is well-known that the Riccati differential equation is concerned with applications in pattern formation in dynamic games, linear systems with Markovian jumps, diffusion problems, river flows, and econometric models. Therefore, many researchers have developed several analytical and numerical methods in solving FRDE problems, since it can be considered as one of the examples of application problems in FDE. Table 12 presents the result on absolute error for solving FRDE problem by using FEAM3 at different intervals, N = 10, 100, 1000 when α = 1.0. Based on the table, the absolute error for solving FRDE by using FEAM3 decreases as N increases. Thus, it gives meaning whereby, as the step size decreases, the approximate solutions approach the exact solutions. For comparison purposes, Table 13 shows the comparison result of solving the FRDE problem between FEAM3, FVIM, and MHPM when α = 1.0 and N = 10. According to the table, FEAM3 managed to give a comparable result as FVIM and MHPM. Furthermore, a graph shows the approximate solution at each point, and t for solving the FRDE problem is also presented in Figure 10, where the approximate solutions do indeed approach the exact solution as α increases. Therefore, this implies that FEAM3 is able to perform well in solving nonlinear FDEs.

Conclusions
This paper has proposed a numerical method known as a fractional explicit Adams method of order 3, FEAM3. The numerical results obtained for each example authenticate that FEAM3 is adequate to preserve accuracy in solving both linear and nonlinear FDEs. The results also validate the convergence analysis where the approximate solutions are indeed converged and approach the exact solution as the step size, h, decreases. Additionally, FEAM3 is proven to be capable of achieving comparable results as the existing methods in each example. Furthermore, this paper also includes solving problems of FDE for the case of single and higher order, where the results have shown that the increment in the value of α yields better accuracy in solving both linear and nonlinear FDE problems. Therefore, it is proven that FEAM3 is competent and reliable to act as an alternative method to solve different kinds of FDEs. Funding: The authors gratefully acknowledge the financial support of Research University Grant (Putra Impact Grant: UPM/800-3/3/1/9629200) from Universiti Putra Malaysia.

Conflicts of Interest:
The authors declare no conflict of interest.