A Convergent Crank-Nicolson Galerkin Scheme for the Benjamin-Ono Equation

In this paper we prove the convergence of a Crank-Nicolson type Galerkin finite element scheme for the initial value problem associated to the Benjamin-Ono equation. The proof is based on a recent result for a similar discrete scheme for the Korteweg-de Vries equation and utilizes a local smoothing effect to bound the $H^{1/2}$-norm of the approximations locally. This enables us to show that the scheme converges strongly in $L^{2}(0,T;L^{2}_{\text{loc}}(\mathbb{R}))$ to a weak solution of the equation for initial data in $L^{2}(\mathbb{R})$ and some $T>0$. Finally we illustrate the method with some numerical examples.


Introduction
In this paper we consider a fully discrete Crank-Nicolson Galerkin scheme for the Cauchy problem associated to the Benjamin-Ono (BO) equation, which reads where T > 0 is fixed, u : R × [0, T ) → R is the unknown, u 0 is the initial data and H is the Hilbert transform defined by Hu(x, ·) := p.v. 1 π R u(x − y, ·) y dy, where p.v. denotes the Cauchy principal value. The equation was derived independently by Benjamin [3] and Ono [26], and serves as a model equation for weakly nonlinear long waves with weak nonlocal dispersion. There has been done much work on the well-posedness of (1) and improvement of regularity restrictions on the initial data, we mention [21], [1] and [28], where in the latter global well-posedness for initial data in H s for s ≥ 1 was proved using a gauge transform resembling the famous Cole-Hopf transform for the viscous Burgers equation. By refining this transform, [20] extended the result to s ≥ 0.
The BO equation is formally completely integrable [13,23], a property shared by the well-known Korteweg-de Vries (KdV) equation. The integrability is closely related to the fact that the BO equation admits an infinite number of conserved quantities [1,23] and a Lax pair [4], hence it can be formulated as a Hamiltonian system. Another feature of the BO equation is the existence of families of explicit solutions called solitons [2], and these localized, solitary wave solutions are a consequence of a delicate balance between dispersion and nonlinear convection. Having to approximate the simultaneous appearance of the two aforementioned effects makes the task of finding reliable numerical methods for the BO equation rather challenging, and we mention some of the various methods which have been proposed. In [29] the authors present a finite difference scheme for which they prove error estimates for smooth solutions of the periodic version of the BO equation, while in [9] they consider operator splitting methods of Godunov and Strang type for (1) and prove corresponding L 2 -convergence rates given sufficiently regular initial data. On the other hand, in [27,6] the authors consider Fourier spectral methods and prove error estimates for sufficiently regular solutions of the periodic equation. A high order hybrid finite element-spectral scheme for the Benjamin equation, for which the BO equation is a special case, is presented in [8] where the authors give experimental convergence rates for the method.
However, we emphasize that the goal of this paper is not to present an efficient, high-order numerical method, but rather to formulate a discrete scheme for which one is able to prove existence of a sequence of approximations that converge locally to a weak solution of the BO equation for low regularity initial data, namely functions in L 2 (R). Thus, the scheme differs from the previously mentioned papers in that it can be used as a constructive proof of the existence of solutions to (1). In this respect, the paper at hand is more in the spirit of [10] where the authors present finite difference schemes for (1) and its periodic counterpart which are proved to converge to to classical solutions given initial data in H 2 .
We will here consider a Crank-Nicolson type Galerkin scheme for finding weak solutions to the BO equation based on a method for the KdV equation due to Dutta and Risebro [12] which can be seen as a generalization to higher order temporal approximations of [11], where the implicit Euler method is used instead of Crank-Nicolson for the temporal discretization. The motivation for using these rather simple time integrators is that they are easier to analyze compared to multi-step integrators such as higher order Runge-Kutta methods. Our strategy for establishing convergence will follow closely that of the above papers, in particular by using a local smoothing effect inherent to the equation, but some more work is required here to treat the dispersive term which contrary to the case of KdV is nonlocal for the BO equation. In fact, it is exactly this added challenging aspect of the nonlocal dispersion which is our main motivation for applying the methods of [11,12] to (1). The smoothing effect of the BO equation is also weaker than that of KdV, which combined with the nonlocal nature of the Hilbert transform makes it natural to consider fractional Sobolev spaces, hence making our estimates more involved than in the case of KdV. We note that the results presented in the current paper is in full based on [15]. For a treatment of convergence rates for the scheme discussed here given sufficiently regular initial data the reader is referred to [14].
The paper is structured as follows. In Section 2 we establish some preliminary technical results, e.g. the local smoothing effect mentioned above for a semidiscretized weak formulation of (1). The fully discrete scheme is presented in Section 3, where we prove existence and uniqueness of its solutions for each time step. The main part of the paper is contained in Section 4, where we first prove that the solutions of the fully discrete scheme in Section 3 exhibit the local smoothing effect from Section 2. Then we move on to our main result, Theorem 4.1, which establishes the existence of a sequence of approximate functions which converge locally in L 2 to a weak solution of (1). Finally, in Section 5 some numerical experiments are presented to illustrate the discrete scheme.

Preliminary estimates
In the following we will give a brief explanation of our strategy. Let us momentarily define a weak solution to the Benjamin-Ono equation (1) to be a function where ·, · denotes the standard L 2 -inner product. We now discretize the equation in time using a Crank-Nicolson method. Let ∆t be the time step size, u n ≈ u(·, n∆t) and u n+1/2 := (u n+1 + u n )/2. For a given u 0 ∈ L 2 (R), define u n to be the solution of for all v ∈ H 2 (R) and n ≥ 0. Assuming that the above equation has a unique solution u n+1 we may choose v = u n+1 + u n in (3) which yields Here we have used that the Hilbert transform is antisymmetric, which is a consequence of the following lemma.
Lemma 2.1. The Hilbert transform is a linear operator with the following properties: (i) (Skew-symmetricity) Assume f ∈ L p (R) for 1 < p < ∞ and g ∈ L q (R) where 1/p + 1/q = 1. Then we have (ii) (Commutes with differentiation) For a differentiable function f we have (iii) (L 2 -isometry) The transform preserves the L 2 -norm, Remark 2.1. Setting p = 2 and g = f in (i) in Lemma 2.1 gives Hf, f = 0, which is the aforementioned antisymmetry property. Also, combining (ii) and (iii) shows that the transform is in fact isometric for all Sobolev norms H k with k ∈ {0} ∪ N.
These properties will be essential in our convergence analysis and for a proof of the above lemma we refer to [18, p. 317]. Note that throughout this paper C will denote various positive constants which exact value is of no importance to the arguments. Likewise, C(R) will denote various positive constants depending on the parameter R and so on. Now, in [12] they introduce a smooth, positive and non-decreasing cut-off function ϕ and use integration by parts to derive a local smoothing effect which bounds u n locally in H 1 -norm. This technique is originally due to Kato [22], and is in fact a consequence of the commutator identity  [17,16], and the identity relevant for us is found in [17, p. 227] and reads where D β , β > 0 denotes the homogeneous fractional derivative defined by and R 1/2 (ϕ) is some remainder operator. Here and in the following we use standard notation for the Fourier transform, defined by for a suitable function f . Note the sign change on the left-hand side of (5) compared to that of Ginibre and Velo which is due to their use of the Hilbert transform with opposite sign. According to equation (40) in [17, p. 228] we have where |||·||| denotes the operator norm in L 2 (R). Using (5) and (6) we are able to bound the H 1/2 -norm of u n locally, where H s (R) is the real ordered Sobolev space of functions u such that (1 + |ξ| 2 ) s/2û ∈ L 2 (R) with the corresponding norm u H s (R) = (1 + |ξ| 2 ) s/2û L 2 (R) . Let us define a smooth cut-off function ϕ ∈ C ∞ (R) satisfying: . Properties (a)-(d) are easily achievable by standard mollifier methods. The existence of functions satisfying property (e) can be motivated as follows. If this were not the case we could have started by defining a nonnegative function h such that h(x) = 1 for |x| < R, h(x) = 0 for |x| ≥ R + 1 and 0 ≤ h(x) ≤ 1. Then h 2 ∈ C ∞ c (R) has the same properties as h and by defining ϕ(x) := 1 + √ ϕ x is smooth and compactly supported by definition. Due to the properties of ϕ, v = ϕu n+1/2 is also an admissible test function in H 2 (R), and we will write w := u n+1/2 to save space. Inserting this in (3) and using the identity which is easily attained from integration by parts, we have We rewrite the second term on the left-hand side as where we have used (5). To show that the last term is bounded we use the following fact: Given f ∈ C N (R) and f (k) ∈ L 1 (R) for 0 ≤ k ≤ N we have for some suitable C depending on N and f (k) L 1 (R) for 0 ≤ k ≤ N . This estimate is fairly standard and easily obtained using properties of Fourier transforms of differentiated functions, see e.g. [18, p. 109], and the straightforward inequal- As ϕ xx belongs to C ∞ c (R) and in particular C 2 c (R), we have ϕ (2+k) ∈ L 1 (R) for k = 0, 1, 2. According to (8) we then have | ϕ xx (ξ)| ≤ C (1+|ξ|) 2 and thus ϕ xx L 1 (R) ≤ 2C. Then it follows from (6) that Next we want to estimate the term stemming from the nonlinearity in terms of R |D 1/2 w| 2 ϕ x dx, and the following results will be of use. From Theorem 6.5 in [7] we have for s ∈ (0, 1) and p ∈ [1, ∞) such that sp < n there exists a positive constant C = C(n, p, s) such that the Sobolev space W s,p (R n ) is continuously embedded in L q (R n ) for any q ∈ [p, p * ] with p * := np/(n − sp), where W s,2 (R n ) = H s (R n ).
Next we have the interpolation inequality stated in Proposition 3.1 of [25]: If Now we turn to the third term on the left-hand side of (7) and estimate it as The second inequality above is an application of (9) with n = 1, p = 2, s = 1 4 , and q = p * = 4, while the third inequality comes from (10) with n = 1 s = 1 4 , s 1 = 0, s 2 = 1 2 , and θ = 1 2 . For the first term in the last line we set h = √ ϕ x and use the commutator identity where S 1/2 (h) is a remainder operator. From Proposition 2.1 in [16] we have where |||·||| denotes the operator norm in L 2 (R). This can be estimated as Noting that h ∈ C ∞ c (R) and in particular C 3 c (R), so that h (k) ∈ L 1 (R) for k = 0, 1, 2, 3 and using (8) we have the estimate Thus, taking the L 2 -norm we obtain for some constant C S depending on ϕ x . Inserting the above estimates in (7) we obtain which again implies where we in the last term have used that the L 2 -norm of w = u n+1/2 is bounded by the norm of u 0 . By dropping the positive second term on the left-hand side, summing from n = 0 to n = m − 1 and utilizing that this is a telescoping sum we obtain Also, first summing and then dropping 1 Together these estimates imply that given initial data u 0 ∈ L 2 (R) we have which shows that the solutions of the Crank-Nicolson temporal discretized equation also exhibit the local smoothing effect of the BO equation, as the above sequence space is a temporally discrete analogue of L 2 [0, . Since this smoothing is the main ingredient of the convergence proof in the case of KdV, we want to show that it is present also in our fully discretized element scheme presented in the next section. When formulating the scheme we follow [12] in using test functions of the form ϕv, where v belongs to some finite element space. This makes (7) hold and leads to a H 1/2 -bound like the one obtained here. The problem with this form of the scheme is that one loses the a priori preservation of the L 2norm that was obtained in (4) by directly choosing the test function u n+1/2 . To overcome this difficulty we make use of a CFL condition combined with a majorizing differential equation.

Formulation of the discrete scheme
Here we formulate the Crank-Nicolson type Galerkin scheme under consideration. First we present some remarks on notation and the discretization of time and space. Then we use the weak formulation of the problem and a Crank-Nicolson temporal discretization to define a sequence of functions approximating the exact solution at each discrete time step. We also define an iteration scheme to solve the implicit equation for each time step and show that this has a solution.
3.1. Notation and discretization. We start by partitioning the real line in equally sized elements in the form of intervals. First define the grid points x j = j∆x for j ∈ Z, where ∆x is the spatial discretization parameter or step length. Then the elements can be written as I j = [x j−1 , x j ]. Now turn to the discretization of the time interval considered. Given a fixed time horizon T > 0 and a temporal discretization parameter ∆t we set t n = n∆t for n ∈ {0, 1, . . . , N }, where N + 1 2 ∆t = T . For convenience we also use the notation t n+1/2 = (t n + t n+1 ) /2. Let ϕ be defined as in the previous section. We define the weighted L 2 -inner product u, v ϕ := u, vϕ , and the associated weighted norm u 2,ϕ = u, u ϕ .
3.2. Galerkin scheme. As always for the finite element method we start by deriving a weak formulation of the problem (1), like the one obtained in (2). Applying the Crank-Nicolson temporal discretization to the weak formulation gives (3). Instead of looking for solutions to this equation in H 2 (R) we will look for solutions belonging to a finite-dimensional subspace S ∆x of this Hilbert space. We define the subspace S ∆x as follows; assuming r ≥ 2 is a fixed integer we denote the space of polynomials on the interval I of degree not exceeding r by P r (I). Our goal is to find an approximation u ∆x to the solution of (1) which for all t ∈ [0, T ] belongs to Now define P to be the L 2 -orthogonal projection onto S ∆x . Then we define the sequence {u n } n≥0 through the following procedure: Given u 0 = Pu 0 , find u n+1 ∈ S ∆x which satisfies for all v ∈ S ∆x and n ∈ {0, 1, . . . , N }. Clearly, (13) is an implicit scheme and consequently one must solve a nonlinear equation to obtain u n+1 from u n . The procedure for solving this equation at each time step is described in the following subsection. Note also that u 0 L 2 (R) ≤ u 0 L 2 (R) , and thus from here on we will always use the L 2 -norm of the initial data u 0 as an upper bound for the norm of the approximation u 0 .
The following inverse inequalities presented in [5, p. 142] will be instrumental in our later estimates.
where the constants C 1 , C 2 > 0 are independent of z and ∆x. Note that the leftmost inequality also holds for z instead of z x .

Solvability for one time step.
To show the existence of a solution u n for each time step we define the iteration scheme which is to hold for all test functions v ∈ S ∆x . The existence of a unique solution w ℓ+1 to (15) is guaranteed by noting that one may consider this a Galerkin scheme for a linear problem involving a bilinear form in the variables w ℓ+1 and v. Using the commutator estimate for the part of the bilinear form involving the Hilbert transform, and choosing ∆t small enough, C∆t ≤ 1 2 say, the bilinear form is coercive with respect to the L 2 -norm, which implies positive definiteness of the resulting matrix system.
We now present a lemma that guarantees the solvability of the implicit scheme (13), and the technique for showing this is due to Simon Laumer (private communication).
Lemma 3.1. Choose a constant L such that 0 < L < 1 and set We consider the iteration scheme (15) and assume that the following CFL condition holds, where C 2 is defined in (14) and λ is given by where ∆t is taken sufficiently small. Then there exists a function u n+1 which solves (13) and lim ℓ→∞ w ℓ = u n+1 . In addition, Proof of Lemma 3.1. We start by rewriting (15) as Hu n x , (ϕv) x . From the above equation one derives For the term involving the Hilbert transform we estimate as before and use the fact that ϕ ≥ 1 to obtain We then estimate the term A 2 by repeatedly applying Young's inequality and using (14), which yields where estimates analogous to the preceding ones lead to Collecting the bounds we have the following inequality for ℓ ≥ 1, Assuming ∆t small enough that 3 2 − C∆t ≥ 1 we obtain w ℓ+1 − w ℓ 2 2,ϕ ≤ 2C 2 λ 2 max{ w ℓ 2 2,ϕ , w ℓ−1 2 2,ϕ , u n 2 2,ϕ } w ℓ − w ℓ−1 2 2,ϕ . (19) We will now bound w 1 , and so by setting ℓ = 0 in (15) we get Choosing v = u n +w 1 2 gives Estimating the term involving the Hilbert transform as before, using Young's inequality and (14) leads to Choosing ∆t small enough that 1 4 − C∆t 2 ≥ 1 8 then gives Now we claim that the following holds for ℓ ≥ 1, The proof follows an induction argument. From (20) and (16) we get and so (21c) and (21b) hold for ℓ = 1. Setting ℓ = 1 in (19) while using (16) we obtain which shows that (21a) holds for ℓ = 1. Now assume that (21a) and (21b) hold for ℓ = 1, . . . , m. One then has thus (21b) holds for all ℓ. This result together with (19) and (16) lead to This shows that (21a) holds for all ℓ as well. Since 0 < L < 1 this shows that {w ℓ } is Cauchy and hence converges, which completes the proof of Lemma 3.1.

Convergence of the scheme
In this section we will prove the convergence of the scheme introduced in the previous section. As mentioned earlier we will use a local smoothing effect of the BO equation to obtain a H 1/2 loc (R) estimate of the approximations. We begin with the following important lemma.
Lemma 4.1. Let λ, K and L be defined as in Lemma 3.1 and let u n be the solution of the scheme (13). Assume also that ∆t satisfies for some Y which only depends on u 0 L 2 (R) . 1 Then there exist a positive time T and a constant C, both depending only on u 0 L 2 (R) such that for all n satisfying n∆t ≤ T the following estimate holds In addition, the approximation u n satisfies the following H 1/2 -estimate ∆t (n+ 1 2 )∆t≤T Proof of Lemma 4.1. Starting with (13) it follows that (7) holds, and with the same estimates as in the associated section and the fact that 1 ≤ ϕ(x) ≤ 2+2R we obtain the inequality R u n+1 2 ϕ dx + 2∆t Dropping the term involving the fractional derivative and writing a n = R (u n ) 2 ϕ dx then gives a n+1 ≤ a n + ∆tf (a n+ 1 2 ), with the function f (a) = C a + a 2 .
It is easily seen that a n+ 1 2 ≤ (a n + a n+1 )/2 and so {a n } solves the Crank-Nicolson method for the differential inequality Let us then consider the following ordinary differential equation where K comes from Lemma 3.1. It is not difficult to show that this differential equation has a unique solution y which blows up at some finite time T ∞ depending only on the initial condition, and so we choose T = T ∞ /2. Note that the solution of this ordinary differential equation is strictly increasing and convex. We now compare this solution with (26) under the assumption that the CFL condition (22) holds with Y := y(T ). We claim that a n ≤ y(t n ) for all n ≥ 0, and argue by induction. As y(0) = a 0 , the claim holds for n = 0. Now assume that it holds for n ∈ {0, 1, . . . , m}. As 0 ≤ a m ≤ y(T ), (22) implies that (16) holds, and thus Lemma 3.1 gives a m+1 ≤ K 2 a m . Therefore we have The convexity of f then gives a m+1 ≤ a m + ∆tf which proves the claim. Since ϕ ≥ 1 we get the L 2 -stability estimate (23), This proves (24) and completes the proof of Lemma 4.1.

4.1.
Bounds on temporal derivative. We will here obtain bounds on the temporal derivative to be used later in the analysis. The following lemma will be of use.

Lemma 4.2.
Let ψ ∈ C ∞ c (−R, R) and ϕ be defined by properties (a)-(e) in Section 2. Then there exists a projection P : In addition, P satisfies the bounds where the constant C is independent of ∆x.
Proof of Lemma 4.2. The proof is a straightforward adaptation of the L 2 -projection results found in the monograph of Ciarlet [5, p. 146].
In our upcoming estimates we also need the following Sobolev inequality.
From the definitions of the dual norms in H −2 (R) and The above relations together with Lemma 4.2 is used to prove the following lemma regarding the boundedness of the temporal derivatives of the approximate solutions.
where D + t u n is the forward time difference operator Proof of Lemma 4.3. Start by rewriting (13) as which holds for all v ∈ S ∆x . Let ψ ∈ C ∞ c (−R, R) and set v = P (ψ), where P is the projection from Lemma 4.2, to get Using (29), (28) and (23) we estimate the first term on the right-hand side as follows The second term can be estimated as where we have used (23), (28) and the L 2 -isometry of the Hilbert transform. Together this gives , R , and the estimate is proven.
If ψ ∈ C ∞ c (R) then P (ψ) ∈ S ∆x . By the exact same arguments as above, but this time on R instead of [−R, R], we get 4.2. Convergence to a weak solution. Prior to stating our theorem of convergence we define the weak solution of the Cauchy problem (1) in the following way.
Definition 4.1. Let Q > 0 and u 0 ∈ L 2 (R). Then u ∈ L 2 0, T ; ). Now we define the approximate solution u ∆x ∈ S ∆x , which will be shown to converge to a weak solution of (1), by the interpolation formula (36) We then have the following convergence theorem, which is the main result of the paper.
Theorem 4.1. Let {u n } n∈N be a sequence of functions defined by the scheme (13) and assume that u 0 L 2 (R) is finite. Assume furthermore that ∆t = O(∆x 2 ). Then there exist a positive time T and a constant C, depending only on R and u 0 L 2 (R) such that where u ∆x is given by (36). Moreover, there exists a sequence {∆x j } ∞ j=1 and a function u ∈ L 2 0, T ; as ∆x j j→∞ −−−→ 0. The function u is a weak solution of the Cauchy problem for (1), which is to say that it satisfies (35) with Q = R − 1.
Remark 4.1. The convergent subsequence u ∆xj in Theorem 4.1 can be used as a constructive proof of existence of solutions to (1), as noted in Section 1. On the other hand, owing to the well-posedness for initial data u 0 ∈ L 2 (R) [20] we can in fact conclude that the whole sequence converges as ∆x → 0.
Next we have where in the last inequality we have used that u 0 due to Sobolev embedding. In the last line we apply (24) to the sum, and the inverse equality (14) combined with the assumption ∆t ≤ C∆x 2 to the second term to conclude that (38) holds. For the third relation, note that Using (32) it is then easy to show that (39) holds. If one instead uses (34) and considers the whole real line R, the same argument leads to Also, as ϕ is a positive and bounded smooth function, (37) and (38) give Based on the bounds (42a), (42b) and (39) we apply the Aubin-Simon compactness lemma [19,Lemma 4.4] to the set {ϕu ∆x } to prove the existence of a sequence {∆x j } j∈N such that ∆x j → 0, and a functionũ such that ϕu ∆xj →ũ strongly in L 2 (0, T ; L 2 ([−R, R])), as j approaches infinity. As ϕ ≥ 1, (43) implies that there exists u such that (40) holds. The strong convergence allows passage to the limit in the nonlinearity. Now it remains to prove that u indeed is a weak solution of (1). In the following we consider the standard L 2 -projection of a function ψ with k + 1 continuous derivatives into the space S ∆x for some k ∈ N 0 , denoted by P. That is, Now let v ∆x = Pv and observe that we may write Now we estimate the preceding terms. From (30), (44) and (39) we obtain From (23), (29) and (44) we get The next terms may be rewritten as Using (30), (32), (14) and (23) we have the estimate Similarly we get Furthermore we have Thus these terms can be estimated like the preceding ones, and as ∆t/ √ ∆x = O(∆x) we obtain Using the L 2 -isometry of the Hilbert transform, (23) and (44) we obtain Finally, from (31), the H 2 -isometry of the Hilbert transform and (41) we have the estimate We combine the preceding estimates to conclude that (45) holds. Also, observe that by passing ∆x → 0 we obtain T )) and integrate by parts to conclude that (35) holds, that is This concludes the proof of Theorem 4.1.

Numerical experiments
Here we present some numerical experiments to illustrate the fully discrete scheme (13). We stress that this is not an exposition of a numerical method for the problem on the real line in the usual sense as these typically are based on applying a numerical scheme for the periodic version of the problem on a sufficiently large domain such that the reference solutions to are close to zero outside of it, see e.g. the experimental sections of [29,10,8]. As this paper concerns the convergence of the discretized real line problem only, we find such an approach irrelevant to the main result presented here. In our illustration we will also consider a discretized domain which is large enough for the reference solutions to be close to zero outside of it, but we use the full line Hilbert transform where we simply neglect the contributions from outside the domain since the reference solution is nearly zero there. We still have to close our system of equations and this we do by imposing a periodic boundary condition by requiring the approximation to take the same value at both ends of the domain. A detailed presentation of our setting will follow.
Inspired by [12] we define S ∆x in the following way. Let f and g be the functions For j ∈ Z we define the basis functions where x j = j∆x. Then {v j } M −M spans a 4M + 2 dimensional subspace of H 2 (R). In the following we define N := 2M , which is the number of elements used in the approximation. We have chosen the weight function ϕ to be 120 + x so that it is a positive, increasing function for which the derivative equals 1 on the domain considered, which in this case is the interval [−100, 100].
In our experiments we have chosen to set ∆t = O(∆x) contrary to the assertion ∆t = O(∆x 2 ) from the theory, as smaller time steps did not lead to significant improvement in the convergence rates of the approximations. In the iteration (15) to obtain u n+1 we chose the stopping condition w ℓ+1 − w ℓ L 2 ≤ 0.002∆x u n L 2 , which typically required 4-7 iterations for the cruder discretisations and 2-3 iterations for the finer ones. On each element, H(v j ) x , (v i ) x was evaluated by applying a seven point Gauss-Legendre (GL) quadrature rule to the principal value integral defining H(v j ) x in eight GL points, followed by applying the eight point GL rule to the inner product. To avoid problems near the singular point, the principal value integral was evaluated using a subtracted Hilbert transform method [24, p. 685].
For t = n∆t, we set u ∆x (x, t) = u n (x, t) = M j=−M u n j v j (x). We have measured the relative error of the numerical approximation compared to the exact solution u, where the L 2 -norms were computed with the trapezoidal rule in the points x j of the finest grid under consideration. As mentioned in the introduction, the BO equation admits an infinite number of conserved quantities [1], and the first three in this hierarchy are typically denoted mass, momentum and energy, respectively. Numerical methods which preserve more of the conserved quantities for completely integrable partial differential equations have been observed to give more accurate approximations than the ones preserving fewer. Therefore we have computed the relative change I n (u ∆x (t)) := (Q n (u ∆x (t)) − Q n (u ∆x (0)))/Q n (u ∆x (0)) for n = 1, 2, 3 to see how well these quantities were conserved in our discretization. We have also included the rate of convergence for E, defined for each intermediate step between element numbers N 1 and N 2 as ln(E(N 1 )) − ln(E(N 2 )) ln(N 2 ) − ln(N 1 ) , where we now see E as a function of the element number N . In our example we consider a solution to (1) presented in [29], namely u s2 (x, t) = 4c 1 c 2 c 1 λ 2 1 + c 2 λ 2 2 + (c 1 + c 2 ) 2 c −1 where λ 1 := x − c 1 t − d 1 and λ 2 := x − c 2 t − d 2 . When c 2 > c 1 and d 1 > d 2 this expression represents a tall soliton overtaking a smaller one while moving to the right. We applied the fully discrete scheme with initial data u 0 (x) = u s2 (x, 0) and parameters c 1 = 0.3, c 2 = 0.6, d 1 = −30 and d 2 = −55. The time step was set to ∆t = 0.5 u 0 −1 L ∞ ∆x and the numerical solutions were computed for t = 90 and t = 180, that is during and after the taller soliton overtaking the smaller one. To approximate the full line problem we set the domain to [−100, 100] with the aforementioned periodic boundary condition. The results are presented in Table 1 and a comparison between the approximation for N = 256 and the exact solution is shown in Figure 1.  The latter shows that the shape and position of the numerical approximation for N = 256 agrees quite well with the exact solution. Still, we observe that for t = 180 there is a visible error in the height of the tallest soliton, which has introduced a small phase error in the approximation. Since the soliton is very narrow this causes a relative error of 30%, as seen from Table 1. This error is actually larger than the error for N = 128 where the approximate solution has pronounced oscillations.  Table 1. Relative L 2 -error, I 1 , I 2 and I 3 at t = 90 and t = 180 for initial data u s2 .
Still, the phase error for N = 128 is much smaller than for N = 256, and since the solitons are so narrow this has a much larger impact on the relative L 2 -norm of the error. Thus an element number of N = 128 or less seems to be too small to get consistent reduction in the errors for increasing N . We see that for N greater than or equal to 256, the relative L 2 -error at both t = 90 and t = 180 decreases, but not in a systematic fashion. Regarding the conserved quantities, we see that these are preserved quite well in the approximation, and increasing the number of elements tends to improve how well they are preserved, most notably for (48) and (49) The lack of systematic reduction for the error is most likely caused by our approximation of the full line problem by restricting to a finite interval whilst still using the full line Hilbert transform. It should also be pointed out that this is a complicated example, as one has to approximate the nonlinear interaction between two passing solitons. Finally, we refer to [14] for a study of theoretical convergence rates for this scheme and a modified scheme aimed at the periodic version of (1), given sufficiently regular initial data.