Comparative Numerical Study of Spline-Based Numerical Techniques for Time Fractional Cattaneo Equation in the Sense of Caputo–Fabrizio

: This study focuses on numerically addressing the time fractional Cattaneo equation involving Caputo–Fabrizio derivative using spline-based numerical techniques. The splines used are the cubic B-splines, trigonometric cubic B-splines and extended cubic B-splines. The space derivative is approximated using B-splines basis functions, Caputo–Fabrizio derivative is discretized, using a ﬁnite difference approach. The techniques are also put through a stability analysis to verify that the errors do not pile up. The proposed scheme’s convergence analysis is also explored. The key advantage of the schemes is that the approximation solution is produced as a smooth piecewise continuous function, allowing us to approximate a solution at any place in the domain of interest. A numerical study is performed using various splines, and the outcomes are compared to demonstrate the efﬁciency of the proposed schemes.


Introduction
The time fractional Cattaneo differential equation (TFCDE) under consideration is [1] ∂v(s, t) ∂t and the boundary conditions, where (s, t) ∈ ∆ = [0, L] × [0, T], 1 < α < 2, g ∈ C[0, T], and f 1 (t), f 2 (t), φ(s), ψ(s) are known functions. Moreover, CF a D α t v(s, t) is the Caputo-Fabrizio derivative given by For mathematical modeling of real-world problems, fractional differential equations are often used. Scientists in a variety of fields are pushed to improve the interpretations of their findings by utilizing the fractional order derivatives, which are particularly useful. In mathematical modeling of many scientific situations, fractional order differential equations provide more accurate information than regular differential equations. Fractional derivatives are used to describe a variety of physical phenomena [2]. This is owing to the fact that fractional operators assess both global and local properties when analyzing system evolution. In addition, integer-order calculus can sometimes contradict the experimental results; therefore, non-integer order derivatives may be preferable [3]. It is difficult to determine the solution to fractional differential equations (FDEs). As a result, a numerical method must be used to obtain the solution to these partial differential equations. To tackle these problems numerically, many approaches have been developed and extended. The existence of solution of FDEs can be seen in [4]. Diethelm et al. presented the predictor-corrector method [5] for the numerical solution of FDEs. Meerschaert and Tadjern [6] developed a finite difference method for a fractional advection-dispersion equation. The homotopy analysis method [7] for the fractional initial value problem was developed by Hashim et al. An eigenvector expansion method for motion containing fractional derivatives was presented by Suarez and Shokooh [8].
When compared to the finite difference approach, other spectral methods, such as the operational matrix method, are particularly popular since they provide good accuracy and take less time to compute. This method works well with fractional ordinary differential equations (ODEs), fractional partial differential equations (PDEs), and variable order PDEs. Jafari et al. [9] gave applications of Legendre wavelets in solving FDEs numerically. The Haar wavelet operational matrix of fractional order integration and its applications in solving fractional order differential equations can be seen in [10]. Chebyshev wavelets [11] were used by Yuanlu for solving a nonlinear fractional order differential equation. Li and Sun [12] developed a generalized block pulse operational matrix method for the solution of FDEs. Obidat [13] used Legendre polynomials to approximate the solution of nonlinear FDEs. Genocchi polynomials [14] were used by Araci to find numerical solutions of FDEs. Grbz and Sezer [15] solved a class of initial and boundary value problems arising in science and engineering using Laguerre polynomials. Caputo and Fabrizio proposed one of the most recent fractional order derivatives. For more applications of this new derivative and the related work, the reader is referred to [1,[16][17][18][19][20][21][22][23][24][25][26][27][28][29].
In comparison to polynomials, the B-splines based collocation methods provide a good approximation rate, are computationally quick, numerically consistent, and have secondorder continuity. To obtain numerical solutions to differential equations, multiple numerical approaches based on various forms of B-splines functions were recently utilized. Inspired by the popularity of spline approaches in finding numerical solutions of fractional partial differential equations, various splines-based numerical techniques have been developed for the numerical solution of the Cattaneo equation involving the Caputo-Fabrizio derivative. The main motivation behind this work is that to the authors' knowledge, this equation has not been solved using the B-splines basis functions. In the current work, B-splines are used to approximate the space derivative, while the Caputo-Fabrizio derivative is approximated using finite differences. Moreover, the presented schemes are tested for stability and convergence analysis.

Numerical Schemes
In this section, the cubic B-splines, extended cubic B-splines and the trigonometric cubic B-splines are used to develop numerical techniques for the numerical solution of time fractional Cattaneo equation (TFCE) (1).

Numerical Scheme Based on Cubic B-Splines
Let τ = T N and h = L M be the step length in space and time direction, respectively. Set t m = mτ, s j = jh, where the positive integers, N and M, are used. The knots s j divide the solution domain ∆ equally into M equal subintervals [s j , s j+1 ], j = 0, 1, . . . , M − 1, where a = s 0 < s 1 < · · · < s M = b. The approximate solution V(s, t) to the exact solution v(s, t) in the following form is acquired by our scheme for solving (1) where C j (t) are unknowns to be found, and B j (s) [30] are cubic B-splines basis (CuBS) functions given by Here, B j−1 (s), B j (s) and B j+1 (s) are survived due to the local support characteristic of the cubic B-splines so that the approximation v m j at the grid point (s j , t m ) at the mth time level is given as The time-dependent unknowns C m j (t) are found using the specified initial and boundary conditions as well as the collocation conditions on B j (s). As a result, the approximation v m j and its required derivatives are where a 1 = 1 6 , a 2 = 4 6 , b 1 = 1 2h , c 1 = 1 h 2 , and c 2 = − 2 h 2 . Let g = {g m : 0 ≤ m ≤ N} be the collection of grid functions on a uniform mesh of the interval [0, T] such that where, and |R Now, we employ the Caputo-Fabrizio fractional derivative and CuBS to establish the numerical scheme for solving (1). Using CuBS and the approximation given in (8), we obtain where R m+1 = O(τ 2 + h 2 ). Thus, by ignoring R m+1 and using the discretization (v m Rearranging the above equation, we obtain where σ = (α − 1 + M 0 τ ) and µ = (α − 1)τ. Using the CuBS approximation (7) in (11), we obtain where η 1 = σa 1 − µc 1 , η 2 = σa 2 − µc 2 , η 3 = σa 1 ,and η 4 = σa 2 . In matrix notation, the above equation is expressed as where the matrices A 1 , A 2 , B 1 , Ψ and G are By combining Equations (12) and (13), we have (M + 3) × (M + 3), a system of linear equations which can be solved uniquely.

Initial State
First of all, it is essential to find the initial vector to initiate the iteration procedure. This vector is obtained from initial conditions as Thus, (M + 3) × (M + 3) a system of linear equations results, and this system can be written in matrix notation as where the matrices A 3 , C 0 and B 2 are and

Numerical Scheme Based on Extended Cubic B-Splines
A cubic B-spline of degree four with a free parameter η is called an extended cubic B-spline. This kind of cubic B-spline was introduced by Han and Liu in 2003. We follow the same notations for the time and space discretizations that we used before. The extended cubic B-spline (ECuBS) basis functions, B 4 j (s, η) are given by where η ∈ [−8, 1]. Here, B 4 j−1 (s), B 4 j (s) and B 4 j+1 (s) are survived due to local support characteristic of the cubic B-splines so that the approximation v m j at the grid point (s j , t m ) at mth time level is given as The time-dependent unknowns C m j (t) are found using the specified initial and boundary conditions as well as the collocation conditions on B j (s). As a result, the approximation v m j and its required derivatives are where ω 1 = 4−η 24 , ω 2 = 8+η 12 , ω 3 = 1 2h , ω 4 = 0, ω 5 = 2+η 2h 2 and ω 6 = − 2+η 2h 2 . By following the same procedure as was done for cubic B-splines and using the ECuBS approximation given in (16), we obtain the following approximation to the solution of (1) where, η 5 = σω 1 − µω 5 , η 6 = σω 2 − µω 5 , η 7 = σω 1 and η 8 = σω 2 . In matrix notation, the above Equation (17) is expressed as where the matrices A 4 , A 5 and B 3 are The above system gives (M + 1) equations in (M + 3) unknowns. For a unique solution, two additional linear equations are necessary. From the boundary conditions, we obtain the required equations as follows , By combining Equations (17) and (18), we have (M + 3) × (M + 3), a system of linear equations, which can be solved uniquely.

Numerical Scheme Based on Trigonometric Cubic B-Splines
We follow the same notations for the time and space discretizations used before. The trigonometric cubic B-spline (TCuBS) basis functions are given by [31] where l(s j ) = sin( j−1 (s), TB 4 j (s) and TB 4 j+1 (s) are survived due to the local support characteristic of the trigonometric cubic B-splines so that the approximation v m j at the grid point (s j , t m ) at mth time level is given as The time-dependent unknowns C m j (t) are found using the specified initial and boundary conditions as well as the collocation conditions on B j (s). As a result, the approximation v m j and its required derivatives are where 2+4 cos(h) . By following the same procedure as was done for cubic B-splines and using the approximation (8) in (1), we obtain where σ = (α − 1 + M 0 τ ) and µ = (α − 1)τ. Using the CuTBS approximation given in (21), we obtain the following approximation to the solution of (1) where η 9 = σζ 1 − µζ 5 , η 10 = σζ 2 − µζ 6 , η 11 = σζ 1 , and η 12 = σζ 2 . In matrix notation, (23) is expressed as where the matrices A 7 , A 8 and B 5 are By combining Equations (23) and (24), we have (M + 3) × (M + 3), a system of linear equations, which can be solved uniquely.

Stability Analysis
This section deals with stability analysis of the scheme based on cubic B-splines. The stability analysis of the schemes based on extended and cubic trigonometric B-splines can be carried out by a similar argument. We use the Fourier method to study the stability analysis of the scheme. LetṼ 0 be the perturbation vector of initial values V 0 andṼ m , 1 ≤ m ≤ N − 1 be the approximate solution of the scheme (12). The error vector δ m is defined as where, Define the grid functions as follows: We can expand δ m (s) into Fourier series as where, and using the Parseval's equality, We can expand δ m j into Fourier series, and because the difference equations are linear, we can analyze the behavior of total error by tracking the behavior of an arbitrary nth component. Based on the above analysis, we can suppose that the solution of (11) has the following form where σ s = 2πl L , I = √ −1. Substituting the above expression into (11), we obtain Using the CuBS approximation given in (7) and Equation (26)  + a 2 + a 1 exp(Iσ s h)), ⇒ d m+1 (σ(a 2 + a 1 (2 cos(σ s h)) − µ(c 2 + c 1 (2 cos(σ s h))) = m ∑ l=1 (M m−l − M m−l+1 )δ t d l (a 2 + a 1 (2 cos(σ s h)) + σd m (a 2 + a 1 (2 cos(σ s h)), which, on further simplification, reduces to where, r = ( c 2 +2c 1 cos(σ s h) a 2 +2a 1 cos(σ s h) ). Proof. We use the mathematical induction for proof. For m = 1, we have from (28), Then, by using Lemmas 1 and 2, we obtain This completes the proof. (12) is unconditionally stable for α ∈ (1, 2).

Theorem 2. The scheme
Proof. By using Theorem 1, Parseval's equality and mτ ≤ T, we obtain which means that the scheme is unconditionally stable.

Convergence Analysis
The convergence of the scheme based on cubic B-splines is presented in this section. The convergence analysis of the extended and cubic trigonometric B-splines based numerical scheme follows accordingly. Let e m From Equation (11) and R m+1 j = O(τ 2 + h 2 ) and noting that e 0 j = 0, we have Define the functions We expand the above functions into Fourier series expansions as where, Applying Parseval's equalities, to the above expression, we have Now, we suppose that e m j = ξ m exp Iσ s jh , R m j = λ m exp Iσ s jh , where σ s = 2πl L . Substituting relations (31) in Equation (29).
The above expression is further simplified as Theorem 3. Let ξ m be the solution of (32), then, there is a positive constant C such that Proof. We use the mathematical induction to prove this claim. For m = 1, we have from (32) Now by using the convergence of the series on the RHS of (30), we know that there exists a constant C 2 such that Theorem 4. The scheme (12) is convergent, and the order of convergence is O(τ 2 + h 2 ).
Proof. By Theorem 3, Equation (30) and mτ ≤ T, we have This completes the proof.

Numerical Findings and Discussion
The efficiency and the validity of the suggested methodologies are confirmed in this part using various test problems by utilizing the L 2 and L ∞ error norms. The numerical results obtained by the proposed schemes are compared. Mathematica 12 was used to obtain the numerical and graphical results.

CuBS TCuBS ECuBS
The corresponding source term is g(s, t) = (1 + t) sin s. The analytic solution of the given problem is v(s, t) = t sin s. In order to achieve the desired numerical results the presented schemes are applied on Example 2. The errors obtained by the schemes are compared with each other in Tables 4-6. For various time stages, a sharp contrast between the exact and approximate solutions is presented in Figure 4. The 2D absolute error profile is plotted in Figure 5. Figure 6 depicts a 3D comparison between the exact and approximate solutions.
The approximate solution using cubic B-spline-based scheme when τ = 0.01 and M = 20 at T = 0.5 and T = 1 for Example 2 are given by respectively.
The corresponding source term is g(s, t) = (1 + t)(1 − s) cos s − 2t sin s. The analytic solution of the given problem is v(s, t) = t(1 − s) cos s. The proposed methodologies are utilized to acquire the numerical results for Example 3.
A comparison of computed errors is provided in Tables 7-9. For various time stages, a close comparison between the exact and approximate solutions is displayed in Figure 7. The 2D error function is plotted in Figure 8. Figure 9 depicts a 3D comparison between the exact and approximate solutions.

Concluding Remarks
The spline-based collocation schemes are developed for the numerical solution of the time fractional Cattaneo differential equation involving the Caputo-Fabrizio time fractional derivative. To begin with, the space derivative involved is approximated using the cubic B-spline. Secondly, using finite differences, the Caputo-Fabrizio derivative is approximated. The stability and convergence analysis of the schemes are also discussed in detail. The splines used are the cubic B-splines, extended cubic B-splines and the trigonometric cubic B-splines. The key advantage is that the approximate solution is obtained as a piecewise continuous function so that approximate solution at any desired position in the domain can be tracked. The efficiency and accuracy of the proposed approaches are confirmed by the experimental findings. The suggested schemes can be applied to a wide range of problems in varied fields of applied sciences.