# An Efficient Explicit Decoupled Group Method for Solving Two–Dimensional Fractional Burgers’ Equation and Its Convergence Analysis

Authors: N. Abdi; H. Aminikhah; A. H. Refahi Sheikhani; J. Alavi; M. Taghipour
Journal: Advances in Mathematical Physics
 (2021)

---

## Abstract

In this paper, the Crank–Nicolson (CN) and rotated four-point fractional explicit decoupled group (EDG) methods are introduced to solve the two-dimensional time–fractional Burgers’ equation. The EDG method is derived by the Taylor expansion and 45° rotation of the Crank–Nicolson method around thex and y axes. The local truncation error of CN and EDG is presented. Also, the stability and convergence of the proposed methods are proved. Some numerical experiments are performed to show the efficiency of the presented methods in terms of accuracy and CPU time.


---

## Body

## 1. Introduction

Fractional calculus is a generalization of the integration and derivation of integer order operators to fractional order that allows us to describe a real system more accurately than integers. Although the fractional order of a real system may be low, it is yet considered as a fractional system. An important feature of fractional calculus is its nonlocality. The fractional derivative (and integrals) of a function is given by a definite integral, thus, it depends on the value of the function over the entire interval [1]. Researchers confirm the existence of interesting phenomena in nature, which cannot be modeled by classical differential equations. To cope with this problem, the nonlocality property of fractional derivative could be a beneficial tool to study our considered system. Fractional calculus has recently been used in various scientific and engineering fields [2–4]. Some fractional calculus applications in modeling and design of control systems are introduced in [5]. The most recent developments and trends in the use of fractional calculus in biomedicine and biology are presented by Ionescu et al. [6]. Based on the fractional calculate, Tang et al. [7] proposed a new four-element creep model; this model accords well with the experimental data of Changshan rock salt. Fractional calculus has an extraordinary potential in signal denoising [8]. Gong et al. discussed the generation conditions of chaotic behavior and proposed the adaptive synchronization control method for a class of fractional-order financial system [9]. Numerous definitions of fractional derivative have been introduced in the literature, amongst are Riemann–Liouville, Caputo, and Caputo–Fabrizio [10]. The Caputo–Fabrizio fractional derivative on the contrary of other derivatives has a nonsingular kernel. Hence, it has been considered by many researchers [11–15]. Since to obtain the exact solution of fractional differential equations is very difficult and sometimes impossible, it is usually approximated by a numerical method such as finite difference method [16–18], finite volume methods [19], and spectral method [20].Burgers’ equation as a nonlinear partial differential equation is widely used in various areas such as fluid mechanics, gas dynamics, and traffic flow which combine nonlinear propagation effects with diffusive one.In this paper, we consider the following two dimensional time-fractional Burgers’ equation:(1)∂αu∂tα+u∂u∂x+∂u∂y=ν∂2u∂x2+∂2u∂y2+fx,y,t,x,y∈Ω,0<t<T,with the initial condition(2)ux,y,0=u0x,y,x,y∈Ω,and the boundary condition(3)ux,y,t=0,x,y∈∂Ω,0≤t≤T,whereν=1/Re, Re is Reynolds number characterizing the strength of viscosity, Ω=0,1×0,1, 0<α<1, and the term ∂αu/∂tα denotes the α order Caputo-Fabrizio fractional derivative of the function ux,y,t defined as:
(4)∂αux,y,t∂tα=Mα1−α∫0t∂ux,y,s∂se−α1−αt−sds,whereMα is a normalization function such that M0=M1=1.We apply the finite difference method to solve Equations (1)–(3). In finite difference, to find the value corresponding to each grid point, it is used natural ordering (indexing the grid of point from left to right and bottom to top by point to point) or group to group. There are several different methods to order the interior mesh point such as natural ordering, diagonal ordering, and alternating diagonal ordering [21]. Different linear systems would be produced by different arrangements of the grid points.We propose two finite difference schemes. The first scheme is given by the Crank–Nicolson difference method (by natural ordering). In this scheme, to obtain a more accurate numerical solution, we should use a smaller mesh size, which requires more storage space and computing time. In order to fix this problem and accelerate the convergence, we use the explicit decoupled group (EDG) method introduced by Abdullah in 1990 [22]. Many studies have been done in reference to the EDG method for examples [23–26]. The EDG method is based on rotating the Crank–Nicolson difference scheme and group ordering of grid points. Applying the EDG method to the Crank–Nicolson difference scheme result in a new scheme in which the dimension of the system is half of the dimension of the system generated by the Crank–Nicolson difference scheme. On this account, half of the grid points are obtained and the rest can be calculated directly. Consequently, the EDG method can be favourably used to reduce the computational cost. In addition, it is worth to notice that we can take advantage of parallel computers to run it.The rest of the paper is organized as follows. In Section2, the Crank–Nicolson difference scheme will be applied to Equations (1)–(3), and also, we give the truncation error. In Section 3, we describe the formulation of the EDG method. The stability of these schemes is discussed in Section 4. In Section 5, we analyze the convergence of these schemes. Numerical examples are carried out to verify the high efficiency of our method in Section 6. Finally, the paper ends with a brief conclusion in Section 7.


## 2. Ù‘The Crank–Nicolson Difference Scheme

For the numerical solution of Equations (1)–(3), we introduce a uniform grid of mesh points xi,yj,tn with xi=iΔx, i=0,1,⋯,I, yj=jΔy, j=0,1,⋯,J, and tn=nΔt, n=0,1,⋯,N.Using the Crank–Nicolson approximation to Equations (1)–(3), we have
(5)∂αu∂tαi,jn+1/2+uuxi,jn+1+uuxi,jn2+uuyi,jn+1+uuyi,jn2=νuxxi,jn+1+uxxi,jn2+uyyi,jn+1+uyyi,jn2+fi,jn+1/2.We use the following linearization technique for nonlinear termuuxn+1 and uuyn+1 [27]
(6)uuxn+1≈un+1uxn+unuxn+1−unuxn,uuyn+1≈un+1uyn+unuyn+1−unuyn.Substituting the above approximation into Equation (5), we yield
(7)∂αu∂tαi,jn+1/2+ui,jn+1uxi,jn+ui,jnuxi,jn+12+ui,jn+1uyi,jn+ui,jnuyi,jn+12=νuxxi,jn+1+uxxi,jn2+uyyi,jn+1+uyyi,jn2+fi,jn+1/2.A discrete approximation to the0CFDtαux,y,t at xi,yj,tn+1/2 can be obtained by the following approximation
(8)∂αuxi,yj,tn+1/2∂tα=11−α∫0tn+12∂uxi,yj,s∂se−α1−αtn+12−sds=11−α∫0tn∂uxi,yj,s∂se−α1−αtn+12−sds+∫tntn+12∂uxi,yj,s∂se−α1−αtn+12−sds=11−α∑k=1n∫tk−1tkui,jk−ui,jk−1Δt+s−tk−12uttxi,yj,cke−α1−αtn+12−sds+∫tntn+12ui,jn+1/2−ui,jnΔt/2+OΔte−α1−αtn+12−sds=11−α∑k=1n∫tk−1tkui,jk−ui,jk−1Δt+s−tk−12uttxi,yj,cke−α1−αtn+12−sds+∫tntn+12ui,jn+1+ui,jn/2−ui,jnΔt/2+OΔte−α1−αtn+12−sds,whereck∈tk−1,tk. Then,
(9)∂αuxi,yj,tn+1/2∂tα=11−α∑k=1n∫tk−1tkui,jk−ui,jk−1Δte−α1−αtn+12−sds+11−α∑k=1n∫tk−1tks−tk−12uttxi,yj,cke−α1−αtn+12−sds+11−α∫tntn+12ui,jn+1−ui,jnΔt+OΔte−α1−αtn+12−sds=1αΔtui,jn+1−ui,jn1−e−α1−α12Δt+∑k=1nui,jk−ui,jk−1e−α1−αn−k+12Δt−e−α1−αn−k+32Δt+R1+R2,which(10)R1=11−α∑k=1n∫tk−1tks−tk−12uttxi,yj,cke−α1−αtn+12−sds,and(11)R2=1α1−e−α1−α12ΔtOΔt.Therefore, we have(12)R1=11−α∑k=1n∫tk−1tks−tk−12uttxi,yj,cke−α1−αtn+12−sds=11−α∑k=1nuttxi,yj,ck∫tk−1tks−tk−12e−α1−αtn+12−sds≤max1≤k≤nuttxi,yj,ckΔt21−α∑k=1n∫tk−1tke−α1−αtn+12−sds=max1≤k≤nuttxi,yj,ckΔt2α∑k=1ne−α1−αn−k+12Δt−e−α1−αn−k+32Δt=max1≤k≤nuttxi,yj,ckΔt2αe−α1−α12Δt−e−α1−αn+12Δt≈max1≤k≤nuttxi,yj,ckΔt2αα1−αn+12−12Δt=max1≤k≤nuttxi,yj,ckΔt21−αnΔt=max1≤k≤nuttxi,yj,ckTΔt21−α=OΔt.By setting(13)wk=e−α1−αk−12Δt−e−α1−αk+12Δt,finally, we obtain(14)∂αuxi,yj,tn+1/2∂tα=1αΔt−wnui,j0+∑k=1n−1wn−k+1−wn−kui,jk+w1ui,jn+ui,jn+1−ui,jn1−e−α1−α12Δt+OΔt.Besides, utilizing the Taylor expansion, we have(15)∂2uxi,yj,tn+1/2∂x2=12ui+1,jn+1−2ui,jn+1+ui−1,jn+1Δx2+ui+1,jn−2ui,jn+ui−1,jnΔx2+OΔt2+Δx2,(16)∂2uxi,yj,tn+1/2∂y2=12ui,j+1n+1−2ui,jn+1+ui,j−1n+1Δy2+ui,j+1n−2ui,jn+ui,j−1nΔy2+OΔt2+Δy2,(17)uxi,yj,tn+12∂uxi,yj,tk+1/2∂x=14Δxui,jnui+1,jn+1−ui−1,jn+1+ui,jn+1ui+1,jn−ui−1,jn+OΔt2+Δx2,(18)uxi,yj,tn+12∂uxi,yj,tn+1/2∂y=14Δyui,jnui,j+1n+1−ui,j−1n+1+ui,jn+1ui,j+1n−ui,j−1n+OΔt2+Δy2.Using the Equations (14)–(18), we derive the following difference scheme which is accurate of the order OΔt+Δx2+Δy2,(19)1αΔt−wnUi,j0+∑k=1n−1wn−k+1−wn−kUi,jk+w1Ui,jn+Ui,jn+1−Ui,jnD0=ν2Ui+1,jn+1−2Ui,jn+1+Ui−1,jn+1Δx2+Ui+1,jn−2Ui,jn+Ui−1,jnΔx2+ν2Ui,j+1n+1−2Ui,jn+1+Ui,j−1n+1Δy2+Ui,j+1n−2Ui,jn+Ui,j−1nΔy2−14ΔxUi,jnUi+1,jn+1−Ui−1,jn+1+Ui,jn+1Ui+1,jn−Ui−1,jn−14ΔyUi,jnUi,j+1n+1−Ui,j−1n+1+Ui,jn+1Ui,j+1n−Ui,j−1n+fxi,yj,tn+12,whereD0=1−e−α/1−α1/2Δt and Ui,jn represents the approximate solution of Equation (1).After simplification, we obtain(20)D0+αΔt4ΔxUi+1,jn−Ui−1,jn+αΔt4ΔyUi,j+1n−Ui,j−1n+αΔtνΔx2+αΔtνΔy2Ui,jn+1−αΔtν2Δx2−Ui,jnαΔt4ΔxUi+1,jn+1−αΔtν2Δx2+Ui,jnαΔt4ΔxUi−1,jn+1−αΔtν2Δy2−Ui,jnαΔt4ΔyUi,j+1n+1−αΔtν2Δy2+Ui,jnαΔt4ΔyUi,j−1n+1=D0−αΔtνΔx2−αΔtνΔy2−w1Ui,jn+αΔtν2Δx2Ui+1,jn+αΔtν2Δx2Ui−1,jn+αΔtν2Δy2Ui,j+1n+αΔtν2Δy2Ui,j−1n+wnUi,j0−∑k=1n−1wn−k+1−wn−kUi,jk+αΔtfi,jn+12.


## 3. Fractional Explicit Decoupled Group Method

Another approximate formula for Equation (1) is obtained by Taylor’s expansion and rotating Equation (20), 45° degrees clockwise around the x−y axis. The Crank–Nicolson rotation formula for Equation (1) is as follows
(21)1αΔt−wnUi,j0+∑k=1n−1wn−k+1−wn−kUi,jk+w1Ui,jn+Ui,jn+1−Ui,jnD0=ν4Ui−1,j+1n+1−2Ui,jn+1+Ui+1,j−1n+1h2+Ui−1,j+1n−2Ui,jn+Ui+1,j−1nh2+ν4Ui+1,j+1n+1−2Ui,jn+1+Ui−1,j−1n+1h2+Ui+1,j+1n−2Ui,jn+Ui−1,j−1nh2−14hUi,jn2Ui+1,j+1n+1+Ui+1,j−1n+1−Ui−1,j−1n+1+Ui−1,j+1n+1+Ui,jn+12Ui+1,j+1n+Ui+1,j−1n−Ui−1,j−1n+Ui−1,j+1n−14hUi,jn2Ui−1,j+1n+1+Ui+1,j+1n+1−Ui−1,j−1n+1+Ui+1,j−1n+1+Ui,jn+12Ui−1,j+1n+Ui+1,j+1n−Ui−1,j−1n+Ui+1,j−1n+fxi,yj,tn+12,whereh=Δx=Δy. Similarly, the above-rotated difference scheme is accurate of order OΔt+Δx2+Δy2. On simplification with rx=αΔt/4h and rxx=αΔtν/2h2, the following equation is obtained
(22)D0+rxUi+1,j+1n−Ui−1,j−1n+2rxxUi,jn+1−rxx2Ui+1,j−1n+1−rxx2Ui−1,j+1n+1−rxx2−rxUi,jnUi+1,j+1n+1−rxx2+rxUi,jnUi−1,j−1n+1=D0−2rxx−w1Ui,jn+rxx2Ui+1,j−1n+rxx2Ui−1,j+1n+rxx2Ui+1,j+1n+rxx2Ui−1,j−1n+wnUi,j0−∑k=1n−1wn−k+1−wn−kUi,jk+αΔtfi,jn+1/2.Utilizing Equation (22) to any group of four points on the solution domain gives a 4×4 system as follows
(23)ai,j−ci,j00−di+1,j+1ai+1,j+10000ai+1,j−b00−bai,j+1Ui,jn+1Ui+1,j+1n+1Ui+1,jn+1Ui,j+1n+1=rhsi,jrhsi+1,j+1rhsi+1,jrhsi,j+1,where(24)ai,j=D0+rxUi+1,j+1n−Ui−1,j−1n+2rxx,ci,j=rxx2−rxUi,jn,di,j=rxx2+rxUi,jn,b=rxx2,with(25)rhsi,jrhsi+1,j+1rhsi+1,jrhsi,j+1=bUi+1,j−1n+1+Ui−1,j+1n+1+di,jUi−1,j−1n+1+D0−2rxx−w1Ui,jn+bUi+1,j−1n+Ui−1,j+1n+Ui+1,j+1n+Ui−1,j−1n+τi,jbUi+2,jn+1+Ui,j+2n+1+ci+1,j+1Ui+2,j+2n+1+D0−2rxx−w1Ui+1,j+1n+bUi+2,jn+Ui,j+2n+Ui+2,j+2n+Ui,jn+τi+1,j+1bUi+2,j−1n+1+ci+1,jUi+2,j+1n+1+di+1,jUi,j−1n+1+D0−2rxx−w1Ui+1,jn+bUi+2,j−1n+Ui,j+1n+Ui+2,j+1n+Ui,j−1n+τi+1,jbUi−1,j+2n+1+ci,j+1Ui+1,j+2n+1+di,j+1Ui−1,jn+1+D0−2rxx−w1Ui+1,jn+bUi+1,jn+Ui−1,j+2n+Ui+1,j+2n+Ui−1,jn+τi,j+1,and(26)τi,jτi+1,j+1τi+1,jτi,j+1=wnUi,j0−∑k=1n−1wn−k+1−wn−kUi,jk+αΔtfi,jn+12wnUi+1,j+10−∑k=1n−1wn−k+1−wn−kUi+1,j+1k+αΔtfi+1,j+1n+12wnUi+1,j0−∑k=1n−1wn−k+1−wn−kUi+1,jk+αΔtfi+1,jn+12wnUi,j+10−∑k=1n−1wn−k+1−wn−kUi,j+1k+αΔtfi,j+1n+12.Equation (23) leads to a decoupled system of 2×2 equations in explicit form
(27)ai,j−ci,j−di+1,j+1ai+1,j+1Ui,jn+1Ui+1,j+1n+1=rhsi,jrhsi+1,j+1,and(28)ai+1,j−b−bai,j+1Ui+1,jn+1Ui,j+1n+1=rhsi+1,jrhsi,j+1.The computational molecule of Equations (27) and (28) is shown in Figure 1.Figure 1
Grid point onx−y plane. Computational molecule of Equation (27) (a) and computational molecule of Equation (28) (b).
(a)(b)From Figure1(a), it can be seen that Equation (27) is executed only by considering the green dots. On the contrary, Equation (28) only runs with red dots. Therefore, the implementation of these two equations is independent of each other, which makes the solution of Equation (1) consume less time.In the EDG method, the grid points are divided into several groups. Each group consists of only two points of the grid (shown in Figure2). We apply one of the Equation (27) or Equation (28) for each group in Figure 2. Therefore, half of the grid points (green dots) are calculated by the rotated finite-difference Equation (22). Before going to the next time level, we obtain other points of the grid (red dots) directly once by taking the Equation (20).Figure 2
Group ordering for EDG method forN=5.


## 4. Stability Analysis

In this section, the stability of the finite difference method is investigated with the Von-Neumann analysis. We first give a lemma aboutwj, which will be used in the stability analysis.Lemma 1 (see [28]).
The coefficientswj in Equation (13) satisfy the following properties.
(29)0≤wj≤CΔt,and
(30)0≤wj−wj+1≤CΔtwj.
To investigate the stability of the difference scheme, the nonlinear termuux+uy in Equations (1)–(3) has been linearized by making the quantity u to a local constant. Thus, the nonlinear term in Equation (1) converts into u^ux+uy, and Equation (1) becomes:
(31)∂αu∂tα+u^∂u∂x+∂u∂y=ν∂2u∂x2+∂2u∂y2+fx,y,t.
LetU~i,jn and U¯i,jn be the approximate solutions of Equations (20) and (22), respectively, and define
(32)ρi,jn=Ui,jn−U~i,jn,(33)ϕi,jn=Ui,jn−U¯i,jn,1≤i≤I,1≤j≤J,1≤n≤N.
Then, by substituting Equation (32) into Equation (20), we have
(34)ναΔtΔx2+ναΔtΔy2+D0ρi,jn+1+αΔtu^4Δx−ναΔt2Δx2ρi+1,jn+1+−αΔtu^4Δx−ναΔt2Δx2ρi−1,jn+1+αΔtu^4Δy−ναΔt2Δy2ρi,j+1n+1+−αΔtu^4Δy−ναΔt2Δy2ρi,j−1n+1=−ναΔtΔx2−ναΔtΔy2−w1+D0ρi,jn+−αΔtu^4Δx+ναΔt2Δx2ρi+1,jn+αΔtu^4Δx+ναΔt2Δx2ρi−1,jn+−αΔtu^4Δy+ναΔt2Δy2ρi,j+1n+αΔtu^4Δy+ναΔt2Δy2ρi,j−1n+wnρi,j0−∑k=1n−1wn−k+1−wn−kρi,jk.
Also by putting Equation (33) in Equation (22), we get
(35)D0+2rxxϕi,jn+1−rxx2ϕi+1,j−1n+1−ϕi−1,j+1n+1−rxx2−u^rxϕi+1,j+1n+1−rxx2+u^rxϕi−1,j−1n+1=D0−w1+2rxxϕi,jn+rxx2ϕi+1,j−1n+ϕi−1,j+1n+rxx2−u^rxϕi+1,j+1n+rxx2+u^rxϕi−1,j−1n+wnϕi,j0−∑k=1n−1wn−k+1−wn−kϕi,jk.
The Fourier series forρnx,y and ϕnx,y is
(36)ρnx,y=∑m2=−∞∞∑m1=−∞∞ξnm1,m2eι2πm1x+m2y,ϕnx,y=∑m2=−∞∞∑m1=−∞∞ηnm1,m2eι2πm1x+m2y,
whereι=−1 and the amplication factors ξn and ηn are defined by
(37)ξnm1,m2=∫01∫01e−ι2πm1τ+m2ερnτ,εdτdε,(38)ηnm1,m2=∫01∫01e−ι2πm1τ+m2εϕnτ,εdτdε.
Introducing the following norm(39)ρn2=∑j=1J−1∑i=1I−1ΔxΔyρi,jn21/2=∫01∫01ρi,jn2dτdε1/2,ϕn2=∑j=1J−1∑i=1I−1ΔxΔyϕi,jn21/2=∫01∫01ϕi,jn2dτdε1/2.
By applying the Parseval’s equality(40)∫01∫01ρnτ,ε2dτdε=∑m2=−∞∞∑m1=−∞∞ξnm1,m22,∫01∫01ϕnτ,ε2dτdε=∑m2=−∞∞∑m1=−∞∞ηnm1,m22,
we obtain(41)ρn22=∑m2=−∞∞∑m1=−∞∞ξnm1,m22,ϕn22=∑m2=−∞∞∑m1=−∞∞ηnm1,m22.
According to the above analysis, we can suppose that the solution of Equations (34) and (35) has the following form
(42)ρi,jn=ξneισxiΔx+σyjΔy,ϕi,jn=ηneισxiΔx+σyjΔy,
whereσx=2m1π,σy=2m2π.Lemma 2.
Ifw1≤D0 in Equation (20) and α∈0,1, then we have
(43)ξn≤ξ0,n=1,2,⋯,N,whereξn is defined in Equation (37).Proof.
Substitutingρi,jn=ξneισxiΔx+σyjΔy into Equation (34), we have
(44)ναΔtΔx2+ναΔtΔy2+D0ξn+1eισxiΔx+σyjΔy+αΔtu^4Δx−ναΔt2Δx2ξn+1eισxi+1Δx+σyjΔy+−αΔtu^4Δx−ναΔt2Δx2ξn+1eισxi−1Δx+σyjΔy+αΔtu^4Δy−ναΔt2Δy2ξn+1eισxiΔx+σyj+1Δy+−αΔtu^4Δy−ναΔt2Δy2ξn+1eισxiΔx+σyj−1Δy=−ναΔtΔx2−ναΔtΔy2−w1+D0ξneισxiΔx+σyjΔy+−αΔtu^4Δx+ναΔt2Δx2ξneισxi+1Δx+σyjΔy+αΔtu^4Δx+ναΔt2Δx2ξneισxi−1Δx+σyjΔy+−αΔtu^4Δy+ναΔt2Δy2ξneισxiΔx+σyj+1Δy+αΔtu^4Δy+ναΔt2Δy2ξneισxiΔx+σyj−1Δy+wnξ0eισxiΔx+σyjΔy−∑k=1n−1wn−k+1−wn−kξkeισxiΔx+σyjΔy,
after simplifications, we can get(45)ξn+1−ναΔt−1Δx2−1Δy2+cosσxΔxΔx2+cosσyΔyΔy2+D0+ιαΔtu^2sinσxΔxΔx+sinσyΔyΔy=ξnναΔt−1Δx2−1Δy2+cosσxΔxΔx2+cosσyΔyΔy2+D0−w1−ιαΔtu^2sinσxΔxΔx+sinσyΔyΔy+wnξ0−∑k=1n−1wn−k+1−wn−kξk.
First, lettingn=0 in Equation (45), we obtain
(46)ξ1=γξ0,
where(47)γ=ναΔtA+D0−ιαΔtu^/2B−ναΔtA+D0+ιαΔtu^/2B,A=−1Δx2−1Δy2+cosσxΔxΔx2+cosσyΔyΔy2,B=sinσxΔxΔx+sinσyΔyΔy.
In the above expression, it is clear that the real part of the numerator is smaller than the real part of the denominator. Thus, the magnitude of the numerator is smaller than the denominator. So we have(48)ξ1≤ξ0.
Now, suppose that we have proved thatξn≤ξ0,n=1,2,⋯,m.
We should prove this forn=m+1. Using Equation (45), we get
(49)ξm+1≤ναΔtA+D0−w1−ιαΔtu^/2B−ναΔtA+D0+ιαΔtu^/2Bξ0+wm−∑k=1m−1wm−k+1−wm−k−ναΔtA+D0+ιαΔtu^/2Bξ0.
Using Lemma1, we have
(50)ξm+1≤ναΔtA+D0−w1−ιαΔtu^/2B+w1−ναΔtA+D0+ιαΔtu^/2Bξ0.Consider the following two cases:Case 3.
IfναΔtA+D0−w1>0, then, we have
(51)ξm+1≤ναΔtA+D0−ναΔtA+D0ξ0≤ξ0.Case 4.
IfναΔtA+D0−w1≤0, then, we have
(52)ξm+1≤2w1−ναΔtA−D0−ναΔtA+D0ξ0.Therefore,(53)2w1−ναΔtA−D0−ναΔtA+D0≤1,That is,w1≤D0, or ∣ξm+1∣≤∣ξ0∣.By mathematical induction, we finish the proof.Theorem 5.
Forα∈0,1, the finite difference scheme Equation (20) is stable if w1≤D0.Proof.
Letw1≤D0, using Lemma 2 and Parseval’s equality, we get
(54)ρn2=∑j=1J−1∑i=1I−1ΔyΔxρi,jn2=ΔyΔx∑j=1J−1∑i=1I−1ξneισxiΔx+σyjΔy2=ΔyΔx∑j=1J−1∑i=1I−1ξn2≤ΔyΔx∑j=1J−1∑i=1I−1ξ02=ΔyΔx∑j=1J−1∑i=1I−1ξ0eισxiΔx+σyjΔy2=ρ02.
So the difference scheme Equation (20) is conditionally stable.Lemma 6.
Ifw1≤D0 in Equation (22) and α∈0,1, then, we have
(55)ηn≤η0,n=1,2,⋯,N,where ηn is defined in Equation (38).Proof.
Substitutingϕi,jn=ηneισxiΔx+σyjΔy into Equation (35), we have
(56)D0+2rxxηn+1eισxiΔx+σyjΔy−rxx2ηn+1eισxi+1Δx+σyj−1Δy−rxx2ηn+1eισxi−1Δx+σyj+1Δy−rxx2−u^rxηn+1eισxi+1Δx+σyj+1Δy−rxx2+u^rxηn+1eισxi−1Δx+σyj−1Δy=D0−w1−2rxxηneισxiΔx+σyjΔy+rxx2ηneισxi+1Δx+σyj−1Δy+rxx2ηneισxi−1Δx+σyj+1Δy+rxx2−u^rxηneισxi+1Δx+σyj+1Δy+rxx2+u^rxηneισxi−1Δx+σyj−1Δy+wnη0eισxiΔx+σyjΔy−∑k=1n−1wn−k+1−wn−kηkeισxiΔx+σyjΔy,
by simple computation and noticing thateιβ+e−ιβ=2cosβ and eιβ−e−ιβ=2ιsinβ, we can get
(57)ηn+1D0+rxx2−A+2ιu^rxB=ηnD0−w1−rxx2−A−2ιu^rxB+wnη0+∑k=1n−1wn−k−wn−k+1ηk,
where(58)A=cosσxΔx−σyΔy+cosσxΔx+σyΔy,B=sinσxΔx+σyΔy.
First, settingn=0 in Equation (57), we obtain
(59)η1=D0−rxx2−A−2ιu^rxBD0+rxx2−A+2ιu^rxBη0=D0−rxx2−A2+2u∧rxB2D0+rxx2−A2+2u∧rxB2η0≤η0.
Now, assume that we have proved thatηn≤η0,n=1,2,⋯,m.
We should prove this forn=m+1. Utilizing Equation (57), we obtain
(60)ηm+1≤D0−w1−rxx2−A−2ιu^rxBD0+rxx2−A+2ιu^rxBη0+wm+∑k=1m−1wm−k−wm−k+1D0+rxx2−A+2ιu^rxBη0.
According to Lemma1, we have
(61)ηm+1≤D0−w1−rxx2−A−2ιu^rxB+w1D0+rxx2−A+2ιu^rxBη0.Consider the following two cases:Case 7.
IfD0−w1−rxx2−A>0, then, we have
(62)ηm+1≤D0−rxx2−AD0+rxx2−Aη0≤η0.Case 8.
IfD0−w1−rxx2−A≤0, then, we have
(63)ηm+1≤2w1−D0+rxx2−AD0+rxx2−Aη0.So that,(64)2w1−D0+rxx2−AD0+rxx2−A≤1,This meansw1≤D0, or ∣ηm+1∣≤∣η0∣.By mathematical induction, the proof is complete.Theorem 9.
Forα∈0,1, the finite difference scheme Equation (22) is stable if w1≤D0.Proof.
Supposew1≤D0, from Lemma 6 and Parseval’s equality, we get
(65)ϕn2=∑j=1J−1∑i=1I−1ΔyΔxϕi,jn2=ΔyΔx∑j=1J−1∑i=1I−1ηneισxiΔx+σyjΔy2=ΔyΔx∑j=1J−1∑i=1I−1ηn2≤ΔyΔx∑j=1J−1∑i=1I−1η02=ΔyΔx∑j=1J−1∑i=1I−1η0eισxiΔx+σyjΔy2=ϕ02.
So the difference scheme Equation (22) is conditionally stable.


## 5. Convergence Analysis

We first introduce some notations and lemmas which will be used in the convergence analysis.(66)δx2vi,jn=vi+1,jn−2vi,jn+vi−1,jnΔx2,δy2vi,jn=vi,j+1n−2vi,jn+vi,j−1nΔy2,Δx0vi,jn=vi+1,jn−vi−1,jn2Δx,Δy0vi,jn=vi,j+1n−vi,j−1n2Δy,v,w=ΔxΔy∑j=1J−1∑i=1I−1vijwij,v=ΔxΔy∑j=1J−1∑i=1I−1vij21/2,c0=maxx,y,t∈0,L×0,L×0,Tux,y,t,∂u∂xx,y,t,∂u∂yx,y,t.It is straightforward to show(67)Δx0v,v=0,Δy0v,v=0.Notice that in this section we supposeC stands for a positive constant independent of Δt,Δx,Δy,i,j, and n, which may take different values at different places.Lemma 10.
(Discrete Gronwall’s inequality [29]).
Supposed be a nonnegetive constant, zn and fn are nonnegative sequences. Let
(68)zn≤d+∑0≤k<nfkzk,n≥0,
then(69)zn≤dexp∑0≤j<nfj,n≥0.Theorem 11.
The Crank–Nicolson scheme Equation (20) is convergent and the order of convergence is OΔt+Δx2+Δy2.Proof.
Letei,jn be the error at xi,yj,tn as defined below
(70)ei,jn=uxi,yj,tn−Ui,jn=ui,jn−Ui,jn,1≤j≤J,1≤i≤I,1≤n≤N.
Substituting Equation (70) in Equation (19), we get the following error equations
(71)−wneij0+∑k=1n−1wn−k+1−wn−keijk+w1eijn+eijn+1−eijnD0=ναΔt2δx2ei,jn+1+ei,jn+δy2ei,jn+1+ei,jn−αΔt2σ4+σ5+αΔtRi,jn+12,1≤j≤J−1,1≤i≤I−1,1≤n≤N,e0,jn=eI,jn=0,1≤n≤N,0≤j≤J,ei,0n=ei,Jn=0,1≤n≤N,0≤i≤I,ei,j0=0,0≤j≤J,1≤i≤I,
where(72)Ri,jn+12=OΔt+Δx2+Δy2,σ4=ei,jnΔx0ui,jn+1+ui,jnΔx0ei,jn+1+ei,jn+1Δx0ui,jn+ui,jn+1Δx0ei,jn,σ5=ei,jnΔy0ui,jn+1+ui,jnΔy0ei,jn+1+ei,jn+1Δy0ui,jn+ui,jn+1Δy0ei,jn.
Multiplying Equation (71) by ΔxΔyei,jn+1+ei,jn and summing up for i from 1 to I−1 and j from 1 to J−1, we obtain
(73)∑k=1n−1wn−k+1−wn−kek+w1en+en+1−enD0,en+1+en=ναΔt2δx2en+1+en,en+1+en+δy2en+1+en,en+1+en−αΔt2σ4+σ5,en+1+en+αΔtRn+12,en+1+en,1≤n≤N.
It is clear thatδx2en+1+en,en+1+en≤0 and δy2en+1+en,en+1+en≤0, so we have
(74)D0en+12≤D0−w1en2+w1en+1,en+∑k=1n−1wn−k−wn−k+1ek,en+1+en+αΔt2σ4+σ5,en+1+en+αΔtRn+12,en+1+en,1≤n≤N.
Now, we estimate the third, fourth, and fifth terms of the right-hand side of Equation (74), respectively
(75)∑k=1n−1wn−k−wn−k+1ek,en+1+en=∑k=1n−1wn−k−wn−k+1ek,en+1+∑k=1n−1wn−k−wn−k+1ek,en.
Using Young inequalityab≤εa2+1/4εb2,a,b∈ℝ, and Lemma 1, we have
(76)∑k=1n−1wn−k−wn−k+1ek,en+1+en≤CΔt2∑k=1n−1ek2+CΔten2+en+12.
For the fourth term, utilizing Young inequality and Equation (67), we obtain
(77)σ4,en+1+en≤c01+ε1+ε2+ε3+ε4+ε5en+12+c01+14ε1+14ε2+14ε3+14ε4+14ε5en2.
Similarly, we get(78)σ5,en+1+en≤c01+ε6+ε7+ε8+ε9+ε10en+12+c01+14ε6+14ε7+14ε8+14ε9+14ε10en2.
Also for the fifth term, we have(79)Rn+12,en+1+en≤ε11+14ε12Rn+122+ε12en+12+14ε11en2.
For the second term on the right-hand side of Equation (74), we get
(80)w1en+1,en≤w1ε13en+12+14ε13en2.
Combining Equations (76)–(80) in Equation (74) and using Lemma 1, we obtain
(81)D0∥en+1∥2≤CΔt2∑k=1n−1ek2+Δtλ1en+12+λ2en2+αΔtε11+14ε12Rn+1/22,1≤n≤N,
where(82)λ1=αc022+ε1+ε2+ε3+ε4+ε5+ε6+ε7+ε8+ε9+ε10+αε12+Cε13+C,λ2=αc022+14ε1+14ε2+14ε3+14ε4+14ε5+14ε6+14ε7+14ε8+14ε9+14ε10+C4ε13+α4ε11+2C.
Therefore, we get(83)C−λ1en+12≤CΔt∑k=1n−1ek2+λ2en2+αε11+14ε12Rn+1/22,1≤n≤N.
In the definition ofλ1, we choose epsilons which C−λ1 be positive, so we have
(84)en+12≤CΔt∑k=1n−1ek2+Cen2+αCRn+1/22,1≤n≤N.
Using Lemma10 and nΔt=T, we get
(85)en+12≤CRn+1/22≤OΔt+Δx2+Δy2,1≤n≤N,
and this completes the proof.Theorem 12.
The EDG method Equation (22) is convergent and the order of convergence is OΔt+Δx2+Δy2.Proof.
The proof is similar to Theorem11.


## 6. Numerical Results

In this section, some numerical examples are considered to demonstrate the efficiency and accuracy of the proposed methods.In numerical examples, we suppose thatux,t, UCNx,t, and UEDGx,t denote the exact, Crank–Nicolson, and the EDG solution, respectively. Also in all Tables, the CN is an abbreviation for Crank–Nicolson Method.The results obtained in this study show that the suggested methods have excellent stability, and they have verified the validity and effectiveness of the presented methods. Notably, we perform all of the computations by MATLAB R2019a software on a 64-bit PC with 2.30 GHz processor and 8 GB memory.Example 13.
Consider Equation (1) with the exact solution
(86)ux,y,t=t31−x221−y22,0≤x,y,t≤1.In Tables1 and 2, the maximum of absolute errors and CPU Times for Crank–Nicolson and EDG methods for α=0.1, 0.3, and α=0.7, 0.9 with T=1, Δt=0.01, Re=10, and different values of N are tabulated, respectively. These results confirm the convergent results. In Figure 3(a), the exact, Crank–Nicolson, and EDG solutions for y=0.5, Re=100, N=49, and Δt=0.01 in different values of t are illustrated. Furthermore in Figure 3(b), the absolute error at T=1 for Re=100, N=49, Δt=0.01, and α=0.5 is portrayed. According to the Figures, we can see that our numerical solutions correspond to the exact solutions.Example 14.
In this example, we assume that the exact solution of Equation (1) is
(87)ux,y,t=x2y2sinπxsinπysint,0≤x,y,t≤1.Table 1
Absolute errors of Example13 at T=1 for Re=10, Δt=0.01, α=0.1, and α=0.3.
Nα=0.1α=0.3CNTimeEDGTimeCNTimeEDGTime94.2593e−030.756.1387e−030.743.9001e−030.226.0143e−030.11199.8379e−040.931.6263e−030.699.0356e−041.011.6124e−030.51491.9672e−049.822.6073e−044.691.8465e−0410.612.6040e−044.73999.1297e−05123.159.7575e−0539.878.8260e−05118.029.3869e−0542.01Table 2
Absolute errors of Example13 at T=1 for Re=10, Δt=0.01, α=0.7, and α=0.9.
Nα=0.7α=0.9CNTimeEDGTimeCNTimeEDGTime92.7474e−030.165.6163e−030.121.8409e−030.175.2940e−030.12196.6801e−040.841.5736e−030.604.8985e−040.861.5707e−030.59491.5163e−0410.262.5934e−045.681.4882e−0410.382.5869e−045.35998.2965e−05119.258.6473e−0541.441.0295e−04121.481.0451e−0442.37Figure 3
The exact, Crank–Nicolson, and EDG solutions in different values oft and y=0.5 (a) and the absolute error of Crank–Nicolson method at T=1 (b) of Example 13 for Δt=0.01, Re=100, N=49, and α=0.5.
(a)(b)The values of initial and boundary conditions andf can be achieved using the exact solution. The exact, Crank–Nicolson, and EDG solutions for y=0.5, Re=50, N=49, and Δt=0.01 in different values of t are shown in Figure 4(a) and in Figure 4(b), the absolute error of Crank–Nicolson method at T=1 for Re=50, N=49, Δt=0.01, and α=0.5 is depicted. Based on the Figures, we can see that the numerical solutions are a good approximation of the exact solutions. In Tables 3 and 4, the maximum of absolute errors and CPU Times for Crank–Nicolson and EDG methods for α=0.1, 0.3, and α=0.7, 0.9 with T=1, Δt=0.01, Re=100, and different values of N are tabulated, respectively. The EDG method generates the numerical solution with almost the same accuracy as the Crank–Nicolson method, but uses less time-consuming in comparison to the Crank–Nicolson method.Example 15.
Assume that the exact solution of Equation (1) is as follows:
(88)ux,y,t=t2x21−x3sinπyexpx+y,0≤x,y,t≤1.Figure 4
The exact, Crank–Nicolson and EDG solutions in different values oft and y=0.5 (a) and the absolute error of Crank–Nicolson method at T=1 (b) of Example 14 for Δt=0.01, Re=50, N=49, and α=0.5.
(a)(b)Table 3
Absolute errors of Example14 at T=1 for Re=100, Δt=0.01, α=0.1, and α=0.3.
Nα=0.1α=0.3CNTimeEDGTimeCNTimeEDGTime94.1318e−030.401.0556e−020.263.7717e−030.539.8104e−030.28198.8813e−041.042.6531e−030.738.1156e−041.172.5458e−030.73491.3364e−0412.144.1507e−044.811.2247e−0415.244.0109e−044.86993.4588e−05122.311.0389e−0440.553.1779e−05130.361.0045e−0440.71Table 4
Absolute errors of Example14 at T=1 for Re=100, Δt=0.01, α=0.7, and α=0.9.
Nα=0.7α=0.9CNTimeEDGTimeCNTimeEDGTime92.9459e−030.578.2550e−030.342.4730e−030.967.5148e−030.26196.3611e−041.392.3322e−030.845.4346e−041.562.2536e−030.84499.7981e−0515.193.7198e−045.698.6200e−0517.263.5924e−045.32992.5774e−05126.659.3273e−0544.772.4146e−05136.809.1092e−0543.01The maximum of absolute errors and CPU Times for Crank–Nicolson and EDG methods forα=0.1, 0.3, 0.7, 0.9 with T=1, Δt=0.01, Re=5, and different values of N are shown in Tables 5 and 6, respectively. In Figure 5(a), the exact, Crank–Nicolson, and EDG solutions for y=0.5, Re=20, N=49, and Δt=0.01 in different values of t are illustrated. Also, the absolute error of Crank–Nicolson method at T=1 for Re=20, N=49, Δt=0.01, and α=0.5 is portrayed in Figure 5(b). The numerical experiments verified our theoretical results once again.Example 16.
Consider Equation (1) with the exact solution
(89)ux,y,t=t2cosπxcosπy,0≤x,y,t≤1.Table 5
Absolute errors of Example15 at T=1 for Re=5, Δt=0.01, α=0.1, and α=0.3.
Nα=0.1α=0.3CNTimeEDGTimeCNTimeEDGTime92.1509−030.517.1351e−030.242.0937e−030.196.9914e−030.12194.8726e−041.091.7356e−030.594.7437e−040.911.7078e−030.56497.3239e−0511.882.6457e−044.867.1213e−0511.802.6032e−045.04991.7989e−05121.066.4817e−0540.011.7399e−05120.956.3688e−0540.27Table 6
Absolute errors of Example15 at T=1 for Re=5, Δt=0.01, α=0.7, and α=0.9.
Nα=0.7α=0.9CNTimeEDGTimeCNTimeEDGTime91.9317e−030.216.5894e−030.121.8282e−030.236.3412e−030.13194.3742e−040.971.6292e−030.604.1260e−041.011.5796e−030.61496.4857e−0512.212.4776e−045.295.9238e−0512.582.3834e−045.35991.5095e−05123.575.9878e−0541.381.2086e−05124.235.5847e−0541.61Figure 5
The exact, Crank–Nicolson and EDG solutions in different values oft and y=0.5 (a) and the absolute error of Crank–Nicolson method at T=1 (b) of Example 15 for Δt=0.01, Re=20, N=49, and α=0.5.
(a)(b)This example does not apply to the initial and boundary conditions of the article, but the results of this example are as good as other examples.In Tables7 and 8, the maximum of absolute errors and the CPU time consumed by our proposed methods for α=0.2, 0.4, and 0.6, 0.8 at T=1, Δt=0.01, Re=2 with different mesh sizes are presented, respectively. Similar to other examples, the EDG method is faster than the Crank–Nicolson method. The exact, Crank–Nicolson, and EDG solutions for y=0.5, Re=10, N=49, α=0.5, and Δt=0.01 in different values of t are displayed in Figure 6(a). In addition, the absolute error of the Crank–Nicolson method at T=1, α=0.5, Re=10, N=49, and Δt=0.01 is displayed in Figure 6(b).Table 7
Absolute errors of Example16 at T=1 for Re=2, Δt=0.01, α=0.2, and α=0.4.
Nα=0.2α=0.4CNTimeEDGTimeCNTimeEDGTime92.5900e−030.212.1227e−020.192.5752e−030.202.1177e−020.18195.9920e−040.835.3687e−030.445.9600e−040.845.3597e−030.59499.3436e−059.418.2492e−045.199.2993e−059.708.2383e−045.12992.5225e−05117.522.0397e−0439.942.5145e−05119.212.0376e−0439.86Table 8
Absolute errors of Example16 at T=1 for Re=2, Δt=0.01, α=0.6, and α=0.8.
Nα=0.6α=0.8CNTimeEDGTimeCNTimeEDGTime92.5546e−030.222.1112e−020.202.5253e−030.242.1039e−020.22195.9150e−040.845.3476e−030.605.8469e−040.885.3318e−030.62499.2340e−059.888.2234e−045.059.0870e−0510.198.2010e−045.51992.5001e−05119.142.0348e−0441.282.4238e−05123.962.0280e−0441.67Figure 6
The exact, Crank–Nicolson, and EDG solutions in different values oft (a) and the absolute error of Crank–Nicolson method at T=1 (b) of Example 16 for Δt=0.01, Re=10, N=49, and α=0.5.
(a)(b)Example 17.
In this example, we assume that the exact solution of Equation (1) is
(90)ux,y,t=t2x−x22y−y22,0≤x,y,t≤1.In Figure7(a), the exact, Crank–Nicolson, and EDG solutions for y=0.5, Re=60, N=49, and Δt=0.01 in different values of t are shown. Besides, the absolute error of Crank–Nicolson is demonstrated in Figure 7(b). Obviously, our schemes are very accurate and quickly converge to the exact solution. The maximum of absolute errors and CPU Times for Crank–Nicolson and EDG methods for α=0.2, 0.4, and 0.6, 0.8 at T=1, Δt=0.01, and Re=100 with various mesh sizes are expressed in Tables 9 and 10, respectively.Figure 7
The exact, Crank–Nicolson and EDG solutions in different values oft (a) and the absolute error of Crank–Nicolson method at T=1 (b) of Example 16 for Δt=0.01, Re=60, N=49, and α=0.5.
(a)(b)Table 9
Absolute errors of Example17 at T=1 for Re=100, Δt=0.01, α=0.2, and α=0.4.
Nα=0.2α=0.4CNTimeEDGTimeCNTimeEDGTime92.3544e−050.371.2078e−040.152.0308e−050.321.1028e−040.14195.3089e−060.893.7501e−050.554.5547e−060.903.5462e−050.73497.8572e−079.706.1862e−064.846.5077e−079.725.9056e−065.23991.9520e−07117.731.5238e−0639.651.9513e−07118.381.4377e−0640.22Table 10
Absolute errors of Example17 at T=1 for Re=100, Δt=0.01, α=0.6, and α=0.8.
Nα=0.6α=0.8CNTimeEDGTimeCNTimeEDGTime91.6919e−050.329.8582e−050.261.3331e−050.378.5506e−050.27193.7515e−060.893.3183e−050.822.8606e−060.943.0637e−050.79494.9131e−079.915.5845e−065.412.7513e−0710.065.1849e−065.50991.9501e−07118.221.3268e−0641.161.9462e−071191.1529e−0641.91


## 7. Conclusion

In this paper, we introduced the Crank–Nicolson method and the EDG method which derived from 45° rotation of the Crank–Nicolson approximation point and Taylor expansion to solve the 2D time-fractional Burgers’ equation with Caputo-Fabrizio derivative. The error analysis and local truncation error of these methods gave in detail. The stability of the proposed numerical methods is analyzed by the Von-Neumann method. The convergence analysis of the CN and EDG methods proved. Some test problems chose to investigate the applicability and practical efficiency. From Tables1–10, the results showed a good agreement with the exact solution, and the EDG method was faster than the CN method. Numerical experiments showed the efficiency of the proposed methods in terms of CPU time and accuracy.



---

*Source: 6669287-2021-03-16.xml*