---
abstract: |
    The aim of this paper is to numerically solve a diffusion differential problem having time derivative of fractional order. To this end we propose a collocation-Galerkin method that uses the fractional splines as approximating functions. The main advantage is in that the derivatives of integer and fractional order of the fractional splines can be expressed in a closed form that involves just the generalized finite difference operator. This allows us to construct an accurate and efficient numerical method. Several numerical tests showing the effectiveness of the proposed method are presented.\
    [**Keywords**]{}: Fractional diffusion problem, Collocation method, Galerkin method, Fractional spline
author:
- 'Laura Pezza[^1], Francesca Pitolli[^2]'
title: 'A fractional spline collocation-Galerkin method for the time-fractional diffusion equation'
---

Introduction. {#sec:intro}
=============

The use of fractional calculus to describe real-world phenomena is becoming increasingly widespread. Integro-differential equations of [*fractional*]{}, [*i.e.*]{} positive real, order are used, for instance, to model wave propagation in porous materials, diffusive phenomena in biological tissue, viscoelastic properties of continuous media [@Hi00; @Ma10; @KST06; @Ta10]. Among the various fields in which fractional models are successfully used, viscoelasticity is one of the more interesting since the memory effect introduced by the time-fractional derivative allows to model anomalous diffusion phenomena in materials that have mechanical properties in between pure elasticity and pure viscosity [@Ma10]. Even if these models are empirical, nevertheless they are shown to be consistent with experimental data.\
The increased interest in fractional models has led to the development of several numerical methods to solve fractional integro-differential equations. Many of the proposed methods generalize to the fractional case numerical methods commonly used for the classical integer case (see, for instance, [@Ba12; @PD14; @ZK14] and references therein). But the nonlocality of the fractional derivative raises the challenge of obtaining numerical solution with high accuracy at a low computational cost. In [@PP16] we proposed a collocation method especially designed for solving differential equations of fractional order in time. The key ingredient of the method is the use of the fractional splines introduced in [@UB00] as approximating functions. Thus, the method takes advantage of the explicit differentiation rule for fractional B-splines that allows us to evaluate accurately the derivatives of both integer and fractional order.\
In the present paper we used the method to solve a diffusion problem having time derivative of fractional order and show that the method is efficient and accurate. More precisely, the [*fractional spline collocation-Galerkin method*]{} here proposed combines the fractional spline collocation method introduced in [@PP16] for the time discretization and a classical spline Galerkin method in space.\
The paper is organized as follows. In Section \[sec:diffeq\], a time-fractional diffusion problem is presented and the definition of fractional derivative is given. Section \[sec:fractBspline\] is devoted to the fractional B-splines and the explicit expression of their fractional derivative is given. The fractional spline approximating space is described in Section \[sec:app\_spaces\], while the fractional spline collocation-Galerkin method is introduced in Section \[sec:Galerkin\]. Finally, in Section \[sec:numtest\] some numerical tests showing the performance of the method are displayed. Some conclusions are drawn in Section \[sec:concl\].

A time-fractional diffusion problem. {#sec:diffeq}
====================================

We consider the [*time-fractional differential diffusion problem*]{} [@Ma10] $$\label{eq:fracdiffeq}
\left \{ \begin{array}{lcc} 
\displaystyle D_t^\gamma \, u(t, x) - \frac{\partial^2}{\partial x^2} \, u(t, x) = f(t, x)\,, & \quad t \in [0, T]\,, &  \quad x \in [0,1] \,,\\ \\
u(0, x) = 0\,, &  & \quad x \in [0,1]\,, \\ \\
u(t, 0) = u(t, 1) = 0\,, & \quad t \in [0, T]\,,    
\end{array} \right.$$ where $ D_t^\gamma u$, $0 < \gamma < 1$, denotes the [*partial fractional derivative*]{} with respect to the time $t$. Usually, in viscoelasticity the fractional derivative is to be understood in the Caputo sense, [*i.e.*]{} $$\label{eq:Capfrac}
D_t^\gamma \, u(t, x) = \frac1{\Gamma(1-\gamma)} \,  \int_0^t \, \frac{u_t(\tau,x)}{(t - \tau)^\gamma} \, d\tau\,, \qquad t\ge 0\,,$$ where $\Gamma$ is the Euler’s gamma function $$\Gamma(\gamma+1)= \int_0^\infty \, s^\gamma \, {\rm e}^{-s} \, ds\,.$$ We notice that due to the homogeneous initial condition for the function $u(t,x)$, solution of the differential problem (\[eq:fracdiffeq\]), the Caputo definition (\[eq:Capfrac\]) coincides with the Riemann-Liouville definition (see [@Po99] for details). One of the advantage of the Riemann-Liouville definition is in that the usual differentiation operator in the Fourier domain can be easily extended to the fractional case, [*i.e.*]{} $${\cal F} \bigl(D_t^\gamma \, f(t) \bigr) = (i\omega)^\gamma {\cal F} (f(t))\,,$$ where ${\cal F}(f)$ denotes the Fourier transform of the function $f(t)$. Thus, analytical Fourier methods usually used in the classical integer case can be extended to the fractional case [@Ma10].

The fractional B-splines and their fractional derivatives. {#sec:fractBspline}
==========================================================

The [*fractional B-splines*]{}, [*i.e.*]{} the B-splines of fractional degree, were introduced in [@UB00] generalizing to the fractional power the classical definition for the polynomial B-splines of integer degree. Thus, the fractional B-spline $B_{\alpha}$ of degree $\alpha$ is defined as $$\label{eq:Balpha}
B_{\alpha}(t) := \frac{{ \Delta}^{\alpha+1} \, t_+^\alpha} {\Gamma(\alpha+1)}\,, \qquad \alpha > -\frac 12\,,$$ where $$\label{eq:fracttruncpow}
t_+^\alpha: = \left \{ \begin{array}{ll}
t^\alpha\,, & \qquad t \ge 0\,, \\ \\
0\,, & \qquad \hbox{otherwise}\,,
\end{array} \right. \qquad \alpha > -1/2\,,$$ is the [*fractional truncated power function*]{}. $\Delta^{\alpha}$ is the [*generalized finite difference operator*]{} $$\label{eq:fracfinitediff}
\Delta^{\alpha} \, f(t) :=  \sum_{k\in \NN} \, (-1)^k \, {\alpha \choose k} \, f(t-\,k)\,, 
\qquad \alpha \in \RR^+\,,$$ where $$\label{eq:binomfrac}
{\alpha \choose k} := \frac{\Gamma(\alpha+1)}{k!\, \Gamma(\alpha-k+1)}\,, \qquad k\in \NN\,, \quad \alpha \in \RR^+\,,$$ are the [*generalized binomial coefficients*]{}. We notice that ’fractional’ actually means ’noninteger’, [*i.e.*]{} $\alpha$ can assume any real value greater than $-1/2$. For real values of $\alpha$, $B_\alpha$ does not have compact support even if it belongs to $L_2(\RR)$. When $\alpha=n$ is a nonnegative integer, Equations (\[eq:Balpha\])-(\[eq:binomfrac\]) are still valid; $\Delta^{n}$ is the usual finite difference operator so that $B_n$ is the classical polynomial B-spline of degree $n$ and compact support $[0,n+1]$ (for details on polynomial B-splines see, for instance, the monograph [@Sc07]).

The fractional B-splines for different values of the parameter $\alpha$ are displayed in Figure \[fig:fractBsplines\] (top left panel). The classical polynomial B-splines are also displayed (dashed lines). The picture shows that the fractional B-splines decay very fast toward infinity so that they can be assumed compactly supported for computational purposes. Moreover, in contrast to the polynomial B-splines, fractional splines are not always positive even if the nonnegative part becomes more and more smaller as $\alpha$ increases.

  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
        ![Top left panel: The fractional B-splines (solid lines) and the polynomial B-splines (dashed lines) for $\alpha$ ranging from 0 to 4. Top right panel: The fractional derivatives of the linear B-spline $B_1$ for $\gamma = 0.25, 0.5, 0.75$. Bottom left panel: The fractional derivatives of the cubic B-spline $B_3$ for $\gamma$ ranging from 0.25 to 2. Bottom right panel: The fractional derivatives of the fractional B-spline $B_{3.5}$ for the $\gamma$ ranging from 0.25 to 2. Ordinary derivatives are displayed as dashed lines. []{data-label="fig:fractBsplines"}](Fig_fract_Bspline.png "fig:"){width="6cm"}         ![Top left panel: The fractional B-splines (solid lines) and the polynomial B-splines (dashed lines) for $\alpha$ ranging from 0 to 4. Top right panel: The fractional derivatives of the linear B-spline $B_1$ for $\gamma = 0.25, 0.5, 0.75$. Bottom left panel: The fractional derivatives of the cubic B-spline $B_3$ for $\gamma$ ranging from 0.25 to 2. Bottom right panel: The fractional derivatives of the fractional B-spline $B_{3.5}$ for the $\gamma$ ranging from 0.25 to 2. Ordinary derivatives are displayed as dashed lines. []{data-label="fig:fractBsplines"}](Fig_fractder_Bspline_linear.png "fig:"){width="6cm"}
   ![Top left panel: The fractional B-splines (solid lines) and the polynomial B-splines (dashed lines) for $\alpha$ ranging from 0 to 4. Top right panel: The fractional derivatives of the linear B-spline $B_1$ for $\gamma = 0.25, 0.5, 0.75$. Bottom left panel: The fractional derivatives of the cubic B-spline $B_3$ for $\gamma$ ranging from 0.25 to 2. Bottom right panel: The fractional derivatives of the fractional B-spline $B_{3.5}$ for the $\gamma$ ranging from 0.25 to 2. Ordinary derivatives are displayed as dashed lines. []{data-label="fig:fractBsplines"}](Fig_fractder_Bspline_cubica.png "fig:"){width="6cm"}   ![Top left panel: The fractional B-splines (solid lines) and the polynomial B-splines (dashed lines) for $\alpha$ ranging from 0 to 4. Top right panel: The fractional derivatives of the linear B-spline $B_1$ for $\gamma = 0.25, 0.5, 0.75$. Bottom left panel: The fractional derivatives of the cubic B-spline $B_3$ for $\gamma$ ranging from 0.25 to 2. Bottom right panel: The fractional derivatives of the fractional B-spline $B_{3.5}$ for the $\gamma$ ranging from 0.25 to 2. Ordinary derivatives are displayed as dashed lines. []{data-label="fig:fractBsplines"}](Fig_fractder_Bspline_alpha3p5.png "fig:"){width="6cm"}
  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

The fractional derivatives of the fractional B-splines can be evaluated explicitly by differentiating (\[eq:Balpha\]) and (\[eq:fracttruncpow\]) in the Caputo sense. This gives the following differentiation rule $$\label{eq:diffrule_tronc}
D^{\gamma}_t \, B_{\alpha} (t)=  \frac{\Delta^{\alpha+1} \, t_+^{\alpha-\gamma}} {\Gamma(\alpha-\gamma+1)}\,, \qquad 0 < \gamma < \alpha + \frac12\,,$$ which holds both for fractional and integer order $\gamma$. In particular, when $\gamma, \alpha$ are nonnegative integers, (\[eq:diffrule\_tronc\]) is the usual differentiation rule for the classical polynomial B-splines [@Sc07]. We observe that since $B_\alpha$ is a causal function with $B_\alpha^{(n)}(0)=0$ for $n\in \NN\backslash\{0\}$, the Caputo fractional derivative coincides with the Riemann-Liouville fractional derivative.\
From (\[eq:diffrule\_tronc\]) and the composition property $\Delta^{\alpha_1} \, \Delta^{\alpha_2} = \Delta^{\alpha_1+\alpha_2}$ it follows [@UB00] $$\label{eq:diffrule_2}
D^{\gamma}_t \, B_{\alpha} = \Delta ^{\gamma} \, B_{\alpha-\gamma}\,,$$ [*i.e.*]{} the fractional derivative of a fractional B-spline of degree $\alpha$ is a fractional spline of degree $\alpha-\gamma$. The fractional derivatives of the classical polynomial B-splines $B_n$ are fractional splines, too. This means that $D^{\gamma}_t \, B_{n}$ is not compactly supported when $\gamma$ is noninteger reflecting the nonlocal behavior of the derivative operator of fractional order.\
In Figure \[fig:fractBsplines\] the fractional derivatives of $B_1$ (top right panel), $B_3$ (bottom left panel) and $B_{3.5}$ (bottom right panel) are displayed for different values of $\gamma$.

The fractional spline approximating spaces. {#sec:app_spaces}
===========================================

A property of the fractional B-splines that is useful for the construction of numerical methods for the solution of differential problems is the [*refinability*]{}. In fact, the fractional B-splines are [*refinable functions*]{}, [*i.e.*]{} they satisfy the [*refinement equation*]{} $$B_\alpha(t) = \sum_{k\in \NN} \, a^{(\alpha)}_{k} \, B_\alpha(2\,t-k)\,, \qquad t \ge 0\,,$$ where the coefficients $$a^{(\alpha)}_{k} := \frac{1}{2^{\alpha}} {\alpha+1 \choose k}\,,\qquad k\in \NN\,,$$ are the [*mask coefficients*]{}. This means that the sequence of nested approximating spaces $$V^{(\alpha)}_j(\RR) = {\rm span} \,\{B_\alpha(2^j\, t -k), k \in \ZZ\}\,, \qquad j \in \ZZ\,,$$ forms a [*multiresolution analysis*]{} of $L_2(\RR)$.

As a consequence, any function $f_j(t)$ belonging to $V^{(\alpha)}_j(\RR)$ can be expressed as $$f_j(t) = \sum_{k\in \ZZ}\, \lambda_{jk} \, B_\alpha(2^j\, t -k)\,,$$ where the coefficient sequence $\{\lambda_{j,k}\} \in \ell_2(\ZZ)$. Moreover, any space $V^{(\alpha)}_j(\RR)$ reproduces polynomials up to degree $\lceil \alpha\rceil$, [*i.e.*]{} $x^d \in V^{(\alpha)}_j(\RR)$, $ 0 \le d \le \lceil \alpha\rceil$, while its approximation order is $\alpha +1$. We recall that the polynomial B-spline $B_n$ reproduces polynomial up to degree $n$ whit approximation order $n+1$ [@UB00].

To solve boundary differential problems we need to construct a multiresolution analysis on a finite interval. For the sake of simplicity in the following we will consider the interval $I=[0,1]$. A simple approach is to restrict the basis $\{B_\alpha(2^j\, t -k)\}$ to the interval $I$, [*i.e.*]{} $$\label{eq:Vj_int}
V^{(\alpha)}_j(I) = {\rm span} \,\{B_\alpha(2^j\, t -k), t\in I, -N \le k \le 2^j-1\}\,, \qquad j_0 \le j\,,$$ where $N$ is a suitable index, chosen in order the significant part of $B_\alpha$ is contained in $[0,N+1]$, and $j_0$ is the starting refinement level. The drawback of this approach is its numerical instability and the difficulty in fulfilling the boundary conditions since there are $2N$ boundary functions, [*i.e.*]{} the translates of $B_\alpha$ having indexes $ -N\le k \le -1$ and $2^j-N\le k \le 2^j-1$, that are non zero at the boundaries. More suitable refinable bases can be obtained by the procedure given in [@GPP04; @GP04]. In particular, for the polynomial B-spline $B_n$ a B-basis $\{\phi_{\alpha,j,k}(t)\}$ with optimal approximation properties can be constructed. The internal functions $\phi_{\alpha,j,k}(t)=B_\alpha(2^j\, t -k)$, $0 \le k \le 2^j-1-n$, remain unchanged while the $2n$ boundary functions fulfill the boundary conditions $$\begin{array}{llcc}
\phi_{\alpha,j,-1}^{(\nu)}(0) = 1\,, & \phi_{\alpha,j,k}^{(\nu)}(0) = 0\,, &\hbox{for} & 0\le \nu \le -k-2\,, -n \le k \le -2\,,\\ \\
\phi_{\alpha,j,2^j-1}^{(\nu)}(1) = 1\,, & \phi_{\alpha,j,2^j+k}^{(\nu)}(1) = 0\,, & \hbox{for} & 0\le \nu \le -k-2\,, -n \le k \le -2\,,\\ \\
\end{array}$$ Thus, the B-basis naturally fulfills Dirichlet boundary conditions.\
As we will show in the next section, the refinability of the fractional spline bases plays a crucial role in the construction of the collocation-Galerkin method.

The fractional spline collocation-Galerkin method. {#sec:Galerkin}
==================================================

In the collocation-Galerkin method here proposed, we look for an approximating function $u_{s,j}(t,x) \in V^{(\beta)}_s([0,T]) \otimes  V^{(\alpha)}_j([0,1])$. Since just the ordinary first spatial derivative of $u_{s,j}$ is involved in the Galerkin method, we can assume $\alpha$ integer and use as basis function for the space $V^{(\alpha)}_j([0,1])$ the refinable B-basis $\{\phi_{\alpha,j,k}\}$, [*i.e.*]{} $$\label{uj}
u_{s,j}(t,x) = \sum_{k \in {\cal Z}_j} \, c_{s,j,k}(t) \, \phi_{\alpha,j,k}(x)\,,$$ where the unknown coefficients $c_{s,j,k}(t)$ belong to $V^{(\beta)}_s([0,T])$. Here, ${\cal Z}_j$ denotes the set of indexes $-n\le k \le 2^j-1$.\
The approximating function $u_{s,j}(t,x)$ solves the variational problem $$\label{varform}
\left \{ 
\begin{array}{ll} 
\displaystyle \left ( D_t^\gamma u_{s,j},\phi_{\alpha,j,k} \right ) -\left ( \frac {\partial^2} {\partial x^2}\,u_{s,j},\phi_{\alpha,j,k} \right ) = \left ( f,\phi_{\alpha,j,k} \right )\,, & \quad k \in {\cal Z}_j\,, \\ \\
u_{s,j}(0, x) = 0\,, & x \in [0,1]\,, \\ \\
u_{s,j}(t, 0) = 0\,, \quad u_{s,j}(t,1) = 0\,, & t \in [0,T]\,,
\end{array} \right.$$ where $(f,g)= \int_0^1 \, f\,g$.\
Now, writing (\[varform\]) in a weak form and using (\[uj\]) we get the system of fractional ordinary differential equations $$\label{fracODE}
\left \{ \begin{array}{ll}
M_j \, D_t^\gamma\,C_{s,j}(t) + L_j\, C_{s,j}(t) = F_j(t)\,, & \qquad t \in [0,T]\,, \\ \\
C_{s,j}(0) = 0\,,
\end{array} \right.$$ where $C_{s,j}(t)=(c_{s,j,k}(t))_{k\in {\cal Z}_j}$ is the unknown vector. The connecting coefficients, i.e. the entries of the mass matrix $M_j = (m_{j,k,i})_{k,i\in{\cal Z}_j}$, of the stiffness matrix $L_j = (\ell_{j,k,i})_{k,i\in{\cal Z}_j}$, and of the load vector $F_j(t)=(f_{j,k}(t))_{k\in {\cal Z}_j}$, are given by $$m_{j,k,i} = \int_0^1\, \phi_{\alpha,j,k}\, \phi_{\alpha,j,i}\,, \qquad \ell_{j,k,i} = \int_0^1 \, \phi'_{\alpha,j,k} \, \phi'_{\alpha,j,i}\,,$$ $$f_{j,k}(t) = \int_0^1\, f(t,\cdot)\, \phi_{\alpha,j,k}\,.$$ The entries of $M_j$ and $L_j$ can be evaluated explicitly using (\[eq:Balpha\]) and (\[eq:diffrule\_tronc\]), respectively, while the entries of $F_j(t)$ can be evaluated by quadrature formulas especially designed for wavelet methods [@CMP15; @GGP00].

To solve the fractional differential system (\[fracODE\]) we use the collocation method introduced in [@PP16]. For an integer value of $T$, let $t_p =  p/2^q$, $0\le p \le 2^q\,T$, where $q$ is a given nonnegative integer, be a set of dyadic nodes in the interval $[0,T]$. Now, assuming $$\label{ck}
c_{s,j,k}(t) = \sum_{r\in {\cal R}_s} \, \lambda_{k,r}\,\chi_{\beta,s,r}(t) \,, \qquad k \in {\cal Z}_j\,,$$ where $\chi_{\beta,s,r}(t)=B_\beta(2^s\,t-r)$ with $B_\beta$ a fractional B-spline of fractional degree $\beta$, and collocating (\[fracODE\]) on the nodes $t_p$, we get the linear system $$\label{colllinearsys}
(M_j\otimes A_s + L_j\otimes G_s) \,\Lambda_{s,j}  =F_j\,,$$ where $\Lambda_{s,j}=(\lambda_{k,r})_{r\in {\cal R}_s,k\in {\cal Z}_j}$ is the unknown vector, $$\begin{array}{ll}
A_s= \bigl( a_{p,r} \bigr)_{p\in {\cal P}_q,r\in {\cal R}_s}\,, & \qquad a_{p,r} = D_t^\gamma \, \chi_{\beta,s,r}(t_p)\,, \\ \\
G_s=\bigl(g_{p,r}\bigr)_{p \in {\cal P}_q,r\in {\cal R}_s}\,, & \qquad g_{p,r} = \chi_{\beta,s,r}(t_p)\,, 
\end{array}$$ are the collocation matrices and $$F_j=(f_{j,k}(t_p))_{k\in{\cal Z}_j,p \in {\cal P}_q}\,,$$ is the constant term. Here, ${\cal R}_s$ denotes the set of indexes $-\infty < r \le 2^s-1$ and ${\cal P}_q$ denotes the set of indexes $0<p\le 2^qT$. Since the fractional B-splines have fast decay, the series (\[ck\]) is well approximated by only few terms and the linear system (\[colllinearsys\]) has in practice finite dimension so that the unknown vector $\Lambda_{s,j}$ can be recovered by solving (\[colllinearsys\]) in the least squares sense.\
We notice that the entries of $G_s$, which involve just the values of $\chi_{\beta,s,r}$ on the dyadic nodes $t_p$, can be evaluated explicitly by (\[eq:Balpha\]). On the other hand, we must pay a special attention to the evaluation of the entries of $A_s$ since they involve the values of the fractional derivative $D_t^\gamma\chi_{\beta,s,r}(t_p)$. As shown in Section \[sec:fractBspline\], they can be evaluated efficiently by the differentiation rule (\[eq:diffrule\_2\]).

In the following theorem we prove that the fractional spline collocation-Galerkin method is convergent. First of all, let us introduce the Sobolev space on bounded interval $$H^\mu(I):= \{v \in L^2(I): \exists \, \tilde v \in H^\mu (\RR) \ \hbox{\rm such that} \ \tilde v|_I=v\}, \quad \mu\geq 0\,,$$ equipped with the norm $$\|v\|_{\mu,I}  = \inf_{\tilde v \in H^\mu(\RR), \tilde v|_I=v} \|\tilde v\|_{\mu,\RR}\,,$$ where $$H^\mu(\RR):= \{v: v\in L^2(\RR) \mbox{ and } (1+|\omega|^2)^{\mu/2}  {\cal F}(v)(\omega) \in L^2(\RR)\}, \quad \mu\geq 0\,,$$ is the usual Sobolev space with the norm $$\| v \| _{\mu,\RR} =\bigl \| (1+|\omega|^2)^{\mu/2}  {\cal F}(v)(\omega) \bigr \| _{0,\RR}\,.$$

\[Convergence\] Let $$H^\mu(I;H^{\tilde \mu}(\Omega)):= \{v(t,x): \| v(t,\cdot)\|_{H^{\tilde \mu}(\Omega)} \in H^\mu(I)\}, \quad \mu, {\tilde \mu} \geq 0\,,$$ equipped with the norm $$\|v\|_{H^\mu(I;H^{\tilde \mu}(\Omega))} := \bigl \| \|v(t,\cdot)\|_{H^{\tilde \mu}(\Omega)} \bigr\|_{\mu,I}\,.$$ Assume $u$ and $f$ in (\[eq:fracdiffeq\]) belong to $H^{\mu}([0,T];H^{\tilde \mu}([0,1]))$, $0\le \mu$, $0\le \tilde \mu$, and $H^{\mu-\gamma}([0,T];$ $H^{{\tilde \mu}-2}([0,1]))$, $0\le \mu-\gamma$, $0\le \tilde \mu-2$, respectively. Then, the fractional spline collocation-Galerkin method is convergent, [*i.e.*]{}, $$\|u-u_{s,j}\|_{H^0([0,T];H^0([0,1]))} \, \to 0  \quad \hbox{as} \quad s,j \to \infty\,.$$ Moreover, for $\gamma \le \mu \le \beta+1$ and $1 \le \tilde \mu \le \alpha +1$ the following error estimate holds: $$\begin{array}{lcl}
\| u-u_{s,j}\|_{H^0([0,T];H^0([0,1]))} &\leq & \left (\eta_1 \, 2^{-j\tilde \mu}  + \eta_2  \, 2^{-s\mu} \right ) \| u\|_{H^\mu([0,T];H^{\tilde \mu}([0,1]))}\,, 
\end{array}$$ where $\eta_1$ and $\eta_2$ are two constants independent of $s$ and $j$.

Let $u_j$ be the exact solution of the variational problem (\[varform\]). Following a classical line of reasoning (cf. [@Th06; @FXY11; @DPS94]) we get $$\begin{array}{l}
\| u-u_{j,s}\|_{H^0([0,T];H^0([0,1]))} \leq \\ \\ 
\rule{2cm}{0cm} \leq \|u-u_{j}\|_{H^0([0,T];H^0([0,1]))} + \| u_j-u_{j,s}\|_{H^0([0,T];H^0([0,1]))}\, \leq 
\\ \\ \rule{2cm}{0cm} \leq \eta_1  \, 2^{-j\tilde \mu}\, \| u\|_{H^0([0,T];H^{\tilde \mu}([0,1]))} + \eta_2 \, 2^{-s\mu} \, \| u\|_{H^\mu([0,T];H^0([0,1]))} \leq 
\\ \\ \rule{2cm}{0cm} \leq  \left ( \eta_1  \, 2^{-j\tilde \mu} + \eta_2 \, 2^{-s\mu} \right ) \, \|u\|_{H^\mu([0,T];H^{\tilde \mu}([0,1]))}\,.
\end{array}$$

Numerical tests. {#sec:numtest}
================

To shown the effectiveness of the fractional spline collocation-Galerkin method we solved the fractional diffusion problem (\[eq:fracdiffeq\]) for two different known terms $f(t,x)$ taken from [@FXY11]. In all the numerical tests we used as approximating space for the Galerkin method the (polynomial) cubic spline space. The B-splines $B_3$, its first derivatives $B_3'$ and the B-basis $\{\phi_{3,3,k}\}$ are displayed in Figure \[fig:Bcubic\]. We notice that since the cubic B-spline is centrally symmetric in the interval $[0,4]$, the B-basis is centrally symmetric, too. All the numerical tests were performed on a laptop using a Python environment. Each test takes a few minutes.

  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Left panel: The cubic B-spline (red line) and its first derivative (blue line). Right panel: The B-basis $\{\phi_{3,3,k}(x)\}$.[]{data-label="fig:Bcubic"}](Fig_Bspline_n3.png "fig:"){width="45.00000%"}   ![Left panel: The cubic B-spline (red line) and its first derivative (blue line). Right panel: The B-basis $\{\phi_{3,3,k}(x)\}$.[]{data-label="fig:Bcubic"}](OptBasis_alpha3.png "fig:"){width="45.00000%"}
  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Example 1
---------

In the first test we solved the time-fractional diffusion equation (\[eq:fracdiffeq\]) in the case when $$f(t,x)=\frac{2}{\Gamma(3-\gamma)}\,t^{2-\gamma}\, \sin(2\pi x)+4\pi^2\,t^2\, \sin(2\pi x)\,.$$ The exact solution is $$u(t,x)=t^2\,\sin(2\pi x).$$ We used the fractional B-spline $B_{3.5}$ as approximating function for the collocation method and solved the problem for $\gamma = 1, 0.75, 0.5, 0.25$. The fractional B-spline $B_{3.5}$, its first derivative and its fractional derivatives are shown in Figure \[fig:fract\_Basis\] along with the fractional basis $\{\chi_{3.5,3,r}\}$. The numerical solution $u_{s,j}(t,x)$ and the error $e_{s,j}(t,x) = u(t,x)-u_{s,j}(t,x)$ for $s=6$ and $j=6$ are displayed in Figure \[fig:numsol\_1\] for $\gamma = 0.5$. In all the numerical tests we set $q = s+1$.

  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Left panel: The fractional B-spline $B_{3.5}$ (green line), its first derivative (red line) and its fractional derivatives of order $\gamma =0.75$ (blue line), 0.5 (cyan line), 0.25 (black line). Right panel: The fractional basis $\{\chi_{3.5,3,r}\}$(right).[]{data-label="fig:fract_Basis"}](Fig_Bspline_n3p5.png "fig:"){width="45.00000%"}   ![Left panel: The fractional B-spline $B_{3.5}$ (green line), its first derivative (red line) and its fractional derivatives of order $\gamma =0.75$ (blue line), 0.5 (cyan line), 0.25 (black line). Right panel: The fractional basis $\{\chi_{3.5,3,r}\}$(right).[]{data-label="fig:fract_Basis"}](FractBasis_alpha3p5.png "fig:"){width="45.00000%"}
  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Example 1. The numerical solution (left panel) and the error (right panel) for $j=6$ and $s=6$ when $\gamma = 0.5$.[]{data-label="fig:numsol_1"}](Fig_NumSol_jGK6_sref7_beta3p5_gamma0p5_ex1.png "fig:"){width="45.00000%"}   ![Example 1. The numerical solution (left panel) and the error (right panel) for $j=6$ and $s=6$ when $\gamma = 0.5$.[]{data-label="fig:numsol_1"}](Fig_Error_jGK6_sref7_beta3p5_gamma0p5_ex1.png "fig:"){width="45.00000%"}
  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

We analyze the behavior of the error as the degree of the fractional B-spline $B_\beta$ increases. Figure \[fig:L2\_error\_1\] shows the $L_2$-norm of the error as a function of $s$ for $\beta$ ranging from 2 to 4; the four panels in the figure refer to different values of the order of the fractional derivative. For these tests we set $j=5$. The figure shows that for $s \le 4$ the error provided by the polynomial spline approximations is lower than the error provided by the fractional spline approximations. Nevertheless, in this latter case the error decreases reaching the same value, or even a lower one, of the polynomial spline error when $s=5$. We notice that for $\gamma=1$ the errors provided by the polynomial spline approximations of different degrees have approximatively the same values while the error provided by the polynomial spline of degree 2 is lower in case of fractional derivatives. In fact, it is well-known that fractional derivatives are better approximated by less smooth functions [@Po99].

  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
    ![Example 1: The $L_2$-norm of the error as a function of $q$ for different values of $\gamma$. Each line corresponds to a spline of different degree: solid lines correspond to the polynomial splines; non solid lines correspond to fractional splines.[]{data-label="fig:L2_error_1"}](Fig_L2Error_gam1_ex1.png "fig:"){width="45.00000%"}    ![Example 1: The $L_2$-norm of the error as a function of $q$ for different values of $\gamma$. Each line corresponds to a spline of different degree: solid lines correspond to the polynomial splines; non solid lines correspond to fractional splines.[]{data-label="fig:L2_error_1"}](Fig_L2Error_gam0p75_ex1.png "fig:"){width="45.00000%"}
   ![Example 1: The $L_2$-norm of the error as a function of $q$ for different values of $\gamma$. Each line corresponds to a spline of different degree: solid lines correspond to the polynomial splines; non solid lines correspond to fractional splines.[]{data-label="fig:L2_error_1"}](Fig_L2Error_gam0p5_ex1.png "fig:"){width="45.00000%"}   ![Example 1: The $L_2$-norm of the error as a function of $q$ for different values of $\gamma$. Each line corresponds to a spline of different degree: solid lines correspond to the polynomial splines; non solid lines correspond to fractional splines.[]{data-label="fig:L2_error_1"}](Fig_L2Error_gam0p25_ex1.png "fig:"){width="45.00000%"}
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Then, we analyze the convergence of the method for increasing values of $j$ and $s$. Table \[tab:conv\_js\_fract\_1\] reports the $L_2$-norm of the error for different values of $j$ and $s$ when using the fractional B-spline $B_{3.5}$ and $\gamma = 0.5$. The number of degrees-of-freedom is also reported. The table shows that the error decreases when $j$ increases and $s$ is held fix. We notice that the error decreases very slightly when $j$ is held fix and $s$ increases since for these values of $s$ we reached the accuracy level we can expect for that value of $j$ (cf. Figures \[fig:L2\_error\_1\]). The higher values of the error for $s=7$ and $j=5,6$ are due to the numerical instabilities of the basis $\{\chi_{3.5,s,r}\}$ which result in a high condition number of the discretization matrix. The error has a similar behavior even in the case when we used the cubic B-spline space as approximating space for the collocation method (cf. Table \[tab:conv\_js\_cubic\_1\]).

                                          3                4                5                6
  -------------------------------- ---------------- ---------------- ---------------- ----------------
   $\sharp V_j^{(\alpha)}([0,1])$         9                17               33               65
                 5                  0.02037 (369)    0.00449 (697)    0.00101 (1353)   0.00025 (2665)
                 6                  0.02067 (657)    0.00417 (1241)   0.00093 (2409)   0.00024 (4745)
                 7                  0.01946 (1233)   0.00381 (2329)   0.00115 (4521)   0.00117 (8905)

  : Example 1: The $L_2$-norm of the error for increasing values of $s$ and $j$ when using the fractional B-spline of degree $\beta=3.5$. The numbers in parenthesis are the degrees-of-freedom. Here, $\gamma =0.5$.

\[tab:conv\_js\_fract\_1\]

                                          3                4                5                6
  -------------------------------- ---------------- ---------------- ---------------- ----------------
   $\sharp V_j^{(\alpha)}([0,1])$         9                17               33               65
                 5                  0.02121 (315)    0.00452 (595)    0.00104 (1155)   0.00025 (2275)
                 6                  0.02109 (603)    0.00443 (1139)   0.00097 (2211)   0.00023 (4355)
                 7                  0.02037 (1179)   0.00399 (2227)   0.00115 (4323)   0.00115 (8515)

  : Example 1: The $L_2$-norm of the error for increasing values of $s$ and $j$ when using the cubic B-spline. The numbers in parenthesis are the degrees-of-freedom. Here, $\gamma =0.5$.

\[tab:conv\_js\_cubic\_1\]

Example 2
---------

In the second test we solved the time-fractional diffusion equation (\[eq:fracdiffeq\]) in the case when $$\begin{array}{lcl}
f(t,x) & = &  \displaystyle \frac{\pi t^{1-\gamma}}{2\Gamma(2-\gamma)} \left( \, {_1F_1}(1,2-\gamma,i\pi\,t) + \,{_1F_1}(1,2-\gamma,-i\pi\,t) \right) \, \sin(\pi\,x)  \\ \\
& + & \pi^2 \, \sin(\pi\,t)  \, \sin(\pi\,x)\,,
\end{array}$$ where $_1F_1(\alpha,\beta,z)$ is the Kummer’s confluent hypergeometric function defined as $$_1F_1(\alpha,\beta, z) =
\frac {\Gamma(\beta)}{\Gamma(\alpha)} \, \sum_{k\in \NN} \, \frac {\Gamma(\alpha+k)}{\Gamma(\beta+k)\, k!} \, z^k\,, \qquad \alpha \in \RR\,, \quad -\beta \notin \NN_0\,,$$ where $\NN_0 = \NN \backslash \{0\}$ (cf. [@AS65 Chapter 13]). In this case the exact solution is $$u(t,x)=\sin(\pi t)\,\sin(\pi x).$$ We performed the same set of numerical tests as in Example 1. The numerical solution $u_{s,j}(t,x)$ and the error $e_{s,j}(t,x)$ for $s=5$ and $j=6$ are displayed in Figure \[fig:numsol\_2\] in the case when $\gamma = 0.5$. Figure \[fig:L2\_error\_2\] shows the $L_2$-norm of the error as a function of $s$ for $\beta$ ranging from 2 to 4 and $j=5$; the four panels in the figure refer to different values of the order of the fractional derivative. Tables \[tab:conv\_js\_fract\_2\]-\[tab:conv\_js\_cubic\_2\] report the $L_2$-norm of the error for different values of $j$ and $s$ and $\beta = 3.5, 3$, respectively. The number of degrees-of-freedom is also reported.\
Figure \[fig:L2\_error\_2\] shows that value of the error is higher than in the previous example but it decreases as $s$ increases showing a very similar behavior as that one in Example 1. The values of the error in Tables \[tab:conv\_js\_fract\_2\]-\[tab:conv\_js\_cubic\_2\] are approximatively the same as in Tables \[tab:conv\_js\_fract\_1\]-\[tab:conv\_js\_cubic\_1\].

  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Example 2: The numerical solution (left panel) and the error (right panel) when $j=6$ and $s=5$. []{data-label="fig:numsol_2"}](Fig_NumSol_jGK6_sref6_beta3p5_gamma0p5_ex2.png "fig:"){width="45.00000%"}   ![Example 2: The numerical solution (left panel) and the error (right panel) when $j=6$ and $s=5$. []{data-label="fig:numsol_2"}](Fig_Error_jGK6_sref6_beta3p5_gamma0p5_ex2.png "fig:"){width="45.00000%"}
  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
    ![Example 2: The $L_2$-norm of the error as a function of $s$ for different values of $\gamma$. Each line corresponds to a spline of different degree: solid lines correspond to the polynomial splines; non solid lines correspond to fractional splines.[]{data-label="fig:L2_error_2"}](Fig_L2Error_gam1_ex2.png "fig:"){width="45.00000%"}    ![Example 2: The $L_2$-norm of the error as a function of $s$ for different values of $\gamma$. Each line corresponds to a spline of different degree: solid lines correspond to the polynomial splines; non solid lines correspond to fractional splines.[]{data-label="fig:L2_error_2"}](Fig_L2Error_gam0p75_ex2.png "fig:"){width="45.00000%"}
   ![Example 2: The $L_2$-norm of the error as a function of $s$ for different values of $\gamma$. Each line corresponds to a spline of different degree: solid lines correspond to the polynomial splines; non solid lines correspond to fractional splines.[]{data-label="fig:L2_error_2"}](Fig_L2Error_gam0p5_ex2.png "fig:"){width="45.00000%"}   ![Example 2: The $L_2$-norm of the error as a function of $s$ for different values of $\gamma$. Each line corresponds to a spline of different degree: solid lines correspond to the polynomial splines; non solid lines correspond to fractional splines.[]{data-label="fig:L2_error_2"}](Fig_L2Error_gam0p25_ex2.png "fig:"){width="45.00000%"}
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

                                          3                4                5                6
  -------------------------------- ---------------- ---------------- ---------------- ----------------
   $\sharp V_j^{(\alpha)}([0,1])$         9                17               33               65
                 5                  0.01938 (369)    0.00429 (697)    0.00111 (1353)   0.00042 (2665)
                 6                  0.01809 (657)    0.00555 (1241)   0.00507 (2409)   0.00523 (4745)
                 7                  0.01811 (1233)   0.01691 (2329)   0.01822 (4521)   0.01858 (8905)

  : Example 2: The $L_2$-norm of the error for increasing values of $s$ and $j$ for the fractional B-spline of degree $\beta=3.5$. The numbers in parenthesis are the degrees-of-freedom. Here, $\gamma=0.5$.

\[tab:conv\_js\_fract\_2\]

                                          3                4                5                6
  -------------------------------- ---------------- ---------------- ---------------- ----------------
   $\sharp V_j^{(\alpha)}([0,1])$         9                17               33               65
                 5                  0.01909 (315)    0.00404 (595)    0.00102 (1155)   0.00063 (2275)
                 6                  0.01810 (603)    0.00546 (1139)   0.00495 (2211)   0.00511 (4355)
                 7                  0.01805 (1179)   0.01671 (2227)   0.01801 (4323)   0.01838 (8515)

  : Example 2: The $L_2$-norm of the error for increasing values of $s$ and $j$ for the cubic B-spline. The numbers in parenthesis are the degrees-of-freedom. Here, $\gamma=0.5$.

\[tab:conv\_js\_cubic\_2\]

Conclusion {#sec:concl}
==========

We proposed a fractional spline collocation-Galerkin method to solve the time-fractional diffusion equation. The novelty of the method is in the use of fractional spline spaces as approximating spaces so that the fractional derivative of the approximating function can be evaluated easily by an explicit differentiation rule that involves the generalized finite difference operator. The numerical tests show that the method has a good accuracy so that it can be effectively used to solve fractional differential problems. The numerical instabilities arising in the fractional basis when $s$ increases can be reduced following the approach in [@GPP04] that allows us to construct stable basis on the interval. Moreover, the ill-conditioning of the linear system (\[colllinearsys\]) can be reduced using iterative methods in Krylov spaces, such as the method proposed in [@CPSV17]. Finally, we notice that following the procedure given in [@GPP04], fractional wavelet bases on finite interval can be constructed so that the proposed method can be generalized to fractional wavelet approximating spaces.

[10]{}

Milton Abramowitz and Irene A. Stegun. , volume 55. Dover Publications, 1965.

Dumitru Baleanu, Kai Diethelm, Enrico Scalas, and Juan J. Trujillo. Fractional calculus. models and numerical methods. , 3:10–16, 2012.

Francesco Calabr[ò]{}, Carla Manni, and Francesca Pitolli. Computation of quadrature rules for integration with respect to refinable functions on assigned nodes. , 90:168–189, 2015.

Daniela Calvetti, Francesca Pitolli, Erkki Somersalo, and Barbara Vantaggi. Bayes meets [K]{}rylov: preconditioning [CGLS]{} for underdetermined systems. , in press.

Wolfgang Dahmen, Siegfried Pr[ö]{}ssdorf, and Reinhold Schneider. Wavelet approximation methods for pseudodifferential equations: I stability and convergence. , 215(1):583–620, 1994.

Neville Ford, Jingyu Xiao, and Yubin Yan. A finite element method for time fractional partial differential equations. , 14(3):454–474, 2011.

Walter Gautschi, Laura Gori, and Francesca Pitolli. Gauss quadrature for refinable weight functions. , 8(3):249–257, 2000.

Laura Gori, Laura Pezza, and Francesca Pitolli. Recent results on wavelet bases on the interval generated by [GP]{} refinable functions. , 51(4):549–563, 2004.

Laura Gori and Francesca Pitolli. Refinable functions and positive operators. , 49(3):381–393, 2004.

Rudolf Hilfer. . World Scientific, 2000.

Francesco Mainardi. . World Scientific, 2010.

Arvet Pedas and Enn Tamme. Numerical solution of nonlinear fractional differential equations by spline collocation methods. , 255:216–230, 2014.

Laura Pezza and Francesca Pitolli. A multiscale collocation method for fractional differential problems. , 147:210–219, 2018.

Igor Podlubny. , volume 198. Academic Press, 1998.

Larry L. Schumaker. . Cambridge University Press, 2007.

Hari Mohan Srivastava and Juan J. Trujillo. . Elsevier, 2006.

Vasily E. Tarasov. . Springer Science & Business Media, 2011.

Vidar Thomée. . Springer-Verlag, 2006.

Michael Unser and Thierry Blu. Fractional splines and wavelets. , 42(1):43–67, 2000.

Mohsen Zayernouri and George Em Karniadakis. Fractional spectral collocation method. , 36(1):A40–A62, 2014.

[^1]: [*Dept. SBAI, University of Roma ”La Sapienza”*]{}, Via A. Scarpa 16, 00161 Roma, Italy. e-mail: [laura.pezza@sbai.uniroma1.it]{}

[^2]: [*Dept. SBAI, University of Roma ”La Sapienza”*]{}, Via A. Scarpa 16, 00161 Roma, Italy. e-mail:
