Approximation of Caputo time-fractional diffusion equation using redefined cubic exponential B-spline collocation technique

Abstract: The purpose of this work is to find the numerical solution of the Caputo timefractional diffusion equation using the modified cubic exponential B-spline (CExpB-spline) collocation technique. First, the CExpB-spline functions are modified and then used to discretize the space derivatives. Three numerical examples are considered for checking the efficiency and accuracy of the method. The obtained results are compared with those reported earlier showing that the present technique gives highly accurate results. Von Neumann stability is carried out which gives the guarantee that the technique is unconditionally stable. The rate of convergence is also obtained. Furthermore, this technique is efficient and requires less storage.


Introduction
Fractional differential equations have been applied in modelling of many fields of physics and processes such as earthquakes, optics, finance, hydrology, traffic flow, fluid mechanics, fractional kinetics, mathematical biology, measurement of visco-elastic material properties, electrical network, electro-chemistry electro-magnetic, signal processing, control theory, acoustics, material sciences [2,5,8,9,17,22,39,[41][42][43][44][45][46]. The nonlocal property is the most important advantage of these equations showing state of a complex system does not depend only on its current state but also depends on its all-previous states. Due to this reason, the fractional calculus has becoming more and more popular.
The fractional diffusion equations (FDEs) are generalized forms of classical diffusion equations. These equations have been solved by many researchers. For example, authors of [6,11,24,26,32] have been applied finite element method, compact FDM, Crank-Nicolson FDM, B-spline based method, implicit difference approximation, respectively, to solve FDEs. Khader [14] considered an efficient method based upon Chebyshev approximations. Murio [20] developed an implicit unconditionally stable method. Sun et al. [25] applied a semi-analytical FEM to solve a class of these equations. Tadjeran et al. [27] applied aforesaid method with spatial extrapolation to solve a class of variable coefficients FDEs. Lin and Xu [15] constructed a method based on FDM and Legendre spectral method while Ghanbari and Atangana [7] presented a method based on the product-integration rule for solving aforesaid equations. Ç elik and Duman [3] examined a Crank-Nicolson method with the Riesz fractional derivative. Zhai and Feng [35] introduced a block-centered FDM for solving these equations.
Additionally, Zhai et al. [34,36] constructed unconditionally stable FDMs to solve the time-space fractional diffusion and three-dimensional time-fractional subdiffusion equations, respectively. Wu and Zhai [33] proposed a high order FDM to solve the 2D time-fractional convection-dominated diffusion equation while Zhai et al. [37,38] proposed a high-order compact FDM and ADI method to solve the 3D fractional convection-diffusion equation. Furthermore, Hanert [10] proposed a flexible numerical method to discretize the space-time FDEs, where pseudo-spectral expansion has been used for discretization of the time derivative and either high order pseudo-spectral or low order FE and FD have been used for space derivative. Murillo and Yuste [19] considered explicit difference scheme, where three-point centered formula has been used to approximate the spatial derivative. Verma et al. [31] studied nonlinear diffusion equations analytically and numerically using the classical Lie symmetry method.
In this paper, we consider a Caputo time-fractional diffusion equation subject to the initial condition the Dirichlet boundary conditions where c D α t u (x, t n ) denotes the Caputo time-fractional derivative [40] as

Discretization of the problem
First, the Caputo time-fractional derivative is discretized. For this, we partition the time domain [0, T ] uniformly as 0 = t 0 < t 1 , ..., < t N = T , where interval ∆t = t n − t n−1 = T N for n = 1, 2, . . . , N. Here N is the number of time mesh. By the definition of Caputo derivative at time t = t n , the time-fractional derivative is approximated as whereb l = (l + 1) 1−α − l 1−α ∀ 0 ≤ l ≤ N and T r n ∆t gives the truncation error ≤c u ∆t 2−α , where the constantc u is only related to u(x, t). Lemma 3.1 The below relations should be hold for the coefficientsb l Proof. For the proof of Lemma 3.1, see [16].
In Eq (2.2), the free parameterκ have been used for obtaining the different forms of CExpB-spline functions. The set of E p m (x) ∀ m = −1, 0, 1, . . . , M + 1 forms a basis over the problem domain. We assume that the approximation u M to the exact u(x, t) at the point (x m , t n ) is expressed in terms of linear combinations of the CExpB-spline functions and unknown time-dependent quantities as follows: where C i (t) are the unknown quantities which we have to evaluate for the approximated solution u M (x, t) at x i , t j . Since each CExpB-spline covers four elements, so each element is covered by four CExpB-splines. So the variation of the u M (x, t), over the element, can be written as: Using Eq (2.4), the u(x, t) and its first two derivatives at the knots in terms of C n m are given as: Now, we modify the CExpB-spline basis functions which generate a new set of CExpB-spline basis functions. The procedure, for modifying the basis functions, is given as [1,12,13]: where {Ê p 0 (x),Ê p 1 (x) . . . . . .Ê p M (x)} forms a basis over the problem domain. Next, the modified form of the approximated solutions, as a linear combination of modified CExpB-spline functions, is given by Next, the initial vector C 0 m , m = 0, 1, . . . , M can be obtained by initial condition and boundary values of its derivatives as: (2.10) Eq (2.10) yields a (M + 1) × (M + 1) system as: Now, discretizing the Eq (1.1) by using (2.1) and Crank-Nicolson method, we have Using Eqs (2.4) and (2.6) in Eq (2.12) and simplifying the terms, we have where For m = 0, from Eq (2.13), we get , n = 1, 2, ..., N. (2.14) For m=M, from Eq (2.13), we get At time step t n , n = 1, 2, ..., N, the Eqs (2.13)-(2.15) can be formulated into a (M + 1) × (M + 1) linear system as follows:

Stability analysis
The von-Neumann stability analysis [4,18,21,23] have been done in this section. Taking f (x, t) = 0 in Eq (1.1) and discretizing, we get Simplifying and rearranging the terms, we get which can be written as Now we consider one Fourier mode out of the full solution C n m = δ n e kiφ as trial solutions at a given point x i , where φ = θh. The θ and h are the mode number and element size respectively, and k = √ −1. Substituting the trial solution in above equation and simplifying the terms, we get Now, using above equation and the values of A, B, D and E in Eq (3.3), we have .

(3.6)
From (3.6), we get |δ| ≤ 1 , and hence the method is unconditionally stable for the discretized system of the Caputo time-fractional diffusion equation.

Computational results and discussions
In this section, we consider numerical examples of the Caputo time-fractional diffusion equation in order to check the accuracy and efficiency of the method. We use the following error norms : The ROC is calculated by where, E h 1 and E h 2 represent the errors at space mesh sizes h 1 and h 2 , respectively and E k 1 and E k 2 represent the errors at time mesh sizes k 1 and k 2 , respectively.

Example 1
First, we consider the Caputo time-fractional diffusion equation We solve the Caputo time-fractional diffusion Eq (4.3) for different values of M and N. Table  1 shows the comparison of the present method with the finite element method [6] and cubic B-spline collocation method [24] for different values of M, fixed ∆t = 0.001 at the time-fractional order α = 0.5 in terms of L 2 and L max error norms. This table shows that the present method gives better results than the results obtained by those available in [6,24]. Moreover, one can notice from Table 1 that the error norms are decreasing as we increase the space as well as time mesh sizes. The analytical and approximate solutions together with the absolute errors are presented in Figure 1 forκ = 11, ∆x = 0.001, and α = 0.5 at ∆t = 0.2, ∆t = 0.1 , and ∆t = 0.05, respectively. Figure 2 shows the comparison between analytical and approximate solutions for M = 64, ∆t = 0.001, and α = 0.5 at different t. From these figures, we notice that there is an excellent agreement between analytical and approximate solutions. Also, the absolute errors are decreasing on increasing the time interval.

Example 2
Next, we consider the Caputo time-fractional diffusion equation with the exact solution We solve the Caputo time-fractional diffusion Eq (4.4) for various values of M. Table 2 shows the maximum absolute error norms L max for different values of time-fractional orders i. e. α = 0.2, 0.5 and 0.8 at fixed time mesh size of N = 50 for different space mesh sizes. One can see from Tables 2 and 3 that the error norms are decreasing as we increase the mesh sizes. As we can notice from this table that the present method is of order O (∆x) 2 , (∆t) 1+α . The approximate solutions for time levels t = 0.1, 0.4, 0.6 and 0.8 are plotted in Figure 3 for M = 32, N = 50 and time-fractional order α = 0.5. One can notice that the amplitude of the approximate solutions is increasing on incresing the time t.

Example 3
Now, we consider the Caputo time-fractional diffusion equation subject to the initial and boundary conditions u (x, 0) = 0, 0 < x < 1, and with the exact solution Now, we solve the Caputo time-fractional diffusion Eq (4.5) for different values of M and N. Table  4 shows the comparison in term of L max error norms obtained by the present method and the method based on the product-integration (PI) rule presented in [7] for ∆x = 1 101 , β = 6, different values of N and α at t = 1. As we can see that the present method gives more accurate results than the method presented in [7]. Also, we notice that the error norms are decreasing as we increase the time mesh sizes. Figure 4 shows the analytical and approximate solutions together with absolute error norms for ∆x = 1 101 , α = 0.85, β = 5, ∆t = 0.01 at T = 0. 35 showing that errors are decreasing on increasing the grid sizes.

Conclusions
A modified CExpB-spline collocation technique has been presented for solving the Caputo timefractional diffusion equation . The modified CExpB-spline collocation technique is used to discretize the space derivatives. The three examples of the Caputo time-fractional diffusion equation have been considered. The obtained results show that the present method gives more accurate results than the results obtained in [6,7,24]. It is observed numerically that the method is second-order accurate in space and (1 + α) order in time. The stability analysis shows that the method is unconditionally stable. Moreover, the implementation of the present method is easy and needs low memory storage which is the advantage. The present method can easily be extended for solving higher dimensional fractional PDEs.