A Galerkin Solution for Burgers ’ Equation Using Cubic B-Spline Finite Elements

and Applied Analysis 3 Table 1 x x −2 x −1 x x 1 x 2 φ 0 1 4 1 0 φ ′ 0 3/h 0 −3/h 0 φ ′′ 0 6/h2 −12/h2 6/h2 0 set of splines {φ−1, φ0, . . . , φN, φN 1} form a basis for functions defined over a, b . We seek the approximation uN x, t to the solution u x, t , which uses these splines as trial functions uN x, t N 1 ∑ −1 δ t φ x , 3.1 where the δ are time-dependent quantities to be determined from 2.3 . Defining a local coordinate system for the finite element x , x 1 with nodes at x and x 1, each cubic B-spline covers four elements 22 ; consequently each element x , x 1 is covered by four splines φ −1, φ , φ 1, φ 2 which are given in terms of a local coordinate system η, where η x − x and 0 ≤ η ≤ h, by ⎪⎨ ⎪⎩ φ −1 φ φ 1 φ 2 1 h3 ⎪⎪⎨ ⎪⎩ ( h − η)3, h3 3h2 ( h − η) 3h(h − η)2 − 3(h − η)3, h3 3h2η 3hη2 − 3η3, η3, 0 ≤ η ≤ h. 3.2 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 φ x and its two principle derivatives vanish outside the interval x −2, x 2 . In Table 1 we list for convenience the values of φ x and its derivatives φ′ x , φ′′ x at the relevant knots. We now identify the finite elements for the problem with the intervals x , x 1 , and the element nodes with the knots x , x 1. Using 3.1 and Table 1, we see that the nodal parameters u are given in terms of the parameters δ by u u x δ −1 4δ δ 1, u 1 u x 1 δ 4δ 1 δ 2, 3.3 and the variation of a function u over the element x , x 1 is given by u 2 ∑ j −1 δjφj . 3.4 4 Abstract and Applied Analysis In addition, we have the valuable property that δ −1 − δ 2 determine also the first and second derivatives at the nodes element boundaries , and these are also continuous u′ 3 h −δ −1 δ 1 , 3.5 u′′ 6 h2 δ −1 − 2δ δ 1 . 3.6 The finite element equations we will set up will not be expressed in terms of the nodal parameters u , u′ 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 x , x 1 , we have the contribution ∫x 1 x v ut εuux dx ν ∫x 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

The Governing Equation
Consider Burgers' equation in the form where ε and ν are positive parameters and the subscripts t and x denote temporal and spatial differentiation, respectively.Appropriate boundary conditions will be chosen from u a, t α, u b, t β, u x a, t u x b, t 0, u xx a, t u xx b, t 0.

2.2
Applying the Galerkin approach 19-21 to 2.1 with weight function v, 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.

The B-Spline Finite Element Solution
The region a, b is partitioned into N finite elements of equal length h by knots x i such that a x 0 < x 1 • • • < x N b, and let φ i x be the cubic B-splines with knots at the points x i .The set of splines {φ −1 , φ 0 , . . ., φ N , φ N 1 } form a basis for functions defined over a, b .We seek the approximation u N x, t to the solution u x, t , 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 x , x 1 with nodes at x and x 1 , each cubic B-spline covers four elements 22 ; consequently each element x , x 1 is covered by four splines φ −1 , φ , φ 1 , φ 2 which are given in terms of a local coordinate system η, where η x − x and 0 ≤ η ≤ h, by

3.2
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 φ x and its two principle derivatives vanish outside the interval x −2 , x 2 .In Table 1 we list for convenience the values of φ x and its derivatives φ x , φ x at the relevant knots.
We now identify the finite elements for the problem with the intervals x , x 1 , and the element nodes with the knots x , x 1 .
Using 3.1 and Table 1, we see that the nodal parameters u are given in terms of the parameters δ by and the variation of a function u over the element x , x 1 is given by u 2 j −1 δ j φ j .

3.4
In addition, we have the valuable property that δ −1 − δ 2 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 u , u 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 x , x 1 , we have the contribution which also depends on the parameters δ e k and will be used in the following theoretical discussions.
The element matrices A e , C e , L e , and B e can be determined algebraically from 3.2 .Then, we obtain where the δ n are the parameters at time nΔt, and Δt is the timestep.Then 3.14 , can be written as

3.18
The matrices A, C are independent of the time; so, they will remain constant throughout the calculations.While the matrix B is dependent on the time, it must therefore be recalculated at each step.Apply the boundary conditions and eliminate δ −1 , δ N 1 from 3.18 which then becomes N 1 × N 1 septadiagonal system and can be solved using septadiagonal algorithm.
The time evolution of the approximate solution u N x, t is determined by the time evolution of the vector δ n .This is found by repeatedly solving the system 3.18 , once the initial vector δ 0 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 B depends on δ 1/2 δ n δ n 1 .The following solution procedure is followed.
i At time t 0, for the initial step of the inner iteration, we approximate B by B * calculated from δ 0 only and obtain a first approximation to δ 1 from 3.18 .We then iterate, using 3.18 with matrix B calculated from δ 1/2 δ 0 δ 1 , for up to 10-times to refine the approximation to δ 1 .
ii At the other timesteps we use for the matrix B, at the first step of inner iteration, B * derived from δ * δ n 1/2 δ n − δ n−1 , to obtain a first approximation to δ n 1 , by solving 3.18 .We then iterate, using 3.18 with matrix B calculated from δ 1/2 δ n δ n 1 , two-or three-times to refine the approximation to δ n 1 .

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 k is the mode number, and h is the element size, is determined for a linearization of the numerical scheme 3.18 .
In this linearization, we assume that the quantity u in the nonlinear term uu x is locally constant.This is equivalent to assuming that the corresponding values of δ j are also constant and equal to d.
A linearized recurrence relationship corresponding to 3.18 is then given by where where g is the growth factor for the mode.Then, the modulus of the growth factor is Therefore the linearized scheme is unconditionally stable.

The Initial State
From the initial condition u x, 0 on the function u x, t , we must determine the initial vector δ 0 so that the time evolution of δ, using 3.18 , can be started.
First rewrite 3.1 for the initial condition as where δ 0 j are unknown parameters to be determined.To do this we require u N x, 0 to satisfy the following constraints: a it must agree with the initial condition u x, 0 at the knots; 3.3 leads to N 1 conditions, that is, u N x , 0 u x , 0 , 0, . . ., N, 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, u x 0 , 0 u x N , 0 0.
The initial vector δ 0 is then determined as the solution of a matrix equation derived from 3.3 -3.6 Mδ 0 b.

5.2
This system can be solved using the Thomas algorithm.

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 L 2 -and L ∞ -error norms.
a Burgers' equation has the analytic solution 23 where t 0 exp 1/8ν .With 6.1 , at time t 1 as initial condition, numerical solutions are determined up to a time of t 4, and use as boundary conditions u 0, t u 1, t 0. We discuss the three cases.
i Figures 1 and 2 shows us the behavior of the numerical solutions for viscosity coefficients ν 0.005 and 0.0005 at times from t 1 to 3.1 and from t 1 to 3.25, 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 L 2 -and L ∞ -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 ν 0.0005, h 0.005, Δt 0.01, at different times t 1.7, 2.5, and 3.25, 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.
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 t 0, over the range 0 ≤ x ≤ 30, with boundary conditions u 0, t 1, u 30, t 0, t > 0. 6.3For ν 1/2, a very weak shock wave develops, when ν 1/8, we obtain a moderate shock wave and when ν 1/24, a strong shock wave is produced.As the value of ν is decreased the propagation front becomes steeper.The L 2 -and L ∞ -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 ν 1/24, at different times from t 0 to 10.
3.7 and identifying the weight function v η with the cubic splines we obtain 2 written in matrix form as A e dδ e dt εδ e T L e δ e νC e δ e 0j φ k dη, 3.10 where i, j, k take only the values − 1, , 1, 2, for this typical element x , x 1 .The matrices A e , C e are therefore 4 × 4 and the matrix L e is 4 × 4 × 4.An associated 4 × 4 matrix B e is given by

Figure 3 : 3 Figure 4 :
Figure 3: The behavior of the numerical solution with different viscosity at different times.

Table 4 :
Comparison of results at different times for ν 0.0005 and 0 ≤ x ≤ 1, with h 0.005 and Δt 0.01, problem a .

Table 5 :
Comparison of errors at different times for ν 0.005 and 0 ≤ x ≤ 1, with h 0.01, problem a .

Table 6 :
Comparison of errors at different times for ν 0.01 and 0 ≤ x ≤ 1, with h 0.01, problem a .

Table 9 :
Comparison of results for various numerical schemes with exact at time t 0.5, for ν 0.01.