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

whose solution $u(x,t)$ is considered in a class of non-periodic functions defined on and $t\ge 0.$

For $\alpha =0$ and $n=1$ the equation (1) is called Benjamin Bona Mahony or the regularized long-wave equation. It has been introduced by the above mentioned mathematicians in the early ’70 (see [1]) as a counterpart of the well known Korteweg-de Vries (KdV) equation

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 Mahony-Burgers (BBM-B 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 long-wave 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 Crank-Nicolson-type 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

where the initial data $g\left(x\right)$ satisfies the conditions

1. 1.

   $g\in W_2^1\left(ℝ\right)\cap C^2\left(ℝ\right),$

2. 2.

Then the problem (3) has a solution $u\in {𝒞}_{\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.,

Additionally, $E_0$ is conserved over time, i.e., $E\left(u\right)=E_0$ where

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(ℝ\right)\cap C^1\left(ℝ\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 ${𝒞}_T$ of continuous functions $u(x,t)$ on $ℝ\times [0,T)$ equipped with the maximum norm and its subspace ${𝒞}_T^{l,m}$ of functions $u$ such that ${\partial }_x^i{\partial }_t^{l}u\in {𝒞}_T$ for $0\le i\le l$ and $0\le j\le m$ equipped with the subsequent norm. The space $W_2^1\left(ℝ\right)$ is the usual Sobolev space of $L_2$ functions over $ℝ$ 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 non-periodic 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 Runge-Kutta (RK) methods of orders $2$, $3$ and $4,$ and linear multistep methods such as Adams-Bashforth (AB for short), Adams-Moulton (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 Lax-stability 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 Runge-Kutta 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 BBM-B 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 time-discretization 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.

In order to solve problem (3) we rewrite the differential equation in the following convenient form

In this case we have solved the following Cauchy’s problem

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.

It is again convenient to write the (1) equation in the following form

The MoL based on the HC for to the initial value problem attached to the equation (14) reads

The corresponding linearized spatial discretization operator for (15) is now

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.

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 Runge-Kutta higher than $2$ or AB of order $3$. The pseudospectrum for $L_S^B$ is provided in Fig. 9.

The numerical solution for BBM-B 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

where the linear operator $P$ belongs to a class of pseudo-differential operators and to BBM equation supplied with a forcing term or with slightly modified asymptotic conditions for the initial front.

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 space-time 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.