INFORMATICA Informatica 1822-88440868-49520868-4952 Vilnius University INFOR540 10.15388/23-INFOR540 Research Article Numerical Approximations of the Riemann–Liouville and Riesz Fractional Integrals https://orcid.org/0000-0002-3695-9569 CiesielskiMariuszmariusz.ciesielski@icis.pcz.pl∗

M. Ciesielski received his MSc degree in computer science and PhD degree in mechanical engineering from the Czestochowa University of Technology, in 2000 and 2005, respectively. He has been working in the Faculty of Mechanical Engineering and Computer Science there since 2000. Currently, he is an assistant professor in the Department of Computer Science. He is an author and a co-author of over 100 scientific papers, 42 of them were published in Web of Science journals. His research interests include the numerical modelling of heat transfer processes, computational methods of mechanics, numerical algorithms, fractional calculus and its applications.

 https://orcid.org/0000-0003-1206-6923 GrodzkiGrzegorzgrzegorz.grodzki@icis.pcz.pl

G. Grodzki received his MSc and PhD degree in mechanic science from the Czestochowa University of Technology, Poland, in 1992 and 2000, respectively. At this time, he was conducting research measurements with use LDA (Laser Doppler Anemometry), and digital analysis and processing of randomly sampled signals. Since 2000, he worked at the Institute of Mathematics and Computer Science. Currently, he is an assistant professor at the Department of Computer Science, Czestochowa University of Technology, Poland. He is an author of over 30 papers in refereed journal and conference papers. He is also the co-author of 2 scientific monographs. His current interests include security informatics systems and computer network, optimization of organization strategy.

 Department of Computer Science, Czestochowa University of Technology, Dabrowskiego 73, 42-200 Czestochowa, Poland Corresponding author. 202422112023351214652023112023 © 2024 Vilnius University2024 Open access article under the CC BY license.

In this paper, the numerical algorithms for calculating the values of the left- and right-sided Riemann–Liouville fractional integrals and the Riesz fractional integral using spline interpolation techniques are derived. The linear, quadratic and three variants of cubic splines are taken into account. The estimation of errors using analytical methods are derived. We show four examples of numerical evaluation of the mentioned fractional integrals and determine the experimental rate of convergence for each derived algorithm. The high-precision calculations are executed using the 128-bit floating-point numbers and arithmetic routines.

 fractional calculus numerical integration numerical algorithms spline interpolation Introduction

The fractional calculus (e.g. Kilbas et al., 2006; Podlubny, 1999; Rubin, 1996) is treated as a branch of mathematical analysis dealing with differential-integral equations, in which the integrals are of the convolutional type and have power-type kernels. The fractional calculus refers to the fractional integration and fractional differentiation. It can be seen in the extensive literature in the past few decades that the fractional calculus has received increasing attention for its applications in science and engineering – it is difficult to list examples of literature without omitting even the most important ones. This work mainly deals with the development and investigation of numerical methods for the fractional integrals.

We recall the following definitions of the Riemann–Liouville fractional integrals. The left-sided fractional integral of order α > 0 of the given function y ( x ) on the interval [ a , b ] is defined as follows (Kilbas et al., 2006; de Oliveira and Machado, 2014) I a + α y ( x ) = 1 Γ ( α ) ∫ a x y ( ξ ) ( x − ξ ) 1 − α d ξ , for x > a , where Γ denotes the Gamma function. Whereas, the right-sided fractional integral is of the form I b − α y ( x ) = 1 Γ ( α ) ∫ x b y ( ξ ) ( ξ − x ) 1 − α d ξ , for x < b . The Riesz fractional integral (known also as the Riesz potential), in general, is defined in the infinite n-dimensional domain (Rubin, 1996; Kilbas et al., 2006). Here, we limit our considerations to its one-dimensional definition in the finite interval of the real axis, and then it can be regarded as the linear combination of the left- and right-sided Riemann–Liouville integrals. The Riesz fractional integral of order α > 0, α ≠ 1 , 3 , 5 , … of the function y ( x ) on the finite interval [ a , b ] is defined as (see: Kilbas et al., 2006) R I [ a , b ] α y ( x ) = 1 2 Γ ( α ) cos ( α π / 2 ) ∫ a b y ( ξ ) | ξ − x | 1 − α d ξ = 1 2 cos ( α π / 2 ) ( I a + α y ( x ) + I b − α y ( x ) ) , for a < x < b .

The above fractional integral operators play an important role in the fractional calculus, especially to transform fractional differential equations into fractional integral equations. Therefore, among others, it is necessary to develop appropriate numerical methods for approximate evaluation of these fractional integrals. In view of the numerous applications of the fractional integral operators, there is great demand for efficient and accurate algorithms for their numerical calculations, especially for the integrand functions that have complicated forms or the explicit analytical forms of the fractional integrals that are so far unknown. Generally, different ways of discretization of the integrand function yield a series of quadrature formulas, whereas different formulas of coefficients in these quadratures give different accuracies. Numerical approximations of the fractional integrals are most often based on polynomial interpolation from which numerical schemes can be derived with a specified accuracy. Among the pioneering works in this field, the book by Oldham and Spanier (1974) can be distinguished. The reviews of different numerical methods (so far known) for fractional integrals and derivatives can be found in Cai and Li (2020), Li and Zeng (2015), Almeida et al. (2015), Blaszczyk and Siedlecki (2014), Baleanu et al. (2012), Malinowska et al. (2015), Blaszczyk et al. (2018), Odibat (2006), Dimitrov (2021), Budak et al. (2023). Due to the wide applications of fractional calculus, it is becoming increasingly important to develop numerical algorithms with high accuracy, fast convergence and less storage memory.

In this work, we focused primarily on the study of numerical algorithms for approximation of fractional integrals (1)–(3) by applying several interpolation methods for the integrand function using different kinds of splines. Next, we derived and presented general numerical integration schemes for the left and right-sided Riemann–Liouville and the Riesz fractional integrals. For each numerical scheme, we estimated the computational error using analytical methods. We tested the quality of the obtained numerical algorithms on three examples by examining numerical errors and determining the experimental order of convergence.

 Numerical Algorithms

The integrand function y ( x ) given in Eqs. (1)–(3) should be defined on interval [ a , b ] and sufficiently smooth. In the proposed approach, function y ( x ) is replaced by same interpolation formulas, e.g. an interpolation polynomial of arbitrary degree or an interpolation spline, as y ( x ) ≅ s ( x ) . Often instead of using high-degree polynomials, one can apply the spline curves (Press et al., 2007; Engeln-Müllges and Uhlig, 1996). The advantage of such an approach is that the spline curves are at most third degree polynomials and hence the oscillations between interpolation points are characterized by smaller amplitude and, after all, the splines are only piecewise continuous.

Spline Interpolations

The spline is a set of piecewise polynomials linked in the set of points. Let us assume that the interval [ a , b ] is divided into N sub-intervals [ x i , x i + 1 ], for i = 0 , 1 , … , N − 1, with the constant step length Δ x = ( b − a ) / N. The coordinates of N + 1 nodal points are the following: a = x 0 < x 1 < ⋯ < x N = b, x i = a + i Δ x. At the given set of points, the values of function y ( x ) are determined as y i = y ( x i ), for i = 0 , 1 , … , N. The piecewise function s ( x ) can be expressed as s ( x ) = s 0 ( x ) , if x ∈ [ x 0 , x 1 ] , s 1 ( x ) , if x ∈ [ x 1 , x 2 ] , … s N − 1 ( x ) , if x ∈ [ x N − 1 , x N ] , where s i ( x ) are polynomials of degree p in each sub-interval [ x i , x i + 1 ], for i = 0 , 1 , … , N − 1, and here we write these polynomials in the form s i ( x ) = ∑ k = 0 p c k , i ( x − x i ) k , where c k , i , for k = 0 , 1 , … , p, are the coefficients of polynomial s i ( x ) in i-th sub-interval [ x i , x i + 1 ]. The way of determination of these coefficients depends on the kind of spline interpolation used. Additionally, the important feature of the spline s ( x ) that interpolates the set of the data points ( x 0 , y 0 ) , ( x 1 , y 1 ) , … , ( x N , y N ) is the relationship s ( x i ) = y i , for i = 0 , 1 , … , N.

Linear Spline Interpolation

The linear interpolation is the simplest form of interpolation where the set of the data points is approximated by a piecewise linear function (the first degree polynomial), which means that the adjacent data points ( x i , y i ) and ( x i + 1 , y i + 1 ), for i = 0 , 1 , … , N − 1, wherein N ⩾ 2, are connected by straight lines. The polynomial of first degree ( p = 1) in i-th sub-interval is constructed as s i ( x ) = c 0 , i + c 1 , i ( x − x i ) .

The linear Lagrange interpolating polynomial in the interval x ∈ [ x i , x i + 1 ], for i = 0 , 1 , … , N − 1, that passes through the points ( x i , y i ) and ( x i + 1 , y i + 1 ) results from the relation y ( x ) = x − x i + 1 x i − x i + 1 y ( x i ) + x − x i x i + 1 − x i y ( x i + 1 ) + y ″ ( x ¯ i ) 2 ! ( x − x i ) ( x − x i + 1 ) = y i + y i + 1 − y i Δ x ( x − x i ) + y ″ ( x ¯ i ) 2 ( x − x i ) ( x − x i + 1 ) = s i ( x ) + Err 1 i ( x ) , and hence, the polynomial coefficients in Eq. (7) are as follows c 0 , i = y i , c 1 , i = y i + 1 − y i Δ x . Moreover, the error term in i-th interval is given as Err 1 i ( x ) = y ″ ( x ¯ i ) 2 ( x − x i ) ( x − x i + 1 ) , where x ¯ i ∈ [ x i , x i + 1 ] is a certain point in the interval. The error term involves y ″ , so the linear spline gives the exact result for any function y ( x ), where the second derivative is equal to zero. In literature, this approach used for integration is known as the trapezoidal rule (Press et al., 2007).

 Quadratic Spline Interpolation

In this kind of approximation, the second degree polynomial ( p = 2) in each i-th sub-interval is used s i ( x ) = c 0 , i + c 1 , i ( x − x i ) + c 2 , i ( x − x i ) 2 , where x ∈ [ x i , x i + 1 ], for i = 0 , 1 , … , N − 1, and N ⩾ 2 must also be even. To determine the coefficients of parabola (quadratic segment), three data points are needed, and hence each parabola spans two sub-intervals.

Here, we suppose that function y ( x ) in the considered double interval [ x i − 1 , x i + 1 ], for i = 1 , 3 , … , N − 1, is expanded in the Taylor polynomial of third order about the central point x i as y ( x ) = y ( x i ) + y ′ ( x i ) ( x − x i ) + y ″ ( x i ) 2 ! ( x − x i ) 2 + y ‴ ( x i ) 3 ! ( x − x i ) 3 + y ( 4 ) ( x ¯ 1 i ) 4 ! ( x − x i ) 4 , where x ¯ 1 i ∈ [ x i − 1 , x i + 1 ] is a certain point in the interval. The values of y ′ ( x i ) and y ″ ( x i ) in Eq. (12) are replaced by well-known relationships y ′ ( x i ) = 1 2 Δ x ( y ( x i + 1 ) − y ( x i − 1 ) ) − ( Δ x ) 2 6 y ‴ ( x ¯ 2 i ) , y ″ ( x i ) = 1 ( Δ x ) 2 ( y ( x i + 1 ) − 2 y ( x i ) + y ( x i − 1 ) ) − ( Δ x ) 2 12 y ( 4 ) ( x ¯ 3 i ) , where x ¯ 2 i ∈ [ x i − 1 , x i + 1 ] and x ¯ 3 i ∈ [ x i − 1 , x i + 1 ], and hence it was obtained y ( x ) = y ( x i ) + ( y ( x i + 1 ) − y ( x i − 1 ) 2 Δ x − ( Δ x ) 2 6 y ‴ ( x ¯ 2 i ) ) ( x − x i ) + 1 2 ! ( y ( x i + 1 ) − 2 y ( x i ) + y ( x i − 1 ) ( Δ x ) 2 − ( Δ x ) 2 12 y ( 4 ) ( x ¯ 3 i ) ) ( x − x i ) 2 + y ‴ ( x i ) 3 ! ( x − x i ) 3 + y ( 4 ) ( x ¯ 1 i ) 4 ! ( x − x i ) 4 , or written in the form y ( x ) = s i ( x ) + Err 2 i ( x ) , for x ∈ [ x i − 1 , x i + 1 ] , where s i ( x ) = y i + y i + 1 − y i − 1 2 Δ x ( x − x i ) + y i + 1 − 2 y i + y i − 1 2 ( Δ x ) 2 ( x − x i ) 2 and Err 2 i ( x ) = 1 6 y ‴ ( x i ) ( x − x i ) 3 − ( Δ x ) 2 6 y ‴ ( x ¯ 2 i ) ( x − x i ) + 1 24 y ( 4 ) ( x ¯ 1 i ) ( x − x i ) 4 − ( Δ x ) 2 24 y ( 4 ) ( x ¯ 3 i ) ( x − x i ) 2 ≅ 1 6 ( ( x − x i ) 3 − ( Δ x ) 2 ( x − x i ) ) y ‴ ( x ¯ i ) + 1 24 ( ( x − x i ) 4 − ( Δ x ) 2 ( x − x i ) 2 ) y ( 4 ) ( x ¯ i ) , wherein x ¯ i ∈ [ x i − 1 , x i + 1 ] is taking into account as common values of x i , x ¯ 1 i , x ¯ 2 i and x ¯ 3 i .

In order to maintain a uniform form of the spline notation (as in Eq. (11)), two cases are considered c k , i = c ¯ ¯ k , i , for i even , c ¯ k , i , for i odd , for k = 0 , 1 , 2 .

If i is odd (i.e. i = 1 , 3 , … , N − 1), then the coefficients c ¯ k , i can be taken directly from Eq. (16) as c ¯ 0 , i = y i , c ¯ 1 , i = y i + 1 − y i − 1 2 Δ x , c ¯ 2 , i = y i − 1 − 2 y i + y i + 1 2 ( Δ x ) 2 , and then the relationships s i ( x i − 1 ) = y i − 1 , s i ( x i ) = y i , s i ( x i + 1 ) = y i + 1 occur.

If i is even (i.e. i = 0 , 2 , … , N − 2), then the coefficients c ¯ ¯ k , i are calculated from the transformation of Eq. (16) by substitution i ← i + 1 on the right side of equation, and hence they have the following forms c ¯ ¯ 0 , i = y i , c ¯ ¯ 1 , i = − 3 y i + 4 y i + 1 − y i + 2 2 Δ x , c ¯ ¯ 2 , i = y i − 2 y i + 1 + y i + 2 2 ( Δ x ) 2 , and the properties s i ( x i ) = y i , s i ( x i + 1 ) = y i + 1 , s i ( x i + 2 ) = y i + 2 hold.

As can be seen in the first case, the quadratic segment is defined in the interval [ x i , x i + 2 ] (for i even), while in the second case it is defined in the interval [ x i − 1 , x i + 1 ] (for i odd), but in fact, both quadratic segments are identical in two adjacent sub-intervals. This approach will simplify the notations of numerical integration later in the work.

 Cubic Spline Interpolation

The cubic splines produce a curve that appears to be seamless and has smooth characteristics. A piecewise continuous curve passes through each of the data points ( x i , y i ), for i = 0 , 1 , … , N, wherein N ⩾ 4, in the given order and the separate polynomial of third degree ( p = 3) (so-called the cubic polynomial segment) in each sub-interval has its own set of coefficients, i.e. s i ( x ) = c 0 , i + c 1 , i ( x − x i ) + c 2 , i ( x − x i ) 2 + c 3 , i ( x − x i ) 3 , where x ∈ [ x i , x i + 1 ], for i = 0 , 1 , … , N − 1. These four coefficients c 0 , i , c 1 , i , c 2 , i and c 3 , i for i-th polynomial, in each of the N sub-intervals should be determined, and hence, in order to define the whole spline, a total of 4 N independent dependencies imposed on the spline are required.

In the case of this kind of spline, the following assumptions are made (Burden et al., 2016):

the polynomials match the data points at both ends of each i-th sub-interval: s i ( x i ) = y i , s i ( x i + 1 ) = y i + 1 , for i = 0 , 1 , … , N − 1, (two conditions for each sub-interval give 2 N dependencies in total).

the spline interpolation must be smooth, therefore the following requirements are assumed for the first and second derivatives of the spline (which means that the slope and the curvature must be equal for each pair of neighbouring polynomials that join at each data point): s i − 1 ′ ( x i ) = s i ′ ( x i ), s i − 1 ″ ( x i ) = s i ″ ( x i ), for i = 1 , 2 , … , N − 1, (two conditions for each internal node give further 2 N − 2 dependencies).

the missing two additional conditions should be determined on the basis of the assumed integer order derivatives of the function y ( x) in the boundary nodes x 0 and x N (colloquially called the end point conditions). Three variants of the cubic spline constructions are considered:

variant 1: s 0 ′ ( x 0 ) = Y 0 ′ , s N − 1 ′ ( x N ) = Y N ′ ,

variant 2: s 0 ″ ( x 0 ) = Y 0 ″ , s N − 1 ″ ( x N ) = Y N ″ ,

variant 3: s 0 ‴ ( x 0 ) = Y 0 ‴ , s N − 1 ‴ ( x N ) = Y N ‴ ,

 where Y 0 ′ , Y 0 ″ , Y 0 ‴ , Y N ′ , Y N ″ and Y N ‴ correspond to the values of the first, second and third order derivatives of the function y ( x ) at the nodes x 0 and x N , respectively. In principle, any combination of derivative orders at both boundary nodes can be considered (e.g. s 0 ‴ ( x 0 ) and s N − 1 ′ ( x N )).

To determine the coefficients in Eq. (21), we used substitutions (see: Burden et al., 2016) in order to reduce the number of coefficients to be determined c 0 , i = y i , c 1 , i = y i + 1 − y i Δ x − Δ x 3 ( c 2 , i + 1 + 2 c 2 , i ) , c 3 , i = 1 3 Δ x ( c 2 , i + 1 − c 2 , i ) , and hence, the cubic spline in the i-th sub-interval (Eq. (21)) can be written in the form s i ( x ) = y i + ( y i + 1 − y i Δ x − Δ x 3 ( c 2 , i + 1 + 2 c 2 , i ) ) ( x − x i ) + c 2 , i ( x − x i ) 2 + 1 3 Δ x ( c 2 , i + 1 − c 2 , i ) ( x − x i ) 3 .

It can now be seen that only N + 1 unknown coefficients c 2 , i need to be determined (the last one with the index N is used to evaluate c 1 , N − 1 and c 3 , N − 1 ). The first, second and third derivatives of the spline s i ( x ) are as follows s i ′ ( x ) = ( y i + 1 − y i Δ x − Δ x 3 ( c 2 , i + 1 + 2 c 2 , i ) ) + 2 c 2 , i ( x − x i ) s i ′ ( x ) = + 1 Δ x ( c 2 , i + 1 − c 2 , i ) ( x − x i ) 2 , s i ″ ( x ) = 2 c 2 , i + 2 Δ x ( c 2 , i + 1 − c 2 , i ) ( x − x i ) , s i ‴ ( x ) = 2 Δ x ( c 2 , i + 1 − c 2 , i ) .

It is easy to check that the conditions: s i ( x i ) = y i , s i ( x i + 1 ) = y i + 1 , for i = 0 , 1 , … , N − 1, and s i − 1 ″ ( x i ) = s i ″ ( x i ), for i = 1 , 2 , … , N − 1, are always present. While, from the conditions s i − 1 ′ ( x i ) = s i ′ ( x i ), for i = 1 , 2 , … , N − 1, as a result, the following dependencies are obtained c 2 , i + 1 + 4 c 2 , i + c 2 , i − 1 = 3 y i + 1 − 2 y i + y i − 1 ( Δ x ) 2 , for i = 1 , 2 , … , N − 1 .

To complete the system of N + 1 linear equations, two dependencies resulting from the end point conditions should be added. If the points x 0 and x N are inserted into the appropriate spline derivatives, then we get s 0 ′ ( x 0 ) = y 1 − y 0 Δ x − Δ x 3 ( c 2 , 1 + 2 c 2 , 0 ) , s 0 ″ ( x 0 ) = 2 c 2 , 0 , s 0 ‴ ( x 0 ) = 2 Δ x ( c 2 , 1 − c 2 , 0 ) , s N − 1 ′ ( x N ) = y N − y N − 1 Δ x + Δ x 3 ( c 2 , N − 1 + 2 c 2 , N ) , s N − 1 ″ ( x N ) = 2 c 2 , N , s N − 1 ‴ ( x N ) = 2 Δ x ( c 2 , N − c 2 , N − 1 ) .

Hence, depending on the assumed variant of the cubic spline construction, the missing end point conditions can be determined in the forms

variant 1: 2 c 2 , 0 + c 2 , 1 = 3 ( y 1 − y 0 ( Δ x ) 2 − Y ′ 0 Δ x ),

c 2 , N − 1 + 2 c 2 , N = 3 ( Y ′ N Δ x − y N − y N − 1 ( Δ x ) 2 ),

variant 2: c 2 , 0 = 1 2 Y 0 ″ ,  c 2 , N = 1 2 Y N ″ ,

variant 3: c 2 , 0 − c 2 , 1 = − Δ x 2 Y 0 ‴ ,  − c 2 , N − 1 + c 2 , N = Δ x 2 Y N ‴ .

The above derived relationships create the linear system of equations that can be written in matrix form γ 0 , 0 γ 0 , 1 0 0 ⋯ 0 0 0 1 4 1 0 0 0 0 0 1 4 1 0 0 0 0 0 1 4 0 0 0 ⋮ ⋱ ⋮ 0 0 0 0 4 1 0 0 0 0 0 1 4 1 0 0 0 0 ⋯ 0 γ N , N − 1 γ N , N · c 2 , 0 c 2 , 1 c 2 , 2 c 2 , 3 ⋮ c 2 , N − 2 c 2 , N − 1 c 2 , N = d 0 d 1 d 2 d 3 ⋮ d N − 2 d N − 1 d N , where d i = 3 y i + 1 − 2 y i + y i − 1 ( Δ x ) 2 , for i = 1 , 2 , … , N − 1 and

for variant 1: γ 0 , 0 = 2,  γ 0 , 1 = 1,  d 0 = 3 ( y 1 − y 0 ( Δ x ) 2 − Y 0 ′ Δ x ),

γ N , N − 1 = 1,  γ N , N = 2,  d N = 3 ( Y N ′ Δ x − y N − y N − 1 ( Δ x ) 2 ),

for variant 2: γ 0 , 0 = 1,  γ 0 , 1 = 0,  d 0 = 1 2 Y 0 ″ ,

γ N , N − 1 = 0,   γ N , N = 1,   d N = 1 2 Y N ″ ,

for variant 3: γ 0 , 0 = 1,  γ 0 , 1 = − 1,  d 0 = − Δ x 2 Y 0 ‴ ,

γ N , N − 1 = − 1,  γ N , N = 1,  d N = Δ x 2 Y N ‴ .

The values of Y 0 ′ , Y 0 ″ , Y 0 ‴ , Y N ′ , Y N ″ and Y N ‴ can be determined directly if the first, second and third order derivatives of function y ( x ) can be derived analytically, i.e. Y 0 ′ = y ′ ( x 0 ), Y N ′ = y ′ ( x N ) and so on. Otherwise, we propose to determine the numerical values for them using the forward/backward finite difference schemes of fourth order accuracy (with the uniform grid spacing) in the forms (Fornberg, 1988) Y 0 ′ ≅ 1 Δ x ( − 25 12 y 0 + 4 y 1 − 3 y 2 + 4 3 y 3 − 1 4 y 4 ) , Y N ′ ≅ 1 Δ x ( 25 12 y N − 4 y N − 1 + 3 y N − 2 − 4 3 y N − 3 + 1 4 y N − 4 ) , Y 0 ″ ≅ 1 ( Δ x ) 2 ( 15 4 y 0 − 77 6 y 1 + 107 6 y 2 − 13 y 3 + 61 12 y 4 − 5 6 y 5 ) , Y N ″ ≅ 1 ( Δ x ) 2 ( 15 4 y N − 77 6 y N − 1 + 107 6 y N − 2 − 13 y N − 3 + 61 12 y N − 4 − 5 6 y N − 5 ) , Y 0 ‴ ≅ 1 ( Δ x ) 3 ( − 49 8 y 0 + 29 y 1 − 461 8 y 2 + 62 y 3 − 307 8 y 4 + 13 y 5 − 15 8 y 6 ) , Y N ‴ ≅ 1 ( Δ x ) 3 ( 49 8 y N − 29 y N − 1 + 461 8 y N − 2 − 62 y N − 3 + 307 8 y N − 4 Y N ‴ ≅ − 13 y N − 5 + 15 8 y N − 6 ) .

By solving the system of equations (27), we obtain the values of polynomial coefficients c 2 , i , for i = 0 , 1 , … , N, in Eq. (21), while the values of remaining coefficients c 0 , i , c 1 , i , c 3 , i , for i = 0 , 1 , … , N − 1, can then be directly calculated from relations (23).

Remarks to the cubic splines regarding the global approximation: in order to determine the polynomial coefficients c k , i in Eq. (21), the value of the interpolation function s ( x ) at an arbitrary point x depends on all data points ( x i , y i ), for i = 0 , 1 , … , N. Here, the variation or perturbation of one arbitrary data point ( x i , y i ) affects the construction of the whole spline function s ( x ).

In numerous literature on the spline interpolation, it can be read that if the exact values of the first derivative of function y ( x ) in both boundary nodes are known, then one can build the spline with exact boundary conditions (colloquially called a ‘clamped spline’) – this corresponds to variant 1. Here, we want to extend this nomenclature to the remaining considered variants. Whereas, if the values of the derivatives of function y ( x ) in the boundary nodes x 0 and x N are unknown, then one can set the so-called natural boundary conditions (i.e. the second derivatives in these nodes take values equal to 0) – colloquially called a ‘natural spline’. The natural spline is generally characterized by higher approximation errors than the clamped spline, and for this reason, it was omitted in this work.

Based on the theorem (e.g. Burden et al., 2016), the interpolation error for the clamped cubic spline can be determined. If y ( x ) ∈ C 4 [ a , b ] and s ( x ) is the clamped cubic spline that interpolates function y ( x ) with respect to the nodes x i , for i = 0 , 1 , … , N, then the following relationship occurs for all x ∈ [ a , b ] Err 3 = | y ( x ) − s ( x ) | ⩽ 5 384 max a ⩽ x ⩽ b | y ( 4 ) ( x ) | ( Δ x ) 4 .

Numerical Integration

According to the previously adopted assumptions, we propose an approach to determine the approximate values of the left- and right-sided fractional integrals (1)–(2) of function y ( x ) in the set of data points x = g ∈ { x M }, x M = a + M Δ x, where M = 1 , … , N for integral (1), and M = 0 , … , N − 1 for integral (2), respectively, using the spline function s ( x ) in the following ways I a + α y ( x ) x = g ≅ I a + α s ( x ) x = g = x M = 1 Γ ( α ) ∫ a g s ( ξ ) ( g − ξ ) 1 − α d ξ = ∑ i = 0 M − 1 1 Γ ( α ) ∫ x i x i + 1 s i ( ξ ) ( x M − ξ ) 1 − α d ξ , for M = 1 , … , N , I b − α y ( x ) x = g ≅ I b − α s ( x ) x = g = x M = 1 Γ ( α ) ∫ g b s ( ξ ) ( ξ − g ) 1 − α d ξ = ∑ i = M N − 1 1 Γ ( α ) ∫ x i x i + 1 s i ( ξ ) ( ξ − x M ) 1 − α d ξ , for M = 0 , … , N − 1 . The fractional integrals containing the spline function s ( x ) are the sum of the local integrals.

If we use the notation: x i = a + i Δ x, then the general form of (6) is of the form s i ( x ) = ∑ k = 0 p c k , i ( x − ( a + i Δ x ) ) k and if we put Eq. (35) into Eq. (33), then we obtain I a + α s ( x ) x = x M = ∑ i = 0 M − 1 1 Γ ( α ) ∫ a + i Δ x a + ( i + 1 ) Δ x ∑ k = 0 p c k , i ( ξ − ( a + i Δ x ) ) k ( a + M Δ x − ξ ) 1 − α d ξ = ∑ i = 0 M − 1 ∑ k = 0 p c k , i ( 1 Γ ( α ) ∫ a + i Δ x a + ( i + 1 ) Δ x ( ξ − ( a + i Δ x ) ) k ( a + M Δ x − ξ ) 1 − α d ξ ) or in the form I a + α s ( x ) x = x M = ∑ i = 0 M − 1 ∑ k = 0 p c k , i J i , M L , k , where the particular integrals J i , M L , k , for k = 0 , 1 , … , 3, with regard to the data point x M in i-th sub-interval ( i < M ), one can find in an analytical way, and they are as follows: J i , M L , k = 1 Γ ( α ) ∫ a + i Δ x a + ( i + 1 ) Δ x ( ξ − ( a + i Δ x ) ) k ( a + M Δ x − ξ ) 1 − α d ξ = ξ = a + ( u + i ) Δ x ( Δ x ) α + k Γ ( α ) ∫ 0 1 u k ( M − i − u ) 1 − α d u J i , M L , k = ( Δ x ) α + k Γ ( α + k + 1 ) ( k ! ( M − i ) α + k − d i , M L , k ( M − i − 1 ) α ) , where d i , M L , k = 1 , if k = 0 , ( M − i − 1 ) + ( α + 1 ) , if k = 1 , 2 ( M − i − 1 ) 2 + 2 ( α + 2 ) ( M − i − 1 ) + ( α + 1 ) ( α + 2 ) , if k = 2 , 6 ( M − i − 1 ) 3 + 6 ( α + 3 ) ( M − i − 1 ) 2 + 3 ( α + 2 ) ( α + 3 ) ( M − i − 1 ) + ( α + 1 ) ( α + 2 ) ( α + 3 ) , if k = 3 .

In a similar manner, we proceed with the second integral (34) and we have I b − α s ( x ) | x = x M = ∑ i = M N − 1 1 Γ ( α ) ∫ a + i Δ x a + ( i + 1 ) Δ x ∑ k = 0 p c k , i ( ξ − ( a + i Δ x ) ) k ( ξ − ( a + M Δ x ) ) 1 − α d ξ = ∑ i = M N − 1 ∑ k = 0 p c k , i ( 1 Γ ( α ) ∫ a + i Δ x a + ( i + 1 ) Δ x ( ξ − ( a + i Δ x ) ) k ( ξ − ( a + M Δ x ) ) 1 − α d ξ ) or I b − α s ( x ) | x = x M = ∑ i = M N − 1 ∑ k = 0 p c k , i J i , M R , k , and the integrals J i , M R , k , for k = 0 , 1 , … , 3, also with respect to data point x M in i-th sub-interval ( i ⩾ M ), take the following analytical forms J i , M R , k = 1 Γ ( α ) ∫ a + i Δ x a + ( i + 1 ) Δ x ( ξ − ( a + i Δ x ) ) k ( ξ − ( a + M Δ x ) ) 1 − α d ξ = ξ = a + ( u + i ) Δ x ( Δ x ) k + α Γ ( α ) ∫ 0 1 u k ( i − M + u ) 1 − α d u J i , M R , k = ( Δ x ) α + k Γ ( α + k + 1 ) ( ( − 1 ) k + 1 k ! ( i − M ) α + k + d i , M R , k ( i − M + 1 ) α ) , where d i , M R , k = 1 , if k = 0 , − ( i − M + 1 ) + ( α + 1 ) , if k = 1 , 2 ( i − M + 1 ) 2 − 2 ( α + 2 ) ( i − M + 1 ) + ( α + 2 ) ( α + 1 ) , if k = 2 , − 6 ( i − M + 1 ) 3 + 6 ( α + 3 ) ( i − M + 1 ) 2 − 3 ( α + 2 ) ( α + 3 ) ( i − M + 1 ) + ( α + 3 ) ( α + 2 ) ( α + 1 ) , if k = 3 .

Using the above numerical integration schemes, one can derive the numerical algorithm for the Riesz fractional integral (3) in the following way R I [ a , b ] α y ( x ) | x = g = 1 2 cos ( α π / 2 ) ( I a + α y ( x ) | x = g + I b − α y ( x ) | x = g ) ≅ 1 2 cos ( α π / 2 ) ( I a + α s ( x ) | x = g = x M + I b − α s ( x ) | x = g = x M ) = 1 2 cos ( α π / 2 ) ( ∑ i = 0 M − 1 ∑ k = 0 p c k , i J i , M L , k + ∑ i = M N − 1 ∑ k = 0 p c k , i J i , M R , k ) = 1 2 cos ( α π / 2 ) ∑ i = 0 N − 1 ∑ k = 0 p c k , i J i , M L , k , if i = 0 , … , M − 1 , J i , M R , k , if i = M , … , N − 1 , for M = 1 , … , N − 1.

 Error Estimate for the Composite Scheme of Integration

To determine the global error of fractional integration, we use the previously determined interpolation errors on the individual intervals for each kind of spline, which should be integrated and summed. Here, we limit to the estimation of errors for the left-sided fractional Riemann–Liouville integral, but for the remaining integral operators, the results will be analogous. Additionally, to simplify calculations, the whole (largest) integration interval [ a , b] (i.e. assuming M = N) is taken into account.

The Case of the Linear Spline

Taking the particular error terms Err 1 i (10), for i = 0 , 1 , … , N − 1, we integrate them over each interval separately using the properties of integral calculus Err = ∑ i = 0 N − 1 1 Γ ( α ) ∫ x i x i + 1 Err 1 i ( ξ ) ( x N − ξ ) 1 − α d ξ = ∑ i = 0 N − 1 1 Γ ( α ) ∫ x i x i + 1 1 2 ( ξ − x i ) ( ξ − x i + 1 ) y ″ ( x ¯ i ) 1 ( x N − ξ ) 1 − α d ξ ⩽ 1 2 max i = 0 , 1 , … , N − 1 | y ″ ( x ¯ i ) | ∑ i = 0 N − 1 1 Γ ( α ) ∫ x i x i + 1 ( ξ − x i ) ( ξ − x i + 1 ) ( x N − ξ ) 1 − α d ξ = 1 2 M 2 θ 2 N ( α , Δ x ) , where M 2 = max i = 0 , 1 , … , N − 1 | y ″ ( x ¯ i ) | , θ 2 N ( α , Δ x ) = ∑ i = 0 N − 1 1 Γ ( α ) ∫ x i x i + 1 ( ξ − x i ) ( ξ − x i + 1 ) ( x N − ξ ) 1 − α d ξ = ξ = x i + u Δ x ∑ i = 0 N − 1 ( Δ x ) 2 + α Γ ( α ) ∫ 0 1 u ( u − 1 ) ( N − i − u ) 1 − α d u . After analytical integration and simplifications we obtain θ 2 N ( α , Δ x ) = ( Δ x ) 2 + α ( − 2 N 3 + α Γ ( 3 + α ) + 1 Γ ( 2 + α ) ( 2 ∑ i = 1 N i 1 + α − N 1 + α ) ) . We approximate the finite sum ∑ i = 1 N i 1 + α by using the Euler-Maclaurin formula of the second order. Then we get θ 2 N ( α , Δ x ) = ( Δ x ) 2 + α ( C 2 ( α ) + 1 6 Γ ( 1 + α ) N α ) , where C 2 ( α ) = − 2 Γ ( 3 + α ) + 1 Γ ( 2 + α ) − 1 6 Γ ( 1 + α ) . It can be noticed that for α = 1, we have C 2 ( 1 ) = 0, θ 2 N ( 1 , Δ x ) = ( Δ x ) 3 N / 6 and we obtain compliance with the classical trapezoidal rule of integration. By substituting Eq. (48) into Eq. (45) and taking into account N = ( b − a ) / Δ x, the final estimation of the composed error is as follows Err ⩽ 1 2 M 2 ( ( Δ x ) 2 + α ( C 2 ( α ) + 1 6 Γ ( 1 + α ) ( b − a Δ x ) α ) ) = M 2 ( C 2 ( α ) 2 ( Δ x ) 2 + α + ( b − a ) α 12 Γ ( 1 + α ) ( Δ x ) 2 ) = M 2 ( b − a ) α 12 Γ ( 1 + α ) ( Δ x ) 2 + O ( ( Δ x ) 2 + α ) .

The Case of the Quadratic Spline

Here, the error terms Err 2 i (17), for i = 1 , 3 , … , N − 1, were determined on double adjacent intervals, hence their fractional integration concerns double intervals Err = ∑ i = 1 , 3 , … , N − 1 1 Γ ( α ) ∫ x i − 1 x i + 1 Err 2 i ( ξ ) ( x N − ξ ) 1 − α d ξ = ∑ i = 1 , 3 , … , N − 1 1 Γ ( α ) ∫ x i − 1 x i + 1 1 6 ( ( ξ − x i ) 3 − ( Δ x ) 2 ( ξ − x i ) ) y ‴ ( x ¯ i ) + 1 24 ( ( ξ − x i ) 4 − ( Δ x ) 2 ( ξ − x i ) 2 ) y ( 4 ) ( x ¯ i ) 1 ( x N − ξ ) 1 − α d ξ ⩽ 1 6 max i = 1 , 3 , … , N − 1 | y ‴ ( x ¯ i ) | × ∑ i = 1 , 3 , … , N − 1 1 Γ ( α ) ∫ x i − 1 x i + 1 ( ξ − x i ) 3 − ( Δ x ) 2 ( ξ − x i ) ( x N − ξ ) 1 − α d ξ + 1 24 max i = 1 , 3 , … , N − 1 | y ( 4 ) ( x ¯ i ) | × ∑ i = 1 , 3 , … , N − 1 1 Γ ( α ) ∫ x i − 1 x i + 1 ( ξ − x i ) 4 − ( Δ x ) 2 ( ξ − x i ) 2 ( x N − ξ ) 1 − α d ξ = 1 6 M 3 θ 3 N ( α , Δ x ) + 1 24 M 4 θ 4 N ( α , Δ x ) , where M 3 = max i = 1 , 3 , … , N − 1 | y ‴ ( x ¯ i ) | , M 4 = max i = 1 , 3 , … , N − 1 | y ( 4 ) ( x ¯ i ) | , θ 3 N ( α , Δ x ) = ∑ i = 1 , 3 , … , N − 1 1 Γ ( α ) ∫ x i − 1 x i + 1 ( ξ − x i ) 3 − ( Δ x ) 2 ( ξ − x i ) ( x N − ξ ) 1 − α d ξ = ξ = x i + u Δ x ∑ i = 1 , 3 , … , N − 1 ( Δ x ) 3 + α Γ ( α ) ∫ − 1 1 u 3 − u ( N − i − u ) 1 − α d u , θ 4 N ( α , Δ x ) = ∑ i = 1 , 3 , … , N − 1 1 Γ ( α ) ∫ x i − 1 x i + 1 ( ξ − x i ) 4 − ( Δ x ) 2 ( ξ − x i ) 2 ( x N − ξ ) 1 − α d ξ = ξ = x i + u Δ x ∑ i = 1 , 3 , … , N − 1 ( Δ x ) 4 + α Γ ( α ) ∫ − 1 1 u 4 − u 2 ( N − i − u ) 1 − α d u .

Further transformations (including using analytic integration and changing the indexing order in the series) lead to θ 3 N ( α , Δ x ) = ( Δ x ) 3 + α 6 N 3 + α Γ ( 4 + α ) − 6 Γ ( 3 + α ) ( 2 3 + α ∑ i = 1 N / 2 i 2 + α − N 2 + a ) + 2 N 1 + α Γ ( 2 + α ) , θ 4 N ( α , Δ x ) = ( Δ x ) 4 + α ( 24 N 4 + α Γ ( 5 + α ) − 24 Γ ( 4 + α ) ( 2 4 + α ∑ i = 1 N / 2 i 3 + α − N 3 + a ) + 10 N 2 + α Γ ( 3 + α ) − 2 Γ ( 2 + α ) ( 2 2 + α ∑ i = 1 N / 2 i 1 + α − N 1 + a ) ) .

Next, to estimate the finite sums of the series, here the Euler-Maclaurin formula of the fourth order is used, and finally we have θ 3 N ( α , Δ x ) = ( Δ x ) 3 + α ( C 3 ( α ) + 2 15 Γ ( α ) N α − 1 ) , θ 4 N ( α , Δ x ) ≅ ( Δ x ) 4 + α ( C 4 ( α ) − 2 15 Γ ( α + 1 ) N α + 2 ( α − 1 ) 45 Γ ( α ) N α − 2 ) , where C 3 ( α ) = − 3 · 2 4 + α Γ ( 4 + α ) − 3 · 2 3 + α Γ ( 3 + α ) + 2 2 + α Γ ( 2 + α ) − 2 α 15 Γ ( α ) , C 4 ( α ) = − 6 ( 2 + α ) 2 5 + α Γ ( 5 + α ) + 5 · 2 3 + α Γ ( 3 + α ) − 2 2 + α Γ ( 2 + α ) + 2 1 + α 15 Γ ( 1 + α ) − ( α − 1 ) 2 α 90 Γ ( α ) .

It should be pointed out that for α = 1, one obtains C 3 ( 1 ) = − 2 / 15, C 4 ( 1 ) = 0 and hence θ 3 N ( 1 , Δ x ) = 0, θ 4 N ( 1 , Δ x ) = − 2 ( Δ x ) 5 N / 15 – which corresponds to the classic Simpson’s rule approximation for integral. By substituting θ 3 N ( α , Δ x ) and θ 4 N ( α , Δ x ) from Eq. (55) into Eq. (51) and also assuming that N = ( b − a ) / Δ x one gets Err ⩽ 1 6 M 3 ( C 3 ( α ) ( Δ x ) 3 + α + 2 ( b − a ) α − 1 15 Γ ( α ) ( Δ x ) 4 ) + 1 24 M 4 C 4 ( α ) ( Δ x ) 4 + α − 2 ( b − a ) α 15 Γ ( α + 1 ) ( Δ x ) 4 + 2 ( α − 1 ) ( b − a ) α − 2 45 Γ ( α ) ( Δ x ) 6 . Taking into account the terms containing Δ x with the lowest powers, we have Err = M 3 C 3 ( α ) 6 ( Δ x ) 3 + α + O ( ( Δ x ) 4 ) , for α < 1 , ( M 3 ( b − a ) α − 1 45 Γ ( α ) − M 4 ( b − a ) α 180 Γ ( α + 1 ) ) ( Δ x ) 4 + O ( ( Δ x ) 3 + α ) , for α ⩾ 1 .

The Case of the Cubic Spline

For this case, we estimate the global error without considering the summation of local errors in each interval. Using the known properties of the classical calculus and the fractional integration of the constant, we obtain Err = ‖ I a + α y ( x ) x = b − I a + α s ( x ) x = b ‖ = ‖ I a + α ( y ( x ) − s ( x ) ) x = b ‖ ⩽ I a + α ‖ y ( x ) − s ( x ) ‖ x = b ⩽ ‖ y ( x ) − s ( x ) ‖ I a + α 1 x = b = ‖ y ( x ) − s ( x ) ‖ ( b − a ) α Γ ( 1 + α ) .

By inserting relationship (32) into Eq. (59) we obtain Err ⩽ 5 M 4 384 ( b − a ) α Γ ( 1 + α ) ( Δ x ) 4 , where M 4 = max a ⩽ x ⩽ b | y ( 4 ) ( x ) | or for the discrete form M 4 = max i = 0 , 1 , … , N − 1 | y ( 4 ) ( x ¯ i ) |, x ¯ i ∈ [ x i , x i + 1 ].

 Examples of Computations

In order to verify the correctness of the proposed numerical schemes, sample calculations have been performed. One of the important characteristics of the numerical integration schemes is the determination of the Experimental Order of Convergence (EOC) (Li and Zeng, 2015; Blaszczyk et al., 2018) for each presented scheme.

If the values of fractional integrals at the given data points can be found in an analytical way (i.e. the exact solutions are known), then one can determine the computational error of the numerical integration scheme obtained on the size grid N as err N = I a + α y ( x ) x = g − Ψ N , where Ψ N denotes the numerical/approximate value of I a + α y ( x ) | x = g and simultaneously the exact value of I a + α s ( x ) | x = g . Similarly, the error err N is defined for the right-sided fractional integral and the Riesz integral.

Next, the estimated EOC can be evaluated as EOC N = log 2 | err N / 2 | | err N | and should be performed for the computations over a range of different grid sizes.

Whereas, if the exact analytical result of the definite integral is not known, then the EOC can be determined from the following relationship EOC N = log 2 | Ψ N − Ψ N / 2 | | Ψ 2 N − Ψ N | .

All of the numerical calculations have been performed with the quadruple floating-point precision. Many compilers, as well as numerous mathematical software, support the double precision calculation with the 15–17 significant decimal digits. In the case of the numerical integrations using the spline interpolation, the effect of rounding errors is significant, especially for the precise calculations of the EOC. For this reason, the C++11 in GCC (MinGW-W64) compiler and the quadmath library have been used to build the application. The 128-bit floating point type __float128 for real variables (this type supports the calculation with the 34 significant decimal digits) and the quadmath functions have been applied in the calculations.

We also recall the following properties of the fractional integration of the power functions ( x − a ) γ and ( b − x ) γ (see: Kilbas et al. (2006) ) I a + α ( x − a ) γ = Γ ( γ + 1 ) Γ ( γ + α + 1 ) ( x − a ) γ + α , for γ > − 1 , α > 0 , I b − α ( b − x ) γ = Γ ( γ + 1 ) Γ ( γ + α + 1 ) ( b − x ) γ + α , for γ > − 1 , α > 0 , which are used to determine analytical solutions in the examples below.

The integrand function is the eighth-degree polynomial of the form y ( x ) = x 8 − 8 x 7 + 26 x 6 − 44 x 5 + 40 x 4 − 15 x 3 − 4 x 2 + 5 x + 1 . Numerical values of I a + α y ( x ) x = g = b for a = 0 and b = 2 have been investigated. For the assumed values of a and b, function y ( x ) takes the values y ( a ) = 1, y ( b ) = 3, y ′ ( a ) = 5, y ′ ( b ) = 1, y ″ ( a ) = − 8, y ″ ( b ) = 4.

The choice of polynomial function (66) gives us the possibility to determine the fractional integral of order α > 0 of function (66) for a = 0 at the point x = g = b in the analytical way using (64), which is expressed by the following formula I 0 + α y ( x ) x = b = Γ ( 9 ) Γ ( 9 + α ) b 8 + α − 8 Γ ( 8 ) Γ ( 8 + α ) b 7 + α + 26 Γ ( 7 ) Γ ( 7 + α ) b 6 + α − 44 Γ ( 6 ) Γ ( 6 + α ) b 5 + α + 40 Γ ( 5 ) Γ ( 5 + α ) b 4 + α − 15 Γ ( 4 ) Γ ( 4 + α ) b 3 + α − 4 Γ ( 3 ) Γ ( 3 + α ) b 2 + α + 5 Γ ( 2 ) Γ ( 2 + α ) b 1 + α + Γ ( 1 ) Γ ( 1 + α ) b α , where Γ ( k ) = ( k − 1 ) !, for k = 1 , 2 , … .

In Table 1, we report the numerical errors err N (determined using Eq. (61)) and the calculated values of the EOC (using Eq. (62)) for the selected sets of α ∈ { 0.4 , 0.7 , 1.4 , 2.7 } and N = 100 , 200 , 400 , 800 , 1600 , 3200 , 6400 , 12800.

Results related to Example 1.

α N Linear spline Quadratic spline Cubic spline
Variant 1 Variant 2 Variant 3
err N EOC err N EOC err N EOC err N EOC err N EOC
0.4 100 3.080E−05 – −3.510E−06 – 2.858E−08 – 1.447E−07 – 3.949E−07 –
200 7.018E−06 2.134 −3.700E−07 3.246 4.080E−09 2.808 8.520E−09 4.086 2.026E−08 4.285
400 1.637E−06 2.100 −3.739E−08 3.307 3.188E−10 3.678 4.996E−10 4.092 1.054E−09 4.265
800 3.880E−07 2.076 −3.687E−09 3.342 2.177E−11 3.873 2.957E−11 4.079 5.576E−11 4.240
1600 9.313E−08 2.059 −3.582E−10 3.363 1.418E−12 3.940 1.769E−12 4.063 3.008E−12 4.213
3200 2.256E−08 2.045 −3.449E−11 3.377 9.063E−14 3.968 1.068E−13 4.050 1.654E−13 4.185
6400 5.505E−09 2.035 −3.302E−12 3.385 5.743E−15 3.980 6.497E−15 4.039 9.271E−15 4.157
12800 1.351E−09 2.027 −3.150E−13 3.390 3.622E−16 3.987 3.977E−16 4.030 5.290E−16 4.131
0.7 100 8.235E−05 – −9.581E−07 – 3.687E−08 – 8.270E−08 – 1.814E−07 –
200 2.047E−05 2.008 −8.226E−08 3.542 3.259E−09 3.500 4.657E−09 4.151 8.352E−09 4.441
400 5.102E−06 2.004 −6.841E−09 3.588 2.237E−10 3.865 2.691E−10 4.113 4.084E−10 4.354
800 1.273E−06 2.002 −5.574E−10 3.617 1.440E−11 3.958 1.596E−11 4.075 2.123E−11 4.265
1600 3.180E−07 2.001 −4.479E−11 3.638 9.089E−13 3.985 9.654E−13 4.048 1.165E−12 4.188
3200 7.946E−08 2.001 −3.563E−12 3.652 5.702E−14 3.994 5.912E−14 4.030 6.671E−14 4.127
6400 1.986E−08 2.001 −2.814E−13 3.662 3.570E−15 3.998 3.648E−15 4.018 3.938E−15 4.082
12800 4.963E−09 2.000 −2.211E−14 3.670 2.233E−16 3.999 2.263E−16 4.011 2.373E−16 4.053
1.4 100 1.984E−04 – −6.312E−08 – 2.960E−08 – 4.681e-08 – 8.388E−08 –
200 4.961E−05 2.000 −3.301E−09 4.257 2.226E−09 3.733 2.65509 4.140 3.790E−09 4.468
400 1.240E−05 2.000 −1.772E−10 4.220 1.453E−10 3.937 1.567E−10 4.083 1.917E−10 4.305
800 3.101E−06 2.000 −9.734E−12 4.186 9.178E−12 3.985 9.500E−12 4.044 1.058E−11 4.179
1600 7.752E−07 2.000 −5.460E−13 4.156 5.751E−13 3.996 5.846E−13 4.022 6.182E−13 4.098
3200 1.938E−07 2.000 −3.120E−14 4.129 3.597E−14 3.999 3.626E−14 4.011 3.730E−14 4.051
6400 4.845E−08 2.000 −1.812E−15 4.106 2.248E−15 4.000 2.257E−15 4.006 2.290E−15 4.026
12800 1.211E−08 2.000 −1.068E−16 4.085 1.405E−16 4.000 1.408E−16 4.003 1.418E−16 4.013
2.7 100 2.740E−04 – −1.357E−07 – 3.425E−08 – 5.644E−08 – 1.042E−07 –
200 6.849E−05 2.000 −8.525E−09 3.993 2.620E−09 3.708 3.186E−09 4.147 4.683E−09 4.476
400 1.712E−05 2.000 −5.335E−10 3.998 1.718E−10 3.931 1.871E−10 4.090 2.341E−10 4.322
800 4.281E−06 2.000 −3.335E−11 4.000 1.087E−11 3.983 1.131E−11 4.049 1.278E−11 4.195
1600 1.070E−06 2.000 −2.085E−12 4.000 6.814E−13 3.996 6.944E−13 4.025 7.405E−13 4.109
3200 2.675E−07 2.000 −1.303E−13 4.000 4.262E−14 3.999 4.302E−14 4.013 4.446E−14 4.058
6400 6.689E−08 2.000 −8.144E−15 4.000 2.664E−15 4.000 2.677E−15 4.006 2.722E−15 4.030
12800 1.672E−08 2.000 −5.090E−16 4.000 1.665E−16 4.000 1.669E−16 4.003 1.683E−16 4.015
α Analytical values of I 0 + α y ( x ) | x = 2 calculated using formula (67)
0.4 3.6979129457596915301988815161146608
0.7 4.0856207593403175492511974048448624
1.4 4.3604818404289140653601695680338754
2.7 2.9484099812828967875285769194034989

In the second example, we take a more complicated integrand function of the form y ( x ) = ( x sin ( 3 x 2 ) + 5 x x + 2 ) · exp ( − ( x − 2 ) 3 2 − 2 x ) + x x 8 3 x x 2 + 1 assuming that the analytic form of the left-sided fractional integral is unknown or has a very complicated form. Here, the approximated values of I a + α y ( x ) | x = g = b for a = 1 and b = 4 have been calculated.

In Table 2, the presentation of numerical results Ψ N (where Ψ N = I a + α s ( x ) | x = b obtained on the grid size N) is limited to 15 decimal digits, but all calculations (including the EOC obtained using formula (63)) have been performed for real numbers with the accuracy of 34 significant decimal digits. Identical sets of parameters α and N, as given in Example 1, have been taken into account.

Results related to Example 2.

α N Linear spline Quadratic spline Cubic spline
Variant 1 Variant 2 Variant 3
Ψ N EOC Ψ N EOC Ψ N EOC Ψ N EOC Ψ N EOC
0.4 100 0.129175293650077 – 0.129159283883400 – 0.129159149778395 – 0.129159260371743 – 0.129159333125011 –
200 0.129163284309064 1.973 0.129159195936989 4.003 0.129159190635184 6.329 0.129159191910429 5.277 0.129159192879209 5.720
400 0.129160226303874 1.981 0.129159190452032 3.971 0.129159190127129 3.442 0.129159190145188 4.784 0.129159190219303 4.307
800 0.129159451516164 1.986 0.129159190102355 3.920 0.129159190080380 4.132 0.129159190081100 4.259 0.129159190084947 4.261
1600 0.129159255887306 1.989 0.129159190079257 3.864 0.129159190077714 4.094 0.129159190077754 4.111 0.129159190077939 4.243
3200 0.129159206615184 1.992 0.129159190077671 3.803 0.129159190077558 4.033 0.129159190077560 4.071 0.129159190077569 4.216
6400 0.129159194228093 1.994 0.129159190077557 3.738 0.129159190077548 4.004 0.129159190077549 4.052 0.129159190077549 4.188
12800 0.129159191118241 – 0.129159190077549 – 0.129159190077548 – 0.129159190077548 – 0.129159190077548 –
0.7 100 0.165128141318417 – 0.165103544764964 – 0.165103545293346 – 0.165103583875713 – 0.165103610220807 –
200 0.165109708612410 1.997 0.165103550450332 3.702 0.165103551597517 3.312 0.165103551927078 5.055 0.165103552154690 5.630
400 0.165105091970249 1.998 0.165103550887115 3.870 0.165103550962706 3.948 0.165103550965874 4.443 0.165103550982253 4.292
800 0.165103936437045 1.999 0.165103550916983 3.997 0.165103550921568 4.076 0.165103550921684 4.115 0.165103550922400 4.208
1600 0.165103647338066 1.999 0.165103550918853 4.092 0.165103550919128 4.035 0.165103550919134 4.045 0.165103550919163 4.153
3200 0.165103575029852 2.000 0.165103550918963 4.189 0.165103550918979 4.011 0.165103550918980 4.025 0.165103550918981 4.106
6400 0.165103556947631 2.000 0.165103550918969 4.322 0.165103550918970 4.003 0.165103550918970 4.015 0.165103550918970 4.071
12800 0.165103552426280 – 0.165103550918969 – 0.165103550918970 – 0.165103550918970 – 0.165103550918970 –
1.4 100 0.261754500404957 – 0.261701311442012 – 0.261701464571557 – 0.261701448885203 – 0.261701442194459 –
200 0.261714688712354 2.001 0.261701418707993 4.009 0.261701427658470 4.410 0.261701427393735 3.863 0.261701427109875 3.649
400 0.261704741164795 2.000 0.261701425371550 4.003 0.261701425921852 4.113 0.261701425916636 3.955 0.261701425907174 3.804
800 0.261702254627535 2.000 0.261701425787252 4.001 0.261701425821510 4.028 0.261701425821390 3.980 0.261701425821085 3.907
1600 0.261701633016401 2.000 0.261701425813219 4.000 0.261701425815358 4.007 0.261701425815355 3.990 0.261701425815345 3.954
3200 0.261701477615193 2.000 0.261701425814842 4.000 0.261701425814975 4.002 0.261701425814975 3.995 0.261701425814975 3.977
6400 0.261701438765002 2.000 0.261701425814943 4.000 0.261701425814951 4.000 0.261701425814951 3.998 0.261701425814951 3.989
12800 0.261701429052462 – 0.261701425814949 – 0.261701425814950 – 0.261701425814950 – 0.261701425814950 –
2.7 100 0.351816614974611 – 0.351709046915229 – 0.351709417655548 – 0.351709374905998 – 0.351709355473402 –
200 0.351736132977994 2.001 0.351709305553739 4.004 0.351709327145813 4.448 0.351709326481777 3.792 0.351709325779111 3.398
400 0.351716024501843 2.000 0.351709321670694 4.001 0.351709322999350 4.115 0.351709322986643 3.945 0.351709322963108 3.786
800 0.351710998133673 2.000 0.351709322677255 4.000 0.351709322760002 4.027 0.351709322759708 3.979 0.351709322758956 3.901
1600 0.351709741588551 2.000 0.351709322740154 4.000 0.351709322745321 4.006 0.351709322745313 3.990 0.351709322745290 3.952
3200 0.351709427455202 2.000 0.351709322744085 4.000 0.351709322744408 4.002 0.351709322744407 3.995 0.351709322744407 3.976
6400 0.351709348922048 2.000 0.351709322744330 4.000 0.351709322744351 4.000 0.351709322744351 3.998 0.351709322744351 3.988
12800 0.351709329288771 – 0.351709322744346 – 0.351709322744347 – 0.351709322744347 – 0.351709322744347 –

The next example concerns the determination of numerical values of the Riesz fractional integral. Here, we take into account the polynomial of fifth degree (higher order than the spline of third order) of the form y ( x ) = x 5 − 13 x 4 + 59 x 3 − 108 x 2 + 67 x + 4 .

The values of the Riesz fractional integral of function (69) we designate in the finite interval [ a , b ] at the point x = g (where a ⩽ g ⩽ b). Sample calculations have been carried out for a = 1, b = 5 and g = 2. The point g = 2 corresponds to the node x M , where M = N / 4, and N must be a multiple of 4. In this example, we take the selected sets of α ∈ { 0.25 , 0.75 , 1.25 , 1.75 } and N = 100 , 200 , 400 , 800 , 1600 , 3200 , 6400 , 12800. In order to find the analytical solution of this integral, we first rewrite polynomial (69) by collecting similar terms that match the expressions ( x − a ) and ( b − x ) y ( x ) = ( x − 1 ) 5 − 8 ( x − 1 ) 4 + 17 ( x − 1 ) 3 + ( x − 1 ) 2 − 19 ( x − 1 ) + 10 = − ( 5 − x ) 5 + 12 ( 5 − x ) 4 − 49 ( 5 − x ) 3 + 77 ( 5 − x ) 2 − 37 ( 5 − x ) + 14 , and then the properties (64) and (65) can be directly applied, respectively. Hence, we get the analytical solution R I [ 1 , 5 ] α y ( x ) | x = g = 1 2 cos ( α π / 2 ) Γ ( 6 ) Γ ( 6 + α ) ( ( g − 1 ) 5 + α − ( 5 − g ) 5 + α ) + Γ ( 5 ) Γ ( 5 + α ) ( − 8 ( g − 1 ) 4 + α + 12 ( 5 − g ) 4 + α ) + Γ ( 4 ) Γ ( 4 + α ) ( 17 ( g − 1 ) 3 + α − 49 ( 5 − g ) 3 + α ) + Γ ( 3 ) Γ ( 3 + α ) ( ( g − 1 ) 2 + α + 77 ( 5 − g ) 2 + α ) + Γ ( 2 ) Γ ( 2 + α ) ( − 19 ( g − 1 ) 1 + α − 37 ( 5 − g ) 1 + α ) + Γ ( 1 ) Γ ( 1 + α ) ( 10 ( g − 1 ) α + 14 ( 5 − g ) α ) .

In Table 3, we present the errors err N (where err N = R I [ 1 , 5 ] α y ( x ) | x = 2 − Ψ N , and Ψ N is numerical value of the Riesz integral obtained on the grid size N) and the calculated values of the EOC (Eq. (62)), respectively, for different kinds of the spline interpolation methods.

Results related to Example 3.

α N Linear spline Quadratic spline Cubic spline
Variant 1 Variant 2 Variant 3
err N EOC err N EOC err N EOC err N EOC err N EOC
0.25 100 −2.957E−03 – −1.384E−06 – −1.318E−07 – −1.346E−07 – −1.421E−07 –
200 −7.766E−04 1.929 −1.581E−07 3.130 −8.981E−09 3.876 −9.055E−09 3.893 −9.290E−09 3.935
400 −2.020E−04 1.943 −9.973E−09 3.987 −5.990E−10 3.906 −6.012E−10 3.913 −6.085E−10 3.932
800 −5.214E−05 1.954 −6.280E−10 3.989 −3.941E−11 3.926 −3.947E−11 3.929 −3.970E−11 3.938
1600 −1.338E−05 1.962 −3.950E−11 3.991 −2.566E−12 3.941 −2.568E−12 3.942 −2.575E−12 3.946
3200 −3.418E−06 1.969 −2.482E−12 3.992 −1.658E−13 3.952 −1.658E−13 3.953 −1.661E−13 3.955
6400 −8.698E−07 1.974 −1.558E−13 3.994 −1.065E−14 3.961 −1.065E−14 3.961 −1.065E−14 3.962
12800 −2.207E−07 1.979 −9.774E−15 3.995 −6.803E−16 3.968 −6.803E−16 3.968 −6.806E−16 3.969
0.75 100 −8.977E−03 – −3.265E−06 – 3.319E−07 – 3.356E−07 – 3.575E−07 –
200 −2.251E−03 1.996 −2.213E−07 3.883 2.050E−08 4.017 2.065E−08 4.022 2.134E−08 4.067
400 −5.637E−04 1.997 −1.372E−08 4.011 1.275E−09 4.008 1.280E−09 4.012 1.301E−09 4.035
800 −1.411E−04 1.999 −8.537E−10 4.007 7.945E−11 4.004 7.962E−11 4.007 8.029E−11 4.019
1600 −3.529E−05 1.999 −5.321E−11 4.004 4.958E−12 4.002 4.963E−12 4.004 4.984E−12 4.010
3200 −8.825E−06 1.999 −3.320E−12 4.002 3.096E−13 4.001 3.098E−13 4.002 3.104E−13 4.005
6400 −2.207E−06 2.000 −2.073E−13 4.001 1.934E−14 4.001 1.935E−14 4.001 1.937E−14 4.003
12800 −5.518E−07 2.000 −1.295E−14 4.001 1.208E−15 4.000 1.209E−15 4.001 1.209E−15 4.001
1.25 100 1.125E−02 – 5.353E−06 – −1.499E−06 – −1.558E−06 – −1.765E−06 –
200 2.812E−03 2.000 3.320E−07 4.011 −9.390E−08 3.997 −9.569E−08 4.026 −1.021E−07 4.111
400 7.029E−04 2.000 2.077E−08 3.998 −5.872E−09 3.999 −5.927E−09 4.013 −6.129E−09 4.059
800 1.757E−04 2.000 1.299E−09 3.999 −3.671E−10 4.000 −3.688E−10 4.007 −3.751E−10 4.030
1600 4.393E−05 2.000 8.120E−11 4.000 −2.294E−11 4.000 −2.300E−11 4.003 −2.319E−11 4.015
3200 1.098E−05 2.000 5.075E−12 4.000 −1.434E−12 4.000 −1.436E−12 4.002 −1.442E−12 4.008
6400 2.746E−06 2.000 3.172E−13 4.000 −8.962E−14 4.000 −8.967E−14 4.001 −8.986E−14 4.004
12800 6.864E−07 2.000 1.983E−14 4.000 −5.601E−15 4.000 −5.603E−15 4.000 −5.609E−15 4.002
1.75 100 6.695E−03 – 5.745E−06 – −1.102E−06 – −1.164E−06 – −1.373E−06 –
200 1.674E−03 2.000 3.587E−07 4.001 −6.914E−08 3.994 −7.099E−08 4.035 −7.753E−08 4.146
400 4.185E−04 2.000 2.242E−08 4.000 −4.326E−09 3.999 −4.382E−09 4.018 −4.586E−09 4.079
800 1.046E−04 2.000 1.401E−09 4.000 −2.704E−10 4.000 −2.722E−10 4.009 −2.786E−10 4.041
1600 2.616E−05 2.000 8.759E−11 4.000 −1.690E−11 4.000 −1.696E−11 4.005 −1.716E−11 4.021
3200 6.540E−06 2.000 5.474E−12 4.000 −1.056E−12 4.000 −1.058E−12 4.002 −1.064E−12 4.011
6400 1.635E−06 2.000 3.422E−13 4.000 −6.603E−14 4.000 −6.608E−14 4.001 −6.628E−14 4.005
12800 4.087E−07 2.000 2.138E−14 4.000 −4.127E−15 4.000 −4.128E−15 4.001 −4.134E−15 4.003
α Analytical values of R I [ 1 , 5 ] α y ( x ) | x = 2 calculated using formula (70)
0.25   6.9563532456344804165421264614628538
0.75   42.4546893190059613381179849166915634
1.25 −64.6142429211655969966421680694892918
1.75 −32.5941704287460581059377804482796869

In Fig. 1, the plot of function (69) (which corresponds to α = 0) and the plots of the Riesz fractional integrals R I [ 1 , 5 ] α y ( x ) (71) for α = 0.25 , 0.5 , 0.75 , 1.25 , 1.5 , 1.75 , 2 are shown. While, the numerical errors err N obtained on the interval [ 1 , 5 ] for N = 100 , 200 , 400 , 800 , 1600 , 3200 , 6400 , 12800 and α = 0.75, for five considered methods are illustrated in the plots in Fig. 2. It can be seen for this considered problem that the error values can be both positive and negative, but errors always decrease as N increases.

Plots of R I [ 1 , 5 ] α y ( x ) for function (68) and α ∈ { 0 , 0.25 , 0.5 , 0.75 , 1.25 , 1.5 , 1.75 , 2 }.

The numerical errors err N = R I [ 1 , 5 ] α y ( x ) − Ψ N for different kinds of the spline interpolation methods, α = 0.75 and N = 100 , 200 , 400 , 800 , 1600 , 3200 , 6400 , 12800.

The last example concerns the comparison of the calculation results of the developed methods with an external approach. Here, we take into account the approximation method developed by Dimitrov (2021). His asymptotic expansion formula for the trapezoidal approximation of the left-sided fractional integral is as follows (the symbols in this formula have been adapted to the symbols used in this work) I a + α y ( x ) ≃ 1 Γ ( α ) ( ( Δ x ) α ∑ k = 1 N − 1 y ( x − k Δ x ) k 1 − α + y ( a ) 2 ( x − a ) 1 − α Δ x − ( α − 1 ) ( x − a ) α − 2 y ( a ) − ( x − a ) α − 1 y ′ ( a ) 12 ( Δ x ) 2 − ζ ( 1 − α ) y ( x ) ( Δ x ) α + ζ ( − α ) y ′ ( x ) ( Δ x ) α + 1 − ζ ( − 1 − α ) y ″ ( x ) 2 ( Δ x ) α + 2 + ζ ( − 2 − α ) y ‴ ( x ) 6 ( Δ x ) α + 3 ) + O ( ( Δ x ) 4 ) , where ζ ( · ) is the Riemann zeta function.

Based on one of the examples presented in Dimitrov (2021), we take function y ( x ) = exp ( λ x ), for which the analytical form of the left-sided fractional integral of order α > 0 is the following I a + α exp ( λ x ) = exp ( λ a ) ( x − a ) α E 1 , 1 + α ( λ ( x − a ) ) , where E α , β ( · ) is two-parameter Mittag-Leffler function (Kilbas et al., 2006) and λ is a constant.

Table 4 presents the errors err N and the EOC’s values calculated for three variants of the cubic splines, whereby the results for the Dimitrov’s method have been taken from his work (see: Dimitrov, 2021). The example calculations have been performed for λ = 1, α = 0.5, a = 0, b = 2 and N = 40 , 80 , 160 , 320 , 640. The analytical value of I 0 + 0.5 exp ( x ) | x = 2 is also given in the table.

Results related to Example 4.

N Δ x Dimitrov’s method Cubic spline
Variant 1 Variant 2 Variant 3
err N EOC err N EOC err N EOC err N EOC
40 0.05 no data – 4.87E−08 – 7.98E−08 – 1.66E−07 –
80 0.025 1.06E−09 4.04725 3.46E−09 3.81833 4.66E−09 4.09663 8.45E−09 4.29274
160 0.0125 6.48E−11 4.03414 2.27E−10 3.92785 2.76E−10 4.07693 4.44E−10 4.25181
320 0.00625 3.98E−12 4.02694 1.45E−11 3.96779 1.66E−11 4.05782 2.40E−11 4.20937
640 0.003125 2.34E−13 4.08360 9.17E−13 3.98378 1.01E−12 4.04244 1.33E−12 4.16883
Analytical value of I 0 + 0.5 exp ( x ) | x = 2 = 7.052852096484309014376129 (calculated using (73))

As can be seen in Table 4, the results for all variants of the cubic splines and for the Dimitrov’s method have the identical EOC = 4, but for the same integration step sizes, the obtained numerical errors are slightly larger in the case of the cubic spline methods. As an explanation for the smaller error values in the Dimitrov’s article, one should delve into his method. Equation (72) contains the derivatives (up to the third order) of the integrand function. Knowledge of the analytical values of these derivatives at points a and x significantly improves the quality of numerical evaluation of the fractional integral in his developed method. His method can be considered a composition of the analytical and numerical approaches. Our computational schemes based on the cubic splines are purely numerical, and also the values of derivatives of the integrand function at the end point conditions are determined only numerically. Unfortunately, in the case of a very complicated form of the integrand function, it is often difficult to find the values of the first, second and third derivatives in an analytical way, and hence it can be problematic to use it for the Dimitrov’s method.

 Conclusions

We have derived the numerical integration formulas for evaluating the fractional order integrals. These formulas are based on the interpolation of the integrand function by the linear, quadratic and cubic splines. The integration method that uses the linear spline (i.e. the piecewise linear function) is treated as a complement and is used to compare the obtained results with the methods that use the quadratic and cubic splines.

Analysing the results presented in all tables and graphs (as well as other performed results, but omitted and not presented in this work), it can be noticed that the calculated numerical values of the fractional order integrals have good agreement with the exact analytical solution (if this solution is known, of course), and the numerical errors err N tend to 0 as the grid size N increases for all derived integration methods. We can see that as N increases, the values of the EOC are stabilized and tend to the specified values, and for different kinds of splines used we obtained: EOC = 2 for the linear spline, EOC = min { 3 + α , 4 } for the quadratic spline, EOC = 4 for all variants of the clamped cubic splines. Here, one can observe the consistency of the obtained EOC values for the particular methods with the orders of error (defined by the O notation) which were estimated analytically. Moreover, analysing the numerical results (by comparing the order of numerical errors, as well as the EOC) for three variants of the cubic splines, we can notice that these numerical methods have a similar order of errors, however, it can be seen that variant 1 generally generates slightly smaller errors. If we compare the order of numerical errors for methods using the quadratic and cubic splines (where the E O C = 4 is the same), we can conclude that the quadratic spline gives worse results than the cubic splines.

From a computational point of view, in the case of the methods that use the clamped cubic splines (all variants), the linear system of equations needs to be solved, which can be computationally time consuming, and thus it has certain disadvantages of the methods as a result. But here, in the case of derived methods, we need only to solve a tri-diagonal system of equations, which can be solved easily and faster, e.g. by using the Thomas algorithm of complexity O(N). This problem does not occur in the case of the linear and quadratic splines where the coefficients of the interpolation spline are determined locally, and these algorithms can be easily parallelized.

Summing up, it can be said that the obtained numerical results for the derived integration methods that use the cubic or quadratic splines give higher values of the EOC than for the method that uses the linear splines. This confirms the efficiency and the applicability of the derived numerical algorithms for the fractional integration of function.

 References Almeida, R., Pooseh, S., Torres, D.F.M. (2015). Computational Methods in the Fractional Calculus of Variations. Imperial College Press, London. Baleanu, D., Diethelm, K., Scalas, E., Trujillo, J.J. (2012). Fractional Calculus: Models and Numerical Methods. World Scientific, Singapore. Blaszczyk, T., Siedlecki, J. (2014). An approximation of the fractional integrals using quadratic interpolation. Journal of Applied Mathematics and Computational Mechanics, 13(4), 13–18. https://doi.org/10.17512/jamcm.2014.4.02. Blaszczyk, T., Siedlecki, J., Ciesielski, M. (2018). Numerical algorithms for approximation of fractional integral operators based on quadratic interpolation. Mathematical Methods in the Applied Sciences, 41(9), 3345–3355. https://doi.org/doi:10.1002/mma.4828. Budak, H., Hezenci, F., Kara, H., Sarikaya, M.Z. (2023). Bounds for the error in approximating a fractional integral by Simpson’s rule. Mathematics, 11(10), 16. https://doi.org/10.3390/math11102282. Burden, R.L., Faires, J.D., Burden, A.M. (2016). Numerical Analysis, 10th edition. Cengage Learning, Boston. Cai, M., Li, C. (2020). Numerical approaches to fractional integrals and derivatives: a review. Mathematics, 8(1), 43. https://doi.org/10.1002/mma.4828. de Oliveira, E.C., Machado, J.T. (2014). A review of definitions for fractional derivatives and integral. Mathematical Problems in Engineering, 2014, 238459. https://doi.org/10.1155/2014/238459. Dimitrov, Y. (2021). Approximations of the fractional integral and numerical solutions of fractional integral equations. Communications on Applied Mathematics and Computation, 3, 545–569. https://doi.org/10.1007/s42967-021-00132-7. Engeln-Müllges, G., Uhlig, F. (1996). Numerical Algorithms with C. Springer, Berlin. Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation, 51(184), 699–706. Kilbas, A.A., Srivastava, H.M., Trujillo, J.J. (2006). Theory and Applications of Fractional Differential Equations. Elsevier, Amsterdam. Li, C., Zeng, F. (2015). Numerical Methods for Fractional Calculus. Chapman and Hall/CRC, New York. Malinowska, A.B., Odzijewicz, T., Torres, D.F.M. (2015). Advanced Methods in the Fractional Calculus of Variations. Springer International Publishing, London. Odibat, Z. (2006). Approximations of fractional integrals and Caputo fractional derivatives. Applied Mathematics and Computation, 178(2), 527–533. https://doi.org/doi:10.1016/j.amc.2005.11.072. Oldham, K.B., Spanier, J. (1974). The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order. Academic Press, San Diego. Podlubny, I. (1999). Fractional Differential Equations. Academic Press, San Diego. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P. (2007). Numerical Recipes: The Art of Scientific Computing, 3rd ed. Cambridge University Press, New York. Rubin, B. (1996). Fractional Integrals and Potentials. Taylor and Francis, London.