# A B-Spline Quasi Interpolation Crank–Nicolson Scheme for Solving the Coupled Burgers Equations with the Caputo–Fabrizio Derivative

Authors: M. Taghipour; H. Aminikhah
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, a Crank–Nicolson finite difference scheme based on cubic B-spline quasi-interpolation has been derived for the solution of the coupled Burgers equations with the Caputo–Fabrizio derivative. The first- and second-order spatial derivatives have been approximated by first and second derivatives of the cubic B-spline quasi-interpolation. The discrete scheme obtained in this way constitutes a system of algebraic equations associated with a bi-pentadiagonal matrix. We show that the proposed scheme is unconditionally stable. Numerical examples are provided to verify the efficiency of the method.


---

## Body

## 1. Introduction

The coupled Burgers equations are coupled partial differential equations which are capable of describing realitic polydispersive supensions. The coupled Burgers equation perdicts an interesting phenomenon, which is called phase shifts [1]. This equation is one of the fundamental models in fluid mechanics and arises in gas dynamics, chromatography, and flood waves in rivers [2]. The coupled viscous Burges equation is given by(1)ut−uxx+ηuux+α1uvx=0,x∈a,b,t∈0,T,vt−vxx+ηvvx+α2uvx=0,x∈a,b,t∈0,T,with the initial conditions(2)ux,0=ϕ1x,vx,0=ϕ2x,and the boundary conditions(3)ua,t=f1a,t,ub,t=f2b,t,va,t=g1a,t,vb,t=g2b,t,where η is a real constant and α and β are arbitrary constants depending on the system parameters such as Peclet number, Stokes velocity of particles due to gravity, and Brownian diffusivity [3].Spline is a special function defined piecewise by polynomials. The spline approximation first appeared in a paper by Schoenberg [4]. Spline interpolation is a form of interpolation where the interpolant is a special type of piecewise polynomial called a spline. Applications of spline function in fractional partial differential equations can be found in [5–15].In this article, we consider the following coupled Burgers equation with time fractional derivative:(4)∂γu∂tγ−∂2u∂x2+ηu∂u∂x+α1∂uv∂x=q1,x∈a,b,t∈0,T,0<γ<1,(5)∂γv∂tγ−∂2v∂x2+ηv∂v∂x+α2∂uv∂x=q2,x∈a,b,t∈0,T,0<γ<1,with initial(6)ux,0=ϕ1x,vx,0=ϕ2x,and the boundary conditions(7)ua,t=f1a,t,ub,t=f2b,t,va,t=g1a,t,vb,t=g2b,t,where γ is order of time fractional derivative. Also, η, α1, and α2 are those ones we said before.D0CFtγux,t=∂γux,t/∂tγ denotes the Caputo–Fabrizio derivative of the function ux,t defined as [16](8)D0CFtγux,t=Mγ1−γ∫0tu′x,se−σt−sds,where Mγ is a normalization function such that M0=M1=1 and σ=γ/1−γ.Recently, the Caputo–Fabrizio derivative has received more attention from researchers due to their description of some physical phenomenon [17–26].For the description of memory and some physical properties of various materials and processes, modeling with fractional derivatives is very appropriate. This is the main benefit of fractional derivatives in comparison with classical integer order models, in which such effects are missed. In recent years, the coupled system of Burgers equations with fractional derivatives has been the focus of attention. For example, in [27], the Adomian decomposition method is directly extended to study the coupled Burgers equations with time and space fractional derivatives. Khan et al. [28] proposed the generalized differential transform method (GDTM) and homotopy perturbation method (HPM) for time fractional Burgers and coupled Burgers equations. The fractional variational iteration method (FVIM) to solve a time and space fractional coupled burgers equations is given by Prakash et al. [29]. In [30], a q-homotopy analysis transform method (q-HATM) for time and space fractional coupled Burgers equations is introduced. Aminikhah and Malekzadeh [31] introduced a new homotopy perturbation method for system of variable coefficient coupled Burgers equations with time fractional derivative. In [32], the Laplace-Adomian decomposition method (LADM), the Laplace-variational iteration method (LVIM), and the reduced differential transform method (RDTM) are proposed to solve the one- and two-dimensional fractional coupled Burgers equations. Albuohimad and Adibi derived a hybrid spectral exponential Chebyshev method (HSECM) for time fractional coupled Burgers equations [33]. In [34], authors investigate the fractional coupled viscous Burgers equation involving Mittag–Leffler kernel. In [35], the generalized two-dimensional differential transform method (DTM) was applied to solve the coupled Burgers equations with space and time fractional derivatives. Ozdemir et al. used the Gegenbauer wavelets-based computational methods to find the approximate solutions of the coupled system of Burgers equations with time fractional derivative [36].Our aim is to propose a Crank–Nicolson finite difference scheme using cubic B-spline quasi-interpolation to solve time fractional coupled viscous Burgers equations. The first- and second-order spatial derivatives have been approximated by first and second derivatives of the cubic B-spline quasi-interpolation. This approximations have not been used for the fractional coupled Burgers equations before.The paper is organized as follows. In Section2, we present some basic definitions and concepts of quasi-interpolants. In Section 3, using the quasi-interpolant and Crank–Nicolson finite difference method, we obtain a numerical scheme. The stability of this method is studied in Section 4. In Section 5, some numerical examples are proposed. Finally, conclusions are given in Section 6.


## 2. Univariate Spline Quasi-Interpolants

In this section, we introduce the basic concepts about B-spline and univariate B-spline quasi-interpolants that we will use in Section3.According to [37], let(9)Pd1≔space of univariate polynomials of degree at mostd,and Ω=a,b be an interval that has been partitioned into subintervals via a set of points Δ=xii=0k+1 with(10)a=x0<x1<⋯<xk<xk+1=b.We define the space of univariate polynomial splines of smoothnessr and degree d with knots Δ as(11)SdrΔ≔s∈CrΩ:s|xi,xi+1∈Pd1,i=0,…,k,where 0≤r<d are given integers. We have(12)n≔dimSdrΔ=kd−r+r+1.For a formal proof of this fact, see Theorem 4.4 of [38].Given0≤r<d and Δ=xii=0k+1, the associated extended partition Δe is defined to be yii=0n+d+1, where n is the dimension of SdrΔ given in (11):(13)a=y1=⋯=yd+1,yn+1=⋯=yn+d+1=b,yd+2≤⋯≤yn=x1,…,x1︷d−r,…,xk,…,xk︷d−r.Given an extended partitionΔe, let(14)Qi1t≔1yi+1−yi,yi≤t<yi+1,0,otherwise,for i=1,…,n+d, and let(15)Qimt≔t−yiQim−1t+yi+m−tQi+1m−1tyi+m−yi,yi≤t<yi+m,0,otherwise,for 2≤m≤d+1 and i=1,…,n+d−m+1. Let(16)Nimt≔yi+m−yiQimt,i=1,…,n+d−m+1.We call these the normalized B-splines of orderm (or degree m−1) associated with the extended partition Δe.In [37], univariate B-spline quasi-interpolants can be defined as a formula of the form(17)Qdfx=∑i=1n+dλifNix,where Nii=1n+d are the B-splines forming a basis of SdrΔ.Quasi-interpolants have been heavily studied in the literature. Some basic ideas and sources for further information can be found in [38]. For a good approximations, we need to make sure it reproduces polynomials, i.e., Qp=p for all p∈Pd1. For each i=1,…,n, we assume that the coefficient λi is a linear functional defined on Ca,b that can be computed from samples of f at some set of points σλi in a,b.According to [39], the error of a quasi-interpolation satisfies(18)fx−Qdfx≤Qdd+1!fd+1∞,Dxhxd+1,x∈Dyd,where Dyd=yd+1,yn+1, Dx is the union of the supports of all B-splines Ni, i∼x and fd+1∞,Dx denote the maximum norm of fd+1 on Dx and hx=maxy∈Dxy−x that ∼ is used to indicate proportionality. If the local mesh ratio is bounded, i.e., if the quotients of the lengths of adjacent knot intervals are ≤ry, then the error of the derivatives on the knot intervals yl,yl+1 can be estimated by(19)fjx−Qdfjx≤cd,ryQdfd+1∞,Dxhxd+1−j,for j≤d.Supposea=t0<⋯<tn=b are equally spaced points in the interval a,b. Let(20)λif≔ft0,i=1,1187ft0+18ft1−9ft2+2ft3,i=2,16−fti−3+8fti−2−fti−1,3≤i≤n+1,1182ftn−3−9ftn−2+18ftn−1+7ftn,i=n+2,ftn,i=n+3.Then, (17) defines a linear operator mapping Ca,b into S32Δ with Qp=p for all cubic polynomials p. For approximate derivatives of f by derivatives of Q3f up to the order h3, we can evaluate the value of f′ and f〞 at xi by Q3f′x=∑i=1n+3λifNi′x and Q3f〞x=∑i=1n+3λifNi〞x. We set Y=f0,f1,…,fnT, Y′=f0′,f1′,…,fn′T, and Y〞=f0〞,f1〞,…,fn〞T, where fj′=Q3f′xj, j=1,…,n, and fj〞=Q3f〞xj, j=1,…,n. By solution of the linear systems,(21)fi′=∑i=1n+3λifNi′x,i=0,1,…,n,fi″=∑i=1n+3λifNi″x,i=0,1,…,n,we obtain(22)Y′=1hD1Y,Y″=1h2D2Y,where D1,D2∈ℝn+1×n+1 and is obtained as follows:(23)D1=−1163−321300…00−13−121−1600…00112−23023−1120…000112−23023−112…00⋮⋮⋮⋮⋮⋮⋮⋮⋮00…112−23023−112000…0112−23023−11200…0016−1121300…00−1332−3116,D2=2−54−100…001−21000…00−1653−353−160…000−1653−353−16…00⋮⋮⋮⋮⋮⋮⋮⋮⋮00…−1653−353−16000…0−1653053−1600…0001−2100…00−14−52.


## 3. Numerical Scheme

We consider a gridxi=a+ih,i=0,1,…,M, with h=xi+1−xi. The step length in time is denoted by τ and tk=kτ,τ=T/N,0≤k≤N. So, the domain a,b×0,T is divided into a uniform grid of mesh points xj,tk. The values of the function u at the grid points are denoted ujk=uxj,tk and Ujk is the approximate solution at the point xj,tk.A discrete approximation to theD0CFtγux,t at xj,tk can be obtained by the following approximation [40]:(24)D0CFtγuxj,tk=Mγ1−γ∫0tk∂uxj,s∂se−σtk−sds=Mγ1−γ∑l=1k∫tk−1tk∂uxj,s∂se−σtk−sds=Mγ1−γ∑l=1k∫tk−1tkujl−ujl−1τe−σtk−sds+R=Mγγτ∑l=1kujl−ujl−1wk+R=Mγγτw0ujk−∑l=1k−1wk−l−1−wk−lujl−wk−1uj0+R,where(25)wj=e−στj−e−στj+1.Theorem 1.
Supposevt∈C20,tk. Let(26)A=∑l=1k∫tl−1tlv′τ−vtl−vtl−1Δte−σtk−τdτ.
Then,(27)A≤Δt2σmax0≤t≤tkv〞tσ212−σ3Δt24+⋯tk.Proof.
Using the Taylor series expansion with integral remainder, we have(28)A=1Δt∑l=1k∫tl−1tl∫tl−1τv〞ss−tl−1ds−∫τtlv〞stl−sdse−σtk−τdτ,=1σΔt∑l=1k∫tl−1tls−tl−1e−σtk−tl−e−σtk−s−tl−se−σtk−s−e−σtk−tl−1v〞sds.
Since,(29)∫tl−1tls−tl−1e−σtk−tl−e−σtk−s−tl−se−σtk−s−e−σtk−tl−1ds,=Δt4e−σΔtk−lσ212−σ3Δt24+⋯≤Δt4σ212−σ3Δt24+⋯,and e−σΔtk−l<1 for l=1,2,…,k and tk=kΔt; hence, the result will be achieved.
Now, using Theorem1, we obtain(30)Mγ1−γ∑l=1k∫tl−1tl∂u∂s|xj,s−ujl−ujl−1τe−σtk−sds,≤Mγ1−γτ2σmax0≤t≤tku〞tσ212−σ3τ24+⋯∑l=1ke−στk−l=Oτ2.
We introduce some lemmas which will be used in numerical scheme and stability analysis.Lemma 1.
(see [41]). Suppose ut∈C30,tk+1; then, we have(31)D0CFtγuk+1/2−Ttγuk+1/2=Oτ2,0≤k≤N−1,where(32)Ttγuk+1/2=1γτw0ujk+ujk+12−∑l=1kwk−l−wk−l+1ujl+ujl−12−wkuj0,0≤k≤N−1.Lemma 2.
(see [42]). For the definition Mj, we have Mj>0 and Mj+1<Mj,∀j≤k.Lemma 3.
(see [42]). Suppose ut∈Cx,t4,40,L×0,T; it holds that 0≤wj≤Cτ, and 0≤wj−wj+1≤Cτwj.Now, we present the numerical scheme for solving (4)–(7) based on the Crank–Nicolson method and cubic B-spline quasi-interpolant. We approximate equations (4) and (5) at the point xj,tk+1/2 by the Crank–Nicolson finite difference approximation:(33)∂γuxj,tk+1/2∂tγ−12uxxxj,tk+1+uxxxj,tk+η2uxj,tk+1uxxj,tk+1+uxj,tkuxxj,tk+α12uvxjk+1+uvxjk=q1xj,tk+1/2,(34)∂αvxj,tk+1/2∂tγ−12vxxxj,tk+1+vxxxj,tk+η2vxj,tk+1vxxj,tk+vxj,tk+1vxxj,tk+α22uvxjk+1+uvxjk=q2xj,tk+1/2.The nonlinear terms in equations (33) and (34) are linearized using the following quasi-linearization [43]:(35)uuxjk+1=ujk+1uxjk+ujkuxjk+1−uuxjk+Oτ2,(36)uvxjk+1=ujk+1vxjk+ujkvxjk+1−uvxjk+Oτ2,(37)uxvjk+1=uxjk+1vjk+uxjkvjk+1−uxvjk+Oτ2,(38)vvxjk+1=vjk+1vxjk+vjkvxjk+1−vvxjk+Oτ2.Now, using (22), (24), Lemma 1, and relations (36)–(38), we have the following difference scheme which is accurate of the order Oτ2+h2:(39)w02+γτη2h∑l=0MDjl1Ulk+α1γτ2h∑l=0MDjl1VlkUjk+1−γτ2h2∑l=0MDjl2Ulk+1+γτη2hUjk∑l=0MDjl1Ulk+1+α1γτ2hVjk∑l=0MDjl1Ulk+1+α1γτ2hVjk+1∑l=0MDjl1Ulk+α1γτ2hUjk∑l=0MDjl1Vlk+1=γτq1jk+1/2−w02Ujk+γτ2h2∑l=0MDjl2Ulk+wkUj0+∑l=0Mwk−l−wk−l+1Ujl+Ujl−12,(40)w02+γτη2h∑l=0MDjl1Vlk+α2γτ2h∑l=0MDjl1UlkVjk+1−γτ2h2∑l=0MDjl2Vlk+1+γτη2hVjk∑l=0MDjl1Vlk+1+α2γτ2hUjk∑l=0MDjl1Vlk+1+α2γτ2hUjk+1∑l=0MDjl1Vlk+α2γτ2hVjk∑l=0MDjl1Ulk+1=γτq2jk+1/2−w02Vjk+γτ2h2∑l=0MDjl2Vlk+wkVj0+∑l=0Mwk−l−wk−l+1Vjl+Vjl−12.Therefore, in each time step we solve the following bi-pentadiagonal linear system of dimension2M−2×2M−2:(41)AMk+1=Fk,where M=U1k+1,…,UM−1k+1,V1k+1,…,VM−1k+1′.


## 4. Stability Analysis

To study the stability analysis of the proposed scheme, we use the Fourier method. In applying the Fourier stability method, the nonlinear terms are temporarily frozen, since the stability analysis is strictly only applicable to linear equations. Thus, we have linearized the nonlinear termsuux and uvx in equation (4) by freezing u and v as a local constants β1 and β2, respectively. We have(42)D0CFtαUk+1/2−12Uxxjk+1+Uxxjk+ηβ12Uxjk+1+Uxjk+α1β22Uxjk+1+Uxjk+α1β12Vxjk+1+Vxjk=q1xj,tk+1/2.Substituting approximations (22) and (32) yield the following difference equation:(43)w02Ujk+1−γτ2h2∑l=0MDjl2Ulk+1+γτ2hηβ1+α1β2∑l=0MDjl1Ulk+1+γτα1β12h∑l=0MDjl1Vlk+1=−w02Ujk+∑l=0Mwl−k−wl−k+1Ujl+Ujl−12+wkUj0+γτ2h2∑l=0MDjl2Ulk−γτ2hηβ1+α1β2∑l=0MDjl1Ulk−γτα1β12h∑l=0MDjl1Vlk+q1jk+1/2.LetU˜jk and V˜jk be the approximate solutions of (43), and define ζjk=Ujk−U˜jk,ξjk=Vjk−V˜jk,1≤j≤M−1,0≤k≤N−1, with corresponding vectors:(44)ζk=ζ1k,ζ2k,…,ζM−1kT,ξk=ξ1k,ξ2k,…,ξM−1kT.So, we have(45)w0ζjk+12−γτ2h2−ζj−2k+16+5ζj−1k+13−3ζjk+1+5ζj+1k+13−ζj+2k+16+γτν2ζj−2k+112−2ζj−1k+13+2ζj+1k+13−ζj+2k+112+γτα1β12h112ξj−2k+1−23ξj−1k+1+23ξj+1k+1−112ξj+2k+1=−w02ζjk+∑l=0Mwl−k−wl−k+1ζjl+ζjl−12+γτ2h2−16ζj−2k+53ζj−1k−3ζjk+53ζj+1k−16ζj+2k−γτν2112ζj−2k−23ζj−1k+23ζj+1k−112ζj+2k−γτα1β12h112ξj−2k−23ξj−1k+23ξj+1k−112ξj+2k+wkζj0,where ν=ηβ1+α1β2.Now, we define the grid functions as follows:(46)ζkx=ζjk,xj−h2<x≤xj+h2,0,0≤x≤h2orL−h2<x≤L,ξkx=ξjk,xj−h2<x≤xj+h2,0,0≤x≤h2orL−h2<x≤L.We expand theζkx and ξkx into the following Fourier series expansions:(47)ζkx=∑l=−∞∞Aklei2πlx/L,ξkx=∑l=−∞∞Bklei2πlx/L,where(48)Akl=1L∫0Lζkxe−i2πlx/Ldx,Bkl=1L∫0Lξkxe−i2πlx/Ldx.Applying the Parseval equality,(49)∫0Lζkx2dx=∑l=−∞∞Akl2,∫0Lξkx2dx=∑l=−∞∞Bkl2,(50)∫0Lζkx2dx=∑j=1M−1hζjk2,∫0Lξkx2dx=∑j=1M−1hξjk2,we have(51)ζk22=∑l=−∞∞Akl2,ξk22=∑l=−∞∞Bkl2.Now, we suppose that(52)ζjk=Akeiσxjh,ξjk=Bkeiσxjh,where σx=2lπ/L. Substituting the above relations into (45) leads to(53)w02+γτ6h2cos2σxh−5γτ3h2cosσxhAk+1+i−γτα1β112hsin2σxh+2γτα1β13hsinσxhBk+1+i−ν6sin2σxh+4ν3sinσxhAk+1=−w02−γτ6h2cos2σxh+5γτ3h2cosσxh−3γτ2h2Ak+wkA0+iν6sin2σxh−4ν3sinσxhAk+iγτα1β112hsin2σxh+2γτα1β13hsinσxhBk+∑l=1kwk−l−wk−l+1Al+Al−12.Set(54)X=w02+γτ6h2cos2σxh−5γτ3h2cosσxh,Y=−ν6sin2σxh+4ν3sinσxh,Z=−γτα1β112hsin2σxh+2γτα1β13hsinσxh.We have(55)X+iYAk+1+iZBk+1=−X+iYAk−iZBk+∑l=1kwk−l−wk−l+1Al+Al−12+wkA0so that(56)Ak+1≤Ak+iZX+iYBk+1+Bk+∑l=1kwk−l−wk−l+1Al+Al−121X+iY+wkA0X+iY.Theorem 2.
IfAk is the solution of equation (55), then there are positive constants Ck such that(57)Ak≤CkA0,k=1,2,…,N−1.Proof.
We use the mathematical induction for proof. Fork=1, we have(58)A1≤A0+iZX+iYB1+B0+w0A0X+iY.
Using the convergence of the series on the right-hand side of equation (49), we know that there exists a positive constant P2 such that(59)Bk≤P2A0,k=0,1,…,N−1.
So,(60)A1≤A0+DA0+EA0≤C1A0.
Now, suppose that(61)Ak≤CkA0,k=1,2,…,N−2.
By equation (56), we have(62)Ak+1≤CkA0+FC2A0+G∑l=1kwk−l−wk−l+1ClA0+ClA02+HA0.
Now, assume that(63)C′=maxC1,C2,…,CN−2.
So,(64)Ak+1≤CkA0+FP2A0+GC′A0+HA0≤Ck+1A0.Remark 1.
From (5), with similar way, there are positive constants Qk such that(65)Bk≤QkB0,k=1,2,…,N−1.Theorem 3.
The finite difference schemes (39) and (40) are unconditionally stable for γ∈0,1.Proof.
According to Theorem2 and Remark 1, using (51), we obtain(66)ζk22=∑l=−∞∞Akl2≤∑l=−∞∞Ck2A0l2=Ck2ζ022,ηk22=∑l=−∞∞Bkl2≤∑l=−∞∞Qk2B0l2=Qk2η022,so that(67)Uk−U˜k22≤CkU0−U˜022,Vk−V˜k22≤QkV0−V˜022,which shows that schemes (39) and (40) are unconditionally stable.


## 5. Numerical Results

In this section, we provide two examples to illustrate efficiency of schemes (39) and (40). All experiments are performed on a Windows 10 (64 bit) Intel(R) Core(TM) i7-7500U CPU 2.70 GHz, 8.0 GB of RAM using MATLAB R2017b. In all examples, we use the error norm(68)eτ,h=eN=h∑j=1MejN21/2,where ejn=ujn−Ujn and ejn=vjn−Vjn. We evaluate the convergence order with the following formula:(69)rτ,h=log2eτ,2heτ,h.Example 1.
Consider the coupled Burgers equation with time fractional derivative with exact solutionsux,t=vx,t=etx21−x2:(70)∂γu∂tγ−∂2u∂x2+ηu∂u∂x+α1∂uv∂x=q1,x∈0,1,t∈0,1,0<γ<1,∂γv∂tγ−∂2v∂x2+ηv∂v∂x+α2∂uv∂x=q2,x∈0,1,t∈0,1,0<γ<1,u0,t=0,v0,t=0,u1,t=0,v1,t=0,ux,0=x21−x2,vx,0=x21−x2.
We derive the functionsq1 and q2 with the help from exact solutions. Numerical solution and exact solution have been demonstrated in Figure 1. Tables 1–4 give the approximation errors and CPU times for the difference schemes. We choose different space step sizes to obtain the numerical results and order of convergence. Figure 2 shows the comparison of numerical solution and exact solution for γ=0.9 at t=1 and contour plot of numerical solution. Also, Figure 3 shows the pointwise errors for ux,t.Figure 1
The plot of exact solution (a) and numerical solution (b) atτ=1/1000 and h=1/160 with γ=0.9 for Example 1.
(a)(b)Table 1
Takeτ=0.001 and η=800andα1=0.01,α2=0.01 experiment order of convergence for ux,t at γ=0.1,0.3,0.5 for Example 1.
hγ=0.1γ=0.3γ=0.5eNOrderCPUeNOrderCPUeNOrderCPU1/101.1161e − 031.6557 s1.1133e − 031.4701 s1.1094e − 031.5101 s1/203.8484e − 042.59412.8074 s1.8406e − 042.59662.7374 s1.8282e − 042.60132.8795 s1/404.7604e − 051.95715.4495 s4.7393e-051.95745.2651 s4.7059e − 051.95795.6347 s1/801.1883e − 052.002210.6861 s1.1833e − 052.001810.6804 s1.1756e − 052.001112.0549 s1/1602.7784e − 062.096622.2362 s2.7680e − 062.001823.96062.7529e − 062.094424.3564 sTable 2
Takeτ=0.001 and η=800andα1=0.01andα2=0.01 experiment order of convergence for ux,t at γ=0.7,0.9 for Example 1.
hγ=0.7γ=0.9eNOrderCPUeNOrderCPU1/101.1018e − 031.6541 s1.0653e − 031.7714 s1/201.8045e − 042.61023.1754 s1.7495e − 042.60623.4511 s1/404.6426e − 051.95866.2369 s4.5085e − 051.95626.7489 s1/801.1614e − 051.999113.0887 s1.1364e − 051.988215.1479 s1/1602.7286e − 062.089625.7009 s2.7292e − 062.057931.5002 sTable 3
Takeτ=0.001 and η=200andα1=0.03andα2=0.03 experiment order of convergence for vx,t at γ=0.1,0.3,0.5 for Example 1.
hγ=0.1γ=0.3γ=0.5eNOrderCPUeNOrderCPUeNOrderCPU1/101.1129e − 031.6208 s1.1110e − 031.3725 s1.1061e − 031.5264 s1/201.8595e − 042.58131.5632 s1.8514e − 042.58522.6589 s1.8385e − 042.58892.5895 s1/404.7628e − 051.96505.2549 s4.7411e − 051.96535.2359 s1.1695e − 052.00885.3654 s1/801.1827e − 052.009710.6584 s1.1775e − 052.009510.2545 s1.1695e − 052.008812.0036 s1/1602.7407e − 062.109522.12358 s2.7296e − 062.109023.8547 s2.7131e − 062.107924.1354 sTable 4
Takeτ=0.001 and η=200andα1=0.03andα2=0.03, experiment order of convergence for vx,t at γ=0.7,0.9 for Example 1.
hγ=0.7γ=0.9eNOrderCPUeNOrderCPU1/101.0983e − 031.6151 s1.0614e − 031.6646 s1/201.8136e − 042.59833.0172 s1.7506e − 042.60003.3234 s1/404.6406e − 051.96656.1189 s4.4958e − 051.96126.7189 s1/801.1544e − 052.007213.5179 s1.1262e − 051.997115.1345 s1/1602.6855e − 062.103925.4015 s2.7119e − 062.075531.4012 sFigure 2
The comparison (a) and absolute error (b) between numerical solution and exact solution withτ=1/1000, h=1/160, and γ=0.9 at t=1 for Example 1.
(a)(b)Figure 3
Pointwise errors forux,t with τ=1/1000, h=1/160, and γ=0.9 (a) and contour plot of numerical solution ux,t (b) for Example 1.
(a)(b)Example 2.
Consider the following coupled Burgers equation with time fractional derivative with exact solutionsux,t=vx,t=x1−xsint:(71)∂γu∂tγ−∂2u∂x2+2u∂u∂x+α1∂uv∂x=q1,x∈0,1,t∈0,T,0<γ<1,∂γv∂tγ−∂2v∂x2+2v∂v∂x+α2∂uv∂x=q2,x∈0,1,t∈0,T,0<γ<1,u0,t=0,v0,t=0,u1,t=0,v1,t=0,ux,0=0,vx,0=0.
We solve this problem with the method developed in this article with several values ofT and γ. Graphs of numerical solution and exact solution at different times have been demonstrated in Figure 4. Figure 4 shows that the proposed method is efficient. Table 5 gives the approximation errors for t=0.5,1,1.5,2,2.5,3,4 with different α1,α2 and CPU times. We choose h=1/160 for space step size to obtain the numerical results. Figure 5 shows the comparison of numerical solution and exact solution for γ=0.1 at t=1. Also, Figure 6 shows the pointwise errors and contour plot of numerical solution for ux,t.Figure 4
Graphs of numerical and exact solution at different times withτ=0.001 for Example 2.Table 5
Errors forux,t and vx,t at different times for Example 2.
Tα1α2γ=0.1γ=0.3γ=0.7euNevNCPUeuNevNCPUeuNevNCPU0.50.010.012.2750e −071.0966e −0819.7143s2.2269e −071.0692e −0819.6825 s2.1196e −072.3957e −0922.1880 s10.030.022.1872e −061.9266e −0819.7329 s2.1574e −061.8640e −0819.9076 s2.0966e −065.8427e −0924.0896 s1.50.0010.058.1009e −082.2678e −0819.6884 s8.1626e −082.1703e −0820.8275 s9.8131e −086.6570e −0925.0396 s20.0010.0016.6249e −082.0562e −0819.8765 s6.8923e −081.9294e −0821.4610 s9.1599e −083.4789e −0925.4790 s2.50.0020.0026.3910e −081.3390e −0819.9632 s7.1535e −081.1951e −0822.1602 s1.0437e −073.7207e −0925.9242 s30.010.013.3597e −082.8784e −0919.9462 s6.9805e −081.4200e-0922.4854 s1.3465e −071.2871e −0826.0649 s40.010.016.2662e −071.7635e −0820.2347 s6.3810e −071.8725e −0823.4227 s6.0534e −072.9270e −0826.4143 sFigure 5
The plot of exact solution (a) and numerical solution (b) atτ=1/1000 and h=1/160 with γ=0.1 for Example 2.
(a)(b)Figure 6
Pointwise errors forvx,t with τ=1/1000, h=1/160, and γ=0.1 (a) and contour plot of numerical solution ux,t (b) for Example 2.
(a)(b)


## 6. Conclusion

In this article, we constructed a Crank–Nicolson finite difference scheme based on cubic B-spline quasi-interpolation to solve time fractional coupled Burgers equations. By the Fourier series method, we proved that this scheme is unconditionally stable. Numerical examples have been carried out to show the convergence orders and applicability of the scheme and error norms are calculated with respect to different space step sizes. From the error tables and graphs of exact and numerical solution, we can say that our method has a good accuracy. For the numerical computations, we have used Matlab.



---

*Source: 8837846-2021-01-23.xml*