Numerical Treatment of Time-Fractional Klein–Gordon Equation Using Redefined Extended Cubic B-Spline Functions

In this article we develop a numerical algorithm based on redefined extended cubic B-spline functions to explore the approximate solution of the time-fractional Klein–Gordon equation. The proposed technique employs the finite difference formulation to discretize the Caputo fractional time derivative of order α ∈ (1, 2] and uses redefined extended cubic B-spline functions to interpolate the solution curve over a spatial grid. A stability analysis of the scheme is conducted, which confirms that the errors do not amplify during execution of the numerical procedure. The derivation of a uniform convergence result reveals that the scheme is O(h2 + Δt2−α) accurate. Some computational experiments are carried out to verify the theoretical results. Numerical simulations comparing the proposed method with existing techniques demonstrate that our scheme yields superior outcomes.

The fractional KGE plays a significant role in quantum mechanics, the study of solitons, and condensed matter physics. Many approaches have been adopted to solve equations of Klein/sine-Gordon type efficiently, including the Adomian decomposition method, the variational iteration method [32][33][34], and the homotopy analysis method [35]; see also the references cited in these works. Jafari et al. proposed using fractional B-splines for approximate solution of fractional differential equations [36]. In Vong and Wang [37,38] space compact difference schemes were applied to one-and twodimensional time-fractional Klein-Gordon-type equations, and stability and convergence of the proposed numerical approaches were established with the aid of an energy method. In Dehghan et al. [39] the authors used a meshless method based on radial basis functions to develop an unconditionally stable numerical scheme for fractional Klein/sine-Gordon equations. The Adomian decomposition method and an iterative method were applied in Jafari [40] to solve Klein-Gordon-type equations involving fractional time derivatives. A fully spectral approach was employed in Chen et al. [41] that uses finite differences for time discretization and Legendre spectral approximation in the spatial direction to construct numerical solutions of non-linear partial differential equations involving fractional derivatives. A sinc-Chebyshev collocation method (SCCM) was developed in Nagy [42] for numerical treatment of the time-fractional nonlinear KGE. Recently, in Kanwal et al. [43], Genocchi polynomials were employed together with the Ritz-Galerkin scheme to solve fractional KGEs and diffusion wave equations. A linearized second-order scheme was introduced in Lyu and Vong [44] to solve non-linear time-fractional Klein-Gordon-type equations. Later on, in Doha et al. [45], a space-time spectral approximation was proposed for solving non-linear variable-order fractional Klein/sine-Gordon differential equations.
In this article we propose using redefined extended cubic Bspline (RECBS) functions for numerical solution of the timefractional KGE. RECBS functions are basically a generalization of typical cubic B-spline functions that involve a free parameter which provides the flexibility to fine-tune the solution curve. We employ the usual finite central difference approach to discretize the Caputo fractional time derivative and use RECBS functions for spatial integration.
This article is organized as follows. The Caputo definition of fractional time derivative and the finite difference formulation for temporal discretization are reviewed in section 2; this section also includes a brief introduction to extended cubic B-spline and RECBS functions and their applications to space discretization. The stability analysis of the proposed algorithm is presented in section 3, and the description of theoretical convergence is given in section 4. The approximate results are reported and discussed in section 5. Finally, concluding remarks are given in section 6.

Extended Cubic B-Spline Functions
Let the spatial domain [a, b] be partitioned into M parts of equal length h = b−a M with boundary points a = x 0 < x 1 < · · · < x M = b, where x m = x 0 + mh for m = 0 : 1 : M. For a sufficiently continuous function v(x, t), there always exists a unique extended cubic B-spline (ECBS) approximation V * (x, t): where the ξ m (t) are to be calculated and the fourth-degree ECBS blending functions S m (x, λ) are defined as [47] S m (x, λ) = 1 24h 4 Here λ, with −n(n − 2) ≤ λ ≤ 1, is a real number responsible for fine-tuning the curve, and n gives the degree of the ECBS used to generate different forms of ECBS functions. The approximate solution (V * ) r m = V * (x m , t r ) and its first two derivatives with respect to the spatial variable x at the rth time step can be expressed in terms of ξ m as [48]

Redefined Extended Cubic B-Spline Functions
In the typical ECBS collocation method, the basis functions S −1 , S 0 , . . . , S M+1 do not vanish at the boundaries of the spatial domain when Dirichlet-type end conditions are imposed. Therefore, we need to redefine them so that the resulting set of basis functions will vanish at the boundaries. For this, a weight function (x, t) is introduced to eliminate ξ −1 and ξ M+1 from Equation (9) in the following manner [49]: where the weight function (x, t) and the redefined ECBS (RECBS) functions are given by and.

Space Discretization
Using Equation (12) in Equation (8) at t = t r+1 , we obtain Discretizing at x = x j , we get Using (12), the last expression takes the form Consequently, we get the following system of M + 1 equations in where To start the numerical procedure, we use the given initial conditions to obtain the set of equations The matrix representation of (19) is  The ξ j values are then substituted into (12) to get V 0 . Now we can use (18) for r = 0, 1, 2, . . . , R − 1. However, for r = 0 the term involving V −1 appears in Equation (18). This issue is resolved by using the following substitution derived from the velocity condition given in (2):

STABILITY ANALYSIS
We use the Fourier method to study the stability of the proposed numerical method. Let ε r m andε r m denote, respectively, the exact and approximate growth factors of the Fourier modes. The error, ̺ r m , is given by where For the sake of simplicity, we shall investigate the stability of the proposed scheme with f = 0. The equation for the round-off error is derived from Equations (8) and (21) as The error equation satisfies the end conditions and We define the grid function as Frontiers in Physics | www.frontiersin.org Now, ̺ r (x) can be written in the form of a Fourier series as follows: where Taking the · 2 norm, we get From Parseval's equality we have b a |̺ r (n)| 2 dx= ∞ −∞ |ε n (m)| 2 , so the above expression can be written as Next, we consider the solution in terms of Fourier series, where ι = √ −1 and ν = 2π n b−a . Using Equation (29) in Equation (22) and then dividing by e iνkh gives (βb 1 + ρ 1 b 1 + ρb 4 )ε r+1 e −ινh + (βb 2 + ρ 1 b 2 + ρb 5 )ε r+1 We know that e iνh + e −iνh = 2 cos(νh), so after collecting like terms, the following useful relation is obtained:   Proof: For r = 0 in (31), we have Suppose that the result is true for r = 1 : 1 : R. Then, from Equation (31) we get (13) is unconditionally stable.

CONVERGENCE OF THE SCHEME
To investigate the convergence of the proposed scheme, we follow the approach in Khalid et al. [50]. Before proceeding, we state the following useful theorems [51,52].   For any knot x m , we have From (11)  Then, for x ∈ [x m−1 , x m ], S m (x, λ) and S m−1 (x, λ) are bounded above by 1 12 (8 + λ). Similarly, S m+1 (x, λ) and S m−2 (x, λ) are bounded above by 1 24 (4 − λ) where h is reasonably small and ̥ > 0 is a constant not depending on h.
Proof: LetṼ(x, t) = M m=0 d m (t)η m (x) be the calculated spline for the approximate solution V(x, t) and the exact solution v(x, t). Let The boundary conditions can be rewritten as   We define η r = max{|η r m | : 0 ≤ m ≤ M},ẽ r m = |ζ r m | and e r = max{|e r m | : 0 ≤ m ≤ M}. For r = 0, Equation (35) transforms into the following relation: Using the initial condition e 0 = 0, we obtain Taking absolute values of η r m and ζ r m and with adequately small h, we haveẽ 1 m ≤ 6̥h 4 βh 2 (λ + 2) + 12(−2 − λ)ρ + ρ 1 h 2 (2 + λ) using the boundary conditions, from which we conclude that where ̥ 1 is independent of the spatial grid spacing. Using the induction technique, we assume thatẽ k m ≤ ̥ k h 2 is true for k = 1 : 1 : r. Let ̥ = max{̥ k : 0 ≤ k ≤ r}; then Equation (35) becomes Again, taking absolute values of η r m and ζ r m , we havẽ  Hence, for all values of n,ẽ r+1 m ≤ ̥h 2 .
Taking the infinity norm and applying Lemma (3.1), we obtain Making use of the triangle inequality, we get Using the inequalities (32) and (38) in (39) Using the above theorem with expression (5), it is easy to conclude that the numerical approach converges unconditionally. Therefore, where ̥ is a constant and α ∈ (1, 2]. Hence, theoretically, the proposed scheme is O(h 2 + t 2−α ) accurate.

NUMERICAL RESULTS AND DISCUSSION
To examine the accuracy of the proposed method, we conduct a numerical study of some test problems. The L ∞ and L 2 error norms are calculated as [53] Also, the experimental order of convergence (EOC) is computed by the following important formula [54]: All numerical computations were performed using Mathematica 9.0.      where The initial/end conditions can be extracted from the analytical exact solution (1 − x) For Example 5.1, the piecewise-defined approximate solution obtained using the proposed method with α = 1.25, 0 ≤ x ≤ 1, n = 100, and t = 0.01 is given by  Table 1. It can easily be seen that our scheme is more accurate than the SCCM [42]. In Table 2 the absolute and relative numerical errors are listed for our method with M = 100, t = 0.001, and α = 1.6 at x = 0.4, 0.6, 0.8 when t = 0.4, 0.8. We can see that the computational results are superior to those obtained from the SCCM [42]. Table 3 compares the absolute errors of the proposed method, the variational iteration method (VIM) [34], and the SCCM [42] under different values of α.
where the forcing term f (x, t) on right-hand side is given by f (x, t) = 1 2 Ŵ(3 + α) sin(πx)t 2 + (1 + π 2 )t 2+α sin(πx) For Example 5.2, the piecewise-defined numerical solution obtained using the proposed method with α = 1.5, 0 ≤ x ≤ 1, n = 100, and t = 0.01 is given by The initial/boundary conditions can be extracted from the analytical exact solution v(x, t) = sin(πx)t 2+α . The absolute numerical errors at different grid points of the RECBS solution for Example 5.2 using t = 0.001 and M = 100 are listed in Table 5. Again it can be observed that our scheme is more accurate than the SCCM [42]. Table 6 reports the absolute and relative errors in our numerical computation with M = 100, t = 0.001, and α = 1.6 at x = 0.4, 0.6, 0.8 when t = 0.4, 0.8. It is clear that the results are better than those obtained by the SCCM [42]. Table 7 compares the absolute errors of the proposed method, VIM [34], and SCCM [42] under different values of α.
The EOC in the spatial direction, using t = 0.001 and α = 1.50, is tabulated in Table 8. The experimental rate of convergence of the proposed scheme is found to be in line with the theoretical prediction. Figure 5 shows the behavior at different time stages of numerical solutions obtained using α = 1.5, M = 100, and t = 0.001. The 3D plots of exact and numerical solutions with α = 1.5 and M = 100 are displayed in Figure 6. The absolute error between the exact and approximate solutions using α = 1.3, M = 100, and t = 0.001 is plotted in Figure 7.

CONCLUSION
In this work we have conducted a numerical investigation of the time-fractional Klein-Gordon equation by applying the redefined extended cubic B-spline collocation method. A finite central difference formulation is employed for temporal discretization, while a set of redefined extended cubic B-spline functions is used to interpolate the solution curve in the spatial direction. The unconditional stability of the proposed scheme is established, and the orders of convergence along the space and time grids are shown to be O(h 2 ) and O( t) 2−α , respectively. The computational outcomes of the proposed algorithm show that the order of convergence agrees with the theoretical results.
The numerical scheme has been tested on different problems, and comparison of the results reveals our method's advantage over VIM [34] and SCCM [42].