An exponential B-spline collocation method for fractional sub-diffusion equation

In this article, we propose an exponential B-spline collocation method to approximate the solution of the fractional sub-diffusion equation of Caputo type. The present method is generated by use of the Gorenflo-Mainardi-Moretti-Paradisi (GMMP) scheme in time and an efficient exponential B-spline based method in space. The unique solvability is rigorously discussed. Its stability is well illustrated via a procedure closely resembling the classic von Neumann approach. The resulting algebraic system is tri-diagonal that can rapidly be solved by the known algebraic solver with low cost and storage. A series of numerical examples are finally carried out and by contrast to the other algorithms available in the literature, numerical results confirm the validity and superiority of our method.


Introduction
The basic concept of anomalous diffusion dates back to Richardson's treatise on atmospheric diffusion in 1926 [31]. It has increasingly got recognition since the late 1960s within transport theory. In contrast to a typical diffusion, such process no longer follows Gaussian statistics, then the classic Fick's law fails to apply. Its most striking character is the temporal power-law pattern dependence of the mean squared displacement [22], i.e., χ 2 (t) ∼ κt α , for sub-diffusion, α < 1, while α > 1 for super-diffusion. Anomalous transport behavior is ubiquitous in physical scenarios and due to its universal mutuality, formidable challenges are introduced. In recent decades, fractional partial differential equations (PDEs) enter public vision that compare favorably with the usual models to characterize such transport motions in heterogeneous aquifer and the medium with fractal geometry [1,26]. An explosive interest has been gained among academic circles to scramble to investigate the theoretical properties, analytic techniques, and numerical algorithms for fractional PDEs [2,5,7,18,21,23,29,36].
As a model problem of the class of fractional PDEs described above, the fractional sub-diffusion equation is considered here subjected to the initial and boundary conditions as u(a, t) = g 1 (t), u(b, t) = g 2 (t), 0 < t ≤ T, (1.3) where 0 < α < 1, κ is the positive viscosity constant, and ϕ(x), g 1 (t), g 2 (t) are the prescribed functions with sufficient smoothness. In Eq. (1.1), the timefractional derivative is defined in Caputo sense, i.e., with the Gamma function Γ(·). There have already been some works dedicated to develop numerical algorithms to solve Eqs. (1.1)-(1.3) apart from a few analytic methods that are not always available for general situations. Zhang and Liu derived an implicit difference scheme and proved that it is unconditional stable [37]. Yuste and Acedo studied an explicit difference scheme based on Grünwald-Letnikov formula [34]. Along the same line, a group of weighted average difference schemes was then obtained [33]. In [4], Cui raised a high-order compact difference scheme and its convergence was detailedly discussed; another similar approach was the compact scheme stated in [30], for the fractional sub-diffusion equation with Neumann boundary condition. In [15], an effective spectral method was constructed by using the common L 1 -formula in time and a Legendre spectral approximation in space. Later, this method was extended to the time-space case [14]. The finite element method was considered by Jiang and Ma [10]. The semi-discrete lump finite element method was studied by Jin et al. for a time-fractional model with a nonsmooth right-hand side [11]. Liu et al. described an implicit RBF meshless approach for the time-fractional diffusion equation [16]. Li et al. suggested an adomian decomposition algorithm for the equations of the same type [13]. In [9], the authors solved such a model by the direct discontinuous Galerkin method with the Caputo derivative discretized by a GMMP scheme. Recently, Luo et al. established a quadratic spline collocation method for the fractional sub-diffusion equation [17], where the convergence under L ∞ -norm was analyzed. Sayevand et al. conducted a cubic B-spline collocation method [32], whose stability was provided as well. In [27], a Sinc-Haar collocation method was proposed, which used the Haar operational matrix to convert the original problem into linear algebraic equations.
In the present work, regarding the current interest in efficient numerical algorithms for fractional PDEs, we showcase a collocation method based on exponential B-spline trial function to solve Eqs. (1.1)-(1.3). The Caputo derivative is tackled by GMMP formula and the spatial derivative is approximated in an exponential spline space via a uniform nodal collocation strategy. A von Neumann like procedure leads to its unconditional stability. Its codes are tested on five numerical examples and studied in contrast with the other algorithms. The obtained method is highly accurate and calls for a lower cost to implement. This may make sense to treat the equations as the model we consider here with a long time range. The outline is as follows. In Section 2, we give a concise description of the exponential B-spline trial basis, which will be useful hereinafter. In Section 3, we construct a fully discrete exponential B-spline method on uniform meshes to discretize the model and prove that it is stable. The initial vector is addressed in Section 4, which we require to start our method. To evaluate its accuracy and advantages, numerical examples are covered in Section 5.

Description of exponential spline functions
In the sequel, let a = x 0 < x 1 < x 2 < · · · < x M −1 < x M = b be an equidistant spatial mesh on the interval [a, b], and for M ∈ N + , denote where p j is the value of function p(x) at mesh knot x j . The exponential splines are a kind of piecewise non-polynomial functions that are known as a generalization of the semi-classical cubic splines. They are recognized as a continuum of interpolants ranging from the cubic splines to the linear cases [19]. Also, like the polynomial splines, a basis of exponential B-splines is admitted and an advisable definition is the one introduced by McCartin [20], each of which is support on finite subsegments. On the above mesh together with another six knots x j , j = −3, −2, −1, M + 1, M + 2, M + 3 beyond [a, b], the mentioned exponential B-splines B j (x), j = −1, 0, . . . , M + 1, are given as follows The values of B j (x) at each knot are given as The values of B j (x) and B j (x) at each knot are given as The set of B j (x) ∈ C 2 (R), j = −1, 0, . . . , M + 1, are linearly independent and form an exponential spline space on [a, b]. The non-negative free p is termed "tension" parameter and p → 0 yields cubic spline whereas p → ∞ corresponds to the linear spline. The cubic spline interpolation causes extraneous inflexion points while the exponential spline interpolation allows to remedy this issue.

GMMP scheme for Caputo derivative
To start with, we recall the Caputo and Riemann-Liouville fractional derivatives. Given a smooth enough f (x, t), the α-th Caputo derivative is defined by and the α-th Riemann-Liouville type derivative is defined by where, m − 1 < α < m, m ∈ N is not less than 1. In common sense, (3.7) owns merits in handling the initial-valued problems, and thereby is utilized in time in most instances.

A fully discrete exponential B-spline based scheme
Define with the unknown weights {α j (t)} M +1 j=−1 yet to be determined by some certain restrictions. Discretizing Eq. (1.1) by using (3.12) in time, we have Let α n j = α j (t n ). On replacing u(x, t) by u N (x, t) and imposing the following collocation and boundary conditions at each nodal point x j , j = 0, 1, . . . , M , we obtain 14) and the boundary sets where m = 0, 1, . . . , n − 1. As a result, using Eqs. (3.15)-(3.16) to remove the unknown variables α n −1 , α n M +1 in Eq.(3.14) when j = 0, M , the above system admits a linear system of algebraic equations of size (M + 1) × (M + 1), as below where in which, m = 0, 1, . . . , n, and d n 0 , d n M are as follows The weights α n depends on α n−k , k = 0, 1, . . . , n, at its previous time levels and is found via a recursive style; once α n is obtained, α n −1 , α n M +1 are obvious due to Eqs. (3.15)-(3.16). On the other side, A is a (M + 1) × (M + 1) tri-diagonal matrix, therefore the system can be performed by the well-known Thomas algorithm, which simply needs the arithmetic operation cost O(M + 1).

Initial state
In order to start Eq. (3.17), an appropriate initial vector α 0 to the system is required. To this end, we employ the initial conditions u N (x j , 0) = ϕ(x j ), j = 0, 1, · · · , M, together with the collocation constraints got via Eq. (1.2) explicitly to determine a unique initial vector α 0 by with the notations In the same fashion, K is a (M + 1) × (M + 1) tri-diagonal matrix, so the solution of Eq. (4.18) can also be computed by Thomas algorithm.

Stability and solvability
In this section, our objective is to prove that Eqs. (3.17)-(4.18) are uniquely solvable and unconditionally stable. Ifα n j , n ≥ 1, is a perturbed solution of Eq. (3.14), we shall study how the perturbation ρ n j = α n j −α n j , which solves the homogeneous counterpart of the equation by evolves over time, where Z 0 j , Z n−k j are the quantities like P 0 j , P n−k j with regard to the perturbation. Since the classic von Neumann method does not work for Eq. (5.19), a fractional procedure is employed to analyze its stability. This extension was recently laid down in [35] applied to discuss a non-uniform implicit difference scheme for fractional diffusion equations.  ) . Then, the lemma is ascribed to s − ph < phc − s. Using the following Taylor's expansions

Proof. In virtue of A, A , one gets
and A − 2 A > 0, which implies A is strictly diagonally dominant, so is K. The stable analysis is proceeded as following.

Numerical experiments
In this part, the proposed exponential B-spline collocation method is tested on a couple of numerical examples, which suffice to gauge its accuracy and realistic performance. The computed errors are measured by L 2 -and L ∞ -norms, i.e., and for every concrete problem, the tension parameter p is optimally selected. The resulting algebraic equations are handled by Thomas algorithm and the numerical results may be compared with the other existent methods.

Conclusion
In this research, an efficient exponential B-spline based collocation method is proposed to tackle the diffusion equation with a time-fractional derivative in Caputo sense discretized by a GMMP scheme. It leads to a linear system of algebraic equations with tri-diagonal coefficient matrix and thereby can be solved speedily by Thomas algorithm. The solvability is strictly evaluated and the stable analysis is proceeded by adopting a fractional von Neumann procedure. Its codes are tested on several given models and the numerical results validate that this method is capable of dealing with these equations very well. The comparison with the other methods manifests its practicability and advantages. Moreover, the derived method is easy and economical to carry out, so it can be served as an alternative choice to model the other fractional problems.