An Invariant-Preserving Scheme for the Viscous Burgers-Poisson System

Abstract

We formulate and analyze a new finite difference scheme for a shallow water model in the form of viscous Burgers-Poisson system with periodic boundary conditions. The proposed scheme belongs to a family of three-level linearized finite difference methods. It is proved to preserve both momentum and energy in the discrete sense. In addition, we proved that the method converges uniformly and has second order of accuracy in space. The analysis given in this work is the first time a pointwise error estimation is done on a second-order finite difference operator applied to the Burgers-Poisson system. We validate our findings by performing various numerical simulations on both viscous and inviscous problems. These numerical examples show the efficacy of the proposed method and confirm the proven theoretical results.

1. Introduction

The study of water waves has been one of the fascinating subjects in mathematical modeling as one can gain insights into ranges of natural phenomena. One of the earliest well-known water models is the Korteweg-De Vries (KdV) equation [1]. Its dimensionless form is given by

$$\begin{array}{c}\hfill u_t+u_{xxx}-6u{u}_x=0,\end{array}$$

(1)

where $t,x$ are temporal and spatial variables, respectively. The variable u represents the fluid velocity along the x-axis. The subscripts indicate the partial derivatives. Thanks to its completely integrable solution, the KdV model has been a subject of active research since its discovery. As a water wave model, the KdV equation has some shortcomings, one of which is its incapability to capture the wave-breaking phenomenon. In [2], Whitham proposed a model of the form

$$\begin{array}{cc}& {\partial }_{t}u+u{\partial }_{x}u={\partial }_x[G\ast u],\end{array}$$

(2)

where $G\ast u$ defines a convolution

$$\begin{array}{cc}\hfill [G\ast u]\left(x\right)& ={\int }_{\mathbb{R}}G(x-y)u\left(y\right)dy.\hfill \end{array}$$

(3)

For what is now known as the Whitham equation, the convolution kernel $G=G_{Wh}$ is given by

$$\begin{array}{cc}\hfill G_{Wh}\left(x\right)& =\frac{1}{2\pi }{\int }_{\mathbb{R}}{\left(\frac{g}{\kappa }tanh\left(\kappa h_0\right)e^{i\kappa x}\right)}^{1/2}dy,\hfill \end{array}$$

(4)

where g is the gravitational constant of acceleration, $h_0$ is the undisturbed water depth, assuming flat ground, and $\kappa$ is the wave number. When the ratio $\frac{g}{\kappa }$ is normalized, the kernel $G_{Wh}$ can be approximated by

$$\begin{array}{cc}\hfill G_{BP}\left(x\right)& =\frac{\pi }{4}{e}^{-\pi \left|x\right|/2}.\hfill \end{array}$$

(5)

The resulting model is called Burgers-Poisson equations [3]

$$\begin{array}{cc}\hfill & u_t+u{u}_x={\varphi }_x,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & {\varphi }_{xx}-\varphi =u.\hfill \end{array}$$

(7)

Unlike the KdV equation, the models derived from (2) using the kernels (4) and (5) capture the peaking and breaking of water waves. To understand the behavior of water waves on the shallow flat bottom better, we are interested in numerical approximation of the system (6) and (7) with periodic boundary condition. Not only will this help us gain more insights into the shallow water waves, but also many other natural phenomena with models of similar forms; for instance, the Two-Species-Euler-Poisson system as mentioned in [3]. We also add the viscous term to the Burgers equation in order to cover the one-dimensional version of the Navier-Stokes-Poisson system as mentioned in [4].

Analytic properties of (6) and (7) are first studied in [3] where the traveling solution is founded. Local existence for the smooth solution and global existence for the weak entropy solution are also confirmed. In [5], classification of group invariant solutions for the system (6) and (7) is determined by using the classical Lie method. In [6], another global existence for the entropy solution is studied under some regularity conditions on the initial data. In an alternate approach to study the system (6) and (7), one can also rewrite it into a single-equation form

$$\begin{array}{cc}\hfill & u_t-u_{xxt}+u_x+u{u}_x=3u_x{u}_{xx}+u{u}_{xxx},\hfill \end{array}$$

(8)

which has terms similar to other wave models such as Camassa-Holm equation [7,8], Benjamin-Ono equation [9,10], KdV equation, Rosenua equation [11,12], and the RLW equation [13,14], to name a few.

One of the common properties for these models is invariant-preserving: certain quantities derived from the solution stay constant or do not increase over time. The solution of (6) and (7) possesses two invariants, namely momentum and energy. This raises challenge when solving the model numerically as one needs to maintain the invariant-preserving property while achieving convergence. There have not been many attempts to solve the Burgers-Poisson model or its variations numerically, and only a few of them meet both of the requirements. In [15], the numerical solution of the Burgers-Poisson Equation (8) is studied by using both finite difference and variational iteration methods. The stability and consistency condition is provided, but the invariant-preserving property and convergence analysis are not explored. In [16,17], both exact and approximate explicit solutions of fractional-order Burgers-Poisson equation are founded by using homotopy perturbation method and also the Adomian decomposition method in the latter work. Numerical solutions and exact solutions are compared to show the efficiency of the presented methods, but convergence and stability analysis is not provided. In [4,18], the Local Discontinuous Galerkin (LDG) method is applied to the inviscous problem and the viscous problem, respectively. The analysis on invariant-preserving and convergence is included in these LDG works. However, the LDG method requires solving a large system of unknowns. Multiple steps are involved when employing the scheme and when interpreting the results. We seek a simpler approach that can be easily applied to the model while still satisfying the invariant-preserving and convergence conditions.

Meeting both requirements is not straightforward, but it is even more challenging to provide rigorous proof on analysis part. In this work, we are able to do both by presenting a collection of FDM schemes solving the system (6) and (7) along with analyses on the convergence and invariant preserving. This is the first instance of uniform error analysis on an FDM operator applied to the system. The operator consists of two finite differences using information from three time levels. Since the model is nonlinear, the coefficient matrices need to be redefined at each time step, but the linearization of the method eliminates the need to use iterations to solve for the unknowns in each time step. The $\theta$ parameter gives rise to a family of schemes, allowing us to compare their performances. The error analysis using induction technique in [19] is for the Dirichlet problem, whereas the present work is a periodic problem. This poses a challenge as some tools in the Dirichlet problem may not be available in the periodic case. We overcome this difficulty by combining the $L_2$ and $H_1$ analysis and choosing appropriate test functions.

To the best of our knowledge, there has not been any FDM work on the viscous Burgers-Poisson system, let alone providing analyses on the invariant-preserving property and convergence. The presented FDM scheme offers an easier alternative to the previous works while still achieving convergence and invariant-preserving properties. The ability to preserve invariant helps retain accuracy after a long time. This claim is based on the numerical evidence in [18] where the numerical solution obtained from the non-preserving method fails to capture the essence of the exact solution in the long run. A similar comparison will be presented in this work. The simplicity of the FDM method and the robustness of the proposed scheme allow a broader range of applications in computational simulations of the models. For instance, various types of initial data can be tested to learn the behavior of the shallow water waves.

The proposed method follows traditional FDM strategies: at a fixed point in the computational domain, each term in the equations is approximated by a divided difference quotient using information from nearby points. An alternate approach is to consider an average value on each sub-domain instead of value on each computational point. To get around the nonlinearity, some linearizing transform may be applied to the problem prior to the numerical step. See, for example, [20] where the Hopf-Cole transform is applied to a high order Burgers equation, then the 2D high order Spectral Volume (SV) is applied to the resulting equation. Stability analysis in SV can be found in [21,22] where Kannan and Wang conduct stability analyses on several viscous formulations for the higher SV method. When higher-order spatial derivative is involved, e.g., the KdV type equations, the high order SV scheme can also be employed. See [23,24] for examples.

The remaining part of this paper is organized as follows. In Section 2, we lay out the background of the study-case problem and discuss some analytic properties of the solution that are essential to the stability and error analyses. In Section 3, we outline basic settings for FDM framework and propose a collection of $\theta$-schemes to approximate the solution of the viscous Burger-Poisson system. In addition, we also prove the invariant properties of the proposed schemes. In Section 4, we present the error analysis to show that the method has second-order accuracy. In Section 5, we test the proposed scheme on various examples of the inviscous and viscous Burgers-Poisson equations to verify the theoretical results. In Section 6, we discuss the numerical results to obtain some insights into the properties of the proposed scheme. Some possible future works are also discussed. Finally, some concluding remarks are made in Section 7.

2. Analytic Property

This work formulates and analyzes a collection of three-level linearized schemes to solve the viscous Burgers-Poisson system

$$\begin{array}{cc}\hfill & u_t+u{u}_x-{\varphi }_x=ϵu_{xx},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & {\varphi }_{xx}-\varphi =u,(x,t)\in (a,b)\times (0,T),\hfill \end{array}$$

(10)

with initial condition

$$\begin{array}{c}\hfill u(x,0)=u_0\left(x\right),x\in (a,b),\end{array}$$

(11)

and periodic boundary conditions

$$\begin{array}{cc}\hfill & u(x+L,t)=u(x,t),\varphi (x+L,t)=\varphi (x,t),t>0,\hfill \end{array}$$

(12)

where $L=b-a$.

The nonlinear term $u{u}_x$ in (9) and other models of similar forms pose a challenge for the design of the scheme. In an early attempt to approximate the term $u{u}_x$, Guo and Sanz-Serna [25] employ a sum of two finite differences to approximate the term in the KdV equation. Following this idea, there have been several attempts to approximate the term $u{u}_x$ by splitting it into a sum of products of derivatives of u. Each term in the sum is then approximated by using implicit central difference with or without some average values. In the early years, two layers of numerical solution from consecutive time steps are involved in the schemes. See [26,27,28,29,30] for examples. It is also common to use three layers of numerical solutions. Having information from three consecutive time steps in the schemes helps eliminate the needs to solve nonlinear systems of unknowns in each time step. See [31,32,33] for examples. In [34], the authors use the term Rosenua-KdV-RLW to refer to all three equations and apply a convex average value of two implicit finite differences to approximate the nonlinear term. This gives a collection of second-order $\theta$-schemes proven to preserve invariant when $\theta =1/3$. We adapt their idea to the present problem. In the more recent works, Sun, Wang, and Zhang apply second-order [35] and fourth-order [19] operators to the viscous Burgers’ equation with zero boundary and obtain point-wise convergence analysis. We adapt their technique in the analysis, but some adjustments are needed because of the different types of boundary conditions.

Before we continue to present the scheme, we show that the solution of (9) and (10) with conditions (11) and (12) preserves the momentum and energy.

Theorem 1. 

Let u be a solution of (9)–(12). Define

$$Q\left(t\right)={\int }_a^{b}u(x,t)dx.$$

If $u_0\in L^2(a,b)$, then

$$Q\left(t\right)=Q\left(0\right),$$

for all $t>0$.

Proof. 

From (9), we have $u_t=-u{u}_x+{\varphi }_x+ϵu_{xx}$. Therefore,

$$\begin{array}{cc}\hfill \frac{d}{dt}Q\left(t\right)& ={\int }_a^b{u}_t(x,t)dx\hfill \\ \hfill & ={\int }_a^b-u(x,t)u_x(x,t)+{\varphi }_x(x,t)+ϵu_{xx}(x,t)dx\hfill \\ \hfill & ={\left[-\frac{1}{2}{u}^2(x,t)+\varphi (x,t)+ϵu_x(x,t)\right]}_a^b=0,\hfill \end{array}$$

because of the periodic boundary condition (12). Therefore, the function $Q\left(t\right)$ does not change over time. That is,

$$Q\left(t\right)=Q\left(0\right),$$

for all $t>0$. □

Theorem 2. 

Let u be a solution of (9)–(12). Define

$$E\left(t\right)={\int }_a^b{u}^2(x,t)dx.$$

If $u_0\in L^2(a,b)$, then

$$E\left(t\right)\le E\left(0\right),$$

for all $t>0$. The inequality becomes equality when $ϵ=0$.

Proof. 

Upon differentiating E, we obtain

$$\begin{array}{cc}\hfill \frac{1}{2}\frac{d}{dt}E\left(t\right)& ={\int }_a^{b}u(x,t)u_t(x,t)dx\hfill \\ \hfill & ={\int }_a^b-u^2(x,t)u_x(x,t)+u(x,t){\varphi }_x(x,t)+ϵu(x,t)u_{xx}(x,t)dx.\hfill \end{array}$$

For the second term, a use of (10) after taking inner product with ${\varphi }_x$ implies

$$\begin{array}{cc}\hfill \frac{1}{2}\frac{d}{dt}E\left(t\right)& ={\int }_a^b-u^2(x,t)u_x(x,t)+{\varphi }_{xx}(x,t){\varphi }_x(x,t)-\varphi (x,t){\varphi }_x(x,t)+ϵu(x,t)u_{xx}(x,t)dx\hfill \\ \hfill & ={\left[-\frac{1}{6}{u}^3(x,t)+\frac{1}{2}{\varphi }_x^2(x,t)-\frac{1}{2}{\varphi }^2(x,t)\right]}_a^b+{\int }_a^bϵu(x,t)u_{xx}(x,t)dx\hfill \\ \hfill & =-{\int }_a^bϵu_x^2(x,t)dx\le 0.\hfill \end{array}$$

Therefore, the function $E\left(t\right)$ does not increase over time, that is,

$$E\left(t\right)\le E\left(0\right),$$

for all $t>0$. If $ϵ=0$, then the proof above shows that $E^{\prime }\left(0\right)=0$ which gives $E\left(t\right)=E\left(0\right)$. This completes the proof. □

3. Finite Difference Method

3.1. Discretization

We discretize the spatial domain $[a,b]$ into the partition $x_i=a+ih,i=0,\dots ,M$ for the step size $h=(b-a)/M$. For a time step $\tau =T/N$, we define $t^n=n\tau$ for $n=0,\dots ,N$. Let $u_i^n$ be an approximation of $u(x_i,t^n)$. The numerical solution ${\mathbf{u}}^n={[u_1^n,\dots ,u_M^n]}^{\top }$ at the time $t^n$ is taken from the solution space $Z_h$ defined as

$$\begin{array}{cc}\hfill Z_h& :=\{\mathbf{u}\in {\mathcal{Z}}_h\mid u_j=u_{M+j},j=0,\dots ,M-1\},\hfill \end{array}$$

where ${\mathcal{Z}}_h$ is the set of discrete functions defined by

$$\begin{array}{cc}\hfill {\mathcal{Z}}_h& :=\{\mathbf{u}={[u_1,\dots ,u_M]}^{\top }\}.\hfill \end{array}$$

For any $\mathbf{u},\mathbf{v}\in {\mathcal{Z}}_h$, we define an inner product:

$$\begin{array}{cc}\hfill {\langle \mathbf{u},\mathbf{v}\rangle }_h& =h\sum _{i=1}^M{u}_i{v}_i\hfill \end{array}$$

which allows us to define the discrete $L_2$ norm

$$\begin{array}{cc}\hfill {\parallel \mathbf{u}\parallel }_h^2& ={\langle \mathbf{u},\mathbf{u}\rangle }_h.\hfill \end{array}$$

(13)

We also define the discrete uniform norm

$$\begin{array}{cc}\hfill {\parallel \mathbf{u}\parallel }_{h,\infty }& =\underset{i=1,\dots ,M}{max}|u_i|,\hfill \end{array}$$

(14)

for later use in the analysis.

With the set up above, we define finite difference methods for approximating some partial derivatives and some averages of u at $(x_i,t^n)$ as follows:

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\tau }^{+}{u}_i^n& =\frac{u_i^{n+1}-u_i^n}{\tau },{\mathcal{D}}_{\tau }^0{u}_i^n=\frac{u_i^{n+1}-u_i^{n-1}}{2\tau },\hfill \\ \hfill u_i^{n+1/2}& =\frac{u_i^{n+1}+u_i^n}{2},{\overline{u}}_i^n=\frac{u_i^{n+1}+u_i^{n-1}}{2},\hfill \\ \hfill D_h^{+}{u}_i^n& =\frac{u_{i+1}^n-u_i^n}{h},D_h^0{u}_i^n=\frac{u_{i+1}^n-u_{i-1}^n}{2h},D_h^{-}{u}_i^n=\frac{u_i^n-u_{i-1}^n}{h}.\hfill \end{array}$$

Note that we can apply these operators to a vector in $Z_h$ by acting on each element of the vector.

Using the boundary conditions in the definition of $Z_h$ and summation by parts, one can easily prove the following results.

Lemma 1. 

Let $\mathbf{u},\mathbf{v}\in Z_h$. The following relations hold.

$$\begin{array}{cc}\hfill {\langle D_h^{+}\mathbf{u},\mathbf{v}\rangle }_h& =-{\langle \mathbf{u},D_h^{-}\mathbf{v}\rangle }_h\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\langle D_h^0\mathbf{u},D_h^0\mathbf{u}\rangle }_h& \le {\langle D_h^{+}\mathbf{u},D_h^{+}\mathbf{u}\rangle }_h\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \parallel D_h^0{\mathbf{u}\parallel }_{h,\infty }& \le \parallel D_h^{+}{\mathbf{u}\parallel }_{h,\infty }.\hfill \end{array}$$

(17)

3.2. Formulation of the FDM

First, we define

$$\begin{array}{cc}\hfill {\mathsf{\Psi }}_{\theta }{(\mathbf{u},\mathbf{v})}_i& =2\theta u_i{D}_h^0{v}_i+(1-\theta )D_h^0{\left(uv\right)}_i,\hfill \end{array}$$

(18)

for any $\mathbf{u},\mathbf{v}\in Z_h$. Using the finite difference framework above, we propose a family of implicit three-level finite difference schemes for solving the Burgers-Poisson system (9)–(10) with conditions (11) and (12) as following: for $\theta \in [0,1]$, find $\mathbf{u},\mathbf{\varphi }\in Z_h$ satisfying

$$\begin{array}{cc}\hfill & {\mathcal{D}}_{\tau }^{+}{u}_i^0+\frac{1}{2}{\mathsf{\Psi }}_{\theta }{({\mathbf{u}}^0,{\mathbf{u}}^{1/2})}_i-D_h^0{\varphi }_i^0=ϵD_h^{+}{D}_h^{-}{u}_i^{1/2}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & D_h^{+}{D}_h^{-}{\varphi }_i^0-{\varphi }_i^0=u_i^{1/2}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & u_i^0=u_0\left(x_i\right),\hfill \end{array}$$

(21)

for $i=1,\dots ,M$, and

$$\begin{array}{cc}\hfill & {\mathcal{D}}_{\tau }^0{u}_i^n+\frac{1}{2}{\mathsf{\Psi }}_{\theta }{({\mathbf{u}}^n,{\overline{\mathbf{u}}}^n)}_i-D_h^0{\varphi }_i^n=ϵD_h^{+}{D}_h^{-}{\overline{u}}_i^n,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & D_h^{+}{D}_h^{-}{\varphi }_i^n-{\varphi }_i^n={\overline{u}}_i^n,\hfill \end{array}$$

(23)

for $i=1,\dots ,M,n=1,\dots ,N-1$.

For the implementation of the three-level scheme, the first two initial steps of the solutions are required. Thus, we use (19)–(21) to compute ${\mathbf{u}}^1$. Then, we use (22) and (23) to compute ${\mathbf{u}}^n$, for $n=2,3,\dots ,N$. Given ${\mathbf{u}}^n$ in each step, we use (20) or (23) to solve for ${\mathbf{\varphi }}^n$ first. Then, we use ${\mathbf{\varphi }}^n$ and ${\mathbf{u}}^n$ to solve for ${\mathbf{u}}^{n+1}$ using (19) or (22).

In the initial step, the system (20), with $i=1,\dots ,M$, can be written in a matrix-vector form

$$\begin{array}{cc}\hfill X_h{\mathbf{\varphi }}^0& ={\mathbf{u}}^{1/2},\mathrm{for}i=1,\dots ,M.\hfill \end{array}$$

(24)

where $X_h$ is a circulant tridiagonal matrix with $(1/h^2,-2/h^2-1,1/h^2)$ on the diagonal positions. It can be shown that $X_h$ is strictly diagonally dominant, thus its inverse exists. This allows us to write the system (19), $i=1,\dots ,M$, as a matrix-vector form

$$\begin{array}{cc}\hfill & {\mathbf{u}}^1+\frac{\tau \theta }{4h}{A}^0{\mathbf{u}}^1+\frac{\tau (1-\theta )}{8h}{B}^0{\mathbf{u}}^1-\left(\tau Y_h\right)\left(\frac{1}{2}{X}_h^{-1}{\mathbf{u}}^1\right)-\frac{\tau ϵ}{2}{Z}_h{\mathbf{u}}^1\hfill \\ \hfill & ={\mathbf{u}}^0-\frac{\tau \theta }{4h}{A}^0{\mathbf{u}}^0-\frac{\tau (1-\theta )}{8h}{B}^0{\mathbf{u}}^0+\left(\tau Y_h\right)\left(\frac{1}{2}{X}_h^{-1}{\mathbf{u}}^0\right)+\frac{\tau ϵ}{2}{Z}_h{\mathbf{u}}^0.\hfill \end{array}$$

(25)

For $n=1,2,\dots ,N-1$, the system (22), $i=1,\dots ,M$, can be written as

$$\begin{array}{cc}\hfill & {\mathbf{u}}^{n+1}+\frac{\tau \theta }{2h}{A}^n{\mathbf{u}}^{n+1}+\frac{\tau (1-\theta )}{4h}{B}^n{\mathbf{u}}^{n+1}-\left(2\tau Y_h\right)\left(\frac{1}{2}{X}_h^{-1}{\mathbf{u}}^{n+1}\right)-\tau ϵZ_h{\mathbf{u}}^{n+1}\hfill \\ \hfill & ={\mathbf{u}}^{n-1}-\frac{\tau \theta }{2h}{A}^n{\mathbf{u}}^{n-1}-\frac{\tau (1-\theta )}{4h}{B}^n{\mathbf{u}}^{n-1}+\left(2\tau Y_h\right)\left(\frac{1}{2}{X}_h^{-1}{\mathbf{u}}^{n-1}\right)+\tau ϵZ_h{\mathbf{u}}^{n-1}.\hfill \end{array}$$

(26)

Here, the coefficients $A^n,B^n,Y_h$, and $Z_h$ are circulant tridiagonal matrices whose nonzero entries on the ith row are given by $(-u_i^n,0,u_i^n)$, $(-u_{i-1}^n,0,u_{i+1}^n)$, $(-1/2h,0,1/2h)$, and $(1/h^2,-2,1/h^2)$, respectively.

3.3. Stability Analysis

Let $u_i^n$ be a numerical solution of (9) and (10) obtained from the proposed schemes (19)–(23). We prove that it preserves the invariants in the discrete sense.

Theorem 3. 

Define

$$\begin{array}{cc}\hfill Q_h^n& =h\sum _{i=1}^M{u}_i^{n+1/2}+\frac{\theta \tau }{2}h\sum _{i=1}^M{u}_i^n{D}_h^0{u}_i^{n+1}.\hfill \end{array}$$

(27)

If ${\mathbf{u}}^n$ is a solution of (19)–(23), then

$$\begin{array}{cc}\hfill Q_h^n& =Q_h^{n-1}=\dots =Q_h^0=h\sum _{i=1}^M{u}_i^0+\frac{\theta \tau }{4}h\sum _{i=1}^M{u}_i^0{D}_h^0{u}_i^1,\hfill \end{array}$$

(28)

for any $\theta \in [0,1]$.

Proof. 

Multiply (22) by $\tau h$ and take summation on $i=1,\dots ,M$ to arrive at

$$\begin{array}{cc}\hfill & \frac{1}{2}h\sum _{i=1}^M\left(u_i^{n+1}-u_i^{n-1}\right)+\theta \tau h\sum _{i=1}^M{u}_i^n{D}_h^0{\overline{u}}_i^n+\frac{(1-\theta )}{2}\tau h\sum _{i=1}^M{D}_h^0\left(u_i^n{\overline{u}}_i^n\right)\hfill \\ \hfill & -\tau h\sum _{i=1}^M{D}_h^0{\varphi }_i^n=ϵ\tau h\sum _{i=1}^M{D}_h^{+}{D}_h^{-}{\overline{u}}_i^n.\hfill \end{array}$$

Since all terms other than the first two vanish, we obtain

$$\begin{array}{cc}\hfill \frac{1}{2}h\sum _{i=1}^M\left(u_i^{n+1}-u_i^{n-1}\right)+\theta \tau h\sum _{i=1}^M{u}_i^n{D}_h^0{\overline{u}}_i^n& =0,\hfill \end{array}$$

and hence,

$$\begin{array}{cc}\hfill \frac{1}{2}h\sum _{i=1}^M\left(u_i^{n+1}-u_i^{n-1}\right)+\frac{\theta \tau }{2}h\sum _{i=1}^M\left(u_i^n{u}_{i+1}^{n+1}-u_i^n{u}_{i+1}^{n+1}+u_i^n{u}_{i-1}^{n-1}-u_i^n{u}_{i-1}^{n-1}\right)& =0.\hfill \end{array}$$

After rewriting,

$$\begin{array}{cc}\hfill \frac{1}{2}h\sum _{i=1}^M\left(u_i^{n+1}-u_i^{n-1}\right)+\frac{\theta \tau }{2}h\sum _{i=1}^M\left(u_i^n{u}_{i+1}^{n+1}-u_i^n{u}_{i+1}^{n+1}+u_{i-1}^n{u}_i^{n-1}-u_{i+1}^n{u}_i^{n-1}\right)& =0,\hfill \end{array}$$

and thus,

$$\begin{array}{cc}\hfill \frac{1}{2}h\sum _{i=1}^M\left(u_i^{n+1}-u_i^{n-1}\right)+\frac{\theta \tau }{2}h\sum _{i=1}^M\left(u_i^n{D}_h^0{u}_i^{n+1}-u_i^{n-1}{D}_h^0{u}_i^n\right)& =0.\hfill \end{array}$$

This shows that $Q_h^n-Q_h^{n-1}=0$, for $n=1,\dots ,N-1$. Similarly, multiply (19) by $\tau h$ and take summation over i to obtain

$$\begin{array}{cc}\hfill & h\sum _{i=1}^M{u}_i^1=h\sum _{i=1}^M{u}_i^0-\frac{\theta \tau }{2}h\sum _{i=1}^M{u}_i^0{D}_h^0{u}_i^1=0.\hfill \end{array}$$

Substitute this into $Q_h^0$ to find that

$$Q_h^0=h\sum _{i=1}^M{u}_i^0+\frac{\theta \tau }{4}h\sum _{i=1}^M{u}_i^0{D}_h^0{u}_i^1,$$

as needed. □

Theorem 4. 

Define

$$\begin{array}{cc}\hfill E_h^n& =\frac{1}{2}\parallel {\mathbf{u}}^{n+1}{\parallel }_h^2+\frac{1}{2}{\parallel {\mathbf{u}}^n\parallel }_h^2.\hfill \end{array}$$

(29)

Let ${\mathbf{u}}^n$ be a solution of (19) and (23). If $\theta ={\displaystyle \frac{1}{3}}$, then

$$\begin{array}{cc}\hfill & E_h^n+2\tau ϵ\sum _{k=1}^n\parallel {D_h^{+}\overline{\mathbf{u}}}^k{\parallel }_h^2+\tau ϵ\parallel D_h^{+}{\mathbf{u}}^{1/2}{\parallel }_h^2={\parallel {\mathbf{u}}^0\parallel }_h^2,\hfill \end{array}$$

(30)

for any $n>0$.

Proof. 

First, consider the summation

$$\begin{array}{cc}\hfill & h\sum _{i=1}^M\left[u_i^n\left(D_h^0{\overline{u}}_i^n\right){\overline{u}}_i^n+D_h^0\left(u_i^n{\overline{u}}_i^n\right){\overline{u}}_i^n\right]\hfill \\ \hfill & =h\sum _{i=1}^M\left[u_i^n{\overline{u}}_{i+1}^n{\overline{u}}_i^n-u_i^n{\overline{u}}_{i-1}^n{\overline{u}}_i^n+u_{i+1}^n{\overline{u}}_{i+1}^n{\overline{u}}_i^n-u_{i-1}^n{\overline{u}}_{i-1}^n{\overline{u}}_i^n\right]\hfill \\ \hfill & =h\sum _{i=1}^M\left[u_i^n{\overline{u}}_{i+1}^n{\overline{u}}_i^n-u_{i-1}^n{\overline{u}}_{i-1}^n{\overline{u}}_i^n\right]+h\sum _{i=1}^M\left[-u_i^n{\overline{u}}_{i-1}^n{\overline{u}}_i^n+u_{i+1}^n{\overline{u}}_{i+1}^n{\overline{u}}_i^n\right]=0.\hfill \end{array}$$

(31)

Using (31), we multiply (22) by $2\tau h{\overline{u}}_i^n$ and take summation over i to obtain

$$\begin{array}{cc}\hfill & \frac{1}{2}h\sum _{i=1}^M\left[{\left(u_i^{n+1}\right)}^2-{\left(u_i^{n-1}\right)}^2\right]+\tau (3\theta -1)h\sum _{i=1}^M{D}_h^0\left(u_i^n{\overline{u}}_i^n\right){\overline{u}}_i^n-2\tau h\sum _{i=1}^M\left(D_h^0{\varphi }_i^n\right){\overline{u}}_i^n\hfill \\ \hfill & +2ϵ\tau h\sum _{i=1}^M{\left(D_h^{+}{u}_i^n\right)}^2=0.\hfill \end{array}$$

(32)

Next, we multiply (23) by $h{D}_h^0{\varphi }_i^n$ and take the summation over i. Together with Lemma 1, this yields

$$\begin{array}{cc}\hfill h\sum _{i=1}^M{\overline{u}}_i^n\left(D_h^0{\varphi }_i^n\right)& =h\sum _{i=1}^M\left(D_h^{+}{D}_h^{-}{\varphi }_i^n\right)\left(D_h^0{\varphi }_i^n\right)-h\sum _{i=1}^M{\varphi }_i^n\left(D_h^0{\varphi }_i^n\right)\hfill \\ \hfill & ={\langle D_h^{+}{D}_h^{-}{\mathbf{\varphi }}^n,D_h^0{\mathbf{\varphi }}^n\rangle }_h-{\langle {\mathbf{\varphi }}^n,D_h^0{\mathbf{\varphi }}^n\rangle }_h=0.\hfill \end{array}$$

(33)

By taking $\theta ={\displaystyle \frac{1}{3}}$, from (32) and (33), we arrive at

$$\begin{array}{cc}\hfill & \frac{1}{2}h\sum _{i=1}^M\left[{\left(u_i^{n+1}\right)}^2-{\left(u_i^{n-1}\right)}^2\right]+2ϵ\tau h\sum _{i=1}^M{\left(D_h^{+}{\overline{u}}_i^n\right)}^2=0.\hfill \end{array}$$

(34)

This shows that

$$\begin{array}{cc}\hfill & E_h^n-E_h^{n-1}+2ϵ\tau h\sum _{i=1}^M{\left(D_h^{+}{\overline{u}}_i^n\right)}^2=0,\hfill \end{array}$$

(35)

for $n=1,\dots ,N-1$.

Using similar idea, one can show from (19)–(21) that

$$\begin{array}{cc}\hfill & \parallel {\mathbf{u}}^1{\parallel }_h^2-{\parallel {\mathbf{u}}^0\parallel }_h^2+ϵ\tau h\sum _{i=1}^M{\left(D_h^{+}{u}_i^{1/2}\right)}^2=0.\hfill \end{array}$$

(36)

Substitute this into the definition of $E_h^0$ to derive the desired result. □

Remark 1. 

As a consequence of Theorem 4, we obtain

$$\parallel {\mathbf{u}}^n{\parallel }_h^2\le \parallel {\mathbf{u}}^{n-1}{\parallel }_h^2\le \dots \le {\parallel {\mathbf{u}}^0\parallel }_h^2.$$

This, in turns, guarantees the existence and uniqueness of the solution of the schemes (19)–(23).

4. Convergence Analysis

In this section, we study the error analysis of the numerical solution obtained from the proposed scheme. The following results are necessary for the error estimation.

Lemma 2 

(Discrete Sobolev’s inequality [36]). If ${\mathbf{u}}^n\in Z_h$, then there exists a constant α, depending only on L, such that

$$\parallel {\mathbf{u}}^n{\parallel }_{h,\infty }\le \alpha \left(\parallel {\mathbf{u}}^n{\parallel }_h+{\parallel D_h^{+}{\mathbf{u}}^n\parallel }_h\right).$$

Let $v_i^n$ and ${\phi }_i^n$ be the exact solutions at the point $(x_i,t^n)$. Define $e_i^n=v_i^n-u_i^n$ and ${\epsilon }_i^n={\phi }_i^n-{\varphi }_i^n$. From (19)–(21), we have that the truncation errors at the time step $t^0$, denoted ${\mathcal{T}}^0$ and ${\mathcal{R}}^0$, satisfy

$$\begin{array}{cc}\hfill & {\mathcal{D}}_{\tau }^{+}{e}_i^0+\frac{1}{2}\left[{\mathsf{\Psi }}_{\theta }{({\mathbf{v}}^0,{\mathbf{v}}^{1/2})}_i-{\mathsf{\Psi }}_{\theta }{({\mathbf{u}}^0,{\mathbf{u}}^{1/2})}_i\right]-D_h^0{\epsilon }_i^0=ϵD_h^{+}{D}_h^{-}{e}_i^{1/2}+{\mathcal{T}}_i^0,\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & D_h^{+}{D}_h^{-}{\epsilon }_i^0-{\epsilon }_i^0=e_i^{1/2}+{\mathcal{R}}_i^0.\hfill \end{array}$$

(38)

Using Taylor series expansion, one can show there exists $c_1$ such that

$$\begin{array}{cc}\hfill & |{\mathcal{T}}_i^0|\le c_1(\tau +h^2),\mathrm{and}\left|{\mathcal{R}}_i^0\right|\le c_1{h}^2.\hfill \end{array}$$

On the other hand, we have from (22) and (23) that the truncation errors at the time step $t^n,n\ge 1$, denoted ${\mathcal{T}}^n$ and ${\mathcal{R}}^n$, satisfy

$$\begin{array}{cc}\hfill & {\mathcal{D}}_{\tau }^0{e}_i^n+\frac{1}{2}\left[{\mathsf{\Psi }}_{\theta }{({\mathbf{v}}^n,{\overline{\mathbf{v}}}^n)}_i-{\mathsf{\Psi }}_{\theta }{({\mathbf{u}}^n,{\overline{\mathbf{u}}}^n)}_i\right]-D_h^0{\epsilon }_i^n=ϵD_h^{+}{D}_h^{-}{\overline{e}}_i^n+{\mathcal{T}}_i^n,\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill & D_h^{+}{D}_h^{-}{\epsilon }_i^n-{\epsilon }_i^n={\overline{e}}_i^n+{\mathcal{R}}_i^n.\hfill \end{array}$$

(40)

Using Taylor series expansion, one can show there exists $c_2$ such that

$$\begin{array}{cc}\hfill & |{\mathcal{T}}_i^n|\le c_2({\tau }^2+h^2),\mathrm{and}\left|{\mathcal{R}}_i^n\right|\le c_2{h}^2.\hfill \end{array}$$

The next theorem establishes a relation between the errors $e_i^n$ and ${\epsilon }_i^n$.

Theorem 5. 

Let $e_i^n$ and ${\epsilon }_i^n$ be solutions of (37)–(40). The following inequalities hold for $n=0,1,\dots ,N$.

$$\begin{array}{cc}\hfill \parallel {\mathbf{\epsilon }}^n{\parallel }_h& \le \parallel {\mathbf{e}}^{\zeta }{\parallel }_h+{\parallel {\mathcal{R}}^n\parallel }_h\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill \parallel D_h^{+}{\mathbf{\epsilon }}^n{\parallel }_h& \le \parallel {\mathbf{e}}^{\zeta }{\parallel }_h+{\parallel {\mathcal{R}}^n\parallel }_h\hfill \end{array}$$

(42)

Here, ${\mathbf{e}}^{\zeta }={\mathbf{e}}^{1/2}$ when $n=0$. Otherwise, ${\mathbf{e}}^{\zeta }={\overline{\mathbf{e}}}^n$.

Proof. 

We have from (38) and (40) that

$$X_h{\mathbf{\epsilon }}^n={\mathbf{e}}^{\zeta }+{\mathcal{R}}^n.$$

Since the matrix $X_h$ has eigenvalues of the form ${\lambda }_i=-\left[1+(2/h)(1-cos(\frac{i\pi }{M+1}))\right]$, $i=1,2,\dots ,M$, we have $\parallel X_h^{-1}{\parallel }_h\le 1$. Therefore,

$$\parallel {\mathbf{\epsilon }}^n{\parallel }_h\le \parallel X_h^{-1}{\parallel }_h\parallel {\mathbf{e}}^{\zeta }+{\mathcal{R}}^n{\parallel }_h\le \parallel {\mathbf{e}}^{\zeta }{\parallel }_h+{\parallel {\mathcal{R}}^n\parallel }_h,$$

which proves (41). Using Lemma (1), we have from (38) and (40) that

$$\begin{array}{cc}\hfill \parallel D_h^{+}{\mathbf{\epsilon }}^n{\parallel }_h^2& ={\langle D_h^{+}{\mathbf{\epsilon }}^n,D_h^{+}{\mathbf{\epsilon }}^n\rangle }_h=-{\langle D_h^{+}{D}_h^{-}{\mathbf{\epsilon }}^n,{\mathbf{\epsilon }}^n\rangle }_h\hfill \\ \hfill & =-{\langle {\mathbf{\epsilon }}^n+{\mathbf{e}}^{\zeta }+{\mathcal{R}}^n,{\mathbf{\epsilon }}^n\rangle }_h\le -{\langle {\mathbf{e}}^{\zeta },{\mathbf{\epsilon }}^n\rangle }_h-{\langle {\mathcal{R}}^n,{\mathbf{\epsilon }}^n\rangle }_h\hfill \\ \hfill & \le \left(\parallel {\mathbf{e}}^{\zeta }{\parallel }_h+{\parallel {\mathcal{R}}^n\parallel }_h\right){\parallel {\mathbf{\epsilon }}^n\parallel }_h.\hfill \end{array}$$

By (41), we get

$$\parallel D_h^{+}{\mathbf{\epsilon }}^n{\parallel }_h^2\le {\left(\parallel {\mathbf{e}}^{\zeta }{\parallel }_h+{\parallel {\mathcal{R}}^n\parallel }_h\right)}^2$$

which proves (42) as needed. □

Theorem 6. 

Define

$$c_3=\underset{x\in [a,b],t\in [0,T]}{max}\left\{\right|v(x,t)|,|v_x(x,t)\left|\right\},$$

$$c_4=\frac{c_3ϵ\theta }{4}+\frac{c_3ϵ(1-\theta )}{8}+\frac{c_3^2{(1-\theta )}^2}{4}+max\left\{\frac{c_3^2{\theta }^2}{2},\frac{3ϵ+1}{2}+\frac{c_3ϵ(1-\theta )}{4}\right\},$$

$$c_5=2(ϵ+8)c_2,c_6=\frac{c_3ϵ(1-\theta )}{4}+2c_3^2{(1-\theta )}^2+\frac{1}{2}max\left\{A_1,A_2\right\}where$$

$$A_1=c_3ϵ\theta +8c_3^2{\theta }^2+8{\theta }^2+8{(1-\theta )}^2$$

$$A_2=3ϵ+4+ϵ(1-\theta )+{(1-\theta )}^2,$$

$$c_7=c_1(2ϵ+6),c_8=2\alpha \sqrt{\left(2\frac{c_5}{ϵ}\right)exp\left(\frac{4c_{6}T}{ϵ}\right)}.$$

If $({\tau }^{3/2}+h^2)\le min\left\{\frac{1}{2\alpha }\sqrt{\frac{ϵ}{c_5}},\frac{1}{c_8}\right\}$ and $\tau \le min\left\{1,\frac{ϵ}{2c_4},\frac{ϵ}{4c_6}\right\}$, then the solution ${\mathbf{u}}^n$ of the scheme (19)–(23) converges in the sense

$$\parallel {\mathbf{e}}^n{\parallel }_{h,\infty }\le C_8({\tau }^{3/2}+h^2).$$

Proof. 

Since ${\mathbf{e}}^0=\mathbf{0}$, we have

$$\begin{array}{cc}\hfill {\mathsf{\Psi }}_{\theta }({\mathbf{v}}^0,{\mathbf{v}}^{1/2})-{\mathsf{\Psi }}_{\theta }({\mathbf{u}}^0,{\mathbf{u}}^{1/2})& ={\mathsf{\Psi }}_{\theta }({\mathbf{v}}^0,{\mathbf{e}}^{1/2})+{\mathsf{\Psi }}_{\theta }({\mathbf{e}}^0,{\mathbf{v}}^{1/2})-{\mathsf{\Psi }}_{\theta }({\mathbf{e}}^0,{\mathbf{e}}^{1/2})\hfill \\ \hfill & ={\mathsf{\Psi }}_{\theta }({\mathbf{v}}^0,{\mathbf{e}}^{1/2}).\hfill \end{array}$$

Thus, the error Equation (37) is reduced to

$$\begin{array}{cc}\hfill & \frac{e_i^1}{\tau }+\frac{1}{4}{\mathsf{\Psi }}_{\theta }{({\mathbf{v}}^0,{\mathbf{e}}^1)}_i-D_h^0{\epsilon }_i^0=\frac{ϵ}{2}{D}_h^{+}{D}_h^{-}{e}_i^1+{\mathcal{T}}_i^0.\hfill \end{array}$$

(43)

Take the inner product of the system (43) where $i=1,\dots ,M$ with $ϵ{\mathbf{e}}^1+\frac{2}{\tau }{\mathbf{e}}^1$ to get

$$\begin{array}{cc}\hfill & {\langle \frac{e^1}{\tau },ϵ{\mathbf{e}}^1+\frac{2}{\tau }{\mathbf{e}}^1\rangle }_h-{\langle \frac{ϵ}{2}{D}_h^{+}{D}_h^{-}{\mathbf{e}}^1,ϵ{\mathbf{e}}^1+\frac{2}{\tau }{\mathbf{e}}^1\rangle }_h\hfill \\ \hfill & ={\langle {\mathcal{T}}^0,ϵ{\mathbf{e}}^1+\frac{2}{\tau }{\mathbf{e}}^1\rangle }_h+{\langle D_h^0{\mathbf{\epsilon }}^0,ϵ{\mathbf{e}}^1+\frac{2}{\tau }{\mathbf{e}}^1\rangle }_h-{\langle \frac{1}{4}{\mathsf{\Psi }}_{\theta }({\mathbf{v}}^0,{\mathbf{e}}^1),ϵ{\mathbf{e}}^1+\frac{2}{\tau }{\mathbf{e}}^1\rangle }_h.\hfill \end{array}$$

(44)

The left-hand side of (44) can be simplified into

$$\begin{array}{cc}\hfill & \frac{ϵ}{\tau }\parallel {\mathbf{e}}^1{\parallel }_h^2+2\parallel \frac{1}{\tau }{\mathbf{e}}^1{\parallel }_h^2+\frac{ϵ}{\tau }\parallel D_h^{+}{\mathbf{e}}^1{\parallel }_h^2+\frac{{ϵ}^2}{2}{\parallel D_h^{+}{\mathbf{e}}^1\parallel }_h^2.\hfill \end{array}$$

(45)

Using Lemma 1, Theorem 5, and Cauchy–Schwarz inequality, the first two terms on the right-hand side of (44) yield

$$\begin{array}{cc}\hfill & {\langle {\mathcal{T}}^0,ϵ{\mathbf{e}}^1+\frac{2}{\tau }{\mathbf{e}}^1\rangle }_h+{\langle D_h^0{\mathbf{\epsilon }}^0,ϵ{\mathbf{e}}^1+\frac{2}{\tau }{\mathbf{e}}^1\rangle }_h\hfill \\ \hfill & \le ϵ\parallel {\mathcal{T}}^0{\parallel }_h\parallel {\mathbf{e}}^1{\parallel }_h+2\parallel {\mathcal{T}}^0{\parallel }_h{\parallel \frac{1}{\tau }{\mathbf{e}}^1\parallel }_h\hfill \\ \hfill & +ϵ\left(\frac{1}{2}\parallel {\mathbf{e}}^1{\parallel }_h+{\parallel {\mathcal{R}}^0\parallel }_h\right)\parallel {\mathbf{e}}^1{\parallel }_h+2\left(\frac{1}{2}\parallel {\mathbf{e}}^1{\parallel }_h+{\parallel {\mathcal{R}}^0\parallel }_h\right){\parallel \frac{1}{\tau }{\mathbf{e}}^1\parallel }_h\hfill \\ \hfill & \le \left(\frac{ϵ}{2}\parallel {\mathcal{T}}^0{\parallel }_h^2+\frac{ϵ}{2}{\parallel {\mathbf{e}}^1\parallel }_h^2\right)+\left(4\parallel {\mathcal{T}}^0{\parallel }_h+\frac{1}{4}{\parallel \frac{1}{\tau }{\mathbf{e}}^1\parallel }_h\right)+\left(ϵ\parallel {\mathbf{e}}^1{\parallel }_h^2+\frac{ϵ}{2}{\parallel {\mathcal{R}}^0\parallel }_h^2\right)\hfill \\ \hfill & +\left(\frac{1}{2}\parallel {\mathbf{e}}^1{\parallel }_h^2+\frac{1}{2}\parallel \frac{1}{\tau }{\mathbf{e}}^1{\parallel }_h^2+4\parallel {\mathcal{R}}^0{\parallel }_h^2+\frac{1}{4}{\parallel \frac{1}{\tau }{\mathbf{e}}^1\parallel }_h^2\right)\hfill \\ \hfill & =\left(\frac{3ϵ}{2}+\frac{1}{2}\right)\parallel {\mathbf{e}}^1{\parallel }_h^2+{\parallel \frac{1}{\tau }{\mathbf{e}}^1\parallel }_h^2+\left(\frac{ϵ}{2}+4\right)\left(\parallel {\mathcal{T}}^0{\parallel }_h^2+{\parallel {\mathcal{R}}^0\parallel }_h^2\right).\hfill \end{array}$$

(46)

As for the last term in (44), note that

$$\begin{array}{cc}\hfill D_h^0{\left(uv\right)}_i& =\frac{1}{2}{u}_{i+1}{D}_h^{+}{v}_i+\left(D_h^0{u}_i\right)v_i+\frac{1}{2}{u}_{i-1}{D}_h^{+}{v}_{i-1}.\hfill \end{array}$$

(47)

Therefore, using Lemma 1, we arrive at

$$\begin{array}{cc}\hfill & -{\langle \frac{1}{4}{\mathsf{\Psi }}_{\theta }({\mathbf{v}}^0,{\mathbf{e}}^1),ϵ{\mathbf{e}}^1+\frac{2}{\tau }{\mathbf{e}}^1\rangle }_h=-\frac{\theta }{2}h\sum _{i=1}^M{v}_i^0{D}_h^0{e}_i^1(ϵe_i^1+\frac{2}{\tau }{e}_i^1)\hfill \\ \hfill & -\frac{1-\theta }{4}h\sum _{i=1}^M\left(\frac{1}{2}{v}_{i+1}^0{D}_h^{+}{e}_i^1+\left(D_h^0{v}_i^0\right)e_i^1+\frac{1}{2}{v}_{i-1}^0{D}_h^{+}{e}_{i-1}^1\right)(ϵe_i^1+\frac{2}{\tau }{e}_i^1)\hfill \\ \hfill & \le \left[\frac{c_3\theta }{2}{\parallel D_h^{+}{\mathbf{e}}^1\parallel }_h+\frac{c_3(1-\theta )}{4}\left(\parallel D_h^{+}{\mathbf{e}}^1{\parallel }_h+{\parallel {\mathbf{e}}^1\parallel }_h\right)\right]\left(ϵ\parallel {\mathbf{e}}^1{\parallel }_h+2{\parallel \frac{1}{\tau }{\mathbf{e}}^1\parallel }_h\right)\hfill \\ \hfill & \le \left(\frac{c_3ϵ\theta }{4}+\frac{1}{2}{c}_3^2{\theta }^2+\frac{c_3ϵ(1-\theta )}{8}+\frac{c_3^2{(1-\theta )}^2}{4}\right){\parallel D_h^{+}{\mathbf{e}}^1\parallel }_h^2\hfill \\ \hfill & +\left(\frac{c_3ϵ\theta }{4}+\frac{3c_3ϵ(1-\theta )}{8}+\frac{c_3^2{(1-\theta )}^2}{4}\right)\parallel {\mathbf{e}}^1{\parallel }_h^2+\frac{1}{2}{\parallel \frac{1}{\tau }{\mathbf{e}}^1\parallel }_h^2.\hfill \end{array}$$

(48)

Substituting (45), (46), and (48) into (44), we obtain

$$\begin{array}{cc}\hfill \frac{ϵ}{\tau }(\parallel {\mathbf{e}}^1{\parallel }_h^2+\parallel D_h^{+}{\mathbf{e}}^1{\parallel }_h^2)& \le c_4(\parallel {\mathbf{e}}^1{\parallel }_h^2+\parallel D_h^{+}{\mathbf{e}}^1{\parallel }_h^2)+\frac{c_5}{2}{(\tau +h^2)}^2\hfill \\ \hfill \left(1-\frac{c_4}{ϵ}\tau \right)(\parallel {\mathbf{e}}^1{\parallel }_h^2+\parallel D_h^{+}{\mathbf{e}}^1{\parallel }_h^2)& \le \frac{c_5}{2ϵ}\tau {(\tau +h^2)}^2.\hfill \end{array}$$

For $\frac{c_4}{ϵ}\tau \le \frac{1}{2}$, we have

$$\begin{array}{c}\hfill \parallel {\mathbf{e}}^1{\parallel }_h^2+{\parallel D_h^{+}{\mathbf{e}}^1\parallel }_h^2\le \frac{c_5}{ϵ}\tau {(\tau +h^2)}^2.\end{array}$$

(49)

From this, we get

$$\begin{array}{cc}\hfill & \parallel {\mathbf{e}}^1\parallel \le \sqrt{\frac{ϵ}{c_5}}({\tau }^{3/2}+h^2)and\parallel D_h^{+}{\mathbf{e}}^1\parallel \le \sqrt{\frac{ϵ}{c_5}}({\tau }^{3/2}+h^2),\hfill \end{array}$$

(50)

if $\tau <1$. From the hypothesis of the theorem, we have $\parallel {\mathbf{e}}^1{\parallel }_h+{\parallel D_h^{+}{\mathbf{e}}^1\parallel }_h\le \frac{1}{\alpha }$. This shows that $\parallel {\mathbf{e}}^1{\parallel }_{h,\infty }\le 1$. For the time step $t^n,n>0$, we use induction on n and proceed in a similar manner. Assume

$$\begin{array}{cc}\hfill & \parallel {\mathbf{e}}^k{\parallel }_h+{\parallel D_h^{+}{\mathbf{e}}^k\parallel }_h\le \frac{1}{\alpha },\hfill \end{array}$$

(51)

that is, $\parallel {\mathbf{e}}^k{\parallel }_{h,\infty }\le 1$, for $k=0,1,\dots ,n$. We shall show that (51) holds for $k=n+1$. First, we take inner product of (39), where $i=1,\dots ,M$, with $2ϵ{\overline{\mathbf{e}}}^n+2{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n$ to get

$$\begin{array}{cc}\hfill & {\langle {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n,2ϵ{\overline{\mathbf{e}}}^n+2{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h-2ϵ{\langle D_h^{+}{D}_h^{-}{\overline{\mathbf{e}}}^n,ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h\hfill \\ \hfill & =2{\langle {\mathcal{T}}^n,ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h+2{\langle D_h^0{\mathbf{\epsilon }}^n,ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h\hfill \\ \hfill & -{\langle {\mathsf{\Psi }}_{\theta }({\mathbf{v}}^n,{\overline{\mathbf{v}}}^n)-{\mathsf{\Psi }}_{\theta }({\mathbf{u}}^n,{\overline{\mathbf{u}}}^n),ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h.\hfill \end{array}$$

(52)

The left-hand side of (52) can be simplified into

$$\begin{array}{cc}\hfill & \frac{ϵ}{2\tau }\left(\parallel {\mathbf{e}}^{n+1}{\parallel }_h^2-{\parallel {\mathbf{e}}^{n-1}\parallel }_h^2\right)+2{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h^2\hfill \\ \hfill & +\frac{ϵ}{2\tau }\left(\parallel D_h^{+}{\mathbf{e}}^{n+1}{\parallel }_h^2-{\parallel D_h^{+}{\mathbf{e}}^{n-1}\parallel }_h^2\right)+2{ϵ}^2{\parallel {\overline{\mathbf{e}}}^n\parallel }_h^2.\hfill \end{array}$$

(53)

Using Lemma 1 and Theorem 5, the first two terms on the right-hand side of (52) satisfy

$$\begin{array}{cc}\hfill 2{\langle {\mathcal{T}}^n,ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h& \le 2ϵ\parallel {\mathcal{T}}^n{\parallel }_h\parallel {\overline{\mathbf{e}}}^n{\parallel }_h+2\parallel {\mathcal{T}}^n{\parallel }_h{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h\hfill \\ \hfill & \le (ϵ+2)\parallel {\mathcal{T}}^n{\parallel }_h^2+ϵ\parallel {\overline{\mathbf{e}}}^n{\parallel }_h^2+\frac{1}{2}{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h^2\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill 2{\langle D_h^0{\mathbf{\epsilon }}^n,ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h& \le 2\parallel D_h^{+}{\mathbf{\epsilon }}^n{\parallel }_h\left(ϵ\parallel {\overline{\mathbf{e}}}^n{\parallel }_h+{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h\right)\hfill \\ \hfill & \le 2\left(\parallel {\overline{\mathbf{e}}}^n{\parallel }_h+{\parallel {\mathcal{R}}^n\parallel }_h\right)\left(ϵ\parallel {\overline{\mathbf{e}}}^n{\parallel }_h+{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h\right)\hfill \\ \hfill & \le \left(3ϵ+4\right)\parallel {\overline{\mathbf{e}}}^n{\parallel }_h^2+(ϵ+4)\parallel {\mathcal{R}}^n{\parallel }_h^2+\frac{1}{2}{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h^2.\hfill \end{array}$$

(55)

As for the last term in (52), note that

$$\begin{array}{cc}\hfill {\mathsf{\Psi }}_{\theta }({\mathbf{v}}^n,{\overline{\mathbf{v}}}^n)-{\mathsf{\Psi }}_{\theta }({\mathbf{u}}^n,{\overline{\mathbf{u}}}^n)& ={\mathsf{\Psi }}_{\theta }({\mathbf{v}}^n,{\overline{\mathbf{e}}}^n)+{\mathsf{\Psi }}_{\theta }({\mathbf{e}}^n,{\overline{\mathbf{v}}}^n)-{\mathsf{\Psi }}_{\theta }({\mathbf{e}}^n,{\overline{\mathbf{e}}}^n).\hfill \end{array}$$

Using (47) and Lemma 1, we have

$$\begin{array}{cc}\hfill & -{\langle {\mathsf{\Psi }}_{\theta }({\mathbf{v}}^n,{\overline{\mathbf{e}}}^n),ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h=-2\theta h\sum _{i=1}^M{v}_i^n{D}_h^0{\overline{e}}_i^n(ϵ{\overline{e}}_i^n+{\mathcal{D}}_{\tau }^0{e}_i^n)\hfill \\ \hfill & -(1-\theta )h\sum _{i=1}^M\left(\frac{1}{2}{v}_{i+1}^n{D}_h^{+}{\overline{e}}_i^n+\left(D_h^0{v}_i^n\right){\overline{e}}_i^n+\frac{1}{2}{v}_{i-1}^n{D}_h^{+}{\overline{e}}_{i-1}^n\right)(ϵ{\overline{e}}_i^n+{\mathcal{D}}_{\tau }^0{e}_i^n)\hfill \\ \hfill & \le \left(2c_3\theta {\parallel D_h^{+}{\overline{\mathbf{e}}}^n\parallel }_h+c_3(1-\theta )\left(\parallel D_h^{+}{\overline{\mathbf{e}}}^n{\parallel }_h+{\parallel {\overline{\mathbf{e}}}^n\parallel }_h\right)\right)\left(ϵ\parallel {\overline{\mathbf{e}}}^n{\parallel }_h+{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h\right)\hfill \\ \hfill & \le \left(c_3ϵ\theta +8c_3^2{\theta }^2+\frac{c_3ϵ(1-\theta )}{2}+4c_3^2{(1-\theta )}^2\right){\parallel D_h^{+}{\overline{\mathbf{e}}}^n\parallel }_h^2\hfill \\ \hfill & +\left(c_3ϵ\theta +\frac{3c_3ϵ(1-\theta )}{2}+4c_3^2{(1-\theta )}^2\right)\parallel {\overline{\mathbf{e}}}^n{\parallel }_h^2+\frac{1}{4}{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h^2.\hfill \end{array}$$

(56)

Similarly, we have

$$\begin{array}{cc}\hfill & -{\langle {\mathsf{\Psi }}_{\theta }({\mathbf{e}}^n,{\overline{\mathbf{v}}}^n),ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h\le \left(c_3ϵ\theta +c_3ϵ(1-\theta )\right){\parallel {\overline{\mathbf{e}}}^n\parallel }_h^2\hfill \\ \hfill & +\left(\frac{c_3ϵ(1-\theta )}{2}+4c_3^2{(1-\theta )}^2\right){\parallel D_h^{+}{\mathbf{e}}^n\parallel }_h^2\hfill \\ \hfill & +\left(c_3ϵ\theta +8c_3^2{\theta }^2+\frac{c_3ϵ(1-\theta )}{2}+4c_3^2{(1-\theta )}^2\right)\parallel {\mathbf{e}}^n{\parallel }_h^2+\frac{1}{4}{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h^2.\hfill \end{array}$$

(57)

Using (47), (51), and Lemma 1, we find that

$$\begin{array}{cc}\hfill & {\langle {\mathsf{\Psi }}_{\theta }({\mathbf{e}}^n,{\overline{\mathbf{e}}}^n),ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h=2\theta h\sum _{i=1}^M{e}_i^n{D}_h^0{\overline{e}}_i^n(ϵ{\overline{e}}_i^n+{\mathcal{D}}_{\tau }^0{e}_i^n)\hfill \\ \hfill & +(1-\theta )h\sum _{i=1}^M\left(\frac{1}{2}{e}_{i+1}^n{D}_h^{+}{\overline{e}}_i^n+\left(D_h^0{e}_i^n\right){\overline{e}}_i^n+\frac{1}{2}{e}_{i-1}^n{D}_h^{+}{\overline{e}}_{i-1}^n\right)(ϵ{\overline{e}}_i^n+{\mathcal{D}}_{\tau }^0{e}_i^n)\hfill \\ \hfill & \le 2\theta \parallel D_h^{+}{\overline{\mathbf{e}}}^n{\parallel }_h\left(ϵ\parallel {\overline{\mathbf{e}}}^n{\parallel }_h+{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h\right)+(1-\theta ){\parallel D_h^{+}{\overline{\mathbf{e}}}^n\parallel }_h\left(ϵ\parallel {\overline{\mathbf{e}}}^n{\parallel }_h+{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h\right)\hfill \\ \hfill & +(1-\theta )\alpha \left(\parallel {\overline{\mathbf{e}}}^n{\parallel }_h+{\parallel D_h^{+}{\overline{\mathbf{e}}}^n\parallel }_h\right)\left(\frac{1}{\alpha }\right)\left(ϵ\parallel {\overline{\mathbf{e}}}^n{\parallel }_h+{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h\right).\hfill \end{array}$$

Rewriting the above terms to arrive at

$$\begin{array}{cc}\hfill {\langle {\mathsf{\Psi }}_{\theta }({\mathbf{e}}^n,{\overline{\mathbf{e}}}^n),ϵ{\overline{\mathbf{e}}}^n+{\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\rangle }_h& \le \left(ϵ\theta +8{\theta }^2+ϵ(1-\theta )+8{(1-\theta )}^2\right){\parallel D_h^{+}{\overline{\mathbf{e}}}^n\parallel }_h^2\hfill \\ \hfill & +\left(ϵ\theta +2ϵ(1-\theta )+{(1-\theta )}^2\right)\parallel {\overline{\mathbf{e}}}^n{\parallel }_h^2+\frac{1}{2}{\parallel {\mathcal{D}}_{\tau }^0{\mathbf{e}}^n\parallel }_h^2.\hfill \end{array}$$

(58)

Substituting (53)–(58) into (52), we obtain

$$\begin{array}{cc}\hfill & \frac{ϵ}{2\tau }\left(\parallel {\mathbf{e}}^{n+1}{\parallel }_h^2-{\parallel {\mathbf{e}}^{n-1}\parallel }_h^2\right)+\frac{ϵ}{2\tau }\left(\parallel D_h^{+}{\mathbf{e}}^{n+1}{\parallel }_h^2+{\parallel D_h^{+}{\mathbf{e}}^{n-1}\parallel }_h^2\right)\hfill \\ \hfill & \le c_6\left(\parallel {\mathbf{e}}^{n+1}{\parallel }_h^2+2\parallel {\mathbf{e}}^n{\parallel }_h^2+{\parallel {\mathbf{e}}^{n-1}\parallel }_h^2\right)\hfill \\ \hfill & +c_6\left(\parallel D_h^{+}{\mathbf{e}}^{n+1}{\parallel }_h^2+2\parallel D_h^{+}{\mathbf{e}}^n{\parallel }_h^2+{\parallel D_h^{+}{\mathbf{e}}^{n-1}\parallel }_h^2\right)+c_7{({\tau }^2+h^2)}^2.\hfill \end{array}$$

(59)

Define

$$S^n=(\parallel {\mathbf{e}}^n{\parallel }_h^2+\parallel D_h^{+}{\mathbf{e}}^n{\parallel }_h^2)+(\parallel {\mathbf{e}}^{n+1}{\parallel }_h^2+\parallel D_h^{+}{\mathbf{e}}^{n+1}{\parallel }_h^2).$$

One can show that (59) leads to

$$\begin{array}{cc}\hfill \frac{ϵ}{2\tau }(S^n-S^{n-1})& \le c_6(S^n+S^{n-1})+c_7{({\tau }^2+h^2)}^2\hfill \\ \hfill \left(1-\frac{2c_6}{ϵ}\tau \right)S^n& \le \left(1+\frac{2c_6}{ϵ}\tau \right)S^{n-1}+\frac{2c_7}{ϵ}\tau {({\tau }^2+h^2)}^2.\hfill \end{array}$$

For $\frac{2c_6}{ϵ}\tau \le \frac{1}{2}$, we arrive at

$$\begin{array}{c}\hfill S^n\le \left(1+\frac{4c_6}{ϵ}\tau \right)S^{n-1}+\frac{4c_7}{ϵ}\tau {({\tau }^2+h^2)}^2.\end{array}$$

(60)

A use of discrete Grönwall inequality shows

$$\begin{array}{cc}\hfill S^n& \le exp\left(\frac{4c_6}{ϵ}n\tau \right)\left(S^0+\frac{4c_7}{ϵ}n\tau {({\tau }^2+h^2)}^2\right)\hfill \\ \hfill & \le exp\left(\frac{4c_{6}T}{ϵ}\right)\left(S^0+\frac{4c_{7}T}{ϵ}{({\tau }^2+h^2)}^2\right).\hfill \end{array}$$

(61)

Take the value of $S_0$ from (50), we derive

$$\begin{array}{cc}\hfill S^n& \le exp\left(\frac{4c_{6}T}{ϵ}\right)\left(2\frac{c_5}{ϵ}{({\tau }^{3/2}+h^2)}^2+\frac{4c_{7}T}{ϵ}{({\tau }^2+h^2)}^2\right)\hfill \\ \hfill & \le \frac{c_8^2}{{\left(2\alpha \right)}^2}{({\tau }^{3/2}+h^2)}^2.\hfill \end{array}$$

(62)

We have from (62) that

$$\parallel {\mathbf{e}}^{n+1}{\parallel }_h\le \frac{c_8}{2\alpha }({\tau }^{3/2}+h^2)\mathrm{and}{\parallel D_h^{+}{\mathbf{e}}^{n+1}\parallel }_h\le \frac{c_8}{2\alpha }({\tau }^{3/2}+h^2).$$

Using Lemma 2, we get $\parallel {\mathbf{e}}^{n+1}{\parallel }_{h,\infty }\le c_8({\tau }^{3/2}+h^2)$ as needed. □

Remark 2. 

1. 

We are losing $1/2$ power in τ because of the initial approximation step. However, using a predictor–corrector method, one may obtain $\mathcal{O}\left({\tau }^2\right)$ in time.

2. 

As a consequence, we obtain the estimates ${\parallel e\parallel }_h$ and ${\parallel e\parallel }_{h,\infty }$.

5. Numerical Results

In this section, we verified the invariant conservation property and order of convergence through various numerical test problems. We tested the proposed scheme on the inviscous problem using exact solutions found in [3]. As for the viscous problem $(ϵ\ne 0)$, the test problem was taken from [4].

In Theorem 4, the optimal value of $\theta$ is $1/3$. The simulations using $\theta =1/3$ was conducted against other values of $\theta$ to compare their performances. The errors were measured by the discrete Euclidean norm ${\parallel ·\parallel }_h$ and the discrete uniform norm ${\parallel ·\parallel }_{h,\infty }$ defined in (13) and (14), respectively. The order of accuracy r with respect to the norm $\parallel ·\parallel$ was computed using the formula

$$r={log}_2\left(\frac{\parallel {\mathbf{e}}_1\parallel }{\parallel {\mathbf{e}}_2\parallel }\right)$$

where ${\mathbf{e}}_1$ is the errors at nodal points resulting from the step size twice larger than that of ${\mathbf{e}}_2$.

5.1. Accuracy Test for the Inviscous Problem

The simulation was conducted on the steady solution, which is periodic on the interval $[-p,p]$, with initial data

$$\begin{array}{cc}\hfill u_0\left(x\right)& =\frac{4}{3}\left(\frac{cosh(x/2)}{cosh(p/2)}-1\right).\hfill \end{array}$$

(63)

Using domain $p=2$, final time $T=1$, and $\tau =h$, we obtained an optimal order of convergence as shown in

Table 1. Errors and orders showing optimal convergence rates for the inviscous problem, with $\theta =0,\frac{1}{3}$.

N	$\mathit{\theta }=0$				$\mathit{\theta }=1/3$
	${\parallel \mathbf{e}\parallel }_{\mathit{h}}$		${\parallel \mathbf{e}\parallel }_{\mathit{h},\infty }$		${\parallel \mathbf{e}\parallel }_{\mathit{h}}$		${\parallel \mathbf{e}\parallel }_{\mathit{h},\infty }$
	Error	Order	Error	Order	Error	Order	Error	Order
40	6.7498 × 10−4		7.8284 × 10−4		4.2413 × 10−4		5.4619 × 10−4
80	1.7452 × 10−4	1.95	2.0982 × 10−4	1.90	1.1021 × 10−4	1.94	1.4711 × 10−4	1.89
160	4.4369 × 10−5	1.98	5.4335 × 10−5	1.95	2.8086 × 10−5	1.97	3.8179 × 10−5	1.95
320	1.1186 × 10−5	1.99	1.3826 × 10−5	1.97	7.0886 × 10−6	1.99	9.7254 × 10−6	1.97
640	2.8083 × 10−6	1.99	3.4874 × 10−6	1.99	1.7806 × 10−6	1.99	2.4543 × 10−6	1.99

and

Table 2. Errors and orders showing optimal convergence rates for the inviscous problem, with $\theta =\frac{2}{3},1$.

N	$\mathit{\theta }=2/3$				$\mathit{\theta }=1$
	${\parallel \mathbf{e}\parallel }_{\mathit{h}}$		${\parallel \mathbf{e}\parallel }_{\mathit{h},\infty }$		${\parallel \mathbf{e}\parallel }_{\mathit{h}}$		${\parallel \mathbf{e}\parallel }_{\mathit{h},\infty }$
	Error	Order	Error	Order	Error	Order	Error	Order
40	1.7901 × 10−4		2.7492 × 10−4		9.8426 × 10−5		7.4600 × 10−5
80	4.7136 × 10−5	1.93	7.6875 × 10−5	1.84	2.4617 × 10−5	2.00	1.8652 × 10−5	2.00
160	1.2084 × 10−5	1.96	2.0333 × 10−5	1.92	6.1550 × 10−6	2.00	4.6661 × 10−6	2.00
320	3.0585 × 10−6	1.98	5.2290 × 10−6	1.96	1.5388 × 10−6	2.00	1.1665 × 10−6	2.00
640	7.6934 × 10−7	1.99	1.3259 × 10−6	1.98	3.8470 × 10−7	2.00	2.9163 × 10−7	2.00

. The performances for each value of $\theta$ are shown in Figure 1 and Figure 2. We see from Figure 1 that all the errors are vanishing for any value of $\theta$. However, as shown in Figure 2, the quantity $\parallel u^n{\parallel }_h$ is stable only when $\theta =1/3$.

5.2. Accuracy Test for the Viscous Problem

For the case $ϵ\ne 0$, we used the non-homogeneous example with $\varphi (x,t)=sech(x-t)$ provided in [4]. The simulation was conducted on the domain $[-20,20]$ with the final time $T=1$ and $\tau =h$. The errors and orders of accuracy are given in

Table 3. Errors and orders showing optimal convergence rates for the viscous problem, with $\theta =0,\frac{1}{3}$.

N	$\mathit{\theta }=0$				$\mathit{\theta }=1/3$
	${\parallel \mathbf{e}\parallel }_{\mathit{h}}$		${\parallel \mathbf{e}\parallel }_{\mathit{h},\infty }$		${\parallel \mathbf{e}\parallel }_{\mathit{h}}$		${\parallel \mathbf{e}\parallel }_{\mathit{h},\infty }$
	Error	Order	Error	Order	Error	Order	Error	Order
80	9.8070 × 10−1		8.8961 × 10−1		1.0900		1.0686
160	1.5003 × 10−1	2.71	1.6312 × 10−1	2.45	1.4306 × 10−1	2.93	1.5557 × 10−1	2.78
320	3.0028 × 10−2	2.32	3.0737 × 10−2	2.41	2.9229 × 10−2	2.29	3.1383 × 10−2	2.31
640	7.0535 × 10−3	2.09	7.3580 × 10−3	2.06	7.0430 × 10−3	2.05	7.6141 × 10−3	2.04
1280	1.7182 × 10−3	2.04	1.8173 × 10−3	2.02	1.7445 × 10−3	2.01	1.8972 × 10−3	2.00

and

Table 4. Errors and orders showing optimal convergence rates for the viscous problem, with $\theta =\frac{2}{3},1$.

N	$\mathit{\theta }=2/3$				$\mathit{\theta }=1$
	${\parallel \mathbf{e}\parallel }_{\mathit{h}}$		${\parallel \mathbf{e}\parallel }_{\mathit{h},\infty }$		${\parallel \mathbf{e}\parallel }_{\mathit{h}}$		${\parallel \mathbf{e}\parallel }_{\mathit{h},\infty }$
	Error	Order	Error	Order	Error	Order	Error	Order
80	1.1878		1.2909		1.2501		1.4776
160	1.5447 × 10−1	2.94	1.7385 × 10−1	2.89	1.7512 × 10−1	2.84	2.0320 × 10−1	2.86
320	3.3077 × 10−2	2.22	3.5948 × 10−2	2.27	3.7731 × 10−2	2.21	4.1052 × 10−2	2.31
640	8.0712 × 10−3	2.03	8.7244 × 10−3	2.04	9.2118 × 10−3	2.03	9.9053 × 10−3	2.05
1280	2.0101 × 10−3	2.01	2.1747 × 10−3	2.00	2.2928 × 10−3	2.01	2.4630 × 10−3	2.01

. Note that the optimal order of convergence can be achieved for the viscous problem as well. In Figure 3, all errors are vanishing for any value of $\theta$.

From the results in Section 5.1 and Section 5.2, we see that the scheme converges uniformly for any value of $\theta$. However, we see in Figure 2 that the quantity $\parallel u^n{\parallel }_h$ is preserved only when $\theta =1/3$. This agrees with the result from Theorem 4. In the subsequent experiment, we show that the scheme also performs better in a long-time simulation when $\theta =1/3$.

5.3. Invariant-Preserving Test

To see how well the scheme performs on a long-time simulation, we tested it on the exact solution with cusp:

$$\begin{array}{cc}\hfill u(x,t)& =\frac{4}{3}\left(e^{\left|x\right|/2}-1\right),\hfill \end{array}$$

(64)

and observed how the sharp profile changed over time. The simulation was done using $M=640$ on the domain $[-20,20]$ with time step $\mathsf{\Delta }t={10}^{-4}$. In Figure 4, we see that when $\theta =1/3$ the scheme preserves the shape of numerical solution, but when we used $\theta =0$, the numerical solution becomes unstable as time increases. This shows that the invariant-preserving method performs better in the long run.

5.4. Asymptotic Test

We tested the invariant-preserving scheme ($\theta =1/3$) on the initial data

$$u_0\left(x\right)=\left\{\begin{array}{cc}0.5,\hfill & 0\le x\le 1,\hfill \\ -x+1.5,\hfill & 1\le x\le 4,\hfill \\ -2.5,\hfill & x\ge 4.\hfill \end{array}\right.$$

The solution is expected to converge to a stationary solution. In Figure 5, we observe a stable pattern formation as analyzed in [3]. Here, we used $h=1/16$.

6. Discussion

Results from Section 5 verified that the numerical solutions converge to the exact solutions with second order of accuracy with respect to the spatial variable. The convergence can be inferred from the error plots in Figure 1 and Figure 3 while the rate of convergence can be inferred from the numeric evidence in Table 1, Table 2, Table 3 and Table 4. Note that the second-order convergence occurs at any value of $\theta$. However, Figure 2 shows that the invariant is preserved only when $\theta =1/3$. This behavior agrees with the results we proved in Section 3. The comparison in Section 5.3 confirmed that the invariant-preserving scheme performs better in the long run. This hints that the invariant-preserving property, if any, should be taken into account when one wants to design a numerical method.

The experiment on the asymptotic problem in Section 5.4 ensured that the proposed scheme is robust and can be applied to a broader range of problems. Possible applications include but are not limited to adaption to other types of boundary conditions, changing to various types of initial data, and adding some extra terms to the equation. All of these possibilities allow us to simulate wider ranges of natural phenomena. On the other hand, one can improve the accuracy of presented. In [19], a fourth-order operator is applied to a viscous Burgers’ equation. The same operator can be applied to the Burgers-Poisson system to obtain a fourth-order method as well. Besides this, the idea of the proposed method can also be generalized to higher dimension. While the simplicity of the finite diffidence method makes the generalization of the implementation straightforward, the analysis part may not be so. One may opt for other approach such as the spectral volume (SV) method [20,21,22] which has high order of accuracy and can be readily applied to higher dimensional settings.

7. Conclusions

In this paper, we proposed a family of implicit finite difference $\theta$-schemes to solve the system of Burgers-Poisson equations. The method used information from three time levels and eliminated the need to solve unknowns using iterative nonlinear solvers. When $\theta =1/3$, we showed that the resulting scheme has a discrete invariant-preserving property. We also proved that the numerical solution converges uniformly and has second order of accuracy with respect to the spatial variable. Finally, the numerical experiments verified that the optimal order of convergence is achieved for any value of $\theta$, but the invariant is stable only when $\theta =1/3$. All the numerical evidences obtained from the simulations agree with the theoretical results.