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

Authors: Shams Ul Arifeen; Sirajul Haq; Farhan Golkarmanesh
Journal: Complexity
 (2022)

---

## Abstract

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.L2 and L∞ 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.


---

## Body

## 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 [2–6]. 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 [11–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 different applications [15, 16]. This work considers multiterm time-fractional DEs as follows:(1)Pα1,α2,…,αm∂tYw,t=∂2Yw,t∂w2+fw,t,w,t∈c,d×0,T,Yw,0=h1w,w∈c,d,Y0,t=v1t,YL,t=v2t,t>0,where h1w, v1t, and v2t are smooth functions and Pα1,α2,…,αm∂t is a differential operator defined as follows:(2)Pα1,α2,…,αm∂t=∑i=1mai∂tαi,ai>0,m∈N+,0<αm<αm−1⋯<α2<α1≤1.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 L1 and L2 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 [33–35]. 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 [40–45] 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Ω=c,d, taking wl∈c,d such that c=w0<w1<w2⋯wN=d and h=wl+1−wl. The B-spline of degree r is denoted by Ψl,rw and defined as follows:(3)Ψl,0w=1,w∈wl,wl+1,0,Otherwise,(4)Ψl,rw=w−wlwl+r−wlΨl,r−1w+wl+r+1−wwl+r+1−wl+1Ψl+1,r−1w,l=01N.The formula in (4) is called the Cox−de Boor recursion formula. The linear B-spline is defined as follows:(5)Ψl,1w=1hw−wl−1,w∈wl−1,wl,wl+1−w,w∈wl,wl+1.The definition of a cubic B-spline is(6)Ψl,3w≡Ψlw==g11hw−wl−2,w∈wl−2,wl−1,g21hw−wl−1,w∈wl−1,wl,g21hwl+1−w,w∈wl,wl+1,g11hwl+2−w,w∈wl+1,wl+2,0,Otherwise,where g1w=w3 and g2w=1+3w+3w2−3w3. Approximating Yw,t by cubic B-spline, we get(7)Yw,t=∑l=−1N+1CltΨlw,where Cl are unknown parameters to be determined. Using transformation ξ=w−wl,0≤ξ≤h in (6) takes the following form:(8)Ψl−1=1h3h−ξ3,Ψl=1h3h3+3h2h−ξ+3hh−ξ2−3h−ξ3,Ψl+1=1h3h3+3h2ξ+3hξ2−3ξ3,Ψl+2=ξh3.All splines apart fromΨl−1,Ψl,Ψl+1, and Ψl+2 are zero over the interval wl,wl+1. Using (8) in terms of (7) as follows:(9)Yξ,t=∑k=l−1l+2CktΨkξ.The values ofΨkw and its derivatives at the nodal points are given by(10)Ψkwl=4,ifk=l,1,ifk−l=1,0,ifk−l≥2,.(11)Ψk′wl=0,ifk=l,−3h,ifk=1−1,−3h,ifk=l+1,0,ifk−l≥2,(12)Ψk″wl=−12h2,ifk=l,6h2,ifk−1=1,0,ifk−1≥2.Using equations (10)–(12), it yields(13)Yl=Cl−1+4Cl+Cl+1,Yl′=3hCl+1−Cl−1,Yl″=6h2Cl+1−2Cl+Cl−1,l=01N.


## 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 0,T be discretized using mesh points ts=sΔt,s=0,1,2,…,n where Δt=T/n is the temporal mesh size and T is the maximum time limit. The time-fractional derivative can be approximated by an L1 approximation formula [46] as follows:(14)∂tαYw,tn+1=1Γ1−α∫0tn+1∂Yw,z∂ztn+1−z−αdz,=1Γ1−α∑p=0n∫tptp+1∂Yw,z∂ztn+1−z−αdz.Using forward difference for time derivative leads to the following:(15)∂tαYw,tn+1=1Γ1−α∑p=0nYp+1−YpΔt+OΔt∫tptp+1tn+1−z−αdz,=1Γ1−α∑p=0nYp+1−YpΔtn+1−p1−α−n−p1−α1−αΔt1−α+OΔt2−α,=Δt−αΓ2−α∑p=0nYn+1−p−Yn−pp+11−α−p1−α+OΔt2−α,=rα∑p=0nsαpYn+1−p−Yn−p+OΔt2−α,where rα=Δt−α/Γ2−α and sαp=p+11−α−p1−α.The same procedure is applied for the discretization of fractional derivative of orderαi,i=1,2,3.


### 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 2n−1. The source term is approximated as follows:(16)∫−11gydy≈∑j=14Tjgyj,where Tj are weight factors and yj are the quadrature or sampling points. The sampling points can be obtained from the roots of Legendre polynomial Pny=1/2nn!dny2−12/dyn and weight points from the integration of Lagrange polynomial. The four-point Gauss quadrature of sampling and weight points are(17)y=∓15±23035and,T=18±3036.Lety=h/2+h/2ξ and substitute in (16); it becomes(18)h2∫0hgh2+h2ξdξ≈h2∑j=14Tjgh2+h2ξj.


### 3.3. Cubic B-Spline Galerkin’s Approximation

Applying Galerkin’s approach to (1) with weight function W leads to weak formulation as follows:(19)∫cdWPα1,α2,…,αm∂tY−∂2Y∂w2−fdw=0,(20)∫cd∑i=1mai∂tαiWY+∂W∂w∂Y∂wdw−WYw|cd−∫cdWfdw=0.Form=3, (20) becomes(21)∫cda1∂tα1WY+a2∂tα2WY+a3∂tα3WY+∂W∂w∂Y∂wdw−WYw|cd=∫cdWfdw.LetW=Ψi,i=l−1,l,l+1,l+2, and using (6) and (9) in (21), one has(22)∑j=l−1l+2∫wlwl+1ΨiΨja1∂α1Cj∂t1α+a2∂α2Cj∂tα2+a3∂α3Cj∂tα3+Ψi′Ψj′Cjdw−ΨiΨj′Cj|wlwl+1,=∫wlwl+1Ψifdw.Using (8) in (22) leads to the following:(23)∑j=l−1l+2∫0hΨiΨjdξ∑k=13al∂αlCj∂tαl+∫0hΨi′Ψj′dξCj−ΨiΨj′Cj|0h,=∫0hΨifξ+wi,tdξ.The matrix form of (23) leads to(24)∑j=13ajA1e∂αjCe∂tαj+A2e−A3eCe=Fe.Where Cl−1,Cl,Cl+1,Cl+2T are the unknowns and A1e,A2e, A3e, and Fe are element matrices and can be expressed as follows:(25)A1e=h1402012960112911889336060933118812916012920,A2e=110h1821−36−321102−87−36−36−8710221−3−362118,A3e=3h10−104−1−411−4−140−101,Fe=∫0hΨl−1fdξ∫0hΨlfdξ∫0hΨl+1fdξ∫0hΨl+2fdξ.The source term is evaluated by Gauss−Legendre quadrature in Section3.2 as follows:(26)Fie≈h2∑j=14TjfjΨi,j,i=l−1,l,l+1,l+2.Assembling all element contributions give the following system of equations:(27)∑i=13aiA1∂αiC∂tαi+A2−A3C,=F.The matricesA1, A2, and A3 are septa-diagonal matrices of the form(28)A1:h1401,120,1191,2416,1191,120,1,A2:110h−3,−72,−45,240,−45,−72,−3,A3:3h0,0,0,0,0,0,0.Using (14) in (28), it becomes(29)∑i=13ai∂αiC∂tαi≈∑i=13airαiCn+1−rαiCn+rαi∑p=1nsαipCn+1−p−Cn−p,where(30)rαi=Δt−αiΓ2−αi,sαip==p+11−αi−p1−αi,i=1,2,3.Applyingθ-weighted scheme [49] and using (29) in (27), it yields(31)B1Cn+1=B2Cn+B3,where(32)B1=∑i=13airαiA1+θA2−θA3,B2=∑i=13airαiA1+1−θ−A2+A3,B3=∑p=1n∑i=13airαisαipCn+1−p−Cn−pA1+θFn+1+1−θFn.The scheme in (31) is derived for three-term time-fractional DEs. Substituting m=1,2 in (20), one can derive the same scheme for one-term and two-term time-fractional DEs. The initial vector C0 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).


## 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 0,T be discretized using mesh points ts=sΔt,s=0,1,2,…,n where Δt=T/n is the temporal mesh size and T is the maximum time limit. The time-fractional derivative can be approximated by an L1 approximation formula [46] as follows:(14)∂tαYw,tn+1=1Γ1−α∫0tn+1∂Yw,z∂ztn+1−z−αdz,=1Γ1−α∑p=0n∫tptp+1∂Yw,z∂ztn+1−z−αdz.Using forward difference for time derivative leads to the following:(15)∂tαYw,tn+1=1Γ1−α∑p=0nYp+1−YpΔt+OΔt∫tptp+1tn+1−z−αdz,=1Γ1−α∑p=0nYp+1−YpΔtn+1−p1−α−n−p1−α1−αΔt1−α+OΔt2−α,=Δt−αΓ2−α∑p=0nYn+1−p−Yn−pp+11−α−p1−α+OΔt2−α,=rα∑p=0nsαpYn+1−p−Yn−p+OΔt2−α,where rα=Δt−α/Γ2−α and sαp=p+11−α−p1−α.The same procedure is applied for the discretization of fractional derivative of orderαi,i=1,2,3.


## 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 2n−1. The source term is approximated as follows:(16)∫−11gydy≈∑j=14Tjgyj,where Tj are weight factors and yj are the quadrature or sampling points. The sampling points can be obtained from the roots of Legendre polynomial Pny=1/2nn!dny2−12/dyn and weight points from the integration of Lagrange polynomial. The four-point Gauss quadrature of sampling and weight points are(17)y=∓15±23035and,T=18±3036.Lety=h/2+h/2ξ and substitute in (16); it becomes(18)h2∫0hgh2+h2ξdξ≈h2∑j=14Tjgh2+h2ξj.


## 3.3. Cubic B-Spline Galerkin’s Approximation

Applying Galerkin’s approach to (1) with weight function W leads to weak formulation as follows:(19)∫cdWPα1,α2,…,αm∂tY−∂2Y∂w2−fdw=0,(20)∫cd∑i=1mai∂tαiWY+∂W∂w∂Y∂wdw−WYw|cd−∫cdWfdw=0.Form=3, (20) becomes(21)∫cda1∂tα1WY+a2∂tα2WY+a3∂tα3WY+∂W∂w∂Y∂wdw−WYw|cd=∫cdWfdw.LetW=Ψi,i=l−1,l,l+1,l+2, and using (6) and (9) in (21), one has(22)∑j=l−1l+2∫wlwl+1ΨiΨja1∂α1Cj∂t1α+a2∂α2Cj∂tα2+a3∂α3Cj∂tα3+Ψi′Ψj′Cjdw−ΨiΨj′Cj|wlwl+1,=∫wlwl+1Ψifdw.Using (8) in (22) leads to the following:(23)∑j=l−1l+2∫0hΨiΨjdξ∑k=13al∂αlCj∂tαl+∫0hΨi′Ψj′dξCj−ΨiΨj′Cj|0h,=∫0hΨifξ+wi,tdξ.The matrix form of (23) leads to(24)∑j=13ajA1e∂αjCe∂tαj+A2e−A3eCe=Fe.Where Cl−1,Cl,Cl+1,Cl+2T are the unknowns and A1e,A2e, A3e, and Fe are element matrices and can be expressed as follows:(25)A1e=h1402012960112911889336060933118812916012920,A2e=110h1821−36−321102−87−36−36−8710221−3−362118,A3e=3h10−104−1−411−4−140−101,Fe=∫0hΨl−1fdξ∫0hΨlfdξ∫0hΨl+1fdξ∫0hΨl+2fdξ.The source term is evaluated by Gauss−Legendre quadrature in Section3.2 as follows:(26)Fie≈h2∑j=14TjfjΨi,j,i=l−1,l,l+1,l+2.Assembling all element contributions give the following system of equations:(27)∑i=13aiA1∂αiC∂tαi+A2−A3C,=F.The matricesA1, A2, and A3 are septa-diagonal matrices of the form(28)A1:h1401,120,1191,2416,1191,120,1,A2:110h−3,−72,−45,240,−45,−72,−3,A3:3h0,0,0,0,0,0,0.Using (14) in (28), it becomes(29)∑i=13ai∂αiC∂tαi≈∑i=13airαiCn+1−rαiCn+rαi∑p=1nsαipCn+1−p−Cn−p,where(30)rαi=Δt−αiΓ2−αi,sαip==p+11−αi−p1−αi,i=1,2,3.Applyingθ-weighted scheme [49] and using (29) in (27), it yields(31)B1Cn+1=B2Cn+B3,where(32)B1=∑i=13airαiA1+θA2−θA3,B2=∑i=13airαiA1+1−θ−A2+A3,B3=∑p=1n∑i=13airαisαipCn+1−p−Cn−pA1+θFn+1+1−θFn.The scheme in (31) is derived for three-term time-fractional DEs. Substituting m=1,2 in (20), one can derive the same scheme for one-term and two-term time-fractional DEs. The initial vector C0 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 as(33)Cqn=Λneiqkh,where k is the mode number and h is the element size. One can express (31) in terms of nodal parameters as(34)∑p=q−3q+3ζp−q+4Cpn+1,=∑p=q−3q+3σp−q+4Cpn.Using (33) in (34), we get the following:(35)∑p=q−3q+3ζp−q+4Λn+1eipkh,=∑p=q−3q+3σp−q+4Λneipkh,(36)Λn+1eiqkhζ1e−3ikh+ζ2e−2ikh+ζ3e−ikh+ζ4+ζ5eikh+ζ6e2kh+ζ7e3ikh,=Λneiqkhσ1e−3ikh+σ2e−2ikh+σ3e−ikh+σ4+σ5eikh+σ6e2ikh+σ7e3ikh+∑p=1nΛn+1−p−Λn−peiqkhμ1e−3ikh+μ2e−2ikh+μ3e−ikh+μ4+μ5eikh+μ6e2kh+μ7e3ikh.Substituting Euler’s formulaeikh=coskh+isinkh in (36), it becomes(37)∑q=−33ζq+4Λn+1cosqkh+isinqkh,=∑q=−33σq+4Λn+μq+4∑p=1nΛn+1−p−Λn−pcosqkh+isinqkh.The values ofζ, σ and μ are(38)ζ1=s1−3s2,ζ2=120s1−72s2,ζ3=1191s1−45s2,ζ4=2416s1+240s2,ζ5=1191s1−45s2,ζ6=120s1−72s2,ζ7=s1−3s2,σ1=s1−3s2,σ2=120s1−72s2,σ3=1191s1−45s2,σ4=2416s1−240s2,σ5=1191s1−45s2,σ6=120s1−72s2,σ7=s1−3s2,μ1=s3,μ2=120s3,μ3=1191s3,μ4=2416s3,μ5=1191s3,μ6=120s3,μ7=s3,where(39)s1=∑i=13airαih140,s2=120h,s3=∑i=13airαisαiph140,β=∑p=1nΛ1−p−Λ−p.Simplifying (37) leads to(40)Λn+1=a−ba+b+aa+bβΛn,a=2416+2382coshk+240cos2hk+2cos3hk,b=240−90coshk−144cos2hk−6cos3hk.Since it is clear that for all values ofa, b, and β, the amplification factor Λ≤1. 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.L∞ and L2 norms are used to check the efficiency and accuracy of the proposed method as follows:(41)L∞=maxjY˜j−Yj,L2=h∑j=1NY˜j−Yj2,where h is the space step size. The convergence rate of the numerical results can also be calculated using the following formula:(42)Convergence Rate=C.R,=log2L∞h,2ΔtL∞h,Δt.


### 5.1. Example 1

We consider single-term time-fractional differential equations [24].(43)∑i=12ai∂αiYw,t∂tαi−∂2Yw,t∂w2,=2t−t2+2t2−α2Γ3−α2ew,with initial and boundary conditions(44)Yw,0=0,w∈0,1,Y0,t=t2,Y1,t=et2,t>0.The exact solution(45)Yw,t=t2ew.The method is discussed in Section3 for multiterm time-fractional differential equations that have been used to find a solution of the problem. The values of the parameters are α1=1,α2=0.85, θ=0.5, and ai=1,i=1,2. 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 Δt=h=0.1 and t=0.1 are shown in Figure 1. The graphs clearly show that as the values of w and t 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.Table 1
Error ofY at different step sizes.
hΔtError[24]C.R1/41/165.8667E −032.0112E −02—1/81/641.7728E −034.7622E −031.72651941/161/2565.0769E −041.1347E −031.80401001/321/10241.4073E −042.7067E −041.8510179Figure 1
Solution profile of approximate, exact, and error ofY at h=Δt=0.1,α1=1α2=0.85 and plot at t=1. (a) Approximate solution. (b) Exact solution. (c) Approximate vs Exact. (d) Error.
(a)(b)(c)(d)


### 5.2. Example 2

We consider two-term time-fractional differential equation [14].(46)∑i=12ai∂αiYw,t∂tαi,=∂2Yw,t∂w2+2t2−α1Γ3−α1+2t2−α2Γ3−α2w−w2+21+t2.The initial and boundary conditions are(47)Yw,0=0,w∈0,1,Y0,t=w−w2,Y1,t=0,t>0.The exact solution is(48)Yw,t=1+t2w−w2.The proposed method is applied to find the results of the problem as discussed in Section3 with domain [0, 1]. Taking ai=1,i=1,2,θ=0.5 and for different values of αi are given in Table 2. L2 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 different values of αi. 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 h=Δt=0.1 and t=1. The graph of exact and approximate solutions matches each other, showing good agreement between approximate and exact solutions.Table 2
Errors norm ofY at different step sizes.
α1=0.95,α2=0.2,t=1α1=0.5,α2=0.2,Δt=1/5×104h=2−10h=2−3h=2−4ΔtL2[14]tL2[14]L2[14]1/101.9115E −037.9200E −030.0013.4744E −041.6500E −024.1255E −054.1400E −031/209.8958E −043.7900E −030.011.5960E −048.0400E −031.6479E −052.0000E −031/405.0937E −041.8200E −031.004.6638E −051.8600E −031.6970E −054.6400E −04Figure 2
Graphical solution of approximate, exact, and error ofY at h=Δt=0.1 and plot at t=1. (a) Approximate solution. (b) Exact solution. (c) Approximate vs Exact. (d) Error.
(a)(b)(c)(d)


### 5.3. Example 3

We consider two-term time-fractional differential equation as follows:(49)∑i=12ai∂αiYw,t∂tαi,=∂2Yw,t∂w2+2t2−α1Γ3−α1+2t2−α2Γ3−α2sin2πw+4π2t2sin2πw.The initial and boundary conditions are(50)Yw,0=0,w∈0,1,Y0,t=0,Y1,t=0,t>0.The exact solution is given by(51)Yw,t=t2sin2πw.The proposed method is implemented with domain0,1 to two-term time-fractional DEs. The values of ai=1,i=1,2,θ=0.5 and the results are obtained for different values of αi which are given in Table 3. L2 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 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 Δt=h=0.05 and t=1are displayed in Figure 3. The exact and approximate graphs match each other, showing good agreement between approximate and exact solutions.Table 3
Errors norm ofY at different step sizes.
α1=0.9,α2=0.5,h=0.1,α1=0.5,α2=0.25,h=0.1ΔtL2L∞C.RρMΔtL2L∞C.RρM1/103.444E −034.705E −03—4.37E −011/104.118E −035.615E −03—6.81E −021/201.955E −032.702E −030.80029.13E −011/202.099E −032.899E −030.95372.13E −011/301.434E −032.001E −030.43338.82E −011/301.409E −031.971E −030.55662.96E −011/401.167E −031.641E −030.28619.05E −011/401.060E −031.501E −030.39303.53E −011/501.004E −031.421E −030.2077—1/508.498E −041.217E −030.30263.97E −01Figure 3
Graphical solution of approximate, exact, and error ofY at h=Δt=0.05 and plot at t=1. (a) Approximate solution. (b) Exact solution. (c) Approximate vs Exact. (d) Error.
(a)(b)(c)(d)Figure 4
Stability plot of example 3 for different values ofΔt.


### 5.4. Example 4

We consider three-term time-fractional differential equation [50].(52)∑i=13ai∂αiYw,t∂tαi−∂2Yw,t∂w2,=2w21−w2a1t2−α1Γ3−α1+a2t2−α2Γ3−α2+a3t2−α3Γ3−α3−26w2−6w+11+t2.The initial and boundary conditions are(53)Yw,0=w21−w2,w∈0,1,Y0,t=0,Y1,t=0,t>0.The exact solution is(54)Yw,t=w21−w21+t2.In example 4, the same methodology has been used with domain [0, 1]. Here the values ofa1=a2=1,a3=0.5,α1=0.7,α2=0.6, and α3=0.4. L2 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 Δt=h=0.1,t=1 are displayed in Figure 5. The exact and approximate plots coincide, indicating they are in good agreement.Table 4
Error norm ofY at different step sizes.
hΔtError[50]C.R1/41/167.3088E −032.1871E −011/101/1006.4504E −043.5611E −023.50221/201/4004.2460E −048.9080E −030.60331/401/16001.0858E −042.2258E −031.96731/801/64002.7468E −055.5614E −041.9829Figure 5
Graphical solution of approximate, exact, and error ofY at h=Δt=0.1 and plot at t=1. (a) Approximate solution. (b) Exact solution. (c) Approximate vs Exact. (d) Error.
(a)(b)(c)(d)


### 5.5. Example 5

We consider three-term time-fractional differential equation [50].(55)∑i=13ai∂αiYw,t∂tαi−∂2Yw,t∂w2,=w1+α11−wΓα2+α3+1a1tα2+α3−α1Γα2+α3−α1+1+a2tα3Γα3+1+a3tα2Γα2+1−1+α1wα1−1α1−2+α1w1+tα2+α3.The initial and boundary conditions are(56)Yw,0=w1+α11−w,w∈0,1,Y0,t=0Y1,t=0,t>0.The exact solution(57)Yw,t=w1+α11−w1+tα2+α3.The solution is computed over the domain [0, 1] and the values of parameters area1=a3=1,a2=0.5, θ=0.5, α1=0.8,α2=0.7, and α3=0.6 . 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 t=1, Δt=0.1, and h=0.1 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.Table 5
Errors norm ofY at different step sizes.
hΔtError[50]C.R1/41/163.1781E −032.4880E −02—1/101/1004.7588E −046.1204E −032.73951/201/4001.3246E −042.0016E −031.84501/401/16003.6478E −056.3605E −041.86051/801/64001.0012E −051.9816E −041.8653Figure 6
Graphical solution of approximate, exact, and error ofY at h=Δt=0.1 and plots at t=1. (a) Approximate solution. (b) Exact solution. (c) Plots for different value of α’s. (d) Error.
(a)(b)(c)(d)


## 5.1. Example 1

We consider single-term time-fractional differential equations [24].(43)∑i=12ai∂αiYw,t∂tαi−∂2Yw,t∂w2,=2t−t2+2t2−α2Γ3−α2ew,with initial and boundary conditions(44)Yw,0=0,w∈0,1,Y0,t=t2,Y1,t=et2,t>0.The exact solution(45)Yw,t=t2ew.The method is discussed in Section3 for multiterm time-fractional differential equations that have been used to find a solution of the problem. The values of the parameters are α1=1,α2=0.85, θ=0.5, and ai=1,i=1,2. 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 Δt=h=0.1 and t=0.1 are shown in Figure 1. The graphs clearly show that as the values of w and t 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.Table 1
Error ofY at different step sizes.
hΔtError[24]C.R1/41/165.8667E −032.0112E −02—1/81/641.7728E −034.7622E −031.72651941/161/2565.0769E −041.1347E −031.80401001/321/10241.4073E −042.7067E −041.8510179Figure 1
Solution profile of approximate, exact, and error ofY at h=Δt=0.1,α1=1α2=0.85 and plot at t=1. (a) Approximate solution. (b) Exact solution. (c) Approximate vs Exact. (d) Error.
(a)(b)(c)(d)


## 5.2. Example 2

We consider two-term time-fractional differential equation [14].(46)∑i=12ai∂αiYw,t∂tαi,=∂2Yw,t∂w2+2t2−α1Γ3−α1+2t2−α2Γ3−α2w−w2+21+t2.The initial and boundary conditions are(47)Yw,0=0,w∈0,1,Y0,t=w−w2,Y1,t=0,t>0.The exact solution is(48)Yw,t=1+t2w−w2.The proposed method is applied to find the results of the problem as discussed in Section3 with domain [0, 1]. Taking ai=1,i=1,2,θ=0.5 and for different values of αi are given in Table 2. L2 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 different values of αi. 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 h=Δt=0.1 and t=1. The graph of exact and approximate solutions matches each other, showing good agreement between approximate and exact solutions.Table 2
Errors norm ofY at different step sizes.
α1=0.95,α2=0.2,t=1α1=0.5,α2=0.2,Δt=1/5×104h=2−10h=2−3h=2−4ΔtL2[14]tL2[14]L2[14]1/101.9115E −037.9200E −030.0013.4744E −041.6500E −024.1255E −054.1400E −031/209.8958E −043.7900E −030.011.5960E −048.0400E −031.6479E −052.0000E −031/405.0937E −041.8200E −031.004.6638E −051.8600E −031.6970E −054.6400E −04Figure 2
Graphical solution of approximate, exact, and error ofY at h=Δt=0.1 and plot at t=1. (a) Approximate solution. (b) Exact solution. (c) Approximate vs Exact. (d) Error.
(a)(b)(c)(d)


## 5.3. Example 3

We consider two-term time-fractional differential equation as follows:(49)∑i=12ai∂αiYw,t∂tαi,=∂2Yw,t∂w2+2t2−α1Γ3−α1+2t2−α2Γ3−α2sin2πw+4π2t2sin2πw.The initial and boundary conditions are(50)Yw,0=0,w∈0,1,Y0,t=0,Y1,t=0,t>0.The exact solution is given by(51)Yw,t=t2sin2πw.The proposed method is implemented with domain0,1 to two-term time-fractional DEs. The values of ai=1,i=1,2,θ=0.5 and the results are obtained for different values of αi which are given in Table 3. L2 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 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 Δt=h=0.05 and t=1are displayed in Figure 3. The exact and approximate graphs match each other, showing good agreement between approximate and exact solutions.Table 3
Errors norm ofY at different step sizes.
α1=0.9,α2=0.5,h=0.1,α1=0.5,α2=0.25,h=0.1ΔtL2L∞C.RρMΔtL2L∞C.RρM1/103.444E −034.705E −03—4.37E −011/104.118E −035.615E −03—6.81E −021/201.955E −032.702E −030.80029.13E −011/202.099E −032.899E −030.95372.13E −011/301.434E −032.001E −030.43338.82E −011/301.409E −031.971E −030.55662.96E −011/401.167E −031.641E −030.28619.05E −011/401.060E −031.501E −030.39303.53E −011/501.004E −031.421E −030.2077—1/508.498E −041.217E −030.30263.97E −01Figure 3
Graphical solution of approximate, exact, and error ofY at h=Δt=0.05 and plot at t=1. (a) Approximate solution. (b) Exact solution. (c) Approximate vs Exact. (d) Error.
(a)(b)(c)(d)Figure 4
Stability plot of example 3 for different values ofΔt.


## 5.4. Example 4

We consider three-term time-fractional differential equation [50].(52)∑i=13ai∂αiYw,t∂tαi−∂2Yw,t∂w2,=2w21−w2a1t2−α1Γ3−α1+a2t2−α2Γ3−α2+a3t2−α3Γ3−α3−26w2−6w+11+t2.The initial and boundary conditions are(53)Yw,0=w21−w2,w∈0,1,Y0,t=0,Y1,t=0,t>0.The exact solution is(54)Yw,t=w21−w21+t2.In example 4, the same methodology has been used with domain [0, 1]. Here the values ofa1=a2=1,a3=0.5,α1=0.7,α2=0.6, and α3=0.4. L2 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 Δt=h=0.1,t=1 are displayed in Figure 5. The exact and approximate plots coincide, indicating they are in good agreement.Table 4
Error norm ofY at different step sizes.
hΔtError[50]C.R1/41/167.3088E −032.1871E −011/101/1006.4504E −043.5611E −023.50221/201/4004.2460E −048.9080E −030.60331/401/16001.0858E −042.2258E −031.96731/801/64002.7468E −055.5614E −041.9829Figure 5
Graphical solution of approximate, exact, and error ofY at h=Δt=0.1 and plot at t=1. (a) Approximate solution. (b) Exact solution. (c) Approximate vs Exact. (d) Error.
(a)(b)(c)(d)


## 5.5. Example 5

We consider three-term time-fractional differential equation [50].(55)∑i=13ai∂αiYw,t∂tαi−∂2Yw,t∂w2,=w1+α11−wΓα2+α3+1a1tα2+α3−α1Γα2+α3−α1+1+a2tα3Γα3+1+a3tα2Γα2+1−1+α1wα1−1α1−2+α1w1+tα2+α3.The initial and boundary conditions are(56)Yw,0=w1+α11−w,w∈0,1,Y0,t=0Y1,t=0,t>0.The exact solution(57)Yw,t=w1+α11−w1+tα2+α3.The solution is computed over the domain [0, 1] and the values of parameters area1=a3=1,a2=0.5, θ=0.5, α1=0.8,α2=0.7, and α3=0.6 . 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 t=1, Δt=0.1, and h=0.1 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.Table 5
Errors norm ofY at different step sizes.
hΔtError[50]C.R1/41/163.1781E −032.4880E −02—1/101/1004.7588E −046.1204E −032.73951/201/4001.3246E −042.0016E −031.84501/401/16003.6478E −056.3605E −041.86051/801/64001.0012E −051.9816E −041.8653Figure 6
Graphical solution of approximate, exact, and error ofY at h=Δt=0.1 and plots at t=1. (a) Approximate solution. (b) Exact solution. (c) Plots for different value of α’s. (d) Error.
(a)(b)(c)(d)


## 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 usingL2 and L∞ 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.



---

*Source: 3160725-2022-11-23.xml*