L
 1
 -Multiscale Galerkin’s Scheme with Multilevel Augmentation Algorithm for Solving Time Fractional Burgers’ Equation

<jats:p>In this paper, we consider the initial boundary value problem of the time fractional Burgers equation. A fully discrete scheme is proposed for the time fractional nonlinear Burgers equation with time discretized by <jats:inline-formula>
                     <math xmlns="http://www.w3.org/1998/Math/MathML" id="M2">
                        <mi>L</mi>
                        <mn>1</mn>
                     </math>
                  </jats:inline-formula>-type formula and space discretized by the multiscale Galerkin method. The optimal convergence orders reach <jats:inline-formula>
                     <math xmlns="http://www.w3.org/1998/Math/MathML" id="M3">
                        <mi mathvariant="script">O</mi>
                        <mfenced open="(" close=")">
                           <mrow>
                              <msup>
                                 <mrow>
                                    <mi>τ</mi>
                                 </mrow>
                                 <mrow>
                                    <mn>2</mn>
                                    <mo>−</mo>
                                    <mi>α</mi>
                                 </mrow>
                              </msup>
                              <mo>+</mo>
                              <msup>
                                 <mrow>
                                    <mi>h</mi>
                                 </mrow>
                                 <mrow>
                                    <mi>r</mi>
                                 </mrow>
                              </msup>
                           </mrow>
                        </mfenced>
                     </math>
                  </jats:inline-formula> in the <jats:inline-formula>
                     <math xmlns="http://www.w3.org/1998/Math/MathML" id="M4">
                        <msup>
                           <mrow>
                              <mi>L</mi>
                           </mrow>
                           <mrow>
                              <mn>2</mn>
                           </mrow>
                        </msup>
                     </math>
                  </jats:inline-formula> norm and <jats:inline-formula>
                     <math xmlns="http://www.w3.org/1998/Math/MathML" id="M5">
                        <mi mathvariant="script">O</mi>
                        <mfenced open="(" close=")">
                           <mrow>
                              <msup>
                                 <mrow>
                                    <mi>τ</mi>
                                 </mrow>
                                 <mrow>
                                    <mn>2</mn>
                                    <mo>−</mo>
                                    <mi>α</mi>
                                 </mrow>
                              </msup>
                              <mo>+</mo>
                              <msup>
                                 <mrow>
                                    <mi>h</mi>
                                 </mrow>
                                 <mrow>
                                    <mi>r</mi>
                                    <mo>−</mo>
                                    <mn>1</mn>
                                 </mrow>
                              </msup>
                           </mrow>
                        </mfenced>
                     </math>
                  </jats:inline-formula> in the <jats:inline-formula>
                     <math xmlns="http://www.w3.org/1998/Math/MathML" id="M6">
                        <msup>
                           <mrow>
                              <mi>H</mi>
                           </mrow>
                           <mrow>
                              <mn>1</mn>
                           </mrow>
                        </msup>
                     </math>
                  </jats:inline-formula> norm, respectively, in which <jats:inline-formula>
                     <math xmlns="http://www.w3.org/1998/Math/MathML" id="M7">
                        <mi>τ</mi>
                     </math>
                  </jats:inline-formula> is the time step size, <jats:inline-formula>
                     <math xmlns="http://www.w3.org/1998/Math/MathML" id="M8">
                        <mi>h</mi>
                     </math>
                  </jats:inline-formula> is the space step size, and <jats:inline-formula>
                     <math xmlns="http://www.w3.org/1998/Math/MathML" id="M9">
                        <mi>r</mi>
                     </math>
                  </jats:inline-formula> is the order of piecewise polynomial space. Then, a fast multilevel augmentation method (MAM) is developed for solving the nonlinear algebraic equations resulting from the fully discrete scheme at each time step. We show that the MAM preserves the optimal convergence orders, and the computational cost is greatly reduced. Numerical experiments are presented to verify the theoretical analysis, and comparisons between MAM and Newton’s method show the efficiency of our algorithm.</jats:p>


Introduction
In this paper, we consider the following time fractional Burgers equation [1][2][3][4][5][6][7]: with the initial and boundary conditions, given by where 0 < α < 1,Ω = fðx, tÞ | 0 ≤ x ≤ 1, 0 < t ≤ Tg,u 0 ðxÞ and f ðx, tÞ are given functions, and the notation c 0 D α t denotes the Caputo fractional partial derivative of order α, defined by in which Γð·Þ represents the Gamma fuction. The time fractional Burgers equation is a kind of nonlinear subdiffusion convection equation occurring in several physical problems such as unidirectional propagation of weakly nonlinear acoustic waves through a gas-filled pipe, propagation of weak shock, compressible turbulence shallow-water waves, shock waves in a viscous medium, waves in bubbly liquids, and electromagnetic waves [1,4,8]. Till now, there have been several analytical techniques developed to solve the time fractional Burgers equation. These methods include the Cole-Hopf transformation, Laplace transform, variable separation method [8], Adomian decomposition method [4], homotopy analysis method [6], and so on.
However, even if some fractional differential equations can be solved, the expressions of their exact solutions are often expressed by special functions, which are difficult to apply in practice. Moreover, due to the nonlocality of the fractional operators, analytical methods do not always work well on most fractional differential equations in real applications. Hence, it is of great importance to develop reliable and efficient numerical methods for solving fractional differential equations. Nowadays, the numerical methods cover the quadratic B-spline Galerkin method [2], cubic B-spline finite element method [3], finite difference methods [1,[9][10][11][12][13], and Fourier pseudospectral schemes [14]. There are also some other numerical methods (see, for example, [8,[15][16][17]).
In this paper, we first present a fully discrete scheme for solving the time fractional Burgers equation with the time approximated by the L1-type formula and the space discretization based on the multiscale Galerkin method. We give rigorous convergence analysis for the fully discrete scheme, which shows that the scheme enjoys the optimal convergence order Oðτ 2−α + h r Þ in the L 2 norm and Oðτ 2−α + h r−1 Þ in the H 1 norm, respectively, where τ,h, and r are the time step size, space step size, and the order of piecewise polynomial space, respectively. Since the time fractional Burgers equation is a nonlinear differential equation, the fully discrete scheme results in a system of nonlinear algebraic equation at each time step. Iteration methods such as the Newton iteration method and the quasi-Newton iteration method are often employed to solve these nonlinear equations. In this case, a large amount of computational effort is demanded to compute and update the Jacobian matrix in each iteration process. The higher accuracy of the approximate solution is required, the larger dimension of the subspace is needed, and the longer computational time is consumed. To overcome this problem, we develop the multilevel augmentation method for solving the fully discrete scheme. The MAM solves a nonlinear equation at a high level consisting of two parts: solving the nonlinear equation only in a fixed initial subspace with the dimension much lower than that of the whole approximate subspace; compensating the error by matrix-vector multiplications at the high level. The MAM reduces the computational costs significantly and leads to a fast solution for the fully discrete scheme. We prove that the MAM preserves the same optimal convergence order as the original fully discrete scheme. The idea of MAM was first introduced in [18] for solving the linear Fredholm integral equations of the second kind. The theoretical setting of MAM was established by Chen et al. in [19] for solving operator equations covering both first kind and second kind equations; they further develop MAM for solving the nonlinear Hammerstein integral equation in [20]. We modified the framework and extended the idea of MAM to solve general nonlinear operator equations of the second kind and applied it to the Sine-Gordon equation in [21]. Readers are referred to [22][23][24][25][26][27] and the references therein for more applications of MAM.
This paper is organized in seven sections. In "Preliminaries," some necessary notations, multiscale orthonormal bases in Sobolev space, and useful lemmas are introduced. In "L1 Scheme for Discretization of Caputo Derivative in Time," we introduce the L1-formula for time discretization. In "Fully Discrete Scheme and Convergence," a fully discrete scheme for time fractional Burgers equation is established, and the convergence analysis are given. The MAM and its convergence analysis are developed in "Multilevel Augmentation Method for Solving the Fully Discrete Scheme." The numerical experiments are provided in "Numerical Experiments" to verify the theoretical estimates. Finally, a conclusion is included in "Conclusion."

Preliminaries
Denote I = ½0, 1. Let ð·, · Þ stand for the inner product on the space L 2 ðIÞ with the L 2 norm k·k 2 . We denote by H 1 0 ðIÞ the Sobolev space of elements u satisfying the homogeneous boundary conditions that uð0Þ = uð1Þ = 0. The inner product and norm of H 1 0 ðIÞ are defined by respectively. Let n be a positive integer, we denote by X n the subspace of H 1 0 ðIÞ whose elements are the piecewise polynomials of order r with knots j/2 n , j − 1 ∈ ℤ 2 n −1 , where the notation ℤ n ≔ f0, 1, 2 ⋯ , n − 1g: Obviously, the sequence of X n is nested, that is which yields the following decomposition: where W n is the orthogonal complement of X n−1 in X n : It is easily concluded from the definition of X n and W n that the dimensions of X n and W n are given by xðnÞ ≔ dim ðX n Þ = ðr − 1Þ2 n − 1 and wðnÞ ≔ dim ðW n Þ = xðnÞ − xðn − 1Þ = ðr − 1Þ2 n−1 ,respectively.
Define two affine mappings on the interval I by φ 0 ðxÞ = x/2 and φ 1 ðxÞ = x + 1/2, x ∈ I, which map the interval ½0, 1 into ½0, 1/2 and ½1/2, 1, respectively. Associated with the two mappings, we introduce two linear operators as follows: Lemma 1 (see [26]). Let w ij ,j ∈ ℤ wðiÞ , be an orthonormal basis of W i ,i ≥ 1: Then the functions fT 0 w ij , T 1 w ij : j ∈ ℤ wðiÞ g form an orthonormal basis for W i+1 : Lemma 1 shows that the space W n can be recursively constructed by the linear operators T 0 and T 1 once W 1 has been given. Therefore, the basis of the space X n can be 2 Journal of Function Spaces constructed by Lemma 1 step by step. For the details of the construction and more, the readers can refer to [26]. Let P n be an orthogonal projection operator from H 1 0 ðIÞ into X n with respect to the inner product h·, · i, that is, for all u ∈ H 1 0 ðIÞ, or The following approximation results on the operator P n will be used later. Throughout this paper, unless stated otherwise, c denotes a generic positive constant whose value may differ in different occurrences.

Fully Discrete Scheme and Convergence
In this section, we present a fully discrete scheme for the time fractional Burgers equation (1), and we derive the error estimates and convergence of the proposed fully discrete scheme. The Galerkin method associated with the multiscale basis introduced in "Preliminaries" is employed to discretize the spatial variable. The fully discrete scheme in weak formulation Let k, m be two fixed positive integers and n ≔ k + m.
Step 1: obtain the approximation of initial value function u 0 n ≔ P n u 0 ðxÞ = ∑ ði,jÞ∈J n c ij w ij ðxÞ, where c ij = hu 0 ðxÞ, w ij i 1 ,ði, jÞ ∈ J n : Step 2: for i = 1 : N (T = Nτ), do the following: 3 Journal of Function Spaces for (1) reads as follows: for each t = t i ,i = 1, 2 ⋯ n, find u i n ∈ X n such that where f i = f ðx, tÞj t=t i and Π n denotes the interpolation operator.
We present an optimal error estimate of the fully discrete scheme (16) in the following theorem.
Theorem 5. Suppose that the problem (1)-(2) has a unique solution u ∈ C 2 ð½0, T ; L 2 ðΩÞÞ ∩ C 1 ð½0, T ; H r ðΩÞÞ: Then where u i = uðx, tÞj t=t i : Proof. Denote e i n ≔ P n u i − u i n ,i = 0, 1, 2 ⋯ N: We conclude from (1) and (16) that e i n satisfies Taking v n = e i n in (18), we have We estimate the terms of the right-hand side of (19) one by one. For the first term in right-hand side of (19), using the Cauchy-Schwarz inequality and Lemma 2, we have To estimate the second term in the right-hand side of (19), we conclude from Lemma 4 that For the last term in the right-hand side of (19), using integration by parts and the Cauchy-Schwartz inequality, we have where ku i k ≤ M and ku i n k ≤ M, due to the smoothness of u and the approximation u n ∈ X n ⊂ H 1 0 ðIÞ: On the other hand Substituting (20)- (23) into (19) and noting that 1 = a 0 > a 1 > ⋯>a i >⋯, we deduce that Choose τ such that τ −α /2Γð2 − αÞ > ððM/2Þ + 1Þ, and denote σ = τ −α /2Γð2 − αÞ,λ = τ −α /2Γð2 − αÞ − ððM/2Þ + 1Þ; then, we have By Gronwall's inequality, we have which, together with Lemma 2 and the initial error estimate, yields that This completes the proof.
Remark 6. If we choose v n = D α τ e i n in (19) and make a similar analysis as the above Theorem 5, we can obtain the optimal convergence order in H 1 norm

Multilevel Augmentation Method for Solving the Fully Discrete Scheme
At each time step, the fully discrete scheme (16) leads to a nonlinear system, which makes the computational cost expensive. We present a fast multilevel augmentation method in this section to solve these nonlinear systems. To this end, we rewrite (16) into where r i n = ∑ i−1 k=1 ða i−k−1 − a i−k Þu k n + a i−1 u 0 n and μ = τ −α /Γð2 − αÞ: Define a nonlinear operator K :X ⟶ X as follows: Similar to the proof of Lemma 3 in [23], we applied the Riesz representation theorem to the right-hand side of (29); there exists a elementf i ∈ X, such that Then, Equation (29) can be reformulated as or equivalently Since Equation (16) has been reformulated as a nonlinear operator equation of the second kind (33), and K has the properties (P1) and (P2) described in [23], then MAM developed in [21] is applicable.

Journal of Function Spaces
We now briefly describe the MAM for solving (33). As we presented in "Preliminaries," the approximation subspace sequence is nested, for a fixed positive integer k,n ≔ k + m, m is any nonnegative integer, and we have the following decomposition: Now, we are in a position to solve (33) with n ≔ k + m, and k is fixed and smaller than n: Firstly, we solve (33) with n ≔ k exactly and obtain u i k : Next, we obtain an approximation of u i k+1 of (33) with n ≔ k + 1: To this end, we decompose : With the help of (34), Equation (33) with n ≔ k + 1 can be rewritten as an equivalent form as Note that Equation (35) becomes The u i k+1 in the right-hand side can be approximated by the previous level solution u i k,0 ≔ u i k : We compute Replace u i,H k+1 in (36) by u i,H k,1 , and solve u i,L k,1 ∈ X k from Let which is an approximation to the solution u i k+1 : This procedure is repeated m times to obtain an approximation u i k,m of the solution u i k+m of (33) with n = k + m. The solution u i k,m is called a multilevel augmentation solution. Since at any step l = 0, 1, ⋯m, we only need to invert the same nonlinear operator P k ðI − KÞ with a fixed small k instead of the nonlinear operator P k+l ðI − KÞ: This means the algorithm has a high computational efficiency. At every time step, the fully discrete scheme (33) is solved by the MAM, and the whole process can be summarized as the following algorithm: Theorem 7. Let u be the exact solution of (1) and u i k,m be the approximation solution obtained by Algorithm 1. Suppose that the solution of Equation (33) belongs to H r ðIÞ for i = 1, 2, ⋯T/τ: Then, there exist a positive integer N such that for all k ≥ N and m ∈ ℕ Proof. As stated in [20,21,23], u i k,m is the solution of the equation: where v ∈ X k+m ,v 2 = ðP k+m − P k Þv: Rearranging the terms, we have Noting that the exact solution u at t = t i satisfies Subtracting (45) from (46), we obtain that for all v ∈ X k+m ,v 2 = ðP k+m − P k Þv Denote ρ i ≔ u i − P k+m u i and e i ≔ P k+m u i − u i k,m , then u i − u i k,m = ρ i + e i : Using these notations and noting that ðρ i x , v x Þ = 0, we derive the error equation as follows: Let M 1 = max fku i k,m k ∞ , i ∈ ℤ N g and M ′ = max fM 1 , Mg, where M is the positive constant appearing in (22); then, ku i k ∞ ≤ M′,ku i k,m k ∞ ≤ M′: We take v = e i in (48) and estimate the terms in the right-hand side of (48).
For the first three terms in the right-hand side of (48), similar to the analysis of (20)- (22), we have By the Cauchy-Schwartz inequality, Young's inequality, and noting that kv 2 k 2 ≤ kvk 2 , we have For the last term in the right-hand side of (48), it follows from integration by parts, the Cauchy-Schwartz inequality, and the Young inequality that On the other hand side, as presented in (23), we have Combining (49)-(54) and ðe i x , e i x Þ = je i j 2 1 , we have Choose τ such that τ −α /2Γð2 − αÞ > 1 + M ′ + ðσ 1 /4Þ, and denote σ = τ −α /2Γð2 − αÞ, μ ≔ σ − ð1 + M ′ + ðσ 1 /4ÞÞ: Then, we have When the exact solution of (33) belongs to H r ðIÞ, there exists a positive integer N, for all k ≥ N and any m ∈ ℕ (see [21,23]): Combining (56) and (57), we conclude from Gronwall's inequality that Noting that ke 0 k 2 ≤ ch r , then This completes the proof. (48) and make a similar analysis as the above Theorem 7, we can obtain the optimal convergence order in H 1 norm:

Numerical Experiments
We present in this section numerical examples to illustrate the efficiency and accuracy of our proposed method. The computer programs are run on a personal computer with 2.5G CPU and 8G memory.
Example 1. We consider the time fractional Burgers equation (1) with the exact solution: The corresponding initial condition and forcing term are Both piecewise linear (r = 2) and quadratic (r = 3) multiscale orthonormal bases introduced in "Preliminaries" are employed in our numerical approximation. The numerical results are reported in Tables 1-5. k and m stand for the numbers of initial level and augmentation level used in the MAM, respectively. xðnÞ denotes the dimension of approximation subspace X n with n = k + m: Table 1 shows the L 2 errors and temporal convergence rates for different α using the MAM with ðk, mÞ = ð3, 7Þ for linear basis and ðk, mÞ = ð2, 5Þ for quadratic basis. It is seen that our numerical scheme has an accuracy of 2 − α, which is in agreement with our theoretical analysis. In the spatial direction, we illustrate the accuracy, convergence order, and computational efficiency of the MAM, with a comparison to those of the direct Newton's method (DNM) for solving the fully discrete scheme (16). The numerical results listed in Tables 2 and  3 are linear basis cases, and Tables 4 and 5 are quadratic basis cases. We can easily see from these tables that both MAM and DNM have the optimal convergence orders in the H 1 norm (1 for the linear case and 2 for the quadratic  9 Journal of Function Spaces case) and in the L 2 norm (2 for the linear case and 3 for the quadratic case). We also observe that MAM and DNM have nearly the same accuracy, while MAM takes significantly less time than DNM. To intuitively show the approximation effect, we plot in Figure 1 the absolute error surface of the approximation solution u 2,5 obtained by MAM.
Example 2. We consider the time fractional Burgers equation (1) with initial condition: The exact solution of this problem is u x, t ð Þ= t 2 sin 2πx ð Þ: ð64Þ The numerical results are presented in Figure 2 and Tables 6-9, where Figure 2 displays the convergence orders in temporal direction with different α, and Tables 6-9 show the accuracy, convergence order, and computing time for the spatial direction. All the numerical results verify our theoretical analysis and also show the efficiency of the proposed algorithm.

Conclusion
In this article, the L1-discretization formula and the multiscale Galerkin method are adopted to discrete the Caputo fractional derivative and spatial variable, respectively, and the multilevel augmentation algorithm is proposed for solving the resulting fully discrete scheme which is a nonlinear system at each time step. The MAM only needs to solve nonlinear systems in a fixed subspace with much lower dimension than that for the whole approximation subspace and compensate the error by multiplications of matrices and vectors at the high level. Therefore, the computational cost is greatly reduced. Numerical experiments are presented to confirm our theoretical results. Compared with the DNM, the proposed MAM has substantial advantages in computing time and is suitable for solving large-scale and high-accuracy problems.

Data Availability
The data of numerical simulation used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.