Finite Integration Method with Shifted Chebyshev Polynomials for Solving Time-Fractional Burgers’ Equations

: The Burgers’ equation is one of the nonlinear partial differential equations that has been studied by many researchers, especially, in terms of the fractional derivatives. In this article, the numerical algorithms are invented to obtain the approximate solutions of time-fractional Burgers’ equations both in one and two dimensions as well as time-fractional coupled Burgers’ equations which their fractional derivatives are described in the Caputo sense. These proposed algorithms are constructed by applying the ﬁnite integration method combined with the shifted Chebyshev polynomials to deal the spatial discretizations and further using the forward difference quotient to handle the temporal discretizations. Moreover, numerical examples demonstrate the ability of the proposed method to produce the decent approximate solutions in terms of accuracy. The rate of convergence and computational cost for each example are also presented.


Introduction
Fractional calculus has received much attention due to the fact that several real-world phenomena can be demonstrated successfully by developing mathematical models using fractional calculus. More specifically, fractional differential equations (FDEs) are the generalized form of integer order differential equations. The applications of the FDEs have been emerging in many fields of science and engineering such as diffusion processes [1], thermal conductivity [2], oscillating dynamical systems [3], rheological models [4], quantum models [5], etc. However, one of the interesting issues for the FDEs is a fractional Burgers' equation. It appears in many areas of applied mathematics and can describe various kinds of phenomena such as mathematical models of turbulence and shock wave traveling, formation, and decay of nonplanar shock waves at the velocity fluctuation of sound, physical processes of unidirectional propagation of weakly nonlinear acoustic waves through a gas-filled pipe, and so on, see [6][7][8]. In order to understand these phenomena as well as further apply them in the practical life, it is important to find their solutions. Some powerful numerical methods had been developed for solving the fractional Burgers' equation, such as finite difference methods (FDM) [9], Adomian decomposition method [10], and finite volume method [11]. Moreover, in 2015, Esen and Tasbozan [12] gave a numerical solution of time fractional Burgers' equation by assuming that the solution u(x, t) can be approximated by a linear combination of products of two functions, one of which involves only x and the other involves only t. Recently, Yokus and kaya [13] used the FDM to find the numerical

Preliminaries
Before embarking into the details of the FIM-SCP for solving time-fractional differential equations, we provide in this section the basic definitions of fractional derivatives and shifted Chebyshev polynomials. The necessary notations and some important facts used throughout this paper are also given. More details on basic results of fractional calculus can be found in [23] and further details of Chebyshev polynomials can be reached in [22]. Definition 1. Let p, µ, and t be real numbers such that t > 0, and C µ = {u(t) | u(t) = t p u 1 (t), where u 1 (t) ∈ C[0, ∞) and p > µ} .
If an integrable function u (t) ∈ C µ , we define the Riemann-Liouville fractional integral operator of order where Γ(·) is the well-known Gamma function.

Definition 2.
The Caputo fractional derivative D α of u(t) ∈ C m −1 , with u(t) ∈ C m µ if and only if u (m) ∈ C µ , is defined by where m ∈ N and t > 0.
(iii) Let {x k } n k=1 be a set of zeros of T * n (x) defined in (2), and define the shifted Chebyshev matrix T by . . . . . . . . .

Improved FIM-SCP
In this section, we improve the technique of Boonklurb et al. [21] to construct the first and higher order integration matrices in one and two dimensions. We note here that Boonklurb et al. used Chebyshev polynomials to construct the integration matrices and obtained numerical algorithms for solving linear differential equations, whereas in this work, we use the shifted Chebyshev polynomials to construct first and higher order shifted Chebyshev integration matrices to obtain numerical algorithms that are applicable to solve time-fractional Burgers' equations on any domain [0, L] rather than [−1, 1].

One-Dimensional Shifted Chebyshev Integration Matrices
Let M ∈ N and L ∈ R + . Define an approximate solution u(x) of a certain PDE by the linear combination of shifted Chebyshev polynomials (1), i.e., Let x k , k ∈ {1, 2, 3, ..., M}, be the grid points generated by the zeros of the shifted Chebyshev polynomial T * M (x) defined in (2). Substituting each x k into (3), then (3) can be expressed as and we let it be denoted by u = Tc. The coefficients {c n } M−1 n=0 can be obtained by computing c = T −1 u. Let U (1) (x k ) denote the single layer integration of u from 0 to x k . Then, for k ∈ {1, 2, 3, ..., M} or in matrix form: We denote the above equation by U (1) = Tc = TT −1 u := Au, where A = TT −1 := [a ki ] M×M is called the "shifted Chebyshev integration matrix" for the improved FIM-SCP in one dimension. Next, let us consider the double layer integration of u from 0 to x k that denoted by U (2) (x k ). We have for k ∈ {1, 2, 3, ..., M}, it can be written in matrix form as U (2) = A 2 u. The m th layer integration of u from 0 to x k , denoted by U (m) (x k ), can be obtained in the similar manner, that is, a ki m · · · a i 1 j u(x j ) for k ∈ {1, 2, 3, ..., M}, or written in the matrix form as U (m) = A m u.

Two-Dimensional Shifted Chebyshev Integration Matrices
Let M, N ∈ N and L 1 , L 2 ∈ R + . Divide the domain [0, L 1 ] × [0, L 2 ] into a mesh with M nodes by N nodes along the horizontal and the vertical directions, respectively. Let x k , where k ∈ {1, 2, 3, ..., M}, be the grid points generated by the shifted Chebyshev nodes of T * M (x) and let y s , where s ∈ {1, 2, 3, ..., N}, be the grid points generated by the shifted Chebyshev nodes of T * N (y). Thus, there are M × N grid points in total. For computation, we index the numbering of grid points along the x-direction by the global numbering system (Figure 1a) and along y-direction by the local numbering system (Figure 1b).
Let U (1) x and U (1) y be the single layer integrations with respect to the variables x and y, respectively. For each fixed y, we have U (1) x (x k , y) in the global numbering system as For k ∈ {1, 2, 3, ..., M}, (4) can be expressed as U (1) x (·, y) = A M u(·, y), where A M = TT −1 is the M × M matrix. Thus, for each y ∈ {y 1 , y 2 , y 3 , ..., y N }, is the shifted Chebyshev integration matrix with respect to x-axis and ⊗ is the Kronecker product defined in [24]. Similarly, for each fixed x, U (1) y (x, y s ) can be expressed in the local numbering system as For s ∈ {1, 2, 3, ..., N}, (5) can be written as U (1) . . .

U
(1) We shall denote the above matrix equation by U We notice that the elements of u and u are the same but different positions in the numbering system. Thus, we can transform U (1) y and u in the local numbering system to the global numbering system by using the permutation matrix P = [p ij ] MN×MN , where each p ij is defined by for all k ∈ {1, 2, 3, ..., M} and s ∈ {1, 2, 3, ..., N}. We obtain that U (1) y = P U (1) y and u = P u. Therefore, we have U (1) y = A y u, where A y = P A y P −1 = P(I M ⊗ A N )P is the shifted Chebyshev integration matrix with respect to y-axis in the global numbering system.

Remark 1 ([21]
). Let m, n ∈ N, the multi-layer integrations in the global numbering system can be represented in the matrix forms as follows, • the m th layer integration with respect to x is U the n th layer integration with respect to y is U the multi-layer integration with respect to both x and y is U (m,n) xy = A m x A n y u.

The Numerical Algorithms for Time-Fractional Burgers' Equations
In this section, we derive the numerical algorithms based on our improved FIM-SCP for solving time-fractional Burgers' equations both in one and two dimensions. The numerical algorithm for solving time-fractional coupled Burgers' equations is also proposed. To demonstrate the effectiveness and the efficiency of our algorithms, some numerical examples are given. Moreover, we find the time convergence rates and CPU times(s) of each example in order to demonstrate the computational cost. We note here that we implemented our numerical algorithms in MatLab R2016a. The experimental computer system is configured as: Intel(R) Core(TM) i7-6700 CPU @ 3.40 GHz. Finally, the graphically numerical solutions of each example are also depicted.

Algorithm for One-Dimensional Time-Fractional Burgers' Equation
Let L and T be positive real numbers and α ∈ (0, 1]. Consider the time-fractional Burgers' equation with a viscosity parameter ν > 0 as follows.
subject to the initial condition and the boundary conditions where f (x, t), φ(x), ψ 1 (t), and ψ 2 (t) are given functions. Let us first linearize (7) by determining the iteration at time t m = m(∆t), where ∆t is the time step and m ∈ N. Then, we have where u m = u(x, t m ) is the numerical solution at the m th iteration. For the Caputo time-fractional derivative term defined in Definition 2, we have Using the first-order forward difference quotient to approximate the derivative term in (11), we get where w j = (∆t) −α Γ(2−α) (j + 1) 1−α − j 1−α . Thus, (10) becomes In order to eliminate the derivative terms in (13), we apply the modified FIM by taking the double layer integration. Then, for each shifted Chebyshev node x k , k ∈ {1, 2, 3, ..., M}, we obtain where d 1 and d 2 are the arbitrary constants of integration. Next, we consider the nonlinear term in (14). By using the technique of integration by parts, we have where . Thus, for k ∈ {1, 2, 3, ..., M}, (15) can be expressed in matrix form as For computational convenience, we reduce the above equation into the matrix form: Consequently, for k ∈ {1, 2, 3, ..., M} by hiring (16) and the idea of Boonklurb et al. [21], we can convert (14) into the matrix form as For the boundary conditions (9), we can change them into the vector forms by using the linear combination of the shifted Chebyshev polynomial at the m th iteration as follows.

Example 1. Consider the time-fractional Burgers' Equation
and the boundary conditions The exact solution given by Esen and Tasbozan [12] is u * (x, t) = t 2 e x . In the numerical test, we choose the kinematic viscosity ν = 1, α = 0.5 and ∆t = 0.00025. Table 1 presents the exact solution u * (x, 1), the numerical solution u(x, 1) by using our FIM-SCP in Algorithm 1, and the solution obtained by the quadratic B-spline finite element Galerkin method (QBS-FEM) proposed by Esen and Tasbozan [12]. The comparison between the absolute errors E a (as the difference in absolute value between the approximate solution and the exact solution) of the two methods shows that our FIM-SCP is more accurate than QBS-FEM for M = 10 and similar accuracy for other M. Algorithm 1 acquires the significant improvement in accuracy with less computational nodal points M and regardless the time steps ∆t and the fractional order derivatives α. With the selection of α = 0.5 and M = 40, Table 2 shows the comparison between the exact solution u * (x, 1) and the numerical solution u(x, 1) using Algorithm 1 for various values of ∆t ∈ {0.05, 0.01, 0.005, 0.001}. Table 3 illustrates the comparison between the exact solution u * (x, 1) and the numerical solution u(x, 1) by our method for ∆t = 0.001, M = 40, and α ∈ {0.1, 0.25, 0.75, 0.9}. Moreover, the convergence rates are estimated by using our FIM-SCP with the discretization points M = 20 and step sizes ∆t = 2 −k for k ∈ {4, 5, 6, 7, 8}. In Table 4, we observe that these time convergence rates for the ∞ norm indeed are almost O(∆t) for the different α ∈ (0, 1). Then, we also find the computational cost in term of CPU time(s) in Table 4. Finally, the graph of our approximate solutions, u(x, t), for different times, t, and the surface plot of the solution under the parameters ν = 1, M = 40, and ∆t = 0.001, are provided in Figure 2.

Algorithm 1 The numerical algorithm for solving one-dimensional time-fractional Burgers' equation
Output: An approximate solution u(x, T). 10:

end for
12:

13:
Find u m by solving the iterative linear system (21).
14: end while 15:    x u * (x, 1) and the boundary conditions The exact solution given by Yokus and Kaya [13] In our numerical test, we choose the kinematic viscosity ν = 1, α = 0.8, M = 50 and ∆t = 0.001. Table 5 presents the exact solution u * (x, 0.02), the numerical solution u(x, 0.02) by using our FIM-SCP in Algorithm 1, and the solution obtained by using the expansion method and the Cole-Hopf transformation (EPM-CHT) proposed by Yokus and Kaya in [13]. The error norms L 2 and L ∞ of this problem between our FIM-SCP and EPM-CHT with α = 0.8 for the various values of nodal grid points M ∈ {5, 10, 20, 25, 50} and step size ∆t = 1/M are illustrated in Table 6. We see that our Algorithm 1 achieves improved accuracy with less computational cost. Furthermore, we estimate the convergence rates of time for this problem by using our FIM-SCP with the discretization nodes M = 20 and step sizes ∆t = 2 −k for k ∈ {4, 5, 6, 7, 8} which are tabulated in Table 7. We observe that these rates of convergence for the ∞ norm indeed are almost linear convergence O(∆t) for the different values α ∈ (0, 1). Then, we also calculate the computational cost in term of CPU time(s) as shown in Table 7. Figure 3a,b depict the numerical solutions u(x, t) at different times t and the surface plot of u(x, t), respectively.

Algorithm for Two-Dimensional Time-Fractional Burgers' Equation
Let L 1 and L 2 be positive real numbers, Ω = (0, L 1 ) × (0, L 2 ), and α ∈ (0, 1]. Consider the two-dimensional time-fractional Burgers' equation with a viscosity ν > 0, subject to the initial condition and the boundary conditions where f , φ, ψ 1 , ψ 2 , ψ 3 , and ψ 4 are given functions. As ∂ ∂x ( u 2 2 ) = u ∂u ∂x and ∂ ∂y ( u 2 2 ) = u ∂u ∂y , we can transform (22) to Let us linearize (25) by imposing the iteration at time t m = m(∆t) for m ∈ N and ∆t is an arbitrary time step. Thus, we have where f m = f (x, y, t m ) and u m = u(x, y, t m ) is the numerical solution at the m th iteration. Next, consider the fractional order derivative in the Caputo sense as defined in Definition 2, by using (12), then (26) becomes where w j = (∆t) −α Γ(2−α) (j + 1) 1−α − j 1−α . The above equation can be transformed to the integral equation by taking twice integrations over both x and y, we have where g 1 (y), g 2 (y), h 1 (x), and h 2 (x) are the arbitrary functions emerged in the process of integration which can be approximated by the shifted Chebyshev polynomial interpolation. For r ∈ {1, 2}, define where h (i) r and g (j) r , for i ∈ {0, 1, 2, ..., M − 1} and j ∈ {0, 1, 2, ..., N − 1}, are the unknown values of these interpolated points. Next, we divide the domain Ω into a mesh with M nodes by N nodes along x-and y-directions, respectively. We denote the nodes along the x-direction by x = {x 1 , x 2 , x 3 , ..., x M } and the nodes along the y-direction by y = {y 1 , y 2 , y 3 , ..., y N }. These nodes along the xand y-directions are the zeros of shifted Chebyshev polynomials T * M (x) and T * N (y), respectively. Thus, the total number of grid points in the system is P = M × N, where each point is an entry in the set of Cartesian product x × y ordering as global type system, i.e., (x i , y i ) ∈ x × y for i ∈ {1, 2, 3, ..., P}. By substituting each node in (27) and hiring A x and A y in Section 3.2, we can change (27) to the matrix form as Simplifying the above equation yields where each parameter contained in (29) can be defined as follows.
From (28), we obtain Φ x and Φ y , where For the boundary conditions (24), we can transform them into the matrix form, similar the idea in [21], by employing the linear combination of the shifted Chebyshev polynomials as follows,

•
Left & Right boundary conditions: For each fixed y ∈ {y 1 , y 2 , y 3 , ..., y N }, then where I M and I N are, respectively, the M × M and N × N identity matrices, T −1 M and T −1 N are, respectively, the M × M and N × N matrices defined in Lemma 1, P is defined in (6), and the other parameters are Finally, we can construct the system of iterative linear equations from Equations (29)-(33) for a total of P + 2(M + N) unknowns, including u m , g 1 , g 2 , h 1 and h 2 , as follows, Thus, the approximate solutions u m can be reached by solving (34) in conjunction with the initial condition (23), that is, Therefore, an arbitrary solution u(x, y, t) at any fixed time t can be estimated from Example 3. Consider the 2D time-fractional Burgers' Equation (22) for (x, y) ∈ Ω = (0, 1) × (0, 1) and t ∈ (0, 1] with the forcing term subject to the both homogeneous of initial and boundary conditions. The analytical solution of this problem is u * (x, y, t) = t(x 2 − x)(y 2 − y). For the numerical test, we pick ν = 100, α = 0.5, ∆t = 0.01, and M = N = 10.
In Table 8, the solutions approximated by our FIM-SCP Algorithm 2 are presented in the space domain Ω for various times t. We test the accuracy of our method by measuring it with the absolute error E a . In addition, we seek the rates of convergence via ∞ norm of our Algorithm 2 with the nodal points M = N = 10 and different step sizes ∆t = 2 −k for k ∈ {4, 5, 6, 7, 8}, we found that these convergence rates approach to the linear convergence O(∆t) as shown in Table 9 together with the CPU times(s). Also, the graphically numerical solutions are provided in Figure 4.

Example 4. Consider the 2D Burgers' Equation
subject to the boundary conditions corresponding to the analytical solution given by Cao et al. [14] is u * (x, y, t) = t 3 (1 − x 2 ) 2 (1 − y 2 ) 2 . By picking the parameters ν = 0.1, α = 0.5, and M = N = 10, the comparison of error norm L 2 between our FIM-SCP via Algorithm 2 and the discontinuous Galerkin method combined with finite different scheme (DGM-FDS) presented by Cao et al. [14] are displayed in Table 10 at time t = 0.1. We can see that our method gives a higher accuracy than the DGM-FDS at the same step size ∆t. Next, we provide the CPU times(s) and time convergence rates based on ∞ norm of our algorithm for this problem in Table 11. Then, we see that they converge to the linear rate O(∆t). Finally, the graphical solutions of this Example 4 are provided in Figure 5.

Algorithm for Time-Fractional Coupled Burgers' Equation
Consider the following coupled Burgers' equation with fractional time derivative for α ∈ (0, 1] subject to the initial conditions and the boundary conditions where f (x, t), g(x, t), φ 1 (x) , φ 2 (x), ϕ 1 (t), ϕ 2 (t) , ϕ 3 (t) , and ϕ 4 (t) are the given functions. The procedure of using our FIM for solving u and v are similar, we only discuss here the details in finding the approximate solution u.
We begin with linearizing the system (35) by taking the an iteration of time t m = m(∆t) for m ∈ N, where ∆t is a time step. We obtain where u m = u(x, t m ) and v m = v(x, t m ) are numerical solutions of u and v in the m th iteration, respectively. Next, let us consider the fractional time derivative for α ∈ (0, 1] in the Caputo sense by using the same procedure as in (12), by taking the double layer integration on both sides, we obtain where w γ j = (∆t) −γ Γ(2−γ) (j + 1) 1−γ − j 1−γ for γ ∈ {α, β}, and d 1 , d 2 , d 3 , and d 4 are arbitrary constants of integration. For the nonlinear terms in (38) and (39), by using the same process as in (15), we let For computational convenience, we express q 1 (x k ) and q 2 (x k ) into matrix forms as where T is defined in (17) and other parameters obtained on (40) and (41) are Consequently, using (40), (41), and the procedure in Section 3.1, we can convert both (38) and (39) into the matrix forms as

Rearranging the above system yields
where I is the M × M identity matrix and other parameters are defined by The boundary conditions (37) are transformed into the vector forms by using the same process as in (19) and (20), that is, where t l = [1, −1, 1, ..., (−1) M−1 ] and t r = [1, 1, 1, ..., 1]. Finally, starting from the initial guesses we can construct the system of the m th iterative linear equations for finding numerical solutions. The approximate solutions of u can be obtained from (42) and (44) while the approximate solutions of v can be reached by using (43) and (45): and    For any fixed t, the approximate solutions of u(x, t) and v(x, t) on the space domain can be obtained by computing u( Example 5. Consider the time-fractional coupled Burgers' Equation (35) for x ∈ (0, 1) and t ∈ (0, 1] with the forcing terms subject to the homogeneous initial conditions and the boundary conditions corresponding to the analytical solution given by Albuohimad and Adibi [25] is u * (x, t) = v * (x, t) = xt 3 . For the numerical test, we choose the kinematic viscosity ν = 1, α = β = 0.5 and M = 40. Table 12 presents the exact solution u * (x, 1) and the numerical solutions u(x, 1) together with v(x, 1) by using our FIM-SCP through Algorithm 3. The accuracy is measured by the absolute error E a . Table 13 displays the comparison of the error norms L ∞ of our approximate solutions and the approximate solutions obtained by using the collocation method with FDM (CM-FDM) introduced by Albuohimad and Adibi in [25]. As can be seen from Table 13, our FIM-SCP Algorithm 3 is more accurate. Next, the time convergence rates based on ∞ and CPU times(s) of this problem that solved by Algorithm 3 are demonstrated in Table 14. Since the approximate solutions u and v are the same, we only present the graphical solution of u in Figure 6. Set m = m + 1.

Example 6.
Consider the time-fractional coupled Burgers' Equation (35) for x ∈ (0, 1) and t ∈ (0, 1] with the forcing terms f (x, t) = Γ(4)t −α Γ(4 − α) + 1 t 3 sin(x) and g(x, t) = Γ(4)t −β Γ(4 − β) subject to the homogeneous initial conditions and the boundary conditions corresponding to the analytical solution given by Albuohimad and Adibi [25] is u * (x, t) = v * (x, t) = t 3 sin(x). For the numerical test, we choose the viscosity ν = 1, α = β = 0.5 and M = 5. Table 15 provides the comparison of error norms L ∞ between our FIM-SCP and the CM-FDM in [25] for various values of ∆t and M, it show that our method is more accurate. Moreover, Table 16 illustrates the rates of convergence and CPU times(s) for M = 20. Figure 7a,b show the numerical solutions u(x, t) at different times t and the surface plot of u(x, t), respectively. Note that we only show the graphical solution of u(x, t) since the approximate solutions u(x, t) and v(x, t) are the same.

Conclusions and Discussion
In this paper, we applied our improved FIM-SCP to develop the decent and accurate numerical algorithms for finding the approximate solutions of time-fractional Burgers' equations both in one-and two-dimensional spatial domains and time-fractional coupled Burgers' equations. Their fractional-order derivatives with respect to time were described in the Caputo sense and estimated by forward difference quotient. According to Example 1, even though, we obtain similar accuracy, however, it can be seen that our method does not require the solution to be separable among the spatial and temporal variables. For Example 2, the results confirm that even with nonlinear FDEs, the FIM-SCP provides better accuracy than FDM. For two dimensions, Example 4 shows that even with the small kinematic viscosity ν, our method can deal with a shock wave solution, which is not globally continuously differentiable as that of the classical Burgers' equation under the same effect of small kinematic viscosity ν. We can also see from Examples 5 and 6 that our proposed method can be extended to solve the time-fractional Burgers' equation and we expect that it will also credibly work with other system of time-fractional nonlinear equation. We notice that our method provides better accuracy even when we use a small number of nodal points. Evidently, when we decrease the time step, it furnishes more accurate results. Also, we illustrated the time convergence rate of our method based on ∞ norm, we observe that it approaches to the linear convergence O(∆t). Finally, we show the computational cost in terms of CPU time(s) for each example. An interesting direction for our future work is to extend our technique to solve space-fractional Burgers' equations and other nonlinear FDEs.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript.

CM-FDM
collocation method with finite difference method DGM-FDS discontinuous Galerkin method with finite different scheme EPM-CHT expansion method with Cole-Hopf transformation FDE fractional differential equation FDM finite difference method FIM finite integration method FIM-SCP finite integration method with shifted Chebyshev polynomial PDE partial differential equation QBS-FEM quadratic B-spline finite element Galerkin method