Dual Approach as Empirical Reliability for Fractional Differential Equations

by   P. B. Dubovski, et al.

Computational methods for fractional differential equations exhibit essential instability. Even a minor modification of the coefficients or other entry data may switch good results to the divergent. The goal of this paper is to suggest the reliable dual approach which fixes this inconsistency. We suggest to use two parallel methods based on the transformation of fractional derivatives through integration by parts or by means of substitution. We introduce the method of substitution and choose the proper discretization scheme that fits the grid points for the by-parts method. The solution is reliable only if both methods produce the same results. As an additional control tool, the Taylor series expansion allows to estimate the approximation errors for fractional derivatives. In order to demonstrate the proposed dual approach, we apply it to linear, quasilinear and semilinear equations and obtain very good precision of the results. The provided examples and counterexamples support the necessity to use the dual approach because either method, used separately, may produce incorrect results. The order of the exactness is close to the exactness of fractional derivatives approximations.



There are no comments yet.


page 1

page 2

page 3

page 4


Substitution Method for Fractional Differential Equations

Numerical solving differential equations with fractional derivatives req...

Shifted Monic Ultraspherical Approximation for solving some of Fractional Orders Differential Equations

The purpose of this paper is to show and explain a new formula that indi...

Fractional integrals, derivatives and integral equations with weighted Takagi-Landsberg functions

In this paper we find fractional Riemann-Liouville derivatives for the T...

Estimating Discretization Error with Preset Orders of Accuracy and Fractional Refinement Ratios

In solution verification, the primary goal is finding an accurate and re...

A Rotating-Grid Upwind Fast Sweeping Scheme for a Class of Hamilton-Jacobi Equations

We present a fast sweeping method for a class of Hamilton-Jacobi equatio...

A modified EM method and its fast implementation for multi-term Riemann-Liouville stochastic fractional differential equations

In this paper, a modified Euler-Maruyama (EM) method is constructed for ...

Fractional-Step Runge–Kutta Methods: Representation and Linear Stability Analysis

Fractional-step methods are a popular and powerful divide-and-conquer ap...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

There are several definitions of fractional derivatives, of which the best known are the Riemann-Liouville and Caputo definitions. Caputo’s definition (1967) of fractional derivative of order is especially well-suited for initial value problems:


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 (by-parts 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


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 by-parts method. Both methods produce almost-identical 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 well-known method of integration by parts (by-parts method), in Section 4 we describe the methods of discretization. Next Section 5 is critical: we provide the examples of linear equations with the by-parts 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 well-known 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 right-hand 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:


and therefore


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


Since the increase in leads to the decrease in , we obtain .
This approach implies that Caputo derivative can be implemented as follows:


where and, hence, .

3 Integration by parts

In this well-known method the singularity is eliminated by integration by parts


where and after the elimination of singularity, the integral can be approximated by using the trapezoidal rule as follows


where since at the last term of the trapezoidal rule disappears. Finally, in place of expression (6), for the by-parts method we obtain


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 by-parts 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:

  1. use finite differences to represent a few consecutive derivatives:

  2. using formula (5) represent fractional derivative ;

  3. multiply each derivative by the corresponding ;

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

  5. 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 by-parts method
The steps to solve equation like (2) using the by-parts 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):

Remaining derivatives in the sum of equation (8):

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


Then, assuming the convergence of series, we obtain


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
Table 1: Calculation of fractional derivatives with step .

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


In Figure 1, , (left) and (right). In case 1a equation (14) looks as follows:


and has exact solution . Here is beta-function. In case 1b, equation (14) becomes


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.

Figure 1: 1a (left), equation (15): both substitution and integration by parts methods match each other and the exact solution. 1b (right), equation (16): both substitution and by-parts methods diverge and don’t match neither each other nor the exact solution.

In Figure 2 we consider slightly different equations


with exact solution . Also, we consider equation


with exact solution . In Figure 2a (equation (17)) the method of substitution fails whereas the by-parts 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.

Figure 2: 2a (left), equation (17): by-parts method matches the solution, but substitution fails. 2b (right), equation (18): by-parts method fails, but the substitution method is valid.

6 Reliability for quasilinear equations

Example 1. Let’s consider equation


with initial condition . The calculated solution is represented in Figure 3.

Figure 3: Solution to equation (19) with step obtained by both methods, which provide almost identical results. The solution is reliable.

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 By-parts 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
Table 2: Solution to equation (19). Step .

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)


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 by-parts 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
Table 3: Solution to equation (20) found by each method with step .

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


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


then our two methods produce quite different outputs, which indicate that the found by either method solution is not reliable (see Figure 4).

Figure 4: ’Solutions’ to equation (22) with step . Red line – substitution method, blue – by-parts method. The solutions are not reliable.

If we slightly modify the entries in (22) the computation becomes more stable. For example, by changing to , we arrive at equation


The solution for equation (23) is presented in Figure 5 using both methods.

Figure 5: Solution for (23) with step . The methods produce almost identical result, which is reliable.

For further analysis, we replace terms in (23) by similar function, which provides exact solution . Thus, we consider the following equation:


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 well-known by-parts 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.


  • [1] R. Albadarneh, M. Zerqat, I. Batiha. Numerical Solutions for linear and non-linear 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 Riemann-Louisville fractional calculus, g-Jacobi functions and F-Gauss 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.797-814
  • [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.