Approximation of the Riesz–Caputo Derivative by Cubic Splines

Abstract

Differential problems with the Riesz derivative in space are widely used to model anomalous diffusion. Although the Riesz–Caputo derivative is more suitable for modeling real phenomena, there are few examples in literature where numerical methods are used to solve such differential problems. In this paper, we propose to approximate the Riesz–Caputo derivative of a given function with a cubic spline. As far as we are aware, this is the first time that cubic splines have been used in the context of the Riesz–Caputo derivative. To show the effectiveness of the proposed numerical method, we present numerical tests in which we compare the analytical solution of several boundary differential problems which have the Riesz–Caputo derivative in space with the numerical solution we obtain by a spline collocation method. The numerical results show that the proposed method is efficient and accurate.

1. Introduction

In the last few decades, differential equations with fractional derivatives have emerged as a main tool for constructing models of real-world phenomena closer to experimental observations. Nowadays, models using fractional derivatives in time and/or in space are commonly used in various fields of sciences, such as biology, physics, mechanics, economics, control theory, just to cite a few (see, for instance, [1,2,3,4,5,6,7,8,9] and references therein). The reason lies in the fact that the fractional derivative is a nonlocal operator and it allows to approach the memory and hereditary properties of physical systems in a clear way.

In the field of geophysics, fractional differential equations are widely used for modeling anomalous diffusion in porous media, a phenomenon in which particles have been observed to spread at a rate which is incompatible with classical Brownian motion. In this context, the Riesz fractional derivative has been shown to be more attractive compared to the left-sided derivative since it is able to take into account contributions from both sides of the domain [2,10]. In fact, the Riesz derivative is a linear combination of the left- and right-sided derivatives, resulting in a fractional derivative centrally symmetric on finite domains. Although in literature the Riesz derivative defined through the left and right Riemann–Liouville derivatives is commonly used to model fractional diffusion, the Riesz–Caputo derivative seems more suitable for avoiding nonphysical issues (see [11]). For this reason, in this paper, we consider boundary value problems with the Riesz–Caputo derivative in space. In particular, the aim of the paper is to construct an approximation of the Riesz–Caputo derivative by cubic splines and use this approximation to solve fractional differential boundary value problems (FBVPs) using the spline collocation method introduced in [12,13]. The literature offers several numerical methods to approximate the Riesz derivative—finite difference methods [14,15,16,17], finite element methods [18,19,20], spectral methods [21,22,23], mesh-free methods [24,25]—while the use of spline approximations is still limited [26,27,28]. The main contribution of this paper is to give the explicit expression of the cubic spline approximation for the Riesz–Caputo derivative of a given function and to show its effectiveness in solving FBVPs using the spline collocation method. The choice of cubic splines, instead of linear or quadratic splines, was made because they yield smooth approximations while maintaining a low computational cost. To make the evaluation of the Riesz–Caputo derivative more efficient, we exploit the symmetry of the spline basis and avoid in this way the computation of the right Caputo derivative of the basis functions.

To the best of our knowledge, this is the first time that cubic splines have been used in the context of Riesz–Caputo derivatives.

The paper is organized as follows. In Section 2.1, we describe the fractional differential problem we are interested in. In Section 2.2, we recall some basic properties of the cubic B-spline and of the cubic spline basis on finite intervals. Section 2.3 is devoted to the approximation of the Riesz–Caputo derivative, while the collocation method is described in Section 2.4. To show the accuracy of the numerical approximation, in Section 3, we solve some FBVPs using the spline collocation method. Finally, in Section 4, we discuss conclusions and possible future developments.

2. Materials and Methods

2.1. Fractional Boundary Value Problems

We consider the one-dimensional boundary value problem

$$\left\{\begin{array}{c}{\displaystyle \frac{d^{\gamma }}{{d\left|x\right|}^{\gamma }}y\left(x\right)+f\left(x\right)y\left(x\right)=g\left(x\right),x\in (0,L),}\hfill \\ \\ {\rho }_{r0}y\left(0\right)+{\zeta }_{r0}y\left(L\right)=c_r,1\le r\le m,\hfill \end{array}\right.$$

(1)

where $0<\gamma <2$, $m=\lceil \gamma \rceil$, f and g are known functions and ${\rho }_{r0},{\zeta }_{r0},c_r$ are known parameters.

Here, $\frac{d^{\gamma }}{{d\left|x\right|}^{\gamma }}y\left(x\right)$ denotes the Riesz–Caputo fractional derivative in space defined as

$$\frac{d^{\gamma }}{{d\left|x\right|}^{\gamma }}y\left(x\right):=\frac{c\left(\gamma \right)}{2}(D_{0,x}^{\gamma }+{(-1)}^m{D}_{x,L}^{\gamma })y\left(x\right),$$

(2)

where the left and right Caputo derivatives are defined by

$$\begin{array}{c}\hfill D_{0,x}^{\gamma }y\left(x\right):=\frac{1}{\Gamma (m-\gamma )}{\int }_0^x\frac{y^{\left(m\right)}\left(\tau \right)}{{(x-\tau )}^{-m+\gamma +1}}d\tau ,\end{array}$$

(3)

$$\begin{array}{c}\hfill D_{x,L}^{\gamma }y\left(x\right):=\frac{{(-1)}^m}{\Gamma (m-\gamma )}{\int }_x^L\frac{y^{\left(m\right)}\left(\tau \right)}{{(\tau -x)}^{-m+\gamma +1}}d\tau ,\end{array}$$

(4)

where $\Gamma$ is the Euler gamma function.

In definition (2), $c\left(\gamma \right)$ is a given function. Common choices are

$$c\left(\gamma \right)=\frac{\pm 1}{cos\left(\frac{\pi }{2}\gamma \right)},c\left(\gamma \right)=1.$$

Here, we prefer the latter definition since in this case, the Riesz–Caputo derivative agrees with the ordinary derivative when $\gamma$ is an integer. Existence and unicity of the solution to (1) were discussed in [29,30]. For further details and properties of the Caputo and Riesz–Caputo derivatives, see [29,31,32,33].

2.2. The Cubic B-Spline Basis on Finite Intervals

Our aim is to approximate the solution to Equation (1) by a cubic spline, i.e., a piece-wise polynomial of degree 3. Cubic spline on the finite interval $[0,L]$ can be represented as a linear combination of a given cubic spline basis [34,35].

Without loss of generality, we assume L to be an integer greater than or equal to 4 and choose integer knots in the interval $[0,L]$:

$$\begin{array}{c}{x}_0=x_1=x_2=x_3=0,\hfill \\ \\ x_j=j-3,4\le j\le M-4,\hfill \\ \\ x_{M-3}=x_{M-2}=x_{M-1}=x_M=L.\hfill \end{array}$$

where $M=L+6$. We observe that the endpoints of the interval are repeated 4 times, meaning that they are knots of multiplicity 4. The internal knots $x_j$, $4\le j\le M-4$, are simple knots.

The cubic B-spline basis

$$\mathcal{N}=\{N_k\left(x\right),0\le k\le L+2\},x\in [0,L],$$

is formed by $L+3$ basis functions, 3 edge functions for each endpoints plus $L-3$ internal functions. The left- and right-edge functions have index $0\le k\le 2$ and $L\le k\le L+2$, respectively, while the internal functions correspond to the values $3\le k\le L-1$.

We recall that the basis functions $N_k$, $0\le k\le L+2$, are non-negative and compactly supported with

$$supp\left(N_k\right)=[max(0,k-3),min(k+1,L)],0\le k\le L+2.$$

Moreover, the basis $\mathcal{N}$ is centrally symmetric, i.e.,

$$N_k\left(x\right)=N_{L+2-k}(L-x),0\le k\le L+2,$$

(5)

and satisfies the endpoint conditions:

$$N_k\left(0\right)={\delta }_{0,k},N_k\left(L\right)={\delta }_{L+2,k},0\le k\le L+2,$$

(6)

where ${\delta }_{j,k}$ denotes the Kronecker delta.

For the convenience of the reader, the explicit expressions of the basis functions are given in Appendix A. For details on the construction of spline bases, see [35,36].

The function system $\mathcal{N}$ can be adapted to any sequence of equidistant knots by dilation:

$${\mathcal{N}}_h=\{N_k\left(\frac{x}{h}\right),0\le k\le M_h\},x\in [0,L],$$

(7)

where $M_h=\frac{L}{h}+2$. Figure 1 shows the basis for $L=1$ and $h=\frac{1}{8}$.

2.3. Approximation of the Riesz–Caputo Derivative

To approximate the Riesz–Caputo derivative of a give function y, we approximate the function by a cubic spline represented by a linear combination of the basis functions ${\mathcal{N}}_h$:

$$y\left(x\right)\approx y_h\left(x\right)=\sum _{k=0}^{M_h}{a}_{k,h}{N}_{k,h}\left(x\right),x\in [0,L].$$

(8)

We recall that cubic spline approximations belong to $C^2[0,L]$ and are exact on polynomials of degree less than or equal to 3 [34,35].

The Riesz–Caputo derivative (2) can be approximated as

$$\frac{d^{\gamma }}{{d\left|x\right|}^{\gamma }}y\left(x\right)\approx \frac{d^{\gamma }}{{d\left|x\right|}^{\gamma }}{y}_h\left(x\right)=\sum _{k=0}^{M_h}{a}_{k,h}\frac{d^{\gamma }}{{d\left|x\right|}^{\gamma }}{N}_{k,h}\left(x\right),x\in [0,L].$$

(9)

Thus, we need to evaluate the Riesz–Caputo derivative of the basis functions:

$$\frac{d^{\gamma }}{{d\left|x\right|}^{\gamma }}{N}_{k,h}\left(x\right):=\frac{1}{2}(D_{0,x}^{\gamma }+{(-1)}^m{D}_{x,L}^{\gamma })N_{k,h}\left(x\right),0\le k\le L+2.$$

To this end, it is sufficient to evaluate the derivative when $h=1$, since for any function f (see [12]),

$$D_{0,x}^{\gamma }f\left(\frac{x}{h}\right)=\frac{1}{h^{\gamma }}{D}_{0,h^{-1}x}^{\gamma }f\left(\frac{x}{h}\right),D_{x,L}^{\gamma }f\left(\frac{x}{h}\right)=\frac{1}{h^{\gamma }}{D}_{h^{-1}x,h^{-1}L}^{\gamma }f\left(\frac{x}{h}\right).$$

The explicit expression of the left Caputo derivative $D_{0,x}^{\gamma }{N}_k$ was given in [13,37] (see Appendix B). The right Caputo derivative can be evaluated by using the symmetry property (5).

Proposition 1.

Let $\Phi =\{{\varphi }_k\left(x\right),0\le k\le M\}$ be a function basis centrally symmetric in the interval $[0,L]$, i.e.,

$${\varphi }_k\left(x\right)={\varphi }_{M-k}(L-x),0\le k\le M,x\in [0,L].$$

(10)

Then, the left and right Caputo derivatives of the basis function Φ are related by the following symmetry properties:

$$D_{x,L}^{\gamma }{\varphi }_k\left(x\right)=D_{0,L-x}^{\gamma }{\varphi }_k(L-x)=D_{0,L-x}^{\gamma }{\varphi }_{M-k}(L-x),0\le k\le M.$$

(11)

Proof. 

For the first equality, using definition (4) and the change of variables $\overline{\tau }=L-\tau$ we obtain

$$\begin{array}{cc}\hfill D_{x,L}^{\gamma }{\varphi }_k\left(x\right)& =\frac{{(-1)}^m}{\Gamma (m-\gamma )}{\int }_x^L\frac{{\varphi }_k^{\left(m\right)}\left(\tau \right)}{{(\tau -x)}^{-m+\gamma +1}}d\tau \hfill \\ \hfill & =\frac{1}{\Gamma (m-\gamma )}{\int }_0^{L-x}\frac{{\varphi }_k^{\left(m\right)}(L-\overline{\tau })}{{(L-x-\overline{\tau })}^{-m+\gamma +1}}d\overline{\tau }\hfill \\ \hfill \\ \hfill & =D_{0,L-x}^{\gamma }{\varphi }_k(L-x).\hfill \end{array}$$

For the second equality, using relation (10) in definition (4) and the same change of variables as above, we obtain

$$\begin{array}{cc}\hfill D_{x,L}^{\gamma }{\varphi }_k\left(x\right)& =\frac{{(-1)}^m}{\Gamma (m-\gamma )}{\int }_x^L\frac{{\varphi }_k^{\left(m\right)}\left(\tau \right)}{{(\tau -x)}^{-m+\gamma +1}}d\tau \hfill \\ \hfill & =\frac{1}{\Gamma (m-\gamma )}{\int }_x^L\frac{{\varphi }_{M-k}^{\left(m\right)}(L-\tau )}{{(\tau -x)}^{-m+\gamma +1}}d\tau \hfill \\ \hfill & =\frac{1}{\Gamma (m-\gamma )}{\int }_0^{L-x}\frac{{\varphi }_{M-k}^{\left(m\right)}\left(\tilde{\tau }\right)}{{(L-x-\tilde{\tau })}^{-m+\gamma +1}}d\tilde{\tau }\hfill \\ \hfill \\ \hfill & =D_{0,L-x}^{\gamma }{\varphi }_{M-k}(L-x).\hfill \end{array}$$

□

Using the symmetry properties (5), Proposition (1) and Equations (A4)–(A12), we obtain the explicit expressions of the Riesz–Caputo derivatives $\frac{d^{\gamma }}{{d\left|x\right|}^{\gamma }}{N}_k$. We observe that the derivatives retain the same symmetry/antisymmetry properties as the ordinary derivatives.

The Riesz–Caputo derivative of the basis functions is plotted in Figure 2 and Figure 3 for different values of the order of the fractional derivative.

2.4. The Collocation Method

To show the effectiveness of the approximations we constructed in Section 2.3, we solve the FBVP (1) by the spline collocation method introduced in [12,13]. First, we approximate the solution to the differential problem by a spline function as in Equation (8). Then, by collocating the differential problem in the equidistant points

$$X_{\delta }=\{x_p=\delta p,0\le p\le M_{\delta }\},M_{\delta }={\delta }^{-1}L,$$

(12)

we obtain the linear system

$$\left\{\begin{array}{cc}{\displaystyle \sum _{k=0}^{M_h}{a}_{k,h}\frac{d^{\gamma }}{{d\left|x\right|}^{\gamma }}{N}_{k,h}\left(x_p\right)+f\left(x_p\right)\sum _{k=0}^{M_h}{a}_{k,h}{N}_{k,h}\left(x_p\right)=g\left(x_p\right),}\hfill & 1\le p\le M_{\delta }-1,\hfill \\ \\ {\displaystyle {\rho }_{r0}\sum _{k=0}^{M_h}{a}_{k,h}{N}_{k,h}\left(x_0\right)+{\zeta }_{r0}\sum _{k=0}^{M_h}{a}_{k,h}{N}_{k,h}\left(x_{M_{\delta }}\right)=c_r,}\hfill & 1\le r\le m,\hfill \end{array}\right.$$

(13)

where the coefficients $\left\{a_{k,h}\right\}$ are the $M_h+1$ unknowns. To avoid the linear system being underdetermined, we choose the step $\delta$ such that $M_{\delta }-1+m\ge M_h+1$.

3. Numerical Results

We will now test the collocation method described in Section 2.4 to solve Equation (1) where we will set without loss of generality $L=1$. For all the numerical tests, the known terms $f\left(x\right)$ and $g\left(x\right)$ are listed in Appendix C.

The linear system we obtain is an overdetermined linear system with $M_{\delta }-1+m$ equations and $M_h+1$ unknowns that can be solved using the least squares method.

In the following examples, we will denote by $\parallel e_h\parallel$ the infinity norm of the approximation error, evaluated as

$$\parallel e_h\parallel =\underset{0\le s\le 4M_{\delta }}{max}|y\left(x_s\right)-y_h\left(x_s\right)|,$$

(14)

where $y\left(x\right)$ is the exact solution to (1), $y_h\left(x\right)$ is the numerical solution (see (8)) and $x_s=\frac{\delta }{4}s,0\le s\le 4M_{\delta }$. The norm in (14) is evaluated through the local error $e\left(x\right)=y\left(x\right)-y_h\left(x\right)$ on a finer grid with step size $\delta /4$, being $\delta$ the step of the collocation points (see (12)). In the following numerical tests, we will set $\delta =h/2$.

3.1. Example 1

In the first example, we consider $0<\gamma <1$. We set the parameters ${\rho }_{10}={\zeta }_{10}=c_1=1$.

The analytical solution is $y\left(x\right)=x$ so that the collocation method is exact. In

Table 1. Example 1: The infinity norm of the error $\parallel e_h\parallel$ for h = 1/8.

$\mathit{\gamma }$	0.25	0.5	0.75
$\parallel e_h\parallel$	$4.36\times {10}^{-14}$	$1.38\times {10}^{-14}$	$1.24\times {10}^{-14}$

, we list the error we obtain for different values of $\gamma$ when $h=1/8$. As expected, the error is in the order of machine precision.

3.2. Example 2

To test the convergence of the collocation method, in the second example, we set the function $g\left(x\right)$ in order to obtain the analytical solution $y\left(x\right)=x^4$. The parameters are set to be ${\rho }_{10}={\zeta }_{10}=c_1=1$ for $0<\gamma <1$, and ${\rho }_{10}={\zeta }_{20}=c_2=1$, ${\rho }_{20}={\zeta }_{10}=c_1=0$, for $1<\gamma <2$.

The behavior of the error for different values of $\gamma$ and h is shown in

Table 2. Example 2: The infinity norm of the error $\parallel e_h\parallel$ for different values of $\gamma$ and h.

$\mathit{\gamma }$	0.3	0.6	0.9	1.2	1.5	1.8
$h=\frac{1}{8}$	$4.53\times {10}^{-5}$	$1.61\times {10}^{-5}$	$1.79\times {10}^{-5}$	$1.59\times {10}^{-4}$	$3.24\times {10}^{-4}$	$2.49\times {10}^{-3}$
$h=\frac{1}{16}$	$2.48\times {10}^{-6}$	$9.18\times {10}^{-7}$	$1.07\times {10}^{-6}$	$1.54\times {10}^{-5}$	$4.91\times {10}^{-5}$	$5.14\times {10}^{-4}$
$h=\frac{1}{32}$	$1.34\times {10}^{-7}$	$5.36\times {10}^{-8}$	$6.39\times {10}^{-8}$	$1.54\times {10}^{-6}$	$7.37\times {10}^{-6}$	$1.05\times {10}^{-4}$
$h=\frac{1}{64}$	$7.29\times {10}^{-9}$	$3.20\times {10}^{-9}$	$3.72\times {10}^{-9}$	$1.59\times {10}^{-7}$	$1.10\times {10}^{-6}$	$2.15\times {10}^{-5}$
$h=\frac{1}{128}$	$1.35\times {10}^{-9}$	$2.64\times {10}^{-10}$	$2.32\times {10}^{-10}$	$1.72\times {10}^{-8}$	$1.64\times {10}^{-7}$	$4.37\times {10}^{-6}$

. The numerical results show that the error increases when $\gamma$ increases. Moreover, the error decreases fast with h when $0<\gamma <1$, while it decreases slower when $1<\gamma <2$.

3.3. Example 3

Here, $0<\gamma <2$. We set the function $g\left(x\right)$ in order to obtain the periodic solution $y\left(x\right)=sin\left(\omega x\right)$ with $\omega =\pi ,5\pi$. For $0<\gamma <1$ in the boundary conditions, we set ${\rho }_{10}={\zeta }_{10}=1$ with $c_1=sin\left(\omega \right)$, while for $1<\gamma <2$, we set ${\rho }_{20}={\zeta }_{10}=c_1=0$, ${\rho }_{10}={\zeta }_{20}=1$, $c_2=sin\left(\omega \right)$.

The analytical and the numerical solutions are shown in Figure 4 together with the error $e\left(x\right)=y\left(x\right)-y_h\left(x\right)$.

In

Table 3. Example 3: The infinity norm of the error $\parallel e_h\parallel$ for different values of $\gamma$ and h when $\omega =\pi$.

$\mathit{\gamma }$	0.4	0.8	1.2	1.6
$h=\frac{1}{4}$	$2.11\times {10}^{-3}$	$1.31\times {10}^{-3}$	$2.82\times {10}^{-3}$	$3.87\times {10}^{-3}$
$h=\frac{1}{8}$	$6.66\times {10}^{-5}$	$4.71\times {10}^{-5}$	$2.27\times {10}^{-4}$	$5.81\times {10}^{-4}$
$h=\frac{1}{16}$	$2.46\times {10}^{-6}$	$2.10\times {10}^{-6}$	$2.63\times {10}^{-5}$	$1.08\times {10}^{-4}$
$h=\frac{1}{32}$	$1.30\times {10}^{-7}$	$1.31\times {10}^{-7}$	$3.50\times {10}^{-6}$	$2.08\times {10}^{-5}$
$h=\frac{1}{64}$	$7.59\times {10}^{-9}$	$8.67\times {10}^{-9}$	$4.88\times {10}^{-7}$	$3.99\times {10}^{-6}$
$h=\frac{1}{128}$	$4.80\times {10}^{-10}$	$6.05\times {10}^{-10}$	$6.90\times {10}^{-8}$	$7.61\times {10}^{-7}$

and

Table 4. Example 3: The infinity norm of the error $\parallel e_h\parallel$ for different values of $\gamma$ and h when $\omega =5\pi$.

$\mathit{\gamma }$	0.4	0.8	1.2	1.6
$h=\frac{1}{8}$	$1.14\times {10}^{-1}$	$6.78\times {10}^{-2}$	$2.70\times {10}^{-1}$	$1.46\times {10}^0$
$h=\frac{1}{16}$	$5.79\times {10}^{-3}$	$3.58\times {10}^{-3}$	$5.81\times {10}^{-3}$	$8.21\times {10}^{-2}$
$h=\frac{1}{32}$	$1.77\times {10}^{-4}$	$1.24\times {10}^{-4}$	$1.04\times {10}^{-3}$	$6.95\times {10}^{-3}$
$h=\frac{1}{64}$	$6.13\times {10}^{-6}$	$6.50\times {10}^{-6}$	$1.59\times {10}^{-4}$	$8.58\times {10}^{-4}$
$h=\frac{1}{128}$	$3.33\times {10}^{-7}$	$4.43\times {10}^{-7}$	$2.31\times {10}^{-5}$	$1.25\times {10}^{-4}$

, the infinity norm of the error for different values of $\gamma$ and h is listed. As expected, the most oscillatory solution, which is the hardest function to approximate, requires a small step h to obtain a small error.

We also show how the left and right contributions of the Riesz–Caputo derivative behave for this periodic function in the case $0<\gamma <1$ (Figure 5) and $1<\gamma <2$ (Figure 6). The plots show that when $0<\gamma <1$, the limit of the Riesz–Caputo derivative for $\gamma \to 1$ approaches the ordinary first derivative (in this case $y^{\prime }\left(x\right)=cos\left(\omega x\right)$); it is the same case for the left and right contributions, the latter with a minus sign due to the factor ${(-1)}^m$. Similarly, when $1<\gamma <2$, the Riesz–Caputo derivative approaches $y^{″}\left(x\right)=-sin\left(\omega x\right)$ for $\gamma \to 2$, and for this case ($m=2$), the left and right contributions appear with the same sign.

3.4. Example 4

In the last example, we set the function $g\left(x\right)$ in order to obtain the exponential solution $y\left(x\right)=e^x$. The parameters are set to be ${\rho }_{10}={\zeta }_{10}=1$ with $c_1=1+e$ for $0<\gamma <1$, while for $1<\gamma <2$, we set ${\rho }_{10}={\zeta }_{20}=c_1=1$, $c_2=e$ and ${\rho }_{20}={\zeta }_{10}=0$.

The behavior of the error for different values of $\gamma$ and h is shown in

Table 5. Example 4: The infinity norm of the error $\parallel e_h\parallel$ for different values of $\gamma$ and h.

$\mathit{\gamma }$	0.25	0.75	1.25	1.75
$h=\frac{1}{4}$	$2.45\times {10}^{-5}$	$1.45\times {10}^{-5}$	$1.11\times {10}^{-4}$	$1.31\times {10}^{-4}$
$h=\frac{1}{8}$	$3.20\times {10}^{-6}$	$1.15\times {10}^{-6}$	$1.07\times {10}^{-5}$	$2.53\times {10}^{-5}$
$h=\frac{1}{16}$	$2.83\times {10}^{-7}$	$7.89\times {10}^{-8}$	$1.12\times {10}^{-6}$	$4.91\times {10}^{-6}$
$h=\frac{1}{32}$	$1.92\times {10}^{-8}$	$5.13\times {10}^{-9}$	$1.17\times {10}^{-7}$	$9.72\times {10}^{-7}$
$h=\frac{1}{64}$	$1.31\times {10}^{-9}$	$3.33\times {10}^{-10}$	$1.25\times {10}^{-8}$	$1.96\times {10}^{-7}$
$h=\frac{1}{128}$	$5.80\times {10}^{-10}$	$4.66\times {10}^{-10}$	$1.82\times {10}^{-9}$	$4.02\times {10}^{-8}$

. The numerical results show that the error increases slightly when $\gamma$ increases, while it decreases fast with h.

4. Discussion and Conclusions

In this paper, we discuss the spline collocation method introduced in [12,13] to solve fractional boundary value problems with the Riesz–Caputo derivative. We propose an efficient method for evaluating the fractional derivatives exploiting the symmetry properties of both the spline basis functions and the Riesz–Caputo operator.

We propose several tests in which the good performance of the numerical method is shown. The errors listed in Table 1, Table 2, Table 3, Table 4 and Table 5 show that we obtain an accurate solution even using a small number of collocation points.

In these examples, we also look at the behavior of the left and right contributions of the Riesz–Caputo operator to better understand their roles, and we notice that while the left fractional derivative is always zero at the left edge, the same holds for the right fractional derivative at the opposite edge.

These preliminary studies suggest that the proposed method is accurate and convergent, with convergence order depending on the fractional order of the derivative, in addition to the smoothness of the solution. Moreover, the method can be easily generalized to fractional advection–diffusion problems (see [38]). We plan to further investigate these points in the near future.