Numerical approximation of solitary waves of the Benjamin equation

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS Available online at www.sciencedirect.com

Mathematics and Computers in Simulation xxx (2012) xxx–xxx

Original article

Numerical approximation of solitary waves of the Benjamin equation V.A. Dougalis a,b , A. Durán c,d,∗ , D.E. Mitsotakis e a Mathematics Department, University of Athens, 15784 Zographou, Greece Institute of Applied and Computational Mathematics, F.O.R.T.H. 71110, Heraklion, Greece c Department of Applied Mathematics, University of Valladolid, Campus Miguel Delibes, P/Belen 15 E-47011 Valladolid, Spain d IMUVA, Instituto de Investigación en Matemáticas, University of Valladolid, Spain Institute for Mathematics and its Applications, University of Minnesota, 114 Lind Hall, 207 Church Street SE, Minneapolis, MN 55455, United States b

e

Received 20 October 2011; received in revised form 18 April 2012; accepted 10 July 2012

Abstract This paper presents several numerical techniques to generate solitary-wave profiles of the Benjamin equation. The formulation and implementation of the methods focus on some specific points of the problem: on the one hand, the approximation of the nonlocal term is accomplished by Fourier techniques, which determine the spatial discretization used in the experiments. On the other hand, in the numerical continuation procedure suggested by the derivation of the model and already discussed in the literature, several algorithms for solving the nonlinear systems are described and implemented: the Petviashvili method, the Preconditioned Conjugate Gradient Newton method and two Squared-Operator methods. A comparative study of these algorithms is made in the case of the Benjamin equation; Newton’s method combined with Preconditioned Conjugate Gradient techniques, emerges as the most efficient. The resulting numerical profiles are shown to have a high accuracy as travelling-wave solutions when they are used as initial conditions in a time-stepping procedure for the Benjamin equation. The paper also explores the generation of multi-pulse solitary waves. © 2012 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Benjamin equation; Numerical computation of solitary waves; Numerical continuation; Spectral methods

1. Introduction The goal of this paper is the description and comparison of some techniques to compute solitary-wave profiles of the Benjamin equation ut + αux + βuux − γHux − δuxxx = 0,

(1)

∗ Corresponding author at: Department of Applied Mathematics, University of Valladolid, Campus Miguel Delibes, P/Belen 15 E-47011 Valladolid, Spain. E-mail addresses: [email protected], [email protected] (V.A. Dougalis), [email protected] (A. Durán), [email protected] (D.E. Mitsotakis).

0378-4754/$36.00 © 2012 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.matcom.2012.07.008

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 2

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

In (1), α, β, γ, δ are nonnegative constants and H = H∂x , where H denotes the Hilbert transform  ∞ f (y) 1 Hf (x) = P.V. dy. π −∞ x − y The Fourier symbol of the operator H is  = |ξ| u(ξ). (Hu)(ξ) Particular cases of (1) are the Korteweg–de Vries (KdV) equation (γ = 0, δ > 0) and the Benjamin–Ono (BO) equation (γ > 0, δ = 0). Eq. (1) arises in the study of unidirectional propagation of long internal waves of small amplitude at the interface of two ideal fluids (the heavier of which has infinite depth) in the presence of surface tension [5–7,2]. Global well-posedness in L2 for the corresponding initial-value problem is proved in [25] (where other results concerning generalized versions of (1) are also obtained). It is well known that the functionals [6]  ∞ C(u) = (2) u dx, −∞

I(u) =

1 2



∞

−∞

u2 dx,

E(u) = αI(u) +



∞

−∞

(3) 

 β 3 1 1 2 u − γuHu + δux dx, 6 2 2

(4)

are preserved by sufficiently smooth, suitably vanishing at infinity solutions of (1). The quantity (4) determines the Hamiltonian formulation of (1) ut = JδE(u),

J = −∂x .

in suitable function spaces and where δE denotes the variational derivative. Solitary-wave solutions u(x, t) = ϕ(x − cs t), cs > 0 of (1) were initially studied, in terms of existence, stability and asymptotic behaviour, by Benjamin in [5–7]. The profile ϕ satisfies, assuming that ϕ → 0 as |x|→ ∞, the equation F (ϕ, γ) = δE(ϕ) − cs δI(ϕ) = (−cs + α)ϕ +

β 2 ϕ − γHϕ − δϕ = 0. 2

(5)

Benjamin argued that (5) has solutions that exhibit an oscillating behaviour in some parts of the spatial domain and a monotonic decay at infinity as 1/|x|2 . A complete theory of existence and stability of solitary-wave solutions for small values of γ is provided in [2]. Explicit formulas for the profiles (except in the KdV and BO cases) are not known. Several techniques of numerical approximation involve numerical continuation algorithms on some parameter in (1). Numerical continuation is applied in [2] to design a code to approximate the solitary waves. The results reveal some properties, such as the symmetry about their crests and (except in the KdV and BO cases) the existence of a finite number of oscillations. In addition, numerical studies in [22] suggest an inelastic interaction of solitary waves and consequently a nonintegrable character of (1) when γ > 0, δ > 0. Note that if we look for solutions of (5) in the form  α − cs 2(α − cs ) ϕ(X) = − ψ (6) X , β δ then the profile ψ satisfies ˜ = ψ − 2γHψ ˜ Q(ψ, γ) − ψ − ψ2 = 0,

(7)

where γ γ˜ = √ . 2 δ(α − cs )

(8)

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

3

˜ is zero, we obtain the well known solitary-wave solutions of This is the notation used in [2]. When γ (and therefore γ) the KdV equation. Thus, we will consider (5) to be dependent on γ as a continuation parameter. The form (7) will be used when necessary. The literature on the numerical computation of travelling-wave solution profiles of nonlinear wave equations includes many different techniques. Some variational methods [4,15] are based on the minimization of some relevant functional; there are also classical procedures, such as shooting methods [37], fixed point iterations with stabilizing factors or renormalization methods [30,28,1], Newton-like methods [9,11,33,34]. The so-called imaginary-time evolution methods (see [35] and references therein) are based on integrating numerically a related evolution equation, for which the profile is a stationary solution and in which t is replaced by it. The squared-operator methods [36,34] iterate a modified equation whose linearization operator at the profile is the square of that of the original equation, guaranteeing the convergence. Finally, ‘iterative cleaning’ procedures [8,22,16] generate a solitary-wave profile by running a timestepping code with some initial condition, isolating the leading pulse and using it as a new initial condition for the code in an iterative way, taking advantage of the resolution of solutions into solitary waves. In this paper we focus on continuation techniques and we apply them to the Benjamin equation. Section 2 describes two directions that the numerical continuation algorithm can follow. One is continuation by differentiation, assuming a smooth dependence of the solitary wave on the continuation parameter. This leads to the so-called Davidenko equation [31,10]. A second way, which will be considered in this paper, is given by incremental methods. Here, the parameter interval is partitioned, and starting from a value of the continuation parameter for which the solitary-wave profile is known, we compute a spatial approximation to the profile at the next value of the parameter by solving iteratively the nonlinear system (5) for the profiles. This iterative process can be implemented with many alternatives, some of which will be described and applied to the Benjamin equation. Specifically, we will consider the Petviashvili method (as an example of a fixed-point iteration technique), the Preconditioned Conjugate Gradient Newton (CGN) method and the Squared-Operator method (SOM), with its modified version MSOM. For each method we will explain its theoretical formulation and implementation in Section 3. This description will be general, independent of the spatial discretization, although the particular spectral method that is used in the experiments will be emphasized. As the classical treatment of the nonlocal term makes use of the Discrete Fourier Transform (DFT), the computations that illustrate the comparisons between the methods are implemented with Fourier collocation discretization in space. The fast convergence of the spatial approximation makes the errors in space negligible compared with the errors in the corresponding continuation-iteration process. This permits us to establish some conclusions about the behavior of the iterative methods. In Section 4 we consider the computational efficiency of the methods when they are applied to approximate solitary waves of the Benjamin equation for various values of the continuation parameter γ ∈ [0, 1). In the tests we compare the number of iterations and the cpu time and we conclude that, among the methods considered, the CGN is the most efficient. Finally, in Section 5 we study the accuracy of the solitary-wave profiles (computed with CGN) as travelling-wave solutions when used as initial data in a time-stepping code for the Benjamin equation. The generation of multipulse solitary waves is also be considered. 2. Continuation algorithms From now on we assume that (5) is discretized, in some way, on a sufficiently long interval (− l, l) with periodic boundary conditions, giving rise to a system of equations Fk (ϕh , γ) ≡ (−cs + α)ϕh,k +

β (Gh (ϕh ))k − γ(Hh ϕh )k − δ(Dh2 ϕh )k = 0, 2

k = 0, . . . , N − 1.

(9)

In (9): • ϕh = (ϕh,0 , . . ., ϕh,N−1 ) denotes the approximation to the values of ϕ at the grid xj = − l + hj, j = 0, . . ., N − 1, where h = 2l/N. It is assumed that N is even and ϕh is extended as a periodic function ϕh = {ϕj }j∈Z , ϕj+N = ϕj , defined on the extended uniform grid xj = −l + hj, j ∈ Z • Gh (ϕ) is an approximation of the nonlinear term ϕ2 . • Hh and Dh2 are, respectively, discretizations of the operators H and ∂xx . Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 4

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

2.1. Continuation by differentiation If we are interested in the computation of the profile ϕ(γ * ), corresponding to some parameter value γ * , a first way to implement the continuation procedure is by differentiation [31]. Moreover, there is an analytic mapping ϕ :  → H 1 (R) ∩ C∞ (R) γ −→ ϕ(γ), satisfying (5) on some interval  = [0, γ * ), cf. [2]. If γ 0 = 0 then ϕ0 = ϕ(γ 0 ) corresponds to the known KdV profile, and ϕ = ϕ(γ) is a solution of the initial-value problem (ivp) dϕ ∂F = − (ϕ, γ), dγ ∂γ ϕ(0) = ϕ0 ,

F  (ϕ, γ)

γ ∈ (0, γ ∗ ),

(10)

where F is the Frechet derivative with respect to ϕ. The system (10) is sometimes called the Davidenko equation [31,10]. Then, an ODE solver on the interval [0, γ * ) can be used to approximate the discrete version of (10): dϕh ∂Fh =− (ϕh , γ), dγ ∂γ ϕh (0) = ϕ0,h , Fh (ϕh , γ)

γ ∈ (0, γ ∗ ),

where Fh = (Fh,0 , . . ., Fh,N−1 ) is defined in (9). The resolution requires the computation and handling of the Jacobian Fh (ϕh , γ); the difficulties in the implementation will be described later, when treating Newton’s method. 2.2. Incremental methods In the present paper we will consider a second type of continuation technique that consists of defining an homotopic path γ0 = 0 < γ1 < . . . < γM = γ ∗ , with hM = γ * /M for some M and intermediate values γ m = mhM , m = 0, . . ., M. The idea is to start from the known profile ϕ0 = ϕ(γ 0 ), and compute ϕ(γ m ) from ϕ(γ m−1 ) using an iterative method in (5) or the spatially discretized version (9) with ϕ(γ m−1 ) (or an extrapolation from previous values) as the starting iteration. This class of techniques is known as ‘incremental methods’. The literature considers many possibilities in the choice of the iteration procedure. In the case of the Benjamin equation and among the possibilities cited in Section 1, we focused on several examples of two groups of methods: Variants of classical techniques such as fixed-point iteration or Newton’s method, and squared-operator methods, based on considering a related equation and properties of the corresponding linearized operator [34]. 3. Description of some incremental methods and their application to the Benjamin equation 3.1. Petviashvili method Several methods are modifications of the classical fixed-point procedure (which typically diverges for this kind of problems [28]). One of them is the Petviashvili method [30]. This was originally implemented to compute solitary waves of the Kadomtsev–Petviashvili I (KPI) equation; its use to compute stationary and solitary-wave solutions has also been proposed in many other papers (cf. e.g. [28,23] and references therein). The method is a version of the fixed-point iteration with a suitable iteration function. In our case, for a general step of the γ-continuation algorithm, Eq. (9) can be written as the system β Sh ϕh ≡ (−cs + α − γHh − δDh2 )ϕh = − (Gh ϕh ). 2 Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

5

If the operator Sh = (−cs + α) − γHh − δDh2 is invertible, then the approximate solitary wave ϕh satisfies β ϕh = − Sh−1 (Gh ϕh ). 2

(11)

This corresponds to the equation satisfied by the exact profile ϕ β Sϕ = ((−cs + α) − γH − δ∂xx ) ϕ = − ϕ2 . 2

(12)

The Fourier transform allows to write (12) in the form (ξ) = − ϕ

 2 ](ξ) β [ϕ . 2 −cs + α − γ|ξ| + δξ 2

The Petviashvili method incorporates a stabilizing factor to ensure convergence in the following way. If we multiply (12) by ϕ and integrate, we obtain by the Parseval identity

∞

∞ 2 ϕ(ξ)|2 dξ −∞ ((−cs + α − γH − δ∂xx )ϕ)ϕ dx −∞ (−cs + α − γ|ξ| + δξ )|

∞ = 1. m(ϕ) = =

2 ∞  2 ](ξ) ϕ(ξ)dξ −(β/2)[ϕ −∞ −(β/2)ϕ ϕ dx −∞

This suggests considering the stabilizing factor for the approximation and to formulate the Petviashvili method as the iteration procedure for (11)   β −1  [ν+1] [ν]

, ν = 0, 1, . . . ϕh = mν − S h G h ϕ h (13) 2 for initial ϕh[0] , mν = m(ϕh[ν] ) and where the power is a free parameter of great importance for the convergence of the method. We note that in the case of convergence, it is expected that mν → 1 as ν→ ∞. This convergence can be analysed by using the theory explained in [28] (see also [23]). In our case if min|ψ(x)| < x∈R

1 , 2

where ψ is given in (6), then the method is locally convergent for ∈ (1, 3), with the maximum rate of convergence attained when = 2. In addition, a technique to improve the rate of convergence of the Petviashvili method (and, in fact, of any iterative method to find solitary waves) is described in [23,24] and called the mode elimination technique. This generates the so-called accelerated Petviashvili method. 3.1.1. Fourier collocation spatial discretization As mentioned in Section 1, the presence of the nonlocal term in (5) suggests the use of Fourier techniques for its approximation. Since Fourier collocation will be used for the experiments in Sections 4 and 5, the discretization (9) will be expressed in Fourier variables. As a result (13) has the form    (Gh ϕh[ν] )p β 

p − (ϕh[ν+1] )p = m 2 (p))) 2 (−c + α − γ H h (p) + δ(D s h

p ∈ Z,

(14)

where ( ϕh )p denotes the pth discrete Fourier coefficient of ϕh , the stabilizing factor is approximated by using 

p = − m

 [ν] 2 (p)))|(ϕ 2 h (p) + δ(D + α − γH h )p | h ,      [ν] [ν] p∈Z (β/2)(Gh ϕh )p ϕh

p∈Z (−cs

p

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 6

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

2 (p) denote, respectively, the Fourier symbol of the approximations H and D2 . In computational terms h (p), D and H h h h (and using FFT) (14) is formulated as     ϕh[ν] . ∗ ϕh[ν] β

p  [ν+1]

p ϕh , −N/2 ≤ p ≤ N/2, ν = 0, 1, . . . , =− m 2 (−c + α − γ|p| + δp2 ) s p where 1. The product ϕh[ν] . ∗ ϕh[ν] is understood component wise. 2. The approximation to the stabilizing factor has the form 



p = m

−N/2≤p≤N/2

    −cs + α − γ|p| + δp2 | ϕh[ν] |2

β  [ν] [ν]  (ϕ[ν] )p ϕ . ∗ ϕh −N/2≤p≤N/2 p h 2 h



p

,

(the double prime in the sums shows that the first and last terms are divided by two). 3. We took Dh2 = D2 (the square of the pseudospectral differentiation operator) and Hh as the discretization of the  p = |p|ϕˆ p . operator H (with the periodic version of the Hilbert Transform) with discrete Fourier coefficient (Hϕ) Let us mention that spectral techniques have been widely used for discretizing terms in nonlocal partial differential equations involving the Hilbert transform; cf. e.g. [2,10,22,27,29], and their references. Recent work in this direction includes [3,13,14]. 3.2. SOM and MSOM methods Two other techniques to compute solitary-wave solutions of the Benjamin equation use a different point of view. They are the so-called SOM (Squared Operator Method) and its variant MSOM (Modified Squared Operator Method). A complete study of these algorithms and many references may be found in e.g. [36]. Consider Eq. (12) written in the form β L0 ϕ ≡ −Sϕ − ϕ2 = 0. 2 Let M be a real-valued, positive definite Hermitian operator. Instead of integrating numerically the system ϕt = ±M −1 L0 ϕ, to obtain approximations to the equilibrium L0 ϕ = 0, the idea consists of integrating ϕt = −M −1 L∗ M −1 L0 ϕ,

L = −S − βϕ,

(15)

where L is the linearization operator of (12) at ϕ and * denotes the L2 adjoint (in our case L* = L). Use of the explicit Euler method leads to the iterative SOM ϕ[ν+1] = ϕ[ν] − (M −1 LM −1 L0 ϕ)|ϕ=ϕ[ν] τ.

(16)

The operator M is introduced as a preconditioning operator, to accelerate the convergence. Typical choices of M are the linear part of L0 (when it is positive definite, or modified to be positive definite), and operators of the form M = c − ∂xx for some constant c. When M = I, the implementation is simpler but the convergence is very slow. The choice of the Euler method is justified in [36] with the argument that other time-stepping schemes require more computations, making the resulting iteration less efficient. This technique is related to the functional minimization method for Hamiltonian equations, proposed and explained in [18]; see the comments in [36]. Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

7

The implementation of (16) requires splitting the second term on the right-hand side into four terms. Taking into account the form of L0 and L, we have  β ϕ[ν+1] = ϕ[ν] − τ M −1 SM −1 Sϕ[ν] + M −1 SM −1 (ϕh[ν] . ∗ ϕh[ν] ) + βM −1 ϕh[ν] M −1 Sϕh[ν] 2  β2 + M −1 ϕh[ν] M −1 (ϕh[ν] . ∗ ϕh[ν] ) . 2 For the Fourier collocation version used in the numerical experiments we took M = c − ∂xx and the last two terms in the expression in brackets above were computed with the suitable combination of FFT/IFFT. The local convergence of (16) is ensured as long as τ is below a certain threshold. Explicitly, this happens when τ < τ max , with τ max = − (2/min ), where min is the minimum eigenvalue of the squared operator L = −M −1 L∗ M −1 L. The selection of τ to optimize the convergence rate is made by minimizing the convergence factor R = max|1 + τ| = max{|1 + max τ|, |1 + min τ|}, 

(17)

for the nonzero eigenvalues  of L. The smallest R occurs when [36] τ = τ∗ = −

2 , max + min

(18)

for which R∗ =

min − max . min + max

(19)

The Modified Squared Operator Method MSOM is a variant of the SOM whose goal is to identify and eliminate the most harmful direction that slows down the convergence. Its description and analysis are also made in [36] and the formulas for the MSOM are  ϕ[ν+1] = ϕ[ν] − M −1 LM −1 L0 ϕ − α[ν] G[ν] , L∗ M −1 L0 ϕ G[ν] |ϕ=ϕ[ν] τ, (20) where α[ν] =

1 MG[ν] , G[ν]

−

1 τ LG[ν] , M −1 LG[ν]

,

and G[ν] = ϕ[ν] or G[ν] = ϕ[ν] − ϕ[ν−1] . Local convergence of (20) is ensured under the conditions for the convergence of the SOM, plus the assumptions LG[0] = / 0 and   2 1 MG[0] , LG[0] . , β= τ < τM = min − , min β − min MG[0] , G[0] The appendix in [36] provides a family of more general squared operator methods, by introducing powers of the acceleration operator M. Compared to (16), the implementation of (20) has the extra cost of computing the second term in the parenthesis. In all cases, and for the Fourier collocation approximation used in this paper, the inner products have been computed via the Parseval identity, the second choice of G[ν] has been used and M = c − ∂xx . The latter gave slightly better results (with suitable values of c) than other alternatives, as, for example, the linear part of L0 . 3.3. CG-Newton method The last technique described in this paper is the so-called Conjugate Gradient Newton method, or CGN method [33]. Given ϕ[ν] , Newton’s method generates the next iterate with the formulas ϕ[ν+1] = ϕ[ν] + ϕ[ν] , Lν ϕ

[ν]

= −L0 ϕ , [ν]

(21) (22)

where Lν is the linearization operator in (15) evaluated at ϕ[ν] . The classical Newton’s method uses direct methods to solve (22). Its implementation, for the computation of travelling-wave profiles, is usually affected by several difficulties, Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 8

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

mainly due to the singular Jacobian operator, the translational invariance of the equation under study, and the possibility of small denominators [2,11,12]. In particular, since (1) is invariant under spatial translations, the approximations to Lν will be singular or almost singular and the previously mentioned computational problems will be present in this case. This can be avoided by imposing the symmetry of the profile at x = 0 [11]. In the case of spectral methods, this symmetry condition can be imposed by using only cosine bases, instead of the complex exponential, and omitting the constant [10]. To this end, one of the Discrete Cosine Transforms (DCT) and its fast version may be used. In the case of (1), the implementation of the DCT is not trivial, since it must take into account the treatment, in the frequency domain, of the differential and nonlocal operators, along with the convolution [27,26,32]. Furthermore, an additional problem is present [9], since, due to the convolution and the nonlocal terms, the whole computation of Lν is required. An alternative to the classical version is given by the resolution of the discrete version of (22) with the Conjugate Gradient method [19], using the fact that L is self-adjoint. The iteration is usually preconditioned by a self-adjoint and positive definite acceleration operator M. Note that in the case of the Benjamin equation the technique must consist then of three processes: the homotopy procedure, the advance of the approximation (21), and the (preconditioned) Conjugate Gradient iteration to obtain the incremental term ϕ[ν] . Taking ϕ(0) = 0 for simplicity, the Preconditioned Conjugate Gradient method for the Newton correction (22) has, in its ith iteration, the following form [33]: R(0) = −L0 ϕ[ν] ; D(0) = M −1 R(0) , R(i) , M −1 R(i) , D(i) , Lν D(i) ϕ(i+1) = ϕ(i) + a(i) D(i) , a(i) =

(23)

R(i+1) = R(i) − a(i) Lν D(i) , R(i+1) , M −1 R(i+1) , R(i) , M −1 R(i) = M −1 R(i+1) + b(i+1) D(i) .

b(i+1) = D(i+1)

The standard L2 inner product is used. It is worth recalling some remarks from [33]. Normally, the Conjugate Gradient method may be applied to self-adjoint, positive definite operators. When this is not satisfied, an alternative version of (23) can be used [33,34]. On the other hand, the zero eigenvalue of L, which makes the solution of (22) not unique and generates the previously mentioned computational problems in Newton’s method, does not seem to cause a trouble here. The symmetry of the profile is ensured in a simple way [33], by taking an even initial guess ϕ[0] . As mentioned before, a resolution with direct methods also has the same problem (the discretization of L is usually a singular, or nearly singular, matrix) and it can be corrected similarly, but the imposition of this assumption is computationally harder [11]. For the corresponding Fourier collocation version of the method used in the next section, all computations are performed in the Fourier space and have, in the ith inner iteration, the form 

 (0) R



p

= (−cs

+ α − γ|p| + δp2 )

 ϕh[ν]

 + p

β  [ν] [ν] ϕ . ∗ ϕh , p 2 h

 (0) ) (R p  (0) ) = , (D p c + p2        (i) (i) LD = (−cs + α − γ|p| + δp2 ) D + β ϕh[ν] D(i) , p







p

p

 (i) |2 /(c + p2 ) −N/2≤p≤N/2 | (R p (i) a =    ,    (i) (i) LD D −N/2≤p≤N/2 p p       (i+1) (i) D(i) p , = ϕ + a(i)  ϕ p

p

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx



 (i+1) (R

p

9

    (i) (i) = (R − a(i) LD , p







p

 (i+1) |2 /(c + p2 ) (R p  =  ,   (i) |2 /(c + p2 ) −N/2≤p≤N/2 | (R p   (i+1) R p   (i+1) ) = (i) ) . + b(i+1) (D (D p p c + p2

b(i+1)

−N/2≤p≤N/2 |

4. Comparison of methods and generation of solitary waves In this and the following section we present a series of numerical experiments in order to compare the performance of the methods considered above when used to generate solitary-wave profiles for the Benjamin equation (1) and (9). The numerical experiments can be divided into two groups. The first, treated in this section, makes general comparisons between the methods by measuring the errors as functions of the number of iterations and of the cpu time. The second group is presented in Section 5 and uses the results of the most efficient method to emerge from the first group of experiments. The corresponding profiles are taken as initial conditions of a numerical evolution method integrating (1) in time. The evolution of the approximate solitary waves is monitored up to T = 300. The goal is to study the accuracy of the computed profiles as travelling-wave solutions of the Benjamin equation. We will also discuss the computation of multi-pulse solitary waves. 4.1. Introductory remarks and parameter selection for the methods The values α = 1, β = 1, δ = 1, cs = 0.75 in (1) have been fixed, so that the continuation parameter is γ. As shown in [2], waves of physical relevance correspond to γ between 0 and 1. (Note that with the previous values, γ coincides with the γ˜ in (8).) The techniques implemented will be: the Petviashvili method [PM], the SOM ([SOM]) and MSOM ([MSOM]) methods, and the Conjugate-Gradient Newton method ([CGN]). The homotopy algorithm has been performed with a stepsize γ = 10−3 . The solitary-wave profile is then computed at values γ j = γ 0 + jγ, j = 1, . . ., M, with γ 0 = 0, for which the solitary wave is known (KdV). The performance of each method will be analysed at the final value of the γ considered and at a general step of the numerical continuation, from γ j to γ j+1 . The error at the νth step of the iteration has been measured by E[ν] =

||ϕh[ν] − ϕh[ν−1] || ||ϕh[ν] ||

,

in the Euclidean norm. For the [CGN], the norm of the residual error EL[ν] = ||L0 ϕh[ν] ||,

(24)

has also been used. Three parameters control the procedure: The iteration stops when a maximum number of iterations is exceeded, or when the errors E[ν] or EL[ν] are smaller than a given tolerance TOL. This general strategy has some variants, depending on the specific implementation of the method, especially in the value of TOL and the way of counting the number of iterations. The study of comparison of the methods requires first some preliminary remarks concerning particular aspects of each of them. For the [PM], the local convergence can be observed by studying the stabilizing factor mν in one continuation step γ j → γ j+1 , Fig. 1 shows, for the values of γ indicated, the behaviour of the error in the stabilizing factor mν as a function of the number of iterations. When γ > 0.1, similar results are observed. It is also noted that the stabilizing factor approaches 1 to machine accuracy quite fast. For instance, for γ = 0.034, the discrepancy |mν − 1| is O(10−13 ) in NITER = 15 iterations and O(10−16 ) when NITER = 25. For the [SOM] and [MSOM], the corresponding discrete version of M = c − ∂xx has been taken as the preconditioning operator. The behaviour of these two methods strongly depends on the parameters c, τ and γ. The main conclusion Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 10

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx −7

γ=0.001

−7

(a) 0.5

x 10

γ=0.1

x 10

(b)

0.5

0

0

−0.5

mν−1

m −1 ν

−0.5

−1

−1

−1.5 −1.5 −2 −2

−2.5 0

−2.5

10

20

30

40

50

60

−3

0

5

10

NITER

15

20

25

30

35

40

45

50

NITER

Fig. 1. Discrepancy mν − 1 of the stabilizing factor of [PM] vs. number of iterations for the continuation steps: (a) 0 → 0.001; (b) 0.099 → 0.1.

is that as γ ↑ 1, the rate of convergence of the methods goes to one very fast. This slows down both methods and sometimes prevents the SOM and MSOM to be competitive to other alternatives, specially for γ close to 1. The techniques described in [36] to estimate the rate of convergence of the SOM can be applied in the case at hand. They are as follows: For some values of c, the eigenvalues of the discretization of the iteration operator L are computed. Then the convergence factor R* (c) is calculated from (19). The dependence of R* on τ is obtained by estimating (17) for some values of the step size. The optimal opt τ practically coincides with that of formula (18). For the solitary-wave profile it is necessary to compute the iteration matrix (see the operator L in (15)). In the following experiments, for each value of γ, we compute the solitary wave given by the [PM] with N = 512 points on the interval [− 64, 64]. Fig. 2 shows, for γ = 0.1, 0.5, 0.9, 0.99, the convergence rate R* (c) as a function of c. Its minimum value, and the corresponding value of c = copt are shown in Table 1. Observe that, as γ tends to 1, the minimum of R* also grows, becoming, for the last γ, very close to one. Fig. 3 shows the rate R, evaluated at the optimal c = copt , as a function of τ. The value of the step size which gives the minimum, τ opt , coincides with the value given by (18) to four digits, since τ = 10−3 . These experiments do not seem to recommend choosing [SOM], at least for γ close to one. Observe that the dependence on γ is also somehow incompatible with the continuation algorithm. For instance, the value of τ opt corresponding to γ = 0.9 gives values of R that are far from the corresponding minimum R* associated with previous values of γ. Hence, it is hard to establish a common τ opt for the whole continuation algorithm. An alternative strategy would be to compute the best value of c and τ for each γ, but this is very costly. The first points to emphasize in the case of [CGN] are the stopping criteria and the number of iterations. Recall that the method contains, in a step γ j → γ j+1 , two loops, corresponding to the advance of the Newton correction and to the application of the CG method. For the first one, the stopping criterion is controlled by a fixed maximum number of iterations and by a small tolerance (TOL = 10−15 in our case) for the residue EL[ν] . For the CG iteration, the technique Table 1 [SOM]: Computed values of copt , R* (copt ) and τ opt . γ

copt

R* (copt )

τ opt

0.1 0.5 0.9 0.95 0.99

0.248 0.259 0.270 0.270 0.270

0.203 0.602 0.979 0.995 0.999

1.0641 1.653 2.274 2.310 2.326

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx γ=0.1

γ=0.5 1

0.8

0.9

R*

R*

1

0.6

0.8

0.4

0.2

11

0.7

0

0.5

1

0.6

1.5

0

0.5

c

1

1.5

1

1.5

c γ=0.99

γ=0.9 1.005

1

1

0.9999

0.995 *

R

R

*

0.9999 0.99

0.9998 0.985 0.9998

0.98 0.975

0

0.5

1

1.5

0.9997

0

0.5

c

c

Fig. 2. [SOM]: Convergence rate R vs. c following formula (19). γ=0.1

γ=0.5

1.6

1

1.4 0.9

1.2 R

R

1 0.8

0.8 0.6

0.7

0.4 0.2

0

0.5

1

1.5

2

0.6

2.5

0

0.5

1

Δτ

1.5

2

2.5

Δτ

γ=0.9

γ=0.99 1.002

1 1.001

0.98

R

R

0.99

1

0.97 0.999 0.96 0.95

0

1

2 Δτ

3

0.998

0

1

2

3

Δτ

Fig. 3. [SOM]: Convergence rate R = R(τ, copt ) vs. τ following formula (17).

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 12

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx 4

(a)

7

x 10

(b) γ=0.1 γ=0.5 γ=0.9 γ=0.95 γ=0.99

6

12

10

8 4

CPU (seconds)

Number of iterations

5

3

6

4

2 2

1

0 0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

0 0.2

0.4

0.6

0.8

c

1

1.2

1.4

1.6

1.8

2

c

Fig. 4. [CGN]: Number of iterations (a) and cpu time (b) vs. c.

described in [33] has been used. This adapts the process to the accuracy of the outer iteration. The error is measured by the norm of the function R(i) in the CG formulas, so that if ||R(i) ||M < CG ||R(0) ||M , the iteration stops. Here || · ||M is the M−1 -weighted L2 norm ||ϕ||M = ϕ, M −1 ϕ 1/2 ,

M = c − ∂xx ,

where · , · denotes the L2 inner product in (− l, l), R(0) = − L0 ϕ[ν] is the residue of the method and CG is an error tolerance parameter, which relates the two iterations. In order to optimize the two procedures, it is recommended in [33] to take CG between 10−1 and 10−3 . We found that CG = 10−2 gives good results in our case. Finally, the comparisons with other methods, shown in Figs. 9 and 10, are made by computing the error with the usual L2 norm of the residue (24). The implementation of two iterative procedures implies the existence of two iterations: one for the Newton correction and one for the CG method. The Newton correction requires one or two iterations, so we have chosen to present data for the total number of inner CG iterations. This number estimates the global computational cost of the method. On the other hand, the [CGN] method also requires the selection of the parameter c of the preconditioning operator. The data of Fig. 4 may be used as an aid to this. The figure shows, for five values of γ, the number of iterations (left) and the cpu time (right) vs. the parameter c. Fig. 4 suggests that the best value of c is likely to be in the interval [0.25, 0.3]. A further refinement suggests that c = 0.275 is a good value; this was used in the numerical experiments in the sequel. 4.2. A comparative study of generation of solitary waves 4.2.1. Solitary-wave profiles Fig. 5 shows some solitary-wave profiles of the Benjamin equation obtained with the [PM] by using a pseudospectral code with N = 2048 points on the interval [− 256, 256], from the known KdV profile (γ = 0) to the one corresponding to γ = 0.99. The generation of the symmetric lobules (whose number increases with γ) is observed (cf. [2]). It is worth mentioning that the other methods gave entirely similar profiles. Some differences in their performance will be pointed out below. 4.2.2. Continuation procedure As noted before, the continuation method was implemented with a stepsize γ = 10−3 , while the tolerance for the procedures was taken as TOL = 10−15 for all the methods except for [SOM], where we took TOL = 5 ×10−15 . Table 2 Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx γ=0

γ=0.1

0

0

−0.2

−0.2

−0.4

−0.4

−0.6

−0.6

−0.8 −300 −200 −100

0 γ=0.5

13

100

200

300

−0.8 −300 −200 −100

0 γ=0.9

100

200

300

0 100 γ=0.99

200

300

200

300

0.1 0 0 −0.2

−0.1

−0.4

−0.2

−300 −200 −100

0 100 γ=0.95

200

300

−300 −200 −100

0.1

0.1

0

0

−0.1

−0.1

−0.2 −300 −200 −100

0

100

200

300

−0.2 −300 −200 −100

0

100

Fig. 5. Solitary-wave profiles for the [PM] on the interval [− 256, 256] and with N = 2048.

shows, for [PM], the total number of iterations and the cpu time required to get to the corresponding value of γ, when the errors in each iterative step, E[ν] or EL[ν] , are less than the tolerance. The last column shows the average number of iterations required in an application of the [PM]; that is, to go from γ j to the next one γ j+1 = γ j + γ. These values are computed using information from the interval of the two most recent values of γ. (For example, the quantity 43.5 for γ = 0.5 is obtained by calculating the mean value of iterations between γ = 0.1 and γ = 0.5; the quantity 32.4 is obtained using the iterations between γ = 0.5 and γ = 0.9, and so on.) This may help to identify the range of γ for which the iterative process is slower. In this sense, it is seen that for the last values of γ (when the oscillations grow in number), the computations seem to be harder. While up to γ = 0.9, the [PM] computes a step γ j → γ j+1 in not more than 50 iterations, this number starts to grow for γ greater than 0.95 and, for instance, the method requires the same cpu time and number of iterations to go from γ = 0 to γ = 0.1 as those to go from γ = 0.95 to γ = 0.99. The global mean value of iterations, up to γ = 0.99, is about 42.95. Table 3 provides the same information but for [SOM]. (Here the computation for γ = 0.9 attained the maximum number of iterations allowed.) Clearly, the convergence is slow. The situation is slightly different in the case of [MSOM]. (Here we used the same values of c and τ as in [SOM].) The improvement, compared to [SOM], can be seen in Table 4. (At γ = 0.99, the maximum number of iterations is attained and the results are not shown.) We observe that, for a smaller tolerance (TOL = 10−15 now), the [MSOM] is able to give results for γ greater than 0.9 in a computationally reasonable time improving those of the [SOM]. Table 2 Number of iterations and cpu time for solitary-wave generation with the [PM]. γ

NITER

cpu (s)

NITER (mean value)

0.1 0.5 0.9 0.95 0.99

5040 22,436 35,374 38,074 42,517

2.03 9.68 15.25 16.44 18.47

50.4 43.5 32.4 54.0 111.1

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 14

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

Table 3 Number of iterations and cpu time for solitary-wave generation with the [SOM]. The computations corresponding to γ = 0.9 attained the maximum number (105 ) of iterations. γ

NITER

cpu (s)

NITER (mean value)

0.1 0.5

8410 51,457

7.64 48.78

84.1 108.3

Table 4 Number of iterations and cpu time for solitary-wave generation with the [MSOM] and TOL = 10−15 . The computations corresponding to γ = 0.99 attained the maximum (105 ) of iterations. γ

NITER

cpu (s)

NITER (mean value)

0.1 0.5 0.9 0.95

3743 21,546 58,075 79,328

3.90 22.54 62.34 85.84

37.4 44.5 91.3 531.3

Table 5 Number of iterations and cpu time for solitary-wave generation with the [CGN]. TOL = 10−15 . γ

NITER

cpu (s)

NITER (mean value)

0.1 0.5 0.9 0.95 0.99

1577 9351 20,842 23,515 27,185

0.29 2.26 4.29 4.76 5.32

15.8 19.4 28.7 66.8 91.8

In the case of the [CGN], with CG = 10−2 , TOL = 10−15 , c = 0.275, the corresponding information about the continuation process appears in Table 5. The best performance of the [CGN] is evident, specially for values of γ closer to one. To verify this we display in Fig. 6 the number of iterations and the cpu time as functions of γ, for the three more competitive methods. The Petviashvili method [PM] continues to be more efficient in the global continuation procedure than [MSOM], both in number of iterations and computational time, although these are comparable for small values of γ. 4

(a) 10

x 10

9

(b) 120

PM MSOM CGN

PM MSOM CGN

100

8

CPU time(in seconds)

Number of iterations

7 6 5 4 3

80

60

40

2 20 1 0 0.1

0.2

0.3

0.4

0.5

γ

0.6

0.7

0.8

0.9

1

0 0.1

0.2

0.3

0.4

0.5

γ

0.6

0.7

0.8

0.9

Fig. 6. Number of iterations (a) and cpu time (b) vs. γ, for [PM], [MSOM] and [CGN].

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

1

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx γ=0.001

γ=0.035

−2

−2 PM SOM MSOM

−6

−4 log10(ERROR)

log10(ERROR)

−4

−8 −10 −12 −14 −16

−6 −8 −10 −12 −14

0

20

40 60 NITER

80

−16

100

0

20

40 60 NITER

80

100

80

100

γ=0.1 −2

−4

−4

−6

−6

log10(ERROR)

log10(ERROR)

γ=0.070 −2

−8 −10 −12 −14 −16

15

−8 −10 −12 −14

0

20

40 60 NITER

80

−16

100

0

20

40 60 NITER

Fig. 7. Logarithm of the error E[ν] vs. number of iterations for the continuation steps: 0 → 0.001, 0.034 → 0.035, 0.069 → 0.070, 0.099 → 0.1.

4.2.3. Iteration errors The four methods can also be compared in terms of the iteration errors. Fig. 7 shows, for small values of γ, the logarithm of the error E[ν] against the number of iterations needed by the [SOM], the [PM] and the [MSOM]. In all cases, the [MSOM] appears as the most efficient. However, this best performance of the [MSOM] does not continue as γ grows. This is observed in Fig. 8. Here, [MSOM] and [PM] are compared with values of γ closer to one and [PM] behaves better. The comparison of the errors associated with the residue reveals that [CGN] has the advantage. Fig. 9 shows, for small values of γ, the L2 norm of the error of the corresponding residue as a function of the number of iterations. This figure shows that the [CGN] attains, for a fixed number of iterations, the smallest error, and needs, for a fixed error γ=0.5

(a) −2

γ=0.535

(b) −2

−4

−4

−6

−6 log10(ERROR)

log10(ERROR)

PM MSOM

−8

−10

−8

−10

−12

−12

−14

−14

−16 0

10

20

30 NITER

40

50

60

−16 0

5

10

15

20

25 30 NITER

35

40

45

50

Fig. 8. Logarithm of the error E[ν] vs. number of iterations for the continuation steps: (a) 0.499 → 0.5; (b) 0.534 → 0.535.

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 16

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx γ=0.001

γ=0.035 0 PM SOM MSOM CGN

−5

−2 log10(ERROR)

log10(ERROR)

0

−10

−4 −6 −8 −10 −12

−15

0

20

40 60 NITER

80

−14

100

0

20

40 60 NITER

0

−2

−2

−4

−4

−6 −8 −10 −12 −14

100

80

100

γ=0.1

0

log10(ERROR)

log10(ERROR)

γ=0.070

80

−6 −8 −10 −12

0

20

40 60 NITER

80

−14

100

0

20

40 60 NITER

Fig. 9. Logarithm of the error EL[ν] vs. number of iterations for the continuation steps: 0 → 0.001, 0.034 → 0.035, 0.069 → 0.070, 0.099 → 0.1.

level, the minimum number of iterations. This persists during the whole continuation procedure, as it can be seen, for example, for some values of γ close to 0.5 and the three methods [PM], [MSOM] and [CGN], in Fig. 10. 5. Generation of multi-pulse solitary waves and time-stepping procedure According to the experiments presented above, it is clear that [CGN] emerges as the most efficient method from the single solitary-wave computations. This method is used to perform the second group of experiments, namely the γ=0.5

(a) −2

γ=0.535

(b) −2 PM MSOM CGN

−4

−4

−6 log10(ERROR)

log10(ERROR)

−6

−8

−8

−10

−10 −12 −12

−14 0

−14

10

20

30 NITER

40

50

60

−16 0

5

10

15

20

25 30 NITER

35

40

45

50

Fig. 10. Logarithm of the error EL[ν] vs. number of iterations for the continuation steps: (a) 0.499 → 0.5; (b) 0.534 → 0.535.

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx γ=0

17

γ=0.5 0.1

0

0 −0.1

−0.2

−0.2 −0.4

−0.3 −0.4

−0.6 −0.5 −0.8 −300 −200 −100

0

100

200

300

−300 −200 −100

γ=0.9 0.2

0

0

−0.2

−0.2

0

100

200

300

100

200

300

γ=0.95

0.2

−0.4 −300 −200 −100

0

100

200

300

−0.4 −300 −200 −100

0

Fig. 11. Two-pulse solitary-wave profiles for the [CGN] on the interval [− 256, 256] and with N = 2048.

generation of multi-pulse solitary waves and the computation of solitary-wave profiles that are posed initial conditions of a temporal integration of the Benjamin equation. 5.1. Multi-pulse solitary waves We started the continuation procedure at γ = 0 using as initial profile a linear superposition of KdV solitary waves of the same amplitude (corresponding to cs = 0.75) and whose centers are located at x = −10, 10 (in the case of a two-pulse) and at x = −10, 0, 10 (in the case of a three-pulse). We used N = 2048 points on the interval [− 256, 256]. Fig. 11 shows the profiles of two-pulse solitary waves that emerged from the [CGN] procedure for different values of γ, while three-pulses are captured in Fig. 12. In these experiments, the performance of the numerical method deteriorated. For instance, we have observed that, in the case of two-pulses, and γ = 0.99 [CGN] was not convergent. Fig. 13 displays a magnification of the three-pulses generated when γ = 0.9, 0.95. 5.2. Temporal evolution of numerical solitary waves A necessary test to evaluate the accuracy of the computed solitary-wave profiles is checking that they are travellingwave solutions of the time-dependent partial differential equation. To this end, we took the solitary-wave profiles generated by [CGN] as initial conditions and integrated the Benjamin equation in time, monitoring the amplitude and speed to assess the accuracy of the travelling-wave numerical solution. The numerical method here uses a Fourier pseudospectral discretization in space, while the temporal integration is carried out with the third-order SDIRK method (see e.g. [21] and references therein), corresponding to the tableau

(25) Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 18

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx γ=0

γ=0.5 0.1

0

0 −0.1

−0.2

−0.2 −0.3

−0.4

−0.4 −0.5

−0.6

−0.6 −0.8 −300 −200 −100

0

100

200

300

−300 −200 −100

γ=0.9

0

100

200

300

100

200

300

γ=0.95

0.1

0.1

0

0

−0.1

−0.1

−0.2

−0.2

−0.3

−0.3

−0.4

−0.4

−0.5 −300 −200 −100

0

100

200

−0.5 −300 −200 −100

300

0

Fig. 12. Three-pulse solitary-wave profiles for the [CGN] on the interval [− 256, 256] and with N = 2048.

and modified with a projection technique [20], to preserve the discrete versions of the invariants C and F, in (2) and (3) respectively, defined as Ch (Z) = h

N−1 

h 2 Zj , 2 N−1

Zj ,

Ih (Z) =

j=0

j=0

(a)

(b)

0.1

0.1

0

0

−0.1

−0.1

−0.2

−0.2

−0.3

−0.3

−0.4

−0.4

−0.5 −60

−40

−20

0

20

40

60

−0.5 −60

−40

−20

0

20

40

60

Fig. 13. Three-pulse solitary-wave profiles for the [CGN] on the interval [− 256, 256] and with N = 2048. Magnification of Fig. 12 for: (a) γ = 0.9; (b) γ = 0.95.

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx (a)

(b)

0.1

(c)

0.1

0

0.1

0

0

−0.1

−0.1

−0.2

−0.2

−0.2

u

u

u

−0.1

−0.3

−0.3

−0.3

−0.4

−0.4

−0.4

−0.5

−0.5

−0.5

−300

−200

−100

0

100

200

300

19

−300

−200

−100

x

0

100

200

300

−300

−200

−100

0

x

100

200

300

x

Fig. 14. Numerical solution from the initial s-w profile for γ = 0.5: (a) t = 0; (b) t = 200; (c) t = 300.

0

0.05 0

−0.05

−0.05

−0.05

u

u

0.05

0

u

0.05

−0.1

−0.1

−0.15

−0.1

−0.15

−0.15

−0.2

−0.2

−0.2

−0.25

−0.25

−0.25

−300

−200

−100

0

x

100

200

300

−300

−200

−100

0

x

100

200

300

−300

−200

−100

0

x

100

200

300

Fig. 15. Numerical solution from the initial s-w profile for γ = 0.9: (a) t = 0; (b) t = 200; (c) t = 300.

for Z = (Z0 , . . ., ZN−1 )T . We observe that preservation of the discrete version of the Hamiltonian (4), Eh = αIh + Gh with ⎛

N−1 

Gh (Z) = h ⎝

j=0

⎞ δ γ β (DZ)2j − (Zj (Hh Z)j ) + (Z)3j ⎠ , 2 2 6

cannot be simultaneously achieved. The reason is that at the exact solitary wave, the quantities δE and δI are dependent, as (5) shows. This also happens in the discrete case. Thus, the numerical projection technique to preserve the discrete versions of the functionals, applied to E and I at the same time [20,21], does not work. The step sizes have been taken as h = 0.25 for the spatial pseudospectral discretization (used for the computation of the initial profiles as well) and t = 6.25 × 10−3 for the temporal discretization. (Smaller values of t did not change the results significantly.) 0.1

0

0

0

u

−0.1

−0.2 −300

u

0.1

u

0.1

−0.1

−200

−100

0

x

100

200

300

−0.2 −300

−0.1

−200

−100

0

100

200

300

−0.2 −300

−200

−100

x

0

100

200

300

x

Fig. 16. Numerical solution from the initial s-w profile for γ = 0.99: (a) t = 0; (b) t = 200; (c) t = 300.

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 20

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx x 10

x 10

0.1 0

2

2

1

1

0

0

−0.1

u

u

u

−0.2 −0.3 −1

−1

−2

−2

−0.4 −0.5

−300

−200

−100

0

100

200

−300

300

−200

−100

0

100

200

300

−300

−200

−100

x

x

0

100

200

300

100

200

300

x

Fig. 17. Magnification of Fig. 14. x 10

0.05

x 10

2

2

1

1

0

0

−0.05

u

u

−0.1

u

0

−0.15 −1

−1

−2

−2

−0.2 −0.25

−300

−200

−100

0

100

200

−300

300

−200

−100

0

100

200

300

−300

−200

−100

x

x

0

x

Fig. 18. Magnification of Fig. 15.

Figs. 14–16 show the evolution of the initial profile, corresponding to γ = 0.5, 0.9, 0.99 respectively, up to time t = 300. More information is provided in Figs. 17–19. They correspond to the previous figures, but here the waves at times t = 200, 300 have been magnified. In all cases, no spurious forward or backward oscillations are observed and the corresponding profiles propagate without any disturbances. In order to assess the accuracy of the numerical travelling-wave solution we study its amplitude and speed. In Table 6, we show the minimum value of the initial profile, obtained by computing the minimum of the corresponding interpolating trigonometric polynomial. Table 6 has been used to generate Fig. 20. Here, for each value of γ (including the KdV case γ = 0), the initial amplitude has been compared with the corresponding amplitude (computed in the same way) of the numerical solution up to t = 300. For all γ, the amplitude error is of order of 10−9 and does not grow with time. A similar study can be made for the speed. The initial profiles were all computed with a speed cs = 0.75, which was the speed of reference to generate the errors observed in Fig. 21. The speed of the travelling wave has been computed x 10

x 10

0.1

−0.1

−0.2 −300

−200

−100

0

x

100

200

300

2

1

1

0

0

u

u

u

0

2

−1

−1

−2

−2

−300

−200

−100

0

100

200

x

300

−300

−200

−100

0

100

200

300

x

Fig. 19. Magnification of Fig. 16.

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx −10

9

−9

x 10

1.2

γ=0 (KdV)

x 10

γ=0.5

1 Amplitude error

Amplitude error

8 7 6 5 4 3

0.8 0.6 0.4 0.2

0

100

200

0

300

0

100

Time

200

300

Time

−10

5

21

−10

x 10

2

x 10

γ=0.99

γ=0.95 Amplitude error

Amplitude error

4 3 2

1

1 0

0

100

200

0

300

0

100

Time

200

300

Time

Fig. 20. Temporal evolution of the error between the amplitude of the initial solitary-wave profile and of the evolving solution. −9

7

−8

x 10

γ=0 (KdV)

6

4 3 2

γ=0.5

0.6 0.4 0.2

1 0

x 10

0.8

5

Error in speed

Error in speed

1

0

100

200

0

300

0

100

Time −9

3.5

300

−9

x 10

2

γ=0.95

3 2.5

Error in speed

Error in speed

200 Time

2 1.5 1

x 10

γ=0.99

1.5

1

0.5

0.5 0

0

100

200 Time

300

0

0

100

200

300

Time

Fig. 21. Temporal evolution of the error between the speed of the initial solitary-wave profile and of the evolving solution.

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 22

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

Table 6 Computed minimum of the initial solitary-wave profiles, for some values of γ. γ

Amplitude

0 0.1 0.5 0.9 0.95 0.99

−0.75 −7.183404555524147e−001 −5.411747962142054e−001 −2.280941079089593e−001 −1.656679129094486e−001 −9.035736815781438e−002

−15

−15

x 10

6

γ=0 (KdV)

4

Energy error

Energy error

6

2

0

0

100

200

γ=0.5

4

2

0

300

x 10

0

100

200

Time −16

−16

x 10

1.2

γ=0.95

5

1

4

0.8

Energy error

Energy error

6

3 2 1 0

300

Time x 10

γ=0.99

0.6 0.4 0.2

0

100

200

0

300

0

100

200

Time

300

Time

Fig. 22. Temporal evolution of the error between Eh of the initial solitary-wave profile and of the numerical solution.

0.2

0.2

0

0

u

u

0

−0.2

−0.4 −300

u

0.2

−0.2

−200

−100

0

x

100

200

300

−0.4 −300

−0.2

−200

−100

0

x

100

200

300

−0.4 −300

−200

−100

0

100

200

300

x

Fig. 23. Numerical evolution from the initial two-pulse solitary-wave profile for γ = 0.9: (a) t = 0; (b) t = 200; (c) t = 300.

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

23

in the standard way [16], by dividing the increment of the spatial translation of the center of the pulse over t at each time step. The errors in the speed are also quite small and remain bounded. This and the previous figure suggest that the computed initial solitary-wave profiles are good approximations to the real ones and that the temporal discretization is accurate. We have previously mentioned that the time integrator is not initially designed to preserve the discrete version of the Hamiltonian. In Fig. 22, the evolution of the error between the values of Eh of the numerical solution and of the initial profile up to t = 300 is displayed. The method virtually preserves the discrete Hamiltonian, which is also a suggestion of good approximation, since the integrator is only constructed to preserve the quadratic quantity Ih [17]. The accuracy of the computed profiles is also reflected in the simulation of the multi-pulse waves. As an illustration, Fig. 23 shows the evolution of a two-pulse solitary-wave profile, for γ = 0.9, at different times. Acknowledgement This research has been supported by MICINN project MTM2010-19510/MTM. References [1] M.J. Ablowitz, Z.H. Musslimani, Spectral renormalization method for computing self-localized solutions to nonlinear systems, Optics Letters 30 (2005) 2140–2142. [2] J.P. Albert, J.L. Bona, J.M. Restrepo, Solitary-wave solutions of the Benjamin equation, SIAM Journal of Applied Mathematics 59 (1999) 2139–2161. [3] D.M. Ambrose, J. Wilkening, Computation of time-periodic solutions of the Benjamin-Ono equation, Journal of Nonlinear Science 20 (2010) 277–308. [4] W. Bao, Q. Du, Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow, SIAM Journal on Scientific Computing 25 (5) (2004) 1674–1697. [5] T.B. Benjamin, Internal waves of permanent form in fluids of great depth, Journal of Fluid Mechanics 29 (1967) 559–592. [6] T.B. Benjamin, A new kind of solitary wave, Journal of Fluid Mechanics 245 (1992) 401–411. [7] T.B. Benjamin, Solitary and periodic waves of a new kind, Philosophical Transactions of the Royal Society of London – Series A: Mathematical, Physical and Engineering Sciences 354 (1996) 1775–1806. [8] J.L. Bona, M. Chen, A Boussinesq system for two-way propagation of nonlinear dispersive waves, Physica D 116 (1998) 191–224. [9] J.P. Boyd, Solitons from sine waves: analytical and numerical methods for non-integrable solitary and cnoidal waves, Physica D 21 (1986) 227–246. [10] J.P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed., Dover Publications, New York, 2000. [11] J.P. Boyd, Deleted residuals, the QR-factored Newton iteration, and other methods for formally overdetermined determinate discretizations of nonlinear eigenproblems for solitary, cnoidal, and shock waves, Journal of Computational Physics 179 (2002) 216–237. [12] J.P. Boyd, Why Newton’s method is hard for travelling waves: small denominators, KAM theory, Arnold’s linear Fourier problem, nonuniqueness, constraints and erratic failure, Mathematics and Computers in Simulation 74 (2007) 72–81. [13] J.P. Boyd, Z. Xu, Comparison of three spectral methods for the Benjamin-Ono equation: Fourier pseudospectral, rational Christov functions and Gaussian radial basis functions, Wave Motion 48 (2011) 702–706. [14] J.P. Boyd, Z. Xu, Numerical and perturbative computations of solitary waves of the Benjamin-Ono equation and higher order nonlinearity using Christov rational basis functions, Journal of Computational Physics 231 (2012) 1216–1229. [15] M. Caliari, A. Ostermann, S. Rainer, M. Thalhammer, A minimisation approach for computing the ground state of Gross-Pitaevskii systems, Journal of Computational Physics 228 (2000) 349–360. [16] V.A. Dougalis, A. Durán, M.A. López-Marcos, D. Mitsotakis, A numerical study of the stability of solitary waves of the Bona-Smith family of Boussinesq systems, Journal of Nonlinear Science 17 (2007) 569–607. [17] J. de Frutos, J.M. Sanz-Serna, Accuracy and conservation properties in numerical integration: the case of the Korteweg–de Vries equation, Numerische Mathematik 75 (1997) 421–445. [18] J.J. Garcia-Ripoll, V.M. Perez-Garcia, Optimizing Schrödinger functionals using Sobolev gradients: applications to quantum mechanics and nonlinear optics, SIAM Journal on Scientific Computing 23 (2001) 1316–1334. [19] G. Golub, C. Van Loan, Matrix Computations, 3rd ed., The Johns Hopkins University Press, Baltimore, 1996. [20] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration, Structure-Preserving Algorithms for Ordinary Differential Equations, Springer-Verlag, New York/Heidelberg/Berlin, 2002. [21] E. Hairer, S.P. Norsett, G. Wanner, Solving Ordinary Differential Equations I. Nonstiff Problems, 2nd ed., Springer Series in Computational Mathematics, vol. 8, Springer-Verlag, New York/Heidelberg/Berlin, 1993. [22] H. Kalisch, J.L. Bona, Models for internal waves in deep water, Discrete and Continuous Dynamical Systems 6 (2000) 1–20. [23] T.I. Lakoba, J. Yang, A generalized Petviashvili method for scalar and vector Hamiltonian equations with arbitrary form of nonlinearity, Journal of Computational Physics 226 (2007) 1668–1692. [24] T.I. Lakoba, J. Yang, A mode elimination technique to improve convergence of iteration methods for finding solitary waves, Journal of Computational Physics 226 (2007) 1693–1709.

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008

+Model MATCOM-3828; No. of Pages 24 24

ARTICLE IN PRESS

V.A. Dougalis et al. / Mathematics and Computers in Simulation xxx (2012) xxx–xxx

[25] F. Linares, M. Scialom, On generalized Benjamin type equations, Discrete and Continuous Dynamical Systems 12 (2005) 161–174. [26] S.A. Martucci, Symmetric convolution and the discrete sine and cosine transforms, IEEE Transactions on Signal Processing 42 (1994) 1038–1051. [27] H. Olkkonen, P. Pesola, J.T. Olkkonen, Computation of Hilbert transform via discrete cosine transform, Journal of Signal and Information Processing 1 (2010) 18–23. [28] D.E. Pelinovsky, Y.A. Stepanyants, Convergence of Petviashvili’s iteration method for numerical approximation of stationary solutions of nonlinear wave equations, SIAM Journal on Numerical Analysis 42 (2004) 1110–1127. [29] B. Pelloni, V.A. Dougalis, Numerical solution of some nonlocal, nonlinear, dispersive wave equations, Journal of Nonlinear Science 10 (2000) 1–22. [30] V.I. Petviashvili, Equation of an extraordinary soliton, Soviet Journal of Plasma Physics 2 (1976) 257–258. [31] W.C. Rheinboldt, Numerical continuation methods: a perspective, Journal of Computational and Applied Mathematics 124 (2000) 229–244. [32] G. Strang, The Discrete Cosine Transform, SIAM Review 41 (1999) 135–147. [33] J. Yang, Newton-conjugate-gradient methods for solitary wave computations, Journal of Computational Physics 228 (2009) 7007–7024. [34] J. Yang, Nonlinear Waves in Integrable and Nonintegrable Systems, SIAM, Philadelphia, 2010. [35] J. Yang, T.I. Lakoba, Accelerated imaginary-time evolution methods for the computation of solitary waves, Studies in Applied Mathematics 120 (2008) 265–292. [36] J. Yang, T.I. Lakoba, Universally-convergent squared-operator iteration methods for solitary waves in general nonlinear wave equations, Studies in Applied Mathematics 118 (2007) 153–197. [37] J. Yang, B.A. Malomed, D.J. Kaup, Embedded solitons in second-harmonic-generating systems, Physical Review Letters 83 (1999) 1958–1961.

Please cite this article in press as: V.A. Dougalis, et al. Numerical approximation of solitary waves of the Benjamin equation, Math. Comput. Simul. (2012), http://dx.doi.org/10.1016/j.matcom.2012.07.008