\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2015 (2015), No. 251, pp. 1--12.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2015/251\hfil Lattice Boltzmann method]
{Lattice Boltzmann method for coupled \\ Burgers equations}

\author[Y. Duan, L. Kong, X. Chen \hfil EJDE-2015/251\hfilneg]
{Yali Duan, Linghua Kong, Xianjin Chen}

\address{Yali Duan (corresponding author) \newline
School of Mathematical Sciences,
University of Science and Technology of China,
Hefei, Anhui 230026, China}
\email{ylduan01@ustc.edu.cn}

\address{Linghua Kong \newline
School of Mathematics and Information Science,
Jiangxi Normal University, Nanchang, Jiangxi, 330022, China}
\email{konglh@mail.ustc.edu.cn}

\address{Xianjin Chen \newline
School of Mathematical Sciences,
University of Science and Technology of China,
Hefei, Anhui 230026, China}
\email{chenxjin@ustc.edu.cn}

\thanks{Submitted May 1, 2015. Published September 25, 2015.}
\subjclass[2010]{65M12, 65M99}
\keywords{Lattice Boltzmann model; equilibrium distribution;
\hfill\break\indent coupled Burgers equations;
Chapman-Enskog expansion}

\begin{abstract}
 In this paper, we propose a lattice Boltzmann model for  coupled  
 Burgers equations (CBEs). With a proper time-space scale and the
 Chapman-Enskog expansion,  the  governing equations are recovered 
 successfully from the lattice Boltzmann  equations, and the resulting 
 local equilibrium distribution functions are  also obtained. 
 The partial derivative $\partial(uv)/\partial x$ in the model
 is treated as a source term  and discretized with a 2nd-order central 
 difference  scheme.  Numerical experiments show that the numerical results
 by the  Lattice Boltzmann Method  (LBM)  either agree well with the 
 corresponding exact solutions or are quite comparable  with those available
 numerical results in the literature.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{example}[theorem]{Example}
\allowdisplaybreaks

\section{Introduction}

 In this article, we are  concerned with the lattice Boltzmann
numerical solutions of the coupled  Burgers equations (CBEs) which
were first derived by Esipov {\cite{Esipov}} to study the model of
polydispersive sedimentation. The coupled Burgers equations can be
written as the  nonlinear partial differential equations
\begin{gather}
\frac{\partial u}{\partial t}+\eta u\frac{\partial u}
{\partial x}+\alpha\frac{\partial (uv)}{\partial x}
-\frac{\partial^2 u}{\partial x^2}=0 , \quad
x\in[a,b],\; t\in[0,T],\label{e1} \\
\frac{\partial v}{\partial t}+\xi v
\frac{\partial v}{\partial x}+\beta\frac{\partial (uv)}{\partial x}
-\frac{\partial^2 v}{\partial x^2}=0 , \quad
x\in[a,b],\; t\in[0,T]\label{e2}
\end{gather}
with initial conditions
\begin{equation}
u(x,0)=\phi(x),\quad v(x,0)=\varphi(x),   \label{e3}
\end{equation}
and boundary conditions
\begin{equation}
u(a,t)=h_1(a,t),\quad u(b,t)=h_2(b,t),\quad
v(a,t)=l_1(a,t),\quad v(b,t)=l_2(b,t),   \label{e4}
\end{equation}
where $\eta$, $\xi$ are real constants, $\alpha$, $\beta$ are
arbitrary constants depending on the system parameters such as the
{P\'{e}clet} number, the Stokes velocity of particles due to
gravity, and the Brownian diffusivity. This system  is a simple
model of sedimentation or evolution of scaled volume concentrations
of two kinds of particles in fluid suspensions or collision under the
effect of gravity.

Note that exact solutions of \eqref{e1}--\eqref{e4}
are  generally not available.
Numerically, \eqref{e1}--\eqref{e4} have been studied by
many authors with different methods, see \cite{Esipov}-\cite{Mittal}.
 Esipov \cite{Esipov} provided some numerical simulations for the CBEs
and compared them with available experimental results.  Kaya \cite{Kaya}
used the decomposition method to successfully obtain some solutions of both
homogenous and inhomogeneous CBEs in the form of a convergent power
series. In order to solve the CBEs, Abdou and Soliman \cite{Abdou} used a variational
iteration method, while Dehghan et al.\ \cite{Dehghan} developed and used the
Adomian-P\'{a}de technique. Khater et al.\ \cite{Khater} used
the Chebyshev spectral collocation method to solve the CBEs and successfully obtained
some approximate solutions.  Rashid et al.\ \cite{Rashid} proposed
a Fourier pseudo-spectral method. Based on a Crank-Nicolson
formulation for time integration and cubic B-spline functions for
space integration, Mittal et al.\ \cite{Mittal} also proposed a useful numerical
scheme for the CBEs.

Recently, the lattice Boltzmann method (LBM) has been developed
as  an alternative and promising numerical scheme for simulating
fluid flows \cite{Benzi,Qian, Chen,Luo00} and for solving various
mathematical physics problems
\cite{Dawson,Yan00,Zhang09,Duan,Shi,Shen1,Lai}. The LBM can be
either regarded as an extension of the lattice gas automaton
\cite{McNamara88} or as a special discrete form of the Boltzmann
equation for the kinetic theory {\cite{He971}}. Unlike conventional
numerical schemes which are on a basis of a discretization of partial differential
equations under macroscopic conservation laws, the LBM is based
on microscopic models and mesoscopic kinetic equations. Thanks to
its advantages in geometrical flexibility, natural parallelity,
simplicity of programming and numerical efficiency, the LBM has made
great success in many fields and been extensively applied to solve
numerous problems. It also shows some potentials to simulate
nonlinear systems, such as reaction-diffusion
equations \cite{Dawson}, convection-diffusion equations {\cite{Shi}},
wave equations \cite{Yan00} and Burgers equations
\cite{Duan,Shen1}. In this paper, a lattice Boltzmann model for the
CBEs is developed. We treat the partial derivative
$\partial(uv)/\partial x$  as a source term  and  then use the second order central
difference scheme to approximate such term. With
a proper time-space scale and the Chapman-Enskog
expansion, the resulting governing equations can be recovered  from the
lattice Boltzmann equations and the resulting local equilibrium distribution
functions can also be obtained. Numerical results are presented
and are compared with those available ones in the literature.
Numerical experiments show that the LBM method is reliable
and efficient for solving the CBEs.

This article is organized as follows. In Section 1, a brief introduction
is given. In Section 2, a lattice Boltzmann model for the CBEs is developed. The
application of the LBM to the coupled Burgers equations is also presented
in the section. In Section 3, three numerical examples are given to
illustrate the efficiency and accuracy of the LBM method
numerically. Some conclusions are made in the final section.

\section{Lattice Boltzmann model}

 To propose a lattice Boltzmann model for the
coupled Burgers equations, we rewrite \eqref{e1}--\eqref{e2} in the
 form
\begin{gather} %\label{CBE11}
\frac{\partial u}{\partial t}+\eta u\frac{\partial u}{\partial x}
-\frac{\partial^2 u}{\partial x^2}=-\alpha\frac{\partial (uv)}{\partial x} ,
 \label{e5} \\
%\label{CBE12}
\frac{\partial v}{\partial t}+\xi v\frac{\partial
v}{\partial x}-\frac{\partial^2 v}{\partial
x^2}=-\beta\frac{\partial (uv)}{\partial x} .  \label{e6}
\end{gather}

 According to the theory of the lattice Boltzmann method, it
consists of two steps:
(1) colliding, which occurs when particles arrive at an interaction
node from different directions and possibly change their moving
directions in the light of scattering rules;
(2) streaming, which occurs when each particle moves to the nearest node
 along its moving direction.
Usually, with the single-relaxation-time and the Bhatnagar-Gross-Krook
(BGK)  approximation \cite{BGK}, these two steps can be combined together
and converted into the lattice Boltzmann equations.
The corresponding evolution equations of the
distribution functions in the model read as
\begin{gather} %\label{BGD1}
 f_i(x+e_i\Delta t,t+\Delta t)-f_i(x,t)= -\frac{\Delta
t}{\tau}(f_i(x,t)-f^{eq}_i(x,t))+\Delta t\frac{F_i(x,t)}{b},
 \label{e7} \\
%\label{BGD2}
g_i(x+e_i\Delta t,t+\Delta t)-g_i(x,t)= -\frac{\Delta
t}{\tau}(g_i(x,t)-g^{eq}_i(x,t))+\Delta t\frac{G_i(x,t)}{b},
 \label{e8}
\end{gather}
where $f_i(x,t)$ and $ g_i(x,t)$ are the distribution functions of
particles; $f_i^{eq}(x,t)$ and $g_i^{eq}(x,t)$  are the local
equilibrium distribution functions of particles; $\Delta x$ and
$\Delta t$ are the lattice space and time increments, respectively;
$c=\Delta x/\Delta t$ is the particle speed and $\tau$ is the
dimensionless relaxation time; $F_i(x,t)$ and $ G_i(x,t)$ depend
only on the
 right hand side of \eqref{e5}--\eqref{e6}, respectively;
 $\{ e_i, i=0,1,\dots,b-1 \}$ is the set of discrete velocity
directions, for the 3-bit model in the one-dimensional case,
$\{e_0, e_1, e_2 \}=\{ 0,c,-c\}$.  The macroscopic variables, $u$ and
$v$ are defined in terms of the distribution functions as
\begin{equation} %\label{u}
u=\sum_i f_i=\sum_i f^{eq}_i,~~~~v=\sum_i g_i=\sum_i g^{eq}_i. \label{e9}
\end{equation}

 To derive the macroscopic equations \eqref{e5}--\eqref{e6}
(or \eqref{e1}--\eqref{e2}) from the lattice BGK
model, the Chapman-Enskog expansion is performed under the
assumption that the Kundsen number
$\epsilon=l/L$ is small, where $l$ is the mean free path and $L$ is the
characteristic length. The Chapman-Enskog expansion is then applied to
$f_i(x,t)$ and $g_i(x,t)$, i.e.,
\begin{gather} % \label{f}
f_i=f_i^{eq}+\sum_{n=1}^\infty \epsilon^nf_i^{(n)}
= f^{eq}_i+\epsilon f^{(1)}_i+\epsilon^2f^{(2)}_i+\dots \label{e10}\\
% \label{g}
g_i=g_i^{eq}+\sum_{n=1}^\infty \epsilon^ng_i^{(n)}
= g^{eq}_i+\epsilon g^{(1)}_i+\epsilon^2g^{(2)}_i+\dots \label{e11}
\end{gather}

With time scale $t_1=\epsilon t, t_2=\epsilon^2 t$ and space
scale $ x_1=\epsilon x$, then the time derivation and the space
derivation can be expanded formally as:
\begin{equation} %\label{derivation}
\partial_t=\epsilon\partial_{t_1}+\epsilon^2\partial_{t_2},\quad \partial_
x=\epsilon \partial_{x_1}.\label{e12}
\end{equation}

In \eqref{e7}--\eqref{e8},  assume further that the second order terms
$F_i(x,t)$ and $G_i(x,t)$ take the form
\begin{equation}
F_i(x,t)=\epsilon^2 \tilde{F}_i(x,t),\quad
G_i(x,t)=\epsilon^2 \tilde{G}_i(x,t)\label{e13}
\end{equation}
for some functions $\tilde{F}_i(x,t),\tilde{G}_i(x,t)$.

For simplicity, we only consider how to recover \eqref{e5} from \eqref{e7},
and \eqref{e6} from \eqref{e8} through a similar derivation.
Performing the Taylor expansion on \eqref{e7} up to terms with order
 of $O(\epsilon^3)$
and using \eqref{e10}, \eqref{e12}--\eqref{e13},
one can get the partial differential equations in order of $\epsilon$,
$\epsilon^2$ as
\begin{gather}
(\partial_{t_1}+ e_i \partial_{x_1})f_i^{eq}
=-\frac{1}{\tau \Delta t}f_i^{(1)}, \label{e14} \\
\partial_{t_2}f_i^{eq}+(\partial_{t_1}+e_i\partial_{x_1})f_i^{(1)}+\frac{\Delta
t}{2}(\partial_{t_1}+e_i\partial_{x_1})^2f_i^{eq}=-\frac{1}{\tau
\Delta t}f_i^{(2)}+\frac{\tilde{F}_i( {\bf x},t)}{b}. \label{e15}
\end{gather}

Taking \eqref{e14} times $\epsilon$ plus \eqref{e15} times $\epsilon^2 $,
summing over $i$ and using \eqref{e9}, \eqref{e12}, \eqref{e14} and
\begin{equation}
 \sum_i f^{(k)}_i=0\quad (k\geq 1), \label{e16}
\end{equation}
we have
\begin{equation}
 \partial_{t}u+\partial_x(\sum_{i=0}^{b-1}e_if_i^{eq})+
 \epsilon^2\Delta t(\frac 12-\tau)\sum_{i=0}^{b-1}(\partial_{t_1}+
e_i\partial_{x_1})^2f_i^{eq}
= \sum_{i=0}^{b-1}\epsilon^2\frac{\tilde{F}_i(x,t)}{b}+O(\epsilon^3).
\label{e17}
\end{equation}
Now, in view of \eqref{e5}, one can set
\begin{gather}
\sum_{i=0}^{b-1} e_if_i^{eq}=\frac{\eta}{2}u^{2}, \label{e18} \\
\epsilon^2\sum_{i=0}^{b-1}(\partial_{t_1}+
e_i\partial_{ x_1})^2f_i^{eq}=\lambda\frac{\partial^2u}{\partial
x^2}, \label{e19} \\
\sum_{i=0}^{b-1}\epsilon^2\frac{\tilde{F}_i(x,t)}{b}=\sum_{i=0}^{b-1}\frac{F_i(
x,t)}{b}=-\alpha\frac{\partial (uv)}{\partial x}, \label{e20}
\end{gather}
where
\begin{equation}
\lambda=\frac{1}{\Delta t (\tau-0.5)}.\label{e21}
\end{equation}
The CBEs with the second
order accuracy of truncation error is hence obtained.

In addition, if summing \eqref{e14} over $i$ and using \eqref{e9}, \eqref{e16}
 and \eqref{e18}, we obtain
\begin{equation}
\partial_{t_1}u+\eta u\partial_{x_1}u=0.\label{e22}
\end{equation}
On the other hand, using \eqref{e14}, \eqref{e19} and \eqref{e22}, we  get the
2nd-order moment equation of the resulting local equilibrium distribution function
\begin{equation}
\sum_{i=0}^{b-1} e_{i} e_{i}f_i^{eq}
=\lambda u+\frac{\eta^2}{3}u^{3}. \label{e23}
\end{equation}
In view of \eqref{e9}, \eqref{e18} and \eqref{e23},  the
equilibrium distribution functions $\{f^{eq}_i\}$ of 3-bit model can be
obtained as
\begin{equation}
\begin{gathered}
f^{eq}_0=u-\frac{\eta^2}{3c^2}u^{3}-\frac{\lambda u}{c^2} ,  \\
f^{eq}_1= \frac{\lambda
u}{2c^2}+\frac{\eta^2}{6c^2}u^{3}+\frac{\eta}{4c}u^{2},  \\
f^{eq}_2= \frac{\lambda
u}{2c^2}+\frac{\eta^2}{6c^2}u^{3}-\frac{\eta}{4c}u^{2}.
\end{gathered}
\label{e24}
\end{equation}
Similarly, the
equilibrium distribution functions for $\{g^{eq}_i\}$ can be achieved as
\begin{equation}
\begin{gathered}
g^{eq}_0=v-\frac{\xi^2}{3c^2}v^{3}-\frac{\lambda v}{c^2} ,  \\
g^{eq}_1= \frac{\lambda v}{2c^2}+\frac{\xi^2}{6c^2}v^{3}+\frac{\xi}{4c}v^{2},  \\
g^{eq}_2= \frac{\lambda v}{2c^2}+\frac{\xi^2}{6c^2}v^{3}-\frac{\xi}{4c}v^{2}.
\end{gathered} \label{e25}
\end{equation}

Next, we use the second order central difference  scheme to approximate
partial derivatives $ \partial (uv)/\partial x$ in \eqref{e5}--\eqref{e6}, i.e.,
the source terms
\begin{gather}
\begin{aligned}
F(x,t)&=-\alpha\frac{\partial (uv)}{\partial x}\\
&=-\alpha u\frac{\partial v}{\partial x}-\alpha v\frac{\partial u}{\partial x}\\
&\approx -\alpha\Big[u \frac{v(x+\Delta x,t)-v(x-\Delta x,t)}{2\Delta x}\\
&\quad +v \frac{u(x+\Delta x,t)-u(x-\Delta x,t)}{2\Delta x})\Big],
\end{aligned}\label{e26}
\\
\begin{aligned}
G(x,t)&=-\beta\frac{\partial (uv)}{\partial x}\\
&=-\beta u\frac{\partial v}{\partial x}-\beta v\frac{\partial u}{\partial x}\\
&\approx -\beta\Big[u \frac{v(x+\Delta x,t)-v(x-\Delta x,t)}{2\Delta x}\\
&\quad +v \frac{u(x+\Delta x,t)-u(x-\Delta x,t)}{2\Delta x})\Big].
\end{aligned}\label{e27}
\end{gather}

Finally, in order not to make large errors on the
boundary which can influence the global evolution, we use the
three-adjacent-point difference scheme (second order accuracy) to approximate partial
derivatives $ v\frac{\partial u}{\partial x}$,
$ u\frac{\partial v}{\partial x}$ on the boundary, i.e.,
\begin{gather}
v\frac{\partial u}{\partial x}|_{x=a}
=v(a,t)\frac{-3u(a,t)+4u(a+\Delta x,t)-u(a+2\Delta x,t)}{2\Delta x},\label{e28}
\\
u\frac{\partial v}{\partial x}|_{x=a}
=u(a,t)\frac{-3v(a,t)+4v(a+\Delta x,t)-v(a+2\Delta x,t)}{2\Delta x},\label{e29}
\\
v\frac{\partial u}{\partial x}|_{x=b}=v(b,t)\frac{u(b-2\Delta
x,t)-4u(b-\Delta x,t)+3u(b,t)}{2\Delta x},\label{e30}
\\
u\frac{\partial v}{\partial x}|_{x=b}=u(b,t)\frac{v(b-2\Delta
x,t)-4v(b-\Delta x,t)+3v(b,t)}{2\Delta x}.\label{e31}
\end{gather}

\section{Numerical results and analysis}

In this section,  numerical results of the 3-bit LBM
 for the coupled Burgers equations are given.
 The efficiency and accuracy of the LBM are demonstrated through three numerical
examples. In particular, we compare the LBM numerical simulations  with the exact
solutions and/or some available numerical results
\cite{Khater,Mittal} in the literature.
Here, we use the discrete $L_{\infty}$-norm
error and the relative discrete $L_2$-norm error as follows:
\begin{gather*}
\|E(u)\|_\infty=\max_{0\leq i\leq N}\{|u^{\rm exact}(x_i,t)-u^{\rm num}(x_i,t)|\},
\\
\|E(u)\|_2= \Big(\sum_{i=0}^N|u^{\rm exact}(x_i,t)-u^{\rm num}(x_i,t)|^2\Big)^{1/2}
\Big/\Big(\sum_{i=0}^N|u^{\rm exact}(x_i,t)|^2\Big)^{1/2},
\end{gather*}
where  $N$ is the number of the lattice, $u^{\rm num}(x_i,t)$ and
$u^{\rm exact}(x_i,t)$ represent the
numerical solution and the exact solution, respectively.

\begin{example} \label{examp1} \rm
To examine the performance of the 3-bit LBM for
solving Burgers equations, we set the parameters $\eta=\xi=-2$ and
$ \alpha=\beta=1$.
With initial condition %given by
$$
u(x,0)=v(x,0)=\sin(x),
$$
Equations \eqref{e1}--\eqref{e2} read as follows
\begin{gather}
\frac{\partial u}{\partial t}-2u\frac{\partial u}{\partial x}
+ \frac{\partial (uv)}{\partial x}-\frac{\partial^2 u}{\partial x^2}=0 ,\label{e32}
\\
\frac{\partial v}{\partial t}-2v\frac{\partial v}{\partial
x}+\frac{\partial (uv)}{\partial x}-\frac{\partial^2 v}{\partial
x^2}=0 ,\label{e33}
\end{gather}
whose exact solution is \cite{Kaya}
$$
u(x,t)=v(x,t)=e^{-t}\sin(x).
$$

 To obtain the LBM numerical solution, we choose $\Delta x=0.04$, $\Delta t=0.001$,
$\lambda=c^2/2$, and the domain $\{x, x\in [-\pi,\pi]\}$.
 As before, $c=\Delta x/\Delta t$ is the particle speed.
To test the computational efficiency and accuracy
of the  LBM, we use both $L_2$-norm and
$L_\infty$-norm errors for $u$, $v $ %at  different time stages
as shown in Table \ref{table1}.
The numerical solution $u(x, t)$ is depicted in Figure \ref{fig1} for
 $t=0.1, 0.5, 1$. Because of  the symmetry of $u$ and $v$,
 it is also true for $v(x, t)$.
\end{example}

 \begin{table}[htb]
\caption{$L_\infty$, $L_2$ errors for $u$, $v $ at different time stages}
\label{table1}
\begin{center}
\renewcommand{\arraystretch}{1.3}
 \begin{tabular}{|c|cc|cc|} \hline
&   \multicolumn{2}{c|}{$u$}    & \multicolumn{2}{c|}{$v$}   \\ \cline{2-5}
 $t$ &  $L_\infty$ &   $L_2 $ &  $L_\infty$ &  $L_2 $ \\ \hline
 0.1 &     8.6936$e$--4 & 6.1750$e$--4& 8.6936$e$--4 & 6.1750$e$--4\\
 0.5 &     5.4963$e$--4& 3.1121$e$--4 & 5.4963$e$--4& 3.1121$e$--4 \\
  1&    2.3350$e$--4 & 1.3417$e$--4& 2.3350$e$--4 & 1.3417$e$--4  \\\hline
\end{tabular}
 \end{center}
\end{table}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig1} % CBE1.eps
\end{center}
\caption{A comparison of numerical and
exact solutions in Example \ref{examp1}}
\label{fig1}
\end{figure}

\begin{example} \label{examp2} \rm
Choosing $\eta=\xi=2$, the CBEs in \eqref{e1}--\eqref{e2} are
\begin{gather}
\frac{\partial u}{\partial t}+2u\frac{\partial u}{\partial x}
+ \alpha\frac{\partial (uv)}{\partial x}-\frac{\partial^2 u}{\partial x^2}=0 ,
\label{e34} \\
\frac{\partial v}{\partial t}+2v\frac{\partial v}{\partial
x}+\beta\frac{\partial (uv)}{\partial x}-\frac{\partial^2
v}{\partial x^2}=0.\label{e35}
\end{gather}
If we use the initial conditions
$$
u(x,0)=a_0(1-\tanh (Dx)),\quad
v(x,0)=a_0[\frac{2\beta-1}{2\alpha-1}-\tanh(Dx)],
$$
then the corresponding exact solution $(u,v)$ for \eqref{e34}--\eqref{e35}
 is \cite{Soliman}
\begin{gather*}
u(x,t)=a_0(1-\tanh (D(x-2Dt))),\\
v(x,t)=a_0[\frac{2\beta-1}{2\alpha-1}-\tanh(D(x-2Dt))],
\end{gather*}
where $D=a_0(4\alpha\beta-1)/(4\alpha-2)$, $a_0$ is an arbitrary
constant. The LBM numerical results of this example are showed in
Tables \ref{table2}--\ref{table3} by using the domain
$[-10,10]$, the $L_{\infty}$ error and the relative $L_2$ error with
 $\Delta t=0.01$, $\Delta x=1$,
$c=\Delta x/\Delta t$, $\lambda=c^2/3$.
In both tables, we compare the LBM
results with those available ones in the literature
\cite{Khater,Mittal}.
As is seen from the tables, the 3-bit lattice
Boltzmann method is quite comparable with those methods used in
\cite{Khater,Mittal}.
\end{example}

\begin{table}[htb]
\caption{Comparison of $L_\infty$, $L_2$ errors for $u(x,t)$ in Example \ref{examp2}.}
\label{table2}
\begin{center} 
\renewcommand{\arraystretch}{1.3} 
\small
\begin{tabular}{|ccc|cc|cc|cc|} \hline
& & &  \multicolumn{2}{c|}{Khater \cite{Khater}}   
& \multicolumn{2}{c|}{Mittal \cite{Mittal}}
& \multicolumn{2}{c|}{LBM}  \\ \cline{4-9}
 $t$ &$\alpha$ & $\beta$
&  $\|E(u)\|_\infty$ &   $\|E(u)\|_2 $ &  $\|E(u)\|_\infty$ &  $\|E(u)\|_2 $
&  $\|E(u)\|_\infty $&  $ \|E(u)\|_2 $  \\ \hline
  0.5   & 0.1   &     0.3   &   4.38$e$--5          & 1.44$e$--3
& 4.167$e$--5     & 6.736$e$--4        &  4.028$e$--5 & 6.589$e$--4\\
 &0.3& 0.03  &  4.58$e$--5 & 6.68$e$--4& 4.590$e$--5
&7.326$e$--4& 4.936$e$--5 & 7.335$e$--4  \\ \hline\hline
1      & 0.1        &  0.3  &  8.66$e$--5 & 1.27$e$--3 & 8.258$e$--5
& 1.325e--3& 8.140$e$--5& 1.298$e$--3  \\
  &0.3 &  0.03   & 9.16$e$--5  & 1.30$e$--3& 9.182$e$--5
&1.452$e$--3& 9.014$e$--5& 1.422$e$--3  \\  \hline
\end{tabular}
 \end{center}
\end{table}

\begin{table}[htb]
\caption{Comparison of $L_\infty$, $L_2$ errors for $v(x,t)$ in Example \ref{examp2}.}
\label{table3}
\begin{center} 
\renewcommand{\arraystretch}{1.3}
\small
 \begin{tabular}{|ccc|cc|cc|cc|} \hline
& & & \multicolumn{2}{c|}{Khater \cite{Khater}}  
& \multicolumn{2}{c|}{Mittal \cite{Mittal}}
& \multicolumn{2}{c|}{LBM}  \\ \cline{4-9}
 $t$ &$\alpha$ &$\beta$ 
&  $\|E(v)\|_\infty$ &   $\|E(v)\|_2 $ &  $\|E(v)\|_\infty$
&  $\|E(v)\|_2 $&  $\|E(v)\|_\infty $&  $ \|E(v)\|_2 $  \\ \hline
  0.5   & 0.1   &     0.3   &   4.99$e$--5         & 5.42$e$--4
& 1.480$e$--4     & 9.057$e$--4        &  2.021$e$--5 & 4.705$e$--4\\
 &0.3& 0.03  &  1.81$e$--4 & 1.20$e$--3
&    5.729$e$--4  &1.591$e$--3       & 1.862$e$--4 & 1.301$e$--3  \\ \hline\hline
1      & 0.1        &  0.3  &  9.92$e$--5 & 1.29$e$--3    & 4.770$e$--5
& 1.252e--3            & 4.029$e$--5& 9.389$e$--4  \\
  &0.3 &  0.03   & 3.62$e$--4  & 2.35$e$--3
& 3.617$e$--4 &2.250$e$--3     & 3.657$e$--4& 2.557$e$--3  \\  \hline
\end{tabular}
 \end{center}
\end{table}

\begin{example} \label{examp3} \rm
Consider the CBEs \eqref{e1}--\eqref{e2} with initial condition \cite{Mittal}
\begin{gather*}
u(x,0)= \begin{cases}
\sin(2\pi x) ,  & 0\leq x\leq 0.5, \\
0,  & 0.5< x\leq 1,
\end{cases}
\\
v(x,0)= \begin{cases}
 0,  & 0\leq x\leq 0.5, \\
-\sin(2\pi x), & 0.5< x\leq 1,
\end{cases}
\end{gather*}
and zero boundary conditions. The numerical solution is computed on the domain
$[0,1]$. In Tables \ref{table4} and \ref{table5}, we present a comparison of the
3-bit LBM numerical results  and those obtained by the cubic B-spline
collocation scheme \cite{Mittal}  in term of
maximum values of $u, v$ for different
times $(i.e., t=0.1, 0.2, 0.3, 0.4)$ and different values of $\alpha, \beta$.
Here, we choose $\eta=\xi=2$,
$N=50, \Delta t=0.0001$ and $\lambda=c^2/2$.
The numerical solutions $(u, v)$  are also depicted in
Figures \ref{fig2}--\ref{fig3} with
$\alpha=\beta=10$ and $\alpha=\beta=100$, respectively.
A sharp decay is clearly displayed in Figure \ref{fig3} when
$\alpha, \beta$ are large. Note that the 3-bit LBM numerical results also
agree well with those in \cite{Mittal}.
\end{example}

\begin{table}[htb]
\caption{Comparison of maximum values of $u(x,t), v(x,t)$ for $\alpha =\beta $ = 10.}
\label{table4}
\begin{center}
\renewcommand{\arraystretch}{1.3}
 \begin{tabular}{|c|cc|cc||cc|cc|} \hline
& \multicolumn{2}{c|}{Mittal \cite{Mittal}} 
& \multicolumn{2}{c|}{LBM} 
& \multicolumn{2}{c|}{Mittal \cite{Mittal}} 
& \multicolumn{2}{c|}{LBM}  \\ \cline{2-9}
$t$\hspace{0.1cm} &  $u^{\rm max}$&  $x$   &  $u^{\rm max}$ &  $x$ &
$v^{\rm max}$ & $x$ & $v^{\rm max}$& $x$
\\\hline
 0.1 & 0.14456   & 0.58  &  0.14464 & 0.58  & 0.14306  & 0.66  &  0.14327 & 0.66\\
0.2 & 0.05237   & 0.54  & 0.05241  & 0.54  & 0.04697  & 0.56  & 0.04704  & 0.56  \\
0.3 & 0.01932   & 0.52  & 0.01934  & 0.52  & 0.01725  & 0.52  & 0.01727 & 0.52  \\
0.4 & 0.00718   & 0.50  & 0.00719  & 0.50  &0.00641   & 0.50  & 0.00642 & 0.50  \\
\hline
\end{tabular}
 \end{center}
\end{table}

\begin{table}[htb]
\caption{Comparison of maximum values of $u(x,t), v(x,t)$ for $\alpha =\beta= 100 $.}
\label{table5}
\begin{center}
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|c|cc|cc||cc|cc|} \hline
& \multicolumn{2}{c|}{Mittal \cite{Mittal}} 
& \multicolumn{2}{c|}{LBM} 
& \multicolumn{2}{c|}{Mittal \cite{Mittal}} 
& \multicolumn{2}{c|}{LBM}  \\ \cline{2-9}
$t$\hspace{0.1cm} & $u^{\rm max}$&  $x$   &   $u^{\rm max}$ &  $x$ &
$v^{\rm max}$ & $x$&$v^{\rm max}$& $x$
\\\hline
 0.1 & 0.04175   & 0.46  &  0.04167 & 0.46  & 0.05065  & 0.76  &  0.05080 & 0.76\\
0.2 & 0.01479   & 0.58  & 0.01478  & 0.58  & 0.01033  & 0.64  & 0.001035  & 0.64  \\
0.3 & 0.00534   & 0.54  & 0.00534  & 0.54  & 0.00350  & 0.56  & 0.00351 & 0.56  \\
0.4 & 0.00198   & 0.52  & 0.00198  & 0.52  &0.00129   & 0.52  & 0.00129 & 0.52  \\  \hline
\end{tabular}
 \end{center}
\end{table}

\begin{figure}[htb]
 \begin{center}
\includegraphics[width=0.5\textwidth]{fig2a}\\ % CBE10u.eps
\includegraphics[width=0.5\textwidth]{fig2b} % CBE10v.eps
\end{center}
\caption{Numerical solution $(u(x,t), v(x,t))$ of Example \ref{examp3} for $\alpha=\beta=10$.}
\label{fig2}
\end{figure}

\begin{figure}[htb]
 \begin{center}
\includegraphics[width=0.5\textwidth]{fig3a} \\ % CBE100u.eps
\includegraphics[width=0.5\textwidth]{fig3b} % CBE100v.eps 
\end{center}
\caption{Numerical solution $(u(x,t), v(x,t))$ of Example \ref{examp3} for $\alpha=\beta=100$.}
\label{fig3}
\end{figure}


To the best of our knowledge, there is no exact solution available in the
literature for general values of $\eta, \xi$. Hence, numerical
experiments are important in study of the behavior of solutions for
larger values of $\eta, \xi$.
With $\alpha=\beta=10$, numerical solutions $(u, v)$ in Example \ref{examp3}
 are plotted  in Figures \ref{fig4}--\ref{fig5} to
illustrate the effect of different values of $\eta, \xi$, i.e.,
$\eta=\xi=20, 200$.
From the figures it can be seen that the solution $(u, v)$ decays
to zero as time $t$ or the value of
$\eta, \xi$ increases.
 It can also be observed that the  LBM is capable
of finding numerical solutions for larger values of $\eta, \xi$.
As a comparison, if
taking the numerical solution on a fine mesh of 500
lattices as the ``exact" solution,
then one can see from Figures \ref{fig4}--\ref{fig5} that the numerical solution
on a much coarser mesh (containing 50 lattices) agrees well with such
``exact" solution.

\begin{figure}[htb]
 \begin{center}
\includegraphics[width=0.48\textwidth]{fig4a}  % 20u.eps 
\includegraphics[width=0.48\textwidth]{fig4b} % 20v.eps 
\end{center}
\caption{Numerical solution $(u, v)$ at different time
levels for $\alpha=\beta=10$ and $\eta=\xi=20$.}
\label{fig4}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.48\textwidth]{fig5a} % 200u.eps
\includegraphics[width=0.48\textwidth]{fig5b} % 200v.eps
\end{center}
\caption{Numerical solution $(u, v)$ at different time
levels for $\alpha=\beta=10$ and $\eta=\xi=200$.}
\label{fig5}
\end{figure}


\subsection*{Conclusions}

   In this article, a 3-bit lattice Boltzmann model for the
CBEs is proposed. The partial derivative $\partial(uv)/\partial x$
is treated as a source term  and discretized with a second-order  central
difference scheme. By choosing an appropriate time-space
scale and applying the Chapman-Enskog expansion, the resulting governing
equations are recovered from the lattice Boltzmann
equations and the resulting local equilibrium distribution functions are also
obtained.
Then, three numerical examples are presented to illustrate the efficiency of the LBM.
Numerical experiments show that LBM results agree well with either
the exact solutions or  other available results.
Certain important behaviors of general CBEs are also verified numerically.
Undoubtedly, it is worth further investigating
the LBM in solving other types of nonlinear partial differential equations in future.


\subsection*{Acknowledgements}
The authors are grateful to anonymous referees for
their valuable comments and suggestions. This work was supported by
the NSF of China (Grant No. 11101399, 10901074) and
 the Fundamental Research Funds for the Central Universities (No. WK0010000022),
 and partially by Anhui Provincial NSF (No. 11040606M05).

\begin{thebibliography}{00}

\bibitem{Abdou} M. A. Abdou, A Soliman;
A Variational iteration method for solving
Burgers' and coupled Burgers' equations, \textit{J. Comp. Appl.
Math.} 181 (2005), 245-251.

\bibitem{Benzi} R. Benzi, S. Succi, M. Vergassola;
The lattice Boltzmann equation: theory and applications, \textit{Phys. Report}
  222(3) (1992), 145-197.

\bibitem{BGK}  P. Bhatnagar, E. Gross, M. Krook;
 A model for collision process in gas. I: Small amplitude processed in
charged and neutral one component system, \textit{Phys. Rev.} 94 (1954) 511-525.

\bibitem{Chen}  S. Chen, G. D. Doolen;
 Lattice Boltzmann method for fluid flows,
\textit{Annu. Rev. Fluid Mech.} 30 (1998), 329-364.

\bibitem{Dehghan} M. Dehghan, A. Hamidi, M. Shakourifar;
The solution of coupled Burgers' equations using Adomian-Pade technique,
\textit{Appl. Math. Comput.} 189 (2007), 1034-1047.

 \bibitem{Duan} Y. L. Duan, R. X. Liu;
Lattice Boltzmann model for two-dimensional unsteady Burgers' equation,
\textit{J. Comput. Appl. Math.} 206 (2007) 432-439.

\bibitem{Esipov} S. E. Esipov;
 Coupled Burgers' equations: a model of polydispersive
sedimentation, \textit{Phys Rev. E} 52 (1995), 3711-3718.

\bibitem{He971} X. Y. He, L. S. Lou;
 Theory of the lattice Boltzmann method: From
the Boltzmann equation to the lattice Boltzmann equation,
\textit{Phy. Rev. E} 56 (1997), 6811-6817.

\bibitem{Kaya} D. Kaya;
 An explicit solution of coupled viscous Burgers' equation
by the decomposition method, \textit{IJMMS } 27(11) (2001), 675-680.

\bibitem{Khater} A. H. Khater, R. S. Temsah, M. M. Hassan;
 A Chebyshev spectral collocation method for solving Burgers-type equations,
 \textit{J. Comput. Appl. Math.} 222(2) (2008),  333-350.

\bibitem{Lai} H. L. Lai, C. F. Ma;
 Lattice Boltzmann method for the generalized Kuramoto-Sivashinsky equation,
 \textit{Physica A} 388 (2009) 1405-1412.

\bibitem{Luo00}  L. Luo;
 The lattice-gas and lattice Boltzmann methods: past,present and future.
  Proceedings of International Conference on Applied Computational
Fluid Dynamics. Beijing,   China, October 17-20, Beijing, (2000), 52-83.

\bibitem{McNamara88} G. R. McNamara, G. Zanetti;
 Use of the Boltzmann equation to simulate lattice-gas automata,
\textit{Phys. Rev. Lett.} 61 (1988), 2332-2335.

\bibitem{Mittal} R. C. Mittal, Geeta Arora;
 Numerical solution of the coupled viscous
Burgers' equation, \textit{Commun Nonlinear Sci. Numer. Simulat.} 16
(2011), 1304-1313.

\bibitem{Dawson} S. Ponce Dawson, S. Chen, G. D. Doolen;
 Lattice Boltzmann computations for reaction-diffusion equations, \textit{J. Chem.
Phys.} 2 (1993), 1514-1523.

\bibitem{Qian} Y. Qian, S. Succi, S. Orszag;
Recent advances in lattice Boltzmann computing,
\textit{Annu. Rev. Comput. Phys.}  3 (1995), 195-242.

\bibitem{Rashid} A. Rashid, Ismail AIB;
 A Fourier Pseudospectral method for solving
coupled viscous Burgers equations, \textit{Comp. Methods. Appl.
Math.} 9(4) (2009), 412-420.

\bibitem{Shen1} Z. J. Shen, G. W. Yuan, L. J. Shen;
 Lattice Boltzmann method for Burgers equation,
\textit{Chinese J. Comp. Phys.}  17 (2000), 166-172.

 \bibitem{Shi} B. C. Shi, Z. L. Guo;
Lattice Boltzmann model for nonlinear convection-diffusion equations,
\textit{Phy. Rev. E}  79 (2009), 016701.

\bibitem{Soliman} A. A. Soliman;
 The modified extended tanh-function method for solving Burgers-type equations,
\textit{Physica A}, 361 (2006), 394-404.

\bibitem{Yan00} G. W. Yan;
 A lattice Boltzmann equation for waves,\textit{J. Comp. Phys.} 161 (2000)
61-69.

\bibitem{Zhang09} J. Y. Zhang, G. W. Yan;
A lattice Boltzmann model for the Korteweg-de Vries equation with
two conservation laws, \textit{Comp. Phys. Commu.} 180 (2009),
1054-1062.

\end{thebibliography}


\end{document}

%%
%% End of file `elsarticle-template-1-num.tex'.
