Numerical Solutions of Conformable Fractional Differential Equations by Taylor and Finite Difference Methods

Bu calismada yeni tanimlanan conformable kesirli turevli denklemler icin guvenilir  ve etkili bir metot turettik. Kesirli Taylor acilimindan ilk once Euler ve Taylor metodunu gelistirdik. Bu Taylor acilimi baslangic noktasindan farkli bir noktada acilmis genellestirilmis Taylor serisiridir. Ongorulen metotlar daha etkili ve hizli oldugunu birinci dereceden kesirli diferansiyel denklemlere ve ikinci dereceden salinimli kesirli diferansiyel denklemlere  uygulayarak gosterdik. Ikinci metodumuz ise kesirli diferansiyel denklemi zayif tekil integral denklemine donusturup, carpim intagrasyon kuralini uygulayarak cozmek olacaktir. Bu yeni tanimda ozel tanimli fonksiyonlar olmadigi icin, metotlar daha dogru sonuc verecek ve bilgisayar programlamasi daha kolay olacaktir. Bu ongorulen metotlarin kararlilik ve yakinsakliklari ispatlanmis olup, teorik sonuclari destekleyen sayisal ornekler verilmistir.


Introduction
Recently, fractional differential equations (FDEs) become more attractive and have been developed in theory and applications in science and mathematics. Applications and the theories of fractional differential equations increasingly get more attention nowadays both in science and engineering. Some applications of FDEs can be founded in chemistry, mechanics, physics, control theory and so on. For more details on the application of FDEs, we refer the reader to the references [1][2][3]. Unlike the ordinary differential equations, the analytic solutions of fractional differential equations may not be available. Thus, efficient and reliable numerical methods for solving fractional differential equations are essential and important. Almost all the definitions of the fractional derivative have been defined globally and in non local sense so that they involve fractional integral equations with weakly singular kernels and some special functions such as Gamma and Mittag-Leffler functions. All these definitions does not obey some standard rules and important properties of ordinary derivative such as chain rule or semi-group property. Recently, much simpler and compatible definition of fractional derivative obeying chain rule and semi-group properties based on the basic limit processing so called the conformable fractional derivative has been given in [4]. This new definition of fractional derivative have came to the many researcher's attention and the fundamental properties and some applications of this new fractional derivative have been studied and developed [5]. Further developments and several application can be found in [6][7][8][9] and references therein.
Khalil et al. [4] define the conformable fractional deriva-tive of order α ∈ (0, 1] of a function f : [a, ∞) → R by An easy consequence of this definition is that if f has the classical derivative, then we have the following relation [4] (T a α f )(t) = (t − a) 1−α f (t), where f (t) is the classical derivative of f . We immediately see that the conformable fractional derivative of a constant function is zero. Some basic properties of this conformable fractional derivative can be found in [4,5] in details. This new definition intuitively is a natural extension of standard derivative to non-integer order. Unlike the existing definitions of fractional derivative, there are no special functions such as the Gamma, Beta and Mittag-Leffler functions that are not easy to evaluate and implement in the solutions. This conformable derivative has the physical interpretation as a modification of the classical derivative in direction and magnitude of physical quantity [10].
For easy presentation, we restrict the case when α ∈ (0, 1]. The results that we will find in this paper can be easily extended to the case α ≥ 1. The other commonly used fractional derivatives are the Riemman-Liouville and Caputo derivative defined by, respectively and The corresponding initial value problem can be written as [11] RL D α a y(t) = f (t, y(t)), α ∈ (n, n + 1], t ∈ (a, T ], y (k) (a) = y (k) 0 , k = 0, 1, 2, . . . , n − 1.
(6) Usually, the initial value problems for fractional differential equations are ill-posed because of singularity of the solution at the initial conditions. In general, the differential equation (6) is converted to the following Volterra integral equation If the solution of (6) is smooth enough, then it also solves (7) and vice versa. The integral equation (7) is singular if α ∈ (0, 1) and this singularity makes the numerical solution inefficient and reqıires some special techniques to approximate the solution of (7). On the other hand, the solution of (7) can be written as an expansion of integer and non-integer powers [12] As we see from the equation (8), the solution is not smooth at t = a and present mixed powers of integer and noninteger . This non-smooth property of the solution makes the numerical solution of the integral equation difficult to approximate. The main difference in between classical derivative and fractional derivative is the non local properties of the fractional calculus. This leads to intense computational methods and high order numerical methods that are very limited in literature. Numerical methods for (6), to our knowledge, are based on Riemann-Liouville or Caputo [12] and the Grünwald-Letnikov approach [13].
Several numerical methods such as finite difference [11], finite element [14] and spectral methods for numerical solution of (6) or (7) have been proposed and developed during the last few decades. In [15], the author discussed the stability of the numerical methods for the equation (7) and gave the the disc of stability of predictor-corrector methods. Diethelm et al. [16] proposed Adams-type predictor-corrector method for equation (7) and if the Caputo derivative of the solution is smooth enough, they gave the error bound for the method. Li and Zeng [17] discussed the finite difference method for fractional differential equations. Recently, in the book [18], finite difference methods and finite element methods have been studied and analyzed for fractional ordinary differential and partial differential equations. Usually, the weakly singular kernel of the Volterra type integral equations makes it difficult to have an efficient and high order numerical method. To overcome this inefficiency, the integral equation t a (t − x) α−1 f (x, y) dx have been approximated by choosing suitable quadrature numerical methods. Product integration rule that is a class of convolution quadrature is one of the methods to numerically solve this kind of integral equation introduced in [19]. Different quadrature rule gives the different numerical method such as fractional Euler and fractional Adams methods. In this article, we use two numerical approaches to the problem (3). We first derive a fractional power series that have been used to obtain fractional Euler scheme and high order Taylor numerical method for the equation (3). To the best of our knowledge, this is the first finite difference method for FDEs in the sense of conformable fractional derivative. The power series obtained are not dependent on the initial point or the point in which the conformable fractional derivative is defined. Our second approach is based on the following Volterra-type integral equation with 851 Ş. TOPRAKSEVEN / Numerical Solutions for Conf. Frac. Diff. Eq. by Taylor and Finite Dif. weakly singular kernel by converting FDEs (3) with the help of the relation (2) This form is significantly different from the one given in (7) because there is no delay argument in this formulation. Unless otherwise stated, we always assume α ∈ (0, 1) in this article and numerical methods we proposed here can be extended to α ∈ (n, n + 1] for any n ∈ N. Numerical methods for FDEs are costly with high computational times and have expensive storage since the number of operations increase at each time step because of the singular kernel in the equivalent Volterra integral equation. However, the newly defined fractional derivative is based on limit process and the number of operation increase linearly. Thus, developing numerical methods based on this newly defined fractional derivative is important and has many advantages compared to the existing methods. This paper is organized as follows. In Section 2 the basic definitions and background for FDEs are reviewed. In Section 3 we establish the existence and uniqueness of the problem (3) and reformulation of FDEs in terms of integral equation. Then we introduce numerical methods for solving Equation (3) with uniform meshes in section 4. In Section 5, we prove the stability and error estimates of these numerical schemes. Finally, various numerical examples are given to show that the numerical results confirm the theoretical findings in Section 6 and last section includes the concluding remarks. Throughout, the notations C and c, with or without a subscript, denote generic constants, which may differ at different occurrences, but are always independent of the mesh size.

Definitions and Background on Fractional Calculus
The left and right conformable fractional αth order derivative of a function f : [a, ∞) → R given by [4] ( The conformable fractional integral operator of order α is given by [4] I a Lemma 2.1. [5] If f : [a, ∞) → R is smooth, then we have the relations for α ∈ (0, 1] and t > a.
Lemma 2.2. [5] If y : [a, ∞) → R is infinitely α-differentiable function, for some α ∈ (0, 1] at a neighborhood of a point a , then y has the fractional power series expansion: (12) Here (T a α y) (k) (a) means the application of the fractional derivative k times.
We use the following Banach Contraction Principle ( see e.g., [20]) Let Y be a nonempty closed subset of a Banach space X such that for each n ≥ 1, there exists a constant c n ≥ 0 such that ∑ ∞ n=0 c n < ∞. If the operator A : Y → Y satisfies A n y − A n z ≤ c n y − z , for all n ∈ N, then A has a unique point y * so that Ay * = y *

Existence and Uniqueness of the Solution
In this section, we present the existence and uniqueness of the solution of equation (3). Similar to the classical ordinary differential equation, we first assume the function f is continuous on some domain for proving the existence of the problem and then to prove the uniqueness, we assume it has a a Lipschitz condition with respect to the second variable; i.e., with some constant L > 0 independent of x, y, and z. Now we ready to establish the following two theorems.  The proof of these two theorem will be based on the following observation.
Proof. With the help of (2), we can re-write the problem (3) as Now, integrating both side from a to t to obtain the desired result (13). Obviously, this set is not empty. Thus, we define the operator A : B δ → B δ given by We show that the operator A has a unique point that solves the equation Ay = y by showing that A n is a contraction operator for n sufficiently large. First, we demonstrate Ay is a continuous function. To do this, we choose a ≤ t 1 ≤ t 2 ≤ T and we write To this end, note that Finally we prove that A n is a contraction operator for n ∈ N. More precisely we have for t ∈ (a, T ) We prove this fact by induction on n. For n = 0, it is obvious. Assume that (15) is true for n = k − 1. Now, we write Using the induction hypothesis, we obtain According to Theorem 2.3, if the sum ∑ ∞ n=0 c n converges with c n = L n (t − a) αn α n n! , then the operator A has a unique point so that the problem (3) has a unique solution. In fact, the fractional calculus (see e.g., [5]) implies that the series ). Thus, the proof is completed.

Numerical Methods
In this section, we will derive the numerical schemes for approximating the problem (3). To do this, we introduce the following notations: For a given positive integer N, let t n = nh + a, n = 0, 1, . . . , N be a uniform meshes of the interval (a, T ] where step size h ( for the sake of simplicity assumed to be constant) is then given by h = t n − t n−1 . Let y(t n ) be approximated by y n at the point t = t n .
Numerical methods proposed in this work for the fractional differential equation (3) will be based on the Taylor expansion by the help of Lemma 2.2. We stress out that Lemma 2.2 allows to obtain a fractional power expansion for a function in terms of its Conformable fractional derivatives evaluated at the initial point a. We now obtain a similar Taylor expansion at any other point a 1 > a, so the expansion can be expressed independently from the initial point a.
Theorem 4.1. If y : [a, ∞) → R is infinitely α-differentiable function, for some α ∈ (0, 1] at a neighborhood of a point a 1 ∈ (a, ∞) , then y has the fractional power series expansion: where R 4 (t, a 1 , a) is the reminder term and and H = t − a and L = a 1 − a. Proof.
Using the power series expansion (12) for y(t), (T a α y)(t), (T a α y) (2) with t replaced by the points a 1 , we have where R j 4 is the reminder of the fractional Taylor series of (T a α y) ( j) (a) at the point a 1 for j = 1, 2, 3, 4. Then every expansions above are substituted in (12) up to order of fourth Taylor series and we obtain the equation (16). Furthermore, the reminder term R 4 (t, a 1 , a) is resulted from the reminder R j 4 of each Taylor expansion for j = 1, 2, 3, 4.

The Construction of Numerical Methods Based on Taylor Expansions
Similar to that done in the standard numerical approach to the Cauchy problems for ODE, we construct numerical schemes to solve FDE (3) based on the application of fractional Taylor power series. The first numerical method that is relatively simpler and easier is numerical integration method having a disadvantage of restricted time stepping. Suppose that y(t) is the solution of the conformable fractional differential equation (3) with α ∈ (0, 1]. By using Theorem 4.1, the fractional Taylor expansion of y(t n+1 ) at t = t n gives that by the FDEs (3), we have Next, we define Fractional Euler method that approximates the solution y(t) at the points t n , n = 0, 1, 2, . . . , N for the values y n = y(t n ) so that the algorithm is then given by Now, we estimate the reminder term R 1 (t n+1 ,t n , a) in terms of the step size h and the parameter α. By the fractional Taylor expansion (12), the reminder can be written as (see e.g., [5]) where b n defined by (18) and all the coefficients A 1 , A 2 and A 3 are bounded constants. Local truncation error at t n+1 can be defined by From (20), the truncation error τ(t n+1 ) = h α so that the error can be given by By a similar argument for the construction of Euler methods, we define high order fractional Taylor methods. These fractional Taylor methods are based on the fractional Taylor expansion (16) and have high order approximations up to desired order. Here, we derive second and fourth order fractional Taylor methods that are frequently used in applications. However, one can derive fractional Taylor methods of any order. A disadvantage of the high order approximations is that they require the fractional derivatives of the function up to the order of the method. Because of the singularity of the solution, the high order Taylor methods may not be available. Now, we derive 2α order fractional Taylor methods for the problem (3). We approximate the solution y(t) at the points t n , n = 0, 1, 2, . . . , N for the values y n = y(t n ) so that the algorithm is then given by where b n is defined by (18). The same idea above gives 4α order fractional Taylor method for approximating the problem (3); For n = 1, 2, . . . , N where b n is defined by (18) and The truncation error for the fractional Taylor methods can be obtained similarly as in the Euler Methods. However, we only give outline of the proof because of very long algebraic expressions. Usually, the reminder for the fractional Taylor series is of order h (N+1)α when the series have order of N. We first find all the reminder of the fractional Taylor series for y(t) up to N order and then substituting all the reminder will be of the same order or greater, thus the reminder of the fractional Taylor series for y(t n+1 ) at t = t n would be of order h (N+1)α . Thus, the local truncation error (21) will be of order h Nα . This outline generalizes the truncation error of classical Taylor method to the fractional Taylor method.
In the next section, we give some numerical results for the numerical methods discussed here. We test three methods described above for FDEs so that the exact solutions are available. We first test the numerical methods for the first order FDEs and secondly we give numerical results for second order FDEs that is converted to first order systems of FDEs so that we compare the obtained results with the previous results based on symplectic schemes given in [15] with Caputo definition.

Numerical Examples
The first problem to be tested is the fractional Cauchy problem given by Example 1.
The exact solution is given by y(t) = y 0 exp(λt α /α). This is a simple linear test problem so that we easily find fractional derivatives up to desired order. The successive fractional derivatives give that (T 0 α y) (k) (t) = λ k y(t) for each k ∈ N. Based on this observation, the numerical methods are given as follows: Fractional Euler Method Fractional Taylor of order 2α Method Fractional Taylor of order 4α Method where b n , 2 b n , 3 b n and 4 b n are defined by (25) and y E 0 = z 0 = w 0 = y 0 .
We plot the solutions of the numerical methods and exact solution for different α and the initial value y 0 = 1 in Figure 1. We take λ = −10. As we can see from the Figure 1, the fractional Taylor of order 4α closes to the exact solution and observe that as α increases, the errors decrease as expected due to the error depends on the fractional power of h. We emphasize here that if one uses the Caputo or Riemann-Lioville definition of fractional derivative, then the Cauchy problem (26) leads to the exact solution of the form y(t) = y 0 E α (λt α ) where E α (z) is the e Mittag-Leffler function [15] that generalizes the standard exponential function exp(z) for complex number z. Therefore, one needs some algorithm to evaluate this special function even in this simple case. However, we have a simple solution to the problem similar to the classical differential equation so that we believe that our methods are faster and more efficient to approximate FDEs.
The exact solution is y(t) = t 4 exp(−t α /α) for α ∈ (0, 1]. We take α = 0.5 and we show the errors and the estimated order of convergence of the fractional Euler and fractional Taylor methods in Table 1  By EOC, we show the estimated order of convergence that is given by the formula that log 2 ( e(h,T ) e(h/2,T ) ) where e(h, T ) is the error with the step-size h at t = T . We note that EOC of the fractional Taylor method of order 2α is O(h) and the fractional Taylor method of order 4α is O(h 2 ) as we expect. Moreover, the fractional Euler method is of order O(h) in this numerical solution because the exact solution is smooth. Next, we show that our methods are more efficient and faster than previous proposed finite difference methods in the literature by demonstrating the results on the second order FDEs given below.
Note that we have two initial conditions that depend on y(t) and the first conformable fractional derivative (T 0 α y)(t) at 855 Ş. TOPRAKSEVEN / Numerical Solutions for Conf. Frac. Diff. Eq. by Taylor and Finite Dif.  t = 0. This FDEs is known as Fractional Oscillator and studied in [21]. Similar to the idea in ODE, we first convert this FDEs to the system of first order FDEs as follows In [21], the authors have used the following numerical scheme with the Caputo definition where Γ(z) is the Gamma function. Note that numerical scheme (30) requires high number of operations that increases as N 2 after N steps since all the previous steps are needed for the next step. This increases the computation time. However, this is not the case for our method since the methods require only the previous step to get next step so the operations increase only linearly with N after N steps. Now, we approximate the system of FDEs (29) by proposed methods defined by Fractional Euler Method Fractional Taylor of order 2α Method Fractional Taylor of order 4α Method T y n+1 = T y n + b n T z n h α α − 2 b n T y n ω 2 h 2α where b n , 2 b n , 3 b n and 4 b n are defined by (25) and y E 0 = y T 0 = T y 0 = y 0 and z E 0 = z T 0 = T z 0 = z 0 . The exact solution of the second order FDEs (28) is given by In Figure 2, we show the exact solution and numerical solutions with their corresponding absolute errors for the problem (28). We take ω = 5,y(0) = 1 with z(0) = 0 and 856 Ş. TOPRAKSEVEN / Numerical Solutions for Conf. Frac. Diff. Eq. by Taylor and Finite Dif.  different values of α. Similar results to the first order FDEs solutions, we expect the errors decrease as α converges to 1. Also, we see that the fractional Taylor method of 4α order is much closer to the exact solution while the fractional Euler and fractional Taylor method of 2α order oscillate frequently as ω gets bigger.  Table 3. The maximum error for the fractional Euler and Taylor methods for for Example 28 at T = 4 with α = 0.9 and ω = 5

Numerical Schemes Based on Weakly Singular Integral Equations
Our second numerical approaches we introduce here depend heavily on the following integral: for n = 0, 1, . . . , N − 1, Now, we approximate f (x, y(x)) by choosing suitable numerical approximation g(x, y(x)) of f (x, y(x)) on the interval [t n ,t n+1 ] and different choice will lead to different numerical scheme for numerical solution of the problem (3) or equivalently (13) as shown below.
1. If we choose g(x, y(x)) ≈ f (t n , y(t n )) then we have an explicit method, i.e., y n+1 is given explicitly in terms of known quantities y n and f (t n , y n ) and we call this method as the fractional forward Euler (FFE) method defined as follows where b n = (n + 1) α − n α .
2. If we choose g(x, y(x)) ≈ f (t n+1 , y n+1 ) then we have an implicit method and we call this fractional backward Euler (FBE) method defined by where b n is defined as above.

857
3. If we choose g(x, y(x)) ≈ P(x), where P(x) is an interpolating polynomial as an approximation of f (x, y(x)), more precisely, if we choose g(x, y(x)) as follows (34) then we have the fractional trapezoid method is given by where a k,n+1 defined by Since the trapezoidal method is implicit and can not be solved directly, we naturally propose the following fractional Adams predictor-corrector method for (13) where a k,n+1 defined by (37).
High-order predictor-corrector methods for fractional differential equation have been proposed in [11]. The authors use the interpolation function for approximating f (x, y(x)) based on first and second degree Lagrange interpolation and spline functions. In [22], more similar predictor-corrector approach to ours have been introduced and analyzed for the time fractional Fokker-Planck equation.

Stability and Error Analysis
In this section, we will examine the stability and error estimates for the methods introduced above. Stability estimates are similar to the stability analysis of standard finite difference method for the fractional Euler methods and based on characteristic equation for the amplification factor y n = y 0 r n . The stability analysis of Adams methods means that a small change in the initial condition will not cause the huge error in the numerical solution. In the sequel, we need the following auxiliary lemmas for stability and error analysis.

Lemmas
Lemma 5.1. If α ∈ (0, 1], then we have for x ∈ [0, 1] and for n ≥ 1 For the inequality (42), by applying the mean value theorem for the function f (x) = ( n+x n+1 ) α , we observe that This leads to get The inequality (43) can be proved similarly so we omit the proof of it. Thus, the proof is now completed.
Lemma 5.3. If α ∈ (0, 1] and y(t) solves the equation (13) and f (t, y(t)) bounded on the domain then there is a constant C α independent of h so that the following inequality holds for small h with t ∈ [a, T − h], Proof. By the help of(13), we write Since f is bounded, there is a constant M so that | f | ≤ M, and we obtain t+h−a ) α , we can appeal the inequality (41) and finally have the desired result Hence, the proof is completed.

Stability Analysis
In this subsection, we prove the stability estimates for the methods we introduced in Section 4 for α ∈ (0, 1]. If α ≥ 1 then standard way of proving the stability analysis based on Gronwall's inequalities can be used. First, we prove the stability analysis of the fractional forward Euler method given by (32). We show the stability of the proposed methods by showing that a small perturbation in the initial conditions does not lead to substantial changes in the numerical solution as time progresses.
Theorem 5.7. If y k for k = 1, 2, . . . , n + 1 are the solutions of the fractional Euler methods (32) and f satisfies the Lipschitz condition with respect to second variable locally with a Lipschitz constant L, then the fractional Euler method (32) is conditionally stable.

859
Ş. TOPRAKSEVEN / Numerical Solutions for Conf. Frac. Diff. Eq. by Taylor and Finite Dif. Now using the Lipschitz condition, Lemma 5.1 and Lemma 5.2 and the fact that h α b n ≤ h α (n + 1) α−1 ≤ T α−1 h for n = 0, 1, . . . , N − 1, we get where M = |ŷ 0 |. By the same argument used in the previous theorem, we can prove the following Theorem 5.8. If y k for k = 1, 2, . . . , n + 1 are the solutions of the fractional Euler methods (33) and f satisfies the Lipschitz condition with respect to second variable locally with a Lipschitz constant L, then the fractional Euler method (33) is stable.
Next, we examine the stability analysis of the fractional Adams method and the technique that is similar to the given in [11] and based on the idea of perturbations in the initial condition does not lead to the larger error in the numerical solution.
Theorem 5.9. If y k for k = 1, 2, . . . , n + 1 are the solutions of the fractional Adams methods (38) and f satisfies the Lipschitz condition with respect to y locally with a Lipschitz constant L, then the fractional Adam's methods (38) are unconditionally stable.
This shows that the Fractional Adam's methods are unconditionally stable. Note that α = 1 than this FFE and FBE methods are reduced to the standard Euler methods and the trapezoidal method and our findings are consistent with the known facts.

Error Estimates
In this section, we will derive error estimates for the numerical methods we proposed in the previous section. We first give error analysis for the FFE method (32) in the following theorem. Theorem 5.10. If y(t) solves (13) and y n solves the FFE (32), and assume that f (t, y(t)) is sufficiently smooth and satisfies the Lipschitz condition with respect to second variable with a constant L then we have the following error estimate for FFE method for α ∈ (0, 1]. Proof. Let the error difference e n is defined as e n = y(t n ) − y n with e 0 = 0. Subtracting (32) from (13) gives the following error equation Therefore, we obtain Appealing Lemma (5.4) and Lemma (42) and using the Lipschitz condition, we obtain Now using the discrete Gronwall's inequality [23], we find that |e n | ≤ Ch α .
The proof is now completed. For more detailed Gronwall's inequality for fractional differential equation, we refer the reader to the paper [23].
With the same argument done in the proof of above theorem, we can prove the following theorem Theorem 5.11. If y(t) solves (3) and y n solves the FBE (33), and assume that f (t, y(t)) is sufficiently smooth and satisfies the Lipschitz condition with respect to second variable with a constant L then we have the following error estimate for FBE method for α ∈ (0, 1].
Next, we will prove an error analysis for the fractional Adam's methods in the following theorem Theorem 5.12. If y(t) solves (13) and y n solves the fractional Adam's methods (FAM) (38) and assume that f (t, y(t)) is sufficiently smooth and satisfies the Lipschitz condition with respect to second variable with a constant L then there is a constant C independent of h so that we have the following error estimate for FAM for α ∈ (0, 1]. Proof. Let e k = y(t k ) − y k be the error at t = t k , then subtracting the equation (13) from (38) and letting e 0 = 0, we find the following error equation To finish the proof, we need to find the bound for T i for i = 1, 2, 3. The first term T 1 is bounded by Lemma 5.6 so that Using the Lipschitz condition on f with a Lipschitz constant L, we have Finally we estimate the last term T 3 using the Lipschitz condition on f with a Lipschitz constant L and a n+1,n+1 = 1 as follows where we used the Lemma 5.5 and C 1 := LMT α α(α + 1) .
Thus, we get Now, the result follows from the well-known Gronwall's inequality (see e.g., [5]) Thus, the proof is now completed.

Numerical Examples
In this section, we will give some numerical examples that verify the theoretical findings that we have established in the previous sections.
In this example, we test the FFE (32) and FAM (38) and the numerical results are shown in Table 4 and Table 5. We measure the error by the maximum norm defined by e N := max 0≤k≤N |y(t k ) − y k |.

Conclusion
In this work, we have proposed two numerical approaches to derive the numerical methods for fractional differential equations and examined the stability and convergence of the fractional Euler method, the fractional Adams method based on newly defined fractional derivative definition. We first find the fractional Taylor series for the conformable fractional derivative. The conformable fractional Taylor series have been given at the starting point called Maclaurin fractional Taylor series. To the our best knowledge, this is is the first paper giving the generalized conformable Taylor series any other point that is different from the initial point. From this generalization, we obtain the fractional Euler and fractional Taylor methods of desired order. Frequently used Taylor methods are the Taylor methods of order 2 and order of 4, so we investigate the fractional Taylor methods of order 2α and of order 4α. We have applied these proposed methods to FDEs with known solutions. Generally, numerical methods for FDEs are based on the integral equations that are inverted from differential equations, and then one uses some techniques to approximate the integral equations to derive a numerical scheme for approximating FDEs. Most of these schemes require large number of operations and computations. However, the fractional Taylor methods need only the previous step to compute the solution at the next step. This shows that the number of operations and computations increases linearly unlike the existing methods in the literature. Thus, the proposed methods in this paper are much easier to apply and simpler to implement since they involve no special functions to evaluate, and thus the methods more accurate when compared to the existing finite difference methods for FDEs. Our second approach for numerical methods for FDEs is based on so called product integral rule. Using this technique, we derive fractional Euler and fractional Adams methods and we prove stability and error analysis for these methods. We established new inequalities in this article. We also gave some numerical examples to verify the theoretical findings. The methods can be extended to consider the numerical methods for solving the fractional partial differential equations that will be the future work.