Sinc Collocation Method for Solving the Benjamin-Ono Equation

Abstract

We propose a simple, though powerful, technique for numerical solutions of the Benjamin-Ono equation. This approach is based on a global collocation method using Sinc basis functions. Some properties of the Sinc collocation method required for our subsequent development are given and utilized to reduce the computation of the Benjamin-Ono equation to a system of ordinary differential equations. The propagation of one soliton and the interaction of two solitons are used to validate our numerical method. The method is easy to implement and yields accurate results.

1. Introduction

It is well known that nonlinear partial differential equations (NPDEs) are widely used to describe complex phenomena in various fields of sciences, such as physics, biology, and chemistry.

In this paper, we consider the initial value problem for the Benjamin-Ono (BO) equation of the form 𝑢_𝑡+𝑢𝑢_𝑥+ℋ𝑢_𝑥𝑥=0,(𝑥,𝑡)∈ℝ², (1) together with the following initial and the boundary conditions: 𝑢(𝑥,0)=𝑢₀(𝑥),𝑢,𝑢_𝑥⟶0as|𝑥|⟶∞. (2) Here 𝑢∶ℝ²→ℝ is a real valued function and ℋ denotes the Hilbert transform defined and described by [1] (amongst others): ℋ𝑓𝑥)=1(𝜋𝑃^∞_−∞𝑓𝑥−𝑦𝑦𝑑𝑦, (3) where 𝑃 denotes the Cauchy principal value. This equation was derived by Benjamin [2] and later by Ono [3] as a model for one-dimensional waves in deep water and it has a close relation with the famous KdV equation which models long waves in shallow water [4].

We recall that the BO equation has an infinite sequence of invariants [5], the first three of which are 𝐼₁=^∞_−∞()𝐼𝑢𝑑𝑥mass₂=^∞_−∞𝑢²(),𝐼𝑑𝑥momentum₃=^∞_−∞𝑢³−3𝑢_𝑥ℋ(𝑢)𝑑𝑥energy. (4) The BO has been shown to admit rational analytical soliton solutions in 𝑥 and 𝑡 [6]. However, for a given arbitrary initial condition, finding analytical solutions of the BO equation becomes an intractable problem. Therefore the use of numerical methods plays an important role in the study of the dynamics of the BO equation.

James and Weideman [7] used the Fourier method which implicitly assumes the periodicity of the boundary conditions and a method based on rational approximating function to compute numerical solutions of the BO. The rational method was shown to have spectacular accuracy for solutions that do not wander too far from the origin while, for long dated solutions, the Fourier method was shown to retain superior accuracy. Miloh et al. [8] proposed an efficient pseudospectral method for the numerical solution of the weakly nonlinear Benjamin-Ono equation for arbitrary initial conditions and suggested a practical new relationship for estimating the number of solitons in terms of arbitrary initial conditions. Thomée and Vsaudeva Murthy [9] used the Crank-Nicolson approximation in time and finite difference approximations in space to solve the BO equation. They treated the nonlinear term in a standard conservative fashion and discretized the Hilbert transform by a quadrature formula which was computed efficiently using the fast Fourier transform. Recently, Boyd and Xu [10] compared three spectral methods for solving the Benjamin-Ono equation, Fourier pseudospectral, rational Christov functions, and Gaussian radial basis functions. They highlighted the advantages and the disadvantages of each numerical method. For instance, the Fourier pseudospectral method is very fast through use of the fast Fourier transform (FFT), but requires domain truncation which is unnecessary for a rational basis. Radial basis functions are slow for a given number of grid points; however, they have a very flexible grid adaptation. To the best of our knowledge, there are no records regarding numerical solutions of the BO equation by using Sinc numerical methods.

In this paper, we propose a Sinc collocation method for solving the BO equation. Sinc numerical methods for solving ordinary and partial differential equations have increasingly become popular and have been extensively studied, in particular, for problems on unbounded domains and those having singularity at boundaries. The first Sinc method was introduced by Stenger [11] to solve two-point boundary value problems for second order differential equations. In the field of nonlinear partial differential equations, the use of Sinc numerical methods has shown tremendous potential. For instance, Al-Khaled [12] showed the utility of the Sinc Galerkin method for solving the Korteweg-de Vries equation. Recently, Mokhtari and Mohammadi [13] presented a meshfree method based on the Sinc collocation method to solve the generalized regularized long wave equation. Excellent overviews of the existing Sinc methods for solving ODEs, PDEs, and integral equations can be found in [14, 15].

The layout of this paper is as follows. We describe the formulation of the Sinc collocation method and the discretization of the Hilbert transform in Section 2. Section 3 is devoted to the Sinc collocation discretization of the BO equation. Furthermore, we examine the stability of the method by using a linearized stability analysis. Numerical results illustrating the merits of the new scheme are given in Section 4. Finally, we present our conclusions in Section 5.

2. Sinc Collocation Methods

In this section, we introduce some useful definitions and pertinent theorems of the Sinc function needed for space discretization of the BO equation. We use ℤ to denote the set of all integers, ℝ the set of all real numbers, and ℂ the set of all complex numbers. Let 𝑢 be a function defined on ℝ and ℎ the step-size. Then the Whittaker function is defined by the following series: 𝐶(𝑢,ℎ)(𝑥)=^∞_𝑗=−∞𝑢_𝑗𝑆_𝑗(𝑥), (5) where 𝑆_𝑗 is given as follows: 𝑆_𝑗(𝑥)=Sinc𝑥−𝑗ℎℎ,𝑗=0,±1,±2,…, (6) and the Sinc function defined on the real line ℝ is given by (𝑥)=()Sincsin𝜋𝑥𝜋𝑥,𝑥≠0,1,𝑥=0. (7) Whenever the series given in (5) converges, 𝑢 is approximated by using a finite number of terms. Therefore, for a positive integer 𝑁, (5) implies 𝐶_𝑁(𝑢,ℎ)(𝑥)=^𝑁_𝑗=−𝑁𝑢_𝑗𝑆_𝑗(𝑥). (8) The approximation of the 𝑛th derivative of the function 𝑢 by the Sinc is given by 𝑢^(𝑛)(𝑥)≈^𝑁_𝑗=−𝑁𝑢𝑑𝑗𝑘^𝑛𝑑𝑥^𝑛𝑆_𝑗(𝑥). (9) The Sinc collocation method requires derivatives of Sinc functions to be evaluated at collocation points. We can obtain explicit close form expressions for these; namely, 𝛿⁽⁰⁾_𝑗𝑘=[𝑆_𝑗(𝑥)]_{𝑥=𝑥𝑘}=𝛿1,𝑗=𝑘0,𝑗≠𝑘._(1)𝑗𝑘=𝑑𝑑𝑥[𝑆_𝑗(𝑥)]_{𝑥=𝑥𝑘}=()0,𝑗=𝑘,−1^𝑗−𝑘ℎ𝛿𝑗−𝑘,𝑗≠𝑘._(2)𝑗𝑘=𝑑²𝑑𝑥²[𝑆_𝑗(𝑥)]_{𝑥=𝑥𝑘}=−𝜋²3(),𝑗=𝑘,−2−1^𝑗−𝑘ℎ²𝑗−𝑘²,𝑗≠𝑘. (10)

Lemma 1. The Hilbert transform of a function 𝑢 can be approximated by ℋ[𝑢𝑥)=1](𝜋𝑃^∞_−∞𝑓𝑥−𝑦𝑦𝑑𝑦≈^𝑁_𝑗=−𝑁𝑢_𝑗𝑇_𝑗(𝑥), (11) where 𝑇_𝑗(𝑥)=1−cos𝜋/ℎ𝑥−𝑗ℎ.𝜋/ℎ𝑥−𝑗ℎ (12) The entries of the Hilbert transform matrix at collocation point 𝑥_𝑘=𝑘ℎ are given by 𝐻_𝑗𝑘=𝑇_𝑗𝑥_𝑘=()0,𝑗=𝑘,1−−1^𝑗−𝑘𝜋𝑗−𝑘,𝑗≠𝑘. (13)

In order to state the convergence theorem of the Sinc collocation method, we introduce the following notation and definitions.

Definition 2. Let 𝒟_𝑑 denote the infinite strip region with width 2𝑑 (𝑑>0) in the complex plane 𝒟_𝑑≡{𝑧∈𝐶∣|Imz|<𝑑}. (14) For all 0<𝜀<1, let 𝒟_𝑑(𝜀) be defined by 𝒟_𝑑(𝜀)=||<1𝑧∈𝐶Rez𝜀,|Imz|<𝑑(1−𝜀). (15) Let 𝐻¹(𝒟_𝑑) be the Hardy space over the region 𝒟_𝑑, that is, the set of functions such that lim_{𝜀→0𝜕𝒟𝑑(𝜀)}|𝑢(𝑧)||𝑑𝑧|<∞. (16)

One then gets the following convergence results.

Theorem 3 (see [15]). Assume, with positive constants 𝛼, 𝛽, and 𝑑, that
(1) 𝑢 belongs to 𝐻¹(𝒟_𝑑);
(2) 𝑢 decays exponentially on the real line; that is, |𝑢(𝑥)|≤𝛼exp−𝛽|𝑥|,∀𝑥∈ℝ. (17) Then we have sup_{−∞<𝑥<∞}𝑢^(𝑛)(𝑥)−^𝑁_𝑗=−𝑁𝑢_𝑗𝑆_𝑗^(𝑛)(𝑥)≤𝐶₁𝑁^(𝑛+1)/2exp−𝜋𝑑𝛽𝑁,sup_{−∞<𝑥<∞}ℋ[𝑢](𝑥)−^𝑁_𝑗=−𝑁𝑢_𝑗𝑇_𝑗(𝑥)≤𝐶₂𝑁^1/2exp−𝜋𝑑𝛽𝑁, (18) for some constants 𝐶₁ and 𝐶₂, where the mesh size ℎ is taken as ℎ=𝜋𝑑𝛽𝑁^1/2. (19)

The proof of the above Theorem 3 is beyond the scope of this paper. We refer the reader to [15] for a detailed discussion. Theorem 3 states that, if 𝑢 is an analytic function on an infinite strip containing the real line and allows specific growth restrictions, then it has exponentially decaying absolute errors in the Sinc approximation. Similarly, the approximation of the derivatives and the Hilbert transform of 𝑢 decay exponentially.

3. Construction of the Method

We consider the initial value problem for the Benjamin-Ono equation of the form 𝑢_𝑡+𝑢𝑢_𝑥+ℋ𝑢_𝑥𝑥=0,(𝑥,𝑡)∈ℝ², (20) together with the following initial: 𝑢(𝑥,0)=𝑢₀(𝑥), (21) and the boundary conditions 𝑢−𝑢_𝑥=0as𝑥⟶−∞,𝑢+𝑢_𝑥=0as𝑥⟶∞. (22) Here 𝑢∶ℝ²→ℝ is a real valued function and ℋ denotes the Hilbert transform defined by ℋ𝑓𝑥)=1(𝜋𝑃^∞_−∞𝑓𝑥−𝑦𝑦𝑑𝑦, (23) where 𝑃 denotes the Cauchy principal value.

In order to propose an unconditionally stable methodology, we discretize the time derivative of the BO equation using a classic finite difference formula by the 𝜃-weighted (0≤𝜃≤1) scheme between two successive levels 𝑛 and 𝑛+1 as 𝑢^𝑛+1−𝑢^𝑛𝛿𝑡+𝜃𝑢𝑢_𝑥^𝑛+1+𝜃ℋ𝑢^𝑛+1_𝑥𝑥+(1−𝜃)𝑢𝑢_𝑥^𝑛+1+(1−𝜃)ℋ𝑢^𝑛+1_𝑥𝑥=0, (24) where 𝑢^𝑛=𝑢(𝑥,𝑡^𝑛), 𝑡^𝑛=𝑡^𝑛−1+𝛿𝑡, and 𝛿𝑡 is the time step size.

The nonlinear term (𝑢𝑢_𝑥)^𝑛+1 in (24) can be approximated by Taylor expansion. We obtain the following formula: 𝑢^𝑛+1≈𝑢^𝑛+𝛿𝑡𝑢_𝑡^𝑛≈𝑢^𝑛𝑢+𝛿𝑡^𝑛+1−𝑢^𝑛𝛿𝑡+𝒪𝛿𝑡²,𝑢_𝑥^𝑛+1≈𝑢^𝑛_𝑥𝑢+𝛿𝑡^𝑥𝑛+1−𝑢^𝑛_𝑥𝛿𝑡+𝒪𝛿𝑡². (25) Therefore 𝑢𝑢_𝑥^𝑛+1≈𝑢𝑢_𝑥^𝑛+𝛿𝑡𝑢_𝑡^𝑛𝑢^𝑛_𝑥+𝑢^𝑛𝑢^𝑛_𝑥𝑡+𝒪𝛿𝑡²=𝑢𝑢_𝑥^𝑛𝑢+𝛿𝑡^𝑥𝑛+1−𝑢^𝑛_𝑥𝑢𝛿𝑡^𝑛𝑥+𝑢^𝑛𝑢_𝑥^𝑛+1−𝑢^𝑛_𝑥𝛿𝑡+𝒪𝛿𝑡²,=𝑢^𝑛𝑢_𝑥^𝑛+1+𝑢^𝑛_𝑥𝑢^𝑛+1−𝑢^𝑛𝑢^𝑛_𝑥+𝒪𝛿𝑡². (26) Equation (24) can be rewritten as 𝑢^𝑛+1+𝛿𝑡𝜃𝑢^𝑛𝑢_𝑥^𝑛+1+𝑢^𝑛_𝑥𝑢^𝑛+1+ℋ𝑢^𝑛+1_𝑥𝑥=𝑢_𝑛−𝛿(1−2𝜃)𝑢^𝑛𝑢^𝑛_𝑥+(1−𝜃)ℋ𝑢^𝑛_𝑥𝑥. (27) Now the space variable is discretized into collocation points 𝑥₁=𝑎,…,𝑥_𝑖=𝑎+(𝑖−1)ℎ,…,𝑥_𝑁,||=𝑏ℎ=𝑏−𝑎,𝑁−1 (28) so that the solution of (20) is interpolated and approximated by means of the Sinc functions as 𝑢𝑥,𝑡^𝑛=𝑢^𝑛(𝑥)=^𝑁_𝑗=1𝑢^𝑛_𝑗𝜉_𝑗(𝑥), (29) where the basis functions 𝜉_𝑖(𝑥) are given by 𝜉_𝑗(𝑥)=𝜔₀(𝑥)𝑆,𝑗=1_𝑗(𝑥)𝜔,𝑗=2,…,𝑁−1,₁(𝑥),𝑗=𝑁. (30) Here 𝑆_𝑗(𝑥) is defined as 𝑆_𝑗(𝑥)=Sinc𝑥−𝑗−1ℎ−𝑎ℎ, (31) and the two cubic interpolants [14] 𝜔₀=(𝑥−𝑎)(𝑏−𝑥)²(𝑏−𝑎)²+(2𝑥+𝑏−3𝑎)(𝑏−𝑥)²(𝑏−𝑎)³,𝜔₁=(−2𝑥+3𝑏−𝑎)(𝑥−𝑎)²(𝑏−𝑎)³+(𝑏−𝑥)(𝑥−𝑎)²(𝑏−𝑎)² (32) are introduced to handle the radiation boundary conditions in (22).

The unknown parameters 𝑢_𝑗 in (29) are to be determined by collocation method. Therefore, for each collocation point 𝑥_𝑖, (29) can be written as 𝑢^𝑛𝑥_𝑖=^𝑁_𝑗=1𝑢^𝑛_𝑗𝜉_𝑗𝑥_𝑖,𝑖=1,…,𝑁. (33) The substitution of (33) into (27) at collocation points 𝑥_𝑖 yields the following equation: ^𝑁_𝑗=1𝑢_𝑗^𝑛+1𝜉_𝑗𝑥_𝑖+𝛿𝑡𝜃^𝑁_𝑗=1𝑢^𝑛_𝑗𝜉_𝑗𝑥_𝑖^𝑁_𝑗=1𝑢_𝑗^𝑛+1𝜉^′_𝑗𝑥_𝑖+^𝑁_𝑗=1𝑢^𝑛_𝑗𝜉^′_𝑗𝑥_𝑖^𝑁_𝑗=1𝑢^𝑛_𝑗𝜉_𝑗𝑥_𝑖+^𝑁_𝑖=1^𝑁_𝑗=1𝑢_𝑗^𝑛+1𝜉^′′𝑥_𝑖𝑇_𝑖𝑥_𝑘=^𝑁_𝑗=1𝑢^𝑛_𝑗𝜉_𝑗𝑥_𝑖−𝛿𝑡(1−2𝜃)^𝑁_𝑗=1𝑢^𝑛_𝑗𝜉_𝑗𝑥_𝑖×^𝑁_𝑗=1𝑢^𝑛_𝑗𝜉^′_𝑗𝑥_𝑖+(1−𝜃)^𝑁_𝑖=1^𝑁_𝑗=1𝑢^𝑛_𝑗𝜉^′′𝑥_𝑖𝑇_𝑖𝑥_𝑘. (34) To obtain the matrix representation of the expression (34), we introduce the following matrix and vector notations: 𝑢^𝑛=𝑢^𝑛₁,𝑢^𝑛₂,…,𝑢^𝑛_𝑁^𝑇,𝐈⁽⁰⁾=𝜉_𝑗𝑥_𝑖𝐈,𝑖,𝑗=1,…,𝑁,₍₁₎=𝜉^′_𝑗𝑥_𝑖𝐈,𝑖,𝑗=1,…,𝑁,₍₂₎=𝜉_𝑗^′′𝑥_𝑖𝑇,𝑖,𝑗=1,…,𝑁,𝐇=_𝑗𝑥_𝑖,𝑖,𝑗=1,…,𝑁. (35) Thus we obtain the following system of 𝑁 linear equations in 𝑁 unknown parameters 𝑢^𝑛+1 which can be expressed in a matrix form 𝐌𝑢^𝑛+1=𝐑𝑢^𝑛, (36) where 𝐌=𝐈⁽⁰⁾+𝛿𝑡𝜃{𝐀+𝐁+𝐂},𝐑=𝐈⁽⁰⁾−𝛿𝑡{(1−2𝜃)𝐃+(1−𝜃)𝐂} (37) in which 𝑢^𝑛_𝑥=𝐈⁽¹⁾𝑢^𝑛,𝐀=𝑢^𝑛∗𝐈⁽¹⁾,𝐁=𝑢^𝑛_𝑥∗𝐈⁽⁰⁾,𝐂=𝐇×𝐈⁽²⁾,𝐃=𝑢^𝑛∗𝑢^𝑛_𝑥∗𝐈⁽⁰⁾, (38) where symbols “∗” and “×” mean the component-wise and the matrix multiplications, respectively.

The stability of the approximation (36) is investigated by using the matrix method. The linearized version of (20) is obtained by assuming the quantity 𝑢 in the nonlinear term 𝑢𝑢_𝑥 as a local constant 𝑢. The error at timestep 𝑛 is defined as 𝑒^𝑛=𝑢^𝑛_exact−𝑢^𝑛_num, (39) where 𝑢^𝑛_exact and 𝑢^𝑛_num are the exact and the numerical solutions, respectively. The error matrix equation of the linearized BO equation is given by [𝐈+𝜃𝛿𝑡𝐉]𝑒^𝑛+1=[𝐈+(1−𝜃)𝛿𝑡𝐉]𝑒^𝑛, (40) where 𝐈=𝐈⁽⁰⁾ and 𝐉=𝑢𝐈⁽¹⁾+𝐇×𝐈⁽²⁾. Now, let 𝐏=[𝐈+𝜃𝛿𝑡𝐉]⁻¹[𝐈+(1−𝜃)𝛿𝑡𝐉]; then (40) becomes 𝑒^𝑛+1=𝐏𝑒^𝑛. (41) The numerical scheme is stable if ‖𝐏‖₂⩽1, or equivalently, 𝜌(𝐏)⩽1; that is, the spectral radius of the matrix 𝐏 has to be smaller than or equal to one. From (40) it can be seen that the stability conditions are possible if all the eigenvalues of the matrix [𝐈+𝜃𝛿𝑡𝐉]⁻¹[𝐈+(1−𝜃)𝛿𝑡𝐉] satisfy the following condition: 𝜆_𝐼−𝛿𝑡(1−𝜃)𝜆_𝐽𝜆_𝐼+𝛿𝑡𝜃𝜆_𝐽⩽1, (42) where 𝜆_𝐼 and 𝜆_𝐽 are the eigenvalues of the matrices 𝐼 and 𝐽, respectively.

For 𝜃=0.5, the relation (42) becomes 𝜆_𝐼−0.5𝛿𝑡𝜆_𝐽𝜆_𝐼+0.5𝛿𝑡𝜆_𝐽⩽1. (43) If both 𝜆_𝐼 and 𝜆_𝐽 are complex eigenvalues, we write 𝜆_𝐼=𝑎_ℎ+𝑖𝑏_ℎ and 𝜆_𝐽=𝑎_𝑘+𝑖𝑏_𝑘. In this case, (43) can be rewritten as 𝑎_ℎ−0.5𝛿𝑡𝑎_𝑘+𝑖𝑏_ℎ−0.5𝛿𝑡𝑏_𝑘𝑎_ℎ+0.5𝛿𝑡𝑎_𝑘+𝑖𝑏_ℎ+0.5𝛿𝑡𝑏_𝑘⩽1. (44) The inequality (44) is satisfied if 𝑎_ℎ𝑎_𝑘+𝑏_ℎ𝑏_𝑘⩾0. For real eigenvalues, the inequality (44) holds true for (𝜆_ℎ⩾0 and 𝜆_ℎ⩾0) or (𝜆_ℎ⩽0 and 𝜆_ℎ⩽0). This shows that scheme (36) is unconditionally stable if 𝑎_ℎ𝑎_𝑘+𝑏_ℎ𝑏_𝑘⩾0 for complex eigenvalues and (𝜆_ℎ⩾0 and 𝜆_ℎ⩾0) or (𝜆_ℎ⩽0 and 𝜆_ℎ⩽0) for real eigenvalues.

For 𝜃=0, the inequality (42) becomes 1−𝛿𝑡𝜆_𝐽𝜆_𝐼⩽1. (45) In other words, 𝛿𝑡⩽2𝜆_𝐽𝜆_𝐼,𝜆_𝐽𝜆_𝐼⩾0. (46) Therefore, when 𝜃=0, the scheme (42) is conditionally stable.

4. Numerical Results

The analytical soliton solutions of the BO equation are rational functions in 𝑥 and 𝑡. For example, one and two soliton solutions are given, respectively, by (see [6, 16, 17]) 𝑢(𝑥,𝑡)=4𝑘₁𝑘²₁𝜃²₁,+1𝑢()=𝑥,𝑡4𝑘₁𝑘₂𝑘₁𝜃²₁+𝑘₂𝜃²₂+𝑘₁⁻¹𝑘₂⁻¹𝑘₁+𝑘₂³𝑘₁−𝑘₂⁻²𝑘₁𝑘₂−𝑘₁+𝑘₂²𝑘₁−𝑘₂⁻²²+𝑘₁𝜃₁+𝑘₂𝜃₂², (47)(48) where 𝜃₁=𝑥−𝑘₁𝑡+𝜙₁, 𝜃₂=𝑥−𝑘₂𝑡+𝜙₂, and 𝑘₁, 𝑘₂, 𝜙₁, 𝜙₂ are arbitrary constants. In this section, we examine the proposed algorithm using different test problems related to the propagation of one soliton and interaction of two soliton solutions. The accuracy of the scheme is measured by using the following error norms: ‖𝐸‖₂=‖𝑢−𝑢‖₂=ℎ^𝑁_𝑗=1𝑢_𝑗−𝑢_𝑗^21/2,‖𝐿‖_∞=‖𝑢−𝑢‖_∞=max_{1⩽𝑗⩽𝑁}𝑢_𝑗−𝑢_𝑗, (49) where 𝑢 and 𝑢 represent the exact and approximate solutions, respectively, and ℎ is the minimum distance between any two points of set points for which the errors are evaluated. In our computational work, we use the collocation points 𝑥₁=𝑎,…,𝑥_𝑖=𝑎+(𝑖−1)ℎ,…,𝑥_𝑁||=𝑏,ℎ=𝑏−𝑎.𝑁−1 (50) The BO equation possesses infinite conservation laws [5]; the first three are given as follows: 𝐼₁=^𝑏_𝑎𝑢𝑑𝑥,𝐼₂=^𝑏_𝑎𝑢²𝐼𝑑𝑥,₃=^𝑏_𝑎𝑢³−3𝑢_𝑥ℋ(𝑢)𝑑𝑥, (51) related to the mass, momentum, and energy. The quantities 𝐼₁, 𝐼₂, and 𝐼₃ are applied to measure the conservation properties of the collocation scheme calculated by 𝐼₁≈ℎ^𝑁_𝑗=1𝑢^𝑛_𝑗,𝐼₂≈ℎ^𝑁_𝑗=1𝑢^𝑛_𝑗²,𝐼₃≈ℎ^𝑁_𝑗=1𝑢^𝑛_𝑗³−3𝑢_𝑥^𝑛_𝑗ℋ𝑢^𝑛_𝑗. (52)

4.1. Propagation of Single Solitons

In this experiment, we consider the propagation of single solitons of the BO equation. The initial and boundary conditions are extracted from (47) and (22), respectively. The values of the parameters used in the numerical experiments are 𝑘₁=0.3, 𝜙₁=0, 𝜃=1/2, the number of collocation points 𝑁=200, and 𝛿𝑡=0.001. The soliton from which the initial condition is extracted (47) moves to the right across the space interval 𝑥∈[−100,100] when the time interval is 𝑡∈[0,110]. The numerical tests for this case are performed using the Sinc basis functions.

The error norms 𝐿_∞, 𝐿₂ and conservation quantities 𝐼₁, 𝐼₂, and 𝐼₃ are computed, which are shown in Table 1. From the numerical results given in Table 1 it is observed that, throughout the simulation, the error norms 𝐿_∞ remain less than 1.9803×10⁻⁴ and 𝐿₂ remains less than 9.3644×10⁻⁴, whereas the changes of the invariant 𝐼₁ are less than 10⁻⁴ and the invariants 𝐼₂ and 𝐼₃ at a given time 𝑡 are equal to those of the initial value; our scheme is satisfactorily conservative. The numerical solutions at different time levels (a) and the absolute error (b) at time 𝑡=110 are shown in Figure 1 for 𝑘₁=0.3, 𝜙₁=0, 𝜃=1/2, the number of collocation points 𝑁=200, and 𝛿𝑡=0.001 in the region 𝑥∈[−100,100].

Table 1 

Invariants and errors for a single soliton of the BO equation. 𝑘₁=0.3, 𝜙₁=0, 𝑁=200, 𝛿𝑡=0.001, 𝜃=1/2, and 𝑥∈[−100,100].

(a)

(b)

Figure 1 

Propagation of soliton solutions (a) and absolute error (b) at 𝑡=110 of the BO equation 𝑘1=0.3, 𝜙1=0, 𝑁=200, 𝛿𝑡=0.001, 𝜃=1/2, and 𝑥∈[−100,100].

In the next experiment, we investigate the convergence of our approach in terms of the number of grid points 𝑁 at time 𝑡=1. To ensure that our scheme is not dominated by the time discretization, we choose 𝛿𝑡=10⁻⁵. We take the parameters such as 𝑘₁=0.4, 𝑘₂=0.1, 𝜙₁=30, 𝜙₂=0, 𝛿𝑡=10⁻⁵, 𝜃=1/2, and 𝑥∈[−100,100].

Figure 2 shows that the error norms 𝐿_∞ and 𝐿₂ converge exponentially for the propagation of one soliton solutions and for the interaction of two soliton solutions as we increase the number of grid points 𝑁.

(a)

(b)

Figure 2 

Errors for the propagation of one soliton solutions (a) and for the interaction of two soliton solutions (b) at 𝑡=1 of the BO equation with 𝑘1=0.4, 𝑘2=0.1, 𝜙1=30, 𝜙2=0, 𝛿𝑡=10−5, 𝜃=1/2, and 𝑥∈[−100,100].

4.2. Interaction of Two Solitons

Our second experiment pertains to the interaction of two soliton solutions of the BO equation having different amplitudes and travelling in the same direction. The initial and boundary conditions are extracted from (47) and (22).

To allow the interaction to occur, the experiment was run from 𝑡=0 to 210 in the region [−100,100]. Figure 3 shows the interaction of two soliton solutions of BO equation for 𝑘₁=0.4, 𝑘₂=0.1, 𝜙₁=30, 𝜙₂=0, 𝑁=200, 𝛿𝑡=0.001, 𝜃=1/2, and 𝑥∈[−100,100]. It can be seen that the faster pulse interacts with and emerges ahead of the lower pulse with the shape and velocity of each soliton retained. Because of the nature of the Sinc discretization, the Hilbert transform and matrix differentiation are explicitly obtained with exponential accuracy. The result in Figure 3 shows numerical interaction of the two soliton solution and its time evolution of the 𝐿₂ error from 𝑡=0 to 𝑡=210. In Figure 3(b), it is shown that the Sinc discretization maintains accurate results comparable to those given in [7, 10] for a long period. Therefore, a number of properties of the Benjamin-Ono equation can be checked numerically, including the exhibition of no phase shifts after interaction of the solitons.

(a)

(b)

Figure 3 

Time evolution of the interaction of two soliton solutions (a) and 𝐿2 error (b) of the BO equation. 𝑘1=0.4, 𝑘2=0.1, 𝜙1=30, 𝜙2=0, 𝑁=200 and 300, 𝛿𝑡=0.001, 𝜃=1/2, and 𝑥∈[−100,100].

Numerical check on the conservation mass, momentum, and energy shows that the three quantities remain constant with respect to time as shown in Table 2. Propagation of the single solitary wave and two soliton integrations are simulated well with the proposed algorithms and conservation invariants do not change much during the computer run. Thus Sinc functions can be used to construct approximate numerical methods of the BO equation.

Table 2 

Invariants and errors for interaction of two solitons of the BO equation. 𝑘₁=0.4, 𝑘₂=0.1, 𝜙₁=30, 𝜙₂=0, 𝑁=200, 𝛿𝑡=0.001, 𝜃=1/2, and 𝑥∈[−100,100].

5. Conclusion

A numerical technique based on the Sinc collocation method has been presented for numerical solutions of the BO equation. The efficiency of the method is tested on the problems of propagation of a single soliton as well as interaction of two solitons. The accuracy of solutions is examined in terms of the 𝐿_∞, 𝐿₂ error norms, and the conservation quantities 𝐼₁, 𝐼₂, and 𝐼₃. Stability analysis is performed by the matrix method. The results obtained for the Sinc collocation are very close to analytical ones. Our algorithm was found to be unconditionally stable and exponentially convergent in space and a reliable numerical method for solving the BO equation. The method is easy to implement and compares favourably to Fourier-based spectral methods without any assumptions on periodicity.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

Edson Pindza is thankful to Brad Welch for the financial support from RidgeCape Capital.