1 Introduction
There are several definitions of fractional derivatives, of which the best known are the RiemannLiouville and Caputo definitions. Caputo’s definition (1967) of fractional derivative of order is especially wellsuited for initial value problems:
(1) 
where and is the gamma function.
It is worth mentioning that Caputo derivative was first derived by Soviet mechanic A.N.Gerasimov in 1948 [4].
The main obstacle to numerical solving FrDE is the presence of the mild singularity in its definition. Integrating the fractional derivative by parts helps to fix this problem. This method (byparts method) has been used, e.g., in [1] and [2], where the Caputo fractional derivative was approximated by finite difference scheme. Many observations demonstrate that fractional derivatives can be calculated with high precision. For benchmarking tests, like in [1] and [2], the numerical results exhibit good convergence to exact solutions. However, if the equations are slightly modified, the convergence fails quite often as we show in Section 5. In these cases, we cannot guarantee the convergence to the correct solution. In Figures 1b and 2 we show simple examples of such a divergence, which may happen even for linear equations.
The goal of this paper is to develop a reliable dual approach to find a numerical solution to fractional differential equations. As sufficiently general samples, we consider quasilinear fractional differential equations of the form
(2) 
During discretization we utilize forward, backward, and central differences with higher precision [7].
To address the convergence issues, we apply two methods of resolving the aforementioned singularity: the method of substitution and the byparts method. Both methods produce almostidentical correct results when the solution is stable and, especially, is known analytically. However, we show in Section 5 that after seemingly minor modifications to the equation, either one or both methods fail to match the analytical solutions. Our analysis leads us to conclude that the computations are reliable if both methods produce the same result. For known functions, our analysis is supported by using the Taylor series.
We demonstrate how the dual approach can be used to find solutions and evaluate their reliability for linear and nonlinear fractional differential equations (2).
In Section 2 we introduce the substitution method, in Section 3 we describe the wellknown method of integration by parts (byparts method), in Section 4 we describe the methods of discretization. Next Section 5 is critical: we provide the examples of linear equations with the byparts and substitution methods. We demonstrate that only if the results of the methods match they yield the reliable solution. As an example of the insufficiency of just one method in the computations of fractional differential equations, we reconsider an example from the wellknown monograph [[8], Section 8.4, Example 3] and make the necessary corrections. Section 6 contains the examples of quasilinear fractional equations with essential nonlinearity in the righthand side.
2 Method of substitution
Let us convert the Caputo fractional derivative into a simpler form by substituting in (1)
and , where and therefore :
This yields:
(3) 
and therefore
(4)  
In the last step we used .
The first natural wish for numerical integration is splitting the interval into equal parts. However, for each fractional derivative the step becomes different (since it depends on ). Hence, if the equation contains different derivatives, we have to deal with different grid points for each derivative, and this way is not applicable for equations with several fractional derivatives.
To fix this difficulty, we calculate at each point between and with step and introduce . Then the grid points fit to all derivatives and functions in the equation, and the use of the trapezoid rule yields
(5) 
Since the increase in leads to the decrease in , we obtain .
This approach implies that Caputo derivative can be implemented as follows:
(6) 
where and, hence, .
3 Integration by parts
In this wellknown method the singularity is eliminated by integration by parts
(7)  
where and after the elimination of singularity, the integral can be approximated by using the trapezoidal rule as follows
(8) 
where since at the last term of the trapezoidal rule disappears. Finally, in place of expression (6), for the byparts method we obtain
(9) 
As we can see from (6) and (9), the principal difference in the approximations of fractional derivatives is the use of st derivative in the byparts method whereas the highest derivatives used in the substitution method, is of the th order.
4 Discretization methods
4.1 Numerical discretization for the substitution method
If equation (2) is of the first order with , then we need to represent first derivatives, apply (5), and utilize the implicit method by setting up a system of algebraic equations. We can calibrate the number of equations for the best performance but the basic idea is the same. Our algorithm provides the second order of accuracy. The steps are as follows:

use finite differences to represent a few consecutive derivatives:
(10a) (10b) (10c) 
using formula (5) represent fractional derivative ;

multiply each derivative by the corresponding ;

sum up all expressions and add to the above generated expression the value of for each point;

equate the produced sum to for each point.
This process generates the system of algebraic equations.
It is important to point out that for the equations of the second order with at least one , in addition to the initial condition , the second initial condition is necessary , which we code using forward approximation . In this case the second derivatives need to be represented like in (11b)(11d) (in addition to (10a)(10c)). For higher derivatives their representations can be found, say, in [7].
4.2 Numerical discretization for the byparts method
The steps to solve equation like (2) using the byparts method, are the same as in the substitution method, but the fractional derivative representation changes. For first and second derivatives of the unknown function need to be used. For equations (7) and (8) FDM utilized are (points prior to point are assumed to be already found):
(11a)  
(11b)  
Remaining derivatives in the sum of equation (8):  
(11c)  
(11d) 
The use of these derivatives in formulas (7) and (8) allows us to represent fractional derivatives in quasilinear equation (2). Then we apply the steps like in substitution method aboveand generate the algebraic system of equations for . The repetition of this process through the whole interval provides the values of function on the entire interval.
4.3 Approximation errors
The above mentioned methods of discretization for both substitution and by parts methods lead to the same level of accuracy.
In fact, the best way to estimate the precision of calculation of fractional derivative is to compare it with the exact expression. If it is not available, then the Taylor series approximation can be used when the solution function is known. In this case, the Caputo derivative for the function can be expressed as follows
(12) 
Then, assuming the convergence of series, we obtain
(13) 
Although we could not find any direct reference to formula (13), we believe that it is known and is a kind of a ”folklore” in the mathematical community.
Now, we are ready to check the accuracy of both methods to calculate the fractional derivative. For example, Table 1 presents the accuracy of calculations :
Taylor exp  Substitution  Abs Err Subst  Int. by Parts  Abs Err by Parts  

0.1  0.2824821555  0.2824821407  0.2824821402  
0.2  0.4344599870  0.4344599557  0.4344599549  
0.3  0.5680457063  0.5680456557  0.5680456546  
0.4  0.6996788619  0.6996787873  0.6996787858  
0.5  0.8392329447  0.8392328384  0.8392328364  
0.6  0.9959906149  0.9959904642  0.9959904614 
It is easy to see that we have the accuracy of order as expected. Similar accuracy was observed for other smooth functions. Of course, usually we do not know the Taylor series representation of the solution to differential equations. In this case we can calculate the difference between solutions as a proxy of the error and compare it with the expected accuracy of the calculation of the fractional derivatives (in the presented schemes the accuracy is ). We resulting order of accuracy of FrDE cannot be expected better than the accuracy of the calculation of the fractional derivatives.
The dual approach works even if the equation contains integer derivatives. In this case, integer derivatives can be approximated by tiny deviation , where , and the accuracy of such the deviation of the order of derivative is almost perfect.
5 Valid and failed convergences for linear equations
The use of only one numerical method to solve fractional equations is not sufficient. Having two significantly different methods allows us to compare their outputs and, if they match, conclude that the found solution is reliable. Four examples, reflected in Figures 1 and 2, correspond to linear equation with exact solution
(14) 
In Figure 1, , (left) and (right). In case 1a equation (14) looks as follows:
(15) 
and has exact solution . Here is betafunction. In case 1b, equation (14) becomes
(16) 
and has exact solution .
In Figure 1a the results for both substitution and by parts methods coincide and provide the correct solution . The consilience of these results justifies the proposed dual approach that guarantees the reliability of computations.
In Figure 1b the results for both substitution and by parts methods differ from the exact solution
and, most important, are far from each other. Thus, the computations are not reliable.
In Figure 2 we consider slightly different equations
(17) 
with exact solution . Also, we consider equation
(18) 
with exact solution . In Figure 2a (equation (17)) the method of substitution fails whereas the byparts method is valid. In Figure 2b (equation (18)) we observe the opposite case: the method of substitution is valid but the method of integration by parts fails.
In all examples the computational step is . As we can see, even though the equations look quite similar, the results of calculations differ: from almost perfect match by both methods to complete mismatch, which demonstrates that only if both methods produce almost identical results, then the computation is reliable. This is the additional justification of the proposed dual approach.
6 Reliability for quasilinear equations
Example 1. Let’s consider equation
(19) 
with initial condition . The calculated solution is represented in Figure 3.
Table 2 shows the results of solving equation (19) by two methods. The residual is calculated as the difference between the right and left sides of the equation.
Solution found  Residual  Solution found  Residual  
using Byparts  using Substitution  
0.1  0.0061330982  0.0061330846  
0.2  0.0212821228  0.0212821387  
0.3  0.0431124645  0.0431125027  
0.4  0.0700231242  0.0700231812  
0.5  0.1007208712  0.1007209446  
0.6  0.1341161406  0.1341162285  
0.7  0.1692773586  0.1692774592  
0.8  0.2054041267  0.2054042388  
0.9  0.2418082833  0.2418084054  
1.0  0.2778991084  0.2778992392 
As we can see in Table 2, the results of both methods are almost identical and residual . However, the true error is definitely bigger than , because on top of the error of the solution we must apply the approximation error for fractional derivatives. Consequently, there is no need for us to be delighted is the calculated residual is small.
The actual error between the computations and the solution will be bigger than .
To support further the above warning on the computational errors, let us solve slightly changed equation (19) with the modified function such that it produces the known result . This value for is chosen because its graph is similar to the above computational results for (19). Since , and and , then we consider equation like (19)
(20) 
with the initial condition . Its exact solution is . In this case we can compare the calculated results by both methods with the known analytic solution:
Exact  Solution found  Solution found  Difference b/w  
using byparts  Error  by Substitution  Error  two methods  
0.1  0.01  0.0100507844  0.0100508224  
0.2  0.04  0.0400897392  0.0400897712  
0.3  0.09  0.0901232321  0.0901232606  
0.4  0.16  0.1601525607  0.1601525867  
0.5  0.25  0.2501785451  0.2501785691  
0.6  0.36  0.3602021066  0.3602021290  
0.7  0.49  0.4902245812  0.4902246024  
0.8  0.64  0.6402483126  0.6402483332  
0.9  0.81  0.8102786790  0.8102786996  
1.0  1.00  1.0003337914  1.0003338137 
As we can see in Table 3, both methods produce almost identical results and, therefore, almost the same errors. The difference between the solutions is . However, the error is about , which is, as expected, more than the precision of calculation of fractional derivatives, but what’s important is that the solutions found by both methods almost match. This indicates that it was found correctly and, consequently, is reliable.
Please note, that differential equations with can be evaluated in the same fashion with similar precision. For example, for equation
(21) 
with initial conditions both methods produce almost identical results.
Example 2. For equations containing fractional derivatives of different orders (between and ), the results may be less reliable. Having two different methods of calculation provides us with a path to analyzing the validity of the results.
If we solve semilinear equation
(22) 
then our two methods produce quite different outputs, which indicate that the found by either method solution is not reliable (see Figure 4).
If we slightly modify the entries in (22) the computation becomes more stable. For example, by changing to , we arrive at equation
(23) 
The solution for equation (23) is presented in Figure 5 using both methods.
For further analysis, we replace terms in (23) by similar function, which provides exact solution . Thus, we consider the following equation:
(24) 
Numerical solution of equation (24) produces with error by both almost coinciding methods. This confirms that solution to (23) is reliable, even though close equation (22) suggests opposite. These examples justify once again the necessity of the proposed dual approach.
7 Conclusion
To analyze the computational methods for fractional differential equations, we introduce the substitution numerical method and, along with wellknown byparts method, suggest the dual approach for the reliability of computations for linear and quasilinear fractional differential equations. We show numerically that in the cases when both methods produce almost identical outputs, the solution is reliable, and its level of accuracy is close to the accuracy of the approximations of fractional derivatives. As a tool to verify the accuracy for fractional derivatives, we use the Taylor series expansion. We demonstrate the validity of the dual approach by solving linear and nonlinear fractional differential equations. Our computational examples and counterexamples demonstrate the necessity of the dual approach.
References
 [1] R. Albadarneh, M. Zerqat, I. Batiha. Numerical Solutions for linear and nonlinear Fractional Differential Equations. International Journal of Pure and Applied Mathematics, 106 No. 3, 2016, pp. 859–871.
 [2] R. Albadarneh, M. Zurigat, I. Batiha. Numerical Solutions for linear Fractional Differential Equations of order using Finite Difference Method (FFDM). J. Math. Computer Sci. 16 (2016), 103–111.
 [3] S.P. Mirevski, L. Boyadjiev, R. Scherer. On the RiemannLouisville fractional calculus, gJacobi functions and FGauss functions. Applied Mathematics and Computation 187, No 1, 2007, pp.315–325
 [4] A.N.Gerasimov. A generalization of the deformation laws and its application to the problems of internal friction. Applied Mathematics and Mechanics 12, 1948, pp. 251–260.
 [5] Z. Tomovski, R. Hilfer, H.M. Srivastava. Fractional and operational calculus with generalized fractional derivative operators and Mittag–Leffler type functions. Integral Transforms and Special Functions, 2010, 21(11), pp.797814
 [6] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo. Theory and applications of Fractional Differential equations. Fakulteit der Exacte Wetenschappen, 2006, Amsterdam, The Netherlands.
 [7] J. Mathews, K. Fink. Numerical Methods Using MATLAB. Fourth Edition.
 [8] I. Podlubny. Fractional Differential Equations. Mathematics in Science and Engineering, Vol. 198, 1999.
 [9] N.V.Zhukovskaya, A.A.Kilbas. Soilving homogeneous fractional differential equations of Euler type. Differential Equations 47, No. 12, 2011, pp. 1714–1725.
Comments
There are no comments yet.