Due to the symmetry feature in nature, fractional differential equations precisely measure and describe biological and physical processes. Multiterm time-fractional order has been introduced to model complex processes in different physical phenomena. This article presents a numerical method based on the cubic B-spline finite element method for the solution of multiterm time-fractional differential equations. The temporal fractional part is defined in the Caputo sense while the B-spline finite element method is employed for space approximation. In addition, the four-point Gauss−Legendre quadrature is applied to evaluate the source term. The stability of the proposed scheme is discussed by the Von Neumann method, which verifies that the scheme is unconditionally stable. and norms are used to check the efficiency 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.

1. Introduction

Isaac Newton and Leibniz independently discovered fractional calculus in the 17th century [1]. Fractional calculus involves integration and differentiation of arbitrary orders. Mathematicians studied fraction calculus as an abstract area containing pure mathematical manipulations with limited applications before 1970. Thenceforth, fractional calculus emerged as an application in physics, dynamical systems, control engineering, economics, bioengineering, fluid mechanics, electrochemistry, continuum, and statistical mechanics [26]. In addition, applications can also be found in signal and image processing, heat and mass transfer, electromagnetics, dynamic system, fluid flow in porous and nonporous media, medicine, viscoelastic materials, and anomalous transport [7, 8]. The 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 differential equations were developed to model complex real-world phenomena that a single term could not accurately describe. For example, a two-term fractional order diffusion model was proposed for the total concentration in solute transport to explicitly distinguish the solute’s mobile and immobile status using fractional dynamics. The kinetic equation with two fractional derivatives of different orders describes subdiffusive motion in velocity fields. Multiterm fractional differential equations can model anomalous diffusion phenomena in complex systems and highly heterogeneous aquifers [1114]. 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 different applications [15, 16]. This work considers multiterm time-fractional DEs as follows:where , , and are smooth functions and is a differential operator defined as follows:

Many authors studied and solved (1) using different techniques. Daftardar and Bhalekar considered multiterm fractional diffusion-wave equations using separation of variables [18]. Luchko applied the Mittag−Leffler function and Fourier variable separation for the solution of (1) given in [17]. Li et al. used Mittag−Leffler function and the expansion method for the initial and boundary value problems [19]. The authors in [20] studied the exact solution of (1) using Luckh’s theorem, Mittag−Leffler function, and the Laplace operator. However, the exact solution of fractional differential equation is difficult to find in general. Even if the exact solutions are found, they contain special functions such as hypergeometric, H-function, Wright, and Mittag−Leffler functions. Therefore, it is necessary to develop an efficient numerical method for solving fractional DEs.

Numerous researchers solved (1) using different numerical methods. Ye et al. studied Riesz Caputo fractional DEs to prove the maximum principle and employed a predictor-corrector method with and schemes [20]. Qiao and Xu used the orthogonal spline collocation method for time-fractional diffusion equations [21]. Shiralashetti and Deshi solved (1) using the Haar wavelet collocation method [22]. Li and co-authors applied mixed finite element (FEM) and finite difference methods for time-fractional diffusion-wave equations and proved stability and convergence [23]. The authors in [24] gave two kinds of implicit difference algorithms for multiterm time-fractional wave equations with nonhomogeneous Dirichlet boundary conditions. Ren and Sun proposed a high-precision difference algorithm for one and two-dimensional time-fractional diffusion equations [25]. Dehghan and co-authors employed finite difference and Galerkin spectral methods for the solution of fractional order diffusion-wave equations [26]. Hussain and Haq solved time-fractional diffusion equations via the meshless method [27].

This study aims to introduce a computational technique based on the B-spline finite element method for the solution of a multiterm time-fractional diffusion equation. FEM is a powerful and established method used to approximate the solution of differential equations (DEs). In 1943, Courant introduced a piecewise polynomial for torsion problems [28]. In 1956, Turner et al. developed the stiffness matrices for trusses, beams, and other elements for engineering analysis of structures [29]. For the first 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 first 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].

Numerous researchers solved time-fractional differential equations using the B-spline finite element method. Esen and Tasbozan employed quadratic and cubic B-spline FEM for fractional Burger, gas dynamics, and Schrodinger’s equations [3335]. Khader and Khadijah solved the fractional Klein Gordon equation using quadratic B-spline FEM [36]. Ucar and Feng used FEM for the solution of fractional diffusion equations [37, 38]. In [39], the authors applied the Petrov Galerkin method for the time-fractional KdV equation. For more details, the reader is referred to [4045] and the references therein.

The primary objective of this study is to develop a hybrid scheme comprising of finite difference and finite element methods to find the numerical solution of (1). A cubic B-spline is used as a test and shape function in a finite element method called the Galerkin approach that generates the symmetric matrices. A Caputo fractional derivative is defined 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. The stability of the scheme is examined via the Fourier method. The distribution of the current work is as follows: The recurrence relation and mathematical formulation of cubic B-spline are presented in Section 2. The solution methodology and stability are discussed in Sections 3 and 4. The performance of the proposed method is examined via test problems in Section 5 while the concluding remarks are given in Section 6.

2. Cubic B-Spline

Let us consider , taking such that and . The B-spline of degree is denoted by and defined as follows:

The formula in (4) is called the Cox−de Boor recursion formula. The linear B-spline is defined as follows:

The definition of a cubic B-spline iswhere and . Approximating by cubic B-spline, we getwhere are unknown parameters to be determined. Using transformation in (6) takes the following form:

All splines apart from , and are zero over the interval . Using (8) in terms of (7) as follows:

The values of and its derivatives at the nodal points are given by

Using equations (10)–(12), it yields

3. Proposed Solution Methodology

This section describes the methodology of the proposed scheme for solving the time-fractional differential equations.

3.1. Time Discretization

This part introduces the formulation of temporal-fractional derivative of differential equations, which is given in (1). For this, let the time interval be discretized using mesh points where is the temporal mesh size and is the maximum time limit. The time-fractional derivative can be approximated by an approximation formula [46] as follows:

Using forward difference for time derivative leads to the following:where and .

The same procedure is applied for the discretization of fractional derivative of order .

3.2. Numerical Integration

Numerical integration plays an important role in applied sciences and engineering. The 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. The quadrature rule yields an exact result of degree . The source term is approximated as follows:where are weight factors and are the quadrature or sampling points. The sampling points can be obtained from the roots of Legendre polynomial and weight points from the integration of Lagrange polynomial. The four-point Gauss quadrature of sampling and weight points are

Let and substitute in (16); it becomes

3.3. Cubic B-Spline Galerkin’s Approximation

Applying Galerkin’s approach to (1) with weight function leads to weak formulation as follows:

For , (20) becomes

Let , and using (6) and (9) in (21), one has

Using (8) in (22) leads to the following:

The matrix form of (23) leads toWhere are the unknowns and , , and are element matrices and can be expressed as follows:

The source term is evaluated by Gauss−Legendre quadrature in Section 3.2 as follows:

Assembling all element contributions give the following system of equations:

The matrices , , and are septa-diagonal matrices of the form

Using (14) in (28), it becomeswhere

Applying -weighted scheme [49] and using (29) in (27), it yieldswhere

The scheme in (31) is derived for three-term time-fractional DEs. Substituting in (20), one can derive the same scheme for one-term and two-term time-fractional DEs. The initial vector is determined from (9)–(13). Using (31) together with boundary conditions, the unknown coefficients can be obtained for arbitrary time levels and hence the solution from (9).

4. Stability Analysis

The stability analysis is based on the Von Neumann concept in which the growth factor of a typical Fourier mode is defined aswhere is the mode number and is the element size. One can express (31) in terms of nodal parameters as

Using (33) in (34), we get the following:

Substituting Euler’s formula in (36), it becomes

The values of , and arewhere

Simplifying (37) leads to

Since it is clear that for all values of , , and , the amplification factor . The scheme is unconditionally stable according to the Fourier stability method. The stability of the scheme is verified in example 3 computationally using the Lax−Richtmyer stability criterion.

5. Examples

In this part, the proposed scheme is applied to solve multiterm time-fractional differential equations. The computed solutions are compared with the exact and available techniques in the literature. and norms are used to check the efficiency and accuracy of the proposed method as follows:where is the space step size. The convergence rate of the numerical results can also be calculated using the following formula:

5.1. Example 1

We consider single-term time-fractional differential equations [24].

with initial and boundary conditions

The exact solution

The method is discussed in Section 3 for multiterm time-fractional differential equations that have been used to find a solution of the problem. The values of the parameters are , , and . The solution is computed in the spatial domain [0, 1]. The results have been computed for different step sizes and are mentioned in Table 1. The error norms are computed to check the efficiency of the scheme. Table 1 clearly shows that the error decreases when the mesh size decreases, which verifies the temporal and spatial convergence. The results are compared with those of the finite difference method in [24]. The 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. The surface plots of approximate, exact, and absolute errors at and are shown in Figure 1. The graphs clearly show that as the values of and increase, the solution also increases exponentially. From the figures, it is obvious that the approximate solution coincides with the exact validity of the proposed scheme.

5.2. Example 2

We consider two-term time-fractional differential equation [14].

The initial and boundary conditions are

The exact solution is

The proposed method is applied to find the results of the problem as discussed in Section 3 with domain [0, 1]. Taking and for different values of are given in Table 2. norms are listed in Table 2 to check the performance and accuracy of the proposed technique. As the values of decrease, accuracy increases for different values of . One can see that the present results are better than those of [14]. The surface plots of exact, approximate, and error are shown in Figure 2 at and . The graph of exact and approximate solutions matches each other, showing good agreement between approximate and exact solutions.

5.3. Example 3

We consider two-term time-fractional differential equation as follows:

The initial and boundary conditions are

The exact solution is given by

The proposed method is implemented with domain to two-term time-fractional DEs. The values of and the results are obtained for different values of which are given in Table 3. and norms are calculated and presented in Table 3 to check the performance and accuracy of the technique. Moreover, the spectral radius of matrix is computed to check the stability of the proposed scheme. It has been noticed that is less than or equal to unity, which verifies the stability of the scheme Figure 3. Furthermore, the convergence rate is given in Table 3 and the order is less than unity. The stability of the graph is illustrated in Figure 4, which clearly shows that the data of spectral radius lies in [0, 1]. The surface plots of exact, approximate, and error at and are displayed in Figure 3. The exact and approximate graphs match each other, showing good agreement between approximate and exact solutions.

5.4. Example 4

We consider three-term time-fractional differential equation [50].

The initial and boundary conditions are

The exact solution is

In example 4, the same methodology has been used with domain [0, 1]. Here the values of , and . 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 finite difference scheme [50]. In addition, the convergence rate is calculated and displayed in Table 4, which is approximate to second order. The graphical solutions of exact, approximate, and error at are displayed in Figure 5. The exact and approximate plots coincide, indicating they are in good agreement.

5.5. Example 5

We consider three-term time-fractional differential equation [50].

The initial and boundary conditions are

The exact solution

The solution is computed over the domain [0, 1] and the values of parameters are , , , and . The error norms are computed for different 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 finite difference method [50]. Moreover, the convergence rate is presented in Table 5, which is greater than or equal to order 2. The graphical solutions of approximate, exact, and error at , , and are illustrated in Figure 6. The 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 figures that the approximate and exact solutions are in good agreement with each other.

6. Conclusion

This work proposed a numerical method based on cubic B-spline FEM for approximating multiterm time-fractional differential equations. The Caputo formula combined with a finite difference 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 differential equations have been solved to demonstrate the performance and accuracy of the proposed method. The efficiency and reliability of the suggested method were examined using and norms. All the obtained results showed good agreement with the exact solutions. The 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

The data used to support the findings of this study are available from the first authors upon request.


A preprint has previously been published [16].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Shams ul Arifeen conceptualized the study, designed the methodology, performed the software analysis, and prepared the original draft. Sirajul Haq supervised the study. Sirajul Haq and Farhan Golkarmanesh reviewed and edited the manuscript and performed the visualization. All the authors read and approved the final manuscript.