Numerical Approximations of the Riemann–Liouville and Riesz Fractional Integrals

Ciesielski, Mariusz; Grodzki, Grzegorz

doi:10.15388/23-INFOR540

Numerical Approximations of the Riemann–Liouville and Riesz Fractional Integrals

Article type: Research Article

Authors: Ciesielski, Mariusz^{; ∗} | Grodzki, Grzegorz

Affiliations: Department of Computer Science, Czestochowa University of Technology, Dabrowskiego 73, 42-200 Czestochowa, Poland

Correspondence: [∗] Corresponding author.

Keywords: fractional calculus, numerical integration, numerical algorithms, spline interpolation

DOI: 10.15388/23-INFOR540

Journal: Informatica, vol. 35, no. 1, pp. 21-46, 2024

Received May 2023

Accepted November 2023

Published: 2024

Get PDF

Abstract

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.

1Introduction

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)

(1)

Ia+αy(x)=1Γ(α)∫axy(ξ)(x−ξ)1−αdξ,forx>a,

where Γ denotes the Gamma function. Whereas, the right-sided fractional integral is of the form

(2)

Ib−αy(x)=1Γ(α)∫xby(ξ)(ξ−x)1−αdξ,forx<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)

(3)

RI[a,b]αy(x)=12Γ(α)cos(απ/2)∫aby(ξ)|ξ−x|1−αdξ=12cos(απ/2)(Ia+αy(x)+Ib−αy(x)),fora<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.

2Numerical 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

(4)

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.

2.1Spline 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 [xi,xi+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=x0<x1<⋯<xN=b, xi=a+iΔx. At the given set of points, the values of function y(x) are determined as yi=y(xi), for i=0,1,…,N. The piecewise function s(x) can be expressed as

(5)

s(x)=s0(x),ifx∈[x0,x1],s1(x),ifx∈[x1,x2],…sN−1(x),ifx∈[xN−1,xN],

where si(x) are polynomials of degree p in each sub-interval [xi,xi+1], for i=0,1,…,N−1, and here we write these polynomials in the form

(6)

si(x)=∑k=0pck,i(x−xi)k,

where ck,i, for k=0,1,…,p, are the coefficients of polynomial si(x) in i-th sub-interval [xi,xi+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 (x0,y0),(x1,y1),…,(xN,yN) is the relationship s(xi)=yi, for i=0,1,…,N.

2.1.1Linear 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 (xi,yi) and (xi+1,yi+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

(7)

si(x)=c0,i+c1,i(x−xi).

The linear Lagrange interpolating polynomial in the interval x∈[xi,xi+1], for i=0,1,…,N−1, that passes through the points (xi,yi) and (xi+1,yi+1) results from the relation

(8)

y(x)=x−xi+1xi−xi+1y(xi)+x−xixi+1−xiy(xi+1)+y″(x¯i)2!(x−xi)(x−xi+1)=yi+yi+1−yiΔx(x−xi)+y″(x¯i)2(x−xi)(x−xi+1)=si(x)+Err1i(x),

and hence, the polynomial coefficients in Eq. (7) are as follows

(9)

c0,i=yi,c1,i=yi+1−yiΔx.

Moreover, the error term in i-th interval is given as

(10)

Err1i(x)=y″(x¯i)2(x−xi)(x−xi+1),

where x¯i∈[xi,xi+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).

2.1.2Quadratic Spline Interpolation

In this kind of approximation, the second degree polynomial (p=2) in each i-th sub-interval is used

(11)

si(x)=c0,i+c1,i(x−xi)+c2,i(x−xi)2,

where x∈[xi,xi+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 [xi−1,xi+1], for i=1,3,…,N−1, is expanded in the Taylor polynomial of third order about the central point xi as

(12)

y(x)=y(xi)+y′(xi)(x−xi)+y″(xi)2!(x−xi)2+y‴(xi)3!(x−xi)3+y(4)(x¯1i)4!(x−xi)4,

where x¯1i∈[xi−1,xi+1] is a certain point in the interval. The values of y′(xi) and y″(xi) in Eq. (12) are replaced by well-known relationships

(13)

y′(xi)=12Δx(y(xi+1)−y(xi−1))−(Δx)26y‴(x¯2i),y″(xi)=1(Δx)2(y(xi+1)−2y(xi)+y(xi−1))−(Δx)212y(4)(x¯3i),

where x¯2i∈[xi−1,xi+1] and x¯3i∈[xi−1,xi+1], and hence it was obtained

(14)

y(x)=y(xi)+(y(xi+1)−y(xi−1)2Δx−(Δx)26y‴(x¯2i))(x−xi)+12!(y(xi+1)−2y(xi)+y(xi−1)(Δx)2−(Δx)212y(4)(x¯3i))(x−xi)2+y‴(xi)3!(x−xi)3+y(4)(x¯1i)4!(x−xi)4,

or written in the form

(15)

y(x)=si(x)+Err2i(x),forx∈[xi−1,xi+1],

where

(16)

si(x)=yi+yi+1−yi−12Δx(x−xi)+yi+1−2yi+yi−12(Δx)2(x−xi)2

and

(17)

Err2i(x)=16y‴(xi)(x−xi)3−(Δx)26y‴(x¯2i)(x−xi)+124y(4)(x¯1i)(x−xi)4−(Δx)224y(4)(x¯3i)(x−xi)2≅16((x−xi)3−(Δx)2(x−xi))y‴(x¯i)+124((x−xi)4−(Δx)2(x−xi)2)y(4)(x¯i),

wherein x¯i∈[xi−1,xi+1] is taking into account as common values of xi, x¯1i, x¯2i and x¯3i.

In order to maintain a uniform form of the spline notation (as in Eq. (11)), two cases are considered

(18)

ck,i=c¯¯k,i,forieven,c¯k,i,foriodd,fork=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

(19)

c¯0,i=yi,c¯1,i=yi+1−yi−12Δx,c¯2,i=yi−1−2yi+yi+12(Δx)2,

and then the relationships si(xi−1)=yi−1, si(xi)=yi, si(xi+1)=yi+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

(20)

c¯¯0,i=yi,c¯¯1,i=−3yi+4yi+1−yi+22Δx,c¯¯2,i=yi−2yi+1+yi+22(Δx)2,

and the properties si(xi)=yi, si(xi+1)=yi+1, si(xi+2)=yi+2 hold.

As can be seen in the first case, the quadratic segment is defined in the interval [xi,xi+2] (for i even), while in the second case it is defined in the interval [xi−1,xi+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.

2.1.3Cubic 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 (xi,yi), 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.

(21)

si(x)=c0,i+c1,i(x−xi)+c2,i(x−xi)2+c3,i(x−xi)3,

where x∈[xi,xi+1], for i=0,1,…,N−1. These four coefficients c0,i, c1,i, c2,i and c3,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 4N 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: si(xi)=yi, si(xi+1)=yi+1, for i=0,1,…,N−1, (two conditions for each sub-interval give 2N 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): si−1′(xi)=si′(xi), si−1″(xi)=si″(xi), for i=1,2,…,N−1, (two conditions for each internal node give further 2N−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 x0 and xN (colloquially called the end point conditions). Three variants of the cubic spline constructions are considered:
- • variant 1: s0′(x0)=Y0′, sN−1′(xN)=YN′,
- • variant 2: s0″(x0)=Y0″, sN−1″(xN)=YN″,
- • variant 3: s0‴(x0)=Y0‴, sN−1‴(xN)=YN‴,
where Y0′, Y0″, Y0‴, YN′, YN″ and YN‴ correspond to the values of the first, second and third order derivatives of the function y(x) at the nodes x0 and xN, respectively. In principle, any combination of derivative orders at both boundary nodes can be considered (e.g. s0‴(x0) and sN−1′(xN)).

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

(22)

c0,i=yi,c1,i=yi+1−yiΔx−Δx3(c2,i+1+2c2,i),c3,i=13Δx(c2,i+1−c2,i),

and hence, the cubic spline in the i-th sub-interval (Eq. (21)) can be written in the form

(23)

si(x)=yi+(yi+1−yiΔx−Δx3(c2,i+1+2c2,i))(x−xi)+c2,i(x−xi)2+13Δx(c2,i+1−c2,i)(x−xi)3.

It can now be seen that only N+1 unknown coefficients c2,i need to be determined (the last one with the index N is used to evaluate c1,N−1 and c3,N−1). The first, second and third derivatives of the spline si(x) are as follows

(24)

si′(x)=(yi+1−yiΔx−Δx3(c2,i+1+2c2,i))+2c2,i(x−xi)si′(x)=+1Δx(c2,i+1−c2,i)(x−xi)2,si″(x)=2c2,i+2Δx(c2,i+1−c2,i)(x−xi),si‴(x)=2Δx(c2,i+1−c2,i).

It is easy to check that the conditions: si(xi)=yi, si(xi+1)=yi+1, for i=0,1,…,N−1, and si−1″(xi)=si″(xi), for i=1,2,…,N−1, are always present. While, from the conditions si−1′(xi)=si′(xi), for i=1,2,…,N−1, as a result, the following dependencies are obtained

(25)

c2,i+1+4c2,i+c2,i−1=3yi+1−2yi+yi−1(Δx)2,fori=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 x0 and xN are inserted into the appropriate spline derivatives, then we get

(26)

s0′(x0)=y1−y0Δx−Δx3(c2,1+2c2,0),s0″(x0)=2c2,0,s0‴(x0)=2Δx(c2,1−c2,0),sN−1′(xN)=yN−yN−1Δx+Δx3(c2,N−1+2c2,N),sN−1″(xN)=2c2,N,sN−1‴(xN)=2Δx(c2,N−c2,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: 2c2,0+c2,1=3(y1−y0(Δx)2−Y′0Δx),
c2,N−1+2c2,N=3(Y′NΔx−yN−yN−1(Δx)2),
• variant 2: c2,0=12Y0″, c2,N=12YN″,
• variant 3: c2,0−c2,1=−Δx2Y0‴, −c2,N−1+c2,N=Δx2YN‴.

The above derived relationships create the linear system of equations that can be written in matrix form

(27)

γ0,0γ0,100⋯000141000001410000014000⋮⋱⋮000041000001410000⋯0γN,N−1γN,N·c2,0c2,1c2,2c2,3⋮c2,N−2c2,N−1c2,N=d0d1d2d3⋮dN−2dN−1dN,

where

(28)

di=3yi+1−2yi+yi−1(Δx)2,fori=1,2,…,N−1

and

• for variant 1: γ0,0=2, γ0,1=1, d0=3(y1−y0(Δx)2−Y0′Δx),
γN,N−1=1, γN,N=2, dN=3(YN′Δx−yN−yN−1(Δx)2),
• for variant 2: γ0,0=1, γ0,1=0, d0=12Y0″,
γN,N−1=0, γN,N=1, dN=12YN″,
• for variant 3: γ0,0=1, γ0,1=−1, d0=−Δx2Y0‴,
γN,N−1=−1, γN,N=1, dN=Δx2YN‴.

The values of Y0′, Y0″, Y0‴, YN′, YN″ and YN‴ can be determined directly if the first, second and third order derivatives of function y(x) can be derived analytically, i.e. Y0′=y′(x0), YN′=y′(xN) 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)

(29)

Y0′≅1Δx(−2512y0+4y1−3y2+43y3−14y4),YN′≅1Δx(2512yN−4yN−1+3yN−2−43yN−3+14yN−4),

(30)

Y0″≅1(Δx)2(154y0−776y1+1076y2−13y3+6112y4−56y5),YN″≅1(Δx)2(154yN−776yN−1+1076yN−2−13yN−3+6112yN−4−56yN−5),

(31)

Y0‴≅1(Δx)3(−498y0+29y1−4618y2+62y3−3078y4+13y5−158y6),YN‴≅1(Δx)3(498yN−29yN−1+4618yN−2−62yN−3+3078yN−4YN‴≅−13yN−5+158yN−6).

By solving the system of equations (27), we obtain the values of polynomial coefficients c2,i, for i=0,1,…,N, in Eq. (21), while the values of remaining coefficients c0,i, c1,i, c3,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 ck,i in Eq. (21), the value of the interpolation function s(x) at an arbitrary point x depends on all data points (xi,yi), for i=0,1,…,N. Here, the variation or perturbation of one arbitrary data point (xi,yi) 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 x0 and xN 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)∈C4[a,b] and s(x) is the clamped cubic spline that interpolates function y(x) with respect to the nodes xi, for i=0,1,…,N, then the following relationship occurs for all x∈[a,b]

(32)

Err3=|y(x)−s(x)|⩽5384maxa⩽x⩽b|y(4)(x)|(Δx)4.

2.2Numerical 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∈{xM}, xM=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

(33)

Ia+αy(x)x=g≅Ia+αs(x)x=g=xM=1Γ(α)∫ags(ξ)(g−ξ)1−αdξ=∑i=0M−11Γ(α)∫xixi+1si(ξ)(xM−ξ)1−αdξ,forM=1,…,N,

(34)

Ib−αy(x)x=g≅Ib−αs(x)x=g=xM=1Γ(α)∫gbs(ξ)(ξ−g)1−αdξ=∑i=MN−11Γ(α)∫xixi+1si(ξ)(ξ−xM)1−αdξ,forM=0,…,N−1.

The fractional integrals containing the spline function s(x) are the sum of the local integrals.

If we use the notation: xi=a+iΔx, then the general form of (6) is of the form

(35)

si(x)=∑k=0pck,i(x−(a+iΔx))k

and if we put Eq. (35) into Eq. (33), then we obtain

(36)

Ia+αs(x)x=xM=∑i=0M−11Γ(α)∫a+iΔxa+(i+1)Δx∑k=0pck,i(ξ−(a+iΔx))k(a+MΔx−ξ)1−αdξ=∑i=0M−1∑k=0pck,i(1Γ(α)∫a+iΔxa+(i+1)Δx(ξ−(a+iΔx))k(a+MΔx−ξ)1−αdξ)

or in the form

(37)

Ia+αs(x)x=xM=∑i=0M−1∑k=0pck,iJi,ML,k,

where the particular integrals Ji,ML,k, for k=0,1,…,3, with regard to the data point xM in i-th sub-interval (i<M), one can find in an analytical way, and they are as follows:

(38)

Ji,ML,k=1Γ(α)∫a+iΔxa+(i+1)Δx(ξ−(a+iΔx))k(a+MΔx−ξ)1−αdξ=ξ=a+(u+i)Δx(Δx)α+kΓ(α)∫01uk(M−i−u)1−αduJi,ML,k=(Δx)α+kΓ(α+k+1)(k!(M−i)α+k−di,ML,k(M−i−1)α),

where

(39)

di,ML,k=1,ifk=0,(M−i−1)+(α+1),ifk=1,2(M−i−1)2+2(α+2)(M−i−1)+(α+1)(α+2),ifk=2,6(M−i−1)3+6(α+3)(M−i−1)2+3(α+2)(α+3)(M−i−1)+(α+1)(α+2)(α+3),ifk=3.

In a similar manner, we proceed with the second integral (34) and we have

(40)

Ib−αs(x)|x=xM=∑i=MN−11Γ(α)∫a+iΔxa+(i+1)Δx∑k=0pck,i(ξ−(a+iΔx))k(ξ−(a+MΔx))1−αdξ=∑i=MN−1∑k=0pck,i(1Γ(α)∫a+iΔxa+(i+1)Δx(ξ−(a+iΔx))k(ξ−(a+MΔx))1−αdξ)

(41)

Ib−αs(x)|x=xM=∑i=MN−1∑k=0pck,iJi,MR,k,

and the integrals Ji,MR,k, for k=0,1,…,3, also with respect to data point xM in i-th sub-interval (i⩾M), take the following analytical forms

(42)

Ji,MR,k=1Γ(α)∫a+iΔxa+(i+1)Δx(ξ−(a+iΔx))k(ξ−(a+MΔx))1−αdξ=ξ=a+(u+i)Δx(Δx)k+αΓ(α)∫01uk(i−M+u)1−αduJi,MR,k=(Δx)α+kΓ(α+k+1)((−1)k+1k!(i−M)α+k+di,MR,k(i−M+1)α),

where

(43)

di,MR,k=1,ifk=0,−(i−M+1)+(α+1),ifk=1,2(i−M+1)2−2(α+2)(i−M+1)+(α+2)(α+1),ifk=2,−6(i−M+1)3+6(α+3)(i−M+1)2−3(α+2)(α+3)(i−M+1)+(α+3)(α+2)(α+1),ifk=3.

Using the above numerical integration schemes, one can derive the numerical algorithm for the Riesz fractional integral (3) in the following way

(44)

for M=1,…,N−1.

2.3Error 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.

2.3.1The Case of the Linear Spline

Taking the particular error terms Err1i (10), for i=0,1,…,N−1, we integrate them over each interval separately using the properties of integral calculus

(45)

Err=∑i=0N−11Γ(α)∫xixi+1Err1i(ξ)(xN−ξ)1−αdξ=∑i=0N−11Γ(α)∫xixi+112(ξ−xi)(ξ−xi+1)y″(x¯i)1(xN−ξ)1−αdξ⩽12maxi=0,1,…,N−1|y″(x¯i)|∑i=0N−11Γ(α)∫xixi+1(ξ−xi)(ξ−xi+1)(xN−ξ)1−αdξ=12M2θ2N(α,Δx),

where

(46)

M2=maxi=0,1,…,N−1|y″(x¯i)|,θ2N(α,Δx)=∑i=0N−11Γ(α)∫xixi+1(ξ−xi)(ξ−xi+1)(xN−ξ)1−αdξ=ξ=xi+uΔx∑i=0N−1(Δx)2+αΓ(α)∫01u(u−1)(N−i−u)1−αdu.

After analytical integration and simplifications we obtain

(47)

θ2N(α,Δx)=(Δx)2+α(−2N3+αΓ(3+α)+1Γ(2+α)(2∑i=1Ni1+α−N1+α)).

We approximate the finite sum ∑i=1Ni1+α by using the Euler-Maclaurin formula of the second order. Then we get

(48)

θ2N(α,Δx)=(Δx)2+α(C2(α)+16Γ(1+α)Nα),

where

(49)

C2(α)=−2Γ(3+α)+1Γ(2+α)−16Γ(1+α).

It can be noticed that for α=1, we have C2(1)=0, θ2N(1,Δx)=(Δx)3N/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

(50)

Err⩽12M2((Δx)2+α(C2(α)+16Γ(1+α)(b−aΔx)α))=M2(C2(α)2(Δx)2+α+(b−a)α12Γ(1+α)(Δx)2)=M2(b−a)α12Γ(1+α)(Δx)2+O((Δx)2+α).

2.3.2The Case of the Quadratic Spline

Here, the error terms Err2i (17), for i=1,3,…,N−1, were determined on double adjacent intervals, hence their fractional integration concerns double intervals

(51)

Err=∑i=1,3,…,N−11Γ(α)∫xi−1xi+1Err2i(ξ)(xN−ξ)1−αdξ=∑i=1,3,…,N−11Γ(α)∫xi−1xi+116((ξ−xi)3−(Δx)2(ξ−xi))y‴(x¯i)+124((ξ−xi)4−(Δx)2(ξ−xi)2)y(4)(x¯i)1(xN−ξ)1−αdξ⩽16maxi=1,3,…,N−1|y‴(x¯i)|×∑i=1,3,…,N−11Γ(α)∫xi−1xi+1(ξ−xi)3−(Δx)2(ξ−xi)(xN−ξ)1−αdξ+124maxi=1,3,…,N−1|y(4)(x¯i)|×∑i=1,3,…,N−11Γ(α)∫xi−1xi+1(ξ−xi)4−(Δx)2(ξ−xi)2(xN−ξ)1−αdξ=16M3θ3N(α,Δx)+124M4θ4N(α,Δx),

where

(52)

M3=maxi=1,3,…,N−1|y‴(x¯i)|,M4=maxi=1,3,…,N−1|y(4)(x¯i)|,

(53)

θ3N(α,Δx)=∑i=1,3,…,N−11Γ(α)∫xi−1xi+1(ξ−xi)3−(Δx)2(ξ−xi)(xN−ξ)1−αdξ=ξ=xi+uΔx∑i=1,3,…,N−1(Δx)3+αΓ(α)∫−11u3−u(N−i−u)1−αdu,θ4N(α,Δx)=∑i=1,3,…,N−11Γ(α)∫xi−1xi+1(ξ−xi)4−(Δx)2(ξ−xi)2(xN−ξ)1−αdξ=ξ=xi+uΔx∑i=1,3,…,N−1(Δx)4+αΓ(α)∫−11u4−u2(N−i−u)1−αdu.

Further transformations (including using analytic integration and changing the indexing order in the series) lead to

(54)

θ3N(α,Δx)=(Δx)3+α6N3+αΓ(4+α)−6Γ(3+α)(23+α∑i=1N/2i2+α−N2+a)+2N1+αΓ(2+α),θ4N(α,Δx)=(Δx)4+α(24N4+αΓ(5+α)−24Γ(4+α)(24+α∑i=1N/2i3+α−N3+a)+10N2+αΓ(3+α)−2Γ(2+α)(22+α∑i=1N/2i1+α−N1+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

(55)

θ3N(α,Δx)=(Δx)3+α(C3(α)+215Γ(α)Nα−1),θ4N(α,Δx)≅(Δx)4+α(C4(α)−215Γ(α+1)Nα+2(α−1)45Γ(α)Nα−2),

where

(56)

C3(α)=−3·24+αΓ(4+α)−3·23+αΓ(3+α)+22+αΓ(2+α)−2α15Γ(α),C4(α)=−6(2+α)25+αΓ(5+α)+5·23+αΓ(3+α)−22+αΓ(2+α)+21+α15Γ(1+α)−(α−1)2α90Γ(α).

It should be pointed out that for α=1, one obtains C3(1)=−2/15, C4(1)=0 and hence θ3N(1,Δx)=0, θ4N(1,Δx)=−2(Δx)5N/15 – which corresponds to the classic Simpson’s rule approximation for integral. By substituting θ3N(α,Δx) and θ4N(α,Δx) from Eq. (55) into Eq. (51) and also assuming that N=(b−a)/Δx one gets

(57)

Err⩽16M3(C3(α)(Δx)3+α+2(b−a)α−115Γ(α)(Δx)4)+124M4C4(α)(Δx)4+α−2(b−a)α15Γ(α+1)(Δx)4+2(α−1)(b−a)α−245Γ(α)(Δx)6.

Taking into account the terms containing Δx with the lowest powers, we have

(58)

Err=M3C3(α)6(Δx)3+α+O((Δx)4),forα<1,(M3(b−a)α−145Γ(α)−M4(b−a)α180Γ(α+1))(Δx)4+O((Δx)3+α),forα⩾1.

2.3.3The 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

(59)

Err=‖Ia+αy(x)x=b−Ia+αs(x)x=b‖=‖Ia+α(y(x)−s(x))x=b‖⩽Ia+α‖y(x)−s(x)‖x=b⩽‖y(x)−s(x)‖Ia+α1x=b=‖y(x)−s(x)‖(b−a)αΓ(1+α).

By inserting relationship (32) into Eq. (59) we obtain

(60)

Err⩽5M4384(b−a)αΓ(1+α)(Δx)4,

where M4=maxa⩽x⩽b|y(4)(x)| or for the discrete form M4=maxi=0,1,…,N−1|y(4)(x¯i)|, x¯i∈[xi,xi+1].

3Examples 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

(61)

errN=Ia+αy(x)x=g−ΨN,

where ΨN denotes the numerical/approximate value of Ia+αy(x)|x=g and simultaneously the exact value of Ia+αs(x)|x=g. Similarly, the error errN is defined for the right-sided fractional integral and the Riesz integral.

Next, the estimated EOC can be evaluated as

(62)

EOCN=log2|errN/2||errN|

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

(63)

EOCN=log2|ΨN−ΨN/2||Ψ2N−Ψ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) )

(64)

Ia+α(x−a)γ=Γ(γ+1)Γ(γ+α+1)(x−a)γ+α,forγ>−1,α>0,

(65)

Ib−α(b−x)γ=Γ(γ+1)Γ(γ+α+1)(b−x)γ+α,forγ>−1,α>0,

which are used to determine analytical solutions in the examples below.

Example 1.

The integrand function is the eighth-degree polynomial of the form

(66)

y(x)=x8−8x7+26x6−44x5+40x4−15x3−4x2+5x+1.

Numerical values of Ia+α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

(67)

I0+αy(x)x=b=Γ(9)Γ(9+α)b8+α−8Γ(8)Γ(8+α)b7+α+26Γ(7)Γ(7+α)b6+α−44Γ(6)Γ(6+α)b5+α+40Γ(5)Γ(5+α)b4+α−15Γ(4)Γ(4+α)b3+α−4Γ(3)Γ(3+α)b2+α+5Γ(2)Γ(2+α)b1+α+Γ(1)Γ(1+α)bα,

where Γ(k)=(k−1)!, for k=1,2,….

In Table 1, we report the numerical errors errN (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.

Table 1

Results related to Example 1.

α	N	Linear spline		Quadratic spline		Cubic spline
		Linear spline		Quadratic spline		Variant 1		Variant 2		Variant 3
		errN	EOC	errN	EOC	errN	EOC	errN	EOC	errN	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 I0+αy(x)\|x=2 calculated using formula (67)
0.4	3.6979129457596915301988815161146608
0.7	4.0856207593403175492511974048448624
1.4	4.3604818404289140653601695680338754
2.7	2.9484099812828967875285769194034989

Example 2.

In the second example, we take a more complicated integrand function of the form

(68)

y(x)=(xsin(3x2)+5xx+2)·exp(−(x−2)32−2x)+xx83xx2+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 Ia+α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=Ia+α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.

Table 2

Results related to Example 2.

α	N	Linear spline		Quadratic spline		Cubic spline
		Linear spline		Quadratic 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	–

Example 3.

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

(69)

y(x)=x5−13x4+59x3−108x2+67x+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 xM, 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)

(70)

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

(71)

RI[1,5]αy(x)|x=g=12cos(απ/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 errN (where errN=RI[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.

Table 3

Results related to Example 3.

α	N	Linear spline		Quadratic spline		Cubic spline
		Linear spline		Quadratic spline		Variant 1		Variant 2		Variant 3
		errN	EOC	errN	EOC	errN	EOC	errN	EOC	errN	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 RI[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 RI[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 errN 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.

Fig. 1

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

Fig. 2

The numerical errors errN=RI[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.

Example 4.

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)

(72)

Ia+αy(x)≃1Γ(α)((Δx)α∑k=1N−1y(x−kΔx)k1−α+y(a)2(x−a)1−αΔx−(α−1)(x−a)α−2y(a)−(x−a)α−1y′(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

(73)

Ia+αexp(λx)=exp(λa)(x−a)αE1,1+α(λ(x−a)),

where Eα,β(·) is two-parameter Mittag-Leffler function (Kilbas et al., 2006) and λ is a constant.

Table 4 presents the errors errN 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 I0+0.5exp(x)|x=2 is also given in the table.

Table 4

Results related to Example 4.

N	Δx	Dimitrov’s method		Cubic spline
		Dimitrov’s method		Variant 1		Variant 2		Variant 3
		errN	EOC	errN	EOC	errN	EOC	errN	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 I0+0.5exp(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.

4Conclusions

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 errN 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 EOC=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

1	Almeida, R., Pooseh, S., Torres, D.F.M. ((2015) ). Computational Methods in the Fractional Calculus of Variations. Imperial College Press, London.
2	Baleanu, D., Diethelm, K., Scalas, E., Trujillo, J.J. ((2012) ). Fractional Calculus: Models and Numerical Methods. World Scientific, Singapore.
3	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.
4	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.
5	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.
6	Burden, R.L., Faires, J.D., Burden, A.M. ((2016) ). Numerical Analysis, 10th edition. Cengage Learning, Boston.
7	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.
8	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.
9	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.
10	Engeln-Müllges, G., Uhlig, F. ((1996) ). Numerical Algorithms with C. Springer, Berlin.
11	Fornberg, B. ((1988) ). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation, 51: (184), 699–706.
12	Kilbas, A.A., Srivastava, H.M., Trujillo, J.J. ((2006) ). Theory and Applications of Fractional Differential Equations. Elsevier, Amsterdam.
13	Li, C., Zeng, F. ((2015) ). Numerical Methods for Fractional Calculus. Chapman and Hall/CRC, New York.
14	Malinowska, A.B., Odzijewicz, T., Torres, D.F.M. ((2015) ). Advanced Methods in the Fractional Calculus of Variations. Springer International Publishing, London.
15	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.
16	Oldham, K.B., Spanier, J. ((1974) ). The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order. Academic Press, San Diego.
17	Podlubny, I. ((1999) ). Fractional Differential Equations. Academic Press, San Diego.
18	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.
19	Rubin, B. ((1996) ). Fractional Integrals and Potentials. Taylor and Francis, London.