Optimal B-Spline Bases for the Numerical Solution of Fractional Differential Problems

Abstract

Efficient numerical methods to solve fractional differential problems are particularly required in order to approximate accurately the nonlocal behavior of the fractional derivative. The aim of the paper is to show how optimal B-spline bases allow us to construct accurate numerical methods that have a low computational cost. First of all, we describe in detail how to construct optimal B-spline bases on bounded intervals and recall their main properties. Then, we give the analytical expression of their derivatives of fractional order and use these bases in the numerical solution of fractional differential problems. Some numerical tests showing the good performances of the bases in solving a time-fractional diffusion problem by a collocation–Galerkin method are also displayed.

1. Introduction

In the last few decades, fractional calculus has proved to be a powerful tool for describing real-world phenomena. Differential problems of fractional order, initially introduced to model anomalous diffusion in viscoelastic materials, are now used in several fields, from physics to population dynamics, from signal processing to control theory [1,2,3,4,5]. For an introduction of fractional calculus, we refer to [6,7,8,9].

As fractional differential models have become widespread, the development of efficient numerical methods to approximate their solution has become of primary interest. In fact, the nonlocality of the fractional derivative poses a severe challenge in its approximation, and, to face this problem, several numerical methods were proposed in the literature. For a review on numerical methods, see, for instance, [10,11,12,13,14,15].

In particular, we are interested in the solution of differential problems on bounded intervals. In this case, it is important to have a function basis available that naturally fulfills boundary and/or initial conditions and does not show numerical instabilities at the boundaries. From this point of view, B-spline bases are especially suitable and their use in the numerical solution of classical differential problems is widely diffused (see, for instance, [16,17,18]). Despite the indubitable success of these methods, the use of B-spline bases for the solution of differential problems of fractional order is not yet very common and limited to few examples. For instance, in [19], the authors used a linear B-spline basis to solve a multi-order fractional differential problem by the operational matrix method, while, in [20], a collocation method based on cubic B-spline wavelets was constructed to solve ordinary differential equations of fractional order.

The aim of this paper is to show how optimal B-spline bases can be profitably used for the solution of fractional differential problems. To this end, we not only describe in detail the construction of optimal B-spline bases of any degree but also give the analytical expression of the basis functions and of their fractional derivatives. These expressions are an efficient tool to construct numerical methods based on spline functions.

The paper is organized as follows. In Section 2, we show how to construct polynomial B-spline bases and recall their main properties. In particular, the construction of cardinal B-splines on the real line through the divided difference operator is described in Section 2.1 and Section 2.2, while, in Section 2.3, we show how to construct optimal B-spline bases on bounded intervals. Their analytical expression in case of equidistant nodes is given in Section 2.4. Section 2.5 is devoted to the evaluation of the fractional derivatives of the B-spline bases. Finally, in Section 2.6, we show how to use optimal bases for the numerical solution of a time-fractional diffusion problem by a collocation–Galerkin method. In Section 3, we use the optimal B-spline basis of degree 3 to solve some test problems; some numerical results are also displayed. Section 4 contains a discussion on the results and highlights possible extensions.

2. Materials and Methods

In this section, we define the B-spline bases through the divided difference operator and recall their main properties. Then, we give the analytical expression of the basis functions and of their fractional derivatives. Finally, B-spline bases are used to solve a fractional differential problem by a collocation–Galerkin method.

2.1. The Divided Difference Operator

Let

$$\mathcal{X}=\{x_1,x_2,\dots ,x_n\},\mathrm{with}{x}_1<x_2<\dots <x_n$$

(1)

be a sequence of simple, i.e., non coincident, knots and let

$$D\left(\begin{array}{cccc}{u}_1\left(x\right),& u_2\left(x\right),& \cdots & u_n\left(x\right)\\ x_1,& x_2,& \cdots & x_n\end{array}\right):=\left|\begin{array}{cccc}{u}_1\left(x_1\right)& u_2\left(x_1\right)& \dots & u_n\left(x_1\right)\\ u_1\left(x_2\right)& u_2\left(x_2\right)& \dots & u_n\left(x_2\right)\\ \dots & \dots & \dots & \dots \\ u_1\left(x_n\right)& u_2\left(x_n\right)& \dots & u_n\left(x_n\right)\end{array}\right|$$

(2)

be the determinant of the collocation matrix

$$M\left(\begin{array}{cccc}{u}_1\left(x\right),& u_2\left(x\right),& \cdots & u_n\left(x\right)\\ x_1,& x_2,& \cdots & x_n\end{array}\right):=\left(\begin{array}{cccc}{u}_1\left(x_1\right)& u_2\left(x_1\right)& \dots & u_n\left(x_1\right)\\ u_1\left(x_2\right)& u_2\left(x_2\right)& \dots & u_n\left(x_2\right)\\ \dots & \dots & \dots & \dots \\ u_1\left(x_n\right)& u_2\left(x_n\right)& \dots & u_n\left(x_n\right)\end{array}\right)$$

(3)

associated with the function system $\mathcal{U}=\{u_1\left(x\right),u_2\left(x\right),\dots ,u_n\left(x\right)\}$ and the knot sequence (1).

The divided difference operator $[x_1,x_2,\dots ,x_n]f$ of a function f over the knots $\mathcal{X}$ is defined as

$$[x_1,x_2,\dots ,x_n]f\left(x\right):=\frac{D\left(\begin{array}{c}1,x,x^2,\dots ,x^{n-2},f\left(x\right)\\ x_1,x_2,x_3,\dots ,x_{n-1},x_n\end{array}\right)}{D\left(\begin{array}{c}1,x,x^2,\dots ,x^{n-2},x^{n-1}\\ x_1,x_2,x_3,\dots ,x_{n-1},x_n\end{array}\right)}.$$

(4)

The formula above can be generalized to the case of coincident knots as follows. Let the knots occur more than once, i.e.,

$$x_1\le x_2\le \dots \le x_n,$$

(5)

and let

$$\mathcal{M}=\{m_1,m_2,\dots ,m_r\}\mathrm{with}\sum _{k=1}^r{m}_k=n$$

(6)

be the multiplicities of the coincident knots, i.e.,

$$\tilde{\mathcal{X}}=\{x_1,x_2,\dots ,x_n\}=\{\underset{m_1}{\underbrace{{\eta }_1,\dots {\eta }_1}},\underset{m_2}{\underbrace{{\eta }_2,\dots {\eta }_2}},\dots ,\underset{m_r}{\underbrace{{\eta }_r,\dots {\eta }_r}}\},$$

(7)

where ${\eta }_1<{\eta }_2<\dots <{\eta }_r$ are the non coincident knots of the sequence $\tilde{\mathcal{X}}$.

Assuming the functions in $\mathcal{U}$ are sufficiently smooth, the collocation matrix associated with the function system $\mathcal{U}$ over the knot sequence (7) is given by

$$\tilde{M}\left(\begin{array}{cccc}{u}_1\left(x\right),& u_2\left(x\right),& \cdots & u_n\left(x\right)\\ x_1,& x_2,& \cdots & x_n\end{array}\right):=\left(\begin{array}{cccc}{u}_1\left({\eta }_1\right)& u_2\left({\eta }_1\right)& \dots & u_n\left({\eta }_1\right)\\ u_1^{\prime }\left({\eta }_1\right)& u_2^{\prime }\left({\eta }_1\right)& \dots & u_n^{\prime }\left({\eta }_1\right)\\ \dots & \dots & \dots & \dots \\ u_1^{(m_1-1)}\left({\eta }_1\right)& u_2^{(m_1-1)}\left({\eta }_1\right)& \dots & u_n^{(m_1-1)}\left({\eta }_1\right)\\ \dots & \dots & \dots & \dots \\ u_1\left({\eta }_r\right)& u_2\left({\eta }_r\right)& \dots & u_n\left({\eta }_r\right)\\ u_1^{\prime }\left({\eta }_r\right)& u_2^{\prime }\left({\eta }_r\right)& \dots & u_n^{\prime }\left({\eta }_r\right)\\ \dots & \dots & \dots & \dots \\ u_1^{(m_r-1)}\left({\eta }_r\right)& u_2^{(m_r-1)}\left({\eta }_r\right)& \dots & u_n^{(m_r-1)}\left({\eta }_r\right)\end{array}\right).$$

(8)

Thus, if the function f has sufficient derivatives, the divided difference operator of f over the sequence of multiple knots $\tilde{\mathcal{X}}$ can be defined by using the collocation matrix (8) in Definition (4), i.e.,

$$[x_1,x_2,\dots ,x_n]f\left(x\right):=\frac{\tilde{D}\left(\begin{array}{c}1,x,x^2,\dots ,x^{n-2},f\left(x\right)\\ x_1,x_2,x_3,\dots ,x_{n-1},x_n\end{array}\right)}{\tilde{D}\left(\begin{array}{c}1,x,x^2,\dots ,x^{n-2},x^{n-1}\\ x_1,x_2,x_3,\dots ,x_{n-1},x_n\end{array}\right)},$$

(9)

where $\tilde{D}$ denotes the determinant of the collocation matrix as defined in Equation (8).

2.2. The Cardinal B-Splines

The polynomial splines are piecewise polynomials having a certain degree of smoothness which is related to the degree of the polynomial pieces. Spline functions can be represented as a linear combination of B-splines that form a local basis for the spline spaces [21].

The cardinal B-splines, i.e., the polynomial B-splines having break points on integer knots, can be defined by applying the divided difference operator (4) to the truncated power function defined as

$$x_{+}^n:=\left\{\begin{array}{ccc}{x}^n,\hfill & \mathrm{for}& x\ge 0,\hfill \\ 0,\hfill & \mathrm{for}& x<0.\hfill \end{array}\right.$$

(10)

Then, the cardinal B-spline $B_n$ of degree n on the integer knots

$$\mathcal{I}=\{0,1,\dots ,n+1\}$$

(11)

is given by

$$\begin{array}{ccc}{B}_n\left(x\right)\hfill & :=& (n+1)[0,1,\dots ,n+1]{(y-x)}_{+}^n\hfill \\ & =& {\displaystyle (n+1)\frac{D\left(\begin{array}{c}1,y,y^2,\dots ,y^n,{(y-x)}_{+}^n\\ 0,1,2,\dots ,n,n+1\end{array}\right)}{D\left(\begin{array}{c}1,y,y^2,\dots ,y^n,y^{n+1}\\ 0,1,2,\dots ,n,n+1\end{array}\right)}=\frac{1}{n!}{\Delta }^{n+1}{x}_{+}^n,}\hfill \end{array}$$

(12)

where ${\Delta }^n$ is the finite difference operator

$${\Delta }^{n}f\left(x\right)=\sum _{k=0}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{n}{k}\right)f(x-k).$$

(13)

From the above definition, it follows that the cardinal B-spline $B_n$ is compactly supported on $[0,n+1]$, positive in $(0,n+1)$ and belongs to $C^{n-1}\left(\mathbb{R}\right)$. Moreover, the system of the integer translates

$${\mathcal{B}}_n=\left\{B_n(x-k),k\in \mathbb{Z}\right\}$$

(14)

forms a function basis that is a partition of unity, i.e.,

$$\sum _{k\in \mathbb{Z}}{B}_n(x-k)=1,\forall x\in \mathbb{R},$$

(15)

and reproduces polynomials up to degree n [21].

The basis ${\mathcal{B}}_n$ can be generalized to any set of equidistant knots ${\mathcal{X}}_h=\{x_k=kh,k\in \mathbb{Z}\}$, where $h>0$ is the space step, by scaling, i.e.,

$${\mathcal{B}}_{nh}=\left\{B_n\left(\frac{x}{h}-k\right),k\in \mathbb{Z}\right\}.$$

(16)

For details and further properties on polynomial B-splines, see, for instance, [21].

It is always possible to construct bases on bounded intervals by restricting the B-spline basis ${\mathcal{B}}_{nh}$ on an interval $[a,b]$. However, in this way, we obtain bases that could be numerically unstable since the functions at the boundaries are obtained by truncation. Moreover, boundary or initial conditions are not easy to satisfy. In the following section, we will show how to construct stable bases on the interval that easily satisfy boundary or initial conditions.

2.3. Optimal B-Spline Bases

Optimal B-spline bases on bounded intervals can be defined by introducing multiple knots at the boundaries of the given interval [21].

Let $L\ge n+1$ be an integer and let $\mathcal{I}=\{0,1,\dots ,L-1,L\}$ be the sequence of the integer knots of the interval $I=[0,L]$. We extend $\mathcal{I}$ to a sequence $\tilde{\mathcal{I}}$ of $L+1+2n$ knots by introducing knots of multiplicity $n+1$ at the boundaries of the interval, i.e.,

$$\tilde{\mathcal{I}}=\{x_0,x_1,\dots ,x_{2n+L}\}\mathrm{with}\left\{\begin{array}{c}{x}_0=x_1=\cdots =x_n=0,\hfill \\ x_k=k-n,n+1\le k\le L+n-1,\hfill \\ x_{L+n}=x_{L+n+1}=\cdots =x_{L+2n}=L.\hfill \end{array}\right.$$

(17)

The functions of the optimal B-spline basis

$${\mathcal{N}}_n=\{N_{kn}\left(x\right),0\le k\le L-1+n\}$$

(18)

of degree n with knots $\tilde{\mathcal{I}}$ are given by the divided difference

$$N_{kn}\left(x\right)=(x_{k+n+1}-x_k)[x_k,x_{k+1},\dots ,x_{k+n+1}]{(y-x)}_{+}^n,0\le k\le L-1+n,$$

(19)

so that each function $N_{kn}$ has compact support with

$$\mathrm{supp}{N}_{kn}=[x_k,x_{k+n+1}]\subset I.$$

(20)

From (19), it follows that the basis ${\mathcal{N}}_n$ has $L-n$ interior functions, i.e., the functions $N_{kn}$, $n\le k\le L-1$ having all simple knots, and $2n$ edge functions, i.e., the functions $N_{kn}$ for $0\le k\le n-1$ and $L\le k\le L-1+n$ that have a multiple knot at the left or right boundary, respectively.

The basis ${\mathcal{N}}_n$ is centrally symmetric, i.e.,

$$N_{kn}\left(x\right)=N_{(L-1+n-k)n}(L-x),0\le k\le L-1+n.$$

(21)

Moreover, the functions $N_{kn}$ satisfy the boundary conditions

$$\left\{\begin{array}{ccc}{N}_{kn}\left(0^{+}\right)={\delta }_{k0},\hfill & & 0\le k\le L-1+n,\hfill \\ \\ D_x^r{N}_{kn}\left(0^{+}\right)=0,\hfill & 1\le r\le k-1,\hfill & \\ & & 1\le k\le n-1,\hfill \\ {(-1)}^r{D}_x^r{N}_{kn}\left(0^{+}\right)>0,\hfill & k\le r\le n,\hfill & \end{array}\right.$$

(22)

$$\left\{\begin{array}{ccc}{N}_{kn}\left(L^{-}\right)={\delta }_{k,r-1},\hfill & & 0\le k\le L-1+n,\hfill \\ \\ D_x^r{N}_{(L-1+n-k)n}\left(L^{-}\right)=0,\hfill & 1\le r\le k-1,\hfill & \\ & & 1\le k\le n-1,\hfill \\ {(-1)}^r{D}_x^r{N}_{(L-1+n-k)n}\left(L^{-}\right)>0,\hfill & k\le r\le n,\hfill & \end{array}\right.$$

(23)

where $D_x^r$ denotes the usual derivative of integer order r.

Finally, the basis ${\mathcal{N}}_n$ forms a partition of unity, i.e.,

$$\sum _{k=0}^{L-1+n}{N}_{kn}\left(x\right)=1,\forall x\in I,$$

(24)

and is stable, i.e., for any sequence $c=\{c_0,c_1,\dots ,c_{L-1+n}\}$ with ${\parallel c\parallel }_{\infty }={max}_{0\le k\le L-1+n}\left|c_k\right|<\infty$ it holds

$$\frac{1}{\kappa }{\parallel c\parallel }_{\infty }\le {∥\sum _{k=0}^{L-1+n}{c}_k{N}_{kn}∥}_{L_{\infty }[0,L]}\le {\parallel c\parallel }_{\infty },$$

(25)

where

$$\kappa :={\left({\displaystyle \underset{{\parallel c\parallel }_{\infty }=1}{min}{∥\sum _{k=0}^{L-1+n}{c}_k{N}_{kn}∥}_{L_{\infty }[0,L]}}\right)}^{-1}$$

(26)

is the condition number of the basis ${\mathcal{N}}_n$. From (25), it follows that $\parallel {\sum }_k\Delta c_k{N}_{k3}{\parallel }_{L_{\infty }[0,L]}\le {\parallel \Delta c\parallel }_{\infty }$ where $\Delta c$ denotes a perturbation of the sequence c. As a consequence, optimal B-spline bases are numerically stable so that numerical errors are not amplified when evaluating spline approximations.

The basis ${\mathcal{N}}_n$ can be generalized to any sequence of equidistant knots on a bounded interval $[a,b]$ by mapping $x\to (x-a)/h$, i.e.,

$${\mathcal{N}}_{nh}=\{N_{knh}\left(x\right),0\le k\le L-1+n\},$$

(27)

where $L=(b-a)/h$. The interior functions are the $L-n$ functions

$$N_{knh}\left(x\right)=B_n\left(\frac{x-a}{h}-(k-n)\right),n\le k\le L-1,$$

(28)

while the $2n$ edge functions are

$$N_{knh}\left(x\right)=N_{kn}\left(\frac{x-a}{h}\right),N_{(L-1+n-k)nh}\left(x\right)=N_{knh}\left(b-x\right),0\le k\le n-1.$$

(29)

All the properties above make the basis ${\mathcal{N}}_n$ optimal when used in approximation problems [22,23].

2.4. Analytical Expression of the Optimal B-Spline Bases

Let ${\mathcal{N}}_n$ be the optimal basis (18) associated with the knot sequence (17). The interior functions $N_{kn}$, $n\le k\le L-1$, are the translates of the B-spline $B_n$ having support $[k-n,k+1]$, i.e.,

$$N_{kn}\left(x\right)=B_n(x-k+n),n\le k\le L-1.$$

(30)

The left edge functions $N_{kn}$, $0\le k\le n-1$, having support on $[0,k+1]$, can be evaluated by using the finite difference operator (9) for coincident knots. For $0\le k\le n-1$, we get

$$\begin{array}{ccc}{N}_{kn}\left(x\right)\hfill & =& (k+1)[x_k,x_{k+1},\dots ,x_{k+n+1}]{(y-x)}_{+}^n\hfill \\ & =& {\displaystyle (k+1)\frac{\tilde{D}\left(\begin{array}{cccccccc}1,& y,& \cdots & \cdots & \cdots & \cdots & y^n,& {(y-x)}_{+}^n\\ {\underbrace{0}}_1,& {\underbrace{0}}_2,& \cdots ,& {\underbrace{0}}_{n+1-k},& 1,& \cdots & k,& k+1\end{array}\right)}{\tilde{D}\left(\begin{array}{cccccccc}1,& y,& \cdots & \cdots & \cdots & \cdots & y^n,& y^{n+1}\\ {\underbrace{0}}_1,& {\underbrace{0}}_2,& \cdots ,& {\underbrace{0}}_{n+1-k},& 1,& \cdots & k,& k+1\end{array}\right)}.}\hfill \end{array}$$

(31)

By Definition (8), we get

$$\tilde{D}\left(\begin{array}{cccccccc}1,& y,& \cdots & \cdots & \cdots & \cdots & y^n,& {(y-x)}_{+}^n\\ \underset{1}{\underbrace{0}},& \underset{2}{\underbrace{0}},& \cdots ,& \underset{n+1-k}{\underbrace{0}},& 1,& \cdots & k,& k+1\end{array}\right)=\left|\begin{array}{cc}{D}_{kn}& 0\\ R_{kn}& T_{kn}\end{array}\right|,$$

(32)

where $D_{kn}$ is the $(n-k+1)$ order diagonal matrix

$$D_{kn}=\mathrm{diag}(1,1,2!,\dots ,(n-k)!),$$

(33)

$R_{kn}$ is the $(k+1)\times (n-k+1)$ dimension matrix

$$R_{kn}=\left(i^{j-1},1\le i\le k+1,1\le j\le n-k+1\right),$$

(34)

and $T_{kn}$ is the $(k+1)$ order collocation matrix

$$T_{kn}=M\left(\begin{array}{ccccc}{y}^{n-k+1},& y^{n-k+2},& \cdots & y^n,& {(y-x)}_{+}^n\\ 1,& 2,& \cdots & k,& k+1\end{array}\right).$$

(35)

Thus, for the numerator in the last equality of (31), we have

$$\tilde{D}\left(\begin{array}{cccccccc}1,& y,& \cdots & \cdots & \cdots & \cdots & y^n,& {(y-x)}_{+}^n\\ \underset{1}{\underbrace{0}},& \underset{2}{\underbrace{0}},& \cdots ,& \underset{n+1-k}{\underbrace{0}},& 1,& \cdots & k,& k+1\end{array}\right)=|T_{nk}\left|\right|D_{kn}|.$$

(36)

For the denominator, a similar calculation gives

$$\tilde{D}\left(\begin{array}{cccccccc}1,& y,& \cdots & \cdots & \cdots & \cdots & y^n,& y^{n+1}\\ \underset{1}{\underbrace{0}},& \underset{2}{\underbrace{0}},& \cdots ,& \underset{n+1-k}{\underbrace{0}},& 1,& \cdots & k,& k+1\end{array}\right)=|P_{nk}\left|\right|D_{kn}|,$$

(37)

where $P_{kn}$ is the $k+1$ order collocation matrix

$$P_{kn}=M\left(\begin{array}{ccccc}{y}^{n-k+1},& y^{n-k+2},& \cdots & y^n,& y^{n+1}\\ 1,& 2,& \cdots & k,& k+1\end{array}\right).$$

(38)

Theorem 1.

Let $\tilde{\mathcal{I}}$ be the sequence of multiple knots (17) on the interval $I=[0,L]$. The analytical expression of the left edge functions is given by

$$N_{kn}\left(x\right)=(k+1)\frac{|T_{nk}|}{|P_{kn}|},0\le k\le n-1,$$

(39)

where $|T_{nk}|$ and $|P_{kn}|$ are the determinants of the collocation matrices (35) and (38), respectively.

Proof. 

The analytical expression (39) follows by using (36) and (37) in the last equality in (31). ☐

The right edge functions $N_{L-1+n-k,n}$, $0\le k\le n-1$, having support $[L-k-1,L]$, can be easily obtained recalling the symmetry property (21).

In the following corollary, we give the analytical expression of the first three left boundary functions.

Corollary 1.

Let $\tilde{\mathcal{I}}$ be the sequence of multiple knots (17) on the interval $I=[0,L]$. For $0\le k\le 2$, the left edge functions $N_{kn}$ are given by

$$N_{0n}\left(x\right)={(1-x)}_{+}^n=\{\begin{array}{ll}{\displaystyle {(1-x)}^n,}\hfill & 0\le x\le 1,\hfill \\ 0,\hfill & \mathit{otherwise},\hfill \end{array}\hspace{1em}\hspace{1em}n\ge 1,$$

(40)

$$N_{1n}\left(x\right)=2\left(\frac{1}{2^n}{(2-x)}_{+}^n-{(1-x)}_{+}^n\right)=\{\begin{array}{ll}{\displaystyle \frac{1}{2^{n-1}}{(2-x)}^n-2{(1-x)}^n,}\hfill & 0\le x\le 1,\hfill \\ {\displaystyle \frac{1}{2^{n-1}}{(2-x)}^n,}\hfill & 1<x\le 2,\hfill \\ 0,\hfill & \mathit{otherwise},\hfill \end{array}\hspace{1em}\hspace{1em}n\ge 2,$$

(41)

$$\begin{array}{ll}{N}_{2n}\left(x\right)& =3\left(\frac{1}{2}\frac{1}{3^{n-1}}{(3-x)}_{+}^n-\frac{1}{2^{n-1}}{(2-x)}_{+}^n+\frac{1}{2}{(1-x)}_{+}^n\right)\hfill \\ & =\left\{\begin{array}{cc}{\displaystyle \frac{1}{2}\frac{1}{3^{n-2}}{(3-x)}^n-\frac{3}{2^{n-1}}{(2-x)}^n+\frac{3}{2}{(1-x)}^n,}\hfill & 0\le x\le 1,\hfill \\ {\displaystyle \frac{1}{2}\frac{1}{3^{n-2}}{(3-x)}^n-\frac{3}{2^{n-1}}{(2-x)}^n,}\hfill & 1<x\le 2,\hfill \\ {\displaystyle \frac{1}{2}\frac{1}{3^{n-2}}{(3-x)}^n,}\hfill & 2<x\le 3,\hfill \\ 0,\hfill & \mathit{otherwise},\hfill \end{array}\right.n\ge 3.\hfill \end{array}$$

(42)

It is interesting to observe that ${\displaystyle N_{0n}\left(x\right)=\frac{B_n(x+n)}{B_n\left(n\right)}}$, $0\le x\le 1$. In particular, when $n=1$ (linear case), there is just one left edge function $N_{01}\left(x\right)=B_1(x+1)$ and the basis ${\mathcal{N}}_1$ coincides with the basis ${\mathcal{B}}_1$ restricted to the interval I. We recall that the basis ${\mathcal{N}}_1$ is the unique interpolatory B-spline basis.

2.5. Fractional Derivatives of the Optimal B-Spline Bases

In this section, we give the explicit expression of the derivatives of fractional order of the functions in the basis ${\mathcal{N}}_n$.

Let $\Gamma$ be the Euler’s gamma function

$$\Gamma \left(\beta \right)={\int }_0^{\infty }{z}^{\beta -1}{\mathrm{e}}^{-z}dz.$$

(43)

Definition 1.

Assuming f is a sufficiently smooth function, the Caputo fractional derivative of f is defined as

$${}^C{D}_x^{\beta }f\left(x\right):=\frac{1}{\Gamma \left(m-\beta \right)}{\int }_0^x\frac{f^{\left(m\right)}\left(z\right)}{{(x-z)}^{\beta -m+1}}dz,x\ge 0,m-1<\beta <m,m\in \mathbb{N}.$$

(44)

Definition 2.

The Riemann–Liouville fractional derivative of a function f is defined as

$${}^{RL}{D}_x^{\beta }f\left(x\right):=\frac{1}{\Gamma \left(m-\beta \right)}\frac{d^m}{d{x}^m}{\int }_0^x\frac{f\left(z\right)}{{(x-z)}^{\beta -m+1}}dz,x\ge 0,m-1<\beta <m,m\in \mathbb{N}.$$

(45)

We recall that the Caputo derivative and the Riemann–Liouville derivative can be obtained from each other by

$${}^{RL}{D}_x^{\beta }f\left(x\right)={}^C{D}_x^{\beta }f\left(x\right)+\sum _{k=0}^{m-1}\frac{x^{k-\beta }}{\Gamma (k-\beta +1)}{f}^{\left(k\right)}\left(0^{+}\right).$$

(46)

If f satisfies homogeneous initial conditions, i.e., $f^{\left(k\right)}\left(0^{+}\right)=0$, $0\le k\le m-1$, the Riemann–Liouville derivative coincides with the Caputo derivative.

To use the basis ${\mathcal{N}}_n$ for the solution of differential problems of fractional order, we need the expression of the fractional derivatives of the functions $N_{kn}$.

The derivatives of the interior functions can be easily evaluated by the differentiation rule

$${}^C{D}_x^{\beta }{B}_n\left(x\right)=\frac{{\Delta }^{n+1}{x}_{+}^{n-\beta }}{\Gamma (n-\beta +1)},x\ge 0,0<\beta <n,$$

(47)

where ${\Delta }^n$ is the finite difference operator (13) [15,24]. Since $B_n$ satisfies homogeneous boundary conditions, the Riemann–Liouville derivative is equal to the Caputo derivative.

Theorem 2.

The Caputo derivative of the left edge functions are given by

$${}^C{D}_x^{\beta }{N}_{kn}\left(x\right)=\frac{(k+1)}{|P_{kn}|}\sum _{r=1}^{k+1}{(-1)}^{r+k+1}{}^C{D}_x^{\beta }{(r-x)}_{+}^n\left|T_{kn}^r\right|,0\le k\le n-1,$$

(48)

where

$$T_{kn}^r=\left(\begin{array}{cccc}1& 1& \dots & 1\\ 2^{n-k+1}& 2^{n-k+2}& \dots & 2^n\\ \dots & \dots & \dots & \dots \\ {(r-1)}^{n-k+1}& {(r-1)}^{n-k+2}& \dots & {(r-1)}^n\\ {(r+1)}^{n-k+1}& {(r+1)}^{n-k+2}& \dots & {(r+1)}^n\\ \dots & \dots & \dots & \dots \\ {(k+1)}^{n-k+1}& {(k+1)}^{n-k+2}& \dots & {(k+1)}^n\end{array}\right).$$

(49)

Proof. 

First of all, we write the explicit expression of the matrix $T_{kn}$, i.e.,

$$T_{kn}=\left(\begin{array}{ccccc}1& 1& \dots & 1& {(1-x)}_{+}^n\\ 2^{n-k+1}& 2^{n-k+2}& \dots & 2^n& {(2-x)}_{+}^n\\ \dots & \dots & \dots & \dots & \dots \\ {(k+1)}^{n-k+1}& {(k+1)}^{n-k+2}& \dots & {(k+1)}^n& {(k+1-x)}_{+}^n\end{array}\right).$$

From (39), we get ${}^C{D}_x^{\beta }{N}_{kn}\left(x\right)=\frac{k+1}{|P_{kn}|}{}^C{D}_x^{\beta }\left|T_{kn}\right|.$ The claim follows by expanding the determinant of $T_{kn}$ along the last column. ☐

Thus, to evaluate the fractional derivatives of the edge functions, we need the fractional derivatives of the translates of the truncated power function. By Definition (44), we get

$$\begin{array}{ccc}{}^C{D}_x^{\beta }{(r-x)}_{+}^n\hfill & =& {\displaystyle \frac{1}{\Gamma (m-\beta )}{\int }_0^x\frac{\frac{d^m}{d{z}^m}{(r-z)}_{+}^n}{{(x-z)}^{\beta -m+1}}dz}\hfill \\ & =& {\displaystyle \frac{1}{\Gamma (m-\beta )}\frac{{(-1)}^{m}n!}{(n-m)!}{\int }_0^x\frac{{(r-z)}_{+}^{n-m}}{{(x-z)}^{\beta -m+1}}dz.}\hfill \end{array}$$

(50)

Integration rules for rational functions give [25]

$$\begin{array}{c}{\int }_0^x\frac{{(r-z)}^{n-m}}{{(x-z)}^{\beta -m+1}}dz=\frac{1}{m-\beta }\frac{r^{n-m}}{x^{\beta -m}}{}_2{F}_1(1,m-n,m+1-\beta ,\frac{x}{r})\hfill \\ {\displaystyle \hspace{1em}\hspace{1em}+\frac{\Gamma (n-m+2)\Gamma (m-\beta )}{\Gamma (n-\beta +1)}{\left(-\frac{x}{r}\right)}^{\beta -m+1}\Theta (r-x){}_2{F}_1(\beta -m+1,\beta -n,\beta -m+1,\frac{x}{r}),}\hfill \end{array}$$

(51)

where $\Theta \left(z\right)$ is the Heaviside function and ${}_2{F}_1(a,b,c,z)$ is the hypergeometric function

$${}_2{F}_1(a,b,c,z)=\sum _{k=0}^{\infty }\frac{{\left(a\right)}_k{\left(b\right)}_k}{{\left(c\right)}_k}\frac{z^k}{k!},$$

(52)

with ${\left(q\right)}_k$ denoting the rising Pochhammer symbol

$${\left(q\right)}_k=\left\{\begin{array}{cc}1,\hfill & k=0,\hfill \\ q(q+1)\cdots (q+k-1),\hfill & k>0.\hfill \end{array}\right.$$

(53)

Since $m-n<0$, ${}_2{F}_1(1,m-n,m+1-\beta ,\frac{x}{r})$ is actually a finite sum, i.e.,

$${}_2{F}_1(1,m-n,m+1-\beta ,\frac{x}{r})=\sum _{k=0}^{n-m}{(-1)}^k\left(\genfrac{}{}{0pt}{}{n-m}{k}\right)\frac{{\left(1\right)}_k}{{(m+1-\beta )}_k}{\left(\frac{x}{r}\right)}^k.$$

(54)

Moreover, ${}_2{F}_1(\beta -m+1,\beta -n,\beta -m+1,\frac{x}{r})$ can be evaluated by a direct calculation, i.e.,

$${}_2{F}_1(\beta -m+1,\beta -n,\beta -m+1,\frac{x}{r})=\sum _{k=0}^{\infty }\frac{{(\beta -n)}_k}{k!}{\left(\frac{x}{r}\right)}^k={\left(\frac{r-x}{r}\right)}^{n-\beta }.$$

(55)

Thus, we get the analytical expression of the Caputo derivative (cf. also [20]).

Theorem 3.

For $0<\beta <n$ with $m-1<\beta <m$ and m a positive integer, the Caputo derivative of the translates of the truncated power function is given by

$$\begin{array}{cc}{}^C{D}_x^{\beta }{(r-x)}_{+}^n\hfill & {\displaystyle =\frac{1}{\Gamma (m+1-\beta )}\frac{{(-1)}^{m}n!}{(n-m)!}\frac{r^{n-m}}{x^{\beta -m}}\sum _{k=0}^{n-m}{(-1)}^k\left(\genfrac{}{}{0pt}{}{n-m}{k}\right)\frac{{\left(1\right)}_k}{{(m+1-\beta )}_k}{\left(\frac{x}{r}\right)}^k}\hfill \\ & {\displaystyle +\frac{{(-1)}^{n+m}n!}{\Gamma (n+1-\beta )}{(x-r)}_{+}^{n-\beta }.}\hfill \end{array}$$

(56)

The Riemann–Liouville derivatives of the functions $N_{kn}$ can be obtained by the Caputo derivatives using the initial conditions (22) in Equation (46).

2.6. The Collocation-Galerkin Method

In this section, we show how to use the optimal basis introduced in the previous sections to solve a fractional differential problem, i.e., the time-fractional diffusion problem (cf. [26,27,28])

$$\left\{\begin{array}{ccc}{\displaystyle D_t^{\beta }u(x,t)-\frac{{\partial }^2}{\partial x^2}u(x,t)=f(x,t),}\hfill & x\in [0,\eta ],& t\in [0,\tau ],\\ u(x,0)=u_0,\hfill & x\in [0,\eta ],& \\ u(0,t)=u(\eta ,t)=0,\hfill & & t\in [0,\tau ],\end{array}\right.$$

(57)

where $D_t^{\beta }u$, $0<\beta <1$, denotes the partial Caputo derivative with respect to the time t.

Let

$${\mathcal{V}}_{nh}^0=\overline{\mathrm{span}}\left\{N_{knh}\left(x\right),1\le k\le L-2+n\right\},$$

(58)

where $L=\eta /h$ is the approximating space generated by the optimal B-spline basis

$${\mathcal{N}}_{nh}^0=\left\{N_{knh}\left(x\right),1\le k\le L-2+n\right\},$$

(59)

of the interval $[0,\eta ]$ on equidistant knots with space step $h>0$ (cf. (27)). We observe that the basis ${\mathcal{N}}_{nh}^0$ naturally satisfies the homogeneous boundary conditions

$$N_{knh}\left(0\right)=N_{knh}\left(\eta \right)=0,1\le k\le L-2+n.$$

(60)

For any h held fix, the spline Galerkin method looks for an approximating function

$$u_h(x,t)=\sum _{k=1}^{L-2+n}{c}_k\left(t\right)N_{knh}\left(x\right)\in {\mathcal{V}}_{nh}^0$$

(61)

that solves the variational problem

$$\left\{\begin{array}{cc}{\displaystyle \left(D_t^{\beta }{u}_h,N_{knh}\right)-\left(\frac{\partial }{\partial x^2}{u}_h,N_{knh}\right)=\left(f,N_{knh}\right),}\hfill & 1\le k\le L-2+n,\hfill \\ u_h(x,0)=u_0,\hfill & x\in [0,\eta ],\hfill \end{array}\right.$$

(62)

where $(f,g)={\int }_0^{\eta }fg$. Writing (62) in a weak form and using (61), we get the system of fractional ordinary differential equations

$$Q{D}_t^{\beta }C\left(t\right)+LC\left(t\right)=F\left(t\right),t\in [0,\tau ],$$

(63)

with initial condition

$$\sum _{k=1}^{L-2+n}{c}_k\left(0\right)N_{knh}\left(x\right)=u_0,x\in [0,\eta ].$$

(64)

Here, $C\left(t\right)={\left(c_k\left(t\right)\right)}_{1\le k\le L-2+n}$ is the unknown function vector, $Q={\left(q_{kj}\right)}_{1\le k,j\le L-2+n}$ is the mass matrix, $L={\left({\ell }_{kj}\right)}_{1\le k,j\le L-2+n}$ is the stiffness matrix and $F={\left(f_k\left(t\right)\right)}_{1\le k\le L-2+n}$ is the load vector whose entries are

$$q_{kj}={\int }_0^{\eta }{N}_{knh}{N}_{jnh},{\ell }_{kj}={\int }_0^{\eta }{N}_{knh}^{\prime }{N}_{jnh}^{\prime },f_k\left(t\right)={\int }_0^{\eta }f(t,·)N_{knh}.$$

The entries of Q and L can be evaluated explicitly by using the integration and differentiation rules for B-splines [21]. The entries of F can be evaluated by quadrature formulas especially designed for Galerkin methods [29].

Now, we introduce the sequence of temporal knots. Let

$$\tilde{\mathcal{T}}=\{t_0,t_1,\dots ,t_{T+m}\}\mathrm{with}\left\{\begin{array}{c}{t}_0=t_1=\cdots =t_m=0,\hfill \\ t_r=(r-m)s,m+1\le r\le T+m,\hfill \end{array}\right.$$

(65)

where $T=\tau /s$ with $s>0$ the time step, be a set of equidistant nodes in the interval $[0,\tau ]$ having a knots of multiplicity $m+1$ at the left boundary. We assume the unknown functions $c_k\left(t\right)$ belong to the spline space

$${\mathcal{V}}_{ms}^1=\mathrm{span}\left\{N_{rms}\left(x\right),0\le r\le T-1+m\right\},$$

(66)

so that

$$c_k\left(t\right)=u_0{N}_{0ms}\left(t\right)+\sum _{r=1}^{T-1+m}{\lambda }_{rk}{N}_{rms}\left(t\right),1\le k\le L-2+m.$$

(67)

We observe that just the functions $N_{rms}$, $0\le r\le n$, are boundary functions while all the other functions in the basis are cardinal B-splines.

To solve the fractional differential system (63), we use the collocation method on equidistant nodes introduced in [28,30]. Let $\mathcal{P}=\{t_p,1\le p\le P\}$ with $P\ge T-1+m$ be a set of non coincident collocation points in the interval $(0,\tau ]$. Then, collocating (63) on the nodes $t_p$, we get the linear system

$$(Q\otimes A+L\otimes B)\Lambda =F,$$

(68)

where $\Lambda ={\left({\lambda }_{rk}\right)}_{1\le r\le T-1+m,1\le k\le L-2+n}$ is the unknown vector,

$$A={\left(a_{pr}\right)}_{1\le p\le P,1\le r\le T-1+m},a_{pr}=D_t^{\beta }{N}_{rms}\left(t_p\right),$$

(69)

$$B={\left(b_{pr}\right)}_{1\le p\le P,1\le r\le T-1+m},b_{pr}=N_{rms}\left(t_p\right),$$

(70)

are collocation matrices, and $F={\left(f_k\left(t_p\right)\right)}_{1\le k\le L-2+n,1\le p\le P}$. When $P=T-1+m$, (68) is a square linear system and the collocation matrix A is non singular if and only if [21],

$$t_p\in \mathrm{supp}{N}_{pms}$$

(71)

while B is always non singular [31]. Otherwise, (68) is an overdetermined linear system that can be solved in the least squares sense.

We notice that we must pay a special attention to the evaluation of the entries of A since they involve the values of the fractional derivative $D_t^{\beta }{N}_{rms}$ in the collocation points. This can be done by the differentiation rules given in Section 2.5.

3. Results

To test the performance of the optimal B-spline basis when used to solve fractional differential problems, we solved the time-fractional diffusion problem (57) by the collocation–Galerkin method described in Section 2.6. In the tests, we used the optimal basis ${\mathcal{N}}_3$ having degree 3 since the cubic B-spline has a small support but is sufficiently smooth. In fact, the cubic B-spline belongs to $C^2\left(I\right)$ and its support has length 4. In particular, the basis ${\mathcal{N}}_3$ has 3 right (left) edge functions with support $[0,k]$ ($[L-4+k,L]$), $1\le k\le 3$, while the interior functions have support $[k,k+4]$, $0\le k\le L-4$. In the following section, we give the analytical expression of the functions $N_{k3}$ and of their fractional derivatives; then, we show some numerical results.

3.1. The Optimal B-Spline Basis of Degree $n=3$

When $n=3$, we get the cubic B-spline basis. The optimal B-spline basis ${\mathcal{N}}_3$ of the interval $[0,L]$, with L an integer greater than 3, is associated with the sequence of integer knots

$$\tilde{\mathcal{I}}=\{x_0,x_1,\dots ,x_{L+6}\}$$

(72)

with $x_0=x_1=x_2=x_3=0$, $x_k=k-3$, $4\le k\le L+2$, $x_{L+3}=x_{L+4}=x_{L+5}=x_{L+6}=L$. ${\mathcal{N}}_3$ has $L-3$ interior functions and six edge functions, three for each boundary.

The interior functions $N_{3k}$, $3\le k\le L-1$, are the translates of the cubic cardinal B-spline, i.e.,

$$N_{3k}\left(x\right)=B_3(x-k+3),3\le k\le L-1,$$

(73)

where

$$B_3\left(x\right)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{6}{x}^3,}\hfill & 0\le x\le 1,\hfill \\ {\displaystyle \frac{1}{6}\left(-3x^3+12x^2-12x+4\right),}\hfill & 1\le x\le 2,\hfill \\ {\displaystyle \frac{1}{6}\left(-3{(4-x)}^3+12{(4-x)}^2-12(4-x)+4\right),}\hfill & 2\le x\le 3,\hfill \\ {\displaystyle \frac{1}{6}{(4-x)}^3,}\hfill & 3\le x\le 4,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(74)

From Corollary 1, it follows that the three left edge functions $N_{kn}$, $0\le k\le 2$, are given by

$$N_{03}\left(x\right)={(1-x)}_{+}^3=\{\begin{array}{cc}{\displaystyle {(1-x)}^3,}\hfill & 0\le x\le 1,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}$$

(75)

$$\begin{array}{cc}{N}_{13}\left(x\right)& =\frac{1}{4}{(2-x)}_{+}^3-2{(1-x)}_{+}^3\hfill \\ & =\{\begin{array}{cc}{\displaystyle \frac{1}{4}{(2-x)}^3-2{(1-x)}^3,}\hfill & 0\le x\le 1,\hfill \\ {\displaystyle \frac{1}{4}{(2-x)}^3,}\hfill & 1<x\le 2,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\hfill \end{array}$$

(76)

$$\begin{array}{cc}{N}_{23}\left(x\right)\hfill & {\displaystyle =\frac{1}{6}{(3-x)}_{+}^3-\frac{3}{4}{(2-x)}_{+}^3+\frac{3}{2}{(1-x)}_{+}^3}\hfill \\ & =\begin{array}{c}\left\{\begin{array}{cc}{\displaystyle \frac{1}{6}{(3-x)}^3-\frac{3}{4}{(2-x)}^3+\frac{3}{2}{(1-x)}^3,}\hfill & 0\le x\le 1,\hfill \\ {\displaystyle \frac{1}{6}{(3-x)}^3-\frac{3}{4}{(2-x)}^3,}\hfill & 1<x\le 2,\hfill \\ {\displaystyle \frac{1}{6}{(3-x)}^3,}\hfill & 2<x\le 3,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}\hfill \end{array}$$

(77)

Now, using Theorems 3, we obtain the analytical expression of the Caputo derivative of fractional order $0<\beta <1$ of the translates of the truncated power function

$${}^C{D}_x^{\beta }{(r-x)}_{+}^3=-\frac{3}{\Gamma (4-\beta )}\left((2x^2-2(3-\beta )xr+{(2-\beta )}_2{r}^2)x^{1-\beta }-2{(x-r)}_{+}^{3-\beta }\right).$$

(78)

Thus, substituting (78) in (48), we get

$${}^C{D}_x^{\beta }{N}_{03}\left(x\right)=\frac{3}{\Gamma (4-\beta )}\left(-\left(2x^2-2(3-\beta )x+{(2-\beta )}_2\right)x^{1-\beta }+2{(x-1)}_{+}^{3-\beta }\right),$$

(79)

$${}^C{D}_x^{\beta }{N}_{13}\left(x\right)=\frac{3}{\Gamma (4-\beta )}\left((\frac{7}{2}{x}^2-3(3-\beta )x+{(2-\beta )}_2)x^{1-\beta }-4{(x-1)}_{+}^{3-\beta }+\frac{1}{2}{(x-2)}_{+}^{3-\beta }\right),$$

(80)

$${}^C{D}_x^{\beta }{N}_{23}\left(x\right)=\frac{3}{\Gamma (4-\beta )}\left((-\frac{11}{6}{x}^2+(3-\beta )x)x^{1-\beta }+3{(x-1)}_{+}^{3-\beta }-\frac{3}{2}{(x-2)}_{+}^{3-\beta }+\frac{1}{3}{(x-3)}_{+}^{3-\beta }\right).$$

(81)

The optimal basis and its fractional derivatives in case of equidistant knots on an interval $[a,b]$ can be evaluated by scaling as shown in (27)–(29). The optimal B-spline basis in the interval $[0,2]$ for $h=0.25$ is shown in Figure 1a. The Caputo derivatives of fractional order $\beta =0.25,0.5,0.75$ are shown in Figure 1b–d.

To give an idea of the condition number of the discretization matrix $(Q\otimes A+L\otimes B)$ (cf. (68)), in

Table 1. The condition number $\mathcal{K}$ (rounded to the nearest integer) and the dimensions $n_{rows}\times n_{cols}$ of the discretization matrix $(Q\otimes A+L\otimes B)$ for $\alpha =0.25$ as a function of the space step h and the time step $\delta$.

	$\mathit{\delta }=0.25$		$\mathit{\delta }=0.125$		$\mathit{\delta }=0.0625$		$\mathit{\delta }=0.03125$
$\mathit{h}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$
0.25	143	36×36	125	72×54	119	144×90	110	288×162
0.125	543	68×68	472	136×102	450	272×170	415	544×306
0.0625	2144	132×132	1864	264×198	1776	528×330	1639	1056×594
0.03125	8550	260×260	7434	520×390	7083	1040×650	6535	2080×1170

,

Table 2. The condition number $\mathcal{K}$ (rounded to the nearest integer) and the dimensions $n_{rows}\times n_{cols}$ of the discretization matrix $(Q\otimes A+L\otimes B)$ for $\alpha =0.5$ as a function of the space step h and the time step $\delta$.

	$\mathit{\delta }=0.25$		$\mathit{\delta }=0.125$		$\mathit{\delta }=0.0625$		$\mathit{\delta }=0.03125$
$\mathit{h}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$
0.25	102	36×36	77	72×54	66	144×90	66	288×162
0.125	386	68×68	291	136×102	234	272×170	179	544×306
0.0625	1525	132×132	1148	264×198	923	528×330	706	1056×594
0.03125	6084	260×260	4579	520×390	3681	1040×650	2817	2080×1170

and

Table 3. The condition number $\mathcal{K}$ (rounded to the nearest integer) and the dimensions $n_{rows}\times n_{cols}$ of the discretization matrix $(Q\otimes A+L\otimes B)$ for $\alpha =0.75$ as a function of the space step h and the time step $\delta$.

	$\mathit{\delta }=0.25$		$\mathit{\delta }=0.125$		$\mathit{\delta }=0.0625$		$\mathit{\delta }=0.03125$
$\mathit{h}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$	$\mathcal{K}$	${\mathit{n}}_{\mathit{rows}}\times {\mathit{n}}_{\mathit{cols}}$
0.25	69	36×36	62	72×54	63	144×90	60	288×162
0.125	253	68×68	154	136×102	100	272×170	67	544×306
0.0625	1001	132×132	610	264×198	394	528×330	244	1056×594
0.03125	3991	260×260	2434	520×390	1573	1040×650	972	2080×1170

, we list the condition number of the matrix when using the cubic B-spline basis for different values of the order of the fractional derivative and different values of the space and time steps.

3.2. Numerical Tests

In this section, we show the numerical results we obtained when solving the differential problem (57) for three different expressions of the known term $f(x,t)$ and different values of the fractional derivative. In the tests, we set $\tau =1$ and $\eta =2$ and chose as collocation nodes the equidistant points of the interval $[0,1]$ with time step $\delta$. All the tests were performed in the Mathematica environment [32].

3.2.1. Test 1

To test the accuracy of the method, we first solved the time-fractional diffusion Equation (57) when

$$f(x,t)=\left(\frac{6}{\Gamma (4-\alpha )}{t}^{-\alpha }x(L-x)-2\right)t^3.$$

(82)

In this case, the exact solution is the bivariate polynomial

$$u(x,t)=x(L-x)t^3.$$

(83)

Since the cubic B-spline basis reproduces polynomials up to degree 3, the approximation error of the collocation–Galerkin method is zero so that we expect the numerical solution coincides with the exact solution.

We evaluated the numerical solution for $\delta =0.25$, $s=0.5$ and $h=0.25$ and different values of the order of the time fractional derivative $\beta$.

Let us define the error by

$$e_{h,\delta }(x,t)=u(x,t)-u_{h,\delta }(x,t),$$

(84)

where $u_{h,\delta }(x,t)$ denotes the numerical solution of the fractional differential system (63) obtained by the collocation method.

The tests show that, as expected, the error is on the order of the machine precision even using a large step size. The numerical solution and the error for $\beta =0.25$ are shown in Figure 2.

3.2.2. Test 2

To test the accuracy of the method when approximating time fractional derivatives, we solved the time-fractional diffusion Equation (57) for

$$f(x,t)=x(L-x)\frac{\pi t^{1-\beta }}{2\Gamma (2-\beta )}\left({}_1{F}_1(1,2-\beta ,i\pi t)+{}_1{F}_1(1,2-\beta ,-i\pi t)\right)+2sin\left(\pi t\right),$$

(85)

where ${}_1{F}_1(a,b,z)$ is the Kummer’s confluent hypergeometric function

$${}_1{F}_1(a,b,z)=\frac{\Gamma \left(b\right)}{\Gamma \left(a\right)}\sum _{k=0}^{\infty }\frac{\Gamma (a+k)}{\Gamma (b+k)k!}{z}^k,a\in \mathbb{R},-b\notin {\mathbb{N}}_0,$$

(86)

where ${\mathbb{N}}_0=\mathbb{N}\setminus \left\{0\right\}$ (cf. [9]). In this case, the exact solution is

$$u(x,t)=x(L-x)sin\left(\pi t\right)$$

(87)

so that the error for the spatial approximation is negligible.

The $L_2$-norm of the error

$$\parallel e_{h,\delta }{\parallel }_2=\sqrt{{\int }_0^{\tau }dt{\int }_0^{\eta }dx{\left|e_{h,\delta }(x,t)\right|}^2},$$

(88)

and the numerical convergence order

$${\rho }_{\delta }={log}_2\left(\frac{\parallel e_{h,\delta }{\parallel }_2}{\parallel e_{h,\delta /2}{\parallel }_2}\right),$$

(89)

are listed in

Table 4. Test 2: The $L_2$-norm of the error and the numerical convergence order ${\rho }_{\delta }$ as a function of the time step $\delta$ for different values of the order $\beta$ of the fractional derivative.

	$\mathit{\beta }=0.25$		$\mathit{\beta }=0.5$		$\mathit{\beta }=0.75$
$\mathit{\delta }$	$\parallel {\mathit{e}}_{\mathit{h},\mathit{\delta }}{\parallel }_{\mathbf{2}}$	${\mathit{\rho }}_{\mathit{\delta }}$	$\parallel {\mathit{e}}_{\mathit{h},\mathit{\delta }}{\parallel }_{\mathbf{2}}$	${\mathit{\rho }}_{\mathit{\delta }}$	$\parallel {\mathit{e}}_{\mathit{h},\mathit{\delta }}{\parallel }_{\mathbf{2}}$	${\mathit{\rho }}_{\mathit{\delta }}$
0.25	$0.42\times {10}^{-2}$		$0.50\times {10}^{-2}$		$0.62\times {10}^{-2}$
0.125	$0.32\times {10}^{-3}$	3.71	$0.32\times {10}^{-3}$	3.97	$0.34\times {10}^{-3}$	4.19
0.0625	$0.17\times {10}^{-4}$	4.23	$0.17\times {10}^{-4}$	4.23	$0.19\times {10}^{-4}$	4.16
0.03125	$0.10\times {10}^{-5}$	4.09	$0.11\times {10}^{-5}$	3.95	$0.12\times {10}^{-5}$	3.98

. The numerical solution in the case when $\beta =0.25$ obtained for $\delta =0.03125$, $s=2\delta$ and $h=0.25$ is shown in Figure 3.

The numerical results show that the method converges with convergence order approximatively equal to 4 and gives a good approximation even using few collocation points. Moreover, the accuracy of the numerical solution does not depend on the order of the fractional derivative.

3.2.3. Test 3

Finally, we solved the time-fractional diffusion equation (57) for

$$f(x,t)=sin\left(\pi x\right)\left(\frac{\pi t^{1-\beta }}{2\Gamma (2-\beta )}\left({}_1{F}_1(1,2-\beta ,i\pi t)+{}_1{F}_1(1,2-\beta ,-i\pi t)\right)+{\pi }^{2}sin\left(\pi t\right)\right).$$

(90)

In this case, the exact solution is

$$u(x,t)=sin\left(\pi x\right)sin\left(\pi t\right).$$

(91)

The $L_2$-norm of the error and the numerical convergence order

$${\rho }_{h,\delta }={log}_2\left(\frac{\parallel e_{h,\delta }{\parallel }_2}{\parallel e_{h/2,\delta /2}{\parallel }_2}\right)$$

(92)

are listed in

Table 5. Test 3: The $L_2$-norm of the error and the numerical convergence order ${\rho }_{h,\delta }$ as a function of the time step $\delta$ and the space step h for different values of the order $\beta$ of the fractional derivative.

	$\mathit{\beta }=0.25$		$\mathit{\beta }=0.5$		$\mathit{\beta }=0.75$
$\mathit{h}=\mathit{\delta }$	$\parallel {\mathit{e}}_{\mathit{h},\mathit{\delta }}{\parallel }_{\mathbf{2}}$	$\mathit{\rho }$	$\parallel {\mathit{e}}_{\mathit{h},\mathit{\delta }}{\parallel }_{\mathbf{2}}$	${\mathit{\rho }}_{\mathit{h},\mathit{\delta }}$	$\parallel {\mathit{e}}_{\mathit{h},\mathit{\delta }}{\parallel }_{\mathbf{2}}$	${\mathit{\rho }}_{\mathit{h},\mathit{\delta }}$
0.25	$0.38\times {10}^{-2}$		$0.41\times {10}^{-2}$		$0.46\times {10}^{-2}$
0.125	$0.31\times {10}^{-3}$	3.62	$0.31\times {10}^{-3}$	3.72	$0.31\times {10}^{-3}$	3.89
0.0625	$0.16\times {10}^{-4}$	4.28	$0.16\times {10}^{-4}$	4.28	$0.17\times {10}^{-4}$	4.19
0.03125	$0.98\times {10}^{-6}$	4.03	$0.99\times {10}^{-6}$	4.01	$0.10\times {10}^{-5}$	4.09

. The numerical solution in the case when $\beta =0.25$ obtained for $\delta =h=0.03125$ and $s=2\delta$ is shown in Figure 4.

Also in this case, the numerical solution has a good accuracy that is not affected by the order of the fractional derivative. The numerical convergence order is approximatively 4.

4. Discussion

In this paper, we showed how optimal B-spline bases have good approximation properties that make them particularly suitable to be used in the solution of fractional differential problems. First of all, boundary and initial conditions can be satisfied easily in view of properties (22)–(23) so that the numerical solution we obtained does not suffer from the numerical instabilities near the boundaries that appear when truncated functions at the boundaries are used (cf. [28]). Moreover, optimal B-spline bases are stable (cf. (25)) meaning that numerical errors are not amplified. Nevertheless, Table 1, Table 2 and Table 3 show that the condition number of the discretization matrix increases as the dimension of the matrix increases. The conditioning gets worse when the order of the fractional derivative becomes smaller. In this case, the linear system (68) can be accurately solved by Krylov methods [33].

As for the convergence order, we observe that several methods proposed in the literature to solve fractional differential problems have convergence order that depends on the order of the fractional derivative (see, for instance, [26,27]). In the case of the cubic B-spline basis, the numerical tests show that the numerical convergence order does not depend on the order of the fractional derivative being close to 4 for any value of $\beta$. This result is in line with the Strang–Fix theory for classical differential problems [34]. We expect that the convergence order for the optimal B-spline basis ${\mathcal{N}}_n$ is $n+1$ so that it can be increased easily by increasing the degree of the basis. In fact, optimal bases of high degree and their fractional derivatives can be evaluated by using the analytical expression given in Section 2.4 and Section 2.5. As far as we know, the analytical expressions (39) and (48) are new and, together with (27)–(29), are an easy tool to evaluate optimal B-spline bases on any set of equidistant knots on bounded intervals. Finally, we observe that the accuracy of the numerical solution we obtained is higher than the accuracy obtained in [26] where a quadrature formula was used to approximate the time fractional derivative and a finite element method was used for the spatial approximation.

The optimal bases we described can be used to solve other fractional differential problems, for instance problems involving the Riesz derivative [2], or can be used in other numerical methods, for instance in the operational matrix method or in wavelet methods. We notice that adaptive wavelet methods can be used to solve differential problems whose solutions are non smooth functions. Optimal bases generating multiresolution analyses and wavelet spaces can be obtained starting from a special family of refinable functions [35,36] so that the procedure described above can be generalized to other refinable bases.