Abstract

Numerical solutions for Burgers’ equation based on the Galerkins’ method using cubic B-splines as both weight and interpolation functions are set up. It is shown that this method is capable of solving Burgers’ equation accurately for values of viscosity ranging from very small to large. Three standard problems are used to validate the proposed algorithm. A linear stability analysis shows that a numerical scheme based on a Cranck-Nicolson approximation in time is unconditionally stable.

1. Introduction

A study of Burgers’ equation is important since it arises in the approximate theory of flow through a shock wave propagating in a viscous fluid [1] and in the modeling of turbulence [2]. Burgers’ equation and the Navier-Stokes equation are similar in the form of their nonlinear terms and in the occurrence of higher order derivatives with small coefficients in both [3]. The applications have been found in field as diverse as number theory, gas dynamics, heat conduction, elasticity, and so forth [4]. The exact solutions of the one-dimensional Burgers’ equation have been surveyed by Benton and Platzman [5]. However, difficulties arise in the numerical solution of Burgers’ equation for small values of the viscosity coefficient, that is large Reynolds numbers, which correspond to steep wave fronts [6]. In these cases, numerical methods are likely to produce results which include large nonphysical oscillations unless the size of the elements in unrealistically small. Many studies have been done on the numerical solutions of Burgers’ equation to deal with solutions for the large Reynolds number values. Some of the earlier numerical studies are documented as follows: a finite element method has been given by Caldwell et al. [7] to solve Burgers’ equation by altering the size of the element at each stage using information from three previous steps. Galerkin and Petrov-Galerkin finite element methods involving a time dependent grid [8, 9], and product approximation methods [10], have been used successfully to obtain accurate numerical solutions even for small viscosity coefficients.

Rubin and Graves Jr. [11] have used the spline function technique and quasilinearization for the numerical solution of Burgers’ equation in one space variable. A cubic spline collocation procedures were developed for Burgers’ equation in the papers [1215]. Soliman [16, 17] has solved Burgers’ equation numerically using a new similarity technique, and analytically using the modified extended tanh-function method. Raslan [18] has solved Burgers’ equation using collocation method based on quadratic B-spline finite element method. In this paper we set up the finite element method using Galerkins’ method with cubic splines as shape functions to obtain a numerical solution for Burgers’ equation. The resulting system will be a system of ordinary differential equations which can be solved using a Crank-Nicolson approximation in time.

2. The Governing Equation

Consider Burgers’ equation in the form where and are positive parameters and the subscripts and denote temporal and spatial differentiation, respectively. Appropriate boundary conditions will be chosen from

Applying the Galerkin approach [1921] to (2.1) with weight function , integrating by parts, and choosing the boundary conditions from (2.2) lead to the weak form where the right-hand side of (2.3) is evaluated only at the boundaries. The conditions on the interpolation functions are now simply that only the functions and their first derivative need to be continuous throughout the region. However, we have chosen to use, as trial functions, the very adaptable cubic splines with their well-known advantages.

3. The B-Spline Finite Element Solution

The region is partitioned into finite elements of equal length by knots such that , and let be the cubic B-splines with knots at the points . The set of splines form a basis for functions defined over . We seek the approximation to the solution , which uses these splines as trial functions where the are time-dependent quantities to be determined from (2.3).

Defining a local coordinate system for the finite element with nodes at and , each cubic B-spline covers four elements [22]; consequently each element is covered by four splines which are given in terms of a local coordinate system , where and , by It is the representation of cubic B-spline over a single element which is most appropriate for the finite element approach where all other splines are zero over this element.

The spline and its two principle derivatives vanish outside the interval . In Table 1 we list for convenience the values of and its derivatives , at the relevant knots.

We now identify the finite elements for the problem with the intervals , and the element nodes with the knots .

Using (3.1) and Table 1, we see that the nodal parameters are given in terms of the parameters by and the variation of a function over the element is given by

In addition, we have the valuable property that determine also the first and second derivatives at the nodes (element boundaries), and these are also continuous

The finite element equations we will set up will not be expressed in terms of the nodal parameters but in terms of element parameters , so we shall not directly determine the nodal values as is the case with the usual finite element formulations, but these can always be recovered using (3.3) and (3.5).

We will now set up the element matrices relevant to (2.3). For a typical element , we have the contribution

Using (3.4) in (3.7) and identifying the weight function with the cubic splines we obtain which can be written in matrix form as where .

The element matrices are given by the integrals where take only the values , for this typical element . The matrices are therefore and the matrix is . An associated matrix is given by which also depends on the parameters and will be used in the following theoretical discussions.

The element matrices , and can be determined algebraically from (3.2). Then, we obtain

Combining contributions from all elements leads to the matrix equation where and assembled from the element matrices in the usual way, are septadiagonal in form. We have assumed homogeneous boundary conditions on or and so the term on the right-hand side of (2.3) is zero. The general row for each matrix has the following form: where , for row . The matrices and are symmetric while the form of the matrix has a somewhat more complex structure.

We will time centre on and write where the are the parameters at time , and is the timestep. Then (3.14), can be written as

The matrices , are independent of the time; so, they will remain constant throughout the calculations. While the matrix is dependent on the time, it must therefore be recalculated at each step. Apply the boundary conditions and eliminate from (3.18) which then becomes septadiagonal system and can be solved using septadiagonal algorithm.

The time evolution of the approximate solution is determined by the time evolution of the vector . This is found by repeatedly solving the system (3.18), once the initial vector has been computed from the initial conditions. The system (3.18) is septadiagonal and so can be solved using the septadiagonal algorithm, but an inner iteration is needed at each timestep to scope with the nonlinear terms in which the matrix depends on . The following solution procedure is followed.(i)At time , for the initial step of the inner iteration, we approximate by calculated from only and obtain a first approximation to from (3.18). We then iterate, using (3.18) with matrix calculated from , for up to 10-times to refine the approximation to .(ii)At the other timesteps we use for the matrix , at the first step of inner iteration, derived from , to obtain a first approximation to , by solving (3.18). We then iterate, using (3.18) with matrix calculated from , two- or three-times to refine the approximation to .

4. The Stability Analysis

An investigation into the stability of the numerical scheme (3.18) is based on the Von Neumann theory in which the growth factor of a typical Fourier mode defined as where is the mode number, and is the element size, is determined for a linearization of the numerical scheme (3.18).

In this linearization, we assume that the quantity in the nonlinear term is locally constant. This is equivalent to assuming that the corresponding values of are also constant and equal to .

A linearized recurrence relationship corresponding to (3.18) is then given by where Substituting (4.1) into (4.2) we obtain where

Let and substitute in (4.4) to give where is the growth factor for the mode. Then, the modulus of the growth factor is Therefore the linearized scheme is unconditionally stable.

5. The Initial State

From the initial condition on the function , we must determine the initial vector so that the time evolution of , using (3.18), can be started.

First rewrite (3.1) for the initial condition as where are unknown parameters to be determined. To do this we require to satisfy the following constraints:(a)it must agree with the initial condition at the knots; (3.3) leads to conditions, that is, ,(b)the first or second derivatives of the approximate initial condition shall agree with those of the exact initial condition at both ends of the range; (3.5) or (3.6) produces two further equations, that is, .

The initial vector is then determined as the solution of a matrix equation derived from (3.3)–(3.6) This system can be solved using the Thomas algorithm.

6. The Test Problems

Numerical solutions of Burgers’ equation for three test problems will now be considered. To measure the accuracy of the numerical algorithm, we compute the difference between the analytic and numerical solutions at each mesh point after specified timesteps and use these to calculate the discrete and error norms.(a) Burgers’ equation has the analytic solution [23] where . With (6.1), at time as initial condition, numerical solutions are determined up to a time of , and use as boundary conditions . We discuss the three cases.(i)Figures 1 and 2 shows us the behavior of the numerical solutions for viscosity coefficients and at times from to and from to , respectively. It is seen that for the smaller viscosity value the propagation front is steeper. These graphs agree with those reported by Nguyen and Reynen [23] to within plotting accuracy. In both cases, if the exact solutions are plotted on the same diagram, the curves are indistinguishable. The and error norms are compared with results presented in the papers [9, 14], which are given in Tables 1, 2, and 3. From these tables, we find that the errors obtained by the present algorithm are smaller than the errors calculated by [9, 14].(ii)The case with , , , at different times , and , respectively. The numerical and the analytical solutions computed by the present approach are given in Table 4, with corresponding previous results [15]. From Table 4, we deduce that the numerical solutions computed by the present algorithm are more accurate than those evaluated by [15], at different times.(iii)In this case, we take , , and viscosity coefficients , . The results given in Tables 5 and 6 show us that the computed errors using the present algorithm are smaller than the corresponding errors obtained by Raslan [18], even, when we chose large timestep.(b) A second analytic solution is [11]

We take as initial condition (6.2), at , over the range , with boundary conditions

For , a very weak shock wave develops, when , we obtain a moderate shock wave and when , a strong shock wave is produced. As the value of is decreased the propagation front becomes steeper. The and error norms for these simulations with the corresponding errors obtained by [9], are given in Tables 7 and 8. For large values of the errors are small and as the value of is decreased the errors tend to increase, but for the values of used here, the errors are still acceptable. From Tables 7 and 8, we deduce that the errors computed by the present algorithm are more accurate than those evaluated by [9] at different times. Figure 3 shows us the behavior of the numerical solutions for viscosity coefficients , at different times from to .(c) A third test example, consider the particular analytic solution of Burgers’ equation [8, 10] where , and are constants.

The comparison can be made with [810, 14, 15], these constants are chosen to have the values , , , and . The initial condition is obtained from (6.4), when . The boundary conditions are taken as

From Figure 4, it can be seen that the solution represents a travelling wave, initially situated at , moving to the right with speed . The present method behaves very well until the wave encounters the right-hand boundary. The numerical results at are compared with the analytic solution, and the results presented in papers [810, 14, 15], known as methods of method of lines (ML), standard Galerkin approach (SGA), product approximation Galerkin method (PAG), cubic B-spline collocation method (CSCM), and a numerical solution of Burgers’ equation using cubic B-spline (CSCM), are given in Table 9. It is seen from Table 9 that the agreement between the numerical and the exact solution appears very satisfactorily.

7. Conclusions

We have seen that the algorithm proposed here using Galerkins’ method with cubic spline shape functions compares well in the accuracy with the other published methods, such as methods of lines, standard Galerkin approach, product approximation Galerkin method, and cubic B-spline collocation method. The and error norms are satisfactorily small during the simulations. A linear stability analysis based on the Von Neumann theory shows that the numerical scheme is unconditionally stable. We conclude that a finite element approach based on Galerkins’ method with cubic spline shape functions is eminently suitable for the computation of solutions to Burgers’ equation. We believe that the approach should be suitable for other applications where the continuity of derivatives is essential.