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

Authors: A. A. Soliman
Journal: Abstract and Applied Analysis
 (2012)

---

## 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.


---

## Body

## 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 [12–15]. 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(2.1)ut+εuux-νuxx=0,a≤x≤b,
where ε and ν are positive parameters and the subscripts t and x denote temporal and spatial differentiation, respectively. Appropriate boundary conditions will be chosen from
(2.2)u(a,t)=α,u(b,t)=β,ux(a,t)=ux(b,t)=0,uxx(a,t)=uxx(b,t)=0.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
(2.3)∫ab‍v(ut+εuux)dx+ν∫ab‍vxuxdx=ν[vux]ab,
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[a,b] is partitioned into N finite elements of equal length h by knots xi such that a=x0<x1⋯<xN=b, and let ϕi(x) be the cubic B-splines with knots at the points xi. The 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
(3.1)uN(x,t)=∑ℓ=-1N+1‍δℓ(t)ϕℓ(x),
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){ϕℓ-1ϕℓϕℓ+1ϕℓ+2=1h3{(h-η)3,h3+3h2(h-η)+3h(h-η)2-3(h-η)3,h3+3h2η+3hη2-3η3,η3,0≤η≤h.
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.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
/
h
2
-
12
/
h
2
6
/
h
2
0We 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
(3.3)uℓ=u(xℓ)=δℓ-1+4δℓ+δℓ+1,uℓ+1=u(xℓ+1)=δℓ+4δℓ+1+δℓ+2,
and the variation of a function u over the element [xℓ,xℓ+1] is given by
(3.4)u=∑j=ℓ-1ℓ+2‍δjϕj.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
(3.5)u′=3h(-δℓ-1+δℓ+1),(3.6)u′′=6h2(δℓ-1-2δℓ+δℓ+1).The finite element equations we will set up will not be expressed in terms of the nodal parametersuℓ,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
(3.7)∫xℓxℓ+1‍v(ut+εuux)dx+ν∫xℓxℓ+1‍vxuxdx=0.Using (3.4) in (3.7) and identifying the weight function v(η) with the cubic splines we obtain
(3.8)∑j=ℓ-1ℓ+2‍[∫0hϕiϕjdη]dδjedt+ε∑j=ℓ-1ℓ+2‍∑k=ℓ-1ℓ+2‍[∫0hϕiϕj′ϕkdη]δkeδje+ν∑j=ℓ-1ℓ+2‍[∫0hϕi′ϕj′dη]δje=0,
which can be written in matrix form as
(3.9)Aedδedt+εδeTLeδe+νCeδe=0,
where δe=(δℓ-1e,δℓe,δℓ+1e,δℓ+2e)T.The element matrices are given by the integrals(3.10)Aije=∫0h‍ϕiϕjdη,Cije=∫0h‍ϕi′ϕj′dη,Lijke=∫0hϕiϕj′ϕkdη,
where i,j,k take only the values ℓ-1,ℓ,ℓ+1,ℓ+2, for this typical element [xℓ,xℓ+1]. The matrices Ae,Ce are therefore 4×4 and the matrix Le is 4×4×4. An associated 4×4 matrix Be is given by
(3.11)Bije=∑k=ℓ-1ℓ+2‍Lijkeδke,
which also depends on the parameters δke and will be used in the following theoretical discussions.The element matricesAe,Ce,Le, and Be can be determined algebraically from (3.2). Then, we obtain
(3.12)Ae=h140(2012960112911889336060933118812916012920),Ce=110h(1821-36-321102-87-36-36-8710221-3-362118),(3.13)B11e=-1840(280,1605,630,5)δe,B12e=-1840(150,1305,792,21)δe,B13e=1840(420,2781,1314,21)δe,B14e=1840(10,129,108,5)δe,B21e=-1840(1605,10830,5349,108)δe,B22e=-1840(1305,17640,17541,1314)δe,B23e=1840(2781,25002,17541,792)δe,B24e=1840(129,3468,5349,630)δe,B31e=-1840(630,5349,3468,129)δe,B32e=-1840(792,17541,25002,2781)δe,B33e=1840(1314,17541,17640,1305)δe,B34e=1840(108,5349,10830,1605)δe,B41e=-1840(5,108,129,10)δe,B42e=-1840(21,1314,2781,420)δe,B43e=1840(21,792,1305,150)δe,B44e=1840(5,630,1605,280)δe.Combining contributions from all elements leads to the matrix equation(3.14)Adδdt+Bδ+νCδ=0,
where
(3.15)δ=(δ-1,δ0,δ1,…,δN+1)T,
and A,B,C assembled from the element matrices Ae,Be,Ce in the usual way, are septadiagonal in form. We have assumed homogeneous boundary conditions on u or u′ and so the term on the right-hand side of (2.3) is zero. The general row for each matrix has the following form:
(3.16)A:h140(1,120,1191,2416,1191,120,1),C:-110h(3,72,45,-240,45,72,3),B:1840(-(5,108,129,10,0,0,0)δ,-(21,1944,8130,3888,129,0,0)δ,-(-21,0,17841,35682,8130,108,0)δ,(5,1944,17841,0,-17841,-1944,-5)δ,(0,108,8130,35682,17841,0,-21)δ,(0,0,129,3888,8130,1944,21)δ,(0,0,0,10,129,108,5)δ),
where δ=(δℓ-3,δℓ-2,δℓ-1,δℓ,δℓ+1,δℓ+2,δℓ+3)T, for row ℓ. The matrices A and C are symmetric while the form of the matrix B has a somewhat more complex structure.We will time centre ont=(n+1/2)Δt and write
(3.17)δ=12(δn+δn+1),dδdt=1Δt(δn+1-δn),
where the δn are the parameters at time nΔt, and Δt is the timestep. Then (3.14), can be written as
(3.18)[A+εΔt2B+νΔt2C]δn+1=[A-εΔt2B-νΔt2C]δn.The matricesA, 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 solutionuN(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 timet=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 matrixB, 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.


## 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
(4.1)δjn=δ^neijkh,
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 quantityu in the nonlinear term uux 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
(4.2)α1δℓ-3n+1+α2δℓ-2n+1+α3δℓ-1n+1+α4δℓn+1+α5δℓ+1n+1+α6δℓ+2n+1+α7δℓ+3n+1=α7δℓ-3n+α6δℓ-2n+α5δℓ-1n+α8δℓn+α3δℓ+1n+α2δℓ+2n+α1δℓ+3n,
where
(4.3)α1=r1-3r2-3r3,α2=120r1-168r2-72r3,α3=1191r1-735r2-45r3,α4=2416r1+240r3,α5=1191r1+735r2+45r3,α6=120r1+168r2+72r3,α7=r1+3r2+3r3,α8=2416r1-240r3,r1=h140,r2=εdΔt20,r3=νΔt20.
Substituting (4.1) into (4.2) we obtain
(4.4)(a+120+ib)δ^n+1=(a-120-ib)δ^n,
where
(4.5)i=-1,a=r1cos3kh+120r1cos2kh+1191r1coskh+1208r1,b=(3r2+3r3)sin3kh+(168r2+72r3)sin2kh+(1470r2+45r3)sinkh.Letδ^n+1=gδ^n and substitute in (4.4) to give
(4.6)g=a-120-iba+120+ib,
where g is the growth factor for the mode. Then, the modulus of the growth factor is
(4.7)|g|=gg-=(a-120)2+b2(a+120)2+b2<1.
Therefore the linearized scheme is unconditionally stable.


## 5. The Initial State

From the initial conditionu(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
(5.1)uN(x,0)=∑j=-1N+1‍δj0ϕj(x),
where δj0 are unknown parameters to be determined. To do this we require uN(x,0) to satisfy the following constraints:(a)
it must agree with the initial conditionu(x,0) at the knots; (3.3) leads to N+1 conditions, that is, uN(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′(x0,0)=u′(xN,0)=0.The initial vectorδ0 is then determined as the solution of a matrix equation derived from (3.3)–(3.6)
(5.2)Mδ0=b.
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 discreteL2- and L∞-error norms.(a)
Burgers’ equation has the analytic solution [23]
(6.1)u(x,t)=x/t1+(t/t0)exp(x2/4νt),0≤x≤1,t≥1,where t0=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)
Figures1 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 L2- 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.(iii)
In this case, we takeh=0.01, Δt=0.01, and viscosity coefficients ν=0.005, ν=0.01. 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]
(6.2)u(x,t)=12[1-tanh{14ν(x-12t-15)}].Table 2
Error normL2×103 for problem (a).
ν
=
0.005
ν
=
0.0005
t
=
1.7
t
=
2.4
t
=
3.1
t
=
1.75
t
=
2.5
t
=
3.25
Present
1.206
0.653
0.434
0.581
0.252
0.179
[14]
1.003
0.598
0.522
0.567
0.308
0.291
[9]
2.281
1.175
0.808
2.023
2.094
0.933Table 3
Error normL∞×103 for problem (a).
ν
=
0.005
ν
=
0.0005
t
=
1.7
t
=
2.4
t
=
3.1
t
=
1.75
t
=
2.5
t
=
3.25
Present
3.868
1.679
1.295
5.371
1.619
0.821
[14]
3.730
2.029
1.537
5.880
2.705
2.291
[9]
7.150
2.921
2.297
14.033
12.615
8.394Table 4
Comparison of results at different times forν=0.0005 and 0≤x≤1, with h=0.005 and Δt=0.01, problem (a).
x
t
=
1.7
t
=
1.7
t
=
1.7
t
=
2.5
t
=
2.5
t
=
2.5
t
=
3.25
t
=
3.25
t
=
3.25
Exact
Numer
Numer
Exact
Numer
Numer
Exact
Numer
Numer
Present
[15]
Present
[15]
Present
[15]
0.1
0.05882
0.05882
0.05883
0.04000
0.04003
0.04000
0.03077
0.03082
0.03077
0.2
0.11765
0.11765
0.11765
0.08000
0.08000
0.08000
0.06154
0.06154
0.06154
0.3
0.17647
0.17647
0.17648
0.12000
0.12000
0.12001
0.09231
0.09231
0.09231
0.4
0.23529
0.23529
0.23531
0.16000
0.16000
0.16001
0.12308
0.12308
0.12308
0.5
0.29412
0.29412
0.29414
0.20000
0.20000
0.20001
0.15385
0.15385
0.15385
0.6
0.35294
0.35287
0.35296
0.24000
0.24000
0.24001
0.18462
0.18462
0.18462
0.7
0.00000
0.00000
0.00000
0.28000
0.28000
0.28001
0.21538
0.21538
0.21539
0.8
0.00000
0.00000
0.00000
0.00977
0.00994
0.00811
0.24615
0.24615
0.24616
0.9
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.12433
0.12401
0.12358Table 5
Comparison of errors at different times forν=0.005 and 0≤x≤1, with h=0.01, problem (a).
Time
L
2-error
L
2-error
L
∞-error
L
∞-error
Present
[18]
Present
[18]
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.01
Δ
t
=
0.001
1.5
0.00019
0.00770
0.00132
0.04210
2.0
0.00018
0.01154
0.00102
0.05166
2.5
0.00017
0.01431
0.00082
0.05623
3.0
0.00016
0.01654
0.00069
0.05882
3.5
0.00019
0.01872
0.00109
0.06046
4.0
0.00035
0.02255
0.00268
0.07691Table 6
Comparison of errors at different times forν=0.01 and 0≤x≤1, with h=0.01, problem (a).
Time
L
2-error
L
2-error
L
∞-error
L
∞-error
Present
[18]
Present
[18]
Δ
t
=
0.01
Δ
t
=
0.001
Δ
t
=
0.01
Δ
t
=
0.001
1.07
0.00019
0.00207
0.00168
0.01660
1.14
0.00021
0.00333
0.00166
0.02413
1.35
0.00023
0.00613
0.00147
0.03667
1.70
0.00022
0.00938
0.00120
0.04677
2.40
0.00020
0.01376
0.00086
0.05537
3.10
0.00021
0.01689
0.00074
0.05906
3.24
0.00022
0.01746
0.00089
0.05955
3.31
0.00023
0.01775
0.00096
0.05977Figure 1
Example (a):v=0.005,h=0.02,Δt=0.1.Figure 2
Example (a):v=0.0005,h=0.005,Δt=0.01.We take as initial condition (6.2), at t=0, over the range 0≤x≤30, with boundary conditions
(6.3)u(0,t)=1,u(30,t)=0,t>0.Forν=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 L2- 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.(c)
A third test example, consider the particular analytic solution of Burgers’ equation [8, 10]
(6.4)u(x,t)=[α+μ+(α-μ)exp(η)]1+exp(η),0≤x≤1,t≥0,where η=α(x-μt-β)/ν,α,μ, and β are constants.Table 7
Comparison ofL2×103 at different times for 0≤x≤30, with h=0.1, and Δt=0.025, problem (b).
t
=
2
t
=
4
t
=
6
t
=
8
t
=
10
ν
Present
0.002
0.003
0.003
0.004
0.004
1/2
[9]
0.006
0.010
0.014
0.022
0.043
1/2
Present
0.032
0.036
0.038
0.038
0.038
1/8
[9]
0.228
0.227
0.227
0.227
0.226
1/8
Present
0.431
0.433
0.433
0.433
0.433
1/24
[9]
3.431
3.426
3.426
3.426
3.429
1/24Table 8
Comparison ofL∞×103 at different times for 0≤x≤30, with h=0.1, and Δt=0.025, problem (b).
t
=
2
t
=
4
t
=
6
t
=
8
t
=
10
ν
Present
0.001
0.002
0.002
0.003
0.003
1/2
[9]
0.004
0.006
0.007
0.017
0.045
1/2
Present
0.044
0.049
0.050
0.051
0.051
1/8
[9]
0.346
0.339
0.340
0.347
0.356
1/8
Present
0.937
0.939
0.939
0.939
0.939
1/24
[9]
4.340
4.336
4.363
4.402
4.432
1/24The behavior of the numerical solution with different viscosity at different times.
(a)
(b)
(c)The comparison can be made with [8–10, 14, 15], these constants are chosen to have the values α=0.4, β=0.125, μ=0.6, and ν=0.01. The initial condition is obtained from (6.4), when t=0. The boundary conditions are taken as
(6.5)u(0,t)=1,(1,t)=0.2,t≥0.From Figure4, it can be seen that the solution represents a travelling wave, initially situated at x=β, moving to the right with speed μ. The present method behaves very well until the wave encounters the right-hand boundary. The numerical results at t=0.5 are compared with the analytic solution, and the results presented in papers [8–10, 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.Table 9
Comparison of results for various numerical schemes with exact at timet=0.5, for ν=0.01.
x
SGA [10]
PAG [10]
ML [9]
CSCM [14]
CSCM [15]
Present
Exact
h
=
1
/
18
h
=
1
/
18
h
=
1
/
36
h
=
1
/
36
h
=
1
/
36
h
=
1
/
36
Δ
t
=
.
001
Δ
t
=
.
001
Δ
t
=
.
025
Δ
t
=
.
025
Δ
t
=
.
025
Δ
t
=
.
025
.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
.056
1.000
1.000
1.000
1.000
1.000
1.000
1.000
.111
1.000
1.000
1.000
1.000
1.000
1.000
1.000
.167
1.000
1.000
1.000
1.000
1.000
1.000
1.000
.222
0.998
0.999
1.000
1.000
1.000
1.000
1.000
.278
0.991
0.997
0.998
0.999
0.999
0.999
0.998
.333
0.970
0.982
0.985
0.985
0.986
0.986
0.980
.389
0.862
0.850
0.851
0.847
0.850
0.850
0.847
.444
0.461
0.444
0.447
0.452
0.448
0.448
0.452
.500
0.159
0.171
0.238
0.238
0.236
0.238
0.238
.556
0.300
0.286
0.205
0.204
0.204
0.204
0.204
.611
0.194
0.197
0.201
0.200
0.200
0.200
0.200
.667
0.213
0.211
0.200
0.200
0.200
0.200
0.200
.722
0.211
0.210
0.200
0.200
0.200
0.200
0.200
.778
0.188
0.190
0.200
0.200
0.200
0.200
0.200
.833
0.201
0.207
0.200
0.200
0.200
0.200
0.200
.889
0.191
0.193
0.200
0.200
0.200
0.200
0.200
.944
0.203
0.202
0.200
0.200
0.200
0.200
0.200
1.00
0.200
0.200
0.200
0.200
0.200
0.200
0.200Figure 4
Example (c): the numerical solution withv=0.01 at different times.


## 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. TheL2- and L∞-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.



---

*Source: 527467-2012-09-26.xml*