Boundary Element Method for the Mixed BBM-KdV Equation Compared to Non Standard Boundary Conditions

Abstract

In this chapter, we are interested in the numerical resolution of the mixed BBM-KdV equation defined in unbounded domain. Boundary Element Method (BEM) are introduced to truncate the equation into a considered bounded domain. BEM uses domain decomposition techniques to construct Boundary Condition (BC) as transmission between the bounded domain and its complementary. We then present a suitable approximation of these BC using Discrete Galerkin Method. Numerical simulations are made to show the efficiency of these BC. We also compare with another method that truncates the equation from unbounded to bounded domain, called Non Standard Boundary Conditions (NSBC) which introduces new variables to catch information at the boundary and compose a system to connect all these variables in the bounded domain. Further discussions are made to highlight the advantages of each method as well as the difficulties encountered in the numerical resolution.

1. Introduction

We consider a combination of two linearized typical dispersive partial differential equations that model solitary waves and all interactions between them, given as follows

$$\left\{\begin{array}{lll}{\partial }_{t}u\left(t,x\right)+\alpha {\partial }_{\mathit{xxx}}^{3}u\left(t,x\right)-\beta {\partial }_{\mathit{txx}}^{3}u\left(t,x\right)+\gamma {\partial }_{x}u\left(t,x\right)=& 0& \left(t,x\right)\in {\mathbb{R}}_{+}^{\ast }\times \mathbb{R}\\ u\left(0,x\right)=& u_0\left(x\right)& x\in \mathbb{R}\\ \underset{\mid x\mid \to \pm \infty }{lim}u\left(t,x\right)=& 0& t\in {\mathbb{R}}_{+}^{\ast }\end{array}\right.$$E1

such that $\alpha$, $\beta$ are dispersion parameters and are positive numbers, while $\gamma \in \mathbb{R}$ is the velocity number. In the case $\alpha =0$, we obtain the BBM equation [1] and when $\beta =0$, we get the KdV equation [2]. Our main purpose is to obtain numerical approximation of Eq. (1) when taken in a bounded domain $\left[0,T\right]\times \left[a,b\right]$ with suitable boundary conditions with no spurious reflections. For this regard, we use two different techniques that are BEM and NSBC.

The Boundary Element Method (BEM), also known as the Boundary Integral Equation Method (BIEM), is an alternative deterministic method that incorporates a mesh located, only, at domain boundaries and therefore attractive for free surface problems. There are two types of BEM, the direct BEM which requires a closed boundary so that the physical variables (e.g. pressure and normal velocity in acoustics) can only be considered from one side of the surface (interior or exterior), while the indirect (IBEM) can consider both sides of the surface and does not need a closed surface. In the first part of this chapter, we use this technique of BEM to derive the BC to the Eq. (1) in the domain [0, T] × [a, b]. More precisely, we are going to introduce the BEM to establish BC satisfied by the Eq. (1) on two interface points a and b by solving the same equation in the complementary domain $\mathbb{R}\backslash \left[a,b\right]$. The BEM has significant advantages over the finite element or difference methods (FEM or FDM), as there is no need for discretizing the domain $\mathbb{R}\backslash \left[a,b\right]$ into elements. It only uses infinite boundary condition and transmission condition to compute the solution at $a$ and $b$ as integral equations. Consequently, this integral equations will be fixed as the boundary conditions of the problem (1) on the bounded domain [0, T] × $\left[a,b\right]$. Therefore, the boundary condition are approximated as Fredholm Integral Equations of second kind.

Despite the meshing effort is limited and the system matrices are smaller, the BEM also has disadvantages over the Finite Element Method or Difference Finite Method. In fact, the BEM matrices are mostly populated with complex coefficients. Furthermore, singularities may arise in the solution. These deteriorate the efficiency of the solution and must be prevented [2].

The outline of this chapter is organized as follows. In section 2, we describe the BEM for the mixed BBM-KdV equation [3]. Next, we discuss the special case of the BBM equation and give the approximation of the resulting equation Finite Difference Method. Section 3 presents briefly another method to derive boundary conditions for BBM equation called NSBC introduced in [4]. Finally in section 4, comparison of both methods is given with numerical experiments to highlight the transparency of both BC obtained in sections 2 and 3.

2. Boundary element method for the mixed BBM-KdV equation

Being in one dimensional space, $\mathbb{R}$, the boundary of any bounded interval reduces to two points. Hence, we use the BEM to find two values that might depend on time. For this regard, we consider a bounded domain ${\mathrm{\Omega }}_T=]0,T[\times \mathrm{\Omega }$ where $\mathrm{\Omega }=]a,b[$ and $a,b,T\in \mathbb{R}\text{such that}a<b,T>0$. Note $\mathrm{\Sigma }=\left\{a,b\right\}$ and ${\mathrm{\Sigma }}_T=]0,T[\times \mathrm{\Sigma }$. we take the decomposition $\mathbb{R}={\mathrm{\Omega }}_g\cup \mathrm{\Omega }\cup {\mathrm{\Omega }}_d$, such that, ${\mathrm{\Omega }}_g=]-\infty ,a]$ and ${\mathrm{\Omega }}_d=[b,+\infty [$. The corresponding equations to (1) using Dirichlet-to-Neumann domain decomposition write

$$\left\{\begin{array}{lll}{\partial }_{t}u\left(t,x\right)+\alpha {\partial }_{\mathit{xxx}}^{3}u\left(t,x\right)-\beta {\partial }_{\mathit{txx}}^{3}u\left(t,x\right)+\gamma {\partial }_{x}u\left(t,x\right)& =0& \text{in}{\mathrm{\Omega }}_T\\ u\left(0,x\right)& =u_0\left(x\right)& \mathrm{at}\mathrm{\Omega }\\ {\partial }_{n}u& ={\partial }_{n}w& \mathrm{at}{\mathrm{\Sigma }}_T\end{array}\right.$$E2

$$\left\{\begin{array}{lll}{\partial }_{t}w\left(t,x\right)+\alpha {\partial }_{\mathit{xxx}}^{3}w\left(t,x\right)-\beta {\partial }_{\mathit{txx}}^{3}w\left(t,x\right)+\gamma {\partial }_{x}w\left(t,x\right)& =0& \text{in}{\mathrm{\Omega }}_{g_T}\cup {\mathrm{\Omega }}_{d_T}\\ w\left(0,x\right)& =0& \text{in}{\mathrm{\Omega }}_g\cup {\mathrm{\Omega }}_d\\ w& =u& \text{in}{\mathrm{\Sigma }}_{g_T}\cup {\mathrm{\Sigma }}_{d_T}\\ \underset{\mid x\mid \to +\infty }{lim}w\left(t,x\right)& =0& \mathrm{at}]0,T[\end{array}\right.$$E3

The main object of this section is to prove the following result.

Lemma 2.1 The solution of the evolution Eq. (3) satisfies the following integral equations

$$\begin{array}{c}w\left(t,a\right)-U_2{\mathcal{L}}^{-1}\left(\frac{{\lambda }_1{\left(s\right)}^2}{s}\right)\ast w_x\left(t,a\right)-U_2{\mathcal{L}}^{-1}\left(\frac{{\lambda }_1\left(s\right)}{s}\right)\ast w_{\mathit{xx}}\left(t,a\right)=0\\ w\left(t,b\right)-{\mathcal{L}}^{-1}\left(\frac{1}{{\lambda }_1{\left(s\right)}^2}\right)\ast w_{\mathit{xx}}\left(t,b\right)=0,\\ w_x\left(t,b\right)-{\mathcal{L}}^{-1}\left(\frac{1}{{\lambda }_1\left(s\right)}\right)\ast w_{\mathit{xx}}\left(t,b\right)=0\end{array}$$E4

where ${\mathcal{L}}^{-1}\left(f\left(s\right)\right)$ stands for the inverse Laplace transform of $f$, $\ast$ denotes the convolution operator and ${\lambda }_1$ a function of the time co-variable $s$.

Proof. We apply the Laplace transformation with respect to the time variable $t$ to the exterior problems (3), recall the Laplace transformation

$$\mathcal{L}\left(w\right)\left(s,x\right)≔\tilde{w}\left(s,x\right)={\int }_0^{+\infty }w\left(t,x\right)e^{-\mathit{ts}}\mathit{dt},$$E5

where $s$ stands for the co-variable of time $t$ and verify $\mathfrak{R}\left(s\right)>0$.

We obtain

$$s\tilde{w}\left(s,x\right)+\alpha {\partial }_{\mathit{xxx}}\tilde{w}-\mathit{\beta s}{\partial }_{\mathit{xx}}\tilde{w}+\gamma {\partial }_x\tilde{w}=0,x\ge b,x\le a,\mathfrak{R}\left(s\right)>0$$E6

which is a cubic ordinary differential equation whose solutions are of the form are given explicitly by

$$\widehat{w}\left(s,x\right)=c_1\left(s\right){\mathrm{e}}^{{\lambda }_1\left(s\right)x}+c_2\left(s\right){\mathrm{e}}^{{\lambda }_2\left(s\right)x}+c_3\left(s\right){\mathrm{e}}^{{\lambda }_3\left(s\right)x},x\in \mathbb{R}\backslash \left[a,b\right]$$E7

where ${\lambda }_1\left(s\right),{\lambda }_2\left(s\right),{\lambda }_3\left(s\right)$ denote the roots of the depressed cubic equation

$${\mathit{\alpha \lambda }}^3-\mathit{\beta s}{\lambda }^2+\mathit{\gamma \lambda }+s=0$$E8

The three solutions are given by

$${\lambda }_k\left(s\right)={\mathrm{j}}^{k-1}\zeta \left(s\right)-\frac{{\mathrm{\Theta }}_1\left(s\right)}{{\mathrm{j}}^{k-1}\zeta \left(s\right)}+{\mathrm{\Theta }}_2\left(s\right),k=\mathrm{1,2,3}$$E9

where the complex $\mathrm{j}$ is given by $\mathrm{j}=exp\left(2\mathit{i\pi }/3\right)$,

$$\begin{array}{ll}\zeta \left(s\right)& ={\left(\frac{\sqrt{4\alpha {\gamma }^3-{\beta }^2{s}^2{\gamma }^2+18\alpha \beta s^2\gamma -4{\beta }^3{s}^4+27{\alpha }^2{s}^2}}{23^{32}{\alpha }^2}-\frac{9\alpha \beta s\gamma -2{\beta }^3{s}^3+27{\alpha }^{2}s}{54{\alpha }^3}\right)}^{13},\\ {\mathrm{\Theta }}_1\left(s\right)& =\frac{3\alpha \gamma -{\beta }^2{s}^2}{9{\alpha }^2},\\ {\mathrm{\Theta }}_2\left(s\right)& =\frac{\mathit{\beta s}}{3\alpha }.\end{array}$$E10

Assume that $\mathfrak{R}\left(s\right)>\frac{2{\beta }^3+9{\mathit{\gamma \alpha }}^2}{27{\alpha }^2}$, then roots of the cubic Eq. (8) possess the following separation property

$$\mathfrak{R}\left({\lambda }_1\right)<\frac{\beta }{3\alpha },\mathfrak{R}\left({\lambda }_2\right)>\frac{\beta }{3\alpha }\text{and}\mathfrak{R}\left({\lambda }_3\right)>\frac{\beta }{3\alpha }.$$E11

In fact, we consider the change of variable $\lambda =z-\frac{\beta }{3\alpha }$. Then the cubic Eq. (8) becomes $z^3+\mathit{pz}+q=0$ such that $p=\frac{3\mathit{\alpha \gamma }-{\beta }^2}{3{\alpha }^2}\text{and}q=\frac{27{\alpha }^{2}s+2{\beta }^3-9\mathit{\alpha \beta \gamma }}{27\alpha }$.

Hence under the condition $\mathfrak{R}\left(q\right)=\frac{27{\alpha }^2\mathfrak{R}\left(s\right)+2{\beta }^3-9\mathit{\alpha \beta \gamma }}{27\alpha }$, it follows that the roots $r_i,i=\mathrm{1,2,3}$ of the equation $z^3+\mathit{pz}+q=0$ satisfy

$$\mathfrak{R}\left(r_1\right)<0,\mathfrak{R}\left(r_2\right)>0,\mathfrak{R}\left(r_3\right)>0$$

Now back to Eq. (7), for $x\ge b$ we have from the infinite condition that the coefficients $c_2$ and $c_3$ must vanishe, hence $\tilde{w}\left(x,s\right)=c_1\left(s\right)e^{r_1\left(s\right)x}$, deriving over $x$ and using the continuity of $w$ in the interface yield

$$\widehat{w}\left(s,b\right)-\frac{1}{{\lambda }_1^2\left(s\right)}{\widehat{w}}_{\mathit{xx}}\left(s,b\right)=0,{\widehat{w}}_x\left(s,b\right)-\frac{1}{{\lambda }_1\left(s\right)}{\widehat{w}}_{\mathit{xx}}\left(s,b\right)=0$$E12

Idem for $x\le a$, we have $c_1=0$ and hence

$${\widehat{w}}_{\mathit{xx}}\left(s,a\right)-\left({\lambda }_2\left(s\right)+{\lambda }_3\left(s\right)\right){\widehat{w}}_x\left(s,a\right)+{\lambda }_2\left(s\right){\lambda }_3\left(s\right)\widehat{w}\left(s,a\right)=0.$$E13

As ${\lambda }_1,{\lambda }_2,$ and ${\lambda }_3$ are roots of the cubic Eq. (8) we obtain immediately

$${\lambda }_1\left(s\right){\lambda }_2\left(s\right){\lambda }_3\left(s\right)=-\frac{s}{\alpha }\text{and}{\lambda }_2\left(s\right)+{\lambda }_3\left(s\right)+{\lambda }_1\left(s\right)=-\frac{\beta }{\alpha }$$E14

Then the Eq. (13) becomes in terms of ${\lambda }_1\left(s\right)$

$${\widehat{w}}_{\mathit{xx}}\left(s,a\right)+\left({\lambda }_1\left(s\right)+\frac{\beta }{s}\right){\widehat{w}}_x\left(s,a\right)-\frac{s}{{\mathit{\alpha \lambda }}_1\left(s\right)}\widehat{w}\left(s,a\right)=0$$E15

Now applying the inverse Laplace transform to Eqs. (8) and (10), we infer

$$\begin{array}{c}w\left(t,a\right)-\alpha {\mathcal{L}}^{-1}\left(\frac{{\lambda }_1^2\left(s\right)s+{\lambda }_1\left(s\right)\beta }{s}\right)\ast w_x\left(t,a\right)-\alpha {\mathcal{L}}^{-1}\left(\frac{{\lambda }_1\left(s\right)}{s}\right)\ast w_{\mathit{xx}}\left(t,a\right)=0\\ w\left(t,b\right)-{\mathcal{L}}^{-1}\left(\frac{1}{{\lambda }_1^2\left(s\right)}\right)\ast w_{\mathit{xx}}\left(t,b\right)=0,w_x\left(t,b\right)-{\mathcal{L}}^{-1}\left(\frac{1}{{\lambda }_1\left(s\right)}\right)\ast w_{\mathit{xx}}\left(t,b\right)=0\end{array}$$E16

Therefore, we get the following result describing the problem in the bounded domain satisfied by the restriction on ${\mathrm{\Omega }}_T$ of the original problem (1).

Theorem 1.1 Let $\alpha ,\beta$ be non negative numbers and $\gamma \in \mathbb{R}$. The restriction of (1) to $\mathrm{\Omega }$ is described by the following Initial Boundary Value Problem (IBVP)

$$\left\{\begin{array}{ll}{\partial }_{t}u+\alpha {\partial }_{\mathit{xxx}}u-\beta {\partial }_{\mathit{txx}}u+\gamma {\partial }_{x}u=0& \text{in}{\mathrm{\Omega }}_T\\ u\left(0,x\right)=u_0\left(x\right)& \mathrm{at}\mathrm{\Omega }\\ {\partial }_{n}u\left(t,x\right)=\mathit{Bu}\left(t,x\right)& \mathrm{at}{\mathrm{\Sigma }}_T\end{array}\right.$$E17

where $B$ is derived on ${\mathrm{\Sigma }}_T$ from equations

$$u\left(t,a\right)-\alpha {\mathcal{L}}^{-1}\left(\frac{{\lambda }_1^2\left(s\right)s+{\mathit{\beta \lambda }}_1\left(s\right)}{s}\right)\ast u_x\left(t,a\right)-\alpha {\mathcal{L}}^{-1}\left(\frac{{\lambda }_1\left(s\right)}{s}\right)\ast u_{\mathit{xx}}\left(t,a\right)=0$$

$$u\left(t,b\right)-{\mathcal{L}}^{-1}\left(\frac{1}{{\lambda }_1^2\left(s\right)}\right)\ast u_{\mathit{xx}}\left(t,b\right)=0,u_x\left(t,b\right)-{\mathcal{L}}^{-1}\left(\frac{1}{{\lambda }_1\left(s\right)}\right)\ast u_{\mathit{xx}}\left(t,b\right)=0$$

We emphasize that those boundary conditions strongly depend on $\alpha$ and $\beta$ through the root ${\lambda }_1\left(s\right)$. Some simplifications can be obtained for particular cases allowing direct evaluation of the inverse Laplace transform. Taking for example the BBM equation (for $\alpha =0$), we can get after applying Laplace transformation to (3),

$$\begin{array}{ll}{\partial }_x\tilde{w}\left(s,b\right)& =\frac{\gamma }{2\beta s}\tilde{w}\left(s,b\right)-\frac{\sqrt{{\gamma }^2+4\beta s^2}}{2\mathit{\beta s}}\tilde{w}\left(s,b\right)\\ & =\frac{\gamma }{2\beta }\frac{\tilde{w}\left(s,b\right)}{s}-\frac{{\gamma }^2}{2\beta }\frac{1}{\sqrt{{\gamma }^2+4\beta s^2}}\frac{\tilde{w}\left(s,b\right)}{s}-2\frac{s}{\sqrt{{\gamma }^2+4\beta s^2}}\tilde{w}\left(s,b\right).\end{array}$$E18

$$\begin{array}{ll}{\partial }_x\tilde{w}\left(s,a\right)& =\frac{\gamma }{2\mathit{\beta s}}\tilde{w}\left(s,a\right)+\frac{\sqrt{{\gamma }^2+4\beta s^2}}{2\mathit{\beta s}}\tilde{w}\left(s,a\right)\\ & =\frac{\gamma }{2\beta }\frac{\tilde{w}\left(s,a\right)}{s}+\frac{{\gamma }^2}{2\beta }\frac{1}{\sqrt{{\gamma }^2+4\beta s^2}}\frac{\tilde{w}\left(s,a\right)}{s}+2\frac{s}{\sqrt{{\gamma }^2+4\beta s^2}}\tilde{w}\left(s,a\right).\end{array}$$E19

In this case, we obtain convolution products with Bessel functions after the Laplace inverse transformation as follows

$${\partial }_{x}w\left(t,b\right)=C_1{I}_t\left(w\right)\left(t\right)-C_2\left(J_0^c\ast I_t\left(w\right)\right)\left(t\right)+C_3\left(J_1^c\ast w\right)\left(t\right),$$

$${\partial }_{x}w\left(t,a\right)=C_1{I}_t\left(w\right)\left(t\right)+C_2\left(J_0^c\ast I_t\left(w\right)\right)\left(t\right)-C_3\left(J_1^c\ast w\right)\left(t\right).$$

where we have used the expressions

$$\sqrt{\frac{c}{c^2+s^2}}=\mathcal{L}\left(H\left(t\right)J_0\left(\mathit{ct}\right)\right)\left(s\right)≔\mathcal{L}\left(H\left(t\right)J_0^c\left(t\right)\right)\left(s\right),\text{and}{\left(J_0^c\right)}^{\prime }=-J_1^c,$$

and the notations $c=\frac{\gamma }{2\sqrt{\beta }},C_1=\frac{\gamma }{2\beta },C_2={\left(\frac{\gamma }{2}\right)}^{\frac{3}{2}}{\beta }^{-\frac{5}{4}},C_3=\sqrt{\frac{2}{\gamma \sqrt{\beta }}}.$.

Recall that the Bessel functions can be defined by the following integrals

$$J_0^c\left(t\right)=\frac{1}{\pi }{\int }_0^{\pi }\mathit{cos}\left(\mathit{\text{ctsinτ}}\right)\mathit{d\tau },J_1^c\left(t\right)=\frac{1}{\pi }{\int }_0^{\pi }\mathit{cos}\left(\mathit{\text{ctsinτ}}-\tau \right)\mathit{d\tau }.$$

From this, we may compute

$${\partial }_{n}w\left(t,b\right)=-C_1{I}_t\left(w\right)\left(t\right)+C_2\left(J_0^c\ast I_t\left(w\right)\right)\left(t\right)-C_3\left(J_1^c\ast w\right)\left(t\right),$$

$${\partial }_{n}w\left(t,a\right)=C_1{I}_t\left(w\right)\left(t\right)+C_2\left(J_0^c\ast I_t\left(w\right)\right)\left(t\right)-C_3\left(J_1^c\ast w\right)\left(t\right).$$

Thus the boundary operator $B$ in (2) writes, in the case $\alpha =0$,

$$\mathit{Bu}\left(t,x\right)≔\left\{\begin{array}{rr}{C}_1{I}_{t}u\left(t\right)+C_2\left(J_0^c\ast I_{t}u\right)\left(t\right)-C_3\left(J_1^c\ast u\right)\left(t\right)& x=a\\ -C_1{I}_{t}u\left(t\right)+C_2\left(J_0^c\ast I_{t}u\right)\left(t\right)-C_3\left(J_1^c\ast u\right)\left(t\right)& x=b\end{array}\right.$$E20

Next, we propose an approximation, always for the case $\alpha =0$, of the BBM equation in ${\mathrm{\Omega }}_T$ supplemented with constructed boundary conditions.

2.1 Numerical approximation

This subsection is devoted to the numerical approximation of the obtained IBVP (17) for $\alpha =0$ and $B$ given in (20). Our strategy is to seek numerical simulations that permits to avoid any boundary reflections and in some way renders the fully discrete scheme unconditionally stable.

Let $N,M$ be integers, we define time step $\mathrm{\Delta }t=\frac{T}{M}$ and spatial step $h=\frac{b-a}{N}$. The grids $t_n=n\mathrm{\Delta }t,0\le n\le M$ and $x_i=a+\mathit{ih},0\le i\le N$ are used to discretize ${\mathrm{\Omega }}_T$. Throughout this paper, we denote $u_i^n$ the considered approximation of $u\left(t_n,x_i\right)$ and the set ${\mathbb{N}}_k^l=\left\{m\in \mathbb{N},k\le m\le l\right\}.$

2.1.1 Approximation of the governing equation

We describe a discretization for the BBM equation by the Crank-Nicholson time scheme as follows

$$\frac{u\left(t_{n+1},x\right)-u\left(t_n,x\right)}{\mathrm{\Delta }t}-\beta \frac{{\partial }_{\mathit{xx}}u\left(t_{n+1},x\right)-{\partial }_{\mathit{xx}}u\left(t_n,x\right)}{\mathrm{\Delta }t}+\gamma \frac{{\partial }_{x}u\left(t_{n+1},x\right)+{\partial }_{x}u\left(t_n,x\right)}{2}=0,$$E21

$$u\left(t_0,x\right)=u_0\left(x\right),n\in {\mathbb{N}}_0^{M-1}.$$

For the space finite difference scheme, we use the approximations

$${\partial }_{x}u\left(t,x\right)\approx \frac{1}{2h}\left(u\left(t,x+h\right)-u(t,x-h)\right),$$

$${\partial }_{\mathit{xx}}u\left(t,x\right)\approx \frac{1}{h^2}\left(u\left(t,x+h\right)-2u(t,x)+u(t,x-h)\right).$$

The fully discretization then writes,

$$\begin{array}{l}\frac{u_i^{n+1}-u_i^n}{\mathrm{\Delta }t}-\beta \frac{\left(u_{i+1}^{n+1}-2u_i^{n+1}+u_{i-1}^{n+1}\right)-\left(u_{i+1}^n-2u_i^n+u_{i-1}^n\right)}{h^2\mathrm{\Delta }t}\\ +\gamma \frac{\left(u_{i+1}^{n+1}-u_{i-1}^{n+1}\right)+\left(u_{i+1}^n-u_{i-1}^n\right)}{4h}=0,\left(i,n\right)\in {\mathbb{N}}_1^{N-1}\times {\mathbb{N}}_0^{M-1}\\ u_i^0=u_0\left(x_i\right),i\in {\mathbb{N}}_0^N.\end{array}$$E22

2.1.2 Approximation of the boundary condition

The constructed boundary conditions (BC) contains time convolutions that are non-local and introduces many difficulties, for example, using a direct implementation leads to long and low accuracy. Several techniques have been used to overcome these problems by trying to localize the BC, see [5, 6, 7, 8] for more details. The resulting localized BC are easy to implement and more efficient but tends to depend sensitively on the initial data. In our case, we utilize the Discrete Galerkin Method. The BC are formulated as Fredholm integral equations of second kind. The basic idea is to write the boundary condition on (20) in the form

$${\partial }_{n}u\left(t,a\right)-{\int }_0^t{\overline{K}}_1\left(t,s\right)u\left(s,a\right)\mathit{ds}=0,$$E23

$${\partial }_{n}u\left(t,b\right)-{\int }_0^t{\overline{K}}_2\left(t,s\right)u\left(s,b\right)\mathit{ds}=0.$$E24

where, the introduced Kernels ${\overline{K}}_1,{\overline{K}}_2$ represent a linear combination of the two Bessel functions of order 0 and 1 at the time $\left(t-s\right)$. After a space discretization we obtain

$$u\left(t,x_0\right)-{\int }_0^t{K}_1\left(t,s\right)u\left(s,x_0\right)\mathit{ds}=u\left(t,x_1\right),$$E25

$$u\left(t,x_N\right)-{\int }_0^t{K}_2\left(t,s\right)u\left(s,x_N\right)\mathit{ds}=u\left(t,x_{N-1}\right),$$E26

where $K_1=-h{\overline{K}}_1$ and $K_2=h{\overline{K}}_2$. Both resulting Eqs. (25) and (26) can be identified to the linear integral equation

$$y\left(t\right)-{\int }_{D}K\left(t,s\right)y\left(s\right)\mathit{d\sigma }\left(s\right)=z\left(t\right),t\in D.$$E27

The Eq. (27) is a Fredholm integral equation of second kind, where $D$ is a closed bounded set in ${\mathbb{R}}^m$, with $m\ge 1$. The approximation of such integral equation could be made by a discrete Galerkin method using the quadrature rule of Gauss-Legendre as presented in [9]. Based on this, the BC can be similarly discretized while considering the domain $D$ as the time interval $\left[0,t\right]$. Precisely, we use the Gauss Legendre Quadrature of order $q$, labeled ${\mathit{GLQ}}_q$ with zeros ${\xi }_j$ and weights $w_j$ being in the interval $\left[-1,1\right]$ for $j\in {\mathbb{N}}_0^q$. Let $i\in {\mathbb{N}}_0^{N-1}$, we introduce the following transformation

$$\begin{array}{llll}{F}_i:& \left[-1,1\right]& \to & \left[t_i,t_{i+1}\right]\\ & \xi & \mapsto & t_i\frac{1-\xi }{2}+t_{i+1}\frac{1+\xi }{2}\end{array}$$E28

The approximation of the BC is now given by

$$\begin{array}{c}u\left(t_{n+1},x_0\right)-{\int }_0^{t_{n+1}}{K}_1\left(t_{n+1},s\right)u\left(s,x_0\right)\mathit{ds}=u\left(t_{n+1},x_1\right)\end{array}$$E29

$$\begin{array}{c}u\left(t_{n+1},x_N\right)-{\int }_0^{t_{n+1}}{K}_2\left(t_{n+1},s\right)u\left(s,x_N\right)\mathit{ds}=u\left(t_{n+1},x_{N-1}\right)\end{array}$$E30

that is

$$u_0^{n+1}-{\int }_0^{t_{n+1}}{K}_1\left(t_{n+1},s\right)u\left(s,x_0\right)\mathit{ds}=u_1^{n+1}$$E31

$$u_N^{n+1}-{\int }_0^{t_{n+1}}{K}_2\left(t_{n+1},s\right)u\left(s,x_N\right)\mathit{ds}=u_{N-1}^{n+1}$$E32

For the seek of simplicity, we rewrite the integral terms of (31) and (32) in the form

$$\begin{array}{l}{\int }_0^{t_{n+1}}{K}_i\left(t_{n+1},s\right)u\left(s,x_k\right)\mathit{ds}=A_i{\int }_0^{t_{n+1}}u\left(s,x_k\right)\mathit{ds}+B_i{\int }_0^{t_{n+1}}{J}_0^c\left(t_{n+1}-s\right)\left({\int }_0^{s}u\left(r,x_k\right)\mathit{dr}\right)\mathit{ds}\\ +D_i{\int }_0^{t_{n+1}}{J}_0^c\left(t_{n+1}-s\right)u\left(s,x_k\right)\mathit{ds}\\ ≔A_i{I}_1\left(t_{n+1},x_k\right)+B_i{I}_2\left(t_{n+1},x_k\right)+D_i{I}_3\left(t_{n+1},x_k\right),\end{array}$$

such that $\left(i,k\right)\in \left\{\left(1,0\right),\left(2,N\right)\right\}$, $A_1=A_2=-{\mathit{hC}}_1,B_2=-B_1={\mathit{hC}}_2$ and $D_1=-D_2={\mathit{hC}}_3$, all constants $C_i$ are defined in (??). Thus, basing on the approach presented in [9], the ${\mathit{GLQ}}_q$ applied to the integrals previously defined is described by the following, for $n\in {\mathbb{N}}_0^{M-1}$ and $k\in \left\{0,N\right\}$,

$$I_1\left(t_{n+1},x_k\right)=\frac{\mathrm{\Delta }t}{2}\sum _{i=0}^n\sum _{j=0}^q{w}_{j}u\left(F_i\left({\xi }_j\right),x_k\right).$$E33

The second integral is more complicated since it involves two composing integrals, using Gauss-Legendre quadrature twice yields

$$I_2\left(t_{n+1},x_k\right)=\frac{\mathrm{\Delta }{t}^2}{4}\sum _{i=0}^n\sum _{j=0}^q{w}_j{J}_0\left(t_{n+1}-F_i\left({\xi }_j\right)\right)\sum _{l=0}^j\sum _{m=0}^q{w}_{m}u\left(F_l\left({\xi }_m\right),x_k\right)$$E34

and the remained integral is approximated by

$$I_3\left(t_{n+1},x_k\right)=\frac{\mathrm{\Delta }t}{2}\sum _{i=0}^n\sum _{j=0}^q{w}_j{J}_1\left(t_{n+1}-F_i\left({\xi }_j\right)\right)u\left(F_i\left({\xi }_j\right),x_k\right).$$E35

From approximations (33)–(35), the numerical solution on the interface of (17) can be given by

$$u_0^{n+1}=u_1^{n+1}+f_1^{n+1},$$E36

$$u_N^{n+1}=u_{N-1}^{n+1}+f_2^{n+1}.$$E37

We accomplish this by simply adding (36) and (37) to the discretization of the interior governing Eq. (22). We obtain an implicit scheme that we illustrate by the following system in matrix form

$$\left\{\begin{array}{ll}\left(I+\tilde{A}+\tilde{B}\right)U^{n+1}& =\left(I+A-B\right)U^n+F^n,n\in {\mathbb{N}}_0^{M-1},\\ u_0^{n+1}& =u_1^{n+1}+f_1^{n+1},\\ u_N^{n+1}& =u_{N-1}^{n+1}+f_2^{n+1}.\end{array}\right.$$E38

with

$$U^n=\left(\begin{array}{c}{u}_1^n\\ ⋮\\ u_{N-1}^n\end{array}\right)I=\left(\begin{array}{ccc}1& & 0\\ & \ddots & \\ 0& & 1\end{array}\right),F^n=\left(\begin{array}{c}\left(C_{d_1}+C_{d_2}\right)\left(f_1^n-u_1^n\right)\\ 0\\ ⋮\\ 0\\ \left(C_{d_1}-C_{d_2}\right)\left(f_2^n-u_{N-1}^n\right)\end{array}\right),$$

$$\begin{array}{l}A=C_{d_1}\left(\begin{array}{ccc}2& -1& 0\\ -1& \ddots & -1\\ 0& -1& 2\end{array}\right),B=C_{d_2}\left(\begin{array}{ccc}0& 1& 0\\ -1& \ddots & 1\\ 0& -1& 0\end{array}\right)\\ \end{array}$$

and

$$\begin{array}{l}\tilde{A}=C_{d_1}\left(\begin{array}{ccc}2-C_{d_1}& -1& 0\\ -1& \ddots & -1\\ 0& -1& 2-C_{d_1}\end{array}\right),\tilde{B}=C_{d_2}\left(\begin{array}{ccc}-C_{d_2}& 1& 0\\ -1& \ddots & 1\\ 0& -1& C_{d_2}\end{array}\right).\\ \end{array}$$

where the discretization constants are $C_{d_1}=\frac{\beta }{h^2},C_{d_2}=\frac{\gamma \mathrm{\Delta }t}{4h}$.

3. Non standard boundary conditions for the BBM equation

In [4], we have presented a new method to derive transparent boundary conditions for the BBM equation. These boundary conditions have the advantage of being local in time but needs an additional function construct the BC which means bigger system to be solved. We recall that the problem designed to be the restriction in ${\mathrm{\Omega }}_T$ of the BBM initial Eq. (1) with $\alpha =0$ is given by

$$\left\{\begin{array}{ll}{\partial }_{t}u-\beta {\partial }_{\mathit{txx}}u+\gamma {\partial }_{x}u=f& \text{in}]0,T]\times ]a,b[\\ {\partial }_{t}v-\beta {\partial }_{\mathit{txx}}v+\gamma {\partial }_{x}v=g& \text{in}]0,T]\times ]a,b[\\ {\partial }_{x}u=v& \mathrm{on}]0,T]\times \left\{a,b\right\}\\ \beta {\partial }_{\mathit{tx}}v-\mathit{\gamma v}={\partial }_{t}u-f& \mathrm{on}]0,T]\times \left\{a,b\right\}\\ u\left(0,x\right)=u_0\left(x\right)& \mathrm{on}\left[a,b\right]\\ v\left(0,x\right)=v_0\left(x\right)& \mathrm{on}\left[a,b\right]\end{array}\right.$$E39

4. Numerical examples

We take an initial condition as solitary wave like function locally supported in $\mathrm{\Omega }$. The evolution of the solutions are plotted in different time steps before, under and after traveling the right boundary of the considered bounded domain. We save a reference solution that is numerically calculated in a broaden domain of $\mathrm{\Omega }$ with Dirichlet boundary condition. We compute infinite error between numerical solutions using both formulations presented in this paper and the reference solution. We denote ${\mathit{GLQ}}_i$ for approached solution with BEM and Gauss Legendre Quadrature in (2) for $i\in \left\{\mathrm{0,1,2}\right\}$, while NSBC refers to numerical solution with non standard boundary conditions given in (3). We define the following errors

$$\parallel \left(u-u_{\mathit{ref}}\right)\left(t,a\right){\parallel }_{\infty }=\underset{t\in \left[0,T\right]}{sup}\mid u\left(t,a\right)-u_{\mathit{ref}}\left(t,a\right)\mid ,$$E40

$$\parallel \left(u-u_{\mathit{ref}}\right)\left(t,b\right){\parallel }_{\infty }=\underset{t\in \left[0,T\right]}{sup}\mid u\left(t,b\right)-u_{\mathit{ref}}\left(t,b\right)\mid ,$$E41

$$\parallel \left(u-u_{\mathit{ref}}\right)\left(t,x\right){\parallel }_{\infty }=\underset{\left(t,x\right)\in \left[0,T\right]\times ]a,b[}{sup}\mid u\left(t,x\right)-u_{\mathit{ref}}\left(t,x\right)\mid .$$E42

Let $\beta =\gamma =1$, $T=10$ and $u_0\left(x\right)={\mathit{sech}}^2\left(x\right)$. The considered initial data is locally supported in the interval $\mathrm{\Omega }=]-10,10[$, since $u_0\left(10\right)=u_0\left(-10\right)\approx {\mathrm{8,25.10}}^{-9}$. We fix $h={10}^{-2}$ and we vary the time step $\mathrm{\Delta }t$. For a better comparison of these methods, we compute CPU time, in seconds, needed for each one to obtain numerical solution.

Table 1 shows that both methods give a good approximation of the restriction to ${\mathrm{\Omega }}_T$ of the reference solution. $\mathit{\text{NSBC}}$ give better approximation than $\mathit{GLQ}$. We can remark a slow convergence of ${\mathit{GLQ}}_i$ with respect to $i$ and also time step. However, NSBC gives a good approximation in as much as $\mathrm{\Delta }t$ goes to zero. Furthermore, ${\mathit{GLQ}}_i$ is more expensive in CPU time when $i$ increases than NSBC due to the presence of non local convolutions in time in the boundary condition.

BC	$\mathit{dt}$	${‖\left(u-u_{\mathit{ref}}\right)(t,a)‖}_{\infty }$	${‖\left(u-u_{\mathit{ref}}\right)(t,b)‖}_{\infty }$	${‖\left(u-u_{\mathit{ref}}\right)(t,x)‖}_{\infty }$	CPU time(s)
${\mathit{GLQ}}_0$	${10}^{-2}$	$1.06\times {10}^{-4}$	$6.33\times {10}^{-2}$	$6.03\times {10}^{-2}$	4
${\mathit{GLQ}}_1$		$1.06\times {10}^{-4}$	$6.3\times {10}^{-2}$	$5.9304\times {10}^{-2}$	9
${\mathit{GLQ}}_2$		$1.06\times {10}^{-4}$	$5.895\times {10}^{-2}$	$5.9187\times {10}^{-2}$	15
$\mathit{\text{NSBC}}$		$7.39\times {10}^{-5}$	$5.8\times {10}^{-3}$	$6.8\times {10}^{-3}$	5
${\mathit{GLQ}}_0$	${10}^{-3}$	$1.06\times {10}^{-4}$	$5.92\times {10}^{-2}$	$5.86\times {10}^{-2}$	39
${\mathit{GLQ}}_1$		$1.06\times {10}^{-4}$	$5.86\times {10}^{-2}$	$5.857\times {10}^{-2}$	62
${\mathit{GLQ}}_2$		$1.06\times {10}^{-4}$	$5.852\times {10}^{-2}$	$5.85\times {10}^{-2}$	80
$\mathit{\text{NSBC}}$		$6.7\times {10}^{-5}$	$1.41\times {10}^{-3}$	$1.4\times {10}^{-3}$	50
${\mathit{GLQ}}_0$	${10}^{-4}$	$1.06\times {10}^{-4}$	$5.855\times {10}^{-2}$	$5.858\times {10}^{-2}$	436
${\mathit{GLQ}}_1$		$1.06\times {10}^{-4}$	$5.853\times {10}^{-2}$	$5.852\times {10}^{-2}$	1040
${\mathit{GLQ}}_2$		$1.06\times {10}^{-4}$	$5.85\times {10}^{-2}$	$5.85\times {10}^{-2}$	3205
$\mathit{\text{NSBC}}$		$5.9\times {10}^{-6}$	$7.35\times {10}^{-4}$	$7.3\times {10}^{-4}$	1015

Table 1.

Infinite errors using different boundary conditions.

We also plot in Figure 1, captions at different times of either reference solution and approximated solutions using NSBC and ${\mathit{GLQ}}_i$ for $i=2$. We can see that NSBC follows the refrence solution better than GLQ especially at last times in the right figure. One remarks that no reflections turn back to the bounded domain when the wave is going out from the right boundary using both methods.

Figure 1.

Reference solution and approximated solutions NSBC and ${\mathit{GLQ}}_i$ for $i=2$ at different times for $\mathrm{\Delta }t={10}^{-3}$.

5. Conclusion

We have compared two methods of deriving and approaching boundary conditions for the BBM equation. We presented the BEM for a general equation that is the mixed BBM-KdV equation and that shows the hardness to put easy implemented BC. Furthermore, being non local in time, BC seems to be low accurate and slowly convergent as presented in numerical example. However, this point opens many possibilities trying to improve the accuracy of such BC whether by improving the approximation of convolution product, that comes from Inverse Laplace transformation, via quadrature or exploring a numerical equivalent to such operation such as $\mathcal{Z}$ transformation. We have proposed an other manner to derive local BC that gives better approximation than non local BC. All these conclusions have been made in one space dimension but nothing can be said about the comparison in higher dimension to decide which method is more adapted, this matter will be our interest in future works.