Two Different Methods for Numerical Solution of the Modified Burgers' Equation

A numerical solution of the modified Burgers' equation (MBE) is obtained by using quartic B-spline subdomain finite element method (SFEM) over which the nonlinear term is locally linearized and using quartic B-spline differential quadrature (QBDQM) method. The accuracy and efficiency of the methods are discussed by computing L 2 and L ∞ error norms. Comparisons are made with those of some earlier papers. The obtained numerical results show that the methods are effective numerical schemes to solve the MBE. A linear stability analysis, based on the von Neumann scheme, shows the SFEM is unconditionally stable. A rate of convergence analysis is also given for the DQM.


Introduction
The one-dimensional Burgers' equation first suggested by Bateman [1] and later treated by Burgers' [2] has the form where V is a positive parameter and the subscripts and denote space and time derivatives, respectively. Burgers' model of turbulence is very important in fluid dynamics model and study of this model and the theory of shock waves has been considered by many authors for both conceptual understanding of a class of physical flows and for testing various numerical methods [3]. Relationship between (1) and both turbulence theory and shock wave theory was presented by Cole [4]. He also obtained an exact solution of the equation. Analytical solutions of the equation were found for restricted values of V which represent the kinematics viscosity of the fluid motion. So the numerical solution of Burgers' equation has been subject of many papers. Various numerical methods have been studied based on finite difference [5,6], Runge-Kutta-Chebyshev method [7,8], group-theoretic methods [9], and finite element methods including Galerkin, Petrov-Galerkin, least squares, and collocation [10][11][12][13].
The modified Burgers' equation (MBE) which we discuss in this paper is based upon Burgers' equation (BE) of the form The equation has the strong nonlinear aspects and has been used in many practical transport problems, for instance, nonlinear waves in a medium with low-frequency pumping or absorption, turbulence transport, wave processes in thermoelastic medium, transport and dispersion of pollutants in rivers and sediment transport, and ion reflection at quasiperpendicular shocks. Recently, some numerical studies of the equation have been presented: Ramadan and El-Danaf [14] obtained the numerical solutions of the MBE using quintic B-spline collocation finite element method. A special lattice Boltzmann model is developed by Duan et al. [15]. Daǧ et al. [16] have developed a Galerkin finite element solution of the equation using quintic B-splines and time-split technique.
A solution based on sextic B-spline collocation method is proposed by Irk [17]. Roshan and Bhamra [18] applied a Petrov-Galerkin method using a linear hat function as the trial function and a cubic B-spline function as the test function. A discontinuous Galerkin method is presented by Zhang et al. [19]. Bratsos [20] has used a finite difference scheme 2 The Scientific World Journal based on fourth-order rational approximants to the matrixexponential term in a two-time level recurrence relation for calculating the numerical solution of the equation.
Recently, DQM has become a very efficient and effective method to obtain the numerical solutions of various types of partial differential equations. In 1972, Bellman et al. [21] first introduced differential quadrature method (DQM) for solving partial differential equations. The main idea behind the method is to find out the weighting coefficients of the functional values at nodal points by using base functions of which derivatives are already known at the same nodal points over the entire region. Various researchers have developed different types of DQMs by utilizing various test functions; Bellman et al. [22] have used Legendre polynomials and spline functions in order to get weighting coefficients. Quan and Chang [23,24] have presented an explicit formulation for determining the weighting coefficients using Lagrange interpolation polynomials. Zhong [25], Guo and Zhong [26], and Zhong and Lan [27] have introduced another efficient DQM as spline based DQM and applied it to different problems. Shu and Wu [28] have considered some of the implicit formulations of weighting coefficients with the help of radial basis functions. Nonlinear Burgers' equation is solved using polynomial based differential quadrature method by Korkmaz and Daǧ [29]. The DQM has many advantages over the classical techniques; mainly, it prevents linearization and perturbation in order to find better solutions of given nonlinear equations. Since QBDQM do not need transforming for solving the equation, the method has been preferred.
In the present work, we have applied a subdomain finite element method and a quartic B-spline differential quadrature method to the MBE. To show the performance and accuracy of the methods and make comparisons of numerical solutions, we have taken different values of V.

Subdomain Finite Element Method (SFEM).
We will consider (2) with the boundary conditions chosen from with the initial condition where 1 and 2 are constants. The quartic B-splines ( ) ( = −2(1) + 1) at the knots which form a basis over the interval [ , ] are defined by the relationships [30] ( ) Our numerical treatment for solving the MBE using the subdomain finite element method with quartic B-splines is to find a global approximation ( , ) to the exact solution ( , ) that can be expressed in the following form: where are time-dependent parameters to be determined from both boundary and weighted residual conditions. The nodal values , , , and at the knots can be obtained from (5) and (6) in the following form: For each element, using the local coordinate transformation where −2 , −1 , , +1 , and +2 act as element parameters and B-splines −2 ( ), −1 , , +1 , and +2 as element shape functions. Applying the subdomain approach to (33) with the weight function we obtain the weak form of (2) Using the transformation (8) into the weak form (12) and then integrating (12) term by term with some manipulation by parts result in where the dot denotes differentiation with respect to , and In (13) using the Crank-Nicolson formula and its time derivative that is discretized by the forward difference approach, respectively, we obtain a recurrence relationship between the two time levels and + 1 relating two unknown parameters +1 and , for = − 2, − 1, . . . , + 2, Obviously, the system (16) consists of equations in the +4 unknowns ( −2 , −1 , . . . , +1 ). To get a unique solution of the system, we need four additional constraints. These are obtained from the boundary conditions (3) and can be used to eliminate −2 , −1 , , and +1 from the system (16) which then becomes a matrix equation for the unknowns = ( 0 , 1 , . . . , −1 ) of the form A lumped value of is obtained from ( + +1 ) 2 /4 as The resulting system can be solved with a variant of Thomas algorithm and we need an inner iteration ( * ) +1 = + (1/2)( +1 − ) at each time step to cope with the nonlinear term . A typical member of the matrix system (16) is rewritten in terms of the nodal parameters as 4 The Scientific World Journal Before the solution process begins iteratively, the initial vector 0 = ( 0 , 1 , . . . , −1 ) must be determined by means of the following requirements: If we eliminate the parameters 0 −2 , 0 −1 , 0 , and 0 +1 from the system (16), we obtain × matrix system of the following form: where is
We can see that 2 1 < 2 2 and taking the modulus of (38) gives | | ≤ 1, so we find that the scheme (20) is unconditionally stable. (QBDQM). DQM can be defined as an approximation to a derivative of a given function by using the linear summation of its values at specific discrete nodal points over the solution domain of a problem. Provided that any given function ( ) is enough smooth over the solution domain, its derivatives with respect to at a nodal point can be approximated by a linear summation of all the functional values in the solution domain, namely,

Quartic B-Spline Differential Quadrature Method
where denotes the order of the derivative, ( ) represent the weighting coefficients of the th order derivative approximation, and denotes the number of nodal points in the solution domain. Here, the index represents the fact that ( ) is the corresponding weighting coefficient of the functional value ( ). We need first-and second-order derivative of the function ( ). So, we will find value of (28) for the = 1, 2. If we consider (28), then it is seen that the fundamental process for approximating the derivatives of any given function through DQM is to find out the corresponding weighting coefficients ( ) . The main idea of the DQM approximation is to find out the corresponding weighting coefficients ( ) by means of a set of base functions spanning the problem domain. While determining the corresponding weighting coefficients different basis may be used. Using the quartic B-splines as test functions in the fundamental DQM equation (28) leads to the equation

First-Order Derivative Approximation.
When DQM methodology is applied, the fundamental equality for determining the corresponding weighting coefficients of the firstorder derivative approximation is obtained as Korkmaz used In this process, the initial step for finding out the corresponding weighting coefficients (1) , , = −2, −1, . . . , + 3, of the first nodal point 1 is to apply the test functions , = −1, 0, . . . , + 1, at the nodal point 1 . After all the test The Scientific World Journal 5 functions are applied, we get the following system of algebraic equation system: ] ] .
The weighting coefficients (1) 1, related to the first grid point are determined by solving the system (31). This system consists of + 6 unknowns and + 3 equations. To have a unique solution, it is required to add three additional equations to the system. These are By using these equations which we obtained by derivations, three unknown terms will be eliminated from the system. Consider To determine the weighting coefficients, (1) , , = −1, 0, . . . , + 1, at grid points , 2 ≤ ≤ − 1, we got the following algebraic equation system: ] , +1 (1) , +2 . . .
For the last grid point of the domain , determine weighting coefficients, (1) , , = −1, 0, . . . , + 1, we got the following algebraic equation system: ]   where (2) , represents the corresponding weighting coefficients of the second-order derivative approximations. Similarly, for finding out the weighting coefficients of the first grid point 1 all test functions , = −1, 0, . . . , + 1, are used and the following algebraic equations system is obtained: ] .
Here, the system (37) consists of + 6 unknowns and + 3 equations. To have a unique solution, it is required to add new equations to the system. These are The Scientific World Journal If we used (38) and (39) we obtain the following equations system: ] ] .

Test Problem and Experimental Results
In this section, we obtained numerical solutions of the MBE by the subdomain finite element method and differential The Scientific World Journal 9  quadrature method. The accuracy of the numerical method is checked using the error norms 2 and ∞ , respectively, All numerical calculations were computed in double precision arithmetic on a Pentium 4 PC using a Fortran compiler. The analytical solution of modified Burgers' equation is given in [33] as where 0 is a constant and 0 < 0 < 1. For our numerical calculation, we take 0 = 0.5. We use the initial condition for  Tables 1, 2, 3, and 4. It is clearly seen that the results obtained by the SFEM are found to be more acceptable. Also, we observe from these tables that if the value of viscosity decreases, the value of the error norms will decrease. When the value of viscosity parameter is smaller we get better results. The behaviors of the numerical solutions for viscosity V = 0.01, 0.005, 0.001, space steps ℎ = 0.02, 0.005, and time steps Δ = 0.01, 0.001 at times = 1, 2, 4, 6, and 8 are shown in Figures 1, 2, and 3. As seen in the figures, the top curve is at = 1 and the bottom curve is at = 8. Also, we have observed from the figures that as the time increases the curve of the numerical solution decays. With smaller viscosity value, numerical solution decay gets faster.

Experimental Results for QBDQM.
We calculate the numerical rates of convergence (ROC) with the help of the following formula: 10 The Scientific World Journal   Here ( ) denotes either the 2 error norm or the ∞ error norm in case of using grid points. Therefore, some further numerical runs for different numbers of space steps have been performed. Ultimately, some computations have been made about the ROC by assuming that these methods are algebraically convergent in space. Namely, we suppose that ( ) ∼ for some < 0 when ( ) denotes the 2 or the ∞ error norms by using subintervals. is seen from Figure 4 when we select the solution domain 0 ≤ ≤ 1 the right part of the shock wave cannot be seen clearly. By using the larger domain like 0 ≤ ≤ 1.3 as seen in Figure 5 solution is got better than narrow domain 0 ≤ ≤ 1 shown in Figure 4. The computed values of the error norms 2 and ∞ are presented at some selected times up to = 10. The obtained results are recorded in Tables 5 and 6. As it is seen from the tables, the error norms 2 and ∞ are sufficiently small and satisfactorily acceptable. We obtain better results if the value of viscosity parameter is smaller, as shown in Table 7. The behaviors of the numerical solutions for viscosity V = 0.01 and 0.001 and time step Δ = 0.01 at times = 1, 3, 5, 7, and 9 are shown in Figures  4-6. It is observed from the figures that the top curve is at = 1 and the bottom curve is at = 9. It is obviously seen that smaller viscosity value V in shock wave with smaller amplitude and propagation front becomes smoother. Also, we have seen from the figures that, as the time increases, the curve of the numerical solution decays. With smaller viscosity value, numerical solution decay gets faster. These numerical solutions graphs also agree with published earlier work [13]. Distributions of the error at time = 10 are drawn for solitary waves, Figures 7 and 8, from which maximum error happens at the right hand boundary when greater value of viscosity V = 0.01 is used, and with smaller value of viscosity V = 0.001, maximum error is recorded around the location where shock wave has the highest amplitude. The 2 and ∞ error norms and numerical rate of convergence analysis for V = 0.001 and Δ = 0.01 and different numbers of grid points are tabulated in Table 8. Table 9 presents a comparison of the values of the error norms obtained by the present methods with those obtained by other methods [13,14,17,18]. It is clearly seen from the table that the error norm 2 obtained by the SFEM is smaller than those given in [13,14,17,18] whereas the error norm ∞ is very close to those given in [14,17,18]. The error norm ∞ is better than the paper [13]. For the QBDQM both 2 and ∞ are almost the same as these papers.

Conclusion
In this paper, SFEM and DQM based on quartic B-splines have been set up to find the numerical solution of the MBE (2). The performance of the schemes has been considered by studying the propagation of a single solitary wave. The efficiency and accuracy of the methods were shown by calculating the error norms 2 and ∞ . Stability analysis of the approximation obtained by the schemes shows that the [ 13,14,17,18] whereas ∞ error norm is found to be very close to values given in [13,14,17,18]. The obtained results show that our methods can be used to produce reasonable accurate numerical solutions of modified Burgers' equation. So these methods are reliable for getting the numerical solutions of the physically important nonlinear problems.

Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.