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

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 the 
 
 x
 
 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.


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][3][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 fractionalorder 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][12][13][14][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][17][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: with the initial condition and the boundary condition 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 uðx, y, tÞ defined as: where MðαÞ is a normalization function such that Mð0Þ = Mð1Þ = 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][24][25][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 Section 2, 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.
Using the Crank-Nicolson approximation to Equations (1)-(3), we have We use the following linearization technique for nonlinear term ðuu x Þ n+1 and ðuu y Þ n+1 [27] uu x ð Þ n+1 ≈ u n+1 u n x + u n u n+1 x − u n u n x , uu y À Á n+1 ≈ u n+1 u n y + u n u n+1 y − u n u n y : Substituting the above approximation into Equation (5), we yield A discrete approximation to the CF 0 D α t uðx, y, tÞ at ðx i , y j , t n+1/2 Þ can be obtained by the following approximation 2 Advances in Mathematical Physics where c k ∈ ðt k−1 , t k Þ. Then, which and Therefore, we have By setting finally, we obtain Besides, utilizing the Taylor expansion, we have ∂ 2 u x i , y j , t n+1/2
After simplification, we obtain

Fractional Explicit Decoupled Group Method
Another approximate formula for Equation (1) is obtained by Taylor's expansion and rotating Equation (20), 45°d egrees clockwise around the x − y axis. The Crank-Nicolson rotation formula for Equation (1) is as follows where h = Δx = Δy. Similarly, the above-rotated difference scheme is accurate of order OðΔt + Δx 2 + Δy 2 Þ. On simplification with r x = αΔt/4h and r xx = αΔtν/2h 2 , the following equation is obtained Utilizing Equation (22) to any group of four points on the solution domain gives a ð4 × 4Þ system as follows where Advances in Mathematical Physics with and Equation (23) leads to a decoupled system of 2 × 2 equations in explicit form and The computational molecule of Equations (27) and (28) is shown in Figure 1.
From Figure 1(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 Figure 2). 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).

Stability Analysis
In this section, the stability of the finite difference method is investigated with the Von-Neumann analysis. We first give a lemma about w j , which will be used in the stability analysis.
Lemma 1 (see [28]). The coefficients w j in Equation (13) satisfy the following properties. and To investigate the stability of the difference scheme, the nonlinear term uðu x + u y Þ 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 x + u y Þ, and Equation (1) becomes: LetŨ n i,j and U n i,j be the approximate solutions of Equations (20) and (22), respectively, and define Then, by substituting Equation (32) into Equation (20), we have Also by putting Equation (33) in Equation (22), we get The Fourier series for ρ n ðx, yÞ and ϕ n ðx, yÞ is where ι = ffiffiffiffiffi ffi −1 p and the amplication factors ξ n and η n are defined by Introducing the following norm

Advances in Mathematical Physics
By applying the Parseval's equality we obtain According to the above analysis, we can suppose that the solution of Equations (34) and (35) has the following form where σ x = 2m 1 π, σ y = 2m 2 π. Lemma 2. If w 1 ≤ D 0 in Equation (20) and α ∈ ð0, 1Þ, then we have where ξ n is defined in Equation (37).
Proof. Substituting ρ n i,j = ξ n e ιðσ x iΔx+σ y jΔyÞ into Equation (34), we have after simplifications, we can get First, letting n = 0 in Equation (45), we obtain where 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 Now, suppose that we have proved that jξ n j ≤ jξ 0 j, n = 1, 2, ⋯, m:

Advances in Mathematical Physics
We should prove this for n = m + 1. Using Equation (45), we get Using Lemma 1, we have Consider the following two cases: Therefore, That is, w 1 ≤ D 0 , or |ξ m+1 | ≤ | ξ 0 | . By mathematical induction, we finish the proof.
Proof. Substituting ϕ n i,j = η n e ιðσ x iΔx+σ y jΔyÞ into Equation (35), we have by simple computation and noticing that e ιβ + e −ιβ = 2 cos ðβÞ and e ιβ − e −ιβ = 2ιsinðβÞ, we can get where First, setting n = 0 in Equation (57), we obtain Now, assume that we have proved that jη n j ≤ jη 0 j, n = 1, 2, ⋯, m: 8 Advances in Mathematical Physics We should prove this for n = m + 1. Utilizing Equation (57), we obtain According to Lemma 1, we have Consider the following two cases: So that, This means w 1 ≤ D 0 , or |η m+1 | ≤ | η 0 | . By mathematical induction, the proof is complete.
Proof. Suppose w 1 ≤ D 0 , from Lemma 6 and Parseval's equality, we get So the difference scheme Equation (22) is conditionally stable.

Convergence Analysis
We first introduce some notations and lemmas which will be used in the convergence analysis.
It is straightforward to show Notice that in this section we suppose C stands for a positive constant independent of Δt, Δx, Δy, i, j, and n, which may take different values at different places.
Proof. Let e n i,j be the error at ðx i , y j , t n Þ as defined below Substituting Equation (70) in Equation (19), we get the following error equations 9 Advances in Mathematical Physics where Multiplying Equation (71) by ΔxΔyðe n+1 i,j + e n i,j Þ and summing up for i from 1 to I − 1 and j from 1 to J − 1, we obtain It is clear that jδ 2 x ðe n+1 + e n Þ, e n+1 + e n j ≤ 0 and jδ 2 y ðe n+1 + e n Þ, e n+1 + e n j ≤ 0, so we have Þ e k , e n+1 + e n + αΔt 2 σ 4 + σ 5 , e n+1 + e n + αΔt R n+ 1 2 , e n+1 + e n , 1 ≤ n ≤ N: Now, we estimate the third, fourth, and fifth terms of the right-hand side of Equation (74), respectively Using Young inequality ab ≤ εa 2 + ð1/4εÞb 2 , a, b ∈ ℝ, and Lemma 1, we have For the fourth term, utilizing Young inequality and Equation (67), we obtain

Theorem 12. The EDG method Equation
Proof. The proof is similar to Theorem 11.

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 that uðx, tÞ, U CN ðx, tÞ, and U EDG ðx, 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.

16
Advances in Mathematical Physics  Example 14. In this example, we assume that the exact solution of Equation (1) is The values of initial and boundary conditions and f 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 Example 15. Assume that the exact solution of Equation (1) is as follows: The 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 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.
Example 17. In this example, we assume that the exact solution of Equation (1) is In Figure 7(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

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 Tables 1-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.

Data Availability
All results have been obtained by conducting the numerical procedure and the ideas can be shared for the researchers.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.