Comparison of Fractional Order Derivatives Computational Accuracy - Right Hand vs Left Hand Definition

Fractional calculus involves among the others application of fractional order derivatives (FOD) and fractional order integrals (FOI), which order can assume a real or a complex value. Thus derivatives and integrals of integer order can be considered as one of special cases of fractional calculus.

Fractional derivatives can be calculated using numerous approaches which include Riemann-Liouville fractional derivative (RL) and Caputo fractional derivative (C) [1]. The latter definition is often preferred for simulations of real world processes and systems in a form of fractional differential equations (FDE) [2, 3] due to its ability to accept initial conditions in terms of integer order values. Riemann-Liouville fractional derivative definition does not have this practical application advantage and thus is mostly applied for theoretical considerations. However, regarding the subject of the paper, the substantial difference between RL and C definitions lies in the method of their computation: RL method involves repeated integration of a function and then differentiation of the resulting function with an integer order; C method attempts to arrive at the same result using the same operations, but in the reverse order. Thus using direction terms we can call RL and C left-hand and right-hand definitions, respectively.

The aim of the paper is to compare computational accuracy of fractional derivatives by applying the left-hand and the right hand definitions.

To fulfill this goal we have divided the paper in the following subsections: first there are presented formulas for fractional derivatives computation. The second subsection explains the terms left-hand and right-hand. The next one contains brief description of problems associated with FOD/FOI computation and presents numerous numerical methods which are applied for fractional derivatives computation. Difficulties associated with FOD/FOI application include computation of reference values required for accuracy assessment. Thus, the next section presents formulas and methods for their computation. The section after that one is devoted to computational accuracy comparison of the right- and left-hand definitions. It includes detailed results for computation of FOD for some functions in various ranges and short comment about it. The paper ends with conclusion which summarize the reasons for applying a particular definition for FOD computation.

The following formulas are applied for FOD/FOI computations. To differentiate between fractional and integer orders, fractional orders are commonly denoted by the Greek letters ν and μ and integer orders by the letters n and m.

Riemann-Liouville fractional integral for ν > 0 is defined as

0RLJt(ν)f(t)=1Γ(ν)∫0tf(τ)(t−τ)1−νdτ.$$\begin{array}{} \displaystyle _{0}^{RL}\textrm{J}_{t}^{\left ( \nu \right )}f\left ( t \right )=\frac{1}{\Gamma\left ( \nu \right )}\int_{0}^{t}\frac{f\left ( \tau \right )}{\left ( t-\tau \right )^{1-\nu}}d\tau. \end{array}$$(1)

Riemann-Liouville fractional derivative for ν > 0, μ > 0 is defined as

0RLDt(ν)f(t)=DnJνf(t)=dndtn[1Γ(n−ν)∫0tf(τ)(t−τ)ν−n+1dτ]$$\begin{array}{} \displaystyle _{0}^{RL}\textrm{D}_{t}^{\left ( \nu \right )}f\left ( t \right )=D^{n} J^{\nu} f\left(t\right)=\frac{d^{n}}{dt^{n}}\left[ \frac{1}{\Gamma\left ( n-\nu \right ) } \int_{0}^{t}\frac{f\left ( \tau \right )}{\left ( t-\tau \right )^{\nu-n+1}}\;d\tau \right] \end{array}$$(2)

and

0RLDt(μ)f(t)=DnJνf(t)=dndtn[1Γ(ν)∫0tf(τ)(t−τ)ν−1dτ],μ=n−ν$$\begin{array}{} \displaystyle _{0}^{RL}\textrm{D}_{t}^{\left ( \mu \right )}f\left ( t \right )=D^{n} J^{\nu} f\left(t\right)=\frac{d^{n}}{dt^{n}}\left[ \frac{1}{\Gamma\left ( \nu \right ) }\int_{0}^{t}\frac{f\left ( \tau \right )}{\left ( t-\tau \right )^{\nu-1}}\;d\tau \right],\; \mu=n-\nu \end{array}$$(3)

Caputo fractional derivative

0CDt(ν)f(t)=1Γ(n−ν)∫0tf(n)(τ)(t−τ)ν−n+1dτ,$$\begin{array}{} \displaystyle _{0}^{C}\textrm{D}_{t}^{\left ( \nu \right )}f\left ( t \right )= \frac{1}{\Gamma\left ( n-\nu \right )}\int_{0}^{t}\frac{f^{\left ( n \right )}\left ( \tau \right )}{\left ( t-\tau \right )^{\nu-n+1}}\;d\tau, \end{array}$$(4)

with conditions f (t) = 0 for t ≤ 0, f (0) = 0, f⁽¹⁾ = f ⁽²⁾ ... f⁽ⁿ⁾ = 0.

formulas (2) or (3) and (4) are related by

0RLDt(ν)f(t)=0CDt(ν)f(t)+∑k=0n−1tk−νΓ(k−ν+1)f(k)(0).$$\begin{array}{} \displaystyle _{0}^{RL}D_{t}^{\left ( \nu \right )}f\left ( t \right )=_{0}^{C}D_{t}^{\left ( \nu \right )}f\left ( t \right )+\sum_{k=0}^{n-1}\frac{t^{k-\nu}}{\Gamma\left ( k-\nu+1 \right )}f^{\left(k\right)}\left ( 0 \right ). \end{array}$$(5)

Placing C formula (4) in the right side of equation (5) enables obtaining an equivalent Riemann-Liouville fractional derivative formula [4].

0RLDt(ν)f(t)=∑k=0n−1tk−νf(k)(0)Γ(k−ν+1)+1Γ(n−ν)∫0tf(n)(τ)(t−τ)ν−n+1dτ.$$\begin{array}{} \displaystyle _{0}^{RL}\textrm{D}_{t}^{\left ( \nu \right )}f\left ( t \right ) =\sum_{k=0}^{n-1}\frac{t^{k-\nu} f^{\left( k \right)}\left ( 0 \right )}{\Gamma\left ( k-\nu+1 \right )} +\frac{1}{\Gamma\left ( n-\nu \right )}\int_{0}^{t}\frac{f^{\left ( n \right )}\left ( \tau \right )}{\left ( t-\tau \right )^{\nu-n+1}}\;d\tau. \end{array}$$(6)

In formulas (1)-(6) ν is a real number such that n − 1 < ν < n, n denotes an integer number n = ⌈ν⌉ and Γ(.) is Euler’s Gamma function.

As it was stated in introduction, the reason for selecting such a naming schema for RL and C formulas is due to the method of their computation: first we select order ν ε R₊. Then we calculate order n such that n − 1 < ν < n.

Having selected both numbers, we can calculate fractional derivative using left-hand definition presented in figure 1 and formulas (2) and (3). Here we first integrate a function with order n − ν and then we calculate n^th derivative of the resulting function.

Fig. 1

By applying right-hand definition presented in figure 2 and formula (4) and (6) we attempt to arrive at the same results in the reverse order: first we calculate n^th order derivative of a function and then we integrate the resulting function with order n − ν.

Fig. 2

Riemann-Liouville fractional order derivative and integral and Caputo fractional order derivative formulas (1)-(6) consist of an ordinary integral. However, the integral is difficult to compute due to presence of a singularity at the end of integration interval that is preceded by rapid values increase. Similarly to the term used for such occurrence in the world of differential equations, it will be called the stiffness.

The stiffness can increase or decrease correspondingly to the values of the exponent in the integrand in formulas (1)-(6) and characteristics of a function. However, it is the stiffness and the singularity presence that determine the accuracy of integral computation.

To compute fractional order derivatives and integrals by applying formulas (1)-(6) there can be applied commonly used methods of numerical integration, e.g. midpoint rule, trapezoidal rule, gauss quadratures of a general use, etc. However, their accuracy and effectiveness are heavily decreased by the presence of the singularity and the stiffness. Additionally, computation of FOD/FOI of the fractional orders near integer orders (e.g. 0 and 1) of any function is usually burdened with several hundred percent relative error. Thus we can say that none of commonly applied methods of integration can be considered as methods of general use for the purpose of computing FOD/FOI with high accuracy independently of order and interval.

To increase computational accuracy of FOD/FOI, either an applied method of numerical integration must be adapted to handle a difficult integrand properly or an integrand must be transformed to fit into the requirements of an applied numerical method of integration. Furthermore replacing standard double precision computer arithmetic with arbitrary precision for numerical integration enables increasing many times the overall accuracy of computation.

In our papers [5, 6] we presented two-way approach for computation accuracy increase of FOD/FOI, which included integrand transformation and adaptation of Gauss-Jacobi quadrature weight function with significant additional accuracy increase obtained by applying the GNU GMP [7] and the GNU MPFR [8] arbitrary precision mathematical libraries with C++ MPFR C++ interface [9] for programming.

For the reason of its high-accuracy abilities, efficiency and versatility, Gauss-Jacobi quadrature with adapted weight function is applied for the purpose of accuracy comparison of the FOD/FOI right- and left-hand definitions.

Application of Gauss quadratures enables solving efficiently problems associated with singularities and infinite limits by approximating an integrand and its difficult features by using polynomials. Wherein, particular polynomial is selected according to its properties to address a specific problem encountered in an integrand. If that is the case, the application of one of the Gauss quadratures assures high accuracy integration with a few sampling points.

Application of a Gauss quadrature (except for Gauss-Legendre quadrature) enables removing a problematic area or feature of an integrand from the integration and computing it by applying the formulas for the weights of a quadrature. However, Gauss quadratures can only be applied to a specific types of integrands, which is conditioned by the linkage to their weight functions.

By applying Gauss-Jacobi quadrature with the weight function p(x)

p(x)≡(1−x)λ(1+x)β,λ,β>−1,$$\begin{array}{} \displaystyle p\left ( x \right )\equiv\left ( 1-x \right )^{\lambda}(1+x)^{\beta},\;\lambda,\beta>-1, \end{array}$$(7)

general formula of an approximate calculation of definite integral

∫abp(x)f(x)dx≅∑k=1nAkf(xk)+R.$$\begin{array}{} \displaystyle \int_{a}^{b} p\left(x\right)f\left(x\right) dx \cong \sum_{k=1}^{n} A_{k} f\left(x_{k}\right)+R. \end{array}$$

can be expressed as

∫−11(1−x)λ(1+x)βf(x)dx≈∑k=1nAkf(xk)+R,$$\begin{array}{} \displaystyle \int_{-1}^{1}\left ( 1-x \right )^{\lambda }(1+x)^{\beta} f\left ( x \right )dx \approx \sum_{k=1}^{n}A_{k} f\left ( x_{k} \right)+R, \end{array}$$(8)

in which x_k are zeros of Jacobi polynomial Jn(λ,β)(xk)$\begin{array}{} \displaystyle J_n^{(\lambda ,\beta )}({x_k}) \end{array}$ of degree n and R is an error.

We can use this property for FOD/FOI computation.

Removing kernel from an integrand

0RLDt(ν)f(t)=1Γ(ν)∫0t(t−τ)ν−1︸ kernelf(τ)dτ,$$\begin{array}{} \displaystyle _{0}^{RL}\textrm{D}_{t}^{\left ( \nu \right )}f\left ( t \right )=\frac{1}{\Gamma\left ( \nu \right )}\int_{0}^{t} \underbrace{\left ( t-\tau \right )^{\nu-1}}_\text{ kernel}\; f\left ( \tau \right ) d\tau, \end{array}$$

substituting λ = ν − 1, β = 0, we obtain

∫−11(1−t)λf(t)dt≈∑k=1nAkf(tk)=∑k=1n2v(1−tk2)[Jn(λ,0)'tk]2f(tk)$$\begin{array}{} \displaystyle \int_{ - 1}^1 {{{\left( {1 - t} \right)}^\lambda }} f\left( t \right)dt \approx \sum\limits_{k = 1}^n {{A_k}} f\left( {{t_k}} \right) = \sum\limits_{k = 1}^n {{{{2^v}\left( {1 - t_k^2} \right){{\left[ {J_n^{\left( {\lambda ,0} \right)'}{t_k}} \right]}^2}} \over f}\left( {{t_k}} \right)} \end{array}$$(9)

where A_k are weights of the quadrature.

Transforming the interval [0, t] into 〈−1,1〉

(t−t02)ν∫−11f(u)(1−u)λdu$$\begin{array}{} \displaystyle \left( \frac{t-t_{0}}{2} \right)^{\nu} \int_{-1}^{1}\frac{f\left ( u \right )}{\left ( 1-u \right )^{\lambda}}du \end{array}$$

where

f(u)=f((t−t02)u+(t+t02)),$$\begin{array}{} \displaystyle f \left ( u \right )=f\left ( \left ( \frac{t-t_{0}}{2} \right) u+\left ( \frac{t+t_0}{2} \right ) \right ), \end{array}$$

we obtain a formula (10) for Gauss-Jacobi quadrature (8) that is perfectly suited for FOD/FOI computation by applying formulas (1)-(6), i.e. kernel in integrand is not integrated but computed by applying the formula for the weights A_k of the quadrature and it does not determine accuracy of integration anymore.

From now on Riemann-Liouville fractional derivative formula (6) assumes the following computational form

t0RLDt(ν)f(t)=∑k=0n−1(t−t0)k−νf(k)(t0)Γ(k−ν+1)+1Γ(n−ν)(t−t02)n−ν∫−11f(n)(u)(1−u)ν−n+1du=∑k=0n−1(t−t0)k−νf(k)(t0)Γ(k−ν+1)+1Γ(n−ν)(t−t02)n−ν∑k=1n2ν(1−uk2)[Jn(ν−n+1,0)'uk]︸kernelf(n)(uk)$$\begin{array}{} ^{RL}_{t_{0}}\textrm{D}_{t}^{\left ( \nu \right )}f\left(t\right)=\sum_{k=0}^{n-1}\frac{ \left(t-t_{0}\right)^{k-\nu} f^{\left( k \right)}\left ( t_{0} \right )}{\Gamma\left ( k-\nu+1 \right )} +\frac{1}{\Gamma\left ( n-\nu \right )} \left ( \frac{t-t_{0}}{2} \right )^{n-\nu}\int_{-1}^{1}\frac{ f^{\left(n\right)} \left(u\right)}{\left ( 1-u \right )^{\nu-n+1}}du \\ = \sum_{k=0}^{n-1}\frac{ \left(t-t_{0}\right)^{k-\nu} f^{\left( k \right)}\left ( t_{0} \right )}{\Gamma\left ( k-\nu+1 \right )} +\frac{1}{\Gamma\left ( n-\nu \right )} \left ( \frac{t-t_{0}}{2} \right )^{n-\nu} \underbrace{\sum_{k=1}^{n} \frac{2^{\nu}}{ \left( 1-u^{2}_{k} \right) \left[ J_{n}^{\left(\nu-n+1,0\right){'}} u_{k} \right]^{2}}}_\text{ kernel}\;f^{\left(n\right)}\;\left(u_{k}\right). \end{array}$$(10)

in which n = ⌈ν⌉.

The ordinarily sufficient method of Gauss quadrature 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 [10]. However, its application becomes difficult if nodes and weights must be customized and computed, e.g. if the integrand does not meet fully or at all the requirements of a quadrature or required integration interval is different than foreseen for a quadrature. In all of the cases nodes and weights must be computed prior the integration.

Customized variants of the Gauss-Jacobi quadrature (9) for computation of FOD/FOI by applying formulas (1)-(6) require computation of polynomial of order n, its derivative, nodes and weights.

Existing approaches for computing the nodes and the weights, some of which have been widely used for many years, suffer from O(n2)$\begin{array}{} \displaystyle {\cal O}\left( {{n^2}} \right) \end{array}$ complexity or error which grows with n. With some optimizations, e.g. application of technique which utilizes asymptotic formulas for accurate initial guesses of the roots and efficient evaluations of the degree n Jacobi polynomial, it can be reduced to O(n)$\begin{array}{} \displaystyle {\cal O}\left( n \right) \end{array}$ as it is presented in [11].

Fractional order derivative calculation by applying formulas (1)-(6) requires computation of n-th derivative of a function, where n is an integer. The following formulas are used for that purpose in the paper [12, 13]:

First derivative 3-points stencil Central Difference

f'(x0)=f1−f−12h+O(h2),$$\begin{array}{} \displaystyle f'({x_0}) = {{{f_1} - {f_{ - 1}}} \over {2h}} + {\cal O}({h^2}), \end{array}$$(11)

where f₁ ≡ f(x₀ + h) and f₋₁ ≡ f(x₀ − h).

Second derivative 3-points stencil Central Difference

f''(x0)=f−1−2f0+f1h2+O(h2)$$\begin{array}{} \displaystyle f''({x_0}) = {{{f_{ - 1}} - 2{f_0} + {f_1}} \over {{h^2}}} + {\cal O}({h^2}) \end{array}$$(12)

Third derivative 3-points stencil Central Difference

f'''(x0)=−f−2−2f−1−2f1+f22h3+O(h2)$$\begin{array}{} \displaystyle f'''({x_0}) = {{ - {f_{ - 2}} - 2{f_{ - 1}} - 2{f_1} + {f_2}} \over {2{h^3}}} + {\cal O}({h^2}) \end{array}$$(13)

The limiting factor to the accuracy of computation by applying the formulas (1)-(6) is a precision which uses the computer to store data that are being supplied to it.

The selection of uniform C++ equipped with the standard mathematical library as a main programming tool is not enough nowadays to take full advantage of available hardware. The application of arbitrary precision computing for increasing accuracy and correctness of numerical calculations and Nvidia CUDA parallelization technology for their effectiveness, are the best examples in this context [14].

Application of arbitrary precision makes it possible for the user to choose precision for calculation and for each variable storing a value. It is not machine-depended or IEEE standard types. With its help we can - among the others - increase general accuracy of mathematical computations. However, its application purpose is above all to increase accuracy of numerical calculations, e.g. by eliminating under- and overflows, increasing accuracy of a polynomial zeros finding and derivative and integral calculating.

The importance of elimination of limited precision in computer calculations was aptly presented by Toshio Fukushima in The Astronomical Journal in 2001 by giving the following example: “In the days of powerful computers, the errors of numerical integration are the main limitation in the research of complex dynamical systems, such as the long-term stability of our solar system and of some exoplanets [...]” and gives an example where using double precision leads to an accumulated round-off error of more than 1 radian for angular position of planets [15].

Double precision computer arithmetic is optimized for speed and has many flaws which negatively influence the accuracy of computations, e.g. limitations of number values which double precision variables can hold or no programmer influence on mathematical operations rounding.

However, it is the lack of clarity in handling of intermediate results which troubles the most, i.e. the floating-point standard [16] only defines that the results must be rounded correctly to the destination’s precision and fails to define the precision of destination variable. This choice is commonly made by a system or a programming language. The user can not influence it in any way. Therefore, the same program can return significantly different results depending on the implementation of the IEEE standard.

Arbitrary precision is applied in conjunction with special libraries which include their own data structures and mathematical functions.

Accuracy and the efficiency assessment are crucial components in the process of programming because it enables examining a program for its proper operation and accuracy.

Common approach to the accuracy assessment requires an exact value.

In case of FOD/FOI computation, an effective accuracy assessment is difficult, sometimes not possible due to general lack of formulas for exact values.

Despite the availability of a handful of analytical formulas for fractional order ν=12$\begin{array}{} \displaystyle \nu=\frac{1}{2} \end{array}$ and some computational formulas, they are accessible for selected types of functions only. Some other formulas are in form of series expansion only.

The following formulas for FOD/FOI are used further in the paper for calculations of their exact values:

Trigonometric functions:

0Dt(ν)sin(t)=t1−ν∑k=0∞(−1)kt2kΓ(2k+2−ν),for0<ν<1,$$\begin{array}{} \displaystyle _{0}\textrm{D}_{t}^{\left ( \nu \right )}\sin \left(t\right)=t^{1-\nu} \sum_{k=0}^{\infty}\frac{\left(-1\right)^{k} t^{2 k}}{\Gamma\left(2 k+2-\nu\right)},\; \text{for}\; 0<\nu<1, \end{array}$$(14)

0Dt(ν)cos(t)=t2−ν∑i=0∞(−1)i+1t2i+1Γ(2i+4−ν),for1<ν<2.$$\begin{array}{} \displaystyle _{0}\textrm{D}_{t}^{\left ( \nu \right )}\cos \left(t\right)=t^{2-\nu} \sum_{i=0}^{\infty}\frac{\left(-1\right)^{i+1} t^{2i+1}}{\Gamma\left(2i+4-\nu\right)},\; \text{for}\; 1<\nu<2. \end{array}$$(15)

Power functions

0Dt(ν)tβ=Γ(β+1)Γ(β−ν+1)tβ−ν.$$\begin{array}{} \displaystyle _{0}\textrm{D}_{t}^{\left ( \nu \right )}t^{\beta}=\frac{\Gamma\left(\beta+1\right)}{\Gamma\left(\beta-\nu+1\right)}t^{\beta-\nu}. \end{array}$$(16)

Exponential function f (t) = e^(at) by applying the Mittag-Leffler function with α = β = 1

0Dt(ν)e(at)≈0Dt(ν)Eα,β(at)=t−ν∑k=0∞(at)kΓ(k+1)Γ(k+1−ν)Γ(νk+1).$$\begin{array}{} \displaystyle _{0}D_{t}^{\left(\nu\right)}e^{\left ( at \right )} \approx\;_{0}D_{t}^{\left(\nu\right)}\, E_{\alpha,\beta}\left ( at \right ) =t^{-\nu} \sum_{k=0}^{\infty}\frac{\left ( at \right )^{k}\Gamma\left(k+1\right)}{\Gamma\left ( k+1-\nu \right )\Gamma\left(\nu k+1\right)}. \end{array}$$(17)

In practice ∞ is replaced by large number N.

In our work we have to be able to assess the accuracy of FOD/FOI computation of any function. Thus we have developed a method of accuracy assessment by applying operators concatenation [17] of fractional order differentiation and integration.

The main idea behind this versatile criteria presented in Figure 3 is the ability to compare computed values of FOD/FOI with values of classical - integer order derivatives and integrals, which can be calculated by applying well known (and accurate) formulas.

Fig. 3

To compare computational accuracy of FOD we apply right- and left-hand definitions schemas presented in Figures 1 and 2 to obtain FOD of order ν = 0.3 and 2.3 for the following functions and ranges:

f(t)=e−t,t∈(0,5),$$\begin{array}{} \displaystyle f\left(t\right)=e^{-t},\;t \in \left(0,5\right), \end{array}$$(18)

f(t)=et,t∈(0,3),$$\begin{array}{} \displaystyle f\left(t\right)=e^{t},\;t \in \left(0,3\right), \end{array}$$(19)

f(t)=t,t∈(0,1),$$\begin{array}{} \displaystyle f\left(t\right)=\sqrt{t},\;t \in \left(0,1\right), \end{array}$$(20)

f(t)=sin(t),t∈(0,2π).$$\begin{array}{} \displaystyle f(t) = \sin (t),\:t \in (0,2\pi ). \end{array}$$(21)

To compare computational accuracy of FOD we use relative error and infinity norm ||e_r||_L∞ to measure similarity of two vectors x and x^$\begin{array}{} \hat x \end{array}$. Vector x contains 100 computed values in a given range and vector x^$\begin{array}{} \displaystyle \hat x \end{array}$ contains reference values assumed, i.e. values computed by applying formulas and methods listed in subsection 5):

‖er‖L∞=maxi‖xi−x^i‖‖xi‖.$$\begin{array}{} \displaystyle {||e_r||}{_{{L_\infty }}} = ma{x_i}\frac{{{||x_i} - {{\hat x}_i||}}}{{{||x_i||}}}. \end{array}$$(22)

Tables 1-4 contain maximum relative error between vectors x and x^$\begin{array}{} \displaystyle \hat x \end{array}$ for 100 values of FOD computed in a given range t₀ − t, for N = 4,8,...,64 order of polynomial applied for the method of integration and d -differentiation step required for integer order derivative computation.

Table 1

RL||er||L∞,h=t−t0100points$\begin{array}{} \displaystyle {\rm{RL||}}{e_r}{\rm{|}}{{\rm{|}}_{L\infty }},h = {{t - {t_0}} \over {100points}} \end{array}$ ,d - differentiation step, n − ν = 0.7, n = 1,ν = 0.3

Table 2

C||er||L∞,h=t−t0100points$\begin{array}{} \displaystyle C||{e_r}|{|_{L\infty }},h = {{t - {t_0}} \over {100{\rm{ }}points}} \end{array}$d - differentiation step, n = 1; n−v = 0:7; v = 0:3

Table 3

RL||er||L∞,h=t−t0100points$\begin{array}{} \displaystyle {\rm{RL||}}{e_r}{\rm{|}}{{\rm{|}}_{L\infty }},h = {{t - {t_0}} \over {100{\rm{ }}points}} \end{array}$d - differentiation step, n− v = 0:7; n = 3; v = 2:3

Table 4

C||er||L∞,h=t−t0100points$\begin{array}{} \displaystyle C||{e_r}||{L_\infty },h = \frac{{t - {t_0}}}{{100points}} \end{array}$d - differentiation step, n = 3; n− v = 0:7; v = 2:3

Results presented in tables 1-4 confirm shape of integrand as the main factor which determines accuracy of FOD computation. The term ”shape of integrand” describes the behavior of a function in a given interval. In this scope: if we say ”dramatical shape changes of integrand”, we mean certain behavior of a function, in which the values of its derivatives (1^st, 2^nd and higher) assume high values.

The integrand included in RL and C formulas can be divided into a kernel and a function, i.e.

0RLDt(ν)f(t)=1Γ(ν)∫0t(t−τ)ν−n︸kernelf(τ)dτ,forn−1<ν<n.$$\begin{array}{} \displaystyle _{0}^{RL}\textrm{D}_{t}^{\left ( \nu \right )}f\left ( t \right )=\frac{1}{\Gamma\left ( \nu \right )}\int_{0}^{t} \underbrace{\left ( t-\tau \right )^{\nu-n}}_\text{kernel}\; f\left ( \tau \right ) d\tau,\; \text{for} \;n-1<\nu<n. \end{array}$$

Fig. 4

The kernel, which includes singularity determines difficulty of FOD/FOI computation. By applying Gauss-Jacobi quadrature for computation, it is possible to exclude the kernel from numerical integration and compute it by applying formulas for the weights of the quadrature (9) which assures significantly higher accuracy. However, results presented in tables 1-4 show that the accuracy can also be decreased by the shape of the function. For this reason an amount of sampling points (N) required for highest accuracy of FOD for a particular function varies from 4 to 64.

Accuracy of integer order numerical differentiation depends mainly on the width of sampling step of a function. As it is easy to see, the width required for the highest accuracy is beyond the capabilities of the standard double precision computer arithmetic. In this scope, an application of high precision libraries for programming is required. How many points are included is of secondary importance. However, 3-points stencil central difference is an optimum for the most cases. More points does not increase the accuracy.

FOD can be computed by using numerous approaches. It includes application of integer order differentiation of fractional order integral (Riemann-Liouville formula) or fractional order integration of integer order derivative (Caputo formula).

The goal of numerical experiment described in the paper was to compare accuracy of FOD computation by applying both approaches.

During the experiment we calculated FOD of selected orders for exponential, power and periodic functions in various ranges.

As it was mentioned in the problem description, computation of FOD with high accuracy is a real challenge for the computer and it can not be succeed by applying standard numerical methods. For this purpose it is required application of Gauss-Jacobi quadrature with adapted weight function for integration part and central differences method of higher order for integer order differentiation. Additionally, the standard double precision C++ programming language limitations were eliminated by applying for programming arbitrary precision with a help of mathematical libraries (GNU MPFR/GNU GMP and MPFR C++).

The resulting conclusion is that high accuracy computation of FOD can be achieved by applying both approaches (RL and C), provided the appropriate numerical method of integration and differentiation is used. However, due to the fact that accuracy of numerical integration and differentiation generally depends on a shape of a function, by applying an appropriate sequence of numerical calculations (RL or C formula) selected according to a shape of a particular function, can result with significantly higher accuracy.

An important advantage of RL formula over C formula is its capability of computing FOD of function which n^th derivative does not exist (see results for computation of FOD of order 2.3 for a power function). In such case, application of C formula gives nan (not a number) result, where in the same situation, RL method results with accurate value of FOD.

Remark: To be able to apply C formula for computation of FOD for any function, it is required to apply equivalent Riemann-Liouville formula (6).