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

In this paper, we develop a Hermite cubic spline collocation method (HCSCM) for solving variable-order nonlinear fractional differential equations, which apply C1-continuous nodal basis functions to an approximate problem. We also verify that the order of convergence of the HCSCM is about O(hmin{4−α,p}) while the interpolating function belongs to Cp(p ≥ 1), where h is the mesh size and α 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 RiemannLiouville, Caputo and Patie-Simon sense as well as two-sided cases.


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.
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 operatorrelated 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 for any x ∈ (c, b), where Here, the RL D α a,x denotes the left Riemann-Liouville fractional derivative of α-order; we will give its definition in the following section. The above lemma shows that lim x→c RL D α a,x u(x) = ∞ if the conditions u (c + ) = u (c − ) and α > 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.

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 α(x) ≡ α is a constant in the formulas.
and the right Riemann-Liouville fractional integral of order α(x) > 0 is defined as where Γ(·) is the Euler's gamma function.
and the right Riemann-Liouville fractional derivative of order α(x) > 0 is defined as where n is the positive integer such that n − 1 < α(x) < n.
and the right Caputo fractional derivative of order α(x) > 0 is defined as where n is the positive integer such that n − 1 < α(x) < n.
The well-known relationship between Riemann-Liouville and the Caputo derivative is as follows: x,x R f (x) and C D α(x) x,x R f (x) exist, then and RL D α(x) x, x, Besides the common-used definitions above, the fractional diffusion operators which limit the order 1 < α(x) ≤ 2 are also considered. A definition was proposed by Patie and Simon in [59] as follows.
and the right Patie-Simon (or mixed Caputo) fractional derivative of order 1 < α(x) < 2 is defined as x, From the above definitions and Lemma 2, hold the following relationships: and RL D α(x) x, and PS D α(x) and PS D α(x) x, Proof. Since Then note that 1 < α(x) < 2, The equality (11) is obtained by dividing factor Γ(2 − α(x)). 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 ∈ [x L , x R ] and Riemann-Liouville fractional calculus and for the Caputo fractional derivative, and for the Patie-Simon fractional derivative of 1 < α(x) < 2, The following operators with top-tilde are useful in HCSCM for x > x R , and for x < x L , and for x > x R , and for Operators CD x L ,x R * are defined similarly.
and when x ∈ [x L , x c ), we have Proof. Since Then the first equality in (22) is obtained by dividing factor Γ(α(x)) and the definitions (1) and (18). Other results can be derived by a similar argument.
If D * α(x) and D α(x) * represent all left-sided and right-sided definitions of the abovementioned, respectively, then the two-sided fractional derivative can be written as

Hermite Cubic Spline Collocation Method(HCSCM)
In the Section, the HCSCM is presented. The key role of HCSCM is the collocation fractional differentiation matrix.

Fractional Differentiation Matrix (FDM) for HCSCM
Let Λ := (x L , x R ), the first step is to divide the interval Λ into N elements, that is, which is defined by the following set of nodal basis functions. It contains 2N + 2 functions as follows. The first two functions as otherwise.
For i = 1, 2, ..., N − 1, otherwise. and the last two functions as Therefore, If u N ∈ V N , then can be expanded as In fact, a choice of this points is the Gauss-type quadrature nodes, 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 α(x)th-order differential operators defined in Section 2, we denote by D α the collocation fractional differentiation matrix, which satis- 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 and D α with Dirichlet boundary conditions is is one of the fractional differential operators defined in Section 2. Typically, the matrix 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 D α 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.

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 Then we have the following relations: here D α(x) be one of the seven: For the tilde operators, the following relations also hold: hereD α(x) can be one of the six: RLD a,b * . The nodal basis functions presented in Section 3, are the so-called shape functions after being transferred by (28), that is, and y ∈ [−1, 1]. The Hermite Spline collocation method will perform all the operators mentioned above on the shape functions (31).

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 conve-nience of analysis, denote D α = RL D α(x) x L ,x and let h i = x i − x i−1 = h, σ i,1 = σ 1 , σ i,2 = σ 2 , i = 1, 2, ..., N. Then x i = x 0 + ih, i = 0, 1, ..., N and the collocation points Let Π N : C 1 (Λ) → V N be the piecewise Hermite cubic interpolation operator, determined uniquely by For a function u(x) ∈ C(Λ), the maximum norm is defined by The following results are related to the interpolation errors [62].
where C is a constant number which do not dependent on N.
If u(x) ∈ C p (Λ)(p ≥ 1), the interpolation error holds (see [63]): where s = min{p, 4}. In the following, the error bound is presented for Assume that x ∈ [x j−1 , x j ] for some j, then the above integration can be split as Under the assumption of u(x) ∈ C 4 (Λ), from Taylor's theorem we have whereζ i ,ζ i,k ∈ (x i−1 , x i ). For the second integral in (34), we have Let σh = x − x j−1 . Note that σ ∈ (0, 1), it is easy to know that Hence, by Lemma 6, for k = 0, 1, 2, 3 and every i < j we have and and Now collecting the inequalities (40)- (42) gives the following result for the case of p = 4. Theorem 1. If u(x) ∈ C p (Λ) and p ≥ 1 an integer number, then it holds the error estimate: where C independent on N.
Proof of Theorem 1. When p = 1, the Taylor's theorem gives 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 ∈ (0, 1), the numerical tests show that the estimate (43) also holds.

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): • Graded mesh (Mesh 2): Note: For the two-sided operator, two-sided graded mesh will be used with an even number N: where N h = N 2 and when q = q 1 = q 2 , the two-sided mesh is symmetric. • Geometric mesh(Mesh 3):

Fractional Helmholtz Equations
To measure the accuracy of the HCSCM when the exact solution is known, we define the errors by where u N (x) and u(x) are numerical and exact solution respectively. Let Λ := (x L , x R ) and 1 < α(x) < 2. In this subsection we apply the HCSCM to the following variable-order fractional Helmholtz equation with homogeneous boundary conditions The HCSCM for (44) is to find u N ∈ V N , such that The above equation leads to the following linear system: is the collocation matrix, x c and u as in (26), f = f (x c ) and D α is the fractional differentiation matrix with respect to the fractional operator D α(x) as in (27).

2.
The variable-order α(x) = 1.1 + 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. It is shown that the order of convergence of the approximation is 4 − α. The error E 0 and CPU time for α = 1.4 are listed in Table 2. 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 Figures 1 and 2. In Figure 1, it is clearly shown that the orders of convergence of approximation is about 4 − α which confirms the estimate in Theorem 1. In Figure 2, the orders of convergence of approximation is about 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.
The error E 0 and the orders of convergence are listed in Table 3. It is shown that the order of convergence of the approximation is α. The error E 0 for uniform mesh and for α = 1.1, 1.5, 1.9 are shown in Figure 3. It is clear that the order of convergence of E 0 is α.
The error E 0 for α = 1.2 with three types of mesh are shown in Figure 4. It is shown that the uniform mesh achieves an α 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      If we take u(x) = (1 − x)(1 + x) α(x)−1 , it means that u(x) ∈ C α−1 (Λ), which has very low regularity with x L = −1, x R = 1 and then the right-hand function 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. It is shown that the order of convergence of the approximation is α − 1. The errors E 0 are plotted in Figures 7-10. Compared the Figure 4 with the Figure 7, we can find that the orders of convergence of E 0 are dropped to α − 1 for the exact solution that belongs to C α−1 (Λ), 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 α = 1.1, 1.5 and 1.9 and for the variable-order α(x) = 1.1 + (x + 1)/2.5 in Figures 9 and 10. The black lines are for the MDSCM with fixed N = 3 and the penalty parameter τ = 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.
3 ) of left Riemann-Liouville with variable-order case.

Fractional Burgers Equations
In this subsection, we try to solve the fractional Burgers equation as subject to homogeneous Dirichlet boundary condition and initial condition u(x, 0) = u 0 (x), where > 0, 1 < α(x, t) < 2 and (x, t) ∈ (−1, 1) × (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, ..., where and ∆t is the time stepsize, M the collocation matrix, D α k the fractional differentiation matrix of order α k = α(x, k∆t) with respect to the fractional operator D α(x,t) = D α(x,t) r as in (25), S the collocation first-order differentiation matrix which defines as and notation . * the entry-to-entry multiplication.
Example 3. In this example, we consider the fractional Burgers Equation (47) with the initial condition u 0 (x) = sin(πx).
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 5:(nonsmooth order) α(x, t) = 1.1 + 4|xt| 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 α = 1.1, 1.5 in Figures 12 and 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 Figures 14-17, which are agree with the results in [40,49].      3 ) for all cases.  3 ) for all cases.

Conclusions
In this paper, a Hermite cubic spline collocation method (HCSCM) are developed for solving variable-order nonlinear fractional differential equations, which apply C 1continuous nodal basis functions to approximate problem. It is verified that the order of convergence of the HCSCM is O(h min{4−α,p} ), while the interpolating function belongs to C p (p ≥ 1), where h is the mesh-size and α 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.