Accuracy Problems of Numerical Calculation of Fractional Order Derivatives and Integrals Applying the Riemann-Liouville/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

1−τ→e1−1/u$$\begin{array}{} \displaystyle 1 - \tau \to {e^{1 - 1/u}} \end{array}$$(8)

in kernel of formula (5) and (6) removes singularity from the integrand and it makes it smoother.

For an example function f(t) = 1(t), t₀ = 0,t = 1 and fractional order of integral ν = 0.5 we calculate

∫01(1−τ)−0.5︸g(τ)dx=|1−τ=e1−1/u−dτ=e1−1/u⋅1u2dτ=(e1−1/u)⋅1u2|=…=…∫01(e1−1/u)0.5⋅1u2︸G(u)du.$$\begin{array}{} \displaystyle \int_{0}^{1} \underbrace{\left ( 1-\tau \right )^{-0.5}}_\text{$g\left(\tau\right)$}dx & = \begin{vmatrix} 1-\tau = e^{1-1/u} \\ -d\tau = e^{1-1/u}\cdot\frac{1}{u^{2}} \\ d\tau = \left ( e^{1-1/u} \right ) \cdot\frac{1}{u^{2}} \end{vmatrix} = \dots \nonumber \\ & = \dots \int_{0}^{1} \underbrace{\left ( e^{1-1/u} \right )^{0.5}\cdot \frac{1}{u^{2}}}_\text{$G\left(u\right)$} du. \end{array}$$(9)

τ=1−e1−1/u,u=11−ln(1−τ).$$\begin{array}{} \displaystyle \begin{array}{*{20}{l}} {\tau = 1 - {e^{1 - 1/u}},}\\ {u = \frac{1}{{1 - ln\left( {1 - \tau } \right)}}.} \end{array} \end{array}$$(10)

Transformation (9) and integration interval change (10) is presented in Fig. 2.

Fig. 2

Transformation (9) is now applied for practical accuracy increase using two periodic functions:

f(t)=sin(t),t∈(0,2π),$$\begin{array}{} \displaystyle f\left( t \right) = \sin \left( t \right),\;t \in \left( {0,2\pi } \right), \end{array}$$(11)

f(t)=1.5cos(2t)+2.2cos(4t),t∈(0,2π).$$\begin{array}{} \displaystyle f\left( t \right) = 1.5\cos \left( {2t} \right) + 2.2\cos \left( {4t} \right),\;t \in \left( {0,2\pi } \right). \end{array}$$(12)

Fractional order integration formula (6) for the function (11) and order (−ν) assumes the following form

0RLD2π(−ν)sin(t)=1Γ(ν)∫02π(e2π−2π/u)ν2πu2sin(2π−e2π−2π/u)du$$\begin{array}{} \displaystyle _0^{RL}{\rm{D}}_{2\pi }^{\left( { - \nu } \right)}\sin \left( t \right) = \frac{1}{{\Gamma \left( \nu \right)}}\int_0^{2\pi } {{{\left( {{e^{2\pi - 2\pi /u}}} \right)}^\nu }} \frac{{2\pi }}{{{u^2}}}\sin \left( {2\pi - {e^{2\pi - 2\pi /u}}} \right)du \end{array}$$

and for the function (12) assumes the form

0RLD2π(−ν)1.5cos(2t)+2.2cos(4t)=1Γ(ν)∫02π(e2π−2π/u)ν2πu2+1.5cos(2(2π−e2π−2π/u))+2.2cos(4(2π−e2π−2π/u))du.$$\begin{array}{} \displaystyle _0^{RL}{\rm{D}}_{2\pi }^{\left( { - \nu } \right)}1.5\cos \left( {2t} \right) + 2.2\cos \left( {4t} \right) = \frac{1}{{\Gamma \left( \nu \right)}}\int _0^{2\pi }{\left( {{e^{2\pi - 2\pi /u}}} \right)^\nu }\frac{{2\pi }}{{{u^2}}}\, + 1.5\cos \left( {2\left( {2\pi - {e^{2\pi - 2\pi /u}}} \right)} \right)\, + 2.2\cos \left( {4\left( {2\pi - {e^{2\pi - 2\pi /u}}} \right)} \right)du. \end{array}$$

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 ν = 0.2,0.5,0.8, there are presented only plots for fractional order integrals of order ν = −0.5 and derivatives of order ν = 0.5.

Fig. 3

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 [28], [3] and [12]. 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 [29].

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.

For fractional derivatives of order (ν), differentiation formula (5) for the function (11) assumes the following form

0RLD2π(ν)sin(t)=(2π−ν)sin(0)Γ(1−ν)+1Γ(1−ν)∫02π(e2π−2π/u)ν2πu2cos(2π−e2π−2π/u)du$$\begin{array}{} \displaystyle _0^{RL}{\rm{D}}_{2\pi }^{\left( \nu \right)}\sin \left( t \right) = \frac{{\left( {2{\pi ^{ - \nu }}} \right)\sin \left( 0 \right)}}{{\Gamma \left( {1 - \nu } \right)}} + \frac{1}{{\Gamma \left( {1 - \nu } \right)}}\int_0^{2\pi } {{{\left( {{e^{2\pi - 2\pi /u}}} \right)}^\nu }} \frac{{2\pi }}{{{u^2}}}\cos \left( {2\pi - {e^{2\pi - 2\pi /u}}} \right)du \end{array}$$

and for the function (12) assumes the form

0RLD2π(ν)1.5cos(2t)+2.2cos(4t)=(2π−ν)1.5cos(0)+2.2cos(0)Γ(1−ν)+1Γ(1−ν)∫02π(e2π−2π/u)ν2πu2(−3sin(2(2π−e2π−2π/u)))−8.8sin(4(2π−e2π−2π/u))du.$$\begin{array}{} \displaystyle _0^{RL}{\rm{D}}_{2\pi }^{\left( \nu \right)}1.5\cos \left( {2t} \right) + 2.2\cos \left( {4t} \right) = \frac{{\left( {2{\pi ^{ - \nu }}} \right)1.5\cos \left( 0 \right) + 2.2\cos \left( 0 \right)}}{{\Gamma \left( {1 - \nu } \right)}} + \frac{1}{{\Gamma \left( {1 - \nu } \right)}}\int _0^{2\pi }{\left( {{e^{2\pi - 2\pi /u}}} \right)^\nu }\frac{{2\pi }}{{{u^2}}}\left( { - 3\sin \left( {2\left( {2\pi - {e^{2\pi - 2\pi /u}}} \right)} \right)} \right) - 8.8\sin \left( {4\left( {2\pi - {e^{2\pi - 2\pi /u}}} \right)} \right)du. \end{array}$$

The accuracy of computations are presented in Figs 4(a) and 4(b).

Fig. 4

Final remark:

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 [30]. It also became a subject of works of Schwartz [31], Stenger [32], Takahasi [33], Mori [34], Haber [35] 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.

Two solutions of this problem, which involve another substitution at the singularity point have been published in [36], [37], [34] and lately in [38].

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.

In the following section, the DE Transformation is applied to calculate fractional orders derivatives and integrands applying formulas (5) and (6).

The DE Transformation implemented for the purpose of fractional orders derivatives and integrands computations is based on the transformation, which was proposed by Schwartz [31] as the Tanh Rule. It can be described as follows:

Let us consider the integral

S=∫abf(x)dx,$$\begin{array}{} \displaystyle S = \int_a^b f \left( x \right)dx, \end{array}$$(13)

where f(x) is integrable on interval (a,b). The function f(x) may have singularity at x = a or x = b or at both ends.

First we apply the following variable transformation

x=ϕ(t),ϕ(−∞)=a,ϕ(∞)=b.$$\begin{array}{} \displaystyle x = \phi \left( t \right),\;\phi \left( { - \infty } \right) = a,\;\phi \left( \infty \right) = b. \end{array}$$(14)

We obtain

S=∫−∞∞f(ϕ(t))ϕ′(t)dt.$$\begin{array}{} \displaystyle S = \int _{ - \infty }^\infty f\left( {\phi \left( t \right)} \right)\phi \prime \left( t \right)dt. \end{array}$$(15)

It is important that ϕ(t) possesses the property such as ϕ′(t) converges at least at single exponential rate as t → ±∞

|ϕ′(t)|→exp(−cexp(|t|)),$$\begin{array}{} \displaystyle \left| {\phi \prime \left( t \right)} \right| \to \exp \left( { - c\exp \left( {\left| t \right|} \right)} \right), \end{array}$$(16)

where c is some constant. After that, it is best to apply the Trapezoidal Rule to the transformed integrand expression [32], [39]

S=h∑n=−∞∞f(ϕ(nh))ϕ′(nh),$$\begin{array}{} \displaystyle S = h\mathop \sum \limits_{n = - \infty }^\infty f\left( {\phi \left( {nh} \right)} \right)\phi' \left( {nh} \right), \end{array}$$

where h is subinterval width. Due to the property (16) truncation of the summation process can be done at some arbitrarily chosen n = −N and n = +N

S=h∑n=−N+Nf(ϕ(nh))ϕ′(nh),$$\begin{array}{} \displaystyle S = h\mathop \sum \limits_{ + N}^{n = - N} f\left( {\phi \left( {nh} \right)} \right)\phi \prime \left( {nh} \right), \end{array}$$(17)

N=−N+N+1,$$\begin{array}{} \displaystyle N = - N + N + 1, \end{array}$$

where N states the amount of subintervals.

Since ϕ′(nh) as well as the whole expression f(ϕ(nh))ϕ′ (nh) converges at least at exponential rate at large |n|, the quadrature formula (17) is called the Double Exponential [34].

Due to truncation of the summation process in (17) at some arbitrary chosen n = −N, n = +N, function f(x) can have singularities at x = a and/or x = b as long as it is integrable over the integration interval.

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 N. The optimal strategy is to get both errors equal [38], [34].

The subinterval width h, that defines the evaluation step and the number of sampling points are key values in such strategy. The source [38] suggests the following value of h for the DE Formula

h~log(2πNω/c)N,$$\begin{array}{} \displaystyle h \sim \frac{{\log (2\pi N\omega /c)}}{N}, \end{array}$$

where c is some constant to be taken, usually 1 or π/2 and ω is the distance to the nearest singularity of the integrand defined in [38], [34].

The following equations (18)-(20) present three different transformations with varied convergence rate [33], [38]. Transformation (18) has the slowest convergence and transformation (20), the fastest respectively:

x=ϕ(t)=tanhtp,ϕ′(t)=ptp−1cosh2tp,p=1,3,5,…,$$\begin{array}{} \displaystyle x = \phi \left( t \right) = \tanh \,{t^p},\,\phi \prime \left( t \right) = \frac{{p{t^{p - 1}}}}{{{{\cosh }^2}{t^p}}},\,p = 1,3,5, \ldots , \end{array}$$(18)

x=ϕ(t)=tanh(π2sinht),ϕ′(t)=π2coshtcosh2(π2sinht),$$\begin{array}{} \displaystyle x = \phi \left( t \right) = \tanh \left( {\frac{\pi }{2}\sinh \,t} \right),\,\phi \prime \left( t \right) = \frac{{\frac{\pi }{2}\,\cosh \,t}}{{{{\cosh }^2}\left( {\frac{\pi }{2}\,\sinh \,t} \right)}}, \end{array}$$(19)

x=ϕ(t)=tanh(π2sinht3),ϕ′(t)=32πt2cosht3cosh2(π2sinht3).$$\begin{array}{} \displaystyle x = \phi \left( t \right) = \tanh \left( {\frac{\pi }{2}\sinh \,{t^3}} \right),\,\phi \prime \left( t \right) = \frac{{\frac{3}{2}\pi \,{t^2}\,\cosh \,{t^3}}}{{{{\cosh }^2}\left( {\frac{\pi }{2}\,\sinh \,{t^3}} \right)}}. \end{array}$$(20)

Applying the transformation (19) to a function according to the formula (17), we obtain the following trapezoidal form

S=h∑i=1Nf(b−a2xi+b+a2)wi,$$\begin{array}{} \displaystyle S = h\mathop \sum \limits^N_{i = 1} f\left( {\frac{{b - a}}{2}\,{x_i} + \frac{{b + a}}{2}} \right){w_i}, \end{array}$$(21)

where

xi=f(tanh(π2sinhti)),$$\begin{array}{} \displaystyle {x_i} = f\left( {\tanh \left( {\frac{\pi }{2}\;\sinh \;{t_i}} \right)} \right), \end{array}$$

ti=−ta+(i−1)h,i=1,2,3,⋯N,h=2taN,$$\begin{array}{} \displaystyle {t_i} = - {t_a} + \left( {i - 1} \right)h,\;i = 1,2,3, \cdots N,h = \frac{{2{t_a}}}{N}, \end{array}$$

wi=π2coshticosh2(π2sinhti)⋅(b−a2)$$\begin{array}{} \displaystyle {w_i} = \frac{{\frac{\pi }{2}\;\cosh \;{t_i}}}{{{{\cosh }^2}\left( {\frac{\pi }{2}\sinh \;{t_i}} \right)}} \cdot \left( {\frac{{b - a}}{2}} \right) \end{array}$$

are the nodes, the sampling points and the weights of the DE Quadrature respectively.

Parameter t_a can be calculated using the following formula [34]

ta=precision’sdigits0.46,$$\begin{array}{} \displaystyle {t_a} = {\rm{precision}}{\rm{s digit}}{{\rm{s}}^{0.46}}, \end{array}$$(22)

According to [34] asymptotic error of the formula (21) in terms of subinterval width h of the Trapezoidal Rule is expressed as

|ΔSh|≈exp(−ch),$$\begin{array}{} \displaystyle \left| {\Delta {S_h}} \right| \approx \exp \left( { - \frac{c}{h}} \right), \end{array}$$

and asymptotic error in terms of the number N of the functions evaluations is

|ΔSh(N)|≈exp(−cNlogN).$$\begin{array}{} \displaystyle \left| {\Delta S_h^{\left( N \right)}} \right| \approx \exp \left( { - c\frac{N}{{\log N}}} \right). \end{array}$$

4.2.1

Programming Issues

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 sec² (csinh t) can overflow. Second, the selection of the parameter t_a > 3 (22) can cause underflow and the occurrence of uncontrolled errors (the computed integral values replaced by random numbers).

Increasing value of the parameter t_a we can integrate further left and right from y axis (see Fig. 5(c)) and it increases the accuracy. However, values of integrated function can assume very small values. Too small for double precision variables to distinguish them from 0. The underflow occurs. This negative process can be eliminated by applying high precision variables. Application of GNU MPFR library makes it possible. The library have also at their disposal many mathematical functions which can be used with the high precision variables.

Fig. 5

For t_a = 4.97 there are required 100-digits variables, for t_a = 5.5 - 200-digits variables, for t_a = 6.5 - 500-digits variables, for t_a = 7.5 - 1500-digits variables and for t_a = 8.7 - 5000-digits variables.

Another parameter which influence accuracy of the DE Quadrature (21) is an amount N of subintervals for integration applying the Trapezoidal Rule.

4.2.2

Accuracy of the DE Transformation for Fractional Derivatives and Integrals Computations

For accuracy assessment of the DE Transformation there are selected three functions:

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

f(t)=sin(t),t∈(0,2π),$$\begin{array}{} \displaystyle f\left( t \right) = \sin \left( t \right),\;t \in \left( {0,2\pi } \right), \end{array}$$(24)

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

and four fractional orders:

Fractional derivative of order: ν = 0.9 and ν = 0.5

Fractional integral of order: ν = −0.1 and ν = −0.0001,

Results are presented in Figs. 6, 7 and 8.

Fig. 6

Fig. 7

Fig. 8

Final remarks:

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 y axis it is possible to integrate. The latter is difficult due to very small values occurrence on both ends of the transformed integration interval which cause underflows applying standard precision in computations.

The accuracy increase up to ∼ 10⁻⁵⁰ can be obtained by careful programming applying variables with increased precision. Still, some orders are out of reach as for example the order ν = −0.0001, calculation of which ends up with 90% relative error.

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 [28], [40] 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

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 \gt - 1. \end{array}$$(26)

General, well known formula of 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\limits_{k = 1}^n {{A_k}} f\left( {{x_k}} \right) + R \end{array}$$

with the weight function (26) assumes the form:

∫−11(1−x)λ(1+x)β⋅f(x)dx≈∑k=1nAkf(xk),$$\begin{array}{} \displaystyle \int_{ - 1}^1 {{{\left( {1 - x} \right)}^\lambda }} {(1 + x)^\beta } \cdot f\left( x \right)dx \approx \sum\limits_{k = 1}^n {{A_k}} f\left( {{x_k}} \right), \end{array}$$(27)

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

The Jacobi Polynomial can be determined by applying Rodrigues formula

Jn(λ,β)(x)=(−1)n2n⋅n!(1−x)−λ(1+x)−βdndxn[(1−x)λ+n(1+x)β+n]$$\begin{array}{} \displaystyle J_n^{\left( {\lambda ,\beta } \right)}\left( x \right) = \frac{{{{\left( { - 1} \right)}^n}}}{{2n \cdot n!}}{\left( {1 - x} \right)^{ - \lambda }}{\left( {1 + x} \right)^{ - \beta }}\frac{{{d^n}}}{{d{x^n}}}\left[ {{{\left( {1 - x} \right)}^{\lambda + n}}{{\left( {1 + x} \right)}^{\beta + n}}} \right] \end{array}$$(28)

and weights

Ak=2λ+β+1Γ(λ+n+1)Γ(β+n+1)n!Γ(λ+β+n+1)1(1−xk2)[Jn(λ,β)′xk]2,$$\begin{array}{} \displaystyle {A_k} = {2^{\lambda + \beta + 1}}\frac{{\Gamma \left( {\lambda + n + 1} \right)\Gamma \left( {\beta + n + 1} \right)}}{{n!\Gamma \left( {\lambda + \beta + n + 1} \right)}}\frac{1}{{\left( {1 - x_k^2} \right){{\left[ {J_n^{\left( {\lambda ,\beta } \right)\prime }{x_k}} \right]}^2}}}, \end{array}$$(29)

in which Jn′(xk)$\begin{array}{} \displaystyle J_n^\prime \left( {{x_k}} \right) \end{array}$ is derivative Jn(λ,β)′(xk)$\begin{array}{} \displaystyle J_n^{\left( {\lambda ,\beta } \right)\prime }\left( {{x_k}} \right) \end{array}$ of Jacobi polynomial of degree n and

R=2λ+β+2n+1λ+β+2n+1Γ(λ+n+1)Γ(β+n+1)Γ(λ+β+n+1)Γ2(λ+β+2n+1)n!2nf2n(ξ),ξ∈〈−1,1〉=O(n2)$$\begin{array}{} \displaystyle R = \frac{{{2^{\lambda + \beta + 2n + 1}}}}{{\lambda + \beta + 2n + 1}}\frac{{\Gamma \left( {\lambda + n + 1} \right)\Gamma \left( {\beta + n + 1} \right)\Gamma \left( {\lambda + \beta + n + 1} \right)}}{{{\Gamma ^2}\left( {\lambda + \beta + 2n + 1} \right)}}\frac{{n!}}{{2n}}{f^{2n}}\left( \xi \right),\,\xi \in \langle - 1,1\rangle = {\scr O}({n^2}) \end{array}$$

is the remainder of integral approximation.

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.

Fig. 9

Fig. 10

The Gauss Jacobi Quadrature is applied for calculations of fractional order integrals ν = −0.0001, ν = −0.1 and derivatives ν = 0.5, ν = 0.9 of exponential function (23), trigonometric function (24) and the function (25).

Next, GJ is evaluated with N = 32 degree of polynomial applied in integral approximation. The value is an outcome from the results of the next part accuracy assessment of the method, i.e. comparison of accuracy performance of GJ for various degree of polynomial. It is presented in Figs. 11 and 12.

Fig. 11

Fig. 12

Now, GJ is evaluated for accuracy of calculations for different fractional orders 0 < ν < 1 and different integration intervals (1(1)6). Results are presented in Figs. 13.

Fig. 13

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.

Fig. 14

Final remarks:

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 N = 32 and 100-digit precision does not extend significantly time complexity of calculations (see Fig. 14).