Abstract
We are concerned with stable spectral collocation solutions to nonperiodic Benjamin Bona Mahony (BBM), modified BBM and Benjamin Bona MahonyBurgers (BBMB) initial value problems on the real axis. The spectral collocation is based alternatively on the scaled Hermite and sinc functions. In order to march in time we use several one step and linear multistep finite difference schemes such that the method of lines (MoL) involved is stable in sense of Lax. The method based on Hermite functions ensures the correct behavior of the solutions at large spatial distances and in long time periods. In order to prove the stability we use the pseudospectra of the linearized spatial discretization operators. The extent at which the energy integral of BBM model is conserved over time is analyzed for Hermite collocation along with various finite difference schemes. This analysis has been fairly useful in optimizing the scaling parameter. The effectiveness of our approach has been confirmed by some numerical experiments.
Authors
CălinIoan Gheorghiu
(Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy)
Keywords
References
See the expanding block below.
Paper coordinates
C.I. Gheorghiu, Stable spectral collocation solutions to a class of Benjamin Bona Mahony initial value problems. Appl. Math. Comp., 273 (2016) 10901099
doi: 10.1016/j.amc.2015.10.078
not available yet.
About this paper
Print ISSN
00963003
Online ISSN
Google Scholar Profile
google scholar link
Paper (preprint) in HTML form
Stable Spectral Collocation Solutions to a Class of Benjamin Bona Mahony Initial Value Problems
Abstract
We are concerned with stable spectral collocation solutions to nonperiodic Benjamin Bona Mahony (BBM), modified BBM and Benjamin Bona MahonyBurgers (BBMB) initial value problems on the real axis. The spectral collocation is based alternatively on the scaled Hermite and sinc functions. In order to march in time we use several one step and linear multistep finite difference schemes such that the method of lines (MoL) involved is stable in sense of Lax. The method based on Hermite functions ensures the correct behavior of the solutions at large spatial distances and in long time periods. In order to prove the stability we use the pseudospectra of the linearized spatial discretization operators. The extent at which the energy integral of BBM model is conserved over time is analyzed for Hermite collocation along with various finite difference schemes. This analysis has been fairly useful in optimizing the scaling parameter. The effectiveness of our approach has been confirmed by some numerical experiments.
Key words: regularized long wave equation; scaled Hermite collocation; method of lines; Lax stability; energy conservation;
1 Introduction
The purpose of this paper is the analysis and practical development of the method of lines (MoL for short) based on spectral collocation and various finite difference schemes as an efficient tool to solve the Cauchy’s problem attached to a class of Benjamin Bona Mahony (BBM for short) equations. The most representative equation reads
$${u}_{t}+{u}_{x}\alpha {u}_{xx}+{u}^{n}{u}_{x}{u}_{xxt}=0,0\le \alpha \le 1,n\in \mathbb{N},$$  (1) 
whose solution $u(x,t)$ is considered in a class of nonperiodic functions defined on $$ and $t\ge 0.$
For $\alpha =0$ and $n=1$ the equation (1) is called Benjamin Bona Mahony or the regularized longwave equation. It has been introduced by the above mentioned mathematicians in the early ’70 (see [1]) as a counterpart of the well known Kortewegde Vries (KdV) equation
$${u}_{t}+{u}_{x}+u{u}_{x}+{u}_{xxx}=0.$$  (2) 
When $\alpha =0$ and $n>1$ the equation (1) will be called modified BBM equation and for $\alpha >0$ and $n=1$ we will refer at (1) as to the Benjamin Bona MahonyBurgers (BBMB for short) one.
In (1), as well as in the formulation of (2), $u=u(x,t)$ represents the vertical displacement of the surface of the liquid (ideal fluid) from its equilibrium position, $t$ is the time and $x$ is the horizontal coordinate which increases in the direction of waves propagation. Both equations are in dimensionless form with the length scale taken to be the unperturbed depth $h$ of the liquid and time scale be ${\left(h/g\right)}^{1/2};$ $g$ is the gravity constant.
The authors of [1] compare in various ways their regularized longwave equation with KdV and argue that the latter has some shortcomings which makes that an unsuitable posed model for long waves (see also [2]). We do not want to get involved into the details of this debate but we observe that the BBM model has received less attention from numerical point of view than the KdV one.
One of the purposes of this work is to try to fill this gap providing reasonable numerical wave solutions without an artificial (empirical) truncation of the domain or periodicity assumptions.
For a model equation for long waves, periodic solutions have been obtained in the early era of modern computational mathematics by Galerkin method in [19]. On a finite interval, an initial(Dirichlet) boundary value problem attached to BBM has been solved by Galerkin method in [13] and by a CrankNicolsontype finite difference method in [14].
A kind of Hermite spectral method has been used in [11] in order to solve KdV, a modified KdV and the nonlinear Schroedinger equations. The so called KdV solitons of Boyd have been detected this way. However, in [3] it is observed that the Hermite functions are the exact eigenfunctions of quantum mechanics harmonic oscillator and are the asymptotic eigenfunctions for Mathieu’s equation, the prolate spheroidal wave equation, the associated Legendre equation, and Laplace’s Tidal equation. This close connection with the physics makes Hermite functions a fairly natural choice of basis functions for various fields of science and engineering.
In a recent laborious study [4] the authors extend the framework of the finite volume methods to dispersive water wave models, in particular to Boussinesq type systems.
In this context, we focus in this paper on the mathematical aspects of the spectral collocation methods, mainly based on Hermite functions, as an efficient tool in modeling the dispersion of long waves. Along with some implicit time schemes, the involved MoL, accurately conserves a quadratic invariant of the flow and resolves correctly the asymptotic behavior of the solutions at large distances.
Specifically, we will be concerned in the first instance with the Cauchy’s problem
$$  (3) 
where the initial data $g\left(x\right)$ satisfies the conditions

1.
$g\in {W}_{2}^{1}\left(\mathbb{R}\right)\cap {C}^{2}\left(\mathbb{R}\right),$

2.
$$
Then the problem (3) has a solution $u\in {\mathcal{C}}_{\mathrm{\infty}}^{2,\mathrm{\infty}}.$ At any $t\in [0,\mathrm{\infty})$ the solution $u$ and its derivative ${u}_{x}$ along with all their derivatives with respect to $t$ are asymptotically null i.e.,
$$u,{u}_{x},{u}_{t},{u}_{xt},\mathrm{\dots}\to 0asx\to \pm \mathrm{\infty}.$$  (4) 
Additionally, ${E}_{0}$ is conserved over time, i.e., $E\left(u\right)={E}_{0}$ where
$$E\left(u\right):={\int}_{\mathrm{\infty}}^{\mathrm{\infty}}\left({u}^{2}+{u}_{x}^{2}\right)\mathit{d}x.$$  (5) 
With respect to the first condition on the initial data $g,$ the authors of [1] observe that in fact an element $g$ of ${L}_{2}\left(\mathbb{R}\right)\cap {C}^{1}\left(\mathbb{R}\right)$ and its first derivative ${g}_{x}$, both converge to $0$ at $\pm \mathrm{\infty}.$ The functional $E$ introduced with (5) is called energy integral although it is not necessarily identifiable with the energy in the original physical problems.
In order to describe the smoothness of solution the authors quoted above introduce the space ${\mathcal{C}}_{T}$ of continuous functions $u(x,t)$ on $\mathbb{R}\times [0,T)$ equipped with the maximum norm and its subspace ${\mathcal{C}}_{T}^{l,m}$ of functions $u$ such that ${\partial}_{x}^{i}{\partial}_{t}^{l}u\in {\mathcal{C}}_{T}$ for $0\le i\le l$ and $0\le j\le m$ equipped with the subsequent norm. The space ${W}_{2}^{1}\left(\mathbb{R}\right)$ is the usual Sobolev space of ${L}_{2}$ functions over $\mathbb{R}$ with generalized first derivative which is also ${L}_{2}$ function.
We will heavily rely on this analytical result. Moreover this fundamental result on the existence, uniqueness and asymptotic behavior of nonperiodic solutions to (3) has suggested a collocation approach of the spatial dependence using Hermite (HC for short) as well as sinc functions (SC for short). In this context fairly accurate differentiation matrices corresponding to these functions are needed. They are available in [5] and in the seminal paper [20].
Actually, the HC method is used in the same spirit in which the Laguerre collocation (LC) has been successfully used in our previous papers [6] and [7]. With LC, as well as with HC, we avoid the frequently used empirical treatment of boundary value problems defined on the infinite domains based on the arbitrary truncation of the domain.
In order to march in time we have used gradually one step explicit methods such as RungeKutta (RK) methods of orders $2$, $3$ and $4,$ and linear multistep methods such as AdamsBashforth (AB for short), AdamsMoulton (AM for short) and leap frog schemes. The regions of stability of all these finite differences formulas along with the spectral properties of Hermite and sinc differentiation matrices of moderate orders $N$, i.e., $64\le N\le 128$ have assured the numerical stability of MoL. The necessary and sufficient conditions of Laxstability formulated in [9] and [16] have been rigorously observed.
A particular attention is paid to the conservation of energy integral over time. The time marching schemes as well as the spatial collocation schemes are evaluated with respect to this aspect. More precisely, the scaling parameter involved in the HC has been adjusted in order to optimize the conservation of this invariant. Meantime this scaling factor has also been chosen in order to enlarge the spatial interval of integration. This way we have succeeded in capturing the evolution of the initial front on fairly long time intervals. Some differences between HC and SC methods have been noticed in this context.
The paper is organized as follows. In Section 2 we solve the BBM problem as well as a modified BBM problem by Hermite and sinc spectral collocation methods along with RungeKutta method of order 4. Then we show the stability of this MoL as well as of the MoL based on the other finite difference schemes mentioned above. The last subsection is devoted to the conservation of the energy invariant. As we do not have a general rule to establish a correct value of the scaling parameter, this conservation process is in fact the key in its optimal choosing. The Section 3 parallelizes, up to some point, the previous one corresponding to the BBMB problem. Actually, in both sections the normality of the differentiation matrices as well as of the linearized spatial discretization operators have been quantified using a scalar measure, i.e., the Henrici’s number and their pseudospectra. Then we make use of the Trefethen’s rule of thumb, which states: a MoL is stable if the eigenvalues of the (linearized) spatial discretization operator, scaled by $\mathrm{\Delta}t$, lie in the stability region of the timediscretization operator (see [18], p. 103). Fortunately, the matrices involved in the problems at hand are not far from normal and thus this rule makes accurate predictions. In some way the rule can be thought of as an analogue of the well known CFL condition. Section 4 reviews our main results.
2 The MoL for BBM initial value problem
In order to solve problem (3) we rewrite the differential equation in the following convenient form
$$\left(1{\partial}_{x}^{2}\right){u}_{t}={\partial}_{x}\left(u+\frac{1}{2}{u}^{2}\right).$$  (6) 
2.1 The HC method
The interpolant of a map $f:\mathbb{R}\to \mathbb{R}$ reads
$${p}_{N1}(x):=\sum _{j=1}^{N}\frac{{e}^{{x}^{2}/2}}{{e}^{{x}_{j}^{2}/2}}{\varphi}_{j}\left(x\right){f}_{j},{f}_{j}:=f\left({x}_{j}\right)$$  (7) 
where
$${\varphi}_{j}\left(x\right):=\frac{{H}_{N}\left(x\right)}{{H}_{N}^{\prime}\left({x}_{j}\right)\left(x{x}_{j}\right)}.$$ 
${H}_{N}\left(x\right)$ is the Hermite polynomial of degree $N$ defined by the three term recurrence formula
$${H}_{n+1}\left(x\right):=\left(x{\alpha}_{n}\right){H}_{n}\left(x\right){\beta}_{n+1}^{2}{H}_{n1}\left(x\right),n=0,1,2,\mathrm{\dots},$$ 
with
$${H}_{0}\left(x\right)=1,{H}_{1}\left(x\right)=0,{\alpha}_{n}=0,{\beta}_{n}=1/2n.$$ 
The nodes ${x}_{1},\mathrm{\dots},{x}_{N}$ are the roots of ${H}_{N}\left(x\right)$ indexed in the ascending order. The Hermite differentiation matrices will be denoted by ${D}_{H}^{\left(l\right)}$ for $l=1,2,\mathrm{\dots}.$ They are provided along with the above roots by the MATLAB code herdif from [20]. Actually, we work with the scaled roots ${x}_{1}/b,\mathrm{\dots},{x}_{N}/b$ and consequently the differentiation matrices correspond to these scaled nodes. The Hermite differentiation process is exact for functions of the form
$${e}^{{x}^{2}{b}^{2}/2}p\left(x\right),$$ 
where $p\left(x\right)$ is a polynomial of degree $N1$ or less, assuming exact arithmetic. The scaling factor $b$ will be exploited in our numerical experiments in some respects.
However, the rate of convergence theory for Hermite functions is now well established for a fixed ($N$ independent) scale factor $b$ (see [3], § 17.4 ”Hermite functions”). This rate depends not only on the location of singularities, as it is true for all spectral expansions, but also upon the rate at which a given function $f\left(x\right)$ decays as $\leftx\right\to \mathrm{\infty}.$
The MoL applied to problem (3) along with HC means to approximate the solution $u(x,t)$ of this problem (using the form (6)) with
$${u}_{N,b}(x,t):=\sum _{j=1}^{N}\frac{{e}^{{x}^{2}/2}}{{e}^{{x}_{j}^{2}/2}}{\varphi}_{j}\left(x\right){U}_{j}\left(t\right).$$  (8) 
Thus, the conditions (4) are satisfied and the problem is transformed into the Cauchy’s problem for the following system of ODE
$$\{\begin{array}{c}(I{D}_{H}^{\left(2\right)}){U}_{t}={D}_{H}^{\left(1\right)}(U+\frac{1}{2}U{.}^{2}),t\in [0,T],\\ U\left(0\right)=g\left(X\right),\end{array}$$  (9) 
where $I$ is the unit matrix of order $N$, $X:={\left({x}_{1}/b\mathrm{\dots}{x}_{N}/b\right)}^{T}$ and
$$U\left(t\right):={\left({U}_{1}\left(t\right)\mathrm{\dots}{U}_{N}\left(t\right)\right)}^{T}.$$ 
The ”mass” matrix $I{D}_{H}^{\left(2\right)}$ is nonsingular and independent of $t$ and $U$, so this system is fairly amenable to be solved by some ode builtin functions of MATLAB. The period character $(.)$ in the r. h. s. of the semidiscretized problem (9) signifies the array power of vector $U.$
In order to assert the Lax stability of this MoL we have to analyze the linearized spatial discretization operator attached to (9). This is in fact the matrix
$${L}_{H}:={D}_{H}^{\left(1\right)}{\left(I{D}_{H}^{\left(2\right)}\right)}^{1}.$$  (10) 
The Frobenius norm conditioning and the Henrici’s number of the matrices ${D}_{H}^{\left(1\right)},{D}_{H}^{\left(2\right)}$ and ${L}_{H}$ are displayed in Table 1.
Remember that for a square matrix $A$ the Henrici’s number is defined by
$$Hen\left(A\right):=\left({\Vert A{A}^{\ast}{A}^{\ast}A\Vert}^{1/2}\right)/\Vert A\Vert ,$$ 
where $\parallel \cdot \parallel $ stands for the Frobenius norm of $A$ and ${A}^{\ast}$ is the conjugate transpose of $A.$ Thus a matrix is normal iff $A{A}^{\ast}{A}^{\ast}A=0$ and $0\le Hen\left(A\right)\le \sqrt[4]{2}$ ( see for instance our contribution [8]).
The sinc differentiation matrices are symmetric or skewsymmetric for even respectively odd values of the cut off parameter $N.$ The Hermite differentiation matrices look worse from this point of view but are still closed to the normality.
Matrices  ${D}^{\left(1\right)}$  ${D}^{\left(2\right)}$  $L$  

Method  Hermite  sinc  Hermite  sinc  Hermite  sinc 
Conditioning  9.44e+01  6.29e+01  6.35e+03  4.09e+03  4.67e+01  7.78e+02 
Henrici  1.39e01  0  2.94e02  0  2.27e02  8.69e01 
The $\epsilon $pseudospectra of ${L}_{H}$ are depicted in Fig. 1 for two values of scaling parameter. The eigenvalues of ${L}_{H}$ lie on the right hand side boundary of the stability region of RK3. For $b=1/3$ the operator ${L}_{H}$ is almost normal. This statement is based on its pseudospectrum from Fig. 1 a) and its Henrici’s number (the bold value in Tab. 1). Even for a larger value of scaling factor, i.e. $b=5$ half of the largest contour of the pseudospectrum, i.e., corresponding to $\epsilon =0.1,$ is situated inside the stability region and the other half lies within a distance of order $O\left(\epsilon \right)+O\left(\mathrm{\Delta}t\right)$ for ”small” $\epsilon $ and $\mathrm{\Delta}t$ approaching $0.$
This is the first reason in working with $$
In accordance with the above quoted Trefethen’s rule of thumb, this assures the numerical stability of MoL based on RK3 and HC.
A similar and even better situation occurs when instead of RK3 a higher order method, namely RK4 is used. Fairly clear and explicit figures concerning the stability regions corresponding to various RK, AB and AM methods are available for instance in [15] p.215.
The solution to BBM (3) when $g\left(x\right):=\mathrm{exp}\left({x}^{2}/2\right)$ is depicted in Fig. 2. The scaling factor has been chosen equal to $1/3.$ In this way we enlarge the span of evolution of initial data. A second reason for a scaling factor less than unity.
In order to implement the RungeKutta method we have used the MATLAB builtin function ode45. In spite of the large integration time no floating point overflow has been observed. The whole computation involves some seconds and 284 time steps on our HP xw8400 workstation with clock speed of 3.2 Ghz.
The evolution of the time steps length is illustrated in Fig. 3.
From physical point of view, it appears from Fig. 2 that, as time proceeds, new solitary waves have emerged from the initial profile. More exactly, the first wave emerged immediately from the initial one and the fourth one still had considerable evolution to undergo before it could reasonably be said to have established its asymptotic form. This development over time is perfectly identical with that described in [2] p. 265. In this paper the initial front has the form $\mathrm{exp}\left({x}^{2}/10\right).$
It is also fairly important to mention that the numerical value of the ”energy” functional $E\left(u\right)$ oscillates with largest errors of order $O\left({10}^{4}\right)$ during the integration time (see Section 2.3). For the above initial front we have the estimation ${E}_{0}=2.65868077.$
We have also carried out numerical experiments with the following initial fronts: $\mathrm{arcsin}\left(\frac{1}{1+{x}^{2}}\right),\mathrm{arcsin}\left(\frac{2x}{1+{x}^{2}}\right),\leftx1\right\mathrm{exp}\left(\leftx\right\right)$ and $\sqrt[3]{{\left(x+1\right)}^{2}}\sqrt[3]{{\left(x1\right)}^{2}}.$ All these functions are continuous ${L}_{2}$ ones but not continuously differentiable. They, along with their first derivatives, approach $0$ as $x\to \pm \mathrm{\infty}.$ Thus they do not satisfy the first requirement of smoothness for the initial data imposed in the existence theorem quoted in Introduction. In spite of this lack of smoothness, for each front above, a solution to problem (9) and consequently to the problem (3) has been produced. As natural, all these solutions inherit the defect of smoothness of their initial fronts.
A MATLAB animation can be easily set up in order to show the breaking up of the initial front into subsequent solitary waves.
2.2 The SC method
In order to apply this method we use equidistant nodes with spacing $h,$ symmetric with respect the origin, namely
$${x}_{k}:=\left(k\frac{N+1}{2}\right)h,k=1,\mathrm{\dots},N.$$ 
The interpolant for a real valued function $f$ reads
$${s}_{N}\left(x\right):=\sum _{j=1}^{N}{\varphi}_{j}\left(x\right){f}_{j},$$ 
where the sinc or ”Whittaker cardinal” functions are
$${\varphi}_{j}\left(x\right)=\frac{\mathrm{sin}\left(\pi \left(x{x}_{j}\right)/h\right)}{\pi \left(x{x}_{j}\right)/h}.$$ 
Like the Hermite method, the sinc method contains a free parameter, namely the step size $h.$ The typical estimate for $h$ is $h=C/\sqrt{N}$ for some constant $C$ independent of $N.$ Unfortunately the information provided in [17] for this constant is rather poor. In this work the author analyzes the application of various methods (Galerkin and collocation) based on sinc functions in solving a large variety of equations. It is interesting to remark that he solves differential boundary value problems only on finite intervals but otherwise with endpoint singularities.
The sinc expansion is the easiest of all pseudospectral methods to program and to understand (see [3], § ”Whittaker Cardinal or ”Sinc” Functions”). Unfortunately, its success is assured only for the so called bandlimited functions. Another weakness, shared with the Hermite functions, is that the sinc series is useful only when it approximates a function $f\left(x\right)$ which decays exponentially as $\leftx\right\to \mathrm{\infty}.$ Unlike the Hermite functions, they do not provide the exact solution to any classical eigenvalue problem.
The sinc differentiation matrices are available in [20]. The linearized spatial discretization operator attached to (9) in this case is defined by
$${L}_{S}:={D}_{S}^{\left(1\right)}{\left(I{D}_{S}^{\left(2\right)}\right)}^{1}.$$  (11) 
From the last two columns of Table 1 it is apparent that ${L}_{H}$ is better conditioned and more normal than ${L}_{S}.$ The pseudospectrum for ${L}_{S}$ is depicted in Fig. 4. In order to advocate the versatility of MoL for BBM type equations we overlapped this pseudospectrum with the stability region AB of order 3.
2.3 Conservation of the energy integral
It is fairly important to observe to what extent various finite difference schemes along with the scaled HC method generate a MoL which conserves the ”energy” ${u}^{2}+{u}_{x}^{2}$ over time, i.e., the integral (5) is indeed an invariant of the physical system.
In Fig. 5 a) for a scaling factor $b=1/3$ and in Fig. 5 b) for $b=1$ the outcomes provided by MATLAB codes ode45, ode23 and ode15s in solving (9) are compared. We mention that the codes ode45, ode23 implement respectively RK4 and RK2 schemes and ode15s implements a variable order solver based on the numerical differentiation formulas (NDF).
The RK4 used along with $b=1/3$ performs the best. Actually, this is the main reason we have worked with $b\le 1$. In our numerical experiments this value of $b$ has been proved to be the most suitable.
We have to emphasize that these errors in conserving the energy, as well as the solution itself strongly depend on the scaling factor. All our numerical experiments revealed the necessity for such a factor to be less than unity. In this way for the outmost scaled Hermite roots we have the estimation
$${x}_{1}={x}_{N}=O\left(\sqrt{N}/b\right)\text{as}N\to \mathrm{\infty}.$$ 
Thus, with $$ we extend the span for the wave propagation and we can follow the wave for larger time intervals.
It is also worth mentioning that our numerical estimations satisfy the inequality (see also Fig. 2)
$$\underset{\begin{array}{c}{x}_{1}\le x\le {x}_{N}\\ t\in [0,T]\end{array}}{sup}\left{u}_{N,b}(x,t)\right\le {\left(E\left(u\right)\right)}^{1/2}\approx 1.6305,$$  (12) 
which is in perfect accordance with the analytical results reported in [1]. The estimation is independent of $N$ and $b$ for $64\le N\le 128$ and $$.
In order to compute the integral $E$ we have used the GaussHermite quadrature rule (see for instance [12], § 3.5 ”Quadrature: Gauss–Hermite Formula”). The term ${u}_{x}^{2}$ in the integrand was estimated using the first order Hermite differentiation matrix ${D}_{H}^{\left(1\right)}$.
Unfortunately, for sinc collocation, due to the distribution of the collocation nodes such strategy is rather impossible. This is the main explanation for the fact that in spite of the well defined MoL based on these functions our numerical experiments restrict to the Hermite case. Roughly speaking, it seems that the evolution of some large waves can not be detected using bandlimited functions.
3 The MoL for the modified BBM initial value problem
In this case we have solved the following Cauchy’s problem
$$\{\begin{array}{c}(I{D}_{H}^{\left(2\right)}){U}_{t}={D}_{H}^{\left(1\right)}(U+\frac{1}{n+1}U{.}^{n+1}),t\in [0,T],\\ U\left(0\right)=g\left(X\right),\end{array}$$  (13) 
with the same notations as in (9).
The linearizations of this problem look identical with (10) and (11) respectively. An usual computation (multiplication with $u$ and integration by parts) shows that $E\left(u\right)$ must be conserved over time independent of $n$ in the BBM equation. For $n=2$ the situation is illustrated in Fig. 6 and comparing with Fig. 5 it is clear that the conservation process is at least as accurate as that in BBM problem.
The solution of the modified BBM problem with this value of $n$ looks fairly similar with that depicted in Fig. 2, i.e., certain classes of initial data evolve into a sequence of solitary waves followed by a dispersive wave train. Thus we skip such a picture. In Fig. 6 we use the same notations as in Fig. 5.
One important remark about the linearized operators is well suited at this moment. The conditioning of ${L}_{H}^{B}$ is increasing with the parameter $\alpha $ and the order of approximation $N$ (see the left panel in Fig. 7). This is a fairly intuitive aspect. With respect to the normality of this operator the situation is reversed, namely the Henrici number of ${L}_{H}^{B}$ slightly decreases when $\alpha $ as well as the order of approximation $N$ increases for a fixed value of the scaling parameter $b$. The latter aspect could partially explain the success of the HC method.
4 The MoL for BBMB initial value problem
It is again convenient to write the (1) equation in the following form
$$\left(1{\partial}_{x}^{2}\right){u}_{t}={\partial}_{x}\left(u+\frac{1}{2}{u}^{2}\alpha {\partial}_{x}u\right).$$  (14) 
The MoL based on the HC for to the initial value problem attached to the equation (14) reads
$$\{\begin{array}{c}(I{D}_{H}^{\left(2\right)}){U}_{t}={D}_{H}^{\left(1\right)}(U+\frac{1}{2}U{.}^{2})+\alpha {D}_{H}^{\left(2\right)}U,t\in [0,T],\\ U\left(0\right)=g\left(X\right).\end{array}$$  (15) 
The corresponding linearized spatial discretization operator for (15) is now
$${L}_{H}^{B}:=\left({D}_{H}^{\left(1\right)}+\alpha {D}_{H}^{\left(2\right)}\right){\left(I{D}_{H}^{\left(2\right)}\right)}^{1}.$$ 
The spectral properties of ${L}_{H}^{B}$ and ${L}_{S}^{B}$ depending on the diffusion parameter are summarized in Table 2. At a first glance it is clear that for $\alpha $ decreasing the properties of the Hermite operator are better than those of the sinc one.
$\alpha $  $1.0$  $0.01$  

Method  Hermite  sinc  Hermite  sinc 
Conditioning  1.64e+01  1.20e+00  4.64e+01  9.29e+01 
Henrici  9.54e02  2.79e02  2.27e02  8.29e01 
In order to validate the MoL for stiff problems we have depicted in Fig. 8 the pseudospectrum of ${L}_{H}^{B}$ corresponding to $\alpha =0.01.$ Thus we can conclude that MoL remains stable.
The same conclusion remains valid when the problem (15) is solved by SC coupled with RungeKutta higher than $2$ or AB of order $3$. The pseudospectrum for ${L}_{S}^{B}$ is provided in Fig. 9.
The numerical solution for BBMB problem (15) with the same initial data $g\left(x\right)$ as above is depicted in Fig. 10.
It is rather difficult to compare Fig. 2 and Fig. 10 in order to grasp the meaning of the diffusive term $\alpha {u}_{xx},\alpha =0.01$ on the wave propagation. However, a breaking up phenomenon of the initial wave persists despite the lost of the solution smoothness. From numerical point of view the stiffness introduced by this term makes the time marching process more slow.
We conclude this section with two remarks. Firstly, the spectra of the discrete operators ${L}_{H},{L}_{S},{L}_{H}^{B}$ and ${L}_{S}^{B}$ are almost rectilinear sets of points lying on the ordinate axis, or very close to it, and well included in the interval $[i,i].$ But this interval is the stability region for the leap frog scheme. Thus, this scheme is a possible candidate for a stable MoL.
Secondly, the above strategy extends directly to the more general equations of the form
$${u}_{t}+{u}_{x}+{u}^{n}{u}_{x}{\left(Pu\right)}_{t}=0,$$ 
where the linear operator $P$ belongs to a class of pseudodifferential operators and to BBM equation supplied with a forcing term or with slightly modified asymptotic conditions for the initial front.
5 Concluding remarks
In a sense, this study has succeeded in illustrating numerically the well established analytical and physical work [1]. Thus we have carried out solutions to a class of BBM Cauchy’s problems which inherit the smoothness properties of the initial front, conserve the energy invariant, are uniformly bounded over time and numerically stable.
We have used neither the domain truncation nor the periodicity hypothesis as is usually done. Instead, we have exploited the advantages offered by the scaling factor attached to HC method. Thus, with an ”optimal” scaling factor the accuracy at which the energy integral is satisfied, at least for MoL based on scaled HC and RK4, has been considerably improved. The spectral properties of the linearized spatial discretization operator have been tuned up in order to justify the stability of the MoL algorithm. Eventually the spacetime domain has been enlarged such that the wave evolution is plainly observable. Additionally, the time integration period has been arbitrarily large which is fairly important from practical point of view.
Consequently, we honestly believe that we have added two important examples which confirm the fact that the Trefethen’s rule of thumb provides accurate predictions at least in the case of quasi normal matrices. We also believe that we have provided a solid argument for the validity of the BBM equation as a physical approximation of the dispersion of long waves.
References
 [1] Benjamin, T.B., Bona, J.L., Mahony, J.J., Model Equations for Long Waves in Nonlinear Dispersive Systems, Philos. Trans. Roy. Soc. London Ser.A., 272(1972) 4778
 [2] Bona, J.L., Pritchard, W.G., Scott, L.R., Comparison of Two Model Equations for Long Waves, AMS Lectures in Applied Mathematics, 20(1983) 235267
 [3] Boyd, J.P., Chebyshev and Fourier Spectral Methods, 2nd Ed., DOVER, 2000
 [4] D., Dutykh, T., Katsaounis, D., Mitsotakis, Finite volume schemes for dispersive wave propagation and runup, J. Comput. Phys., 230 (2011), 30353061
 [5] Funaro, D., Polynomial Approximation of Differential Equations, Springer Verlag, Berlin, 1992
 [6] Gheorghiu, C.I.: Laguerre collocation solutions to boundary layer type problems. Numer. Algor. 64, 385401 (2013)
 [7] Gheorghiu, C.I.: Pseudospectral Solutions to some Singular Nonlinear BVPs. Applications in Nonlinear Mechanics. Numer. Algor. doi: 10.1007/s110750149834z
 [8] Gheorghiu, C.I., Spectral Methods for NonStandard Eigenvalue Problems. Fluid and Structural Mechanics and Beyond, Springer Briefs in Mathematics, 2014
 [9] Higham, D.J., Trefethen, L.N., Stifness of ODEs, BIT, 33(1993) 285303
 [10] Kassam, AK., Trefethen, L.N., FourthOrder Time Stepping for Stiff PDEs, SIAM J. Sci. Comp., 26 (2005) 1214–1233.
 [11] Marshall, H.G., Boyd, J.P., Solitons in a Continuously Stratified Equatorial Ocean, J. Phys. Oceanogr., 17 (1987) 10161031
 [12] Olver, F. W. J., Lozier, D. M., Boisvert, R. F., Clark, C. W., eds. , NIST Handbook of Mathematical Functions, Cambridge University Press, 2010
 [13] Omrani, K., The convergence of fully discrete Galerkin approximations for the BenjaminBonaMahony (BBM) equation, Appl. Math. Comput., 180 (2006) 614621
 [14] Omrani, K., Ayadi, M., Finite Difference Discretization of the BenjaminBonaMahonyBurgers Equation, Numer. Methods Partial Differential Eq., 24 (2008) 239248
 [15] Quarteroni, A., Saleri, F., Scientific Computing with MATLAB and Octave, 2nd Ed., Springer, 2006
 [16] Reddy, S.C., Trefethen, L.N., Stability of the method of lines, Numer. Math., 62(1992) 235267
 [17] Stenger, F., Summary of Sinc numerical methods, J. Comput. Appl. Math., 121(2000) 379420
 [18] Trefethen, L.N., Spectral Methods in MATLAB, SIAM, 2000
 [19] Wahlbin, L., Error estimates for a Galerkin method for a class of model equations for long waves, Numer. Math. 23 (1975) 289303
 [20] Weideman, J.A.C., Reddy, S. C., A MATLAB Differentiation Matrix Suite, ACM Trans. on Math. Software, 26 (2000) 465–519
[1] T.B. Benjamin, J.L. Bona, J.J. Mahony, Model equations for long waves in nonlinear dispersive systems, Philos. Trans. Roy. Soc. London Ser.A. 272 (1972) 47–78.
[2] J.L. Bona, W.G. Pritchard, L.R. Scott, Comparison of two model equations for long waves, AMS Lect. Appl. Math. 20 (1983) 235–267.
[3] J.P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd, DOVER, 2000.
[4] D. Dutykh, T. Katsaounis, D. Mitsotakis, Finite volume schemes for dispersive wave propagation and runup, J. Comput. Phys. 230 (2011) 3035–3061.
[5] D. Funaro, Polynomial Approximation of Differential Equations, Springer Verlag, Berlin, 1992.
[6] C.I. Gheorghiu, Laguerre collocation solutions to boundary layer type problems, Numer. Algor. 64 (2013) 385–401.
[7] C.I. Gheorghiu, Pseudospectral solutions to some singular nonlinear BVPs, Appl. Nonlinear Mech. Numer. Algorithm. (2015), doi:10.1007/s110750149834z
[8] C.I. Gheorghiu, Spectral Methods for NonStandard Eigenvalue Problems. Fluid and Structural Mechanics and Beyond, Springer Briefs in Mathematics, Springer, Cham Heidelberg New York Dordrecht London, 2014.
[9] D.J. Higham, L.N. Trefethen, Stifness of ODEs, BIT 33 (1993) 285–303.
[10] A.K. Kassam, L.N. Trefethen, Fourthorder time stepping for stiff PDEs, SIAM J. Sci. Comp. 26 (2005) 1214–1233.
[11] H.G. Marshall, J.P. Boyd, Solitons in a continuously stratified equatorial ocean, J. Phys. Oceanogr. 17 (1987) 10161031.
[12] F.W.J. Olver, D.M. Lozier, R.F. Boisvert, C.W. Clark (Eds.), NIST Handbook of Mathematical Functions, Cambridge University Press, 2010.
[13] K. Omrani, The convergence of fully discrete Galerkin approximations for the BenjaminBonaMahony (BBM) equation, Appl. Math. Comput. 180 (2006) 614–621.
[14] K. Omrani, M. Ayadi, Finite difference discretization of the BenjaminBonaMahonyBurgers Equation, Numer. Methods Partial Differential Eq. 24 (2008) 239–248.
[15] A. Quarteroni, F. Saleri, Scientific Computing with MATLAB and Octave, 2nd Ed., Springer, SpringerVerlag Berlin Heidelberg, 2006.
[16] S.C. Reddy, L.N. Trefethen, Stability of the method of lines, Numer. Math. 62 (1992) 235–267.
[17] F. Stenger, Summary of Sinc numerical methods, J. Comput. Appl. Math. 121 (2000) 379–420.
[18] L.N. Trefethen, Spectral methods in MATLAB, SIAM (2000).
[19] L. Wahlbin, Error estimates for a Galerkin method for a class of model equations for long waves, Numer. Math. 23 (1975) 289–303.
[20] J.A.C. Weideman, S.C. Reddy, A MATLAB differentiation matrix suite, ACM Trans. Math. Software 26 (2000) 465–519.