Computational Study of Multiterm Time-Fractional Differential Equation Using Cubic B-Spline Finite Element Method

Due to the symmetry feature in nature, fractional diferential equations precisely measure and describe biological and physical processes. Multiterm time-fractional order has been introduced to model complex processes in diferent physical phenomena. Tis article presents a numerical method based on the cubic B-spline fnite element method for the solution of multiterm time-fractional diferential equations. Te temporal fractional part is defned in the Caputo sense while the B-spline fnite element method is employed for space approximation. In addition, the four-point Gauss − Legendre quadrature is applied to evaluate the source term. Te stability of the proposed scheme is discussed by the Von Neumann method, which verifes that the scheme is unconditionally stable. L 2 and L ∞ norms are used to check the efciency and accuracy of the proposed scheme. Computed results are compared with the exact and available methods in the literature, which show the betterment of the proposed method.


Introduction
Isaac Newton and Leibniz independently discovered fractional calculus in the 17th century [1]. Fractional calculus involves integration and diferentiation of arbitrary orders. Mathematicians studied fraction calculus as an abstract area containing pure mathematical manipulations with limited applications before 1970. Tenceforth, fractional calculus emerged as an application in physics, dynamical systems, control engineering, economics, bioengineering, fuid mechanics, electrochemistry, continuum, and statistical mechanics [2][3][4][5][6]. In addition, applications can also be found in signal and image processing, heat and mass transfer, electromagnetics, dynamic system, fuid fow in porous and nonporous media, medicine, viscoelastic materials, and anomalous transport [7,8]. Te noteworthy property of fractional order is nonlocality, which provides an excellent description of the memory and hereditary properties of various physical processes [9,10].
Multiterm time-fractional diferential equations were developed to model complex real-world phenomena that a single term could not accurately describe. For example, a twoterm fractional order difusion model was proposed for the total concentration in solute transport to explicitly distinguish the solute's mobile and immobile status using fractional dynamics. Te kinetic equation with two fractional derivatives of diferent orders describes subdifusive motion in velocity felds. Multiterm fractional diferential equations can model anomalous difusion phenomena in complex systems and highly heterogeneous aquifers [11][12][13][14]. Furthermore, the proposed model is used to discretize distributed-order derivatives in DEs. Hence, studies of multiterm fractional order DEs have become essential and valuable for diferent applications [15,16]. Tis work considers multiterm time-fractional DEs as follows: where h 1 (w), v 1 (t), and v 2 (t) are smooth functions and P α 1 ,α 2 ,...,α m (z t ) is a diferential operator defned as follows: Many authors studied and solved (1) using diferent techniques. Daftardar and Bhalekar considered multiterm fractional difusion-wave equations using separation of variables [18]. Luchko applied the Mittag− Lefer function and Fourier variable separation for the solution of (1) given in [17]. Li et al. used Mittag− Lefer function and the expansion method for the initial and boundary value problems [19]. Te authors in [20] studied the exact solution of (1) using Luckh's theorem, Mittag− Lefer function, and the Laplace operator. However, the exact solution of fractional diferential equation is difcult to fnd in general. Even if the exact solutions are found, they contain special functions such as hypergeometric, H-function, Wright, and Mittag− Lefer functions. Terefore, it is necessary to develop an efcient numerical method for solving fractional DEs.
Numerous researchers solved (1) using diferent numerical methods. Ye et al. studied Riesz Caputo fractional DEs to prove the maximum principle and employed a predictor-corrector method with L 1 and L 2 schemes [20]. Qiao and Xu used the orthogonal spline collocation method for time-fractional difusion equations [21]. Shiralashetti and Deshi solved (1) using the Haar wavelet collocation method [22]. Li and co-authors applied mixed fnite element (FEM) and fnite diference methods for time-fractional difusion-wave equations and proved stability and convergence [23]. Te authors in [24] gave two kinds of implicit diference algorithms for multiterm time-fractional wave equations with nonhomogeneous Dirichlet boundary conditions. Ren and Sun proposed a high-precision diference algorithm for one and two-dimensional time-fractional difusion equations [25]. Dehghan and co-authors employed fnite diference and Galerkin spectral methods for the solution of fractional order difusion-wave equations [26].
Tis study aims to introduce a computational technique based on the B-spline fnite element method for the solution of a multiterm time-fractional difusion equation. FEM is a powerful and established method used to approximate the solution of diferential equations (DEs). In 1943, Courant introduced a piecewise polynomial for torsion problems [28]. In 1956, Turner et al. developed the stifness matrices for trusses, beams, and other elements for engineering analysis of structures [29]. For the frst time in 1960, the terminology "Finite Element Method" was used by Clough in his paper on plane elasticity [30]. After the 1960s, FEM gained popularity and was used by mathematicians and engineers to solve complex problems. Besides the FEM, splines are employed for numerical solutions of DEs and TFDEs with piecewise polynomials. Isaac Schoenberg frst introduced the concept of splines in 1946. Carl de Boor worked with Schoenberg and introduced a recursive formula for splines. A notable property of B-spline is its local compact support system that generates sparse matrices; that is why the researchers primarily use B-spline as a basis function in FEM [31,32].

Complexity
Te primary objective of this study is to develop a hybrid scheme comprising of fnite diference and fnite element methods to fnd the numerical solution of (1). A cubic B-spline is used as a test and shape function in a fnite element method called the Galerkin approach that generates the symmetric matrices. A Caputo fractional derivative is defned for the temporal fractional part, and the four-point Gauss-Legendre quadrature is used to evaluate the numerical integration of functions to obtain better accuracy. Te stability of the scheme is examined via the Fourier method. Te distribution of the current work is as follows: Te recurrence relation and mathematical formulation of cubic B-spline are presented in Section 2. Te solution methodology and stability are discussed in Sections 3 and 4. Te performance of the proposed method is examined via test problems in Section 5 while the concluding remarks are given in Section 6.

Cubic B-Spline
Let us consider Ω � [c, d], taking w l ∈ [c, d] such that c � w 0 < w 1 < w 2 · · · w N � d and h � w l+1 − w l . Te B-spline of degree r is denoted by Ψl, r(w) and defned as follows: Te formula in (4) is called the Cox− de Boor recursion formula. Te linear B-spline is defned as follows: Te defnition of a cubic B-spline is where C l are unknown parameters to be determined. Using (6) takes the following form: All splines apart from Ψ l− 1 , Ψ l , Ψ l+1 , and Ψ l+2 are zero over the interval [w l , w l+1 ]. Using (8) in terms of (7) as follows: Te values of Ψ k (w) and its derivatives at the nodal points are given by Using equations (10)- (12), it yields

Proposed Solution Methodology
Tis section describes the methodology of the proposed scheme for solving the time-fractional diferential equations.

Time Discretization.
Tis part introduces the formulation of temporal-fractional derivative of diferential equations, which is given in (1). For this, let the time interval [0, T] be discretized using mesh points t s � sΔt, s � 0, 1, 2, . . . , n where Δt � T/n is the temporal mesh size and T is the maximum time limit. Te timefractional derivative can be approximated by an L 1 approximation formula [46] as follows: Using forward diference for time derivative leads to the following: where Te same procedure is applied for the discretization of fractional derivative of order α i , i � 1, 2, 3.

Numerical Integration.
Numerical integration plays an important role in applied sciences and engineering. Te numerical integration is primarily evaluated using Newton cotes and Gauss quadrature [47,48]. Newton cotes quadrature such as Trapezoidal, Simpson's rule, and their composite forms are appropriate for uniform nodes. Gaussian quadrature is preferable for nonuniform nodes to obtain accuracy even for fewer nodal points. Te quadrature rule yields an exact result of degree 2n − 1. Te source term is approximated as follows: where T j are weight factors and y j are the quadrature or sampling points. Te sampling points can be obtained from the roots of Legendre polynomial P n (y) � 1/(2 n n!)(d n (y 2 − 1) 2 /(dy n )) and weight points from the integration of Lagrange polynomial. Te four-point Gauss quadrature of sampling and weight points are and, Let y � (h/2) + (h/2)ξ and substitute in (16); it becomes

Complexity
Te scheme in (31) is derived for three-term timefractional DEs. Substituting m � 1, 2 in (20), one can derive the same scheme for one-term and two-term time-fractional DEs. Te initial vector C 0 is determined from (9)- (13). Using (31) together with boundary conditions, the unknown coefcients can be obtained for arbitrary time levels and hence the solution from (9).
Substituting Euler's formula e ikh � cos kh + i sin kh in (36), it becomes  Table 3: Errors norm of Y at diferent step sizes.

Complexity
Te values of ζ, σ and μ are Simplifying (37) leads to Since it is clear that for all values of a, b, and β, the amplifcation factor |Λ| ≤ 1. Te scheme is unconditionally stable according to the Fourier stability method. Te stability of the scheme is verifed in example 3 computationally using the Lax− Richtmyer stability criterion.

Examples
In this part, the proposed scheme is applied to solve multiterm time-fractional diferential equations. Te computed solutions are compared with the exact and available techniques in the literature. L ∞ and L 2 norms are used to check the efciency and accuracy of the proposed method as follows: where h is the space step size. Te convergence rate of the numerical results can also be calculated using the following formula: Convergence Rate � C.R, with initial and boundary conditions Te exact solution Te method is discussed in Section 3 for multiterm timefractional diferential equations that have been used to fnd a solution of the problem. Te values of the parameters are α 1 � 1, α 2 � 0.85, θ � 0.5, and a i � 1, i � 1, 2. Te solution is computed in the spatial domain [0, 1]. Te results have been computed for diferent step sizes and are mentioned in Table 1. Te error norms are computed to check the efciency of the scheme. Table 1 clearly shows that the error decreases when the mesh size decreases, which verifes the temporal and spatial convergence. Te results are compared with those of the fnite diference method in [24]. Te table shows that the results obtained using the present method are better. In addition, the convergence rate is calculated and presented in Table 1, which is closest to order 2. Te surface plots of approximate, exact, and absolute errors at Δt � h � 0.1 and t � 0.1 are shown in Figure 1. Te graphs clearly show that as the values of w and t increase, the solution also increases exponentially. From the fgures, it is obvious that the approximate solution coincides with the exact validity of the proposed scheme.
Te initial and boundary conditions are Te exact solution is Te proposed method is applied to fnd the results of the problem as discussed in Section 3 with domain [0, 1]. Taking a i � 1, i � 1, 2, θ � 0.5 and for diferent values of α i are given in Table 2. L 2 norms are listed in Table 2 to check the performance and accuracy of the proposed technique. As the values of Δt decrease, accuracy increases for diferent values of α i . One can see that the present results are better than those of [14]. Te surface plots of exact, approximate, and error are shown in Figure 2 at h � Δt � 0.1 and t � 1. Te graph of exact and approximate solutions matches each other, showing good agreement between approximate and exact solutions.

Example 3.
We consider two-term time-fractional differential equation as follows: Te initial and boundary conditions are Te exact solution is given by Te proposed method is implemented with domain [0, 1] to two-term time-fractional DEs. Te values of a i � 1, i � 1, 2, θ � 0.5 and the results are obtained for diferent values of α i which are given in Table 3. L 2 and L ∞ norms are calculated and presented in Table 3 to check the performance and accuracy of the technique. Moreover, the spectral radius of matrix M is computed to check the stability of the proposed scheme. It has been noticed that ρ(M) is less than or equal to unity, which verifes the stability of the scheme Figure 3. Furthermore, the convergence rate is given in Table 3 and the order is less than unity. Te stability of the graph is illustrated in Figure 4, which clearly shows that the data of spectral radius lies in [0, 1]. Te surface plots of exact, approximate, and error at Δt � h � 0.05 and t � 1are displayed in Figure 3. Te exact and approximate graphs match each other, showing good agreement between approximate and exact solutions.
Te initial and boundary conditions are Te exact solution is In example 4, the same methodology has been used with domain [0, 1]. Here the values of a 1 � a 2 � 1, a 3 � 0.5, α 1 � 0.7, α 2 � 0.6, and α 3 � 0.4. L 2 norms are given in Table 4 to check the performance and accuracy of the proposed technique. It can be seen that as the time step decreases, the accuracy increases. Table 4 shows that the present method gives better results than the implicit fnite diference scheme [50]. In addition, the convergence rate is calculated and displayed in Table 4, which is approximate to second order. Te graphical solutions of exact, approximate, and error at Δt � h � 0.1, t � 1 are displayed in Figure 5. Te exact and approximate plots coincide, indicating they are in good agreement.
Te initial and boundary conditions are Te exact solution Te solution is computed over the domain [0, 1] and the values of parameters are a 1 � a 3 � 1, a 2 � 0.5, θ � 0.5, α 1 � 0.8, α 2 � 0.7, and α 3 � 0.6 . Te error norms are computed for diferent step sizes and are displayed in Table 5. It is clear from the table that the error norms decrease as the time step decreases. Table 5 shows that the results obtained from the present method give better accuracy as compared to the implicit fnite diference method [50]. Moreover, the  convergence rate is presented in Table 5, which is greater than or equal to order 2. Te graphical solutions of approximate, exact, and error at t � 1, Δt � 0.1, and h � 0.1 are illustrated in Figure 6. Te solution plot for various values of α's is displayed in (6c), and one can notice that the peaks of solution increase by decreasing the values of α's. It is clear from the fgures that the approximate and exact solutions are in good agreement with each other.

Conclusion
Tis work proposed a numerical method based on cubic B-spline FEM for approximating multiterm time-fractional diferential equations. Te Caputo formula combined with a fnite diference has been employed for the temporal fractional part, and Gauss− Legendre quadrature was used to evaluate the source term. Five test problems including single, double, and three-term time-fractional diferential equations have been solved to demonstrate the performance and accuracy of the proposed method. Te efciency and reliability of the suggested method were examined using L 2 and L ∞ norms. All the obtained results showed good agreement with the exact solutions. Te stability of the scheme has been explained through Fourier analysis, and it was found that the scheme is unconditionally stable. Comparing approximate solutions with the exact and existing methods in the literature shows that the proposed method has a better outcome.

Data Availability
Te data used to support the fndings of this study are available from the frst authors upon request.

Disclosure
A preprint has previously been published [16].