Soliton solutions to the fifth-order Korteweg - de Vries equation and their applications to surface and internal water waves
K. R. Khusnutdinova Y. A. Stepanyants M. R. Tranter
physics.flu-dyn nlin.PS







 

 

 

headings









Figures/


thefnmarkfootnotetext







C[1]>m#1


















||


[1]



equation1


    ε∂[1]^#1 ( ){}[1]  #1   [ ]1/2  ⇒  [2] cosh ^#1 #2 [2] sech ^#1 #2 [2] tanh ^#1 #2 arccosh[1] arccosh #1 [1]#1[1]#1[1]#1_-∞^∞[1] d #1[2]d #1/d #2[3]d^#1 #2d #3^#1[1] ∂ #1[2]∂ #1/∂ #2[3]∂^#1 #2∂ #3^#1*theoremTheorem[section]
        definitionDefinition[section]
        exaExample[section]
        lemLemma[section]
        corCorollary[section]
        proProposition[section]
        
        
        ^1) Department of Mathematical
        Sciences,
        Loughborough University, 
        Loughborough LE11 3TU, United Kingdom; 
        
        ^2) Faculty of Health, Engineering and Sciences, University of
        Southern Queensland, 
        Toowoomba, QLD, 4350, Australia and 
        ^3) Department of Applied Mathematics, Nizhny Novgorod State
        Technical University n.a. R.E. Alekseev, 
        Nizhny Novgorod,
        603950, Russia.
        
        
        We study solitary wave solutions of the fifth-order Korteweg–de
        Vries equation which contains, besides the traditional quadratic
        nonlinearity and third-order dispersion, additional terms
        including cubic nonlinearity and fifth-order linear dispersion, as
        well as two nonlinear dispersive terms. An exact solitary wave
        solution to this equation is derived and the dependence of its
        amplitude, width and speed on the parameters of the governing
        equation are studied. It is shown that the derived solution can
        represent either an embedded or regular soliton depending on the
        equation parameters. The nonlinear dispersive terms can
        drastically influence the existence of solitary waves, their
        nature (regular or embedded), profile, polarity, and stability
        with respect to small perturbations. We show, in particular, that
        in some cases embedded solitons can be stable even with respect to
        interactions with regular solitons. The results obtained are
        applicable to surface and internal waves in fluids, as well as to
        waves in other media (plasma, solid waveguides, elastic media with
        microstructure, etc.).
        47.35.Bb, 47.35.-i, 47.35.Fg
        
        Soliton solutions to the fifth-order Korteweg–de Vries
        equation 
        
        and their applications to surface and internal water waves
        
            K.R. Khusnutdinova^1), Y.A. Stepanyants^2, 3), and M.R. Tranter^1)Corresponding author:
        Yury.Stepanyants@usq.edu.au
        
            
        ===============================================================================================================================
        
        
        
        § INTRODUCTION
        
        
        The Korteweg–de Vries (KdV) equation
        
            ut + c ux + α u ux + β3ux = 0,
        
        is the well-known model for the description of weakly-nonlinear
        long waves in media with small dispersion (see, for instance,
        <cit.>). It
        is widely used in the theory of long internal waves where it
        describes astonishingly well the main properties of nonlinear
        waves, even when their amplitudes are not small (see, for
        instance, the reviews <cit.>). This
        is the simplest model that combines the typical effects of
        nonlinearity and dispersion, and provides stationary solutions
        describing both periodic and solitary waves. The KdV equation is
        completely integrable and possesses many remarkable properties,
        which can be found in the references cited above.
        
        At the same time, the KdV model cannot provide a detailed
        description of many important features of nonlinear waves observed
        in laboratory experiments, such as the non-monotonic dependence of
        solitary wave speed on amplitude, or the table-top shape of
        large-amplitude solitary waves <cit.>. To
        capture such properties, the first natural step is a
        straightforward extension of the KdV model by retaining the
        next-order nonlinear and dispersive terms in the asymptotic
        expansion of the solutions to primitive equations, for example the
        Euler equations with boundary conditions appropriate for
        oceanographic applications in the case of the ocean gravity waves.
        A rather general form of the extended KdV equation has been
        derived by many authors (in application to surface and internal
        waves as shown in Figure <ref> see, e.g.,
        <cit.>):
        
            ut + α u ux + β3ux + εα_1 u^2 ux + γ_1 u 3ux + γ_2 ux2ux + β_1 5ux = 0.
        
        This equation, written in the coordinate frame moving with the
        speed c, combines the quadratic (∼α) and cubic (∼α_1) nonlinear terms, linear dispersion of the
        3^rd (∼β) and 5^th (∼β_1) orders, and also higher-order nonlinear dispersion terms
        with coefficients γ_1 and γ_2; the parameter
        ε≪ 1 is presumed to be small.
        
        Particular cases of the fifth-order KdV equation, where some
        coefficients are zero, were also derived for plasma waves
        <cit.>, electromagnetic waves in discrete transmission
        lines <cit.>, gravity-capillary water waves
        <cit.>, waves in a floating ice sheet (see
        Ref. <cit.> and references therein).
        
        In general equation (<ref>) is not integrable, but for
        particular choices of coefficients it reduces to one of a set of
        equations that are completely integrable. These are the Gardner
        equation <cit.> (when β_1 = γ_1
        = γ_2 = 0) or its particular case the standard KdV/mKdV
        equation (when either α_1 = 0 or α = 0), as well as
        the Sawada–Kotera and Kaup–Kupershmidt equations (when α =
        β = 0) <cit.>. A comprehensive discussion of
        equation (<ref>) and its properties can be found in
        <cit.>.
        
        The coefficients of equation (<ref>) for surface gravity
        waves are
        
            c = √(gh),   α = 3c/2h,    α_1 = -3c/8h^2,   β = ch^2/6, β_1 = 19ch^4/360,   γ_1 = 5ch/12,   γ_2 = 23ch/24.
        
        For the gravity-capillary surface waves, as well as for internal
        waves in a two-layer fluid, the coefficients are presented in
        Appendix <ref>. All notations are shown in Figure 1.
        
        Unlike the KdV equation, the higher-order model (<ref>) is
        not a Hamiltonian equation and does not preserve the energy, in
        general. However, in the particular cases when it reduces to
        completely integrable models, it clearly becomes Hamiltonian.
        Besides those cases, there is one more particular case of
        α_1 = 0 and γ_2 = 2γ_1 when equation
        (<ref>) becomes Hamiltonian, but nonintegrable
        <cit.>. In the meantime, with the help of a
         near-identity transformation, it can be mapped 
        approximately into one of a number of Hamiltonian equations
        <cit.>. In particular
        the asymptotic near-identity transformation
        
            ũ = u + ε (a u^2 + b u_xx)
        
        maps equation (<ref>) into itself up to terms of o(ε), but with the new coefficients
        
            α̃_̃1̃ = α_1 - a α,   β̃_̃1̃ = β_1,   γ̃_̃1̃ = γ_1, γ̃_̃2̃ = γ_2 - 6 a β + 2 b α,
        
        where a and b are arbitrary parameters. If we choose the parameters a = γ_2 - 2 γ_1 /6 β and b = 0, then equation (<ref>) can be presented in the Hamiltonian form
        
            u_t = xδ H/δ u,
        
        where the Hamiltonian is H = ∫ℋx with the
        density
        
            ℋ = -1/6α u^3 + 1/2β u_x^2  - ε1/12α̃_̃1̃ u^4 + 1/2β̃_̃1̃ u_xx^2 - 1/2γ̃_̃1̃ u u_x^2.
        
        The Hamiltonian form provides conservation of the “mass”I_1 =
        ∫ u x, “wave energy”I_2 = ∫ u^2/2 x, and Hamiltonian I_3 = ∫ℋx. These
        conserved quantities are very useful in the development of
        asymptotic methods and perturbation techniques, as well as helping
        to control the accuracy of numerical schemes. Notice however, as
        has been shown in <cit.>, the formal Hamiltonians which
        follow from the approximate evolution equations are not usually
        the genuine Hamiltonians that can be derived from the primitive
        equations for small-amplitude wave perturbations and which agree
        with the physical energy conservation. Even in the classical KdV
        equation (<ref>) the Hamiltonian does not represent the
        genuine wave energy. To this end the “correct” KdV equation with
        the genuine Hamiltonian was derived in <cit.>; the
        corresponding equation is a particular case of equation
        (<ref>) with β_1 = 0.
        
        If the leading-order evolution equation is integrable, then the
        underlying physical system is said to be asymptotically integrable
        up to O(ε). It turns out that for the higher-order KdV
        equation (<ref>) with special choices of coefficients, it is
        possible to extend the asymptotic integrability, even up to
        O(ε^2)<cit.>. Indeed,
        the generic KdV equation is asymptotically reducible to the
        integrable equation by the nonlocal near-identity transformation
        
            ũ = u + ε a_1 u^2 + b_1 u_xx + c_1 u_x∫_x_0^x u x̃ + d_1 x u_t,
        
        where a_1, b_1, c_1, d_1 are arbitrary constants. This transformation can reduce (<ref>) either to the next member of  the KdV hierarchy, or even to the classical KdV equation, with accuracy up to O(ε^2). In particular, to transfer equation (<ref>) to the KdV equation one should choose coefficients in (<ref>) of the form
        
            a_1 = -18 β^2 α_1 + 2 α^2 β_1 + 3 αβγ_1/9 αβ^2,       b_1 = -6 β^2 α_1 - α^2 β_1 + αβγ_2/2 α^2 β, 
               c_1 = 4 αβ_1 - 2 βγ_1/9 β^2,       d_1 = - β_1/3 β^2.
        
        There are numerous other near-identity transformations; some of them are of special interest because they do not contain a secular term ∼ d_1 as in the formula (<ref>). In particular, the appropriately modified near-identity transformations reducing the higher-order KdV equation to the KdV equation have been successfully used to obtain particular solutions for the higher-order KdV equation from the known solutions of the KdV equation (e.g., the two-soliton solution extending the relevant KdV solution <cit.> and the undular bore solution <cit.>).
        
        In some particular cases when equation (<ref>) is non-integrable it possesses, nevertheless, stationary solitary-type solutions, which can be constructed either numerically or sometimes even analytically. One of the best known cases is the Kawahara equation which follows from equation (<ref>) when α_1 = γ_1 = γ_2 = 0<cit.>. This equation contains a rich family of solitary solutions including solitons with monotonic <cit.> and oscillatory tails <cit.>.
        
        Recently one more particular case of equation (<ref>) was considered, the so–called Gardner–Kawahara equation <cit.>, when only the nonlinear dispersive terms are absent (γ_1 = γ_2 = 0). Such a situation may occur, for example, in a two-layer fluid with surface tension between the layers. Solitary solutions for that equation were constructed numerically <cit.> and it was shown that among them there are “fat solitons”, similar to those seen in the Gardner equation, and solitons with oscillatory tails, such as in the Kawahara equation, as well as their combination – fat solitons with oscillatory tails.
        
        An analytical solitary wave solution of the non-integrable
        equation (<ref>) was found in <cit.> for a
        special set of coefficients, when it is not reduced to the
        Sawada–Kotera or Kaup–Kupershmidt equations. The obtained
        solution does not contain any free parameters and represents an
        example of a so-called “embedded soliton”. Embedded solitons
        co-exist with linear waves propagating with the same speed (they
        are “embedded” into the continuous spectrum of linear waves,
        whereas regular solitons can be called “gap” solitons, because
        they exist when the soliton speed belongs to a gap in the phase
        speed spectrum of a corresponding linearized system
        <cit.>). The term “embedded soliton” was introduced in
        Ref. <cit.>, although such solitons were known since 1974
        when the first analytical example was obtained in Ref.
        <cit.> and their stability was proven numerically in
        Ref. <cit.> (some information about that can be found
        also in Ref. <cit.>). As currently known
        <cit.>, embedded solitons can be both stable and unstable
        with respect to perturbations of small or even big amplitudes
        depending on the particular governing equation or set of
        equations.
        
        Nevertheless, the general problem of existence of solitary wave
        solutions of equation (<ref>) with arbitrary coefficients
        remains open so far, and this circumstance motivated our study.
        Besides the pure academic interest the problem is also topical in
        application to surface and internal waves of large amplitude in
        the ocean. Equation (<ref>) can be considered as the model
        equation capable of describing typical features of large-amplitude
        solitary waves with good accuracy (see, e.g.,
        <cit.> where it was shown that even its
        reduced version, the Gardner equation, provides solutions similar
        to those which can be constructed within the fully nonlinear Euler
        equations). We study the role of higher-order nonlinear dispersive
        terms (∼γ_1, γ_2) and their influence on the shape
        and polarity of solitary wave solutions. By means of the
        Petviashvili numerical method <cit.>
        we construct stationary solutions of equation (<ref>) and
        categorise them in terms of dimensionless parameters. We then
        numerically model non-stationary solutions using a pseudospectral
        scheme similar to that used in <cit.>. We show that these solutions demonstrate
        soliton-like properties in the course of their interaction with
        only minor inelastic effect. We also found an exact analytical
        solution to this equation in the general case without a
        restriction on its coefficients. The solution represents either
        the embedded or a regular (gap) soliton, depending on parameters.
        
        
        
        § DIMENSIONLESS FORM OF THE FIFTH-ORDER KDV EQUATION AND ITS GENERAL PROPERTIES
        
        
        To minimise the number of parameters in the problem let us present equation (<ref>) in dimensionless form using the change of variables
        
            τ = sα^3 t/εα_1√(s/εα_1 β),   ξ = α x √(s/εα_1 β),   υ = εα_1/α s u,
        
        where s = signα_1 β i.e. s = 1 if
        α_1 β > 0, and s = -1 if α_1 β < 0. After
        that the main equation (<ref>) can be presented in the
        divergent form:
        
            υτ + ξυ^2/2 + s υ^3/3 + 2υξ + B 4υξ + G_1/22υ^2ξ + G_2 -
            3G_1/2υξ^2  = 0,
        
        where
        
            B = α^2 s/α_1 ββ_1/β,    G_1 = s/α_1 βαγ_1,    G_2 = s/α_1 βαγ_2,
        
        (notice that s/(α_1 β) > 0). The divergent form
        immediately provides the “mass conservation” integral I_1 as
        defined before. Multiplying this equation by υ and
        integrating either over the entire axis ξ for solitary waves
        or over a period for periodic waves, we derive the “energy
        balance equation”
            I_2τ≡τ∫υ^2/2ξ =  G_1 - G_2/2∫υξ^3 ξ.
        
        This relationship denotes a distinguishing feature of the model
        (<ref>) which describes, in general, either time decay or
        growth of “wave energy”I_2 due to the presence of the
        nonlinear dispersive terms. Such conditionally defined “wave
        energy” is conserved either in the trivial case when G_1 = G_2 =
        0 or in the special case when G_2 = 2G_1<cit.>.
        For stationary waves described by even functions υ(ξ),
        the right-hand side of equation (<ref>) is zero, and their
        “energy” is also conserved for any values of G_1 and G_2.
        (Note that in the case of surface gravity waves, G_1 - G_2/2 =
        -3/2.)
        
        For waves of infinitesimal amplitude, υ→ 0,
        equation (<ref>) can be linearised. Looking for a solution
        in the form υ∼expi ω̃τ - κξ, we obtain the dispersion relation ω̃κ and phase speed V_ph(κ), in the coordinate
        frame moving with speed c, of the form
        
            ω̃ = -κ^3 + B κ^5,    V_ph(κ) ≡ω̃/κ = -κ^2 + B κ^4.
        
        The plot of the phase speed is shown in Figure
        <ref> for three typical values of the parameter
        B. Note that for surface water waves B = 57/5 > 0, therefore
        qualitatively the dispersion curve is similar to curve B = 1.
        The same is true for oceanic internal waves, as follows from the
        expressions for the coefficients of (<ref>) (see Appendix
        <ref>). In some cases for internal waves in laboratory
        tanks the coefficient B can be zero or negative.
        
        
        In the case of B ≤ 0 the phase speed is a monotonic function
        of κ, whereas for B > 0 it has a minimum, V_min =
        -1/(4B), at the point κ_c =  2 B ^-1/2 (see Figure
        <ref>). The concept of phase speed is very useful
        in understanding the process of interaction of a moving source
        (e.g. a solitary wave) with linear waves. In particular, if the
        speed of a source V_s is such that there is no resonance with
        any linear wave i.e. there is no intersection of the dashed line
        in Figure <ref> with the dispersion curve (e.g.
        when V_s < V_min), then the source does not lose
        energy for the wave excitation. Otherwise, in the case of a
        resonance (see the intersection of the dashed line with the curve
        for B = 1) the source, in general, can experience energy losses
        for the generation of a linear wave and, as a result, it gradually
        decelerates. Without external compensation of energy losses, such
        a source usually cannot move steadily. However, there are several
        examples of embedded solitons which can steadily propagate with
        the same speed as a linear wave, but not exciting it effectively.
        The physics of such a phenomenon has not been well understood yet;
        it requires further study which is beyond the scope of this paper.
        Therefore, the no-resonance condition V_s≠ V_ph can
        provide only an indication of when a solitary wave can most likely
        be expected.
        
        
        
        § STATIONARY SOLUTIONS OF THE FIFTH-ORDER KDV EQUATION
        
        
        Consider now stationary solutions to equation (<ref>) in the form of travelling waves depending only on one variable, ζ = ξ - V τ, where V is the wave speed. In this case equation (<ref>) can be reduced to an ODE and integrated once (the constant of integration is set to zero for solitary waves):
        
            B 4υζ + 2υζ - V υ + υ^2/2 + s υ^3/3 + G_1/22υ^2ζ + G_2 - 3G_1/2υζ^2 = 0.
        
        There are five independent parameters in this equation: B,
        G_1, G_2, s, and V, which determine the structure of a
        solitary wave. This equation actually splits into two independent
        equations with different properties; one equation with negative
        cubic term (s = -1), and another one with positive cubic term
        (s = 1). We will derive a particular soliton solution to this
        equation for both s=1 and s=-1. The case of s = 0 with a
        specific link between the parameters, G_2 = 2G_1, has been
        studied in <cit.>.
        
        As a first step, let us consider solitary wave asymptotics at
        plus/minus infinity. Assuming that soliton solutions decay at
        infinity, let us linearise equation (<ref>) (simply omit
        all nonlinear terms) and seek a solution of the remaining linear
        equation in the form υ∼expμζ.
        Substituting this trial solution into the linearised equation
        (<ref>), we obtain an algebraic equation for μ of the
        form (cf. <cit.>):
        
            B μ^4 + μ^2 - V = 0.
        
        The roots of this bi-quadratic equation are
        
            μ_1,2 = ±√(-1 + √(1 + 4BV)/2B),   μ_3,4 = ±√(-1 - √(1 + 4BV)/2B).
        
        Let us analyse the roots in detail (an alternative analysis of the
        roots in the plane of different parameters can be found in
        <cit.>). Firstly we consider the case when the
        parameter B is negative. For negative V we have 4BV
        > 0 and √(1 + 4BV) > 1, therefore the roots μ_1,2 are
        purely imaginary, and the roots μ_3,4 are real. Solutions
        corresponding to purely imaginary roots are not decaying and
        cannot represent solitary waves with zero amplitude at infinity.
        If V = 0, then μ_1,2 = 0 and the corresponding solutions
        do not decay at infinity.
        
        If 0 < V < -1/(4B), then we have √(1 + 4BV) > 0, and all four roots μ_1,2,3,4 are real. In this case soliton solutions are possible, with exponentially decaying asymptotics at infinity. Finally, if V > -1/(4B), then √(1 + 4BV) is complex; all roots are complex-conjugate in pairs μ_1,2 = ± (p_1 + iq_1), μ_3,4 = ± (p_2 + iq_2). Due to the presence of the real parts of the roots, p_1,2, there can be soliton solutions with oscillatory asymptotics. The decay rate of a solitary wave in the far field is determined by the root with the smallest value of |p_1,2|.
        
        Assume now that B is positive. It follows from a similar analysis of roots as was done for negative B that, for V < V_min≡ -1/(4B), solitary waves with oscillatory asymptotics are possible. For V_min < V < 0, the roots are purely imaginary, therefore no solitons with zero asymptotics can exist in this case. For V > 0 the roots
        μ_1,2 are purely real, then embedded solitons with exponential asymptotics are possible.
        
        In the particular case of B = 0, equation (<ref>) has two real roots corresponding to soliton solutions, provided that V > 0. These findings can be summarised with the help of a schematic diagram, as shown in Figure <ref>.
        
        It should be noted that the analysis of roots only predicts possible asymptotics of solitons provided that they exist, but it does not guarantee their existence. In particular, if B = 0, then soliton solutions with monotonically decaying exponential asymptotics could exist for any V > 0 (see case (b) in the diagram), but in fact they exist only for 0 < V < 1/6 (see, e.g., <cit.>).
        
        
         §.§ Solitary wave solution of the fifth-order KdV equation
        
        Let us consider a trial solution to (<ref>) in the form of a sech^2 solitary wave (in a similar way to what was done in <cit.>), taking the form
        
            υζ = A 2K ζ,
        
        where A is the soliton amplitude and K is the parameter
        determining its half-width Δ = 1/K (similar solution was
        constructed in <cit.> for the particular case of
        s = 0 and G_2 = 2G_1). By substitution of this solution into
        equation (<ref>) we obtain
        
            C_2 2K ζ + C_4 4K ζ + C_6 6K ζ = 0,
        
        where
        
            C_2    = 16BK^4 + 4K^2 - V, 
            
            C_4    = -120BK^4 + [2(G_1 + G_2)A - 6]K^2 + A/2, 
            
            C_6    = 120BK^4 + s A^2/3 - 2(G_2 + 2G_1)AK^2
        
        (cf. <cit.> where a slightly different approach was used). Equating the coefficients C_2 and C_4 to zero, we obtain
        
            V    = 4K^2  1 + 4BK^2  = 4 Δ^2 + 4B /Δ^4, 
            
            A    = 12K^2 1 + 20BK^2/1 + 4  G_1 + G_2  K^2 = 12/Δ^2Δ^2 + 20B/Δ^2 + 4  G_1 + G_2
            .
        
        Eliminating A from equation (<ref>) with the help of equation (<ref>), we obtain the quadratic equation for R ≡ K^2 of the form
        
            80B  G_1^2 + G_1 G_2 - 10Bs  R^2 + 4  2G_1^2 + 3 G_1 G_ 2 + G_2^2 - 5 BG_2 - 20Bs  R - 5B + 2G_1 + G_2 - 2s = 0.
        
        This equation has two roots, in general, and the corresponding expression for the half-width of a solitary wave is determined by
        
            Δ_1,2^2 = 40B  G_1 G_2 + G_1^2 - 10Bs /5B  G_2 + 4s  - 3G_1 G_2 - 2G_1^2 - G_2^2 ± 5B - G_1 - G_2 √( 2G_1 + G_2 ^2 - 40Bs).
        
        Thus, we see that a solitary wave solution in the form of (<ref>) does exist for a certain set of parameters B, G_1 and G_2. One of the obvious restrictions on the set of parameters is (we have two expressions as s can be positive or negative)
        
            2G_1 + G_2 ^2 - 40 B s ≥ 0  B ≥ - 2G_1 + G_2 ^2/40for s = -1, or B ≤ 2G_1 + G_2 ^2/40for s = 1.
        
        Below we consider a few particular cases and analyse the corresponding soliton solutions.
        
        
        
         §.§ The Gardner–Kawahara equation (G_1 = G_2 = 0)
        
        In this case the expressions for Δ, V and A simplify considerably. The half-width is given by
        
            Δ^2 = -40Bs/2s ±√(-10Bs).
        
        The soliton derived from (<ref>) can either be an embedded soliton or a regular soliton, dependent upon its speed as defined by (<ref>) and the value of B. Due to the presence of s in (<ref>) we consider two sub-cases: s = -1 and s = 1.
        
        
        
          §.§.§ The Gardner–Kawahara equation with s = -1
        
        In this case the only meaningful solution to (<ref>) is
        
            Δ^2 = 40B/-2 + √(10B),
        
        with B > 2/5 for a real solution. Then from equation
        (<ref>) it follows that V > 0. With B > 0 the phase speed
        dependence of linear waves on the wavenumber is shown in Figure
        <ref> by the upper line. Therefore, solution
        (<ref>) with V > 0 is the embedded soliton moving in
        resonance with a linear wave (however gap solitons can co-exist
        with the embedded soliton as will be shown below). Figure
        <ref> shows the dependence of the soliton parameters
        on B. The soliton amplitude monotonically increases with B,
        whereas its width and speed non-monotonically depend on this
        parameter. The minimum value Δ_min = √(32) is
        attained at B = 8/5, and maximum speed V_max = 5/32 occurs
        at B = 128/45. Soliton profiles for these two cases will be
        shown below in Figure <ref> in comparison with regular
        solitons, numerically constructed for the same values of B. The
        asymptotics of embedded solitons are in agreement with the
        prediction from the analysis of roots of the characteristic
        equation (<ref>) - see Figure <ref>b.
        
        
          §.§.§ The Gardner–Kawahara equation with s = 1
        
        In this case the soliton half-width takes the form
        
            Δ^2 = -40B/2 ±√(-10B), where B < 0.
        
        If the positive sign is taken in front of the square root, then the right-hand side of the expression is positive for all negative B, whereas if the negative sign is taken in front of the square root then the expression in the right-hand side is positive only under the restriction B > -2/5 (see Figure <ref> (a)).
        
        
        Then, from Figure <ref> it follows that for B <
        0 the solution (<ref>) represents a regular soliton if
        V > 0, and an embedded soliton if V < 0. The soliton velocity
        (<ref>) with Δ^2 given by equation (<ref>)
        is
        
            V = 4/Δ^4Δ^2 + 4B  = -5B + 8 ± 3 √(-10B)/50B.
        
        The analysis of this expression shows that, if the positive sign
        is chosen in front of the square root, then solution
        (<ref>) represents a regular soliton with V > 0, if
        -32/5 < B < 0, and an embedded soliton with V < 0, if B <
        -32/5. If the sign in front of the square root is negative, then
        solution (<ref>) represents a regular soliton with V >
        0 within the interval -2/5 < B < 0, whereas beyond this
        interval soliton solutions do not exist (see Figure
        <ref>b).
        
        As follows from equation (<ref>), solitons with positive sign in front of the square root in equation (<ref>) have positive polarity, and solitons with negative sign have negative polarity (see Figure <ref>c). In particular, when G_1 = G_2 = 0, we have
        
            A_1    = -360B/Δ_1^4 = -9/20B 2 - 5B + 2√(-10B), 
            
            A_2    = -360/Δ_2^4 = -9/20B^2 2 - 5B - 2√(-10B).
        
        Vertical dashed lines in Figure <ref> at B = -2/5
        show the limiting value of this parameter below which the negative
        polarity solitons cannot exist. Other vertical dashed lines in
        Figure <ref> at B = -32/5 show the boundary
        between the regular solitons with V > 0 and B
        > -32/5, and embedded solitons with V < 0 and B < -32/5.
        
        The horizontal dashed line in frame b) shows the limiting value of the embedded soliton speed V_lim = -1/10, when B → -∞, and the horizontal dashed line in frame c) shows the limiting value of the embedded soliton amplitude A_lim = 9/4, when B → -∞ (the width of the embedded soliton slowly
        increases when B → -∞ as Δ∼ 2√(-10B)).
        
        Thus, in the interval -2/5 < B < 0, two types of regular
        solitons can coexist with different widths, speeds, amplitudes,
        and polarities. When B < -2/5 only one soliton of positive
        polarity can exist which smoothly transfers into the embedded
        soliton when the parameter B passes through the threshold value
        B = -32/5 where the velocity vanishes. Figure <ref>
        illustrates the profiles of regular solitons of positive and
        negative amplitudes for three values of B in the indicated
        interval.
        
        
         §.§ A particular case of surface gravity waves (s = -1)
        
        
        Using the expressions for the coefficients of equation
        (<ref>) for gravity waves (see equation (<ref>)), we
        obtain the dimensionless parameters B = 57/5, G_1 = 10, G_2 =
        23 (notice that in this case G_2 = 2.3G_1, which is very close
        to the case G_2 = 2G_1 when the energy is conserved – see
        equation (<ref>)). For this set of parameters there is only
        one real root of equation (<ref>)
            Δ_1 = √(8436/√(2305) - 14)≈ 15.75.
        
        After that we find the amplitude and speed of a soliton in the form
        
            A = 3/148(51 - √(2305)) ≈ 6.06 ·
            10^-2,    V = 157√(2305) - 89/390165≈ 1.91
            · 10^-2.
        
        We see again that V > 0 and B > 0, therefore the soliton is in
        resonance with one of the linear waves as shown in Figure
        <ref> by the upper line. Hence, solution
        (<ref>) represents again the embedded soliton with
        exponential asymptotics in agreement with the prediction following
        from the analysis of roots of characteristic equation
        (<ref>) - see Figure <ref>b).
        
        In dimensional variables the amplitude, width and speed of the surface embedded soliton are
        
            A_d = 4hA,   Δ_d = h Δ/6,    V_d = 6V √(gh),
        
        where index d pertains to the dimensional variables, and h is
        the fluid depth. Setting h = 10 m, we obtain A_d ≈ 2.42m, Δ_d ≈ 26.25  m, V_d ≈ 1.14m/s (the total speed is c + V_d ≈ 9.9 m/s +
        1.14 m/s≈ 11.04 m/s). Thus, according to
        this solution, a soliton of small amplitude can exist on the
        surface of the water. In the meantime, the classical KdV theory
        predicts the existence of a soliton which has amplitude A_e
        ≈ 5.37 m and speed V_e ≈ 2.66 m/s
        at the same width Δ_d ≈ 26.25 m.
        
        The linear wave of infinitesimal amplitude propagating with the
        same speed as the embedded soliton (V_ph = c + V_d) has
        dimensional wavelength
        
            λ_d = π√(2β/V_d√(1 + 4 β_1 V_d/β^2) - 1 ).
        
        Substituting the coefficients of equation (<ref>) for pure gravity waves (see equation (<ref>) in Appendix <ref>) and h = 10 m, we obtain λ_d ≈ 32.5.
        
        As mentioned above, the embedded solitons can be stable even with
        respect to big perturbations. We will demonstrate in section
        <ref> that they can survive even after strong interactions
        with regular solitons. The problem of general evolution of
        embedded solitons under the action of small perturbations caused
        by medium inhomogeneity or energy dissipation is still an open
        problem and worth a further study. Some preliminary results can be
        found in <cit.>.
        
        
        
         §.§ A particular case when G_2 = 3G_1
        
        One more special case of G_2 = 3G_1 is worth considering because, in this case, the basic equation (<ref>) simplifies. In this simplified form, we can use the numerical code described in Appendix <ref>, based on the Petviashvili method and adapted for the solution of equation (<ref>) without the last term. This allows us to find the solutions numerically and compare the results obtained with the analytical solutions derived here. This case will also allow us to understand the role of nonlinear dispersion in the energy balance equation (<ref>). According to that equation, the energy is conserved on even solutions, and therefore stationary solutions in the form of solitary waves may exist. However the question is, what happens when solitary waves interact? We consider this problem in Section <ref> by direct numerical solution of the non-stationary equation (<ref>).
        
        In the case of G_2 = 3G_1 equation (<ref>) for the soliton half-width reduces to
        
            Δ_1,2^2 = 80B  2G_1^2 - 5Bs /5B  3G_1 + 4s  - 20G_1^2 + σ 5B - 4G_1 √(5  5G_1^2 - 8Bs ),
        
        where σ = ± 1. We plot the parameter plane (B, G_1) for the half-width, amplitude and speed in Figure <ref> for s = -1 and σ = ± 1, and similarly in Figure <ref> for s = 1 and σ = ± 1. Soliton solutions in the form of equation (<ref>) cannot exist if 5G_1^2 - 8Bs < 0, that is B < -5 G_1^2/8 for s = -1 and B > 5 G_1^2/8 for s = 1 (denoted as region 1a in Figures <ref> and <ref>) or if Δ_1,2^2 < 0 (region 1b). The regions where Δ_1,2^2 > 0 are designated by region 2. For the plot of soliton speed and amplitude, red regions indicate a positive quantity and blue regions indicate a  negative quantity (note that, in the figures for Δ, we assume the positive square root is taken). The difference between the cases for σ = -1 and σ = 1 is shown in the first line of each figure, where the size of region 1b differs in each case and therefore the range of values for which solitons can exist is different in each case.
        
        
        Thus, we see that nonlinear dispersion can affect the existence of
        soliton solutions, their nature (regular or embedded solitons),
        polarity, and, apparently, stability, as discussed in section
        <ref>.
        
        
        
        § NUMERICAL SOLUTIONS FOR STATIONARY SOLITONS
        
        Soliton solutions derived in the previous section represent particular cases of the wide family of stationary solutions containing a class of solitary waves. Solitary wave solutions of equation (<ref>) (or in the stationary case equation (<ref>)) can be constructed, in general, by means of one of the well-known numerical methods, e.g., the Petviashvili method <cit.> or Yang–Lacoba method <cit.>. As mentioned in the Introduction, in some particular cases soliton solutions were found analytically, in other cases numerically. In the particular case of the Gardner–Kawahara equation, when γ_1 = γ_2 = 0 in equation (<ref>) (or G_1 = G_2 = 0 in equation (<ref>)), soliton solutions were constructed and studied in <cit.>.
        
        Here we will construct a family of soliton solutions using the Petviashvili method for some particular cases to compare the numerical solutions with the analytical solutions derived in the previous section. The numerical method in application to equation (<ref>) is described in Appendix <ref>. We consider two particular cases when G_1 =
        G_2 = 0 and when G_2 = 3G_1. In the latter case we will see the influence of nonlinear dispersion on the shape and polarity of solitary waves.
        
        
        
         §.§ The Gardner–Kawahara equation (G_1 = G_2 = 0)
        
        
        We showed in Section <ref> that embedded solitons of the form (<ref>) can exist in the Gardner–Kawahara equation under certain restrictions on the value of B, and these restrictions are different for s = -1 and s = 1. We therefore split the following analysis into two subsections.
        
        
        
          §.§.§ Numerical solutions in the case of negative cubic term (s = -1)
        
        In this case soliton solutions in the form of equation
        (<ref>) can exist only for B > 2/5 and represent the
        embedded soliton. However, soliton solutions in different forms
        can exist both for positive and negative B. For positive B
        such solutions can be constructed numerically for V < V_min≡ -1/(4B) (see Figure <ref>). In Figure
        <ref> we present two families of numerical solutions. In
        panel (a) one can see the analytical solution for the embedded
        soliton (<ref>) (line 1) when B = 8/5 and the soliton
        has the minimum width Δ_min = √(32) and speed V =
        3/20 (see Subsection <ref>). Lines 2, 3, and 4
        show numerically constructed regular solitons for the same value
        of B and V = -0.25, -0.5, and -1 respectively.
        
        
        In panel (b) one can see the analytical solution for the embedded
        soliton (<ref>) (line 1) when B = 128/45 and the
        soliton has the maximal speed V = 5/32 and width Δ =
        16√(2/15) (see Subsection <ref>). Lines 2, 3,
        and 4 show numerically constructed regular solitons for the same
        value of B and V = -0.25, -0.5, and -1 respectively.
        
        When B < 2/5 analytical solutions of equation (<ref>) are
        unknown, however they can be constructed numerically. We managed
        to construct such solutions within a relatively narrow range of
        parameters. In particular, when B = 1/5, solitary-type solutions
        appear in the form of oscillatory wave trains – see Figure
        <ref>a (similar solutions were constructed
        numerically in <cit.>). When V approaches -1.25
        from below, solitons look like envelope solitons of the non-linear
        Schrödinger (NLS) equation (see, e.g. <cit.>). Then, when V
        decreases, the soliton shape smoothly changes and represents a
        negative polarity soliton with oscillatory tails as shown by line
        2 in Figure <ref>a. We did not manage to construct
        numerical solutions for V < -1.7.
        
        
        When B < 0 regular solitons exist for V > 0 and have
        bell-shaped profiles as shown in Figure <ref>b for
        B = -1 and several values of V. However, they exist only
        within a relatively narrow range of V between 0 and 0.17.
        Similar situation occurs for other negative values of B.
        
        
        
          §.§.§ Numerical solutions in the case of positive cubic term (s = 1)
        
        
        In this case solutions in the form of (<ref>) can exist
        only for negative B; they can be in the form of either regular
        solitons or embedded solitons (see Subsection
        <ref>). However soliton solutions in different
        forms can exist, in principal, both for positive and negative B.
        Nevertheless we did not manage to construct soliton solutions
        numerically for positive B. This requires further investigation.
        
        In the case of negative B there are two families of regular
        solitons of positive and negative polarities. Typical examples are
        shown in Figure <ref> for B = -0.1. Line 1 represents
        the analytical solution (<ref>) with V = 2.1 and line 2
        is another analytical solution with V = 0.9 as per equation
        (<ref>). However,  we failed to reproduce these solutions
        numerically by means of the Petviashvili method. Lines 3 and 5 in
        Figure <ref> correspond to solutions found using the
        numerical scheme for V = 2.1 and V = 0.9 respectively (in
        comparison to the analytical solutions of line 1 and line 2). All
        these solitons have exponentially decaying asymptotics at
        infinity.
        
        To explain the difference in the solutions, we refer to equation
        (<ref>). As follows from the analysis of this equation,
        all its roots are real when 0 < V < -1/(4B) (which is the case
        for Figure <ref>). This means that exponential decay
        of a solitary wave can be controlled by one of two real roots at
        plus infinity and one of the other two roots at minus infinity.
        Apparently, soliton solutions obtained analytically and
        numerically correspond to different roots of characteristic
        equation (<ref>).
        
        Then, it follows from numerical solutions that when the soliton
        speed increases, so does the amplitude, but the soliton profile
        becomes non-monotonic; this is again in agreement with the
        qualitative analysis shown in Figure <ref>a for B
        < 0. Lines 4 and 6 illustrate such solutions for V = 100.
        
        
        We did not manage to construct solitons of negative polarity with B = -0.1 and V < 1.95 even when the starting solution for the iteration scheme (see Appendix <ref>) was chosen in the form of the analytical solution (<ref>). However, for V ≥ 1.95, solitons of negative polarity were constructed numerically and are shown in Figure <ref>. They represent almost mirror reflections of solitons of positive polarity. A similar situation occurs in the case of the Gardner equation with positive cubic nonlinearity, when solitons of positive polarity can exist for all amplitudes from 0 to infinity, whereas solitons of negative polarity cannot exist, if their amplitudes are less then some critical value. Below the critical amplitude breathers can exist instead (see, e.g., <cit.> and references therein).
        
        
        
         §.§ The particular case of equation (<ref>) with G_2 = 3G_1
        
        
        Another particular case we study is G_2 = 3G_1. The basic
        equation (<ref>) in this case contains four independent
        parameters B, G_1, s, and V, which  determine the
        structure of solitary waves. This case is convenient from the
        numerical point of view as the iteration scheme is simpler than
        for other cases of G_1 and G_2. Another point of interest in
        this case is the fact that wave energy is not conserved in general
        – see equation (<ref>) (whereas for even solutions the
        energy is conserved), therefore this allows us to understand the
        role of nonlinear dispersion. In the parameter plane (G_1,B)
        there are zones where regular solitons of different polarity can
        exist. The boundaries of these zones are fairly complicated and we
        demonstrate some typical soliton solutions belonging to different
        zones, for the cases when s = -1 and s = 1 as before.
        
        
        
          §.§.§ Numerical solutions in the case of negative cubic term (s = -1)
        
          * B = -1
        
        Consider first the case when B = -1 in equation (<ref>).
        In this case regular solitons can exist for V > 0. In Figure
        <ref> we illustrate the structure of a solitary
        wave when G_1 varies. Line 1 in both panels show the reference
        case when G_1 = 0. If this parameter becomes negative, the
        solitons become wider, their amplitudes slightly decrease, and the
        tops become flatter (see lines 2 and 3 in panel a). However, when
        G_1 approaches -9, the soliton profiles become non-monotonic
        with oscillations on the top (see line 4 in panel a) (similar
        solutions were constructed numerically in
        <cit.>). For smaller values of G_1 < -9 we were
        unable to obtain numerical solutions (apparently, they do not
        exist for such a set of parameters).
        
        When G_1 becomes positive and increases the solitons become narrower, their amplitudes slightly increase first, but then they monotonically decrease when G_1 exceeds 4 (see lines 1, 2, and 3 in panel b). For relatively large values of G_1 soliton tails become non-monotonic with negative minima on the profile (see line 4 in panel b).
        
          * B = 1
        
        When B=1 solutions in the form of regular solitons can exist only
        for negative V < V_min≡ -1/(4B) (see after equation
        (<ref>)). In Figure <ref>a we show the
        structure of solitary waves for V = -0.5 and different values of
        G_1. Solitons in this case have negative polarity and
        oscillating tails. Line 1 corresponds to the case when G_1 = 0.
        If this parameter becomes negative, the solitons become wider and
        their amplitudes increase (see line 2 in panel a). For G_1 < -3
        we were unable to construct numerical solutions (apparently, they
        do not exist for such a set of parameters). When G_1 becomes
        positive and increases, the solitons become narrower and smaller
        (see lines 3 and 4 in panel a).
        
        
        In Figure <ref>b we show the structure of solitary waves for V = -5 and different values of G_1. Solitons in this case can have both negative and positive polarity; they can have slightly oscillating tails or non-monotonic aperiodic tails. Line 1 corresponds to the case when G_1 = 0. If this parameter becomes negative and varies from zero to G_1 = -1.88, the solitons remain qualitatively the same, but become wider and their amplitudes slightly increase as shown in panel b) by line 2. When G_1 further decreases and becomes less than -1.89, the solitons abruptly change their polarity, become taller and narrower with well-pronounced negative minima (see line 3). Further increases in G_1 result in soliton profiles that remain qualitatively the same, but their amplitudes decrease.
        
        When G_1 becomes positive and increases, soliton profiles remain similar to the shape of the soliton for G_1 = 0, but the solitons become narrower and of smaller amplitude (see, e.g., line 4 in panel b).
        
        
        
          §.§.§ Numerical solutions in the case of positive cubic term (s = 1)
        
        
        In this case the analytical solution (<ref>) for B > 0
        can exist only for V > 0 representing an embedded soliton (see
        Figure <ref> and Figure <ref>),
        whereas solutions in the form of regular solitons, apparently,
        cannot exist; we were unable to construct such solutions
        numerically.
        
        For B < 0 the analytical solution (<ref>) exists in the
        form of a regular soliton for V > 0 and embedded soliton for V
        < 0 (see Figure <ref> and Figure
        <ref>). Embedded solitons cannot be reproduced by
        means of the Petviashvili method (see Appendix <ref>),
        therefore we constructed only regular solitons numerically for the
        particular case of B = -1 and V > 0. The structure of these
        solitons depends on their speed and is illustrated by Figure
        <ref> for two particular values of V. In panel (a)
        we show a few typical soliton profiles for V = 0.1 and several
        values of G_1. Line 1 shows the reference case when G_1 = 0.
        If G_1 becomes negative, the solitons become wider and their
        amplitudes slightly decrease (see line 2). For G_1 < -12 we did
        not obtain soliton solutions.
        
        When G_1 becomes positive, the soliton amplitude slightly
        increases first (see line 3), but then it decreases, and solitons
        become narrower (see line 4). For sufficiently large G_1
        solitons profiles become non-monotonic, so that negative minima
        appear in the profiles (see line 5).
        
        
        In panel (b) we show other typical soliton profiles for V = 1 and several values of G_1. Line 1 shows the reference case when G_1 = 0. If G_1 becomes negative, the solitons become wider and their amplitudes slightly decrease (see line 2). When G_1 passes through some critical value between -4.63 and -4.64 the soliton polarity abruptly alters from positive to negative (see line 3). Then, when G_1 further decreases, soliton profiles remain qualitatively similar to what is shown by line 3, but their amplitudes become smaller (see line 4). When G_1 becomes positive soliton profiles remain qualitatively similar to line 1, but their amplitudes gradually decrease and well-pronounced minima appear on both sides of the crests (see line 5).
        
        The solutions constructed numerically do not reproduce the analytical solution (<ref>). The reason is the same as discussed in Section <ref> i.e.  solutions with different asymptotics can coexist for the same set of parameters, and the numerical scheme, apparently, converges to only one of them which is different from solution (<ref>).
        
        
        
         §.§ Stationary solutions for two particular cases of equation (<ref>)
        
        In this subsection we briefly consider particular exact solutions of equation (<ref>). In the first case we set s = G_1 = G_2 = 0 to reduce the general equation (<ref>) to the generalised Kawahara equation containing both third- and fifth-order derivatives. The exact solution to this equation was obtained for the first time by  Yamamoto & Takizawa <cit.> (for further references see also <cit.>)
        
            υζ = -105/169B4x - V τ/Δ,
        
        where V = -36/169B, Δ = √(-52B) and B < 0. There
        are no free parameters in this solution; the amplitude, speed and
        width of Yamamoto–Takizawa (YT) soliton (<ref>) are
        determined by the coefficients of the generalised Kawahara
        equation. The soliton moves with a positive speed, whereas linear
        waves propagate with negative phase speeds, therefore it is a
        regular soliton; its profile is shown in Figure <ref>
        by line 1 for B = -1. It was easily reproduced numerically with
        the help of the Petviashvili method when the speed was chosen in
        accordance with the formula V = -36/169, B ≈ 0.213.
        
        By means of Petviashvili's method we constructed a family of
        soliton solutions with the fixed value of parameter B = -1. All
        solitary solutions of this family are qualitatively similar to the
        YT soliton. In particular, line 2 shows the numerical solution for
        V = 0.25, and line 3 for V = 0.15. It was discovered that the
        soliton amplitude decreases as the speed decreases. At small
        amplitudes the soliton profile becomes indistinguishable from the
        profile of the KdV sech^2-soliton of the same amplitude. This is
        illustrated by Figure <ref> where the numerically
        obtained line 3 practically coincides with line 4, which
        represents the KdV soliton of the same amplitude. Apparently
        within this equation there is a continuous family of solitary wave
        solutions whose profiles depend on their amplitude, and the YT
        soliton is just one particular of the representatives of this
        family.
        
        
        Because the generalised Kawahara equation is non-integrable, one
        can expect that soliton interactions are inelastic, i.e. in the
        process of soliton collisions they radiate small-amplitude
        trailing waves and, as a result, change their parameters. This
        will be confirmed in Section <ref>.
        
        In the second case we present the soliton solution to the
        Kaup–Kupershmidt equation <cit.>. This equation is a particular case of equation
        (<ref>) with the following coefficients: α = β =
        0, εα_1 = 180, εγ_1 = 30,
        εγ_2 = 75, and εβ_1 = 1. With
        such a set of coefficients the equation is completely integrable,
        and its soliton solution has a slightly unusual form:
        
        
            υ(ξ, τ) = 3A 4 2ζ / Δ -
            4ζ / Δ/ 2 +  2ζ / Δ^2,
        
        where A is the soliton amplitude (a free parameter), Δ = 1/√(6A) is the soliton width, ζ = ξ - V τ, and V = (6A)^2 is the soliton speed.
        
        The profile of the Kaup–Kupershmidt (KK) soliton of a unit
        amplitude is shown in Figure <ref>. As one can see, its
        speed is positive, whereas waves of infinitesimal  amplitude
        within the Kaup–Kupershmidt equation have negative phase speeds
        (in the moving coordinate frame). Therefore, the KK soliton is a
        regular soliton too.
        
        
        In contrast to the previous case, the Kaup–Kupershmidt equation
        is completely integrable, therefore soliton interactions are
        elastic, and to a certain extant, trivial. This means that after
        interaction solitons completely restore their original parameters
        and do not radiate small amplitude perturbations. Therefore,
        solitons remain the same as they were before interaction, and the
        only traces of interaction are their shifts in space and time,
        exactly as in the interaction of elastic particles.
        
        
        
        § NUMERICAL STUDY OF SOLITON INTERACTIONS
        
        It was anticipated in <cit.> that equation (<ref>) would be studied numerically to confirm existence and robustness of solitary wave solutions, investigate whether they emerge evolutionary from the arbitrary initial pulse-type perturbations, how do they behave under collisions, whether “they emerge unscathed as true solitons, or is there a small, but nonzero nonelastic effect”. Since that time such investigation was not carried out to the best of our knowledge. Here we will try to illuminate these issues.
        
        To solve equation (<ref>) numerically, we apply a
        pseudospectral technique similar to that used in <cit.>. The equation is solved in the
        Fourier space using a 4 order Runge–Kutta method for time
        stepping, while the nonlinear terms are calculated in the real
        space and transformed back to the Fourier space for use in the
        Runge–Kutta scheme. To remove the aliasing effects, we use the
        truncation 2/3-rule by Orszag in Boyd <cit.>. See Appendix
        <ref> for the description of the numerical scheme.
        
        To generate solitons for a given set of parameters, an initial
        pulse was taken as the initial condition. This pulse was taken in
        the form a standard KdV soliton i.e. of the form
        
            u(x,0) = P 2√(P/12)x/L,
        
        where P corresponds to the amplitude and L is a factor used to
        change the width of this initial pulse. In each case, we can
        control the number of solitons produced, and their amplitudes, via
        the parameters P and L.
        
        To calculate the interaction of regular or embedded solitons, we
        generate regular solitons of the required amplitude using the
        method described above. Once the solitons generated from the pulse
        are well separated, we extract them from the solution. To study
        the interaction of these regular solitons with other regular
        solitons or embedded solitons, we generate an initial condition
        using the extracted soliton and either another extracted soliton
        (for the interaction of regular solitons) or with the embedded
        soliton found analytically, so that their interactions could be
        studied. This extraction was performed so that any radiation
        emerging from the initial pulse would not interfere with the
        collision. In each of the numerical cases considered below, the
        value of P and L are stated for each regular soliton generated
        via this method. Furthermore, we state if the solitons used in the
        proceeding calculations are found analytically or generated from a
        pulse.
        
        With the help of this numerical method we studied interactions of solitary waves with different parameters and different coefficients of the governing equation (<ref>). First of
        all we found that the embedded soliton propagates in all cases, with minimal loss of energy. We have calculated the change of “wave energy”I_2 = ∫ (u^2/2) x (see Introduction) as
        
            Δ I_2 = I_2(t) - I_2(0)/I_2(0),
        
        where I_2(0) is the initial wave energy and I_2(t) is the wave energy at time t.
        
        Using this numerical method with periodic boundary conditions, we
        obtained Δ I_2 = 1.2 · 10^-13 and 2.4 ·
        10^-13 for the embedded solitons in the cases when nonlinear
        dispersion is present, whereas for the regular solitons in these
        cases we obtained Δ I_2 = 1.2 · 10^-5 and 5.2 ·
        10^-6. In the case of the regular solitons, as they were
        generated by a pulse-like initial condition, fast moving radiation
        was generated that would re-enter the domain and interfere with
        the main wave structure. To diminish this effect, we applied “a
        sponge layer” to the solution domain to absorb this radiation and
        prevent it re-entering the solution domain, as detailed in
        Appendix <ref>. This accounts for the lower accuracy in
        the wave energy conservation. It is worth noting that when the
        embedded soliton is perturbed, the energy is no longer conserved
        and therefore the solution eventually breaks down, except in the
        case when there is no nonlinear dispersion (see Section
        <ref> below).
        
        For the regular solitons in all cases, they steadily propagate
        without loss of energy even in the cases when G_2  2G_1. As
        has been mentioned above (see the paragraph after equation
        (<ref>)), in the case of propagation of a stationary wave
        described by an even function, the right-hand side of equation
        (<ref>) vanishes, and wave energy I_2 is conserved.
        However, an interesting question arises about the energy
        conservation in the process of soliton interaction when
        G_2  2G_1. Below we present the results of our numerical
        study of equation (<ref>) with negative cubic
        nonlinearity in the cases when (i) G_1 = G_2 = 0 and
        (ii) G_2 = 3G_1. Equation (<ref>) with positive cubic
        nonlinearity can be studied in a similar way, but such exercises
        require much more computational resources, because soliton
        amplitudes are limited in this case and can be very close to each
        other, whereas their speeds are relatively small, therefore
        soliton interactions take a very long time to compute.
        
        
        
         §.§ The Gardner–Kawahara equation (G_1 = G_2 = 0)
        
        
        As shown in subsection <ref>, there are families of
        regular solitons for positive and negative B, some of which can
        co-exist with the embedded soliton (<ref>)– see Figures
        <ref> and <ref>. Here we present (i) an
        example of pulse disintegration into a number of regular solitons
        (Figure <ref>); (ii) interaction of regular solitons
        (Figure <ref>); and (iii) interaction of regular and
        embedded solitons (Figure <ref>). Figure
        <ref> illustrates that regular solitons with
        non-monotonic profiles asymptotically appear from pulse-type
        initial perturbations in the process of its disintegration.
        Apparently, such solitons can form bound states, the bi-solitons,
        triple-solitons, or multi-solitons, i.e., stationary moving
        formations consisting of two or more binding solitons
        <cit.>,
        but we did not study this phenomenon in our paper.
        
        
        The interaction of two solitons as shown in Figure <ref> demonstrates that the solitons survive after the collision, but a residual small wave packet is generated in the trailing wave field. This clearly indicates that the soliton collision is inelastic.
        
        The most fascinating is Figure <ref> which
        demonstrates (seemingly for the first time) that the embedded
        soliton can survive after  interaction with a regular soliton. The
        interaction is obviously inelastic, so that small disturbances
        appear both in front and behind the embedded soliton. Thus, we see
        that it survives even after collision with a regular soliton, and
        for much longer times (up to t = 700) the embedded soliton keeps
        its identity.
        
        
         §.§ Interactions of solitary waves when G_2 = 3G_1
        
        
        Following the same steps as in subsection <ref>, we consider the cases when the parameters are (i) s = -1, B = -1, G_1 = 4, and G_2 = 12 (see Figure <ref>b) and (ii) s = -1, B = 1, G_1 = -1.88, and G_2 = -5.64 (see Figure <ref>b). First of all we observed steady propagation of solitons in both of these cases. Then we observed the emergence of a number of solitons from an initial pulse with larger amplitude and width; this is illustrated in Figure <ref> for case (i) and Figure <ref> for case (ii).
        
        
        
        Finally we studied the interaction of these solitons. When
        colliding two regular solitons, we observed that after the
        interaction both the solitons survive, but some portion of their
        energy converts into a small wave packet generated in the trailing
        wave field; this is shown in Figure <ref> for case (i)
        and Figure <ref> for case (ii). This is typical for
        the inelastic interaction.
        
        
        
        For the collision of a regular soliton with an embedded soliton,
        in case (i) we observed that only the regular soliton survives the
        collision and an intense wave packet is generated in the trailing
        wave field (see Figure <ref>).
        
        In contrast to that, in case (ii), the regular and embedded
        solitons both survive after the collision and a wave packet is
        generated in front of the embedded soliton (see Figure
        <ref>). In both cases we conclude that the collision
        is inelastic. It is worth noting in this case that our numerics
        were not stable beyond the time considered in the calculation; it
        may be stable for a smaller spatial discretisation, however the
        corresponding time discretisation becomes very small and the
        calculations take a very long time.
        
        
        
         §.§ The generalised Kawahara equation
        
        As we defined in Section <ref> above, the generalised
        Kawahara equation is a particular case of equation (<ref>)
        with the coefficients s = G_1 = G_2 = 0. We analysed the
        stationary solitary solutions, one of which is the YT soliton
        (<ref>). Here we show that solitary waves emerge from
        pulse-type initial perturbations within the generalised Kawahara
        equation. As the initial condition the sech^2 pulse was chosen.
        Figure <ref> illustrates an example of pulse
        disintegration onto two solitary waves accompanied by small
        residual wave train (not visible in the figure).
        
        The soliton interaction is demonstrated in Figure
        <ref>. The initial condition was chosen as the YT
        soliton (1) and the numerically obtained solitary wave (2). We see
        that solitary waves survive the collision and appear after that
        with almost the same amplitudes. However, a small wave train
        appears in the result of interaction, evidence that the
        interaction is inelastic.
        
        
        § CONCLUSION
        
        In this paper we have studied the properties of soliton solutions of the fifth-order KdV equation (<ref>) which is used to describe surface and internal gravity waves, as well as appearing in other applied areas. Using the changes of independent and dependent variables this equation has been reduced to the dimensionless form (<ref>) with the minimum number of independent parameters, only three for positive and negative cubic nonlinearity. In the theory of nonlinear internal waves both these cases of positive and negative cubic nonlinearity can occur depending on the stratification of the fluid <cit.>. In some particular cases equation (<ref>) reduces to well known equations, among which there are Gardner <cit.>, Kawahara <cit.> equations and their generalisation <cit.>, Sawada–Kotera and Kaup–Kupershmidt equations <cit.>.
        
        Equation (<ref>) is non-integrable, in general, and it
        does not provide conservation of the “wave energy". However, it
        permits the existence of solitary wave solutions both with
        monotonic and non-monotonic profiles and even with oscillatory
        asymptotics, which suggests the existence of more complicated
        structures in the form of stationary multi-solitons
        <cit.>. Following
        <cit.>, we have derived an exact soliton solution in
        the general case and have shown that in some cases this solution
        can represent a regular soliton, whereas in others it represents
        an embedded soliton whose speed coincides with the speed of some
        linear wave. Our numerical simulations have confirmed that such
        solitons can propagate without loss of energy, if there are no
        perturbations in the form of other waves or medium inhomogeneity,
        or dissipation. We have identified the areas on the plane of
        parameters where the derived soliton solution exists as the
        embedded or regular soliton and found numerically a number of
        soliton solutions of equation (<ref>) with different sets
        of governing parameters.
        
        In the last section we have studied numerically the emergence of
        solitons from arbitrary pulse-type initial conditions and have
        demonstrated that regular solitons are generated from the initial
        conditions. However in the course of interaction such solitons
        produce small-amplitude trailing waves which evidences that the
        interaction has inelastic character. Moreover, the solitons with
        non-monotonic profiles can form stationary or non-stationary bound
        states similar to those experimentally observed in
        <cit.>.
        
        One particularly interesting observation emerging from our study
        is the stability of embedded solitons with respect to interactions
        with other waves (e.g., the regular solitons, like in this paper).
        This sheds some additional light on the stability problem and
        nature of embedded solitons and indicates that they could be
        observed in natural and laboratory environments. In particular, in
        the absence of nonlinear dispersion, the embedded soliton survives
        the interaction with no visible loss of amplitude. On the other
        hand, the results of our numerical investigations suggest that in
        some cases embedded solitons, apparently, can transfer to 
        radiating solitons under the influence of external perturbations
        (the radiating solitons are quasi-stationary long-living solitary
        waves permanently radiating from one side small-amplitude linear
        waves – see, e.g., <cit.> and
        references therein). We do not touch this interesting possibility
        in the present paper, but it can be a theme for further study.
        
        Finally, it would be interesting to extend the present study to
        the study of the respective ring wave counterparts (see <cit.> and references therein). The  solitary wave
        solutions studied in our present paper provide meaningful
        “initial conditions" for numerical experimentation with the
        amplitude equation describing ring waves.
        
         The authors are thankful to T. Marchant for
        useful discussions and valuable comments. This research was
        initiated within the framework of the Scheme 2 grant of the London
        Mathematical Society (LMS), in November–December 2015. K.R.K. and
        Y.A.S. are grateful to the LMS for the support. Y.A.S.
        acknowledges the funding of this study from the State task program
        in the sphere of scientific activity of the Ministry of Education
        and Science of the Russian Federation (Project No.
        5.1246.2017/4.6).  K.R.K. acknowledges the support received during
        her stay at ESI, Vienna in 2017, where some parts of this work
        were finalised. M.R.T acknowledges the support of the
        Engineering and Physical Sciences Research Council (EPSRC).
        
        
        § COEFFICIENTS OF THE HIGHER-ORDER KDV EQUATION FOR WATER WAVES
        
        
        The coefficients of equation (<ref>) for gravity surface waves are given in (<ref>). Here we present the coefficients for internal waves in two-layer fluid derived in <cit.>. All notations are shown in Figure <ref>.
        
        
            c    = √((ρ_2 - ρ_1) gh_1h_2/ρ_1h_2 + ρ_2h_1),   α = 3c/2h_1h_2ρ_2h_1^2 - ρ_1h^2_2/ρ_1h_2 + ρ_2h_1,   β = h_1h_2/6 c ρ_1h_1 + ρ_2h_2 - 3 σ/c^2/ρ_1h_2 + ρ_2h_1, 
            α_1    = -3c/8h_1^2h_2^2ρ_2^2h_1^4 + ρ_1^2h_2^4 + 2ρ_1ρ_2h_1h_2 4h_1^2 + 7h_1h_2 + 4h_2^2/ρ_1h_2 + ρ_2h_1^2, 
            β_1    = ch_1h_2/90ρ_1 ρ_2  h_1^4 + h_2^4 + 15h_1^2h_2^2/2  +  19/4  h_1h_2 ρ_1^2h_1^2 + ρ_2^2h_2^2 - S_1/ρ_1h_2 + ρ_2h_1^2, 
            γ_1    = c/125h_1h_2 ρ_2^2h_1 - ρ_1^2h_2  + ρ_1 ρ_2  h_1 - h_2  7h_1^2 + 9h_1h_2 + 7h_2^2  - 2S_2/ρ_1h_2 +ρ_2h_1^2, 
            γ_2    = c/2423h_1h_2 ρ_2^2h_1 - ρ_1^2h_2  + ρ_1ρ_2  h_1 - h_2  31h_1^2 + 39h_1h_2 + 31h_2^2  + 2S_2/5/ρ_1h_2 + ρ_2h_1^2,
        
        where
        
            S_1 = 5h_1h_23σ/2c^2ρ_1h_1 + ρ_2h_2 + 3σ/2c^2,    S_2 = 3σ/2c^2ρ_2h_1^2 - ρ_1h_2^2.
        
        
        In particular, when ρ_1 = 0, we obtain the coefficients of
        equation (<ref>) for surface gravity-capillary waves on a
        thin liquid layer:
        
            c = √(gh_2),   α = 3c/2h_2,    α_1 = -3c/8h_2^2,   β = ch_2^2/6 1 - 3 σ/c^2 ρ_2h_2, 
            β_1 = ch_2^4/1819/20 - 3σ/2c^2 ρ_2h_2 1 + 3 σ/2c^2 ρ_2h_2, 
            γ_1 = 5ch_2/12 1 - 3 σ/5c^2 ρ_2h_2,   γ_2 = 23ch_2/24 1 + 15 σ/23c^2 ρ_2h_2.
        
        
        § PETVIASHVILI'S METHOD
        
        
        Let us make a Fourier transform of equation (<ref>) with respect to the variable ζ denoting the Fourier image of the function υ(ζ) by F̂(υ):
        
            B κ^4 - κ^2 - V F̂(υ) = -s/3F̂(υ^3) - 1/2 1 - G_1 κ^2 F̂(υ^2) - 1/2 G_2 - 3G_1 F̂υ'^2 ,
        
        where κ is the parameter of the Fourier transform (the dimensionless wave number).
        
        If we multiply equation (<ref>) by F̂(υ) and integrate it with respect to κ from minus to plus infinity, we obtain the equality
        
            ∫_-∞^+∞ B κ^4 - κ^2 - V F̂(υ) ^2 dκ = -1/2∫_-∞^+∞2s/3F̂(υ^3) +  1 - G_1κ^2 F̂(υ^2) +  G_2 - 3G_1 F̂υ' ^2 F̂(υ) dκ.
        
        If υ(ζ) is an exact solution of equation (<ref>) and F̂(υ) is its Fourier image satisfying equation (<ref>), then it follows from equation (<ref>) that the quantity M, dubbed the stabilising factor and defined below, should be equal to one:
        
            M [υ] = -2 ∫_-∞^+∞ B κ^4 - κ^2 - V F̂(υ) ^2 dκ/∫_-∞^+∞ (2s/3) F̂(υ^3) +  1 - G_1 κ^2 F̂(υ^2) +  G_2 - 3G_1 F̂υ' ^2 F̂(υ) dκ.
        
        However, in general, if υ(ζ) is not a solution of equation (<ref>), then M[υ] is some functional of υ. In the spirit of the Petviashvili method, let us construct the iteration scheme (for details see <cit.>):
        
            F̂(υ_n+1) = -1/2 M^r [υ_n] (2s/3) F̂(υ_n^3) +  1 - G_1 κ^2 F̂(υ_n^2) +  G_2 - 3G_1 F̂υ_n'^2 /B κ^4 - κ^2 - V,
        
        where the factor M is used to provide a convergence of the iterative scheme (otherwise the scheme is not converging), and the exponent r should be taken in the range r = [3/2, 2]. As has been shown in <cit.>, r = 3/2 provides the fastest convergence for pure cubic nonlinearity, whereas r = 2 provides the fastest convergence for pure quadratic nonlinearity. In our calculations we chose r = 7/4 which provided the fastest convergence to the stationary solution for mixed quadratic and cubic nonlinearity.
        
        The convergence is controlled by the closeness of the parameter M to unity. Starting from the arbitrary pulse-type function υ_0(ζ), we conducted calculations with the given parameters B, G1, G2, and V on the basis of the iteration scheme (<ref>) until the parameter M was close to 1, up to small quantity ϵ, i.e., until |M - 1| ≤ϵ (in our calculations it was set to ϵ = 10^-6).
        
        To avoid a singularity in equation (<ref>), the speed of a
        solitary wave should be chosen in such a way that the
        fourth-degree polynomial in the denominator of equation (<ref>)
        does not have real roots. This corresponds to the case when there
        is no resonance between the solitary wave and a linear wave, i.e.,
        V  V_ph(κ)– see equation (<ref>). Under this
        condition only regular solitons can be constructed by means of
        this method, not the embedded soliton.
        
        
        
        § PSEUDOSPECTRAL SCHEME FOR THE FIFTH-ORDER KDV EQUATION
        
        
        The numerical scheme used for the interaction of regular and
        embedded solitons is as follows. We implement a pseudospectral
        scheme using a 4 order Runge–Kutta method for time
        stepping. The time stepping is performed in the Fourier space and
        the nonlinear terms are calculated in the real space and
        transformed back to the Fourier space for use in the method.
        Firstly we consider a solution in the domain [-L, L] and
        transform it to the domain [0, 2π] via the transform
        x̃ = S x  + π with S = π/L. Writing (<ref>)
        in the divergent form, we obtain (omitting tildes)
        
            u_t + S x [αu^2/2 + α_1 u^3/3 + β S^2 2ux + β_1 S^4 4ux + γ_1/2 S^2 2u^2x + γ_2 - 3γ_1/2 (S ux )^2 ]= 0.
        
        The terms u^2, u^3 and (ux)^2 are calculated in the real domain before transforming back to the Fourier space for use in the time-stepping algorithm.
        
        Let us discretise the solution interval by N nodes where N is a power of 2, so we have spacing Δ x = 2 π/N (in our calculation we used N = 2^13 = 8192, so that the spacial resolution was Δξ = 0.15). We use the Discrete Fourier Transform (DFT)
        
            û (k, t) = 1/√(N)∑_j=0^N-1 u (x_j, t) e^-ikx_j,    -N/2≤ k ≤N/2 - 1,
        
        where x_j = j Δ x, and k is an integer representing the discretised (and scaled) wavenumber. The inverse transform is
        
            u (x, t) = 1/√(N)∑_k=-N/2^N/2-1û (k, t) e^ikx_j,    j=0, 1, …, N-1.
        
        We make use of the Fast Fourier Transform (FFT) algorithm to implement these transforms effectively. We introduce the following notation for the last term in square brackets of equation (<ref>) to simplify the expression: ẑ = ℱℱ^-1 i k S û^2, where ℱ and ℱ^-1 represent the forward and inverse Fourier transforms respectively. Applying these transforms to (<ref>) we obtain
        
            û_t = F(û) ≡ -α i k S/2u^2 + α_1 i k S/3u^3 - β i k^3 S^3 û + β_1 i k^5 S^5 û - γ_1/2 i k^3 S^3 u^2 + γ_2 - 3γ1/2 i k S ẑ.
        
        To solve the ODE (<ref>) numerically, we use a 4 order Runge–Kutta method for time stepping. Let us assume that the solution at time t is given by û_j = û (x, j Δ t), where Δ t is the time step of integration. The solution at time t = (j+1)Δ t is given by
        
            û_j+1 = û_j + 1/6 a_j + 2 b_j + 2 c_j + d_j,
        
        where complex quantities a, b, c, and d are defined as
        
            a_j = Δ t F û_j,       b_j = Δ t F û_j + 1/2 a_j, 
               c_j =  Δ t F û_j + 1/2 b_j,       d_j = Δ t F û_j + c_j.
        
        The nonlinear terms u^2 and u^3 were evaluated in the real space and then were transformed to the Fourier space for use in equation (<ref>).
        
        Due to periodic boundary conditions (which is the intrinsic
        feature of the pseudospectral method), the radiated waves can
        re-enter the region of interest and interfere with the main wave
        structures. To alleviate this, we have introduced a damping region
        (“sponge layer”) at each end of the domain to prevent waves
        re-entering. Within the sponge layer we introduce in the left-hand
        side of equation (<ref>) a linear decay term ν 
        r(x)  u, where ν is the coefficient of artificial viscosity
        and
        
            r(x) = 1/2 2 +  tanh D  x - 3L/4 - tanh D  x + 3L/4.
        
        The coefficient D was chosen such that damping occurs only beyond the region of interest and does not affect the main wave structures. The spatially nonuniform decay term was treated numerically in the same way as the nonlinear terms.
        
        99Whitham74
        G. B. Whitham, Linear and nonlinear waves (Wiley
        Interscience Publ., John Wiley and Sons, N. Y., 1974).
        
        Karpman75
        V. I. Karpman, Nonlinear waves in dispersive media (Oxford:
        Pergamon Press, 1975).
        
        Lamb80
        G. L. Lamb, Elements of soliton theory (New York et al.:
        Johh Wiley & Sons, 1980).
        
        Ablowitz81
        M. J. Ablowitz, H. Segur, Solitons and the inverse scattering
        transform (Philadelphia: SIAM, 1981).
        
        Dodd82
        R. K. Dodd, J. C. Eilbeck, J. D. Gibbon, H. C. Morris, 
        Solitons and nonlinear wave equations (London et al.: Academic
        Press, 1982).
        
        Newell85
        A. C. Newell, Solitons in mathematics and physics
        (University of Arizona: SIAM, 1985).
        
        Ostrovsky89
        L. A. Ostrovsky, Yu. A. Stepanyants, Do internal solitons
        exist in the ocean? Rev. Geophys. 27(3), 293–310
        (1989).
        
        Apel07
        J. Apel, L. A. Ostrovsky, Y. A. Stepanyants, J. F. Lynch,
        Internal solitons in the ocean and their effect on
        underwater sound, J. Acoust. Soc. Am. 121(2), 695–722
        (2007).
        
        Ostrovsky05
        L. A. Ostrovsky, Y. A. Stepanyants, Internal solitons in
        laboratory experiments: Comparison with theoretical models, Chaos
        15, 037111 (2005).
        
        Michallet   Barthelemy98
        H. Michallet, E. Barthélemy, Experimental study of interfacial
        solitary waves, J. Fluid Mech. 366, 159–177 (1998).
        
        Benney66
        D. J. Benney, Long non-linear waves in fluid flows, J. Math. and
        Phys. 45(1–4), 52–63 (1966).
        
        Lee74
        Ch.-Y. Lee, R. C. Beardsley, The generation of long nonlinear
        internal waves in a weakly stratified shear flows, J. Geophys.
        Res. 79(3), 453–457 (1974).
        
        Koop81
        C. Koop, G. Butler, An investigation of internal solitary waves in
        a two-fluid system, J. Fluid Mech. 112, 225–251 (1981).
        
        Olver84
        P. J. Olver, Hamiltonian perturbation theory and water
        waves, Contemp. Math. 28, 231–249 (1984).
        
        MarchantSmyth90
        T. R. Marchant, N. F. Smyth, The extended Korteweg–de Vries
        equation and the resonant flow of a fluid over topography, J.
        Fluid Mech. 221, 263–288 (1990).
        
        Lamb96
        K. Lamb, The evolution of internal wave undular bores: comparisons
        of a fully nonlinear numerical model with weakly nonlinear theory,
        J. Phys. Oceanography 26, 2712–2734 (1996).
        
        Grimshaw02
        R. Grimshaw, E. Pelinovsky, O. Poloukhina, Higher-order
        Korteweg–de Vries models for internal solitary waves in a
        stratified shear flow with a free surface, Nonlin. Proc. Geophys.
        9, 221–235 (2002).
        
        Giniyatullin14
        A. R. Giniyatullin, A. A. Kurkin, O. E. Kurkina,
        Y. A. Stepanyants, Generalised Korteweg–de Vries equation for
        internal waves in two-layer fluid, Fundamental and Applied
        Hydrophysics (in Russian) 7(4), 16–28 (2014).
        
        KarRozRut14
        A. Karczewska, P. Rozmej, Ł. Rutkowski, A new nonlinear
        equation in the shallow water wave problem, Phys. Scr.
        89, 054026 (2014).
        
        Karczewska14
        A. Karczewska, P. Rozmej, E. Infeld, Shallow-water soliton
        dynamics beyond the Korteweg–de Vries equation, Phys. Rev. E
        90, 012907 (2014).
        
        Kakutani69
        T. Kakutani, H. Ono,  Weak non-linear hydromagnetic waves in a
        cold collisionless plasma, J. Phys. Soc. Japan  26,
        1305–1318 (1969).
        
        Gorshkov76
        K. A. Gorshkov, L. A. Ostrovsky, V. V. Papko, Interactions
        and bound states of solitons as classical particles, Sov. Phys.
        JETP 44(2), 306–311 (1976).
        
        Abramyan85
        L. A. Abramyan, Yu. A. Stepanyants,  The structure of
        two-dimensional solitons in media with anomalously small
        dispersion, Sov. Phys. JETP  61(5), 963–966 (1985).
        
        Hunter88
        J. K. Hunter, J. Scheurle,  Existence of perturbed solitary
        wave solutions to a model equation for water waves, Physica D
         32 253–268, (1988).
        
        Guyenne14
        P. Guyenne, E. I. Părău, Finite-depth effects on
        solitary waves in a floating ice sheet, J. Fluids Struct.
        49, 242–262 (2014).
        
        Slyunyaev99
        A. Slyunyaev, E. Pelinovskii, Dynamics of large-amplitude
        solitons, JETP 89, 173–181 (1999).
        
        Slyunyaev01
        A. V. Slyunyaev, Dynamics of localized waves with large
        amplitude in a weakly dispersive medium with quadratic and
        positive cubic nonlinearity, JETP 92, 529–534 (2001).
        
        Kichenassamy92
        S. Kichenassamy, P. J. Olver, Existence and nonexistence
        of solitary wave solutions to higher-order model evolution
        equations, SIAM J. Math. Anal. 23(5), 1141–1166,
        (1992).
        
        Champneys97
        A. R. Champneys, M. D. Groves, A global investigation of
        solitary-wave solutions to a two-parameter model for water waves,
        J. Fluid Mech. 342, 199–229 (1997).
        
        Yang01
        J. Yang, Dynamics of embedded solitons in the extended KdV
        equations, Stud. Appl. Math. 106, 337–365 (2001).
        
        Kodama85
        Y. Kodama, Normal forms for weakly dispersive wave equations,
        Phys. Lett. A 112, 193–196 (1985).
        
        Fokas95
        A. S. Fokas, On a class of physically important integrable
        equations, Physica D 87, 145–150 (1995).
        
        FokasLiu96
        A. Fokas, Q. M. Liu, Asymptotic integrability of water waves,
        Phys. Rev. Lett. 77, 2347–2357 (1996).
        
        LiZhiSib97
        Li Zhi, N. R. Sibgatullin, An improved theory of long
        waves on the water surface, J. Appl. Math. Mech. 61(2),
        177–182 (1997).
        
        Kraenkel98
        R. A. Kraenkel, First-order perturbed Korteweg–de Vries solitons,
        Phys. Rev. E. 57(4), 4775–4777 (1998).
        
        MarchantSmyth96
        T. R. Marchant, N. F. Smyth, Soliton interaction for the extended
        Korteweg–de Vries equation, IMA J. Appl. Math. 56, 157–176 (1996).
        
        Marchant99
        T. R. Marchant, Asymptotic solitons of the extended Korteweg–de
        Vries equation, Phys. Rev. E 59, 3745–3748 (1999).
        
        MarchantSmyth2006
        T. R. Marchant, N. F. Smyth, An undular bore solution for the
        higher-order Korteweg–de Vries equation, J. Phys. A: Math. Gen.
        39, L563–L569 (2006).
        
        Kawahara72
        T. Kawahara, Oscillatory solitary waves in dispersive media, J.
        Phys. Soc. Japan 33, 260–264 (1972).
        
        Yamamoto81
        Y. Yamamoto, E. I. Takizawa, On a solution of nonlinear
        time-evolution equation of fifth order, J. Phys. Soc. Japan
        50, 1421–1422 (1981).
        
        Gorshkov79
        K. A. Gorshkov, L. A. Ostrovsky, V. V. Papko, A. S. Pikovsky,
        On the existence of stationary multisolitons, Phys. Lett.
        A 73(3–4), 177-179 (1979).
        
        Kurkina15
        O. Kurkina, N. Singh, Y. Stepanyants, Structure of
        internal solitary waves in two-layer fluid at near-critical
        situation, Comm. Nonlin. Sci. Num. Simulation 22(5),
        1235–1242 (2015).
        
        YMKC01
        J. Yang, B. A. Malomed, D. J. Kaup, A. R. Champneys,
        Embedded solitons: a new type of solitary wave, Math.
        Comp. Simulation  56, 585-600 (2001).
        
        Yang99
        J. Yang, B. A. Malomed, D. J. Kaup, Embedded solitons in
        second-harmonic generating systems, Phys. Rev. Lett. 83,
        1958–1961 (1999).
        
        Nishikawa74
        K. Nishikawa, H. Hojo, K. Mima, H. Ikezi, Coupled
        nonlinear electron-plasma and ion-acoustic waves, Phys. Rev.
        Lett. 33(3), 148–151 (1974).
        
        Gorshkov83
        K. A. Gorshkov, V. A. Mironov, A. M. Sergeev, Coupled
        stationary soliton formations, In: “Nonlinear Waves.
        Self-Organization”, ed. A. V. Gaponov-Grekov and
        M. I. Rabinovich, Nauka, Moscow, 112–128 (1983) (in Russian).
        
        Petviashvili92
        V. I. Petviashvili, O. V. Pokhotelov, Solitary waves in
        plasmas and in the atmosphere, (Gordon and Breach, Philadelphia,
        1992).
        
        Yang10
        J. Yang, Nonlinear waves in integrable and nonintegrable
        systems, SIAM, Philadelphia (2010).
        
        Pelinovsky04
        D. E. Pelinovsky, Y. A. Stepanyants, Convergence of Petviashvili's
        iteration method for numerical approximation of stationary
        solutions of nonlinear wave equations, SIAM J. Numerical Analysis
        42(3), 1110–1127 (2004).
        
        Grimshaw08
        R. Grimshaw, K. Helfrich, Long-time solutions of the
        Ostrovsky equation, Stud. Appl. Math. 121, 71–88 (2008).
        
        Alias13
        A. Alias, R. H. J. Grimshaw, K. R. Khusnutdinova, On strongly
        interacting internal waves in a rotating ocean and coupled Ostrovsky
        equations, Chaos 23, 023121 (2013).
        
        Alias14
        A. Alias, R. H. J. Grimshaw, K. R. Khusnutdinova, Coupled
        Ostrovsky equations for internal waves in a shear flow, Physics
        of Fluids 26, 126603 (2014).
        
        Khusnutdinova17
        K. R. Khusnutdinova, M. R. Tranter, On radiating solitary
        waves in bi–layers with delamination and coupled Ostrovsky
        equations, Chaos 27, 013112 (2017).
        
        Yang08
        J. Yang, T. I. Lakoba, Accelerated imaginary-time
        evolution methods for the computation of solitary waves, Studies
        in Appl. Maths. 120, 265–292 (2008).
        
        OstrEtAl15
        L. A. Ostrovsky, E. N. Pelinovsky, V. I. Shrira,
        Y. A. Stepanyants, Beyond the KDV: Post-explosion development,
        Chaos 25(9), 097620 (2015).
        
        Boyd01
        J. P. Boyd, Chebyshev and Fourier Spectral Methods,
        (Dover, Mineola, New York, 2001).
        
        Obregon98
        M. A. Obregon, Yu. A. Stepanyants, Oblique
        magneto-acoustic solitons in rotating plasma, Phys. Lett. A
        249(4), 315–323 (1998).
        
        Talipova99
        T. G. Talipova, E. N. Pelinovsky, K. Lamb, R. Grimshaw, P.
        Holloway, Cubic nonlinearity effects in the propagation of
        intense internal waves, Doklady Academii Nauk 365(6),
        824–827 (1999) (in Russian). (Engl. transl.: Doklady Earth Sci.
        365(2), 241–244 (1999).)
        
        RuvKurk14
        E. A. Ruvinskaya, O. E. Kurkina, A. A Kurkin, Dynamics of
        nonlinear internal gravitational waves in stratified fluids
        (Nizhny Novgorod State Technical Universitry, Nizhny Novgorod,
        2014, in Russian).
        
        Grimshaw17
        R. H. J. Grimshaw, K. R. Khusnutdinova, K. R. Moore,
        Radiating solitary waves in coupled Boussinesq equations,
        IMA J. Appl. Math.  82(4) 802–820 (2017).
        
        J
        R.S. Johnson, Ring waves on the surface of shear flows: a
        linear and nonlinear theory, J. Fluid Mech. 215 145–160
        (1990).
        
        L
        V.D. Lipovskii, On the nonlinear internal wave theory in
        fluid of finite depth, Izv. Akad. Nauk SSSR, Ser. Fiz. Atm.
        Okeana 21 864–871 (1985).
        
        KZ_JFM
        K.R. Khusnutdinova, X. Zhang, Long ring waves in a
        stratified fluid over a shear flow, J. Fluid Mech. 794
        17–44 (2016).
        
        KZ_PhysD
        K.R. Khusnutdinova, X. Zhang, Nonlinear ring waves in a
        two-layer fluid, Physica D 333 208–221 (2016).
