Hermite Cubic Spline Collocation Method for Nonlinear Fractional Differential Equations with Variable-Order

Abstract

In this paper, we develop a Hermite cubic spline collocation method (HCSCM) for solving variable-order nonlinear fractional differential equations, which apply $C^1$-continuous nodal basis functions to an approximate problem. We also verify that the order of convergence of the HCSCM is about $O\left(h^{min\{4-\alpha ,p\}}\right)$ while the interpolating function belongs to $C^p(p\ge 1)$, where h is the mesh size and $\alpha$ the order of the fractional derivative. Many numerical tests are performed to confirm the effectiveness of the HCSCM for fractional differential equations, which include Helmholtz equations and the fractional Burgers equation of constant-order and variable-order with Riemann-Liouville, Caputo and Patie-Simon sense as well as two-sided cases.

1. Introduction

As a powerful tool for modeling a broad range of non-classical phenomena, fractional calculus has already gained much attention from various science and engineering fields during recent decades. For models of anomalous transport processes and diffusion, there are a lot of fractional partial differential equations proposed in publications [1,2] as well as for the modeling of frequency dependent damping behavior such as in viscoelastic, continuum and statistical mechanics, solid mechanic, economics [3,4], and so on. For modelling the energy supply-demand system, the Caputo-Fabrizio fractional derivative is applied and leads to an interesting fractional energy supply-demand equation [5]. With extensive applications of fractional calculus operators, many fractional differential equations (FDEs) are presented.

Meanwhile, there is increasing demand for a robust method to produce a high accuracy solution to FDEs. Publications on numerical methods for FDEs are largely substantial. A considerable number of them are based on finite difference, see [6,7,8,9,10,11,12,13,14,15,16,17] and the references therein. There are many works based on finite element methods, see [18,19,20,21,22,23]. Methods based on spectral/pseudo-spectral or collocation methods, even the spectral element method, can be seen in [24,25,26,27,28,29,30,31,32,33].

The main challenge of approximating FDEs is the precision deterioration caused by singularity of fractional derivatives [34]. For the spectral-type methods, which are one of the most popular numerical methods due to their high accuracy [35,36,37,38], the singularity of endpoints may damage “the spectral accuracy”. Spectral or collocation methods using fractional polynomials rather than polynomials as basis functions provide a promising way to develop an efficient algorithm for numerically solving FDEs and even fractional operator-related problems. There are theoretical and practical efforts involved in publications such as [29,39,40,41,42,43,44,45,46,47,48].

The demand of flexibility may lead researchers to pursue a multi-domain method or domain decomposition method. A multi-domain spectral collocation method (MDSCM) is suggested to numerically solve FDEs [49]. Authors make use of piecewise continuous ($C^0$) nodal basis functions to approximate the problem. However, a piecewise continuous function may have infinite derivatives of fractional-order. Let us introduce the following result from [49]:

Lemma 1.

Let $\alpha \in (0,2)$ be a constant, and $a,b,c\in \mathbb{R}$ such that $a<c<b$. If $u\in C^2[a,c]\cap C^2[c,b]\cap C[a,b]$ and $u^{\prime }\left(c^{-}\right)$ and $u^{\prime }\left(c^{+}\right)$ exist, then

$${}_{RL}{D}_{a,x}^{\alpha }u\left(x\right)=\sum _{k=0}^1\frac{{(x-a)}^{k-\alpha }}{\Gamma (k+1-\alpha )}{u}^{\left(k\right)}\left(a\right)+\frac{{(x-c)}^{1-\alpha }}{\Gamma (2-\alpha )}[u^{\prime }\left(c^{+}\right)-u^{\prime }\left(c^{-}\right)]+s\left(x\right),$$

for any $x\in (c,b)$, where

$$s\left(x\right)=\left\{\begin{array}{cc}\frac{1}{\Gamma (3-\alpha )}\frac{d}{dx}\left[{\int }_a^x{u}^{″}\left(\tau \right){(x-\tau )}^{2-\alpha }d\tau \right],\hfill & \mathrm{if}\alpha \in (0,1),\hfill \\ \frac{1}{\Gamma (4-\alpha )}\frac{d^2}{d{x}^2}\left[{\int }_a^x{u}^{″}\left(\tau \right){(x-\tau )}^{3-\alpha }d\tau \right],\hfill & \mathrm{if}\alpha \in (1,2).\hfill \end{array}\right.$$

Here, the ${}_{RL}{D}_{a,x}^{\alpha }$ denotes the left Riemann-Liouville fractional derivative of $\alpha$-order; we will give its definition in the following section. The above lemma shows that ${lim}_{x\to c}{}_{RL}{D}_{a,x}^{\alpha }u\left(x\right)=\infty$ if the conditions $u^{\prime }\left(c^{+}\right)\ne u^{\prime }\left(c^{-}\right)$ and $\alpha >1$ are satisfied. To overcome this drawback, $C^1$-continuous nodal basis functions are needed. It is well-known that spline functions are a special class of piecewise polynomials, which provide continuous differentiable solutions over the whole spatial domain with great accuracy. One promising candidate as a $C^1$-continuous nodal basis function is the Hermite cubic spline function.

Spline collocation methods are successfully applied to numerical approximation of differential equations (see [50,51,52] and references therein). However, there are a few publications devoted to the spline collocation method for FDEs. Recently, Liu et al. [53] presented an interesting result of stability and convergence of quadratic spline collocation method for time-dependent fractional diffusion equations. Majeed et al. [54] applied the cubic B-spline collocation method to solve time fractional Burgers’ and Fisher’s equations. Khalid et al. [55] presented a non-polynomial quintic spline collocation method to solve fourth-order fractional boundary value problems involving products terms. Emadifar et al. [56] explored exponential spline interpolation with multiple parameters to find solutions of fractional boundary value problem and conducted the convergence analysis for this technique.

In this paper, our aim is to develop a Hermite cubic spline collocation method (HCSCM) for solving variable-order nonlinear fractional differential equations, which makes use of $C^1$-continuous nodal basis functions to approximate a problem. In particular, the collocation fractional differentiation matrix is derived for fractional derivatives in various senses including Riemann-Liouville, Caputo, Patie-Simon. The main contributions of this work are as follows:

• A set of $C^1$ nodal basis functions are constructed and the corresponding collocation fractional differentiation matrix is derived for the discretization.
• Making use of the Hermite cubic spline collocation method, numerical solution could be found for variable-order nonlinear fractional differential equations. The order of convergence of the HCSCM is also analysed for the left Riemann-Liouville case.
• The effectiveness of the HCSCM is confirmed by solving fractional Helmholtz equations of constant-order and variable-order. With application the HCSCM to the fractional Burgers equation, the numerical fractional diffusion is simulated with different senses.

The paper is organized as follows: in the next Section, some definitions and properties are reviewed for later discussion. The Hermite cubic spline collocation method (HCSCM) is presented in Section 3. The key part is to set up the collocation fractional differentiation matrix. In Section 4, the order of convergence of the HCSCM approximation is analyzed for the left Riemann-Liouville case. Several numerical tests are presented in Section 5. This includes applying HCSCM to fractional Helmholtz equations and fractional Burgers equations. Finally, we conclude in Section 6.

2. Preliminaries

In this Section, some definitions of fractional calculus are reviewed for subsequent discussions. The most common-used definitions of fractional derivatives are possibly the Riemann-Liouville’s and the Caputo’s, found in various publications, such as ([57,58]). The following definitions are variable-order versions, which provide constant-order definitions when $\alpha \left(x\right)\equiv \alpha$ is a constant in the formulas.

Definition 1.

For a function $f\left(x\right),x\in [x_L,x_R]$, the left Riemann-Liouville fractional integral of order $\alpha \left(x\right)>0$ is defined as

$${}_{x_L}{I}_x^{\alpha \left(x\right)}f\left(x\right):=\frac{1}{\Gamma \left(\alpha \right(x\left)\right)}{\int }_{x_L}^x{(x-s)}^{\alpha \left(x\right)-1}f\left(s\right)ds,$$

(1)

and the right Riemann-Liouville fractional integral of order $\alpha \left(x\right)>0$ is defined as

$${}_x{I}_{x_R}^{\alpha \left(x\right)}f\left(x\right):=\frac{1}{\Gamma \left(\alpha \right(x\left)\right)}{\int }_x^{x_R}{(s-x)}^{\alpha \left(x\right)-1}f\left(s\right)ds,$$

(2)

where $\Gamma (·)$ is the Euler’s gamma function.

Definition 2.

For a function $f\left(x\right),x\in [x_L,x_R]$, the left Riemann-Liouville fractional derivative of order $\alpha \left(x\right)>0$ is defined as

$${}_{RL}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right):=\frac{1}{\Gamma (n-\alpha (x\left)\right)}{\left[\frac{d^n}{d{\xi }^n}{\int }_{x_L}^{\xi }{(\xi -s)}^{n-\alpha \left(x\right)-1}f\left(s\right)ds\right]}_{\xi =x},$$

(3)

and the right Riemann-Liouville fractional derivative of order $\alpha \left(x\right)>0$ is defined as

$${}_{RL}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right):=\frac{{(-1)}^n}{\Gamma (n-\alpha (x\left)\right)}{\left[\frac{d^n}{d{\xi }^n}{\int }_{\xi }^{x_R}{(s-\xi )}^{n-\alpha \left(x\right)-1}f\left(s\right)ds\right]}_{\xi =x},$$

(4)

where n is the positive integer such that $n-1<\alpha \left(x\right)<n$.

Definition 3.

For a function $f\left(x\right),x\in [x_L,x_R]$, the left Caputo fractional derivative of order $\alpha \left(x\right)>0$ is defined as

$${}_C{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right):=\frac{1}{\Gamma (n-\alpha (x\left)\right)}{\int }_{x_L}^x{(x-s)}^{n-\alpha \left(x\right)-1}{f}^{\left(n\right)}\left(s\right)ds,$$

(5)

and the right Caputo fractional derivative of order $\alpha \left(x\right)>0$ is defined as

$${}_C{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right):=\frac{{(-1)}^n}{\Gamma (n-\alpha (x\left)\right)}{\int }_x^{x_R}{(s-x)}^{n-\alpha \left(x\right)-1}{f}^{\left(n\right)}\left(s\right)ds,$$

(6)

where n is the positive integer such that $n-1<\alpha \left(x\right)<n$.

The well-known relationship between Riemann-Liouville and the Caputo derivative is as follows:

Lemma 2.

If ${}_{RL}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right),{}_C{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right),{}_{RL}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)$ and ${}_C{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)$ exist, then

$${}_{RL}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)={}_C{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)+\sum _{k=0}^{n-1}\frac{f^{\left(k\right)}\left(x_L\right)}{\Gamma (k+1-\alpha (x\left)\right)}{(x-x_L)}^{k-\alpha \left(x\right)},$$

(7)

and

$${}_{RL}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)={}_C{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)+\sum _{k=0}^{n-1}\frac{{(-1)}^{n-j}{f}^{\left(k\right)}\left(x_R\right)}{\Gamma (k+1-\alpha (x\left)\right)}{(x_R-x)}^{k-\alpha \left(x\right)}.$$

(8)

Besides the common-used definitions above, the fractional diffusion operators which limit the order $1<\alpha \left(x\right)\le 2$ are also considered. A definition was proposed by Patie and Simon in [59] as follows.

Definition 4.

For a function $f\left(x\right),x\in [x_L,x_R]$, the left Patie-Simon (or mixed Caputo) fractional derivative of order $1<\alpha \left(x\right)<2$ is defined as

$${}_{PS}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right):=\frac{1}{\Gamma (2-\alpha (x\left)\right)}{\left[\frac{d}{d\xi }{\int }_{x_L}^{\xi }{(\xi -s)}^{1-\alpha \left(x\right)}{f}^{{}^{\prime }}\left(s\right)ds\right]}_{\xi =x},$$

(9)

and the right Patie-Simon (or mixed Caputo) fractional derivative of order $1<\alpha \left(x\right)<2$ is defined as

$${}_{PS}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right):=\frac{1}{\Gamma (2-\alpha (x\left)\right)}{\left[\frac{d}{d\xi }{\int }_{\xi }^{x_R}{(s-\xi )}^{1-\alpha \left(x\right)}{f}^{{}^{\prime }}\left(s\right)ds\right]}_{\xi =x}.$$

(10)

From the above definitions and Lemma 2, hold the following relationships:

Lemma 3.

If $1<\alpha \left(x\right)<2$ and ${}_{RL}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)$, ${}_C{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)$, ${}_{PS}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)$, ${}_{RL}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)$, ${}_C{D}_{x,x_R}^{\alpha \left(x\right)}$ and ${}_{PS}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)$ exist, then

$${}_{RL}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)={}_{PS}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)+\frac{f\left(x_L\right)}{\Gamma (1-\alpha (x\left)\right)}{(x-x_L)}^{-\alpha \left(x\right)},$$

(11)

and

$${}_{RL}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)={}_{PS}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)+\frac{f\left(x_R\right)}{\Gamma (1-\alpha (x\left)\right)}{(x_R-x)}^{-\alpha \left(x\right)},$$

(12)

and

$${}_{PS}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)={}_C{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)+\frac{f^{\prime }\left(x_L\right)}{\Gamma (2-\alpha (x\left)\right)}{(x-x_L)}^{1-\alpha \left(x\right)},$$

(13)

and

$${}_{PS}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)={}_C{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)+\frac{f^{\prime }\left(x_R\right)}{\Gamma (2-\alpha (x\left)\right)}{(x_R-x)}^{1-\alpha \left(x\right)}.$$

(14)

Proof. 

Since

$${\int }_{x_L}^{\xi }{(\xi -s)}^{1-\alpha \left(x\right)}f\left(s\right)ds=\frac{{(\xi -x_L)}^{2-\alpha \left(x\right)}}{2-\alpha \left(x\right)}f\left(x_L\right)+\frac{1}{2-\alpha \left(x\right)}{\int }_{x_L}^{\xi }{(\xi -s)}^{2-\alpha \left(x\right)}{f}^{\prime }\left(s\right)ds.$$

Then note that $1<\alpha \left(x\right)<2$,

$$\begin{array}{ccc}& & \frac{d^2}{d{\xi }^2}\left[{\int }_{x_L}^{\xi }{(\xi -s)}^{1-\alpha \left(x\right)}f\left(s\right)ds\right]\hfill \\ & =& (1-\alpha \left(x\right)){(\xi -x_L)}^{-\alpha \left(x\right)}f\left(x_L\right)+\frac{d}{d\xi }{\int }_{x_L}^{\xi }{(\xi -s)}^{1-\alpha \left(x\right)}{f}^{\prime }\left(s\right)ds].\hfill \end{array}$$

The equality (11) is obtained by dividing factor $\Gamma (2-\alpha (x\left)\right)$. Other results can be derived by a similar argument. □

There exist the following well-known properties:

Lemma 4.

Let m be an integer number, the following properties hold for $x\in [x_L,x_R]$ and Riemann-Liouville fractional calculus

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_{x_L}{I}_x^{\alpha \left(x\right)}{(x-x_L)}^m& =\frac{m!}{\Gamma \left(m+\alpha \left(x\right)+1\right)}{(x-x_L)}^{m+\alpha \left(x\right)},\hfill \\ \hfill {}_x{I}_{x_R}^{\alpha \left(x\right)}{(x_R-x)}^m& =\frac{m!}{\Gamma \left(m+\alpha \left(x\right)+1\right)}{(x_R-x)}^{m+\alpha \left(x\right)},\hfill \\ \hfill {}_{RL}{D}_{x_L,x}^{\alpha \left(x\right)}{(x-x_L)}^m& =\frac{m!}{\Gamma \left(m-\alpha \left(x\right)+1\right)}{(x-x_L)}^{m-\alpha \left(x\right)},\hfill \\ \hfill {}_{RL}{D}_{x,x_R}^{\alpha \left(x\right)}{(x_R-x)}^m& =\frac{m!}{\Gamma \left(m-\alpha \left(x\right)+1\right)}{(x_R-x)}^{m-\alpha \left(x\right)},\hfill \end{array}\end{array}$$

(15)

and for the Caputo fractional derivative,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_C{D}_{x_L,x}^{\alpha \left(x\right)}{(x-x_L)}^m& =\left\{\begin{array}{cc}\frac{m!}{\Gamma \left(m-\alpha \left(x\right)+1\right)}{(x-x_L)}^{m-\alpha \left(x\right)},\hfill & ifm>\alpha \left(x\right),\hfill \\ 0,\hfill & ifm<\alpha \left(x\right),\hfill \end{array}\right.\hfill \\ \hfill {}_C{D}_{x,x_R}^{\alpha \left(x\right)}{(x_R-x)}^m& =\left\{\begin{array}{cc}\frac{m!}{\Gamma \left(m-\alpha \left(x\right)+1\right)}{(x_R-x)}^{m-\alpha \left(x\right)},\hfill & ifm>\alpha \left(x\right),\hfill \\ 0,\hfill & ifm<\alpha \left(x\right),\hfill \end{array}\right.\hfill \end{array}\end{array}$$

(16)

and for the Patie-Simon fractional derivative of $1<\alpha \left(x\right)<2$,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_{PS}{D}_{x_L,x}^{\alpha \left(x\right)}{(x-x_L)}^m& =\left\{\begin{array}{cc}\frac{m!}{\Gamma \left(m-\alpha \left(x\right)+1\right)}{(x-x_L)}^{m-\alpha \left(x\right)},\hfill & ifm>0,\hfill \\ 0,\hfill & ifm=0,\hfill \end{array}\right.\hfill \\ \hfill {}_{PS}{D}_{x,x_R}^{\alpha \left(x\right)}{(x_R-x)}^m& =\left\{\begin{array}{cc}\frac{m!}{\Gamma \left(m-\alpha \left(x\right)+1\right)}{(x_R-x)}^{m-\alpha \left(x\right)},\hfill & ifm>0,\hfill \\ 0,\hfill & ifm=0.\hfill \end{array}\right.\hfill \end{array}\end{array}$$

(17)

The following operators with top-tilde are useful in HCSCM for $x>x_R$,

$${}_{x_L\ast }{\tilde{I}}_{x_R}^{\alpha \left(x\right)}f\left(x\right):=\frac{1}{\Gamma \left(\alpha \right(x\left)\right)}{\int }_{x_L}^{x_R}{(x-s)}^{\alpha \left(x\right)-1}f\left(s\right)ds,$$

(18)

and for $x<x_L$,

$${}_{x_L}{\tilde{I}}_{x_R\ast }^{\alpha \left(x\right)}f\left(x\right):=\frac{1}{\Gamma \left(\alpha \right(x\left)\right)}{\int }_{x_L}^{x_R}{(s-x)}^{\alpha \left(x\right)-1}f\left(s\right)ds,$$

(19)

and for $x>x_R$,

$${}_{RL}{\tilde{D}}_{x_L\ast ,x_R}^{\alpha \left(x\right)}f\left(x\right):=\frac{1}{\Gamma (n-\alpha (x\left)\right)}{\left[\frac{d^n}{d{\xi }^n}{\int }_{x_L}^{x_R}{(\xi -s)}^{n-\alpha \left(x\right)-1}f\left(s\right)ds\right]}_{\xi =x},$$

(20)

and for $x<x_L$,

$${}_{RL}{\tilde{D}}_{x_L,x_R\ast }^{\alpha \left(x\right)}f\left(x\right):=\frac{{(-1)}^n}{\Gamma (n-\alpha (x\left)\right)}{\left[\frac{d^n}{d{\xi }^n}{\int }_{x_L}^{x_R}{(s-\xi )}^{n-\alpha \left(x\right)-1}f\left(s\right)ds\right]}_{\xi =x}.$$

(21)

Operators ${}_C{\tilde{D}}_{x_L\ast ,x_R}^{\alpha \left(x\right)}$, ${}_C{\tilde{D}}_{x_L,x_R\ast }^{\alpha \left(x\right)}$, ${}_{PS}{\tilde{D}}_{x_L\ast ,x_R}^{\alpha \left(x\right)}$, ${}_{PS}{\tilde{D}}_{x_L,x_R\ast }^{\alpha \left(x\right)}$ are defined similarly.

Lemma 5.

Let $x_L<x_c<x_R$ and $x\in (x_c,x_R]$, then

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {{}_{x_L}I}_x^{\alpha \left(x\right)}f\left(x\right)& ={{}_{x_L\ast }\tilde{I}}_{x_c}^{\alpha \left(x\right)}f\left(x\right)+{{}_{x_c}I}_x^{\alpha \left(x\right)}f\left(x\right),\hfill \\ \hfill {}_{RL}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)& ={{}_{RL}\tilde{D}}_{x_L\ast ,x_c}^{\alpha \left(x\right)}f\left(x\right)+{{}_{RL}D}_{x_c,x}^{\alpha \left(x\right)}f\left(x\right),\hfill \\ \hfill {}_C{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)& ={{}_C\tilde{D}}_{x_L\ast ,x_c}^{\alpha \left(x\right)}f\left(x\right)+{{}_{C}D}_{x_c,x}^{\alpha \left(x\right)}f\left(x\right),\hfill \\ \hfill {}_{PS}{D}_{x_L,x}^{\alpha \left(x\right)}f\left(x\right)& ={{}_{PS}\tilde{D}}_{x_L\ast ,x_c}^{\alpha \left(x\right)}f\left(x\right)+{{}_{PS}D}_{x_c,x}^{\alpha \left(x\right)}f\left(x\right),\hfill \end{array}\end{array}$$

(22)

and when $x\in [x_L,x_c)$, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {{}_{x}I}_{x_R}^{\alpha \left(x\right)}f\left(x\right)& ={{}_{x}I}_{x_c}^{\alpha \left(x\right)}f\left(x\right)+{{}_{x_c}\tilde{I}}_{x_R\ast }^{\alpha \left(x\right)}f\left(x\right),\hfill \\ \hfill {}_{RL}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)& ={{}_{RL}D}_{x,x_c}^{\alpha \left(x\right)}f\left(x\right)+{{}_{RL}\tilde{D}}_{x_c,x_R\ast }^{\alpha \left(x\right)}f\left(x\right),\hfill \\ \hfill {}_C{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)& ={{}_{C}D}_{x,x_c}^{\alpha \left(x\right)}f\left(x\right)+{{}_C\tilde{D}}_{x_c,x_R\ast }^{\alpha \left(x\right)}f\left(x\right),\hfill \\ \hfill {}_{PS}{D}_{x,x_R}^{\alpha \left(x\right)}f\left(x\right)& ={{}_{PS}D}_{x,x_c}^{\alpha \left(x\right)}f\left(x\right)+{{}_{PS}\tilde{D}}_{x_c,x_R\ast }^{\alpha \left(x\right)}f\left(x\right).\hfill \end{array}\end{array}$$

(23)

Proof. 

Since

$${\int }_{x_L}^x{(x-s)}^{\alpha \left(x\right)-1}f\left(s\right)ds={\int }_{x_L}^{x_c}{(x-s)}^{\alpha \left(x\right)-1}f\left(s\right)ds+{\int }_{x_c}^x{(x-s)}^{\alpha \left(x\right)-1}ds,$$

Then the first equality in (22) is obtained by dividing factor $\Gamma \left(\alpha \right(x\left)\right)$ and the definitions (1) and (18). Other results can be derived by a similar argument. □

If ${\mathcal{D}}^{\ast \alpha \left(x\right)}$ and ${\mathcal{D}}^{\alpha \left(x\right)\ast }$ represent all left-sided and right-sided definitions of the above-mentioned, respectively, then the two-sided fractional derivative can be written as

$${\mathbb{D}}_r^{\alpha \left(x\right)}:=r{\mathcal{D}}^{\ast \alpha \left(x\right)}+(1-r){\mathcal{D}}^{\alpha \left(x\right)\ast },0\le r\le 1.$$

(24)

3. Hermite Cubic Spline Collocation Method (HCSCM)

In the Section, the HCSCM is presented. The key role of HCSCM is the collocation fractional differentiation matrix.

3.1. Fractional Differentiation Matrix (FDM) for HCSCM

Let $\Lambda :=(x_L,x_R)$, the first step is to divide the interval $\Lambda$ into N elements, that is,

$$x_L=x_0<x_1<\dots <x_N=x_R.$$

Denote $I_i=[x_{i-1},x_i],i=1,2,\dots ,N$ the $i-$th element and $h_i=x_i-x_{i-1}$ the length of $I_i$. Let ${\mathbb{P}}_N^I$ be the collection of all algebraic polynomials defined on interval I with degree at most N. The piecewise Hermite cubic polynomial space is

$${\mathbb{V}}_N=\{v\in C^1(\Lambda ):v{|}_{I_i}\in {\mathbb{P}}_3^{I_i},i=1,2,\dots ,N\}.$$

which is defined by the following set of nodal basis functions. It contains $2N+2$ functions as follows. The first two functions as

$${\phi }_0\left(x\right)=\left\{\begin{array}{cc}\left(1+2\frac{x-x_0}{h_1}\right){\left(1-\frac{x-x_0}{h_1}\right)}^2,\hfill & ifx\in I_1,\hfill \\ 0,\hfill & otherwise,\hfill \end{array}\right.$$

and

$${\varphi }_0\left(x\right)=\left\{\begin{array}{cc}\left(\frac{x-x_0}{h_1}\right){\left(1-\frac{x-x_0}{h_1}\right)}^2{h}_1,\hfill & ifx\in I_1,\hfill \\ 0,\hfill & otherwise.\hfill \end{array}\right.$$

For $i=1,2,\dots ,N-1,$

$${\phi }_i\left(x\right)=\left\{\begin{array}{cc}\left(3-2\frac{x-x_{i-1}}{h_i}\right){\left(\frac{x-x_{i-1}}{h_i}\right)}^2,\hfill & ifx\in I_i,\hfill \\ \left(1+2\frac{x-x_i}{h_{i+1}}\right){\left(1-\frac{x-x_i}{h_{i+1}}\right)}^2,\hfill & ifx\in I_{i+1},\hfill \\ 0,\hfill & otherwise,\hfill \end{array}\right.$$

and

$${\varphi }_i\left(x\right)=\left\{\begin{array}{cc}-\left(1-\frac{x-x_{i-1}}{h_i}\right){\left(\frac{x-x_{i-1}}{h_i}\right)}^2{h}_i,\hfill & ifx\in I_i,\hfill \\ \left(\frac{x-x_i}{h_{i+1}}\right){\left(1-\frac{x-x_i}{h_{i+1}}\right)}^2{h}_{i+1},\hfill & ifx\in I_{i+1},\hfill \\ 0,\hfill & otherwise.\hfill \end{array}\right.$$

and the last two functions as

$${\phi }_N\left(x\right)=\left\{\begin{array}{cc}\left(3-2\frac{x-x_{N-1}}{h_N}\right){\left(\frac{x-x_{N-1}}{h_N}\right)}^2,\hfill & ifx\in I_N,\hfill \\ 0,\hfill & otherwise,\hfill \end{array}\right.$$

and

$${\varphi }_N\left(x\right)=\left\{\begin{array}{cc}-\left(1-\frac{x-x_{N-1}}{h_N}\right){\left(\frac{x-x_{N-1}}{h_N}\right)}^2{h}_N,\hfill & ifx\in I_N,\hfill \\ 0,\hfill & otherwise.\hfill \end{array}\right.$$

Therefore,

$${\mathbb{V}}_N=\mathrm{span}\{{\phi }_i,{\varphi }_i,i=0,1,\dots ,N\}.$$

If $u_N\in {\mathbb{V}}_N$, then can be expanded as

$$u_N\left(x\right)=\sum _{i=0}^N\left(u_N\left(x_i\right){\phi }_i\left(x\right)+u_N^{\prime }\left(x_i\right){\varphi }_i\left(x\right)\right).$$

In each element $I_i$, $x_{i,1}^c,x_{i,2}^c\in I_i$ are the collocation points, where

$$x_{i,1}^c=x_{i-1}+{\sigma }_{i,1}{h}_i,x_{i,2}^c=x_{i-1}+{\sigma }_{i,2}{h}_i,i=1,2,\dots ,N,$$

and $0\le {\sigma }_{i,1}<{\sigma }_{i,2}\le 1$. In fact, a choice of this points is the Gauss-type quadrature nodes, ${\sigma }_{i,1}=(1-{\sigma }_{i,2})=\frac{\sqrt{3}}{3}$, which is named by orthogonal spline collocation. However, the stable collocation points may not be symmetric in the view of [60,61].

As a collocation approximation to the $\alpha \left(x\right)$th-order differential operators defined in Section 2, we denote by ${\mathbf{D}}^{\alpha }$ the collocation fractional differentiation matrix, which satisfies

$${\left({\mathbf{D}}^{\alpha }{u}_N\right)}_l={\mathcal{D}}^{\alpha \left(x_{ij}^c\right)}{u}_N(x_{ij}^c),j=1,2;i=1,2,\dots ,N.$$

(25)

The structure of the collocation fractional differentiation matrix (FDM) may differ with the ordering of the collocation points and the unknowns. In natural ordering $l=2(i-1)+j$, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{u}& ={[u_0^{\prime },u_1,u_1^{\prime },\dots ,u_{N-1},u_{N-1}^{\prime },u_N^{\prime }]}^T,\hfill \\ \hfill {\mathbf{x}}^c& ={[x_{11}^c,x_{12}^c,x_{21}^c,x_{22}^c,\dots ,x_{N1}^c,x_{N2}^c]}^T.\hfill \end{array}\end{array}$$

(26)

and ${\mathbf{D}}^{\alpha }$ with Dirichlet boundary conditions is

$${\mathbf{D}}^{\alpha }=\left[\begin{array}{ccccc}\mathcal{D}{\varphi }_0\left(x_{11}^c\right)& \mathcal{D}{\phi }_1\left(x_{11}^c\right)& \mathcal{D}{\varphi }_1\left(x_{11}^c\right)& \dots & \mathcal{D}{\varphi }_N\left(x_{11}^c\right)\\ \mathcal{D}{\varphi }_0\left(x_{12}^c\right)& \mathcal{D}{\phi }_1\left(x_{12}^c\right)& \mathcal{D}{\varphi }_1\left(x_{12}^c\right)& \dots & \mathcal{D}{\varphi }_N\left(x_{12}^c\right)\\ ⋮& ⋮& ⋮& & ⋮\\ \mathcal{D}{\varphi }_0\left(x_{N1}^c\right)& \mathcal{D}{\phi }_1\left(x_{N1}^c\right)& \mathcal{D}{\varphi }_1\left(x_{N1}^c\right)& \dots & \mathcal{D}{\varphi }_N\left(x_{N1}^c\right)\\ \mathcal{D}{\varphi }_0\left(x_{N2}^c\right)& \mathcal{D}{\phi }_1\left(x_{N2}^c\right)& \mathcal{D}{\varphi }_1\left(x_{N2}^c\right)& \dots & \mathcal{D}{\varphi }_N\left(x_{N2}^c\right)\end{array}\right],$$

(27)

here $\mathcal{D}={\mathcal{D}}^{\alpha (x_{ij}^c)}$ is one of the fractional differential operators defined in Section 2. Typically, the matrix $\mathcal{D}$ is block-triangular for left and right fractional operators.

Remark 1.

According to the Lemma 2, Lemma 3 and the special nodal basis functions, the collocation FDM ${\mathbf{D}}^{\alpha }$ of the Riemann-Liouville operators is equal to the corresponding FDM of the Patie-Simon operators. For the Caputo operators, the corresponding FDM is different only from the first or last column.

3.2. Computing the Entries of FDM

For ease of computing, operations are shifted from an arbitrary interval $[a,b]$ to the reference interval $[-1,1]$. Let the linear transformation

$$x=\frac{h}{2}\left(y+1\right)+a,\mathrm{or}y=\frac{2}{h}\left(x-a\right)-1,h:=b-a$$

(28)

shift functions $f\left(x\right),\alpha \left(x\right)$ defined on the interval $[a,b]$ to $\widehat{f}\left(y\right),\widehat{\alpha }\left(y\right)$ on the reference interval $[-1,1]$. Then we have the following relations:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_a{I}_x^{\alpha \left(x\right)}f\left(x\right)& ={\left(\frac{h}{2}\right)}^{\widehat{\alpha }\left(y\right)}{{}_{-1}{I}_y}^{\widehat{\alpha }\left(y\right)}\widehat{f}\left(y\right),\hfill \\ \hfill {}_x{I}_b^{\alpha \left(x\right)}f\left(x\right)& ={\left(\frac{h}{2}\right)}^{\widehat{\alpha }\left(y\right)}{{}_y{I}_1}^{\widehat{\alpha }\left(y\right)}\widehat{f}\left(y\right),\hfill \\ \hfill {\mathcal{D}}^{\alpha \left(x\right)}f\left(x\right)& ={\left(\frac{2}{h}\right)}^{\widehat{\alpha }\left(y\right)}{\mathcal{D}}^{\widehat{\alpha }\left(y\right)}\widehat{f}\left(y\right),\hfill \end{array}\end{array}$$

(29)

here ${\mathcal{D}}^{\alpha \left(x\right)}$ be one of the seven: ${{}_{RL}D}_{a,x}^{\alpha \left(x\right)}, {{}_{RL}D}_{x,b}^{\alpha \left(x\right)}$, ${{}_{C}D}_{a,x}^{\alpha \left(x\right)}, {{}_{C}D}_{x,b}^{\alpha \left(x\right)}$, ${{}_{PS}D}_{a,x}^{\alpha \left(x\right)}$, ${{}_{PS}D}_{x,b}^{\alpha \left(x\right)}$, ${\mathbb{D}}_r^{\alpha \left(x\right)}.$

For the tilde operators, the following relations also hold:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_{a\ast }{\tilde{I}}_b^{\alpha \left(x\right)}f\left(x\right)& ={\left(\frac{h}{2}\right)}^{\widehat{\alpha }\left(y\right)}{{}_{-1\ast }\tilde{I}}_1^{\widehat{\alpha }\left(y\right)}\widehat{f}\left(y\right),\hfill \\ \hfill {}_a{\tilde{I}}_{b\ast }^{\alpha \left(x\right)}f\left(x\right)& ={\left(\frac{h}{2}\right)}^{\widehat{\alpha }\left(y\right)}{{}_{-1}{\tilde{I}}_{1\ast }}^{\widehat{\alpha }\left(y\right)}\widehat{f}\left(y\right),\hfill \\ \hfill {\tilde{\mathcal{D}}}^{\alpha \left(x\right)}f\left(x\right)& ={\left(\frac{2}{h}\right)}^{\widehat{\alpha }\left(y\right)}{\tilde{\mathcal{D}}}^{\widehat{\alpha }\left(y\right)}\widehat{f}\left(y\right),\hfill \end{array}\end{array}$$

(30)

here ${\tilde{\mathcal{D}}}^{\alpha \left(x\right)}$ can be one of the six: ${}_{RL}{\tilde{D}}_{a\ast ,b}^{\alpha \left(x\right)}$, ${}_{RL}{\tilde{D}}_{a,b\ast }^{\alpha \left(x\right)}$, ${}_C{\tilde{D}}_{a\ast ,b}^{\alpha \left(x\right)}$, ${}_C{\tilde{D}}_{a,b\ast }^{\alpha \left(x\right)}$, ${}_{PS}{\tilde{D}}_{a\ast ,b}^{\alpha \left(x\right)}$, ${}_{PS}{\tilde{D}}_{a,b\ast }^{\alpha \left(x\right)}$.

The nodal basis functions presented in Section 3, are the so-called shape functions after being transferred by (28), that is,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\xi }_1\left(y\right):& =(2+y){\left(1-\frac{y+1}{2}\right)}^2,\hfill \\ \hfill {\xi }_2\left(y\right):& =\left(\frac{y+1}{2}\right){\left(1-\frac{y+1}{2}\right)}^2,\left(except factor{h}_i\right),\hfill \\ \hfill {\xi }_3\left(y\right):& =\left(2-y\right){\left(\frac{y+1}{2}\right)}^2,\hfill \\ \hfill {\xi }_4\left(y\right):& =-\left(1-\frac{y+1}{2}\right){\left(\frac{y+1}{2}\right)}^2,\left(except factor{h}_i\right),\hfill \end{array}\end{array}$$

(31)

and $y\in [-1,1]$. The Hermite Spline collocation method will perform all the operators mentioned above on the shape functions (31).

4. Order of Convergence of the Approximation with HCSCM

In this Section, the order of convergence of the approximation with HCSCM is analysed. Typically, the left Riemann-Liouville fractional derivative is considered. For convenience of analysis, denote $D^{\alpha }={}_{RL}{D}_{x_L,x}^{\alpha \left(x\right)}$ and let $h_i=x_i-x_{i-1}=h,{\sigma }_{i,1}={\sigma }_1,{\sigma }_{i,2}={\sigma }_2,i=1,2,\dots ,N$. Then $x_i=x_0+ih,i=0,1,\dots ,N$ and the collocation points

$$x_{i,1}^c=x_0+(i-1+{\sigma }_1)h,x_{i,2}^c=x_0+(i-1+{\sigma }_2)h,i=1,2,\dots ,N.$$

Let ${\Pi }_N:C^1(\Lambda )\to {\mathbb{V}}_N$ be the piecewise Hermite cubic interpolation operator, determined uniquely by

$${\Pi }_{N}f\left(x_i\right)=f\left(x_i\right),\frac{d}{dx}\left({\Pi }_{N}f\right)\left(x_i\right)=f^{\prime }\left(x_i\right),i=0,1,\dots ,N,$$

for every $f\in C^1(\Lambda )$.

For a function $u\left(x\right)\in C(\Lambda )$, the maximum norm is defined by

$${\parallel u\parallel }_{\infty }=\underset{x\in \Lambda }{max}\left|u\left(x\right)\right|.$$

The following results are related to the interpolation errors [62].

Lemma 6.

Let $u\left(x\right)\in C^4(\Lambda )$. Then

$${∥\frac{d^j}{d{x}^j}\left(u-{\Pi }_{N}u\right)∥}_{\infty }\le C{h}^{4-j}{\parallel u^{\left(4\right)}\parallel }_{\infty },0\le j\le 3,$$

(32)

where C is a constant number which do not dependent on N.

If $u\left(x\right)\in C^p(\Lambda )(p\ge 1)$, the interpolation error holds (see [63]):

$${∥\frac{d^j}{d{x}^j}\left(u-{\Pi }_{N}u\right)∥}_{\infty }\le C{h}^{s-j}{\parallel u^{\left(s\right)}\parallel }_{\infty },0\le j\le p.$$

(33)

where $s=min\{p,4\}$.

In the following, the error bound is presented for $\parallel D^{\alpha }({\Pi }_{N}u-u){\parallel }_{\infty }$ with constant-order $\alpha \in (1,2)$. Let $\tau \left(x\right)={\Pi }_{N}u-u$, we have

$$D^{\alpha }\tau \left(x\right)=\frac{1}{\Gamma (2-\alpha )}\frac{d^2}{d{x}^2}{\int }_{x_L}^x\frac{\tau \left(s\right)}{{(x-s)}^{\alpha -1}}ds.$$

Assume that $x\in [x_{j-1},x_j]$ for some j, then the above integration can be split as

$${\int }_{x_L}^x\frac{\tau \left(s\right)}{{(x-s)}^{\alpha -1}}ds=\sum _{i=1}^{j-1}{\int }_{x_{i-1}}^{x_i}\frac{\tau \left(s\right)}{{(x-s)}^{\alpha -1}}ds+{\int }_{x_{j-1}}^x\frac{\tau \left(s\right)}{{(x-s)}^{\alpha -1}}ds.$$

(34)

Let ${\overline{x}}_i=(x_{i-1}+x_i)/2,i=1,2,\dots ,N.$ Under the assumption of $u\left(x\right)\in C^4(\Lambda )$, from Taylor’s theorem we have

$$\begin{array}{cc}\hfill & \frac{d^2}{d{x}^2}{\int }_{x_{i-1}}^{x_i}\frac{\tau \left(s\right)}{{(x-s)}^{\alpha -1}}ds=\sum _{k=0}^3\frac{{\tau }^{\left(k\right)}\left({\overline{x}}_i\right)}{k!}\frac{d^2}{d{x}^2}{\int }_{x_{i-1}}^{x_i}\frac{{(s-{\overline{x}}_i)}^k}{{(x-s)}^{\alpha -1}}ds\hfill \\ \hfill & +\frac{1}{24}\frac{d^2}{d{x}^2}{\int }_{x_{i-1}}^{x_i}\frac{u^{\left(4\right)}\left({\zeta }_i\right){(s-{\overline{x}}_i)}^4}{{(x-s)}^{\alpha -1}}ds=:\sum _{k=0}^3{J}_k\left(x\right)+J_4\left(x\right),\hfill \end{array}$$

(35)

where ${\zeta }_i\in (x_{i-1},x_i)$. Now from the Mean Value Theorem for integrals for $k=0,1,2,3$

$$J_k\left(x\right)=\frac{{\tau }^{\left(k\right)}\left({\overline{x}}_i\right)}{k!}{({\widehat{\zeta }}_{i,k}-{\overline{x}}_i)}^k(\alpha -1)[{(x-x_i)}^{-\alpha }-{(x-x_{i-1})}^{-\alpha }],$$

(36)

and

$$J_4\left(x\right)=\frac{u^{\left(4\right)}({\tilde{\zeta }}_i){({\widehat{\zeta }}_{i,4}-{\overline{x}}_i)}^4}{24}(\alpha -1)\left[{(x-x_i)}^{-\alpha }-{(x-x_{i-1})}^{-\alpha }\right],$$

(37)

where ${\tilde{\zeta }}_i,{\widehat{\zeta }}_{i,k}\in (x_{i-1},x_i).$ For the second integral in (34), we have

$$\begin{array}{cc}\hfill & \frac{d^2}{d{x}^2}{\int }_{x_{j-1}}^x\frac{\tau \left(s\right)}{{(x-s)}^{\alpha -1}}ds=\tau \left({\overline{x}}_j\right)(1-\alpha ){(x-x_{j-1})}^{-\alpha }\hfill \\ \hfill & +{\tau }^{\prime }\left({\overline{x}}_j\right)\left[-\frac{h}{2}(1-\alpha ){(x-x_{j-1})}^{-\alpha }+{(x-x_{j-1})}^{1-\alpha }\right]\hfill \\ \hfill & +\frac{{\tau }^{″}\left({\overline{x}}_j\right)}{2}[\frac{h^2}{4}(1-\alpha ){(x-x_{j-1})}^{-\alpha }-h{(x-x_{j-1})}^{1-\alpha }\hfill \\ \hfill & +\frac{2}{2-\alpha }{(x-x_{j-1})}^{2-\alpha }]\hfill \\ \hfill & +\frac{{\tau }^{‴}\left({\overline{x}}_j\right)}{6}[-\frac{h^3}{8}(1-\alpha ){(x-x_{j-1})}^{-\alpha }+\frac{3h^2}{4}{(x-x_{j-1})}^{1-\alpha }\hfill \\ \hfill & +\frac{3h}{2-\alpha }{(x-x_{j-1})}^{2-\alpha }+\frac{6}{(2-\alpha )(3-\alpha )}{(x-x_{j-1})}^{3-\alpha }]\hfill \\ \hfill & +\frac{u^{\left(4\right)}\left({\tilde{\zeta }}_j\right){({\widehat{\zeta }}_j-{\overline{x}}_j)}^4}{24}(1-\alpha ){(x-x_{j-1})}^{-\alpha }=:J\left(x\right).\hfill \end{array}$$

(38)

Let $\sigma h=x-x_{j-1}$. Note that $\sigma \in (0,1)$, it is easy to know that

$${\sigma }^{-\alpha }>{\sigma }^{-\alpha }-{(\sigma +1)}^{-\alpha }>{(\sigma +1)}^{-\alpha }-{(\sigma +2)}^{-\alpha }>\dots$$

(39)

Hence, by Lemma 6, for $k=0,1,2,3$ and every $i<j$ we have

$$|J_k\left(x\right)|\le \frac{|{\tau }^{\left(k\right)}\left({\overline{x}}_i\right)|}{k!}{h}^{k-\alpha }(\alpha -1)\left[{\sigma }^{-\alpha }-{(\sigma +1)}^{-\alpha }\right]\le C{h}^{4-\alpha }{\parallel u^{\left(4\right)}\parallel }_{\infty },$$

(40)

and

$$|J_4\left(x\right)|\le \frac{|u^{(k)}({\tilde{\zeta }}_i)|}{4!}{h}^{4-\alpha }(\alpha -1)\left[{\sigma }^{-\alpha }-{(\sigma +1)}^{-\alpha }\right]\le C{h}^{4-\alpha }{\parallel u^{(4)}\parallel }_{\infty },$$

(41)

and

$$\begin{array}{cc}\hfill & \left|J\left(x\right)\right|\le \left|\tau \left({\overline{x}}_j\right)h^{-\alpha }\left[(1-\alpha ){\sigma }^{-\alpha }\right]\right|\hfill \\ \hfill +& \left|{\tau }^{\prime }\left({\overline{x}}_j\right)h^{1-\alpha }\left[-\frac{(1-\alpha )}{2}{\sigma }^{-\alpha }+{\sigma }^{1-\alpha }\right]\right|\hfill \\ \hfill +& \left|{\tau }^{″}\left({\overline{x}}_j\right)h^{2-\alpha }\left[\frac{1-\alpha }{8}{\sigma }^{-\alpha }-\frac{1}{2}{\sigma }^{1-\alpha }+\frac{1}{2-\alpha }{\sigma }^{2-\alpha }\right]\right|\hfill \\ \hfill +& \left|{\tau }^{‴}\left({\overline{x}}_j\right)h^{3-\alpha }\left[\frac{\alpha -1}{48}{\sigma }^{-\alpha }+\frac{{\sigma }^{1-\alpha }}{8}+\frac{{\sigma }^{2-\alpha }}{2(2-\alpha )}+\frac{{\sigma }^{3-\alpha }}{(2-\alpha )(3-\alpha )}\right]\right|\hfill \\ \hfill +& \left|{(\widehat{{\zeta }_j}-{\overline{x}}_j)}^4{h}^{-\alpha }\left[\frac{u^{\left(4\right)}\left({\tilde{\zeta }}_j\right)}{24}(1-\alpha ){\sigma }^{-\alpha }\right]\right|\hfill \\ \hfill \le & C{h}^{4-\alpha }{\parallel u^{\left(4\right)}\parallel }_{\infty }.\hfill \end{array}$$

(42)

Now collecting the inequalities (40)–(42) gives the following result for the case of $p=4$.

Theorem 1.

If $u\left(x\right)\in C^p(\Lambda )$ and $p\ge 1$ an integer number, then it holds the error estimate:

$$\parallel D^{\alpha }({\Pi }_{N}u-u){\parallel }_{\infty }\le C{h}^{min\{p,4-\alpha \}}{\parallel u^{\left(p\right)}\parallel }_{\infty },$$

(43)

where C independent on N.

Proof of Theorem 1.

When $p=1$, the Taylor’s theorem gives

$$\frac{d^2}{d{x}^2}{\int }_{x_{i-1}}^{x_i}\frac{\tau \left(s\right)}{{(x-s)}^{\alpha -1}}dx=[\tau ({\overline{x}}_i)+{\tau }^{\prime }({\tilde{\zeta }}_i)({\widehat{\zeta }}_i-{\overline{x}}_i)](\alpha -1)[{(x-x_i)}^{-\alpha }-{(x-x_{i-1})}^{-\alpha }],$$

and

$$\frac{d^2}{d{x}^2}{\int }_{x_{j-1}}^x\frac{\tau (s)}{{(x-s)}^{\alpha -1}}dx=[\tau ({\overline{x}}_j)+{\tau }^{\prime }({\tilde{\zeta }}_j)({\widehat{\zeta }}_j-{\overline{x}}_j)](\alpha -1)[{(x-x_{j-1})}^{-\alpha }].$$

So, by (33), the estimate (43) follows. For the cases $p=2$ and $p=3$, the estimate (43) can be obtained by a similar argument. □

Remark 2.

For a real number $p\in (0,1)$, the numerical tests show that the estimate (43) also holds.

5. Applications to Fractional Differential Equations

In this Section, some numerical examples are presented to demonstrate the efficiency of our approximation method. The following three types of meshes are used in numerical tests:

• Uniform mesh (Mesh 1):

  $$x_j=x_L+\frac{(x_R-x_L)j}{N},j=0,1,\dots ,N.$$
• Graded mesh (Mesh 2):

  $$x_j=x_L+(x_R-x_L){\left(\frac{j}{N}\right)}^q,q>1,j=0,1,\dots ,N.$$
  Note: For the two-sided operator, two-sided graded mesh will be used with an even number N:

  $$\begin{array}{ccc}\hfill x_j& =& x_L+\frac{(x_R-x_L)}{2}{\left(\frac{j}{N_h}\right)}^{q_1},j=0,1,\dots ,N_h\hfill \\ \hfill x_j& =& x_R-\frac{(x_R-x_L)}{2}{\left(\frac{N-j}{N_h}\right)}^{q_2},j=N_h+1,N_h+2,\dots N,\hfill \end{array}$$

  where $N_h=\frac{N}{2}$ and when $q=q_1=q_2$, the two-sided mesh is symmetric.
• Geometric mesh (Mesh 3):

  $$x_0=x_L,x_j=x_L+(x_R-x_L)\ast q^{N-j},0<q<1,j=1,2,\dots ,N.$$

5.1. Fractional Helmholtz Equations

To measure the accuracy of the HCSCM when the exact solution is known, we define the errors by

$$E_0=\underset{x\in \{x_1,x_2,\dots ,x_{N-1}\}}{max}\left\{\right|u_N\left(x\right)-u\left(x\right)\left|\right\},$$

where $u_N\left(x\right)$ and $u\left(x\right)$ are numerical and exact solution respectively. Let $\Lambda :=(x_L,x_R)$ and $1<\alpha \left(x\right)<2$. In this subsection we apply the HCSCM to the following variable-order fractional Helmholtz equation with homogeneous boundary conditions

$${\lambda }^{2}u\left(x\right)-{\mathcal{D}}^{\alpha \left(x\right)}u\left(x\right)=f\left(x\right),x\in \Lambda ,u\left(x_L\right)=u\left(x_R\right)=0.$$

(44)

The HCSCM for (44) is to find $u_N\in {\mathbb{V}}_N$, such that

$${\lambda }^2{u}_N\left(x\right)-{\mathcal{D}}^{\alpha \left(x\right)}{u}_N\left(x\right)=f\left(x\right),x\in {\mathbf{x}}^c,u_N\left(x_L\right)=u_N\left(x_R\right)=0.$$

(45)

The above equation leads to the following linear system:

$$\left({\lambda }^2\mathbf{M}-{\mathbf{D}}^{\alpha }\right)\mathbf{u}=\mathbf{f}$$

(46)

where $\mathbf{M}=[{\varphi }_0\left({\mathbf{x}}^c\right),{\phi }_1\left({\mathbf{x}}^c\right),{\varphi }_1\left({\mathbf{x}}^c\right),\dots ,{\phi }_{N-1}\left({\mathbf{x}}^c\right),{\varphi }_{N-1}\left({\mathbf{x}}^c\right),{\varphi }_N\left({\mathbf{x}}^c\right)]$ is the collocation matrix, ${\mathbf{x}}^{\mathbf{c}}$ and $\mathbf{u}$ as in (26), $\mathbf{f}=f\left({\mathbf{x}}^{\mathbf{c}}\right)$ and ${\mathbf{D}}^{\alpha }$ is the fractional differentiation matrix with respect to the fractional operator ${\mathcal{D}}^{\alpha \left(x\right)}$ as in (27).

Example 1.

Our first test of HCSCM is to consider the problem (44) with exact solution $u\left(x\right)=sin\left(\pi x\right)$ at $[x_L,x_R]=[-1,1]$.

The right-hand side function $f\left(x\right)={\lambda }^{2}u\left(x\right)-{\mathcal{D}}^{\alpha \left(x\right)}u\left(x\right)$ in which the fractional derivative term is approximated by when ${\mathcal{D}}^{\alpha \left(x\right)}={{}_{RL}D}_{-1,x}^{\alpha \left(x\right)}$

$${{}_{RL}D}_{-1,x}^{\alpha \left(x\right)}sin\left(\pi x\right)=\sum _{k=0}^L{(-1)}^{k+1}\frac{{\pi }^{2k+1}{(x+1)}^{2k+1-\alpha \left(x\right)}}{\Gamma (2k+2-\alpha (x\left)\right)},$$

and when ${\mathcal{D}}^{\alpha \left(x\right)}={{}_{RL}D}_{x,1}^{\alpha \left(x\right)}$

$${{}_{RL}D}_{x,1}^{\alpha \left(x\right)}sin\left(\pi x\right)=\sum _{k=0}^L{(-1)}^k\frac{{\pi }^{2k+1}{(1-x)}^{2k+1-\alpha \left(x\right)}}{\Gamma (2k+2-\alpha (x\left)\right)},$$

with $L=50$, respectively.

For $\alpha \left(x\right)$, we consider the following two cases:

• The constant-order $\alpha =1.1,1.2,1.4,1.5,1.6,1.8,1.9.$
• The variable-order $\alpha \left(x\right)=1.1+\frac{x+1}{2.5}.$

The aim of this example is to test the accuracy of the proposed method for the smooth solution. In this example, the uniform mesh is used. The error $E_0$ and the orders of convergence are listed in

Table 1. Error $E_0$ and the order of convergence (OC), for Example 1 with Mesh 1: $\alpha \left(x\right)=1.2,1.4,1.6,1.8,({\sigma }_1,{\sigma }_2)=(0.2,0.8),\lambda =0$.

N	$\mathit{\alpha }\left(\mathit{x}\right)=1.2$	OC	$\mathit{\alpha }\left(\mathit{x}\right)=1.4$	OC	$\mathit{\alpha }\left(\mathit{x}\right)=1.6$	OC	$\mathit{\alpha }\left(\mathit{x}\right)=1.8$	OC
20	1.1797 $\times {10}^{-4}$	-	1.1326 $\times {10}^{-4}$	-	1.8116 $\times {10}^{-4}$	-	2.9010 $\times {10}^{-4}$	-
40	1.6776 $\times {10}^{-5}$	2.81	1.7641 $\times {10}^{-5}$	2.68	3.3277 $\times {10}^{-5}$	2.45	6.2240 $\times {10}^{-5}$	2.22
80	2.4257 $\times {10}^{-6}$	2.79	2.8890 $\times {10}^{-6}$	2.61	6.2306 $\times {10}^{-6}$	2.42	1.3406 $\times {10}^{-5}$	2.22
120	7.8179 $\times {10}^{-7}$	2.79	1.0058 $\times {10}^{-6}$	2.60	2.3442 $\times {10}^{-6}$	2.41	5.4714 $\times {10}^{-6}$	2.21
160	3.5026 $\times {10}^{-7}$	2.79	4.7588 $\times {10}^{-7}$	2.60	1.1740 $\times {10}^{-6}$	2.40	2.8984 $\times {10}^{-6}$	2.21
200	1.8588 $\times {10}^{-7}$	2.84	2.6617 $\times {10}^{-7}$	2.60	6.8516 $\times {10}^{-7}$	2.41	1.7712 $\times {10}^{-6}$	2.21
240	1.1199 $\times {10}^{-7}$	2.78	1.6566 $\times {10}^{-7}$	2.60	4.4209 $\times {10}^{-7}$	2.40	1.1876 $\times {10}^{-6}$	2.19

. It is shown that the order of convergence of the approximation is $4-\alpha$.

The error $E_0$ and CPU time for $\alpha =1.4$ are listed in

Table 2. Error $E_0$ and the CPU time for Example 1 with Mesh 1: $\alpha =1.4,({\sigma }_1,{\sigma }_2)=(0.2,0.8),\lambda =0$.

N	${\mathit{E}}_0$	CPU Time (s)
10	3.6097 $\times {10}^{-3}$	0.018
50	8.8401 $\times {10}^{-5}$	0.165
100	1.5063 $\times {10}^{-5}$	0.479
150	5.2966 $\times {10}^{-6}$	1.046
200	2.5185 $\times {10}^{-6}$	1.831
250	1.4119 $\times {10}^{-6}$	2.721
300	8.8569 $\times {10}^{-7}$	3.397
500	2.3094 $\times {10}^{-7}$	7.363
1000	1.0710 $\times {10}^{-7}$	23.029

. Similar results can be obtained for other cases. All the computations are performed by Matlab R2020a on pc with AMD PRO A10-8770 R7, 10 COMPUTE CORES 4C+6G 3.50GHz. The Matlab route inv is used to solve the linear system (46) in our numerical tests. Other faster solver such as LU decomposition, iteration-type methods and so forth might be used to improve the efficiency.

The error $E_0$ is given in Figure 1 and Figure 2. In Figure 1, it is clearly shown that the orders of convergence of approximation is about $4-\alpha$ which confirms the estimate in Theorem 1. In Figure 2, the orders of convergence of approximation is about

$$4-\underset{-1\le x\le 1}{max}\alpha \left(x\right)$$

for the variable-order case.

Example 2.

Our second test of HCSCM is to consider the problem (44) with an exact solution that has low regularity.

When we take $u\left(x\right)=(1-x){(1+x)}^{\alpha \left(x\right)}$ with $x_L=-1,x_R=1$, the right-hand functions are same for left Riemann-Liouville, left Caputo and left Patie-Simon cases, that is,

$$f\left(x\right)={\lambda }^{2}u\left(x\right)-2\Gamma (1+\alpha \left(x\right))+\Gamma (2+\alpha \left(x\right))(1+x).$$

In fact, here we have $u\left(x\right)\in C^{\alpha }(\Lambda )$.

The error $E_0$ and the orders of convergence are listed in

Table 3. Error $E_0$ and the order of convergence (OC) for Example 2 with Mesh 1: $\alpha \left(x\right)=1.2,1.4,1.6,1.8,({\sigma }_1,{\sigma }_2)=(0.2,0.8),\lambda =0$.

N	$\mathit{\alpha }\left(\mathit{x}\right)=1.2$	OC	$\mathit{\alpha }\left(\mathit{x}\right)=1.4$	OC	$\mathit{\alpha }\left(\mathit{x}\right)=1.6$	OC	$\mathit{\alpha }\left(\mathit{x}\right)=1.8$	OC
20	3.6965 $\times {10}^{-3}$	-	3.5347 $\times {10}^{-3}$	-	1.9174 $\times {10}^{-3}$	-	5.9525 $\times {10}^{-4}$	-
40	1.6497 $\times {10}^{-3}$	1.16	1.3360 $\times {10}^{-3}$	1.40	6.3033 $\times {10}^{-4}$	1.60	1.6984 $\times {10}^{-4}$	1.81
80	7.2079 $\times {10}^{-4}$	1.19	5.0527 $\times {10}^{-4}$	1.40	2.0733 $\times {10}^{-4}$	1.60	4.8498 $\times {10}^{-5}$	1.81
120	4.4317 $\times {10}^{-4}$	1.20	2.8620 $\times {10}^{-4}$	1.40	1.0824 $\times {10}^{-4}$	1.60	2.3318 $\times {10}^{-5}$	1.81
160	3.1379 $\times {10}^{-4}$	1.20	1.9124 $\times {10}^{-4}$	1.40	6.8268 $\times {10}^{-5}$	1.60	1.3874 $\times {10}^{-5}$	1.80
200	2.4007 $\times {10}^{-4}$	1.20	1.3990 $\times {10}^{-4}$	1.40	4.7751 $\times {10}^{-5}$	1.60	9.2759 $\times {10}^{-6}$	1.80
240	1.9289 $\times {10}^{-4}$	1.20	1.0836 $\times {10}^{-4}$	1.40	3.5659 $\times {10}^{-5}$	1.60	6.6766 $\times {10}^{-6}$	1.80

. It is shown that the order of convergence of the approximation is $\alpha$.

The error $E_0$ for uniform mesh and for $\alpha =1.1,1.5,1.9$ are shown in Figure 3. It is clear that the order of convergence of $E_0$ is $\alpha$.

The error $E_0$ for $\alpha =1.2$ with three types of mesh are shown in Figure 4. It is shown that the uniform mesh achieves an $\alpha$ order of convergence of $E_0$ and the graded mesh improves significantly the order of convergence. We can also observe that the geometric mesh might achieve “higher accuracy” (see the dotted line with squares Figure 3), although the precisions are damaged for large N. The errors $E_0$ for the Caputo case and Patie-Simon case are plotted in Figure 5 and Figure 6.

If we take $u\left(x\right)=(1-x){(1+x)}^{\alpha \left(x\right)-1}$, it means that $u\left(x\right)\in C^{\alpha -1}(\Lambda )$, which has very low regularity with $x_L=-1,x_R=1$ and then the right-hand function

$$f\left(x\right)={\lambda }^{2}u\left(x\right)+\Gamma (1+\alpha \left(x\right))$$

for left Riemann-Liouville and left Patie-Simon cases (but not for left Caputo case).

The error $E_0$ and the orders of convergence are listed in

Table 4. Error $E_0$ and the order of convergence (OC) for Example 2 with Mesh 1: $\alpha \left(x\right)=1.2,1.4,1.6,1.8,({\sigma }_1,{\sigma }_2)=(0.2,0.8),\lambda =0$.

N	$\mathit{\alpha }\left(\mathit{x}\right)=1.2$	OC	$\mathit{\alpha }\left(\mathit{x}\right)=1.4$	OC	$\mathit{\alpha }\left(\mathit{x}\right)=1.6$	OC	$\mathit{\alpha }\left(\mathit{x}\right)=1.8$	OC
20	1.1624 $\times {10}^{+0}$	-	4.0781 $\times {10}^{-1}$	-	1.0764 $\times {10}^{-1}$	-	1.8809 $\times {10}^{-2}$	-
40	1.0453 $\times {10}^{+0}$	0.15	3.1153 $\times {10}^{-1}$	0.39	7.1834 $\times {10}^{-2}$	0.58	1.0960 $\times {10}^{-2}$	0.78
80	9.1739 $\times {10}^{-1}$	0.19	2.3697 $\times {10}^{-1}$	0.39	4.7653 $\times {10}^{-2}$	0.59	6.3385 $\times {10}^{-3}$	0.79
120	8.4738 $\times {10}^{-1}$	0.20	2.0173 $\times {10}^{-1}$	0.40	3.7430 $\times {10}^{-2}$	0.60	4.5930 $\times {10}^{-3}$	0.79
160	8.0061 $\times {10}^{-1}$	0.20	1.7991 $\times {10}^{-1}$	0.40	3.1524 $\times {10}^{-2}$	0.60	3.6528 $\times {10}^{-3}$	0.80
200	7.6601 $\times {10}^{-1}$	0.20	1.6460 $\times {10}^{-1}$	0.40	2.7588 $\times {10}^{-2}$	0.60	3.0577 $\times {10}^{-3}$	0.80
240	7.3879 $\times {10}^{-1}$	0.20	1.5306 $\times {10}^{-1}$	0.40	2.4738 $\times {10}^{-2}$	0.60	2.6439 $\times {10}^{-3}$	0.80

. It is shown that the order of convergence of the approximation is $\alpha -1$.

The errors $E_0$ are plotted in Figure 7, Figure 8, Figure 9 and Figure 10. Compared the Figure 4 with the Figure 7, we can find that the orders of convergence of $E_0$ are dropped to $\alpha -1$ for the exact solution that belongs to $C^{\alpha -1}(\Lambda )$, which agree with the results in Theorem 1. It is also observed that the orders of convergence of $E_0$ are improved by making use of the graded mesh and the geometric mesh similarly.

The HCSCM is comparable to the MDSCM [49] since both of them are applied piecewise polynomial to approximation problems. The numerical errors are compared by using the HCSCM and the MDSCM. The maximum errors with the degree of freedom are plotted for the constant-order $\alpha =1.1,1.5$ and $1.9$ and for the variable-order $\alpha \left(x\right)=1.1+(x+1)/2.5$ in Figure 9 and Figure 10. The black lines are for the MDSCM with fixed $N=3$ and the penalty parameter $\tau =$ 100,000. By the choice of $N=3$ of the MDSCM, the degree of piecewise polynomial in the HCSCM is the same as ones in the MDSCM. Both the uniform meshes are applied for two methods. It is shown that the accuracy of the HCSCM is better than those of the MDSCM [49] with h-refinement but the orders of convergence are almost same.

5.2. Fractional Burgers Equations

In this subsection, we try to solve the fractional Burgers equation as

$${\partial }_{t}u(x,t)+u(x,t){\partial }_{x}u(x,t)=ϵ{\mathbb{D}}_r^{\alpha (x,t)}u(x,t),$$

(47)

subject to homogeneous Dirichlet boundary condition and initial condition $u(x,0)=u_0\left(x\right)$, where $ϵ>0,1<\alpha (x,t)<2$ and $(x,t)\in (-1,1)\times (0,1]$.

The two-step Crank-Nicolson/leapfrog scheme is first employed for time stepping, then the HCSCM is applied to the resulting equations. Thus, the full discretization scheme reads as: for $k=1,2,\dots$,

$$\left\{\begin{array}{c}\left(\mathbf{M}-\Delta tϵ{\mathbf{D}}^{{\alpha }^{k+1}}\right){\mathbf{u}}^{k+1}=\mathbf{g},\hfill \\ \mathbf{M}{\mathbf{u}}^1=\left(\mathbf{M}+\Delta tϵ{\mathbf{D}}^{{\alpha }^0}\right){\mathbf{u}}^0-\Delta t\left(\mathbf{M}{\mathbf{u}}^0\right).\ast \left(\mathbf{S}{\mathbf{u}}^0\right)\hfill \\ \mathbf{M}{\mathbf{u}}^0=u_0\left(\mathbf{x}\right),\hfill \end{array}\right.$$

(48)

where

$$\mathbf{g}=\left(\mathbf{M}+\Delta tϵ{\mathbf{D}}^{{\alpha }^{k-1}}\right){\mathbf{u}}^{k-1}-2\Delta t\left(\mathbf{M}{\mathbf{u}}^k\right).\ast \left(\mathbf{S}{\mathbf{u}}^k\right),$$

and $\Delta t$ is the time stepsize, $\mathbf{M}$ the collocation matrix, ${\mathbf{D}}^{{\alpha }^k}$ the fractional differentiation matrix of order ${\alpha }^k=\alpha (x,k\Delta t)$ with respect to the fractional operator ${\mathcal{D}}^{\alpha (x,t)}={\mathbb{D}}_r^{\alpha (x,t)}$ as in (25), $\mathbf{S}$ the collocation first-order differentiation matrix which defines as

$$\mathbf{S}=[{\varphi }_0^{\prime }\left({\mathbf{x}}^c\right),{\phi }_1^{\prime }\left({\mathbf{x}}^c\right),{\varphi }_1^{\prime }\left({\mathbf{x}}^c\right),\dots ,{\phi }_{N-1}^{\prime }\left({\mathbf{x}}^c\right),{\varphi }_{N-1}^{\prime }\left({\mathbf{x}}^c\right),{\varphi }_N^{\prime }\left({\mathbf{x}}^c\right)],$$

and notation$.\ast$ the entry-to-entry multiplication.

Example 3.

In this example, we consider the fractional Burgers Equation (47) with the initial condition $u_0\left(x\right)=sin\left(\pi x\right)$.

Our first test is the numerical solutions to the Equation (47) of the left Riemann-Liouville fractional derivative. The following five cases of fractional order are considered in [40,49]:

• Case 1: (constant-order) $\alpha (x,t)=1.1,1.2,1.3,1.5,1.8$;
• Case 2: (monotonic increasing-order) $\alpha (x,t)=1+\frac{5+4x}{10}$;
• Case 3: (monotonic decreasing-order) $\alpha (x,t)=1+\frac{5-4x}{10}$;
• Case 4: (nonsmooth order) $\alpha (x,t)=1.1+\frac{4}{5}|sin\left(10\pi (x-t)\right)|$;
• Case 5: (nonsmooth order) $\alpha (x,t)=1.1+\frac{4\left|xt\right|}{5}$.

In Figure 11, the numerical solutions at $t=1$ is plotted for constant-order cases (Case 1). The obtained numerical result is the same as the one in Fig4.6 ([49]). We also compare some results by the multi-domain spectral collocation method(MDSCM) with those by the presented method(HCSCM) for $\alpha =1.1,1.5$ in Figure 12 and Figure 13. It is shown that by the HCSCM one get smoother numerical solution near the left boundary $x=-1$ than that by the MDSCM.

The numerical solutions at $t=1$ for variable-order cases(Case 2–5) are plotted in Figure 14, Figure 15, Figure 16 and Figure 17, which are agree with the results in [40,49].

The numerical solutions are also computed to the fractional Burgers equation with two-sided operators. The numerical solutions at $t=1$ for various $\alpha$’s and r’s are plotted in Figure 18, Figure 19, Figure 20, Figure 21, Figure 22 and Figure 23.

6. Conclusions

In this paper, a Hermite cubic spline collocation method (HCSCM) are developed for solving variable-order nonlinear fractional differential equations, which apply $C^1$-continuous nodal basis functions to approximate problem. It is verified that the order of convergence of the HCSCM is $O\left(h^{min\{4-\alpha ,p\}}\right)$, while the interpolating function belongs to $C^p(p\ge 1)$, where h is the mesh-size and $\alpha$ the order of the fractional derivative. The effectiveness of the HCSCM is demonstrated by solving fractional Helmholtz equations of constant-order and variable-order, and solving the fractional Burgers equation. The numerical fractional diffusions are compared with different senses.

The HCSCM can be applied to fractional-order differential equations on a two or three dimensional Descartes product domain by nodal basis tensor. Through adjusting the location of collocation points, the stability of the HCSCM can be observed numerically. Our future work will focus on the stability and error analysis of the HCSCM for some FDEs.