[picture] | Computer Modeling in Engineering & Sciences | [picture]

DOI: 10.32604/cmes.2021.012720

ARTICLE

Redefined Extended Cubic B-Spline Functions for Numerical Solution of Time-Fractional Telegraph Equation

Muhammad Amin¹, Muhammad Abbas²,*, Dumitru Baleanu³,⁴,⁵, Muhammad Kashif Iqbal⁶ and Muhammad Bilal Riaz⁷

1 Department of Mathematics, National College of Business Administration & Economics, Lahore, 54660, Pakistan
2 Department of Mathematics, University of Sargodha, Sargodha, 40100, Pakistan
3 Department of Mathematics, Faculty of Arts and Sciences, Cankaya University, Ankara, 06530, Turkey
4 Department of Medical Research, China Medical University Hospital, China Medical University, Taichung, 40402, Taiwan
5 Institute of Space-Sciences, Bucharest, 077125, Romania
6 Department of Mathematics, Government College University, Faisalabad, 38000, Pakistan
7 Department of Mathematics, University of Management and Technology, Lahore, 54700, Pakistan
Corresponding Author: Muhammad Abbas. Email: muhammad.abbas@uos.edu.pk
Received: 10 July 2020; Accepted: 21 December 2020

Abstract: This work is concerned with the application of a redefined set of extended uniform cubic B-spline (RECBS) functions for the numerical treatment of time-fractional Telegraph equation. The presented technique engages finite difference formulation for discretizing the Caputo time-fractional derivatives and RECBS functions to interpolate the solution curve along the spatial grid. Stability analysis of the scheme is provided to ensure that the errors do not amplify during the execution of the numerical procedure. The derivation of uniform convergence has also been presented. Some computational experiments are executed to verify the theoretical considerations. Numerical results are compared with the existing schemes and it is concluded that the present scheme returns superior outcomes on the topic.

Keywords: Extended cubic B-spline; redefined extended cubic B-spline; time fractional telegraph equation; caputo fractional derivative; finite difference method; convergence

1  Introduction

In recent years, fractional calculus has gained a remarkable importance. Fractional derivatives and integrals have manifold applications in science and engineering such as fluid mechanics, chemical physics, electricity, control theory, biomedical, epidemic diseases, hydrology, electro-chemistry, probability theory, signal processing, heat conduction and diffusion problems [1–7]. Many researchers developed fractional-order models to describe real-world problems and studied their analytical and numerical solutions [8–11]. These models involve different types of fractional derivative operators [12–15]. The fractional telegraph equation is one of the fundamental mathematical models arising in the study of electrical signals in transmission line and wave phenomena [16–18]. Basically, it belongs to the family of hyperbolic partial differential equations. Several numerical and analytical techniques have been proposed for solving these type of equations. In [19], the authors employed Adomian decomposition method for solving time and space fractional telegraph equations. Dehghan et al. [20] proposed variational method to explore the series solution to multi space telegraph equation. The authors in [21], employed Homotopy analysis method to explore the analytical solution of telegraph equation involving fractional time derivative. Later on, Hayat et al. [22] used Homotopy perturbation technique to study time fractional telegraph equation. They handled TFTE for both brownian and standard motion. In [23], Wei et al. applied fully discrete local discontinuous Galerkin finite element method to solve fractional telegraph equation. Hosseini et al. [24] studied the numerical solution of fractional telegraph equation by means of radial basis functions. Srivastava et al. [25] employed reduced differential transformation method for second order hyperbolic time fractional telegraph equation in one dimensional space. Wang et al. [26] analyzed an H1-Galerkin mixed finite element method for the numerical solution of time fractional telegraph equation. Modanli et al. [27] solved fractional order telegraph equation by means of Theta method. Xu et al. [28] applied Legendre wavelets direct method for solving fractional order telegraph equation. In [29], Wang et al. utilized spectral Galerkin approximation to study the approximate solution of TFTE. Kamran et al. [30] studied the numerical solution of TFTE by means of a Localized kernal-based approach. Here, in this work, we consider the following fractional order telegraph equation.

∂α∂tαv(s,t)+∂β∂tβv(s,t)+γ1v(s,t)-γ2∂2∂s2v(s,t)=f(s,t),(s,t)∈(0,L)×(0,T),(1)

v(s,0)=ψ1(s),vt(s,0)=ψ2(s),s∈[0,L],(2)

v(0,t)=ϕ1(t),v(L,t)=ϕ2(t),t∈[0,T],(3)

where ψj(s) and ϕj(t) (j=1,2) are given and ∂α∂tαv(s,t), ∂β∂tβv(s,t) represent the Caputo fractional derivatives of order α and β, respectively. It is worth mentioning that in (1), α∈(1,2] and β∈(0,1]. However, this work is restricted to the class of problems involving α=β+1 and α=2β.

In this paper, we have studied the application of a redefined form of extended cubic B-spline (ECBS) functions for the numerical treatment of time-fractional Telegraph equation (TFTE). These functions are generalized forms of cubic B-spline functions involving one free shape parameter which provides the flexibility to modify the shape of the solution curve [31]. Although, the degree of the piecewise polynomials is enhanced by one and the continuity of RECBS remains of order three. A finite-difference formula is used for the discretization of the Caputo time-fractional derivative. Usually, in collocation techniques, the Dirichlet’s type end conditions are imposed where the basis of spline functions vanish, but the typical ECBS functions do not vanish at boundaries. We have employed RECBS functions for spatial discretization, as these basis functions die out on the boundaries where the Dirichlet’s types of conditions are specified. The present approach is novel for the approximate solution of fractional PDEs and as far as we are aware, it has never been employed for this purpose before.

The manuscript is composed as: Section 2 describes the redefined extended cubic B-spline functions. In Section 3, the numerical method has been explained. In Section 4, the stability analysis of proposed method is presented. In Section 5, we have derived the results for theoretical convergence. The approximate results and discussion are reported in Section 6. Finally, the concluding remarks have been given in Section 7.

2  Redefined Extended Cubic B-Spline Functions

Suppose the spatial domain [a,b] be portioned into M parts of equal length h=b-aM such that a=s0<s1<⋯<sM=b, where sm = s0 +mh, m = 0:1:M. We assume the ECBS approximation V*(s,t) for a sufficiently smooth function v(s,t) as

V*(s,t)=∑m=-1M+1ξm(t)λm(s,κ),(4)

where ξm(t) are real constants and λm(s,κ) are ECBS functions [32]:

λm(s,κ)=124h4{4(1-κ)h(s-sm-2)3+3κ(s-sm-2)4,ifs∈[sm-2,sm-1)(4-κ)h4+12h3(s-sm-1)+6h2(2+κ)(s-sm-1)2-12h(s-sm-1)3-3κ(s-sm-1)4,ifs∈[sm-1,sm)(4-κ)h4-12h3(s-sm+1)-6h2(2+κ)(s-sm+1)2+12h(s-sm+1)3+3κ(s-sm-1)4,ifs∈[sm,sm+1)-4h(1-κ)(s-sm+2)3-3κ(s-sm+2)4,ifs∈[sm+1,sm+2)0,otherwise(5)

where -8≤κ≤1 is responsible for fine tuning the shape of the curve. The approximate solution (V*)mr=V*(sm,tr) and its first two derivatives with respect to space variable s, at mth knot and rth time step, in terms of ξm can be expressed as

{(V*)mr=b1ξm-1r+b2ξmr+b1ξm+1r,(Vs*)mr=-b3ξm-1r+b3ξm+1r,(Vss*)mr=b4ξm-1r+b5ξmr+b4ξm+1r,(6)

where b1=4-κ24, b2=16+2κ24, b3=12h, b4=2+κ2h2, b5=-4-2κ2h2. The ECBS functions λ-1,λ0,…,λM+1 do not vanish at the boundaries when Dirichlet type end conditions are imposed. Therefore, we redefine these functions in such a manner that the resulting basis vanish at the boundaries [33]. We eliminate ξ-1r and ξM+1r from Eq. (4) as

V(s,t)=Φ(s,t)+∑m=0Mξmr(t)λ̃m(s,κ),(7)

where the weight function Φ(s,t) and redefined ECBS (RECBS) functions are given by

Φ(s,t)=λ-1(s,κ)λ-1(s0,κ)ϕ1(t)+λM+1(s,κ)λM+1(sM,κ)ϕ2(t),(8)

{λ̃m(s,κ)=λm(s,κ)-λm(s0,κ)λ-1(s0,κ)λ-1(s,κ),m=0,1,λ̃m(s,κ)=λm(s,κ),m=2:1:M-2,λ̃m(s,κ)=λm(s,κ)-λm(sM,κ)λM+1(sM,κ)λM+1(s,κ),m=M-1,M.(9)

3  Numerical Technique

We divide the time domain [0,T] into R subintervals [tr,tr+1] s.t. tr=rΔt, r=0,1,2,…,R and Δt=TR. The Caputo’s time fractional derivative at t = tr+1, for α∈(1,2], can be discretized as

∂α∂tαv(s,tr+1)=∫t0tr+1∂2v(s,w)∂w2(tr+1-w)-α+1Γ(2-α)dw.=1Γ(2-α)∑j=0r∫tjtj+1∂2v(s,w)∂w2(tr+1-w)-α+1dw=1Γ(2-α)∑j=0rv(s,tj+1)-2v(s,tj)+v(s,tj-1)Δt2∫tjtj+1(tr+1-w)-α+1dw+(Eα)Δtr+1=1Γ(2-α)∑j=0rv(s,tj+1)-2v(s,tj)+v(s,tj-1)Δt2∫tr-jtr-j+1(υ)-α+1dυ+(Eα)Δtr+1=1Γ(2-α)∑j=0rv(s,tr-j+1)-2v(s,tr-j)+v(s,tr-j-1)Δt2∫tjtj+1(υ)-α+1dυ+(Eα)Δtr+1=1Γ(3-α)∑j=0rv(s,tr-j+1)-2v(s,tr-j)+v(s,tr-j-1)Δtα[(j+1)2-α-j2-α]+(Eα)Δtr+1=1Γ(3-α)∑j=0rpjv(s,tr-j+1)-2v(s,tr-j)+v(s,tr-j-1)Δtα+(Eα)Δtr+1,(10)

where pj=(j+1)2-α-(j)2-α, υ=(tr+1-w) and (Eα)Δtr+1 is the truncation error.

Also

|(Eα)Δtr+1|≤ρ1(Δt)2-α,(11)

where ρ1 is constant and

•   pj∈ℤ+, ∀j

•   1=p0>p1>p2>p3>⋯>pr, pr→0 as r→∞

•   (2p0-p1)+∑j=1r-1(-pj+1+2pj-pj-1)+(2pr-pr-1)-pr=1

Similarly,

∂β∂tβv(s,tr+1)=∫t0tr+1∂v(s,w)∂w(tr+1-w)-βΓ(1-β)dw=1Γ(1-β)∑j=0r∫tjtj+1∂v(s,w)∂w(tr+1-w)-βdw=1Γ(1-β)∑j=0rv(s,tj+1)-v(s,tj)Δt∫tjtj+1(tr+1-w)-βdw+(Eβ)Δtr+1=1Γ(1-β)∑j=0rv(s,tj+1)-v(s,tj)Δt∫tr-jtr-j+1(υ)-βdυ+(Eβ)Δtr+1=1Γ(1-β)∑j=0rv(s,tr-j+1)-v(s,tr-j)Δt∫tjtj+1(υ)-βdυ+(Eβ)Δtr+1=1Γ(2-β)∑j=0rv(s,tr-j+1)-v(s,tr-j)Δtβ[(j+1)1-β-j1-β]+(Eβ)Δtr+1=1Γ(2-β)∑j=0rqjv(s,tr-j+1)-v(s,tr-j)Δtβ+(Eβ)Δtr+1.(12)

where qj=(j+1)1-β-(j)1-β, υ=(tr+1-w) and (Eβ)Δtr+1 is the truncation error.

Also

|(Eβ)Δtr+1|≤ρ2(Δt)1-β,(13)

where ρ2 is constant and

•   qj∈ℤ+, ∀j

•   1=q0>q1>q2>q3>⋯>qr,qr→0 as r→∞

•   ∑j=0r(qj-qj+1)+qr+1=(q0-q1)+∑j=1r-1(qj-qj+1)+qr=1

Substituting (10) and (12) in (1) at t = tr+1, we get

∑j=0rpjv(s,tr-j+1)-2v(s,tr-j)+v(s,tr-j-1)ΔtαΓ(3-α)+∑j=0rqjv(s,tr-j+1)-v(s,tr-j)ΔtβΓ(2-β)+γ1v(s,tr+1)-γ2∂2∂s2v(s,tr+1)=f(s,tr+1),r=0,1,2,…,R.(14)

Using theta-weighted scheme for θ=1, Eq. (14) takes the following form

α1∑j=0rpj(vr-j+1-2vr-j+vr-j-1)+β1∑j=0rqj(vr-j+1-vr-j)+γ1vr+1-γ2(vss)r+1=fr+1,r=0,1,2,…,R,(15)

where α1=1ΔtαΓ(3-α), β1=1ΔtβΓ(2-β), v(s, tr+1) = vr+1.

For r = 0, v−1 appears in Eq. (15). We use the initial conditions and substitute v-1=v0-Δtψ2(s) to get the following equation

(α1+β1+γ1)v1-γ2(vss)1=(α1+β1)v0+α1Δtψ2(s)+f1.(16)

For r=1,2,…,R, Eq. (15) is reshaped as

(α1+β1+γ1)vr+1-γ2(vss)r+1=(2α1+β1)vr-α1∑j=1rpj(vr-j+1-2vr-j+vr-j-1)-β1∑j=1rqj(vr-j+1-vr-j)-α1vr-1+fr+1.(17)

Now, we discretize the spatial domain [a, b] by M + 1 equally spaced knots a=s0,s1,s2,…,sM=b such that sm = s0 +mh, m=0,1,…,M and assume that the RECBS approximation V(s, t) for the exact solution v(s, t) is given by

V(s,t)=Φ(s,t)+∑m=0Mξmr(t)λ̃m(s,κ),(18)

where Φ(s,t) and λ̃m(s,κ) are defined in (8) and (9), respectively.

Solution at t = t1

The initial solution is given in (2). However, the control points ξi at t = t1 are required to start the main scheme (17). For this purpose, (18) is substituted into (16) to get the following system of equations

(α1+β1+γ1)[Φi1+∑m=i-1i+1ξm1λ̃m(si,κ)]-γ2[(Φss)i1+∑m=i-1i+1ξm1(λ̃m)ss(si,κ)]=(α1+β1)vi0+α1Δtψ2(si)+fi1,i=0,1,…,M.(19)

Solving (19), we get [ξ01,ξ11,…,ξM1]T and substitute these control points into (18) to obtain the approximate solution at t = t1

Solution at t=tr+1,r=1,2,…,R

Using (18) in Eq. (17), we obtain

(α1+β1+γ1)[Φir+1+∑m=i-1i+1ξmr+1λ̃m(si,κ)]-γ2[(Φss)ir+1+∑m=i-1i+1ξmr+1(λ̃m)ss(si,κ)]=(2α1+β1)[Φir+∑m=i-1i+1ξmrλ̃m(si,κ)]-α1∑j=1rpj[Φir-j+1-2Φir-j+Φir-j-1+∑m=i-1i+1(ξmr-j+1-2ξmr-j+ξmr-j-1)λ̃m(si,κ)]-β1∑j=1rqj[Φir-j+1-Φir-j+∑m=i-1i+1(ξmr-j+1-ξmr-j)λ̃m(si,κ)]-α1[Φir-1+∑m=i-1i+1ξmr-1λ̃m(si,κ)]+fir+1,i=0,1,2,…,M.(20)

Eq. (20) represents a set of (M + 1) equations involving (M + 1) unknowns. This system of equations is solved to for ξir+1 and their values are plugged into (18) to get the required solution at (r + 1)th time level.

4  Stability

We apply Fourier method to study the stability of our numerical method. Let εmr and ε̃mr denote the Fourier growth factor and its approximate value. We introduce the error term ϱmr as

ϱmr=εmr-ε̃mr,m=1:1:M-1,r=0:1:R,(21)

where ϱr=[ϱ1r,ϱ2r,…,ϱM-1r]T. Using (21) in (20), the error equation at (r + 1)st time level is given by

(α1+β1+γ1)[b1ϱm-1r+1+b2ϱmr+1+b1ϱm+1r+1]-γ2[b4ϱm-1r+1+b5ϱmr+1+b4ϱm+1r+1]=(2α1+β1)[b1ϱm-1r+b2ϱmr+b1ϱm+1r]-α1∑j=1rpj[b1(ϱm-1r-j+1-2ϱm-1r-j+ϱm-1r-j-1)+b2(ϱmr-j+1-2ϱmr-j+ϱmr-j-1)+b1(ϱm+1r-j+1-2ϱm+1r-j+ϱm+1r-j-1)]-β1∑j=1rqj[b1(ϱm-1r-j+1-ϱm-1r-j)+b2(ϱmr-j+1-ϱmr-j)+b1(ϱm+1r-j+1-ϱm+1r-j)]-α1[b1ϱm-1r-1+b2ϱmr-1+b1ϱm+1r-1],m=1,2,…,M-1.(22)

If ϱmr=εreινmh, where ι=-1 and ν=2πmb-a, then (22) is reshaped as

[(α1+β1+γ1)(2b1 cosνh+b2)-γ2(2b4 cosνh+b5)]εr+1=γ4[2b1 cosνh+b2]εr-α1(2b1 cosνh+b2)∑j=1rpj[εr-j+1-2εr-j+εr-j-1]-β1(2b1 cosνh+b2)∑j=1rqj[εr-j+1-εr-j]-α1[2b1 cosνh+b2]εr-1.(23)

After simplifying (23), we get the following result

εr+1=1η[(1+η1)εr-η1∑j=1rpj[εr-j+1-2εr-j+εr-j-1]-η2∑j=1rqj[εr-j+1-εr-j]-η1εr-1],(24)

where η=1+η3+12η4(2+κ)sin2(νh/2)h2[6+(4-κ)sin2(νh/2)]≥1, η1=α1α1+β1, η2=β1α1+β1, η3=γ1α1+β1 and η4=γ2α1+β1.

For r = 0, the expression (23) takes the following form

|ε1|=1η|(1+η1)ε0|≤(1+η1)|ε0|,∵η≥1.

Now, assuming |εr|≤(1+η1)|ε0| for r > 1, we use (24) to proceed as

|εr+1|=1η[(1+η1)|εr|-η1∑j=1rpj[|εr-j+1|-2|εr-j|+εr-j-1]-η2∑j=1rqj[|εr-j+1|-|εr-j|]-η1|εr-1|]≤(1+η1)2|ε0|-η1(1+η1)∑j=1rpj[|ε0|-2|ε0|+ε0]-η2(1+η1)∑j=1rqj[|ε0|-|ε0|]-η1(1+η1)|ε0|=(1+η1)2|ε0|-η1(1+η1)|ε0|=(1+η1)[1+η1-η1]|ε0|⇒|εr+1|≤(1+η1)|ε0|,∀r.

Consequently, following [34], we have |ϱr|=(1+η1)|ϱ0|,r=0,1,…,R.

Hence, the scheme stable.

5  Convergence

Let Ṽ(s,t)=∑m=0Mdm(t)λ̃m(s) be the computed ECBS for the numerical solution V(s, t) and the analytical solution v(s, t) subject to the interpolating conditions LṼ(sm,t)=f̃(sm,t),m=0,1,…,M. Now, the problem (1) in terms of difference equation L(Ṽ(sm,t)-V(sm,t)), at t = tr, is given by

(α1b1+β1b1+γ1b1+γ2b4)ζm-1r+1+(α1b2+β1b2+γ1b2+γ2b5)ζmr+1+(α1b1+β1b1+γ1b1+γ2b4)ζm+1r+1=(2α1+β1)(b1ζm-1r+b2ζmr+b1ζm+1r)-α1∑k=1rpk[b1(ζm-1r-k+1-2ζm-1r-k+ζm-1r-k-1)+b2(ζmr-k+1-2ζmr-k+ζmr-k-1)+b1(ζm+1r-k+1-2ζm+1r-k+ζm+1r-k-1)]-β1∑k=1rqk[b1(ζm-1r-k+1-ζm-1r-k)+b2(ζmr-k+1-ζmr-k)+b1(ζm+1r-k+1-ζm+1r-k)]-α1(b1ζm-1r-1+b2ζmr-1+b1ζm+1r-1)+fmr+1,m=0,1,…,M,(25)

where ζmr=ξmr-dmr and lmr=h2[fmr-f̃mr].

The boundary conditions can be rewritten as

b1ζm-1r+1+b2ζmr+1+b1ζm+1r+1=0,m=0,M.

Moreover, following [34], we have

∥Dj(v(s,t)-Ṽ(s,t))∥∞≤Ϝjh4-j,j=0,1,2,(26)

Therefore, |lmr|=h2|fmr-f̃mr|≤Ϝa4, where Ϝ does not depend on mesh spacing.

Now, we introduce lr=maxm=0M{|lmr|}, ẽmr=|ζmr| and ẽr=m= 0maxM{|emr|}. For r = 0, Eq. (25) transforms into following relation

(α1b2+β1b2+γ1b2+γ2b5)ζm1=(α1b1+γ2b4)(ζm+11-ζm-11)+(β1b1+γ1b1)(ζm+11-ζm-11)+1h2lm1.

Involving the absolute values of lmr and ζmr, we obtain

ẽm1≤6Ϝh42α1h2(2+κ)+12(2+κ)γ1+6γ2h

Hence, employing the end constraints, we get ẽ1≤Ϝ1h2, where Ϝ1 is independent of spatial grid spacing.

Now, assuming that ẽmk≤Ϝrh2 for r > 1, we set Ϝ=maxj=0r{Ϝj} and plug in absolute values of lmr and ζmr in Eq. (25)

ẽmr+1≤6Ϝh2(α1+β1+γ1)h2(2+κ)-12(2+κ)γ2l[(2α1+β1)(b1ζm-1r+b2ζmr+b1ζm+1r)-α1(b1ζm-1r-1+b2ζmr-1+b1ζm+1r-1)-(α1∑k=0r-1(pj-1-2pj+pj+1)+β1∑j=0r-1(pj-1-pj))Ϝh2+Ϝh2r].

Utilizing the boundary conditions, we obtain ẽmr+1≤Ϝh2.

Hence, the last result is true for all r. Using the result ∑m=0M|λm(s,κ)|≤1.75 [34], we get

∥Ṽ(s,t)-V(s,t)∥∞=∥∑m=0M(dm(t)-ξm(t))λm(s,κ)∥∞≤1.75Ϝh2,(27)

Consequently, using (26) and (27), we get

∥v(s,t)-V(s,t)∥∞≤∥v(s,t)-Ṽ(s,t)∥∞+∥Ṽ(s,t)-V(s,t)∥∞≤Ϝ0h4+1.75Ϝh2.

Hence, in the light of above discussion together with (11) and (13), we conclude that the scheme is O(h2) accurate in spatial direction. However, (11) and (12) imply that the truncation error in temporal direction is O(Δt2-α+Δt1-β). This work is restricted to the class of problems involving α=β+1 and α=2β. Therefore, theoretically the scheme is O(Δt2-α) when α=1+β and O(Δt1-α/2) when α=2β.

6  Numerical Results

To investigate the accuracy of presented technique, some numerical experiments are presented. For this purpose, following error norms have been used

L∞=m= 0maxM|Vm-vm|,L2=h∑m=0M|Vm-vm|2,

Also, the experimental order of convergence (EOC) is computed by following important formula [35]:

EOC=1log2log[L∞(2m)L∞(m)]

Example 6.1. As the first experiment, we take the following multi term TFTE [29]

∂α∂tαv(s,t)+∂β∂tβv(s,t)+v(s,t)-∂2∂s2v(s,t)=f(s,t),(s,t)∈[-1,1]×[0,T],v(s,0)=0,vt(s,0)=0,v(-1,t)=v(1,t)=0, where α=β+1 andf(s,t)=2[t1-βΓ(2-β)+t2-βΓ(3-β)+t22+π2t22]sin(πs).

The exact solution of the problem is v(s,t)=t2 sin(πs).

The absolute error and temporal order of convergence for Example 6.1 along temporal direction using M = 24 and different values of β are reported in Tab. 1. It can easily be seen that our results are more accurate than the scheme based on generalized finite difference method (GFDM) [29]. In Tab. 2, we have computed the absolute errors by setting M = 24, 28 and Δt=0.1 corresponding to different grid points in spatial direction. Tab. 3 gives spatial order of convergence (EOC) subject to β=0.6 and Δt=0.1. The experimental rate of convergence of the current method is found to be in line with the theoretical appraisal. Fig. 1 shows the physical behaviour of approximate solutions at different time levels when β=0.1, M = 24 and Δt=0.1. The 3D visuals of exact and numerical solutions with β=0.1, M = 24 and Δt=0.1 are shown in Fig. 2, whereas, Fig. 3 depicts the absolute error between the exact and approximate solutions using β=0.1, M = 36 and Δt=0.1.

Table 1: Experimental order of convergence (EOC) for Example 6.1 when M = 24 using different values of β

Table 2: Absolute errors for Example 6.1 when Δt=0.1 using different values of β

Table 3: Experimental order of convergence (EOC) for Example 6.1, when β=0.6 and Δt=0.1

Figure 1: Exact and numerical solution for Example 6.1 at different time levels when Δt=0.1, β=0.1 and M = 24

Figure 2: Exact and approximate solution for Example 6.1 with M = 24, Δt=0.1 and β=0.1. (a) 3D plot for exact solution. (b) 3D plot for approximate solution

Figure 3: Absolute error for Example 6.1 when M = 36, β=0.1 and Δt=0.1

The piecewise defined approximate solution for Example 6.1 using proposed algorithm, when β=0.50, -1≤s≤1, M = 20, Δt=0.01, is given by

V(s)={1.48811+s(10.1756+s(12.1734+(2.91397-0.571918s)s)),if s∈[-1.00,-0.90]0.736724+s(6.87798+s(6.74703+(-1.05394-1.65977s)s)),if s∈[-0.90,-0.80]0.307397+s(4.79421+s(2.95783+(-4.11339-2.58515s)s)),if s∈[-0.80,-0.70]0.0996004+s(3.67304+s(0.697273+(-6.1311-3.25748s)s)),if s∈[-0.70,-0.60]⋮⋮-0.0000571467+s(3.13881+s(-0.0497796+(-5.58514-1.65977s)s)),if s∈[-0.20,-0.10]-2.33841×10-15+s(3.14161+s(2.84217×10-14+(-5.20164-0.571918s)s)),if s∈[-0.10,0.00]⋮⋮-0.0996004+s(3.67304+s(-0.697273+s(-6.1311+3.25748s))),if s∈[0.60,0.70]-0.307397+s(4.79421+s(-2.95783+s(-4.11339+2.58515s))),if s∈[0.70,0.80]-0.736724+s(6.87798+s(-6.74703+s(-1.05394+1.65977s))),if s∈[0.80,0.90]-1.48811+s(10.1756+s(-12.1734+(2.91397+0.571918s)s)),if s∈[0.90,1.00]

Example 6.2. Consider the TFTE [27]

∂α∂tαv(s,t)+∂β∂tβv(s,t)-∂2∂s2v(s,t)+v(s,t)=f(s,t),(s,t)∈[0,π]×[0,T],v(s,0)= sin(s),vt(s,0)=0,v(0,t)=v(π,t)=0, where α=2, 0<β≤1 andf(s,t)= sin(s)(6t+6t3-βΓ(4-β)+2(t3+1)).

The analytical solution to this problem is sin(s)(t3+1).

The approximate analytical solution for Example 6.2 using proposed method, when β=0.50, 0≤s≤π, M = 20 and Δt=0.01 is given by

V(s)={-6.76345×10-33+s(1.99275+s(2.22585×10-15+(0.837877-5.58309s)s)),if s∈[0.00,π20]-0.0338266+s(2.68154+s(-4.92919+(14.7627-16.6118s)s)),if s∈[π20,π10]-0.565727+s(8.09009+s(-24.2414+(41.9261-27.2315s)s)),if s∈[π10,3π20]-3.18079+s(25.7792+s(-66.1976+(80.9809-37.1806s)s)),if s∈[3π20,π5]⋮⋮-218.272+s(658.404+s(-739.358+(369.326-69.1931s)s)),if s∈[2π5,9π20]-350.039+s(942.959+s(-947.621+(423.37-70.9399s)s)),if s∈[9π20,π2]⋮⋮-1686.35+s(2603.73+s(-1504.72+(386.245-37.1806s)s)),if s∈[4π5,17π20]-1567.02+s(2280.23+s(-1241.68+(300.275-27.2315s)s)),if s∈[17π20,9π10]-1200.66+s(1651.46+s(-849.506+(193.987-16.6118s)s)),if s∈[9π10,19π20]-511.604+s(665.642+s(-322.721+(69.3213-5.58309s)s)),if s∈[19π20,π]

The absolute numerical errors in RECBS solution for Example 6.2 setting Δt=h at different values of β are listed in Tab. 4. It is clear that our results have better agreement with the exact solution in comparison to the theta-method (TM) [27]. Fig. 4 shows the physical behaviour of approximate solutions at different time levels when β=0.5, M = 40 and Δt=0.025. The 3D visuals of exact and numerical solutions with β=0.5, M = 40 and Δt=0.025 are shown in Fig. 5. Whereas, Fig. 6 depicts the absolute error between the exact and approximate solutions using β=0.75, M = 40 and Δt=0.025.

Table 4: Absolute error norms for Example 6.2 using different values of M and β

Figure 4: Exact and numerical solution for Example 6.2 at different time levels when Δt=0.025, β=0.5 and M = 40

Example 6.3 Consider the multi term TFTE [30]

∂α∂tαv(s,t)+∂β∂tβv(s,t)-∂2∂s2v(s,t)+v(s,t)=f(s,t),(s,t)∈[0,1]×[0,T],v(s,0)=∂∂tv(s,0)=s2-s,v(0,t)=v(1,t)=0. where β=α-1f(s,t)=(s2-s)[t2-αΓ(3-α)+t]-2t.

Figure 5: Exact and approximate solution for Example 6.2 with Δt=0.025, β=0.50 and M = 40. (a) 3D plot for exact solution. (b) 3D plot for approximate solution

Figure 6: Absolute error for Example 6.2 when M = 40, β=0.75 and Δt=0.025

The exact solution is (s2 − s)t.

The numerical solution for Example 6.3, when α=1.50, 0≤s≤1, M = 20, Δt=0.01 is given by

V(s)={0.+s(-1.+s(1.+(3.75255×10-12-3.824×10-11s)s)),if s∈[0.00,0.05]-9.92262×10-16+s(-1.+s(1.+(1.15108×10-11-3.824×10-11s)s)),if s∈[0.05,0.10]-8.65974*10-15+s(-1.+s(1.+(1.90994*10-11-3.824*10-11s)s)),if s∈[0.10,0.15]-3.4639*10-14+s(-1.+s(1.+(2.67164*10-11-3.824*10-11s)s)),if s∈[0.15,0.20]⋮⋮-1.23634×10-12+s(-1.+s(1.+(6.51426×10-11-3.824×10-11)s)),if s∈[0.40,0.45]-1.93268×10-12+s(-1.+s(1.+(7.25322×10-11-3.824×10-11s)s)),if s∈[0.45,0.50]⋮⋮-1.77351×10-11+s(-1.+s(1.+(1.26306×10-10-3.824*10-11s)s)),if s∈[0.80,0.85]-2.24532×10-11+s(-1.+s(1.+(1.33952×10-10-3.824*10-11s)s)),if s∈[0.85,0.90]-2.79385×10-11+s(-1.+s(1.+(1.41483×10-10-3.824*10-11s)s)),if s∈[0.90,0.95]-3.4575×10-11+s(-1.+s(1.+(1.49157×10-10-3.824*10-11s)s)),if s∈[0.95,1.00]

The absolute numerical errors in RECBS solution to Example 6.3 using Δt=0.1, α=1.95 corresponding to different grid points are listed in Tab. 5. It is observed that our results are better than the localized kernel–based method (LKBM) [30]. Fig. 7 shows the physical behaviour of approximate solutions at different time levels when α=1.5, M = 100 and Δt=0.01. The 3D visuals of exact and numerical solutions with α=1.5, M = 100 and Δt=0.01 are shown in Fig. 8. Whereas, Fig. 9 depicts the absolute error between the exact and approximate solutions using α=1.5, M = 100 and Δt=0.01.

Table 5: Absolute error for Example 6.3 when Δt=0.1 for different values of s and α=1.95

Figure 7: Exact and numerical solution for Example 6.3 at different time levels when Δt=0.01, M = 100 and α=1.5

Figure 8: Exact and approximate solution for Example 6.3 with M = 100, Δt=0.01 and α=1.50. (a) 3D plot for exact solution. (b) 3D plot for approximate solution

Figure 9: Absolute error for Example 6.3 when M = 100, α=1.50 and Δt=0.01

Example 6.4

∂α∂tαv(s,t)+∂β∂tβv(s,t)-∂2∂s2v(s,t)=f(s,t),(s,t)∈[0,1]×[0,T],v(s,0)=∂∂tv(s,0)=0,v(0,t)=v(1,t)=0. where α=2β,f(s,t)=1/2t2[8π2tβ+(1+2t-βΓ(3-β)(Γ(3+β)))]sin(2πs).

The analytical solution is t2+β sin(2πs). The piecewise defined approximate solution for Example 6.4, when β=0.6, 0≤s≤1, M = 20, Δt=0.01 is given by

V(s)={-8.71974×10-18+6.28324s-41.6278s3+9.36996s4,if s∈[0.00,0.05]0.0000705375+s(6.27678+s(0.218315+s(-44.8655+27.1927s))),if s∈[0.05,0.10]0.00096495+s(6.23478+s(0.941467+s(-50.3082+42.3536s))),if s∈[0.10,0.15]0.00365344+s(6.14384+s(2.04345+s(-56.0616+53.3686s))),if s∈[0.15,0.20]⋮⋮-0.715608+s(13.5577+s(-26.2909+s(-9.51985+27.1927s))),if s∈[0.40,0.45]-1.47623+s(20.2526+s(-48.3868+s(22.8879+9.36996s))),if s∈[0.45,0.50]⋮⋮-5.49802+s(55.5206+s(-154.071+(157.413-53.3686s)s)),if s∈[0.80,0.85]0.77739+s(26.6075+s(-104.138+(119.106-42.3536s)s)),if s∈[0.85,0.90]11.1777+s(-19.1124+s(-28.7779+(63.9052-27.1927s)s)),if s∈[0.90,0.95]25.9746+s(-81.1203+s(68.6637+(-4.14797-9.36996s)s)),if s∈[0.95,1.00]

The comparison of L2 − norm for Example 6.4 using h = 5, Δt=1R, (R = 20, 40, 80) is reported in Tab. 6. It is found that our proposed algorithm has better accuracy when compared to ECBSM [36]. Tab. 7 shows the comparison of the calculated values of the order of convergence with proposed method for different values of spatial grid points M using β=0.75 and Δt=0.01. Fig. 10 shows the physical behaviour of numerical solutions at different time levels when β=0.5, M = 40 and Δt=0.01. The 3D visuals of exact and numerical solutions with β=0.5, M = 40 and Δt=0.01 are shown in Fig. 11, whereas, Fig. 12 depicts the absolute error between the exact and approximate solutions using β=0.5, M = 40 and Δt=0.01. Fig. 13 represents the behaviour of solution curve for different values of β.

Table 6: Absolute error norms for Example 6.4 using different values of M and β

Table 7: Experimental order of convergence for Example 6.4 using different values of M and Δt=0.01

Figure 10: Exact and numerical solution for Example 6.4 at different time levels when Δt=0.01, β=0.5 and M = 40

Figure 11: Exact and approximate solution for Example 6.4 with M = 20, Δt=0.01 and β=0.5. (a) 3D plot for exact solution. (b) 3D plot for approximate solution

Figure 12: Absolute error for Example 6.4 when M = 20, β=0.5 and Δt=0.01

Figure 13: Exact and numerical solutions for Example 6.4 with different values of β

7  Conclusion

This work is concluded with following remarks:

1.    An efficient algorithm based on a redefined set extended basis splines is proposed for numerical solution of multi-term time-fractional telegraph equation.

2.    The fractional time derivatives have been considered in the Caputo sense.

3.    The finite difference formulae have been used to discretize time-fractional derivatives while the discretization of spatial derivatives has been achieved by means of redefined extended B-spline functions.

4.    The spatial discretization used in this manuscript is superior to the other existing methods because the proposed method give continuous approximation with high accuracy to the solution curve of the unknown function and its derivatives at each and every point of the range of integration.

5.    The stability of presented algorithm has been proved along temporal grid.

6.    The theoretical results show that the accuracy of presented numerical approach in spatial direction is of order O(h2) whereas in time direction it is O(Δt2-α) when α=1+β and O(Δt1-α/2) when α=2β.

7.    The numerical rate of convergence is in the line with theoretical results.

8.    The comparison of error norms reveals that in terms of accuracy and straightforward implementation, the proposed algorithm performs better than the methods in [27,29,30,36].

Acknowledgement: We thank Dr. Nauman Khalid, Govt Post Graduate College, Faisalabad, Pakistan for his assistance in proofreading the manuscript.

Funding Statement: The author(s) received no specific funding for this study.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

 1.  Aguilar, J., Córdova-Fraga, T., Tórres-Jiménez, J., Escobar-Jiménez, R., Olivares-Peregrino, V. et al. (2016). Nonlocal transport processes and the fractional cattaneo-vernotte equation. Mathematical Problems in Engineering, 2016(1), 1–15. DOI 10.1155/2016/7845874.

 2.  Gómez-Aguilar, J., Escobar-Jiménez, R., López-López, M., Alvarado-Martnez, V. (2016). Atangana-baleanu fractional derivative applied to electromagnetic waves in dielectric media. Journal of Electromagnetic Waves and Applications, 30(15), 1937–1952. DOI 10.1080/09205071.2016.1225521.

 3.  Atangana, A. (2015). On the stability and convergence of the time-fractional variable order telegraph equation. Journal of Computational Physics, 293(3), 104–114. DOI 10.1016/j.jcp.2014.12.043.

 4.  Touchent, K. A., Belgacem, F. B. M. (2015). Nonlinear fractional partial differential equations systems solutions through a hybrid homotopy perturbation sumudu transform method. Nonlinear Studies, 22(4), 591–600.

 5.  Atangana, A., Bonyah, E. (2019). Fractional stochastic modeling: New approach to capture more heterogeneity. Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(1), 13118. DOI 10.1063/1.5072790.

 6.  Riaz, M. B., Atangana, A., Iftikhar, N. (2020). Heat and mass transfer in maxwell fluid in view of local and non-local differential operators. Journal of Thermal Analysis and Calorimetry, 96(1), 1–17. DOI 10.1007/s10973-020-09383-7.

 7.  Atangana, A., Bonyah, E., Elsadany, A. (2020). A fractional order optimal 4D chaotic financial model with mittag-leffler law. Chinese Journal of Physics, 65(2), 38–53. DOI 10.1016/j.cjph.2020.02.003.

 8.  Atangana, A., Hammouch, Z. (2019). Fractional calculus with power law: The cradle of our ancestors. European Physical Journal Plus, 134(9), 429. DOI 10.1140/epjp/i2019-12777-8.

 9.  Owolabi, K. M., Hammouch, Z. (2019). Mathematical modeling and analysis of two-variable system with noninteger-order derivative. Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(1), 13145. DOI 10.1063/1.5086909.

10. Owolabi, K. M., Hammouch, Z. (2019). Spatiotemporal patterns in the Belousov–Zhabotinskii reaction systems with atangana-baleanu fractional order derivative. Physica A: Statistical Mechanics and its Applications, 523(303), 1072–1090. DOI 10.1016/j.physa.2019.04.017.

11. Ghalib, M. M., Zafar, A. A., Hammouch, Z., Riaz, M. B., Shabbir, K. (2020). Analytical results on the unsteady rotational flow of fractional-order non-newtonian fluids with shear stress on the boundary. Discrete & Continuous Dynamical Systems-S, 13(3), 683–693. DOI 10.3934/dcdss.2020037.

12. Uçar, S., Uçar, E., Özdemir, N., Hammouch, Z. (2019). Mathematical analysis and numerical simulation for a smoking model with Atangana–Baleanu derivative. Chaos, Solitons & Fractals, 118, 300–306. DOI 10.1016/j.chaos.2018.12.003.

13. Ullah, S., Khan, M. A., Farooq, M., Hammouch, Z., Baleanu, D. (2019). A fractional model for the dynamics of tuberculosis infection using caputo-fabrizio derivative. Discrete & Continuous Dynamical Systems-S, 13(3), 975–993. DOI 10.3934/dcdss.2020057.

14. Asif, N., Hammouch, Z., Riaz, M., Bulut, H. (2018). Analytical solution of a maxwell fluid with slip effects in view of the caputo-fabrizio derivative. European Physical Journal Plus, 133(7), 272. DOI 10.1140/epjp/i2018-12098-6.

15. Singh, J., Kumar, D., Hammouch, Z., Atangana, A. (2018). A fractional epidemiological model for computer viruses pertaining to a new fractional derivative. Applied Mathematics and Computation, 316(30), 504–515. DOI 10.1016/j.amc.2017.08.048.

16. Weston, V., He, S. (1993). Wave splitting of the telegraph equation in R3 and its application to inverse scattering. Inverse Problems, 9(6), 789–812. DOI 10.1088/0266-5611/9/6/013.

17. Dehghan, M., Shokri, A. (2008). A numerical method for solving the hyperbolic telegraph equation. Numerical Methods for Partial Differential Equations: An International Journal, 24(4), 1080–1093. DOI 10.1002/num.20306.

18. El-Azab, M., El-Gamel, M. (2007). A numerical algorithm for the solution of telegraph equations. Applied Mathematics and Computation, 190(1), 757–764. DOI 10.1016/j.amc.2007.01.091.

19. Momani, S. (2005). Analytic and approximate solutions of the space-and time-fractional telegraph equations. Applied Mathematics and Computation, 170(2), 1126–1134. DOI 10.1016/j.amc.2005.01.009.

20. Dehghan, M., Yousefi, S., Lotfi, A. (2011). The use of He’s variational iteration method for solving the telegraph and fractional telegraph equations. International Journal for Numerical Methods in Biomedical Engineering, 27(2), 219–231. DOI 10.1002/cnm.1293.

21. Das, S., Vishal, K., Gupta, P., Yildirim, A. (2011). An approximate analytical solution of time-fractional telegraph equation. Applied Mathematics and Computation, 217(18), 7405–7411. DOI 10.1016/j.amc.2011.02.030.

22. Hayat, U., Mohyud-Din, S. (2013). Homotopy perturbation technique for time fractional telegraph equations. International Journal of Modern Theoretical Physics, 2(1), 33–41. DOI 10.1080/00207160902874653.

23. Wei, L., Dai, H., Zhang, D., Si, Z. (2014). Fully discrete local discontinuous galerkin method for solving the fractional telegraph equation. Calcolo, 51(1), 175–192. DOI 10.1007/s10092-013-0084-6.

24. Hosseini, V. R., Chen, W., Avazzadeh, Z. (2014). Numerical solution of fractional telegraph equation by using radial basis functions. Engineering Analysis with Boundary Elements, 38, 31–39. DOI 10.1016/j.enganabound.2013.10.009.

25. Srivastava, V. K., Awasthi, M. K., Tamsir, M. (2013). Rdtm solution of caputo time fractional-order hyperbolic telegraph equation. AIP Advances, 3(3), 32142. DOI 10.1063/1.4799548.

26. Wang, J., Zhao, M., Zhang, M., Liu, Y., Li, H. (2014). Numerical analysis of an H1-Galerkin mixed finite element method for time fractional telegraph equation. Scientific World Journal, 2014, 14. DOI 10.1155/2014/371413.

27. Modanli, M., Akgül, A. (2017). Numerical solution of fractional telegraph differential equations by theta-method. European Physical Journal Special Topics, 226(16–18), 3693–3703. DOI 10.1140/epjst/e2018-00088-6.

28. Xu, X., Xu, D. (2018). Legendre wavelets direct method for the numerical solution of time-fractional order telegraph equations. Mediterranean Journal of Mathematics, 15(1), 27. DOI 10.1007/s00009-018-1074-3.

29. Wang, Y., Mei, L. (2017). Generalized finite difference/spectral galerkin approximations for the time-fractional telegraph equation. Advances in Difference Equations, 2017(1), 281. DOI 10.1186/s13662-017-1348-2.

30. Uddin, M., Ali, A. (2018). On the approximation of time-fractional telegraph equations using localized kernel-based method. Advances in Difference Equations, 2018(1), 305. DOI 10.1186/s13662-018-1775-8.

31. Xu, G., Wang, G. Z. (2008). Extended cubic uniform B-spline and α-B-spline. Acta Automatica Sinica, 34(8), 980–984. DOI 10.1016/S1874-1029(08)60047-6.

32. Mohyud-Din, S. T., Akram, T., Abbas, M., Ismail, A. I., Ali, N. H. (2018). A fully implicit finite difference scheme based on extended cubic B-splines for time fractional advection-diffusion equation. Advances in Difference Equations, 2018(1), 109. DOI 10.1186/s13662-018-1537-7.

33. Benson, D. A., Schumer, R., Meerschaert, M. M., Wheatcraft, S. W. (2001). Fractional dispersion, lévy motion, and the made tracer tests. Transport in Porous Media, 42(1–2), 211–240. DOI 10.1023/A:1006733002131.

34. Khalid, N., Abbas, M., Iqbal, M. K., Baleanu, D. (2019). A numerical algorithm based on modified extended B-spline functions for solving time-fractional diffusion wave equation involving reaction and damping terms. Advances in Difference Equations, 2019(1), 378. DOI 10.1186/s13662-019-2318-7.

35. Wasim, I., Abbas, M., Amin, M. (2018). Hybrid B-spline collocation method for solving the generalized burgers-fisher and burgers-huxley equations. Mathematical Problems in Engineering, 2018(10), 1–18. DOI 10.1155/2018/6143934.

36. Akram, T., Abbas, M., Ismail, A. I., Ali, N. H. M., Baleanu, D. (2019). Extended cubic B-splines in the numerical solution of time fractional telegraph equation. Advances in Difference Equations, 2019(1), 365. DOI 10.1186/s13662-019-2296-9.

[picture] | This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.