In the paper there are presented and evaluated for effectiveness three methods of accuracy increase of fractional order derivatives and integrals computations for application with the Riemann-Liouville/Caputo formulas. They are based on the ideas of either transforming difficult integrand in the formulas to high-accuracy computations requirements of a applied method of numerical integration or adapting a numerical method of integration to handle with high-accuracy a difficult feature in the integrand. Additional accuracy gain is obtained by incorporating increased precision into computations. The actual accuracy improvement by applying presented methods is compared with the capabilities of wide range of available methods of integration.
Keywords and phrases
- Numerical Integration
- Fractional Order Derivatives and Integrals
Fractional calculus may be described as an extension of the concept of a derivative operator from integer order
Beginnings of fractional calculus indeed date to the same time as classical calculus. However, the question of consequences of “taking square root of 1
The most rapid interest grow of this new topic can be observed in the last 40 years. The first conference “International Conference of Fractional Calculus and Its Applications” entirely dedicated to it took place in 1974 at the University of New Heaven, Connecticut, USA . The first monograph dedicated to fractional calculus was published the same year .
From that time on, the topic won interest of scientists from various areas of mathematics, physics, chemistry, electrical engineering, economy, biological sciences. Initial works had a pure theoretical character. However, lately, can be noticed intensification of research on practical application of this tool in physics, technical, biological and economical sciences. In technical sciences, application of fractional calculus is mostly applied for areas of electrical engineering, electronics and control systems as well as image analysis and processing.
Fractional calculus requires extensive application of numerical methods, e.g. numerical integration. However, commonly applied numerical methods are not suitable for that purpose in their existing form.
Effectiveness of fractional derivatives and integrals practical applications can be increased by improving their computational accuracy. Thus, increasing computational accuracy of fractional order derivatives and integrals became a very important task to study for many computer scientists. Increasing it as high as possible and by using as little as possible hardware resources and time is a matter of aptly selected method, applied hardware, programming language and programming techniques in the process.
Growing popularity of computer as a tool for scientific research led to essential advancements of existing programming languages. C++ can be considered as best example in this context. C++ is a general purpose programming language, which enable hardware speed optimization. It is a flexible language with great versatility and it has huge function library which can be used for developing system software, e.g. operating systems, compilers, editors and data bases.
To be able to solve difficult numerical problem according to a set goal, a scientist has to make some crucial decisions regarding applied for that purpose hardware, programming tools and techniques. They include appropriate selection of computer programming languages because their abilities to process numerical data determine how accurate the solution will be.
The selection of uniform C++ equipped with standard mathematical library as a main programming tool is not enough nowadays to solve numerical problems with satisfactory accuracy.
Double precision arithmetic commonly applied in scientific numerical calculations has many flaws which influences negatively the accuracy of computations, e.g. limitations of number values which double precision variables can hold or no influence of programmer on mathematical operations rounding.
Arbitrary precision mathematical library GNU MPFR , built upon GNU GMP library  application enable any computation to be performed with any precision limited only by hardware, with either real or complex numbers. The library supports proper rounding of all mathematical operations as wel, which is crucial because most of inaccuracies, errors and limitations in scientific computations is caused by finite precision arithmetic and a lack of proper rounding of mathematical operations. The library has also at their disposal many ready to use functions, for example gamma function and reciprocal gamma function, which are extensively used for fractional order derivatives and integrals computations.
Fractional derivative or integral is determined in [
The following definitions of fractional order differentiation and integration are the scope of the following paper:
in which Caputo (C) fractional derivative 
Caputo (C) fractional derivative 
Fractional order derivatives and integrals computational formulas (3), (5) and (6) consist of ordinary integrals. However, the integrals are difficult to compute due to presence of singularity at the end of integration interval preceded by rapid values increase (stiffness).
Each integrand included in the formulas can be divided into a kernel and a function, i.e.
Plot of the kernel of the integrand is presented in Fig. 1. Gray area denotes integral to compute.
Stiffness of the integrand can increase or decrease according to a value of power and a character of a function.
This characteristics of the kernel determines accuracy of computations.
To compute fractional order derivatives and integrals with help of formulas (3), (5) and (6) there can be applied commonly used methods of numerical integration, e.g. Midpoint Rule, Trapezoidal Rule as well as Gauss-Legendre and Gauss-Kronrod Quadratures (more informations on the methods can be found in  and ).
Due to the singularity problem in the integrand there can be applied some techniques advised in the literature , , , : (i) Increasing of subintervals amount, (ii) Exclusion of singularity from integration, (iii)
Unfortunately excluding the singularity from the integration decreases in this case accuracy even more because area which contributes to integral value at most proceeds it directly.
Stretching the difficult section even many times and using an appropriately increased amount of subintervals does not effect accuracy satisfyingly. The end section of the integrand is too stiff for the listed methods of numerical integration.
Reciprocal function integration can be applied to simple, one-valued functions only. Therefore, it may not be applied as a general solution of the problem.
Another solution can be an application of numerical methods which have been developed or modified particularly for application with (3), (5) and (6). There exist only a few such methods. They are based on: (i) Extrapolation , (ii) Trapezoidal Rule , (iii) Modified Trapezoidal Rule for Caputo derivatives only , (iii) Mollification techniques , (iv) Chebyshev polynomials (for fractional integrals only) .
However, the only one algorithm which can be considered as significant is the one based on the Trapezoidal Rule.
Diethelm removed analytically the singularity from the integrand by applying integrating by parts and calculating some quadrature coefficients. Resulting accuracy of fractional order derivatives and integrals computations can be additionally improved by application of the Richardson Extrapolation.
Diethelm’s method has been redesigned by Odibat for application for Caputo fractional derivatives only . His modifications are of minor importance and did not increased overall accuracy capabilities of the method.
Maximum of available computational accuracy of fractional order derivatives and integrals calculations applying formulas (5) and (6) and presently available methods of numerical integration is presented in the tables 1-3.
GL NCm Diet Odiba 8 8.11 10.68 0.0004 0.0024 15 4.25 7.81 0.0004 0.0023 21 3.02 6.6 0.0004 0.0023 61 1.03 3.88 0.0004 0.0022 300 0.21 1.75 0.0002 0.0021 600 0.11 1.24 0.0002 0.0021 1000 0.06 0.96 0.0002 0.0021 GL NCm 8 61.87 13.52 15 29.32 5.88 21 20.17 4.12 61 6.54 1.81 300 1.3 0.76 600 0.65 0.53 1000 0.39 0.41 GL NCm Diet Odiba 8 130.03 45.42 17.67 5.78 15 66.14 33.64 7.95 3.46 21 46.09 28.77 5.18 2.94 61 16.1 17.84 1.53 2.94 300 5.02 9.35 0.62 2.43 600 3.67 7.29 0.56 2.35 1000 3.14 6.18 0.55 2.35
Applied methods include the Grünwald-Letnikov method (GL) (7), Newton-Cotes Midpoint Rule (NCm) and modified Trapezoidal Rule by Diethelm (Diet) and Odibat (Odiba). There are not included any results for Gauss Quadratures because they are not suitable for this kind of calculations (relative error exceeding 90%). The same applies to Diethelm’s and Odibat’s method applied for fractional derivative of the exponential function.
The programs for this part has been created by the author using algorithms described in the freely available papers and literature of the subject.
Computational accuracy measure assumed in the following section is percentage relative error. Computer programs are written using C++ with double precision and complied applying GNU GCC compiler.
Computer configuration include Intel i7 processor, 8GB of RAM and Linux Ubuntu 14.04 LTS 64-bit operating system.
The method by Diethelm and its modification by Odibat enable to calculate fractional derivatives and integrals of some orders and of some functions with accuracy limited at best to 3-4 decimal places.
Therefore, the conclusion is that none of the presently available methods can be considered as a method of general use which assures computations of fractional derivatives and integrals with high accuracy independently from order, interval and type of the function.
Conducted survey suggests that to calculate fractional order derivatives and integrals with higher accuracy than presently is possible either an applied method of numerical integration must be adapted to handle with high-accuracy an integrand with the difficult feature or an integrand must be transformed to fit into a numerical method’s of integration requirements for high accuracy of computations.
In the next three subsections there are presented various approaches to the subject of accuracy increase of fractional order derivatives and integrals computations for application with the Riemann-Liouville and Caputo formulas.
A method of transforming an integrand is analytical integrand transformation via independent variable substitution. It has to be conducted prior numerical integration.
Its rules are well known from analytical integration. However, this time instead of transforming an integrand into one which is analytically easier integrable, the goal is to obtain an integrand with a shape which fits into requirements of a selected method of numerical integration for high accuracy of computations.
Application of variable change
Application of variable change
For an example function
Transformation (9) is now applied for practical accuracy increase using two periodic functions:
and for the function (12) assumes the form
The accuracy of computations are presented in Figs 3(a) and 3(b). Due to similar performance of the methods for all three evaluated orders
For accuracy increase evaluation there are applied well known quadratures. They are considered as highly accurate and efficient if an integrand fulfills their requirements. The methods include: Newton-Cotes Midpoint Rule (NCm), Gauss-Legendre Quadrature (GLeg), Gauss-Kronrod Quadrature (GKr), Gauss-Laguerre Quadrature (GLag), Gauss-Jacobi Quadrature (GJ).
The Grünwald-Letnikov formula (GL) (7) is applied for comparison purposes. Although it operates strictly on a integrated function’s values and can not be tested on transformed integrand, its accuracy can be considered as a point of reference.
Exact values required for error computation are calculated using analytical formulas to be found in ,  and . If a formula is not to be found, there is applied a novel approach to the accuracy assessment of fractional order derivatives and integrals using generalized fractional operators concatenation rules .
In the plots with results, accuracy of integration applying unchanged integrands are denoted with prefix RL, e.g. RL GKr and with modifications applied to the integrand - with prefix mRL, e.g. mRL GKr.
and for the function (12) assumes the form
Analytical independent variable substitution in the integrand is a very effective method of accuracy increase. An average four times accuracy increase of the transformed integrand integration proves it. Application of the technique improves efficiency as well; the increased accuracy is obtained by using half of the sampling points required by the same numerical integration method applied to the unchanged integrand.
How high the accuracy increase actually is depends on the shape of the transformed integrand and how much it agrees with an “optimal” shape for high accuracy integration by applying particular method of integration.
Nevertheless, “manual” transformation is arduous and time consuming, which can hinder some applications.
Transformation of an integrand via independent variable substitution can also be applied in form of a numerical quadrature. Such approach enables to omit arduous analytical transformation of each integrand and automate the whole transformation process.
An efficient and universal method of an integrand transformation is known as the Double Exponential Transformation (DE). The DE Transformation is a joint of hyperbolic function independent variable substitution in the integrand and the Trapezoidal Rule.
Such transformed form will be referred as the Double Exponential Quadrature (the DE Quadrature).
The general idea of this kind of variable transformation was proposed for the first time by Korbov . It also became a subject of works of Schwartz , Stenger , Takahasi , Mori , Haber  and others, who gave it different names.
The DE transformation never achieved a wider popularity due to poor publicity and some drawbacks. The drawbacks are mainly associated with the limitations of the single and double precision arithmetic applied in the initial programming at the time of development, e.g. limited range of applied single and double precision variables which can cause under- and over-flow when calculating weights required for the quadrature.
Unfortunately, application of the solutions decrease overall accuracy abilities of the method.
The problem can also be solved by applying arbitrary precision variables and arbitrary precision mathematical library GNU MPFR, which improves accuracy of the DE Transformation.
The DE Transformation implemented for the purpose of fractional orders derivatives and integrands computations is based on the transformation, which was proposed by Schwartz  as the Tanh Rule. It can be described as follows:
Let us consider the integral
First we apply the following variable transformation
It is important that
Due to truncation of the summation process in (17) at some arbitrary chosen
Two kinds of errors should be taken into consideration when implementing the DE Formula: discretization error, due to the use of the Trapezoidal Rule to approximate an integral and truncation error, because of truncation of infinite sum at some
The subinterval width
The following equations (18)-(20) present three different transformations with varied convergence rate , . Transformation (18) has the slowest convergence and transformation (20), the fastest respectively:
are the nodes, the sampling points and the weights of the DE Quadrature respectively.
and asymptotic error in terms of the number
Calculation of the DE Quadrature weights from the formula (21) can cross the double precision variables limits. First, the denominator in the weight function formula sec2 (
Increasing value of the parameter
Another parameter which influence accuracy of the DE Quadrature (21) is an amount N of subintervals for integration applying the Trapezoidal Rule.
For accuracy assessment of the DE Transformation there are selected three functions:
and four fractional orders:
Fractional derivative of order: Fractional integral of order:
Fractional derivative of order:
Fractional integral of order:
Independent variable substitution in integrand in a form of a numerical quadrature has even more advantages than its “manual” counterpart presented in the previous subsection: it can be applied automatically during the integration, the nodes and the weights of the DE Quadrature can be computed only once and applied many times. Their computation is not difficult.
How high the accuracy increase actually can be depends on the shape of the transformed integrand and how much it agrees with an “optimal” shape for high accuracy integration by applying the Trapezoidal Rule and how “far” left and right from the
The accuracy increase up to ∼ 10−50 can be obtained by careful programming applying variables with increased precision. Still, some orders are out of reach as for example the order
The reason for the unsatisfactory accuracy for the orders near 0 and 1 is the decreasing value of power in the denominator in the Riemann-Liouville/Caputo formulas and limitations of the applied substitution expression.
Using the DE Transformation with the double precision, the accuracy is limited to two decimal places.
Despite the obvious advantages of arbitrary precision application, the user must be aware that they extend immensely time required for computations if amount of precision digits applied in calculations is greater that 100 (see Fig. 8).
There exists a wide range of numerical methods for solutions of singularities and infinite limits issues in numerical integration, e.g. the Gauss Quadratures. They involve application of polynomials in an integral approximation. The polynomials’ properties are selected for the specific integration problem and thereby they assure high accuracy results. The problematic area for integration is computed differently, which assures high accuracy.
Generally, the Gauss Quadratures can be applied to a specific types of integrands only. It is conditioned by the linkage to their weight functions (with exclusion of the Gauss-Legendre Quadrature with the weight function equals 1).
The Gauss Quadratures are commonly applied to compute most difficult integrals, e.g. the improper integrals and integrals with end point singularities.
The ordinarily sufficient method of the Gauss Quadratures application is reduced to the use of tabulated values of nodes and weights. They are freely available ,  and easy to incorporate into a computer program. However, the integration process can become difficult if nodes and weights must be computed, an integrand does not meet exactly or at all requirements of an quadrature’s weight function or an integration interval is different than assumed.
The first case requires the application of mathematical formulas for nodes and weights.
The second issue solution involves manual extraction of a weight function from an integrand or a weight function adjustment to a problematic integrand.
Changing integration interval from default to arbitrary selected involves a conversion if a quadrature’s requirements permit it.
Next there are presented some important details to the Gauss-Jacobi Quadrature numerical application requirements as well as the formulas for the nodes and the weights computations.
A weight function which eliminates definite integration range endpoints singularities is Jacobi weight function
General, well known formula of approximate calculation of definite integral
with the weight function (26) assumes the form:
The Jacobi Polynomial can be determined by applying Rodrigues formula
is the remainder of integral approximation.
and relation 
Transforming the interval [
To calculate fractional order derivatives we proceed the similar way
The Riemann-Liouville fractional order derivative formula (5) assumes the following form
The Gauss-Jacobi Quadrature (27) requires computation of polynomial of order
Existing approaches to nodes and weights computing, some of which have been widely used for many years, suffer from 𝒪(
Accuracy of the method is conducted in the similar conditions as in the case of DE Transformation. First, there is assessed accuracy for double precision (denoted as 16-digits precision) and high precision calculations. Wherein, it is only required to 100-digit variables increase in high precision calculations for GJ to achieve the highest possible accuracy for the method. Hence the lack of results for other precision options present in the DE Transformation accuracy assessment. Results are presented in Figs. 9 and 10.
The Gauss Jacobi Quadrature is applied for calculations of fractional order integrals
Next, GJ is evaluated with
Now, GJ is evaluated for accuracy of calculations for different fractional orders 0 <
Fig. 14 presents comparison of rough time complexity of calculations for GJ, DE Transformation (DE) and Grünwald-Letnikov method (GL) applied with double and 100-digits precision.
Approximation of fractional order derivatives and integrals applying the Gauss-Jacobi Quadrature is the most effective method of accuracy increase from all presented in the paper. The integrand is left intact and instead, the weights are changed according to its properties; the computation of the difficult kernel in the Riemann-Liuoville/Caputo formulas is removed from the integration (it is now calculated using the weight function). Programming the quadrature applying variables with increased precision ascertain constant high accuracy of computations, independently from integrated function, fractional order and interval. Wherein difficult fractional orders of integrals near 0 and fractional orders of derivatives near 1, which are out of reach for any other methods can also be calculated with similar high accuracy.
The degree N of polynomial applied in integral approximation and an amount of digits of precision used for calculation are the only factors which influence calculation accuracy.
Another advantage of the method is that application of high order polynomial
The problems associated with singularities are the most difficult to solve in numerical integration. However, application of a proper transformation to the integrand in the Riemann-Liouville/Caputo formulas or a correct selection of parameters for the weight function of the Gauss-Jacobi Quadrature make possible to calculate fractional order derivatives and integrals with high-accuracy. The accuracy gain can be improved by application of arbitrary precision in computations.
Analytical integrand transformation applying appropriate in context of a selected method of integration substitution expression is very effective and can increase accuracy of numerical integration significantly. However, the method is arduous because it requires manual transformation of each integrand.
Transformation of the integrand via independent variable substitution applied in a form of a numerical quadrature enables to omit arduous analytical transformation of an integrand and automate the whole substitution process.
Application of the DE Quadrature turned the integration process into automatic one. It is effective under the condition of finding and applying substitution expression which transforms an integrand into one which matches optimal requirements of the Trapezoidal Rule.
Approximation of fractional order derivatives and integrals by applying the Gauss-Jacobi Quadrature is the most effective method which assures high-accuracy results up to ∼ 10−120 mark independently from fractional order, integration interval and a type of a function.
D(1/2)f(t), f(t) = t, t ∈ (0,1), relative error in %
D(1/2)f(t), f(t) = e−t, t ∈ (0,5), relative error in %
D(1/2)f(t), f(t) = sin(t), t ∈ (0,2π), relative error in %