Abstract
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 and 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.
1. Introduction
The one-dimensional Burgers’ equation first suggested by Bateman [1] and later treated by Burgers’ [2] has the form where 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 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–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 quasi-perpendicular 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 based on fourth-order rational approximants to the matrix-exponential 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 .
2. Numerical Methods
To implement the numerical schemes, the interval is splitted up into uniformly sized intervals by the nodes , , such that , where .
2.1. Subdomain Finite Element Method (SFEM)
We will consider (2) with the boundary conditions chosen from with the initial condition where and are constants. The quartic B-splines () 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 a typical finite interval is mapped into the interval . Therefore, the quartic B-spline shape functions over the element can be defined as
All other splines, apart from , , , , and , are zero over the element . So, the approximation equation (6) over this element can be written in terms of basis functions given in (9) as where , , , , and act as element parameters and B-splines , , , , and 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 relating two unknown parameters and , for , , where
Obviously, the system (16) consists of equations in the unknowns . 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 , , , and from the system (16) which then becomes a matrix equation for the unknowns of the form A lumped value of is obtained from as The resulting system can be solved with a variant of Thomas algorithm and we need an inner iteration 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 where
Before the solution process begins iteratively, the initial vector must be determined by means of the following requirements:
If we eliminate the parameters , , , and from the system (16), we obtain matrix system of the following form: where is , and . This system is solved by using a variant of Thomas algorithm.
2.2. Linear Stability Analysis
We have investigated stability of the scheme by using the von Neumann method. In order to apply the stability analysis, the MBE needs to be linearized by assuming that the quantity in the nonlinear term is locally constant. The growth factor of a typical Fourier mode is defined as where is mode number and is the element size. Substituting (37) into the scheme (20), we have where We can see that and taking the modulus of (38) gives , so we find that the scheme (20) is unconditionally stable.
2.3. Quartic B-Spline Differential Quadrature Method (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, 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 . 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
2.4. First-Order Derivative Approximation
When DQM methodology is applied, the fundamental equality for determining the corresponding weighting coefficients of the first-order derivative approximation is obtained as Korkmaz used [31] In this process, the initial step for finding out the corresponding weighting coefficients , , of the first nodal point is to apply the test functions , , at the nodal point . After all the test functions are applied, we get the following system of algebraic equation system:
The weighting coefficients related to the first grid point are determined by solving the system (31). This system consists of unknowns and 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, , , at grid points , , we got the following algebraic equation system:
For the last grid point of the domain , determine weighting coefficients, , , we got the following algebraic equation system:
2.5. Second-Order Derivative Approximation
The general form of DQM approximation to the problem on the solution domain is where represents the corresponding weighting coefficients of the second-order derivative approximations. Similarly, for finding out the weighting coefficients of the first grid point all test functions , , are used and the following algebraic equations system is obtained:
Here, the system (37) consists of unknowns and equations. To have a unique solution, it is required to add new equations to the system. These are If we used (38) and (39) we obtain the following equations system: Quartic B-splines have not got fourth-order derivations at the grid points so we cannot eliminate the unknown term by the one more derivation of (39). We will use first-order weighting coefficients instead of second-order weighting coefficients which are introduced by Shu [32] After we use (41), system (42) is obtained. To determine the weighting coefficients , , at grid points , , we got the following algebraic system: where equals .
For the last grid point of the domain with the same idea, determine weighting coefficients , , we got the following algebraic equation system: where equals .
3. Test Problem and Experimental Results
In this section, we obtained numerical solutions of the MBE by the subdomain finite element method and differential quadrature method. The accuracy of the numerical method is checked using the error norms 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 is a constant and . For our numerical calculation, we take . We use the initial condition for (46), evaluating at , and the boundary conditions are taken as and .
3.1. Experimental Results for FEM
For the numerical simulation, we have chosen the various viscosity parameters , space steps , and time steps over the interval . The computed values of the error norms and are presented at some selected times up to . The obtained results are tabulated in 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 , space steps , and time steps at times , and are shown in Figures 1, 2, and 3. As seen in the figures, the top curve is at and the bottom curve is at . 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.
3.2. Experimental Results for QBDQM
We calculate the numerical rates of convergence (ROC) with the help of the following formula:
Here denotes either the 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 when denotes the or the error norms by using subintervals.
For the numerical treatment, we have taken the different viscosity parameters and time step over the intervals and . As it is seen from Figure 4 when we select the solution domain the right part of the shock wave cannot be seen clearly. By using the larger domain like as seen in Figure 5 solution is got better than narrow domain shown in Figure 4. The computed values of the error norms and are presented at some selected times up to . The obtained results are recorded in Tables 5 and 6. As it is seen from the tables, the error norms 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 and and time step at times , and are shown in Figures 4–6. It is observed from the figures that the top curve is at and the bottom curve is at . It is obviously seen that smaller viscosity value 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 are drawn for solitary waves, Figures 7 and 8, from which maximum error happens at the right hand boundary when greater value of viscosity is used, and with smaller value of viscosity , maximum error is recorded around the location where shock wave has the highest amplitude. The and error norms and numerical rate of convergence analysis for and 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 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 and are almost the same as these papers.
4. 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 and . Stability analysis of the approximation obtained by the schemes shows that the methods are unconditionally stable. An obvious conclusion can be drawn from the numerical results that for the SFEM error norm is found to be better than the methods cited in [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.