# A Comparison of Finite Difference and Finite Volume Methods with Numerical Simulations: Burgers Equation Model

Authors: Ali Hasan Ali; Ahmed Shawki Jaber; Mustafa T. Yaseen; Mohammed Rasheed; Omer Bazighifan; Taher A. Nofal
Journal: Complexity
 (2022)

---

## Abstract

In this paper, we present an intensive investigation of the finite volume method (FVM) compared to the finite difference methods (FDMs). In order to show the main difference in the way of approaching the solution, we take the Burgers equation and the Buckley–Leverett equation as examples to simulate the previously mentioned methods. On the one hand, we simulate the results of the finite difference methods using the schemes of Lax–Friedrichs and Lax–Wendroff. On the other hand, we apply Godunov’s scheme to simulate the results of the finite volume method. Moreover, we show how starting with a variational formulation of the problem, the finite element technique provides piecewise formulations of functions defined by a collection of grid data points, while the finite difference technique begins with a differential formulation of the problem and continues to discretize the derivatives. Finally, some graphical and numerical comparisons are provided to illustrate and corroborate the differences between these two main methods.


---

## Body

## 1. Introduction

In finite difference methods (FDMs), the differential equation is approximated by a set of algebraic equations (see [1, 11, 12]). These approximations are usually improved with the use of more terms and the use of smaller grid spacing (also known as more nodes). Starting with a variational formulation of the problem, the finite element technique provides piecewise formulations of functions defined by a collection of grid data points. The finite difference technique begins with a differential formulation of the problem and continues to discretize the derivatives [15]. Like any approximation, finite difference schemes have their limitations such as oscillations and diffusion especially at points of discontinuity. Another issue that FDM runs into is that mass is only conserved when the grid spacing goes to zero and they struggle with irregular geometry [2, 10]. Usually, a FDM's grid looks like the following (Figure 1).Figure 1
An FDM's grid.However, in some geometries, it may not be possible to align the grid with perpendicular grid lines which means one might not be able to approximate some derivatives in terms of the values along a constant “i” line. Because of these shortcomings, some have turned to using control volume formulation or finite volume methods (FVMs). This started back in 1960 when it was called the integration method, and it turned out to have many applications in oscillation theory and diffusive convection [4–8]. Most recently, FVMs have been used in computational fluid dynamics because of their ability to be used on unstructured meshes [3]. In general, in the FVM, the differential equation is integrated over a control volume and then discretized using some approximations. Another aspect of the huge advantages of using the FVM is that material (like mass or energy) is conserved because it uses integration over small volumes, and since the flux, which is the flow of material in space with time, of the material at the common side is represented by the same expression, this implies that global conservation is also ensured. In fact, both FDM and FVM have a huge number of applications in many engineering and natural science fields. Several important applications of the FDM can be found in computer vision, image processing, solving diffusion convection, and thermal problems (see [9, 13, 14]). On the other hand, the applications of FVM appear in many fields, and its software can be used in several branches, such as solid mechanics, thermal and electrical analysis, structural analysis, oscillation theory, and mechanical engineering design (see [16–21]). The structure of this paper is organized as follows. In Section 2, we provide the definition and the mechanism of the finite volume method (FVM). In Section 3, we consider a one-dimensional application of the FVM. In Section 4, we use an inviscid Burgers equation to apply the proposed comparison of the mentioned numerical methods. Finally, several simulations for the Burgers equation and one simulation for the Buckley–Leverett equation with the main results are presented in Section 5.


## 2. Finite Volume Method (FVM)

To understand the finite volume method, we need to start with the conservation form of the transport equation, which is given as(1)∂T∂t+∇⋅u⟶T=∇⋅α∇T.This partial differential equation is used on the infinitesimal volumes described in Figure2.Figure 2
A PDE on an infinitesimal volume.Knowing that the infinitesimal discrete volumes are unaffordable, they would need to be of some larger finite size. Now we derive the conservation form for a finite volumeδV bounded by a surface, which we will denote as δS. Applying the integral gives us the following:(2)∫δV∂T∂tdV+∫δV∇⋅u⟶TdV=∫δV∇⋅α∇TdV.If we assume the volume to be fixed in space, we can interchange the order of integration in space and differentiation in time:(3)ddt∫δVTdV+∫δV∇⋅u⟶TdV=∫δV∇⋅α∇TdV.The understanding of the first integral on the left is the time rate of change of theT inside volume δV. The Gauss divergence theorem is necessary for the next step. The theorem says that the outward flux of a vector field through some closed surface is equal to the volume integral of the divergence over the region inside the surface. So, basically, we are going to change our equation from volume integrals to surface integrals.(4)ddt∫δVTdV+∫δSn⟶⋅u⟶TdV=∫δSn⟶⋅α∇TdV,where n⟶ represents the outward part that is normal to the surface. The surface integral on the left describes the “in and out” advection flux across the surface of the volume. The one on the other side describes the diffusive transport of T across the surface. The equation can be summarized as the rate of change of the T value in δV which is equal to the rate of transport T through the surface by advection and diffusive fluxes. Notice that we have not given the surface a definite size or shape, so this equation could be applied to whatever we would like. Note the order of the derivatives as they appear in the differential equation (1) versus those of the integral form (4). In the differential equation, we have a first-order advection term and a second-order diffusion term, but in the integral form, there is a 0-th order advection term and a first-order diffusion term. This is an important change because it helps with problems that have discontinuities of the spatial derivative. Examples of those types are the problems that have shock or hydraulic jump.Now we introduce the averageT in δV which we will refer to as(5)T¯=1δV∫δVTdV.The integral conservation law can now be written as a time evolution equation:(6)δVdT¯dt+∫δSn⟶⋅u⟶TdV=∫δSn⟶⋅α∇TdV.It is important to note that equations (4) and (6) are exact solutions and no approximation has yet been introduced. The approximation will come from(1)
Temporal integration of the equations.(2)
Calculation of the fluxes in space and time.However, the first step is to divide the domain into a discrete number of cells which we will denote asδVj where the cell average is known. After that, we calculate the advection and diffusion fluxes through function reconstruction, followed by evaluation of the integrals, and finally the temporal integration with a forward marching time step.


## 3. Application of the FVM

We consider here only the 1-D case. Therefore, we first break up the domain as shown in Figure3.Figure 3
A sample domain of a 1-D case.The cell centres are denoted by the integer numbers and the cell edges are denoted by the fractions. In this case, the cells are line segments, and the cell “volumes” are reduced to just cell widthsδx (which we will assume as uniform throughout). This means that the flux integrals reduce to the evaluation of the term at the cell edges, which means the conservation law can be rewritten as(7)dTj¯dt=−Fj+1/2−Fj−1/2δx+Dj+1/2−Dj−1/2δx,where F is the advection flux and D is the diffusive flux at the cell edge xj+1/2.The other assumption that we will make is that we are using a piecewise constant approximation as shown in Figure4.Figure 4
A piecewise constant approximation.This means that we have the following approximation:T=a0.We can now write(8)∫xj+1/2xj+1/2a0dx=δxjTj¯,because we only have one unknown coefficient. Basically, the integral of T over a cell gives us δxjTj¯. This means that the solution is a0=Tj¯, which finishes our function reconstruction. Note that at each cell edge, there are two different values. For the advection term, we should use the upstream edge (what goes into the next term has to come from the previous term). This leads to(9)Tj+1/2=T¯j,uj+1/2≥0,T¯j+1,uj+1/2<0.The other implication of this is that piecewise constant functions cannot be used to compute the diffusive term because they have a zero derivative in space. This will not pose a problem for us later because we will be considering Burgers and the Buckley–Leverett equations which do not have a diffusive term.


## 4. Numerical Methods

Take a look at the inviscid Burgers equation(10)ut+u⋅ux=0.Applying the Euler forward in time and backward in space discretization, one would come up with(11)Ujn+1=Ujn−khUjnUjn−Uj−1n.This is sometimes referred to as the nonconservative upwind scheme. This method works well on nice smooth solutions but will not converge at discontinuous points as the grid is refined (see Figure5).Figure 5
A nonconservative upwind method.To avoid this issue, we apply the conservation form(12)Ujn+1=Ujn−khFUj−pn,Uj−p+1n,…,Uj+qn−FUj−p−1n,Uj−pn,…,Uj+q−1n,where F is the numerical flux function. In general, methods that follow this type of scheme are called conservative methods. So, all the methods to follow were implemented using this type of scheme.For the comparisons, we used the general conservation law:(13)ut+fux=0,where for the Burgers equation, we use(14)fu=12u2,and for the Buckley–Leverett equation, we use(15)fu=u2u2+1/41−u2.The basic setup for the two different kinds of methods is as follows.


### 4.1. Finite Difference Methods

(1)
Set up grid (or nodes).(2)
Approximate and replace all the derivatives by their corresponding finite difference approximations.(3)
Solve the resulting algebraic system of equations.


### 4.2. Finite Volume Methods

(1)
Grid generation.(2)
Discretize the integral conservation equation on control volumes.(3)
Solve the resultant discrete integral/flux equations.


### 4.3. Upwind Conservative Method

In fact, this method is first-order accurate in both time and space; however, it is only stable in the interval0≤ak/h≤1, where k is the size of the time step and h is the size of the space step.Then, applying the standard finite difference upwind definition, we would get(16)Ujn+1=Ujn−khfUjn−fUj−1n.Then, applying this to the Burgers equation, we get(17)Ujn+1=Ujn−kh12Ujn2−12Uj−1n2.Moreover, we apply this method to the Buckley–Leverett equations to get(18)Ujn+1=Ujn−khUjn2Ujn2+1/41−Ujn2−Uj−1n2Uj−1n2+1/41−Uj−1n2.


### 4.4. Lax–Friedrichs Method

Like the upwind conservative method, the Lax–Friedrichs method is first-order accurate in time and space. In addition, it has a stability requirement ofak/h≤1. The Lax–Friedrichs method takes the form(19)Ujn+1=12Uj−1n+Uj+1n−k2hfUj+1n−fUj−1n.Then, applying this to the Burgers equation, we get(20)Ujn+1=12Uj−1n+Uj+1n−kh12Uj+1n2−12Uj−1n2.For the Buckley–Leverett equation, we would have(21)Ujn+1=12Uj−1n+Uj+1n−k2hUj+1n2Uj+1n2+1/41−Uj+1n2−Uj−1n2Uj−1n2+1/41−Uj−1n2.


### 4.5. Lax–Wendroff Method

The Lax–Wendroff method is second-order accurate in both time and space, and it has a known stability ofak/h≤1. For the nonlinear conservation laws, it takes the form(22)Ujn+1=Uj+1n−k2hfUj+1n−fUj−1n+k22h2Aj+1/2fUj+1n−fUjn−Aj−1/2fUjn−fUj−1n.


### 4.6. Godunov’s Scheme

To compare the previous finite difference schemes with a finite volume scheme, we choose to use Godunov's scheme. Although this method is only first-order accurate, it is the basis of other high-order finite volume schemes. To get Godunov's method, we start by lettingUjn be a numerical solution on the n-th layer. After that, we need to define a function U^nx,t for the time interval tn<t<tn+1. At t=tn, we have the following:(23)U^nx,t=Ujn,xj−h2<x<xj+h2,forj=2,…,n−1.We also need to defineU^x,t to be the solution of the collection of Riemann problems on that interval. Then, the solution to the next value, Ujn+1, is simply defined by averaging U^x,tn+1 over the interval xj−h/2<x<xj+h/2. When everything is said and done, this idea reduces to the following:(24)Ujn+1=Ujn−khFUjn,Uj+1n−FUj−1n,Ujn.When we apply this to the Burgers equation, the numerical fluxF is(25)FU,V=u∗22,and when we apply this to the Buckley–Leverett equation,(26)FU,V=u∗2u∗2+1/41−u∗2,where u∗ is defined as follows.IfU≥V, then(27)u∗=U,U+V2>0,V,else.IfU<V, then(28)u∗=U,U>0,V,V<0,0,U≤0≤V.


## 4.1. Finite Difference Methods

(1)
Set up grid (or nodes).(2)
Approximate and replace all the derivatives by their corresponding finite difference approximations.(3)
Solve the resulting algebraic system of equations.


## 4.2. Finite Volume Methods

(1)
Grid generation.(2)
Discretize the integral conservation equation on control volumes.(3)
Solve the resultant discrete integral/flux equations.


## 4.3. Upwind Conservative Method

In fact, this method is first-order accurate in both time and space; however, it is only stable in the interval0≤ak/h≤1, where k is the size of the time step and h is the size of the space step.Then, applying the standard finite difference upwind definition, we would get(16)Ujn+1=Ujn−khfUjn−fUj−1n.Then, applying this to the Burgers equation, we get(17)Ujn+1=Ujn−kh12Ujn2−12Uj−1n2.Moreover, we apply this method to the Buckley–Leverett equations to get(18)Ujn+1=Ujn−khUjn2Ujn2+1/41−Ujn2−Uj−1n2Uj−1n2+1/41−Uj−1n2.


## 4.4. Lax–Friedrichs Method

Like the upwind conservative method, the Lax–Friedrichs method is first-order accurate in time and space. In addition, it has a stability requirement ofak/h≤1. The Lax–Friedrichs method takes the form(19)Ujn+1=12Uj−1n+Uj+1n−k2hfUj+1n−fUj−1n.Then, applying this to the Burgers equation, we get(20)Ujn+1=12Uj−1n+Uj+1n−kh12Uj+1n2−12Uj−1n2.For the Buckley–Leverett equation, we would have(21)Ujn+1=12Uj−1n+Uj+1n−k2hUj+1n2Uj+1n2+1/41−Uj+1n2−Uj−1n2Uj−1n2+1/41−Uj−1n2.


## 4.5. Lax–Wendroff Method

The Lax–Wendroff method is second-order accurate in both time and space, and it has a known stability ofak/h≤1. For the nonlinear conservation laws, it takes the form(22)Ujn+1=Uj+1n−k2hfUj+1n−fUj−1n+k22h2Aj+1/2fUj+1n−fUjn−Aj−1/2fUjn−fUj−1n.


## 4.6. Godunov’s Scheme

To compare the previous finite difference schemes with a finite volume scheme, we choose to use Godunov's scheme. Although this method is only first-order accurate, it is the basis of other high-order finite volume schemes. To get Godunov's method, we start by lettingUjn be a numerical solution on the n-th layer. After that, we need to define a function U^nx,t for the time interval tn<t<tn+1. At t=tn, we have the following:(23)U^nx,t=Ujn,xj−h2<x<xj+h2,forj=2,…,n−1.We also need to defineU^x,t to be the solution of the collection of Riemann problems on that interval. Then, the solution to the next value, Ujn+1, is simply defined by averaging U^x,tn+1 over the interval xj−h/2<x<xj+h/2. When everything is said and done, this idea reduces to the following:(24)Ujn+1=Ujn−khFUjn,Uj+1n−FUj−1n,Ujn.When we apply this to the Burgers equation, the numerical fluxF is(25)FU,V=u∗22,and when we apply this to the Buckley–Leverett equation,(26)FU,V=u∗2u∗2+1/41−u∗2,where u∗ is defined as follows.IfU≥V, then(27)u∗=U,U+V2>0,V,else.IfU<V, then(28)u∗=U,U>0,V,V<0,0,U≤0≤V.


## 5. Results

In this section, several simulations for the Burgers equation and one simulation for the Buckley–Leverett equation are carried out. After we ran each simulation, we compared the results from the finite difference methods and the finite volume method.


### 5.1. Simulation 1: Single Shock

The equation and the solutions of the single shock simulation are given by equation (29) and shown in Figure 6 using the methods mentioned in Section 4.(29)ux,0=1,x<0.1,0,x>0.1.Figure 6
Graphs for the single shock simulation.To satisfy the stability condition for the finite schemes, we usedh=1/60 and k=0.005.We can see from the results that the finite difference methods had issues at the sharp turns (the derivative does not exist). Notice that the Lax–Friedrichs method shows diffusion at the corners while the Lax–Wendroff method produces oscillations. The Godunov method handles the shock much better, but still it is not accurate. In fact, we also wanted to point out that for the Burgers equation, the classical upwind scheme and the Godunov scheme produce the same result (see Figure7). This is basically because f′u>0.Figure 7
Comparing the Godunov method to the upwind method.


### 5.2. Simulation 2: Traveling Pulse

The equation and the solutions of the traveling pulse simulation are given by equation (30) and shown in Figure 8 using the methods mentioned in Section 4.(30)ux,0=0,x<−1,12,−1<x<0,0,x>0.Figure 8
Graphs for the traveling pulse simulation.In the simulation above, we used a smallerh value so that we could better see what happens at the corners; h=0.025 and k=0.01. We obtain the same results that we saw before with the single shock. It is interesting to note that the Lax–Wendroff method does not have oscillations at the top left corner x=0, but it does have oscillations at the other places where the derivative does not exist. Note that we did not include the upwind scheme because we knew that it would produce the same results as the Godunov scheme.


### 5.3. Simulation 3: Buckley–Leverett Equations

The equation and the solutions of Buckley–Leverett simulation are given by equation (31) and shown in Figure 9 using the methods mentioned in Section 4.(31)ux,0=12,x=0,0,x>0.Figure 9
Graphs for the Buckley–Leverett simulation.This produced the same results that we saw in the Burgers equation simulations. The Lax–Wendroff method does not have an oscillation by itself, but it definitely has its issues. To see the difference between the Godunov scheme and the upwind scheme, we had to use a different flux (the Buckley–Leverett flux) as shown in Figure10.Figure 10
Comparing Godunov and upwind methods.


## 5.1. Simulation 1: Single Shock

The equation and the solutions of the single shock simulation are given by equation (29) and shown in Figure 6 using the methods mentioned in Section 4.(29)ux,0=1,x<0.1,0,x>0.1.Figure 6
Graphs for the single shock simulation.To satisfy the stability condition for the finite schemes, we usedh=1/60 and k=0.005.We can see from the results that the finite difference methods had issues at the sharp turns (the derivative does not exist). Notice that the Lax–Friedrichs method shows diffusion at the corners while the Lax–Wendroff method produces oscillations. The Godunov method handles the shock much better, but still it is not accurate. In fact, we also wanted to point out that for the Burgers equation, the classical upwind scheme and the Godunov scheme produce the same result (see Figure7). This is basically because f′u>0.Figure 7
Comparing the Godunov method to the upwind method.


## 5.2. Simulation 2: Traveling Pulse

The equation and the solutions of the traveling pulse simulation are given by equation (30) and shown in Figure 8 using the methods mentioned in Section 4.(30)ux,0=0,x<−1,12,−1<x<0,0,x>0.Figure 8
Graphs for the traveling pulse simulation.In the simulation above, we used a smallerh value so that we could better see what happens at the corners; h=0.025 and k=0.01. We obtain the same results that we saw before with the single shock. It is interesting to note that the Lax–Wendroff method does not have oscillations at the top left corner x=0, but it does have oscillations at the other places where the derivative does not exist. Note that we did not include the upwind scheme because we knew that it would produce the same results as the Godunov scheme.


## 5.3. Simulation 3: Buckley–Leverett Equations

The equation and the solutions of Buckley–Leverett simulation are given by equation (31) and shown in Figure 9 using the methods mentioned in Section 4.(31)ux,0=12,x=0,0,x>0.Figure 9
Graphs for the Buckley–Leverett simulation.This produced the same results that we saw in the Burgers equation simulations. The Lax–Wendroff method does not have an oscillation by itself, but it definitely has its issues. To see the difference between the Godunov scheme and the upwind scheme, we had to use a different flux (the Buckley–Leverett flux) as shown in Figure10.Figure 10
Comparing Godunov and upwind methods.


## 6. Conclusion

Although we used the conservation equations to explain and illustrate the finite volume method, it can be used on any type of differential equation. We choose the conservation equations since usually finite volume methods are used on fluid flow, which is conservative. We can conclude by our comparison that the main difference of these methods can be seen in the way of approaching the solutions, as starting with a variational formulation of the problem, the finite element technique provides piecewise formulations of functions defined by a collection of grid data points, while the finite difference technique begins with a differential formulation of the problem and continues to discretize the derivatives. Moreover, the finite volume method even with the simplifications that we gave produced very good results for the simulations. In fact, it did not appear to have any difficulty handling the shocks and discontinuous derivatives, which is why FVMs are used in modelling fluid flow.



---

*Source: 9367638-2022-06-27.xml*