Numerical Solutions of the Improved Boussinesq Equation with Stokes Damping

Quanxiang Wang

College of Science, Nanjing Agricultural University, Nanjing, China.

DOI: 10.4236/am.2021.124015                

Abstract

In this paper, we propose new finite volume element schemes to numerically solve the improved Boussinesq equation with Stokes damping. The new schemes can inherit characteristic properties of the conservation of mass and the decrease of total energy from the improved Boussinesq equation with Stokes damping. Numerical experiments illustrate that the proposed schemes are second-order accuracy in space and time.

Keywords

Energy Dissipation, Improved Boussinesq Equation, Finite Volume Element, Stokes Damping

1. Introduction

The Boussinesq equation was first presented in the 1870s by Joseph Boussinesq to model the propagation of shallow water waves in multiple directions. It is a fourth-order nonlinear partial differential equation and belongs to KdV family, which has been widely used in math-physical field related to nonlinear wave phenomena, such as ion-sound in plasma, nonlinear lattice waves [1]. The Boussinesq equation is also used to study the disturbance and repair of solitary waves in blood vessels with aneurysm [2] [3]. A general form of the Boussinesq equation is

$u_{tt}=u_{xx}+q{u}_{xxxx}+{\left(u^2\right)}_{xx}\mathrm{,}$ (1)

where $u=u\left(x\mathrm{,}t\right)$ is a sufficiently differentiable function, which for $q=-1$ gives the good Boussinesq or well-posed equation, while for $q=1$ the bad or ill-posed Boussinesq equation.

In this study, we will take $q=1$. If the term $u_{xxxx}$ in Equation (1) is replaced with $u_{xxtt}$, it gives the so-called improved Boussinesq equation

$u_{tt}=u_{xx}+u_{xxtt}+{\left(u^2\right)}_{xx}.$ (2)

Researchers do not observe any instability for the above improved Boussinesq equation [4]. The improved Boussinesq equation (IBq) has been investigated theoretically and numerically by many people.

Yang [5] studied the existence, both locally and globally in time, the uniqueness of solutions and non-existence of global solutions to the initial boundary value problem of a generalized IBq equation. Iskandar and Jain [6] studied the dynamical behavior of IBq equation numerically. El-Zoheiry [7] designed a three-level iterative scheme based on the compact implicit method for solving the improved Boussinesq equation. Bratsos [8] proposed a second order implicit finite difference scheme by rational approximations and constructed a predictor-corrector scheme [9] to avoid the solution of a nonlinear system. Lin et al. [10] used the finite element method with linear B-spline basis functions to discretize the improved Boussinesq equation in space and derived a second-order system. Shokri and Dehghan [11] used collocation and approximated the solution by radial basis functions based on the third-order time discretization to solve the equation. Irk and Dăg [12] obtained numerical simulations of the improved Boussinesq equation by using two finite difference schemes and two finite element methods, based on the second- and third-order time discretization. Zhang and Lu [13] derived quadratic finite volume element schemes to solve the improved Boussinesq equation.

So far the main attention was paid to the improved Boussinesq equation in the absence of damping. However, the damping will change the amplitude of the waves and their shape. Yan et al. [14] studied the high-order modified Boussinesq equation with damping term. In [15], Arevalo et al. investigated the dynamics of a lattice soliton in the presence of damping. By using a multiple-scale perturbation expansion up to second-order, they derived a general expression for the first-order velocity correction. Chen et al. [16] studied the Cauchy problem for the generalized IBq equation with hydrodynamical damping.

The improved Boussinesq equation of interest is here

$u_{tt}=u_{xx}+u_{xxtt}+{\left(u^2\right)}_{xx}-\gamma u_t,\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}x\in \left[a,b\right]$ (3)

where $\gamma$ is the positive constant. According to [15], the damping in the equation corresponds to the Stokes damping. For the improved Boussinesq equation with Stokes damping, we adopt the following periodic boundary conditions

$\frac{{\partial }^l}{\partial x^l}{\frac{{\partial }^{j}u}{\partial t^j}|}_{x=a}=\frac{{\partial }^l}{\partial x^l}{\frac{{\partial }^{j}u}{\partial t^j}|}_{x=b},\left(l=0,1;j=0,1\right).$ (4)

As the analytic solutions for the Equation (3) are not available, it is necessary to design efficient numerical schemes to perform numerical simulation. To investigate the influence of damping on the solution, the proposed schemes should inherit characteristic properties of the system. The finite volume element method has been one of the most commonly used numerical methods for solving partial differential equations [17]. One main attractive property of the method is that, the mass conservation law is maintained, which is fairly desirable for fluid and underground fluid computations. So it has been extensively used in computational fluid dynamics [18] - [23].

In this study, we will propose new finite volume element schemes for the improved Boussinesq equation with Stokes damping. The schemes can preserve mass conservation and energy dissipation of the improved Boussinesq equation. Furthermore, we will study numerically the rate of decay of solutions of the considered equation using our numerical schemes.

The paper is organized as follows: Our numerical method is described in Section 2, which will give good results for deterministic cases. In Section 3, we present some numerical examples to illustrate the accuracy and conservation property of the proposed numerical schemes. In Section 4, the numerical results for the improved Boussinesq equation with damping are given. The conclusions are summarized in the last section.

2. Numerical Method

In this section, new finite volume element schemes will be derived to solve the considered Equation (3).

2.1. Energy Dissipation Property

In order to see the energy dissipation property, it is more convenient to work with the following system representation

$\{\begin{array}{l}{u}_t=v_x\mathrm{,}\\ v_t-v_{xxt}=\frac{\partial }{\partial x}\left(\frac{\delta G}{\delta u}\right)-\gamma v\mathrm{,}\end{array}$ (5)

where $\frac{\delta G}{\delta u}$ is the Euler-Lagrange variational derivative [24] [25] defined by

$\frac{\delta G}{\delta u}=\frac{\partial G}{\partial u}-\frac{\partial }{\partial x}\left(\frac{\delta G}{\delta u_x}\right),\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}G=\frac{u^2}{2}+\frac{u^3}{3}.$

Then, the improved Boussinesq equation with Stokes damping has the following mass conservation property.

Proposition 1. Let u be the solution of Equation (5). Then the total mass ${\displaystyle {\int }_a^b}u\text{d}x$ is independent of t. Namely,

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_a^b}u\text{d}x=0.$ (6)

The conservation of mass can be easily proved with periodic boundary conditions.

Proposition 2. Let u be the solution of Equation (5). Then the total energy ${\displaystyle {\int }_a^b}\tilde{G}\left(u\mathrm{,}{u}_t\mathrm{,}v\right)\text{d}x$ is independent of t. Namely,

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_a^b}\tilde{G}\left(u\mathrm{,}{u}_t\mathrm{,}v\right)\text{d}x\le \mathrm{0,}$ (7)

where

$\tilde{G}\left(u\mathrm{,}{u}_t\mathrm{,}v\right)=\frac{v^2}{2}+\frac{{\left(u_t\right)}^2}{2}+G\left(u\right)\mathrm{.}$

Proof. According to the chain rule, we have

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_a^b}\tilde{G}\text{d}x={\displaystyle {\int }_a^b}\left(v{v}_t+u_t{u}_{tt}+\frac{\delta G}{\delta u}{u}_t\right)\text{d}x$

By the first equality of Equation (5), we can get

$={\displaystyle {\int }_a^b}\left(v{v}_t+v_x{v}_{xt}+\frac{\delta G}{\delta u}{v}_x\right)\text{d}x$

Using the integration by part formula, we have

$={\displaystyle {\int }_a^b}v\left(v_t-v_{xxt}-\frac{\partial }{\partial x}\left(\frac{\delta G}{\delta u}\right)\right)\text{d}x$

At last, according to the second equality of Equation (5), we can complete the proof

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_a^b}\tilde{G}dx=-\gamma {\displaystyle {\int }_a^b}{v}^{\text{2}}\text{d}x\le 0.$

2.2. Trial and Test Function Spaces

The region of interest $\left[a\mathrm{,}b\right]$ can be decomposed into a grid $T_h$ with nodes

$a=x_0<x_1<x_2<\cdots <x_{N-1}<x_N=b.$ (8)

where $x_i-x_{i-1}=h$ and $h=\left(b-a\right)/N$. Denote $T_h=\left\{I_i:I_i=\left[x_{i-1},x_i\right],i=1,2,\cdots ,N\right\}$. Accordingly, we place a dual grid $T_h^{\mathrm{<mo>\; *\; </mo>}}$ with nodes

$a=x_0<x_{1/2}<x_{3/2}<\cdots <x_{N-3/2}<x_{N-1/2}<x_N=b,$ (9)

where $x_{i-1/2}=x_i-h/2,i=1,2,\cdots ,N$ Denote $I_0^{\ast }=\left[x_0,x_{1/2}\right]$, $I_i^{\ast }=\left[x_{i-1/2}\mathrm{,}{x}_{i+1/2}\right]$ and $I_N^{\ast }=\left[x_{N-1/2}\mathrm{,}{x}_N\right]$.

Select the trial function space $U_h$ as the linear element space with respect to $T_h$. The basis function of the node $x_i$ is

${\varphi }_i\left(x\right)=\{\begin{array}{l}1-\left|x-x_i\right|/h,x_{i-1}\le x\le x_{i+1},\\ 0,\text{elsewhere}.\end{array}$ (10)

Then numerical solution U for Equation (5) can be uniquely written as $U={\displaystyle \underset{i=1}{\overset{N}{\sum }}}{U}_i{\varphi }_i\left(x\right)$, where $U_i=U\left(x_i\mathrm{,}t\right)$. So in the element $I_i$, we have

$U=U_{i-1}\left(1-\mu \right)+U_i\mu ,$ (11)

$U^{\prime }=\left(U_i-U_{i-1}\right)/h,$ (12)

where $\mu =\left(x-x_{i-1}\right)/h$.

The test function space $V_h$ corresponding to $T_h^{\mathrm{<mo>\; *\; </mo>}}$ is taken as the piecewise constant function space. The test function of the nodes $x_j$ is

${\psi }_j\left(x\right)=\{\begin{array}{l}\mathrm{1,}{x}_{j-1/2}\le x\le x_{j+1/2}\mathrm{,}\\ \mathrm{0,}\text{elsewhere}\mathrm{,}\end{array}$ (13)

In addition, we take $\Delta t$ as the time step and the numerical solution can be denoted by $U_k^n\simeq u\left(kh\mathrm{,}n\Delta t\right)$.

Recall the definition $\tilde{G}\left(u\mathrm{,}{u}_t\mathrm{,}v\right)=\frac{v^2}{2}+\frac{{\left(u_t\right)}^2}{2}+G\left(u\right)$. According to the definition of Euler-Lagrange derivative, we have $\frac{\delta \tilde{G}}{\delta u}=\frac{\delta G}{\delta u}$, $\frac{\delta \tilde{G}}{\delta u_t}=u_t$ and $\frac{\delta \tilde{G}}{\delta v}=v$. So the Equation (5) can be rewritten as

$\{\begin{array}{l}{u}_t=\frac{\partial }{\partial x}\left(\frac{\delta \tilde{G}}{\delta v}\right)\mathrm{,}\\ \left(1-{\partial }_x^2\right)v_t=\frac{\partial }{\partial x}\left(\frac{\delta \tilde{G}}{\delta u}\right)-\gamma v\mathrm{,}\end{array}$ (14)

where ${\partial }_x^2=\frac{{\partial }^2}{\partial x^2}$.

Using any $\psi \in V_h$ and $\tilde{\psi }\in V_h$ to multiply both sides of Equation (14) and integrating on the interval $\left[a\mathrm{,}b\right]$, we can obtain the following semi-discrete scheme

$\{\begin{array}{l}\left(u_t\mathrm{,}\psi \right)=\left(\frac{\partial }{\partial x}\left(\frac{\delta \tilde{G}}{\delta v}\right)\mathrm{,}\psi \right)\mathrm{,}\\ \left(\left(1-{\partial }_x^2\right)v_t\mathrm{,}\tilde{\psi }\right)=\left(\frac{\partial }{\partial x}\left(\frac{\delta \tilde{G}}{\delta u}\right)-\gamma v\mathrm{,}\tilde{\psi }\right)\mathrm{,}\end{array}$ (15)

where $\left(\cdot \mathrm{,}\cdot \right)$ denotes the inner product of $L^2\left(\left[a\mathrm{,}b\right]\right)$.

Now we discretize the above semi-discrete scheme in time. To simplify the notation, we will use the difference operator ${\delta }_n^{\langle 1\rangle }{u}^n=\frac{u^{n+1}-u^{n-1}}{2\Delta t}$. Additionally, we define the discrete version of Euler-Lagrange derivative $\frac{\delta \tilde{G}}{\delta u}=\frac{u^2}{2}+\frac{u^3}{3}$, $\frac{\delta \tilde{G}}{\delta u_t}=u_t$ and $\frac{\delta \tilde{G}}{\delta v}=v$ by

$\frac{\delta {\tilde{G}}_d}{\delta \left(u^{n+1}\mathrm{,}{u}^n\mathrm{,}{u}^{n-1}\right)}=u^n+\frac{u^n\left(u^{n+1}+u^n+u^{n-1}\right)}{3}\mathrm{,}$ (16)

$\frac{\delta {\tilde{G}}_d}{\delta \left({\delta }_n^{\langle 1\rangle }{u}^n\right)}={\delta }_n^{\langle 1\rangle }{u}^n,$ (17)

$\frac{\delta {\tilde{G}}_d}{\delta \left(v^n\right)}=v^n\mathrm{.}$ (18)

At last, we can get the following energy dissipation finite volume element schemes

$\{\begin{array}{l}\left({\delta }_n^{\langle 1\rangle }{u}^n\mathrm{,}\psi \right)=\left(\frac{\partial }{\partial x}\left(\frac{\delta {\tilde{G}}_d}{\delta \left(v^n\right)}\right)\mathrm{,}\psi \right)\mathrm{,}\\ \left(\left(1-{\partial }_x^2\right){\delta }_n^{\langle 1\rangle }{v}^n\mathrm{,}\tilde{\psi }\right)=\left(\frac{\partial }{\partial x}\left(\frac{\delta {\tilde{G}}_d}{\delta \left(u^{n+1}\mathrm{,}{u}^n\mathrm{,}{u}^{n-1}\right)}\right)-\gamma v^n\mathrm{,}\tilde{\psi }\right)\mathrm{.}\end{array}$ (19)

Under the periodic boundary conditions, the Equation (19) has the following discrete mass conservation law and energy dissipation property.

Theorem 1. (Discrete mass conservation law). Under the discrete periodic boundary condition, then numerical solution by Equation (19) conserves mass

${\displaystyle {\int }_a^b}{u}^n\text{d}x=const\mathrm{.}$ (20)

Theorem 2. (Discrete energy dissipation property). Under the discrete periodic boundary condition, then numerical solution by Equation (19) conserves energy

$\frac{1}{\Delta t}\left({\displaystyle {\int }_a^b}{\tilde{G}}_d^{n+1}\text{d}x-{\displaystyle {\int }_a^b}{\tilde{G}}_d^n\text{d}x\right)\le 0.$ (21)

Proof.

$\begin{array}{l}\frac{1}{\Delta t}{\displaystyle {\int }_a^b}\left({\tilde{G}}_d^{n+1}-{\tilde{G}}_d^n\right)\text{d}x\\ =\left(\frac{\delta {\tilde{G}}_d}{\delta \left(u^{n+1},u^n,u^{n-1}\right)},{\delta }_n^{\langle 1\rangle }{u}^n\right)+\left(\frac{\delta {\tilde{G}}_d}{\delta \left({\delta }_n^{\langle 1\rangle }{u}^n\right)},\frac{{\delta }_n^{\langle 1\rangle }{u}^{n+1}-{\delta }_n^{\langle 1\rangle }{u}^{n-1}}{2\Delta t}\right)+\left(\frac{\delta {\tilde{G}}_d}{\delta \left(v^n\right)},{\delta }_n^{\langle 1\rangle }{v}^n\right)\\ =\left(\frac{\delta {\tilde{G}}_d}{\delta \left(u^{n+1},u^n,u^{n-1}\right)},\frac{\partial }{\partial x}\left(\frac{\delta {\tilde{G}}_d}{\delta \left(v^n\right)}\right)\right)+\left(\frac{\partial }{\partial x}\left(\frac{\delta {\tilde{G}}_d}{\delta \left(v^n\right)}\right),\frac{v_x^{n+1}-v_x^{n-1}}{2\Delta t}\right)+\left(\frac{\delta {\tilde{G}}_d}{\delta \left(v^n\right)},{\delta }_n^{\langle 1\rangle }{v}^n\right)\\ =-\left(\frac{\partial }{\partial x}\left(\frac{\delta {\tilde{G}}_d}{\delta \left(u^{n+1},u^n,u^{n-1}\right)}\right),\frac{\delta {\tilde{G}}_d}{\delta \left(v^n\right)}\right)+\left(\left(1-{\partial }_x^2\right){\delta }_n^{\langle 1\rangle }{v}^n,\frac{\delta {\tilde{G}}_d}{\delta \left(v^n\right)}\right)\\ =-\left(\gamma v^n,\frac{\delta {\tilde{G}}_d}{\delta \left(v^n\right)}\right)=-\gamma {\displaystyle {\int }_a^b}{\left(v^n\right)}^2\text{d}x\le 0\end{array}$

2.3. Finite Volume Element Schemes

For comparison, we employ the standard finite volume element scheme constructed based on the

$\{\begin{array}{l}{u}_t=v_x\mathrm{,}\\ \left(1-{\partial }_x^2\right)v_t={\left(u+u^2\right)}_x-\gamma v\mathrm{.}\end{array}$ (22)

Then we can obtain the following finite volume element schemes

$\{\begin{array}{l}\left({\delta }_n^{\langle 1\rangle }{u}^n\mathrm{,}\psi \right)=\left(v_x^n\mathrm{,}\psi \right)\mathrm{,}\\ \left(\left(1-{\partial }_x^2\right){\delta }_n^{\langle 1\rangle }{v}^n\mathrm{,}\tilde{\psi }\right)=\left({\left(u^n+{\left(u^n\right)}^2\right)}_x-\gamma v^n\mathrm{,}\tilde{\psi }\right)\mathrm{.}\end{array}$ (23)

Theorem 3. (Discrete mass conservation law). Under the discrete periodic boundary condition, then numerical solution by Equation (23) conserves mass

${\displaystyle {\int }_a^b}{u}^n\text{d}x=const\mathrm{.}$ (24)

3. Validation of the Numerical Schemes

In this section, we present some numerical examples to illustrate the accuracy and conservation property of the proposed numerical schemes.

3.1. Accuracy

The initial conditions are taken as

$u\left(x\mathrm{,0}\right)=\alpha {\text{sech}}^2\left(\sqrt{\frac{\alpha }{6}}\frac{x-x_0}{\beta }\right)\mathrm{,}$ (25)

$u_t\left(x,0\right)=2\alpha \sqrt{\frac{\alpha }{6}}{\text{sech}}^2\left(\sqrt{\frac{\alpha }{6}}\frac{x-x_0}{\beta }\right)\tanh \left(\sqrt{\frac{\alpha }{6}}\frac{x-x_0}{\beta }\right),$ (26)

where $\alpha$ is amplitude of the solitary wave which is initially centered at $x=x_0$ travelling with velocity $\beta$. In addition, we will take $x_0=0$, $\alpha =0.5$, $\gamma =0.04$ and the region $\left[-\mathrm{20,80}\right]$ in the following.

As the exact solution for the considered examples is not available, we will compute maximum absolute error for the considered examples using the half mesh principle. The maximum absolute error is defined as follows,

$E\left(h\mathrm{,}\Delta t\right)=\underset{i}{\max }\left|u_h^{\Delta t}\left(i\mathrm{,}n\right)-u_{h/2}^{\Delta t/2}\left(2i\mathrm{,2}n\right)\right|\mathrm{.}$ (27)

The maximum absolute errors Equation (3) and solitary wave with the above parameters have been computed at t= 1. They are displayed in Table 1 for the new finite volume element schemes (NFVEM). Examination of the table shows that the error measures of the finite volume element scheme diminish approximately quadratically as the space step size and time step size are simultaneously halved. The error table illustrates that the new numerical schemes have second-order accuracy in space and time direction.

Furthermore, we compute the problem using standard finite volume element schemes (FVEM). The obtained maximum absolute errors are displayed in Table 2. From the table, we can find that the standard finite volume element schemes also have second-order accuracy in space and time direction. The maximum absolute errors computed by the two numerical schemes are consistent. In other words, the proposed new schemes do not reduce the accuracy of the finite volume method.

3.2. Conservation Properties

Now we consider the conservation properties of the new finite volume element schemes and standard finite volume element schemes. The initial conditions and parameters are the same as the above except $\gamma$.

Figure 1 shows the evolutions of discrete mass of the proposed new finite volume element schemes and standard finite volume element schemes with $\gamma ={10}^{-4}$.

Table 1. The maximum absolute error for the NFVEM at t = 1.

Table 2. The maximum absolute error for the FVEM at t = 1.

Figure 1. The evolutions of the discrete mass: (a) FVEM; (b) NFVEM. $h=2\Delta t=0.1$.

Figure 2. The evolutions of the discrete energy: (a) FVEM; (b) NFVEM.

From the figure, we can find that the discrete mass is well conserved. Actually it is conserved up to the machine accuracy, which well agrees with the discrete mass conservation law. In other words, the Stokes damping will not change the mass of system.

Figure 2 presents the evolutions of discrete energy of the proposed new finite volume element schemes and standard finite volume element schemes with $\gamma ={10}^{-4}$. It is easy to find that the discrete energy obtained by standard finite volume element schemes decreases with oscillations. But the discrete energy obtained by new finite volume element schemes will decrease continuously. So our new finite volume element schemes can preserve the energy dissipation property of the improved Boussinesq equation with Stokes damping. From this section, we can conclude that the proposed numerical schemes are more suitable for solving improved Boussinesq equation with Stokes damping.

4. Conclusion

In this paper, we have proposed a new finite volume element method for solving the improved Boussinesq equation with Stokes damping. The proposed method is second-order accurate in time and space. Furthermore, the new method can preserve the mass conservation and energy dissipation property.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Bogolubsky, I. (1977) Some Examples of Inelastic Soliton Interaction. Computer Physics Communications, 13, 149-155. https://doi.org/10.1016/0010-4655(77)90009-1
[2]	Dehghan, M. and Shakeri, F. (2008) Use of He’s Homotopy Perturbation Method for Solving a Partial Differential Equation Arising in Modeling of Flow in Porous Media. Journal of Porous Media, 11, 765-778. https://doi.org/10.1016/j.cnsns.2009.02.021
[3]	Biswas, A., Milovic, D. and Ranasinghe, A. (2009) Solitary Waves of Boussinesq Equation in a Power Law Media. Communications in Nonlinear Science and Numerical Simulation, 14, 3738-3742. https://doi.org/10.1016/j.cnsns.2009.02.021
[4]	Christiansen, P.L., Muto, V. and Rionero, S. (1992) Solitary Wave Solutions to a System of Boussinesq-Like Equations. Chaos, Solitons and Fractals, 2, 45-50. https://doi.org/10.1016/0960-0779(92)90046-P
[5]	Yang, Z. (1998) Existence and Non-Existence of Global Solutions to a Generalized Modification of the Improved Boussinesq Equation. Mathematical Methods in the Applied Sciences, 21, 1467-1477. https://doi.org/10.1002/(SICI)1099-1476(19981110)21:16%3C1467::AID-MMA968%3E3.0.CO; 2-K
[6]	Iskandar, L. and Jain, P.C. (1980) Numerical Solutions of the Improved Boussinesq Equation. Proceedings of the Indian Academy of Sciences-Mathematical Sciences, 89, 171-181. https://doi.org/10.1007/BF02861996
[7]	El-Zoheiry, H. (2002) Numerical Study of the Improved Boussinesq Equation. Chaos, Solitons and Fractals, 14, 377-384. https://doi.org/10.1016/S0960-0779(00)00271-X
[8]	Bratsos, A.G. (2007) A Second-Order Numerical Scheme for the Improved Boussinesq Equation. Physics Letters A, 370, 145-147. https://doi.org/10.1016/j.physleta.2007.05.050
[9]	Bratsos, A.G. (2009) A Predictor-Corrector Scheme for the Improved Boussinesq Equation. Chaos, Solitons & Fractals, 40, 2083-2094. https://doi.org/10.1016/j.chaos.2007.09.083
[10]	Lin, Q., Wu, Y., Loxton, R. and Lai, S. (2009) Linear B-Spline Finite Element Method for the Improved Boussinesq Equation. Journal of Computational and Applied Mathematics, 224, 658-667. https://doi.org/10.1016/j.cam.2008.05.049
[11]	Shokri, A. and Dehghan, M. (2010) A Not-a-Knot Meshless Method Using Radial Basis Functions and Predictor—Corrector Scheme to the Numerical Solution of Improved Boussinesq Equation. Computer Physics Communications, 181, 1990-2000. https://doi.org/10.1016/j.cpc.2010.08.035
[12]	Irk, D. and Dağ, İ. (2009) Numerical Simulations of the Improved Boussinesq Equation. Numerical Methods for Partial Differential Equations, 26, 1316-1327. https://doi.org/10.1002/num.20492
[13]	Zhang, Z. and Lu, F. (2012) Quadratic Finite Volume Element Method for the Improved Boussinesq Equation. Journal of Mathematical Physics, 53, Article ID: 013505. https://doi.org/10.1063/1.3672197
[14]	Yan, Z., Xie, F. and Zhang, H. (2001) Symmetry Reduction, Integrability and Solitary Wave Solutions to High-Order Modified Boussinesq Equation with Damping term. Communications in Theoretical Physics, 36, 1-6. https://doi.org/10.1088/0253-6102/36/1/1
[15]	Arevalo, E., Gaididei, Y. and Mertens, F.G. (2002) Soliton Dynamics in Damped and Forced Boussinesq Equations. The European Physical Journal B, 27, 63-74.
[16]	Chen, G., Rui, W. and Chen, X. (2011) Cauchy Problem for a Damped Generalized IMBq Equation. Journal of Mathematical Physics, 52, Article ID: 053504. https://doi.org/10.1063/1.3577956
[17]	Li, R., Chen, Z. and Wu, W. (1999) Generalized Difference Methods for Differential Equation: Numerical Analysis of Finite Volume Methods. Marcel Dekker, New York.
[18]	Zhang, Z. (2008) Error Estimates of Finite Volume Element Method for the Pollution in Groundwater Flow. Numerical Methods for Partial Differential Equations, 25, 259-274. https://doi.org/10.1002/num.20340
[19]	Li, J., Chen, Z. and He, Y. (2012) A Stabilized Multi-Level Method for Non-Singular Finite Volume Solutions of the Stationary 3D Navier-Stokes Equations. Numerische Mathematik, 122, 279-304. https://doi.org/10.1007/s00211-012-0462-z
[20]	Luo, Z., Li, H., Sun, P., An, J. and Navon, I. (2013) A Reduced-Order Finite Volume Element Formulation Based on POD Method and Numerical Simulation for Two-Dimensional Solute Transport Problems. Mathematics and Computers in Simulation, 89, 50-68. https://doi.org/10.1016/j.matcom.2012.11.012
[21]	Wang, Q., Zhang, Z. and Li, Z. (2013) A Fourier Finite Volume Element Method for Solving Two-Dimensional Quasi-Geostrophic Equations on a Sphere. Applied Numerical Mathematics, 71, 1-13. https://doi.org/10.1016/j.apnum.2013.03.007
[22]	Wang, Q. and Zhang, Z. (2014) High-Order Upwind Finite Volume Element Schemes for Modeling of Neuronal Firing. International Journal of Computer Mathematics, 91, 625-640. https://doi.org/10.1080/00207160.2013.801463
[23]	Lin, Z., Ricardo, R. and Tian, C. (2014) Finite Volume Element Approximation of an Inhomogeneous Brusselator Model with Cross-Diffusion. Journal of Computational Physics, 256, 806-823. https://doi.org/10.1016/j.jcp.2013.09.009
[24]	Furihata, D. (2001) A Stable and Conservative Finite Difference Scheme for the Cahn-Hilliard Equation. Numerische Mathematik, 87, 675-699. https://doi.org/10.1007/PL00005429
[25]	Matsuo, T. (2007) New Conservative Schemes with Discrete Variational Derivatives for Nonlinear Wave Equations. Journal of Computational and Applied Mathematics, 203, 32-56. https://doi.org/10.1016/j.cam.2006.03.009