# A Multilevel Finite Difference Scheme for One-Dimensional Burgers Equation Derived from the Lattice Boltzmann Method

Authors: Qiaojie Li; Zhoushun Zheng; Shuang Wang; Jiankang Liu
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

An explicit finite difference scheme for one-dimensional Burgers equation is derived from the lattice Boltzmann method. The system of the lattice Boltzmann equations for the distribution of the fictitious particles is rewritten as a three-level finite difference equation. The scheme is monotonic and satisfies maximum value principle; therefore, the stability is proved. Numerical solutions have been compared with the exact solutions reported in previous studies. TheL2, L∞ and Root-Mean-Square (RMS) errors in the solutions show that the scheme is accurate and effective.


---

## Body

## 1. Introduction

The lattice Boltzmann method (LBM) has been introduced as a new computational tool for the study of fluid dynamics and systems governed by partial differential equations. It has made a rapid development in theory and application over the last couple of decades since its inception [1–4]. This method can be either regarded as an extension of the lattice gas automaton [5] or as a special discrete form of the Boltzmann equation for kinetic theory [6]. The lattice Boltzmann models can also be used as partial differential equation (PDE) solvers. By choosing appropriate collision operator or equilibrium distribution, the lattice Boltzmann model is able to recover the PDE of interest. Recently, it has been developed to simulate linear and nonlinear PDE such as Laplace equation [7], Poisson equation [8, 9], the shallow water equation [10], Burgers equation [11], Korteweg-de Vires equation [12], Wave equation [13, 14], reaction-diffusion equation [15, 16], and convection-diffusion equation [17, 18].The numerical schemes based on the LBM are given as a system of two-level explicit difference equations composed of the distribution functions of fictitious particles for each direction in which the particles move. For one-dimensional advection-diffusion problems, Ancona [19] showed that the LB schemes with the velocity model D1Q2 which includes two velocities with speed 1 in opposite directions to each other can be rewritten as the DuFort-Frankel scheme [20] which is a second-order three-level difference scheme. This shows that the accuracy of the LB schemes based on the model D1Q2 is identical to that of the DuFort-Frankel scheme. Suga [21] have proposed a four-level explicit finite difference scheme for 1D diffusion equation which is derived from the lattice Boltzmann method with rest particles. The consistency analysis of the scheme shows that the two parameters which appear in the scheme, the relaxation parameter and the amount of rest particles, can be determined such that the scheme has the truncation error of fourth order. In spite of the vast and successful applications, the numerical stability of the method has not been well understood. For certain specific class of lattice Boltzmann methods, for example, solving for linear and nonlinear convective-diffusive equation, there are some convergence and stability results given by Elton et al. [22].Many works have been developed on lattice Boltzmann method to the Burgers equation in one or higher dimension [23–25]. In those papers, the standard lattice Boltzmann method was used and the macroscopic quantities were computed by the distribution function. However, those models are suffered from the stability. In this paper, we derive a three-level difference scheme for 1D Burgers equation based on the model D1Q2 from the LB schemes. It is generally recognized that LBM is a finite difference scheme of Boltzmann equation that has higher-order discretization error. We develop this method with the point of view above, but, at the same time, we also regard the LBM with BGK model as finite difference method for macroscopic equation. We find such LB scheme is a three-level finite difference one, which is monotonic and satisfies maximum value principle; therefore, we complete the proof of stability.The rest of the paper is organized as follows. Section2 describes the LB scheme with the velocity model D1Q2 and derives the three-level finite difference scheme which is equivalent to the LB scheme. A stability analysis of the scheme is given in Section 3. In Section 4, numerical solutions are compared with exact solutions reported in previous studies. And the conclusions are given in the end.


## 2. The Three-Level Finite Difference Scheme for 1D Burgers Equation Based on the LB Schemes

The one-dimensional Burgers equation take the following form:(2.1)∂u∂t+u∂u∂x=ν∂2u∂x2,
with the initial condition u(x,0)=u0(x). Here, the viscous coefficient ν=1/Re, Re is the Reynolds number. Historically, (2.1) was first introduced by Bateman [26] who gave its steady solutions. It was later treated by Burgers [27] as a mathematical model for turbulence and after whom such an equation is widely referred to as Burgers equation. For a small value of ν, Burgers equation behaves merely as hyperbolic partial differential equation and the problem becomes very difficult to solve as a steep shock-like wave fronts developed.


### 2.1. The Lattice Boltzmann Scheme

According to the theory of the LBM, it consists of two steps:(1) streaming, where each particle moves to the nearest node in the direction of its velocity; (2) colliding, which occurs when particles arriving at a node interact and possibly change their velocity directions according to scattering rules. Fictitious particles are introduced at each of the mesh points x=jΔx(j=…,-2,-1,0,1,2,…), and they move with the velocity ci determined by the D1Q2 model from x to the neighboring mesh point which was shown in Figure 1. The lattice Boltzmann schemes are established on grids with two directions(2.2)[c1,c-1]=[-c,c],
where c=Δx/Δt is the speed in the system. Let fi(x,t) denote the distribution function of the particles moving with velocity ci. So the time evolution of the distribution function fi(x,t) is given by the following lattice Boltzmann equation (LBE) based on the Bhatnagar-Gross-Krook (BGK) model:(2.3)fi(x+ciΔt,t+Δt)=fi(x,t)-1τ(fi(x,t)-fieq(x,t)),
where fieq(x,t) is the local equilibrium distribution function of particles and τ is the dimensionless relaxation time which controls the rate of approach to equilibrium. The change in the distribution function produced by the collision of particles is approximated by the second term on the right-hand side of (2.3). The macroscopic velocity u(x,t) is defined in terms of the distribution function as(2.4)u(x,t)=∑ifi(x,t)=∑ifieq(x,t).Figure 1
D1Q2 model with two velocities in one dimension.In this paper,fieq(x,t) are determined as to satisfy (2.4) and the following conditions:(2.5)∑icifieq(x,t)=u2(x,t)2,∑icicifieq(x,t)=c2u(x,t).Solving these equations determines the equilibrium distribution functions(2.6)f1eq(x,t)=u(x,t)2+u2(x,t)Δt4Δx,f-1eq(x,t)=u(x,t)2-u2(x,t)Δt4Δx.
Applying the Chapman-Enskog expansion [24] yields the above Burgers equation (2.1) from the LBE and the equilibrium distribution functions given by (2.3) and (2.6), respectively. The viscosity ν is defined by ν=(τ-1/2)(Δx2/Δt).


### 2.2. The Multilevel Finite Difference Scheme

Now, we letfi,jn denote fi(jΔx,nΔt) and let ujn denote u(jΔx,nΔt). We note that the subscript i,j combines information about the channel or direction of propagation (i=1,-1) and location (j denotes a grid node). Using the equilibrium distribution function (2.6), the lattice Boltzmann equation (2.3) can be rewritten by classical finite different notation(2.7)f1,j+1n+1=(1-1τ)f1,jn+12τujn+Δt4τΔx(ujn)2,(2.8)f-1,j-1n+1=(1-1τ)f-1,jn+12τujn-Δt4τΔx(ujn)2.According to (2.4), the macroscopic velocity can be computed by(2.9)ujn+1=f1,jn+1+f-1,jn+1=(1-1τ)(f1,j-1n+f-1,j+1n)+12τ(uj-1n+uj+1n)+Δt4τΔx((uj-1n)2-(uj+1n)2)=H(f1,j-1n,f-1,j+1n,uj-1n,uj+1n).In addition,(2.10)f1,j-1n+f-1,j+1n=(uj-1n-f-1,j-1n)+(uj+1n-f1,j+1n)=uj-1n+uj+1n-(f1,j+1n+f-1,j-1n),
while(2.11)f1,j+1n+f-1,j-1n=(1-1τ)f1,jn-1+12τujn-1+Δt4τΔx(ujn-1)2+(1-1τ)f-1,jn-1+12τujn-1-Δt4τΔx(ujn-1)2=(1-1τ)(f1,jn-1+f-1,jn-1)+1τujn-1=(1-1τ)ujn-1+1τujn-1=ujn-1.Then, (2.10) becomes(2.12)f1,j-1n+f-1,j+1n=uj-1n+uj+1n-ujn-1.Substitute (2.12) to (2.9), we finally obtain the following three-level explicit finite difference scheme(2.13)ujn+1=(1-1τ)(uj-1n+uj+1n-ujn-1)+12τ(uj-1n+uj+1n)+Δt4τΔx((uj-1n)2-(uj+1n)2).


## 2.1. The Lattice Boltzmann Scheme

According to the theory of the LBM, it consists of two steps:(1) streaming, where each particle moves to the nearest node in the direction of its velocity; (2) colliding, which occurs when particles arriving at a node interact and possibly change their velocity directions according to scattering rules. Fictitious particles are introduced at each of the mesh points x=jΔx(j=…,-2,-1,0,1,2,…), and they move with the velocity ci determined by the D1Q2 model from x to the neighboring mesh point which was shown in Figure 1. The lattice Boltzmann schemes are established on grids with two directions(2.2)[c1,c-1]=[-c,c],
where c=Δx/Δt is the speed in the system. Let fi(x,t) denote the distribution function of the particles moving with velocity ci. So the time evolution of the distribution function fi(x,t) is given by the following lattice Boltzmann equation (LBE) based on the Bhatnagar-Gross-Krook (BGK) model:(2.3)fi(x+ciΔt,t+Δt)=fi(x,t)-1τ(fi(x,t)-fieq(x,t)),
where fieq(x,t) is the local equilibrium distribution function of particles and τ is the dimensionless relaxation time which controls the rate of approach to equilibrium. The change in the distribution function produced by the collision of particles is approximated by the second term on the right-hand side of (2.3). The macroscopic velocity u(x,t) is defined in terms of the distribution function as(2.4)u(x,t)=∑ifi(x,t)=∑ifieq(x,t).Figure 1
D1Q2 model with two velocities in one dimension.In this paper,fieq(x,t) are determined as to satisfy (2.4) and the following conditions:(2.5)∑icifieq(x,t)=u2(x,t)2,∑icicifieq(x,t)=c2u(x,t).Solving these equations determines the equilibrium distribution functions(2.6)f1eq(x,t)=u(x,t)2+u2(x,t)Δt4Δx,f-1eq(x,t)=u(x,t)2-u2(x,t)Δt4Δx.
Applying the Chapman-Enskog expansion [24] yields the above Burgers equation (2.1) from the LBE and the equilibrium distribution functions given by (2.3) and (2.6), respectively. The viscosity ν is defined by ν=(τ-1/2)(Δx2/Δt).


## 2.2. The Multilevel Finite Difference Scheme

Now, we letfi,jn denote fi(jΔx,nΔt) and let ujn denote u(jΔx,nΔt). We note that the subscript i,j combines information about the channel or direction of propagation (i=1,-1) and location (j denotes a grid node). Using the equilibrium distribution function (2.6), the lattice Boltzmann equation (2.3) can be rewritten by classical finite different notation(2.7)f1,j+1n+1=(1-1τ)f1,jn+12τujn+Δt4τΔx(ujn)2,(2.8)f-1,j-1n+1=(1-1τ)f-1,jn+12τujn-Δt4τΔx(ujn)2.According to (2.4), the macroscopic velocity can be computed by(2.9)ujn+1=f1,jn+1+f-1,jn+1=(1-1τ)(f1,j-1n+f-1,j+1n)+12τ(uj-1n+uj+1n)+Δt4τΔx((uj-1n)2-(uj+1n)2)=H(f1,j-1n,f-1,j+1n,uj-1n,uj+1n).In addition,(2.10)f1,j-1n+f-1,j+1n=(uj-1n-f-1,j-1n)+(uj+1n-f1,j+1n)=uj-1n+uj+1n-(f1,j+1n+f-1,j-1n),
while(2.11)f1,j+1n+f-1,j-1n=(1-1τ)f1,jn-1+12τujn-1+Δt4τΔx(ujn-1)2+(1-1τ)f-1,jn-1+12τujn-1-Δt4τΔx(ujn-1)2=(1-1τ)(f1,jn-1+f-1,jn-1)+1τujn-1=(1-1τ)ujn-1+1τujn-1=ujn-1.Then, (2.10) becomes(2.12)f1,j-1n+f-1,j+1n=uj-1n+uj+1n-ujn-1.Substitute (2.12) to (2.9), we finally obtain the following three-level explicit finite difference scheme(2.13)ujn+1=(1-1τ)(uj-1n+uj+1n-ujn-1)+12τ(uj-1n+uj+1n)+Δt4τΔx((uj-1n)2-(uj+1n)2).


## 3. Stability Analysis

In this section, assumed the initial valueu0(x) is bounded and smooth enough, we will prove the multilevel finite difference scheme is stable in L1⋂L∞ space. Suppose(3.1)u0(x)∈L1,|u0(x)|≤1.It is not difficult to see that, if|ujn|≤1 and(3.2)τ≥1,ΔtΔx≤1,
then the scheme (2.9) is monotonic increase. τ≥1 means(3.3)νΔtΔx2≥12.Now, we will point out that the solution of the scheme (2.13) satisfies the maximum value principle.Lemma 3.1 (maximum value principle).
If initial value|u0(x)|≤1 and the restrictions (3.2) hold, then, for all j∈Z, there are
(3.4)min⁡lul0≤ujn+1≤max⁡lul0,n≥0.Proof.
It is known that if we takef1,j0=uj0/2,f-1,j0=uj0/2, and uLn=maxjujn,uSn=minjujnj∈Z, then, for all j,k∈Z,
(3.5)f1,j1+f-1,k1=H(f1,j-10,f-1,k+10,uj-10,uk+10)=H(uj-102,uk+102,uj-10,uk+10)≤H(uL02,uL02,uL0,uL0)=(1-1τ)(uL02+uL02)+12τ(uL0+uL0)+Δt4τΔx((uL0)2-(uL0)2)=uL0,
and similarly
(3.6)f1,j1+f-1,k1=H(uj-102,uk+102,uj-10,uk+10)≥H(uS02,uS02,uS0,uS0)=uS0.
If we supposeuS0≤f1,jn+f-1,kn≤uL0 is also correct. Particularly j=k, we have uS0≤ujn≤uL0, then
(3.7)f1,jn+1+f-1,kn+1=H(f1,j-1n,f-1,k+1n,uj-1n,uk+1n)≤H(f1,j-1n,f-1,k+1n,uL0,uL0)=(1-1τ)(f1,j-1n+f-1,k+1n)+1τuL0≤uL0.
Similarly, we get(3.8)f1,jn+1+f-1,kn+1≥uS0.
Let j=k, we can get
(3.9)min⁡lul0≤ujn+1≤max⁡lul0,n≥0.
Assume thatũ(x,t) is another solution of (2.1) with subject to initial condition ũ(x,0)=ũ0(x), and the initial condition satisfies |ũ0(x)|≤1. Using the same scheme (2.13) and same restriction condition (3.2), we have the following.Lemma 3.2.
If the conditions of Lemma3.1 are fulfilled, there are inequalities
(3.10)∑jmax⁡(ujn+1,ũjn+1)≤∑jmax⁡(uj0,ũj0),∑jmin⁡(ujn+1,ũjn+1)≥∑jmin⁡(uj0,ũj0).
Denote thatuΔxn={ujn,j∈Z} is the discrete solution of LBE (2.7)–(2.9) at time nΔt, and ∥uΔxn∥L1=∑j|ujn|Δx is the L1 norm of discrete function uΔxn. Then, the solution is stable in the meaning of L1.Theorem 3.3.
IfuΔxn,ũΔxn are the solutions of (2.13), uΔx0,ũΔx0∈L1(R2) with subject to the corresponding initial conditions (3.1) and restrictions (3.2), then there are
(3.11)‖uΔxn-ũΔxn‖L1≤‖uΔx0-ũΔx0‖L1,(3.12)‖uΔxn‖L1≤‖uΔx0‖L1.Proof.
Consider(3.13)|ujn+1-ũjn+1|=max⁡(ujn+1,ũjn+1)-min⁡(ujn+1,ũjn+1).
Summing the absolute value to all j, by Lemma 3.2, we have
(3.14)∑j|ujn+1-ũjn+1|=∑jmax⁡(ujn+1,ũjn+1)-∑jmin⁡(ujn+1,ũjn+1)≤∑jmax⁡(uj0,ũj0)-∑jmin⁡(uj0,ũj0)=∑j|uj0-ũj0|.
If we let ũΔx(x,t)=0 in (3.11), we can get (3.12).Remark 3.4.
The restriction (3.2) is sufficient but not necessary.


## 4. Numerical Experiments

Example 4.1.
We investigate the accuracy of the scheme by solving (2.1) on the domain (t,x)∈(0,T]×[0,1]. The initial condition is u(x,0)=sin(2πx),0≤x≤1, and the homogenous boundary condition is u(0,t)=u(1,t)=0. In this case, the exact Fourier solution is given by [28]
(4.1)u(x,t)=2πν∑n=1∞anexp⁡(-n2π2νt)nsin⁡(nπx)a0+∑n=1∞anexp⁡(-n2π2νt)cos⁡(nπx),
where
(4.2)a0=∫01exp⁡(-(2πν)-1(1-cos⁡(πx)))dx,an=2∫01exp⁡(-(2πν)-1(1-cos⁡(πx)))cos⁡(nπx)dx,n=1,2,….
In comparison with the analytical solutions, the efficiency of proposed model is validated. The following error norms are used to measure the accuracy:(1)
L2-error
(4.3)‖e‖L2=(∑i=1nei2)1/2,(2)
L∞-error
(4.4)‖e‖L∞=Max⁡|ei|,1≤i≤n,(3)
The root mean square (RMS) error(4.5)‖e‖RMS=(∑i=1nei2n)1/2.
The numerical solutions of (2.1), which are computed by using different step size at time T=0.1 for ν=1, are given in Table 1. The above error norms are given in Table 2 for different mesh size.
From Table2, we find that the accuracy measured in L2,L∞ and RMS norm errors increases as the step size decrease. The numerical solutions are in the symmetric pattern as the exact solutions are. Table 3 and Figure 1 show a comparison between numerical and exact solutions at different times for ν=0.005. The curves for distribution of absolute errors at different times are also shown in Figure 2. It is known that the Fourier solutions for ν≤0.001 fail to converge because of the slow convergence of the infinite series [28]. The numerical solution cures for ν=0.001 at different time are drawn in Figure 3, which shows the correct physical behavior.Table 1
Comparison of the LB finite difference solutions with exact solution atT=0.1 for ν=1 with τ=1.
xNumerical solutionExact solutionN=10N=20N=1000.10.008470.010590.011290.011320.20.013700.017150.018280.018330.30.013710.017160.018300.018350.40.008480.010610.011320.011350.50.000000.000000.000000.000000.6−0.00848−0.01061−0.01132−0.011350.7−0.01371−0.01716−0.01830−0.018350.8−0.01370−0.01715−0.01829−0.018330.9−0.00847−0.01059−0.01129−0.01132Table 2
Error norms forν=1 at T=0.1 with different step size.
N∥e∥L2∥e∥L∞∥e∥RMS101.089E-022.789E-031.125E-04204.640E-031.190E-035.000E-051003.631E-039.296E-043.756E-05Table 3
Comparison of the LB finite difference solutions with exact solution forν=0.005 with dx=0.005,dt=0.003, and τ=1.1 at different times.
xt1.42.02.6NumericalExactNumericalExactNumericalExact0.10.063030.063940.045670.046210.035810.036180.20.119750.127840.091330.092410.071620.072340.30.189020.191680.136940.138540.107170.108260.40.250910.254340.178090.180220.133670.135210.50.000000.000000.000000.000000.000000.000000.6−0.25091−0.25434−0.17809−0.18022−0.13367−0.135210.7−0.18902−0.19168−0.13694−0.13854−0.10717−0.108260.8−0.12605−0.12784−0.09133−0.09241−0.07162−0.072340.9−0.06303−0.06394−0.04567−0.04621−0.03581−0.03618Numerical solutions (a) and distribution of absolute errors (b) forν=0.005 at different times with dx=0.005,τ=1.1, and dt=0.003.
(a)
Numerical solutions(b)
Absolute errorsFigure 3
Numerical solutions forν=0.001, at different times with dx=0.001,τ=1 and dt=0.0005.Example 4.2.
Consider Burgers equation with the following forms:(4.6)∂u∂t+u∂u∂x=1Re∂2u∂x2,12≤x≤32,t>0,u(x,0)=1Re(x+tan(x2)),12≤x≤32,u(12,t)=1Re+t[12+tan(Re4(Re+t))],t>0,u(32,t)=1Re+t[32+tan(3Re4(Re+t))],t>0.
It possesses the exact solution [23]
(4.7)u(x,t)=1Re+t[x+tan(xRe2(Re+t))].
In the computation, we compare the result with theD1Q2 and D1Q3 lattice Boltzmann model whose equilibrium distribution functions are taken as
(4.8)f1eq(x,t)=u(x,t)2+u2(x,t)4c,f2eq(x,t)=u(x,t)2-u2(x,t)4c,f0eq(x,t)=23u(x,t),f1eq(x,t)=u(x,t)6+u2(x,t)4c,f2eq(x,t)=u(x,t)6-u2(x,t)4c.
LetRe=500, we give the results of our model, and exact solution as Figure 4 at t=0.4. Table 4 shows the results of the D1Q2, D1Q3, our model and the exact solution at different lattice at time t=0.4. The global relative errors
(4.9)GRE=∑i|uE(xi,t)-uN(xi,t)|∑i|uN(xi,t)|,
which are used to measure the accuracy are presented in Table 5.
From Figure4 and Table 4, we find that the D1Q2,D1Q3, and our model are all in excellent agreement with the exact solutions. The accuracy of the multilevel finite difference model is even higher than the D1Q2 and D1Q3 model. It should be pointed out that in order to attain better accuracy, the LB model requires a relatively small time step Δt but the multilevel finite difference model does not have this restriction.Table 4
Comparison of the results with D1Q2, D1Q3, our model, and exact solution.
xD1Q2 modelD1Q3 modelOur modelExact solution0.50.0015000.0015000.0015000.0015000.60.0017950.0017870.0017880.0017860.70.0021120.0020990.0021030.0020960.80.0024310.0024140.0024200.0024110.90.0027550.0027340.0027420.0027311.00.0030860.0030600.0030690.0030561.10.0034250.0033930.0034020.0033891.20.0037730.0037350.0037420.0037291.30.0041310.0040870.0040920.0040801.40.0045110.0044830.0044510.0044211.50.0050000.0050000.0050000.005000Table 5
Global relative errors with different models.
D1Q2 modelD1Q3 modelOur modelGRE3.2383E-031.7094E-035.8823E-04Figure 4
Comparison of the exact solution and our model. Parameters are:Re=500,dx=0.01,dt=0.002,τ=1.


## 5. Conclusion

In the current study, a three-level explicit finite difference scheme for 1D Burgers equation is derived by rewriting the LB scheme. Furthermore, it is proved that the scheme is conditionally stable. The efficiency and accuracy of the proposed scheme are validated through detail numerical simulation. It can be found that the numerical solutions are in excellent agreement with the analytical solutions. In order to derive LB scheme in a higher dimension, the LBM with the multispeed velocity model will be useful, in which different free parameters will be assigned for different values of the speed. Application of our method to 2D and 3D equations is left for future work.



---

*Source: 925920-2012-05-29.xml*