Two-pulse solutions in the fifth-order KdV equation : rigorous theory and numerical approximations

We revisit existence and stability of two-pulse solutions in the fifth-order Korteweg--de Vries (KdV) equation with two new results. First, we modify the Petviashvili method of successive iterations for numerical (spectral) approximations of pulses and prove convergence of iterations in a neighborhood of two-pulse solutions. Second, we prove structural stability of embedded eigenvalues of negative Krein signature in a linearized KdV equation. Combined with stability analysis in Pontryagin spaces, this result completes the proof of spectral stability of the corresponding two-pulse solutions. Eigenvalues of the linearized problem are approximated numerically in exponentially weighted spaces where embedded eigenvalues are isolated from the continuous spectrum. Approximations of eigenvalues and full numerical simulations of the fifth-order KdV equation confirm stability of two-pulse solutions related to the minima of the effective interaction potential and instability of two-pulse solutions related to the maxima points.


Introduction
One-pulse solutions (solitons) are commonly met in many nonlinear evolution equations where dispersive terms (represented by unbounded differential operators) and nonlinear terms (represented by power functions) are taken in a certain balance. Typical examples of such nonlinear evolution equations with one-pulse solutions are given by the NLS (nonlinear Schrödinger) equation, the Klein-Gordon (nonlinear wave) equation and the KdV (Korteweg-de Vries) equation, as well as their countless generalizations.
One-pulse solutions are the only stationary (traveling) localized solutions of the simplest nonlinear evolution equations. However, uniqueness is not a generic property and bound states of spatially separated pulses can represent other stationary (traveling) localized solutions of the same evolution equation. For instance, two-pulse, three-pulse, and generally N -pulse solutions exist in nonlinear evolution equations with a higher-order dispersion (represented by a higher-order differential operator). The prototypical example of such situation is the fifth-order KdV equation in the form, x ∈ R, t ∈ R + , (1.1) where u : R × R + → R and all coefficients of the nonlinear PDE are normalized by a scaling transformation. See Bridges & Derks [BD02] for a review of history and applications of the fifth-order KdV equation (1.1) to magneto-acoustic waves in plasma and capillary-gravity water waves.
Traveling localized solutions u(x, t) = φ(x − ct) of the fifth-order KdV equation (1.1) satisfies the fourth-order ODE where z = x − ct is the traveling coordinate and one integration of the fifth-order ODE in z is performed subject to zero boundary conditions on φ(z) and its derivatives as |z| → ∞. Existence of localized solutions (homoclinic orbits) to the fourth-order ODE (1.2) was considered by methods of the dynamical system theory. See Champneys [C98] for a review of various results on existence of homoclinic orbits in the ODE (1.2).
In particular, it is proved with the variational method by Buffoni & Sere [BS96] and Groves [G98] (see references to earlier works in [C98]) that the fourth-order ODE (1.2) has the one-pulse solution φ(z) for c > 0, which is the only localized solution of the ODE (1.2) for 0 < c < 1 4 up to the translation φ(z − s) for any s ∈ R. The analytical expression for the one-pulse solution is only available for c = 36 169 < 1 4 with φ(z) = 105 338 sech 4 z 2 √ 13 . (1.3) For c > 1 4 , the fourth-order ODE (1.2) has infinitely many multi-pulse solutions in addition to the onepulse solution [BS96,G98]. The multi-pulse solutions look like multiple copies of the one-pulse solutions separated by finitely many oscillations close to the zero equilibrium φ = 0. Stability and evolution of multi-pulse solutions are beyond the framework of the fourth-order ODE (1.2) and these questions were considered by two equivalent theories in the recent past.
The pioneer work of Gorshkov & Ostrovsky explains multi-pulse solutions of the fifth-order KdV equation (1.1) from the effective interaction potential computed from the one-pulse solution [GO81,GOP84]. When the interaction potential has an alternating sequence of maxima and minima (which corresponds to the case when the one-pulse solution φ(z) has oscillatory decaying tails at infinity), an infinite countable sequence of two-pulse solutions emerge with the distance between the pulses near the sequence of extremal points. Three-pulse solutions can be constructed as a bi-infinite countable sequence of three one-pulse solutions where each pair of two adjacent pulses is located approximately at a distance defined by the two-pulse solution. Similarly, N -pulse solutions can be formed by a (N − 1)infinite countable sequence of N copies of one-pulse solutions. The perturbative procedure in [GO81] has the advantages that both the linear and nonlinear stability of multi-pulse solutions can be predicted from analysis of the approximate ODE system derived for distances between the individual pulses. Numerical evidences of validity of this procedure in the context of the fifth-order KdV equation are reported in [BC97].
problem for the fifth-order KdV equation takes the form ∂ z Hv = λv, z ∈ R, (1.5) where v : R → C is an eigenfunction for a small perturbation of φ(z) in the reference frame z = x − ct and λ ∈ C is an eigenvalue. We say that the eigenvalue λ is unstable if Re(λ) > 0 and v ∈ H 2 (R). We say that the eigenvalue λ is of negative Krein signature if Re(λ) = 0, Im(λ) > 0, v ∈ H 2 (R) and (Hv, v) < 0.
Our interest to this well-studied problem is revived by the recent progress in the spectral theory of nonself-adjoint operators arising from linearizations of nonlinear evolution equations [CP06]. These operators can be defined as self-adjoint operators into Pontryagin space where they have a finite-dimensional negative invariant subspace. Two physically relevant problems for the fifth-order KdV equation (1.1) have been solved recently by using the formalism of operators in Pontryagin spaces. First, convergence of the numerical iteration method (known as the Petviashvili's method) for one-pulse solutions of the ODE (1.2) was proved by the contraction mapping principle in a weighted Hilbert space (which is equivalent to Pontryagin space with zero index) [PS04]. Second, eigenvalues of the spectral stability problem in a linearization of the fifth-order KdV equation (1.1) were characterized in Pontryagin space with a non-zero index defined by the finite number of negative eigenvalues of H [KP05,CP06].
Both recent works rise some open problems when the methods are applied to the N -pulse solutions in the fifth-order KdV equation (1.1), even in the case of two-pulse solutions (N = 2). The successive iterations of the Petviashvili's method do not converge for two-pulse solutions. The iterative sequence with two pulses leads either to a single pulse or to a spurious solution with two pulses located at an arbitrary distance (see Remark 6.5 in [PS04]). This numerical problem arises due to the presence of small and negative eigenvalues of H. A modification of the Petviashvili's method is needed to suppress these eigenvalues similarly to the work of Demanet & Schlag [DS06] where the zero eigenvalue associated to the translational invariance of the three-dimensional NLS equation is suppressed. We shall present the modification of the iterative Petviashvili's method in this article.
Our numerical method complements the recent work of Yang & Lakoba [YL06] where the steepest descent method is developed to suppress any unstable eigenvalues of the iterative method and to enforce its convergence. We note that other techniques for numerical approximation of multi-pulse solutions in the fifth-order KdV equation have been developed by Champneys and coworkers [CS93,OCK03] on the basis of a delicate combination of the numerical shooting method and the continuation techniques.
Open questions also arise when spectral stability of multi-pulse solutions is considered within the linear eigenvalue problem (1.5). By either the Gorshkov-Ostrovsky perturbative procedure or the Sandstede-Lin reduction method, the small eigenvalues of the Jacobian operator H result in small eigenvalues of the linearized operator ∂ z H, which are either pairs of real eigenvalues (one of which is unstable) or pairs of purely imaginary eigenvalues of negative Krein signature (which are neutrally stable but potentially unstable). Both cases are possible in the fifth-order KdV equation in agreement with the count of unstable eigenvalues in Pontryagin spaces (see Theorem 6 in [CP06]). (Similar count of unstable eigenvalues and eigenvalues of negative Krein signatures was developed for the NLS equations in recent papers [KKS04,P05].) Since the real eigenvalues are isolated from the continuous spectrum of the eigenvalue problem (1.5), they are structurally stable and persist with respect to parameter continuations. However, the purely imaginary eigenvalues are embedded into the continuous spectrum of the eigenvalue problem (1.5) and their destiny remains unclear within the reduction methods. It is well known for the NLS-type and Klein-Gordon-type equations that embedded eigenvalues are structurally unstable to the parameter continuations [G90]. If a certain Fermi golden rule related to the perturbation term is nonzero, the embedded eigenvalues of negative Krein signature bifurcate off the imaginary axis to complex eigenvalues inducing instabilities of pulse solutions [CPV05]. (The embedded eigenvalues of positive Krein signature simply disappear upon a generic perturbation [CPV05].) This bifurcation does not contradict the count of unstable eigenvalues [KKS04,P05] and it is indeed observed in numerical approximations of various pulse solutions of the coupled NLS equations [PY05].
From a physical point of view, we would expect that the time evolution of an energetically stable superposition of stable one-pulse solutions remains stable. (Stability of one-pulse solutions in the fifthorder KdV equation (1.1) was established with the variational theory [L99] and the multi-symplectic Evans function method [BD02,BDG02].) According to the Gorshkov-Ostrovsky perturbative procedure, dynamics of well-separated pulses is represented by the Newton particle's law which describes nonlinear stability of oscillations near the minima of the effective interaction potential [GOP84]. Therefore, we would rather expect (on the contrary to embedded eigenvalues in the linearized NLS and Klein-Gordon equations) that the embedded eigenvalues of negative Krein signature are structurally stable in the linear eigenvalue problem (1.5) and persist beyond the leading order of the perturbative procedure. (Multipulse solutions of the NLS and Klein-Gordon equations with well-separated individual pulses are always linearly stable since the small purely imaginary eigenvalues of the Lyapunov-Schmidt reductions are isolated from the continuous spectrum of the corresponding linearized problems [YSJ00].) Since the count of unstable eigenvalues in [CP06] does not allow us to prove structural stability of embedded eigenvalues of negative Krein signature, we address this problem separately by using different analytical and numerical techniques. In particular, we present an analytical proof of persistence (structural stability) of embedded eigenvalues of negative Krein signature in the linearized problem (1.5). We also apply the Fourier spectral method and illustrate the linearized stability of the corresponding two-pulse solutions numerically. Our analytical and numerical methods are based on the construction of exponentially weighted spaces for the linear eigenvalue problem (1.5). (See [PW94] for analysis of exponentially weighted spaces in the context of the generalized KdV equation.) Our results complement the recent work of Sandstede [S06] and Bridges et al. [BCD06], where two-pulse solutions of the fifthorder KdV equation (1.1) were also addressed. Additionally, Lyapunov-Schmidt reductions of two-pulse solutions of the coupled KdV equations were recently considered by Scheel & Wright [SW06].
This article is structured as follows. Section 2 contains summary of available results on existence and stability of one-pulse and two-pulse solutions of the fifth-order KdV equation (1.1). Section 3 presents a modification of the iterative Petviashvili method for convergent numerical approximations of the twopulse solutions in the fourth-order ODE (1.2). Section 4 develops the proof of structural stability of embedded eigenvalues in the eigenvalue problem (1.5) and numerical approximations of unstable and stable eigenvalues in an exponentially weighted space. Section 5 describes full numerical simulations of the fifth-order KdV equation (1.1) to study nonlinear dynamics of two-pulse solutions.
2 Review of available results on two-pulse solutions When c < 0, one pair of roots κ is purely imaginary and the other pair is purely real. When 0 < c < 1 4 , two pairs of roots κ are real-valued. When c > 1 4 , the four complex-valued roots κ are located symmetric about the axes. We will use notations k 0 = Im(κ) > 0 and κ 0 = Re(κ) > 0 for a complex root of (2.1) in the first quadrant for c > 1 4 . The following two theorems summarize known results on existence of one-pulse and two-pulse solutions of the ODE (1.2).
(ii) The Jacobian operator H associated with the one-pulse solution φ(z) has exactly one negative eigenvalue with an even eigenfunction and a simple kernel with the odd eigenfunction φ ′ (z).
Proof. (i) Existence of a symmetric solution φ(z) in H 2 (R) follows by the mountain-pass lemma and concentration-compactness principle (see Theorem 8 in [G98] and Theorem 2.3 in [L99]). The equivalence between weak solutions of the variational theory and strong solutions of the ODE (1.2) is established in Lemma 1 of [G98] and Lemma 2.4 of [L99]. The exponential decay of φ(z) follows by the Stable Manifold Theorem in Appendix A of [BS96]. Global uniqueness of the symmetric solution φ(z) for 0 < c < 1 4 (module the space translation in z) is proved in [AT92]. Local uniqueness of the symmetric one-pulse solution φ(z) for c ≥ 1 4 follows from the fact that the Jacobian operator H has a simple kernel due to the translational invariance (statement (ii)). Finally, the smoothness of the function φ(z) is proved from the ODE (1.2) by the bootstrapping principle [CCP06].
(ii) The Jacobian operator H coincides with the Hessian of the energy functional under a fixed value of the momentum functional [BSS87]. Uniqueness of the negative eigenvalue of H follows from the fact that the one-pulse solution φ(z) is a ground state (a global minimizer) of the constrained manifold (see Proposition 16 in [G98]). The ground state for operators with even derivatives and symmetric potential functions is always given by the even eigenfunction. The kernel of H is simple due to the duality principle in Theorem 4.1 of [BS96]. If it is not simple, then global two-dimensional stable and unstable manifolds coincide and the time for a homoclinic orbit to go from the local unstable manifold to the local stable manifold is uniformly bounded. However, it was shown in [B95] that a sequence of homoclinic solutions {u n } n∈N exists such that the time between local manifolds grows linearly in n. Hence, the kernel of H is simple with the odd eigenfunction φ ′ (z) due to the space translation.
(iii) Smoothness of the map φ(z) from c > 0 to H 2 (R) is a standard assumption (see Assumption 5.1 in [L99]). The positive value of P ′ (c) and the absence of eigenvalues of ∂ z H with Re(λ) > 0 follow from the stability of one-pulse solutions in Main Theorem of [BSS87], Theorems 4.1 and 5.3 of [L99] and Theorem 8.1 of [BD02]. Indeed, the only negative eigenvalue of H becomes a positive eigenvalue in the constrained subspace where the corresponding eigenfunction is orthogonal to φ(z) [SS90]. The two-dimensional algebraic kernel of ∂ z H follows from the derivatives of the ODE (1.2) in z and c: The algebraic kernel of ∂ z H is exactly two-dimensional under the condition P ′ (c) = 0 [PW92].
Theorem 2.2 (Two-pulse solutions) Consider the ODE (1.2) with c > 1 4 . There exists an infinite countable set of two-pulse solutions φ(z) such that φ ∈ H 2 (R) ∩ C 5 (R), φ(−z) = φ(z), φ(z) → 0 exponentially as |x| → ∞, and φ(z) resembles two copies of the one-pulse solutions described in Theorem 2.1 which are separated by small-amplitude oscillatory tails. The members of the set are distinguished by the distance L between individual pulses which takes the discrete values {L n } n∈N . Moreover, for any small δ > 0 there exists γ > 0 such that Proof. Existence of an infinite sequence of geometrically distinct two-pulse solutions with the distances distributed by (2.3) follows by the variational theory in Theorem 1.1 of [BS96] under the assumption that the single-pulse solution φ(z) is isolated (up to the space translations). This assumption is satisfied by Theorem 2.1(ii).
The following theorem describes an asymptotic construction of the two-pulse solutions, which is used in the rest of our paper.
Theorem 2.3 (Lyapunov-Schmidt reductions for two-pulse solutions) Let c > 1 4 and let Φ(z) denote the one-pulse solution described by Theorem 2.1. Let L = 2s be the distance between two copies of the one-pulse solutions of the ODE (1.2) in the decomposition where ϕ(x) is a remainder term. Let W (L) be C m (R + ) function for any m ≥ 0 defined by There exists an infinite countable set of extrema of W (L), which is denoted by {L n } n∈N .
(i) Assume that W ′′ (L n ) = 0 for a given n ∈ N. There exists a unique symmetric two-pulse solution φ(z) ≡ φ n (z) described by Theorem 2.2, such that the distance L between individual pulses satisfies the bound |L − L n | ≤ C n e −κ 0 L (2.6) for some C n > 0.
(ii) The Jacobian H associated with the n-th two-pulse solution φ n (z) has exactly two finite negative eigenvalues with even and odd eigenfunctions, a simple kernel with the odd eigenfunction φ ′ n (z) and a small eigenvalue µ with an even eigenfunction, such that In particular, the small eigenvalue µ is negative when W ′′ (L n ) > 0 and positive when W ′′ (L n ) < 0.
(iii) There exists a pair of small eigenvalues λ of the linearized operator ∂ z H associated with the n-th two-pulse solution φ n (z), such that In particular, the pair is real when W ′′ (L n ) < 0 and purely imaginary (up to the leading order) with negative Krein signature when W ′′ (L n ) > 0.
Proof. When the tails of the one-pulse solution Φ(z) are decaying and oscillatory (i.e. when c > 1 4 ), the smooth function W (L) in (2.5) is decaying and oscillatory in L and an infinite set of extrema {L n } n∈N exists. Let us pick L n for a fixed value of n ∈ N such that W ′ (L n ) = 0 and W ′′ (L n ) = 0.
(i) When the decomposition (2.4) is substituted into the ODE (1.2), one can find the ODE for ϕ(z): (2.9) Let ǫ = O(e −κ 0 L ) be a small parameter that measures the L ∞ -norm of the overlapping term Φ(z − s)Φ(z + s) in the sense that for each ǫ > 0 there exist constants C 0 > 0 and s 0 > 0 such that (2.10) Denote L = 2s and ǫΨ(z, L) = 2Φ(z)Φ(z + L) and rewrite the ODE (2.9) forφ(z) = ϕ(z + s): The ODE (2.11) is defined in function space H 2 (R), where the Jacobian for the one-pulse solution has a simple kernel with the odd eigenfunction Φ ′ (z) by Theorem 2.1(ii). By the Lyapunov-Schmidt reduction method (see [GS85]), there exists a unique solutionφ =φ(z, where the first term has the order of O(ǫ) and the second term has the order of O(ǫ 2 ). The statement follows by the Implicit Function Theorem applied to the scalar equation F (L, ǫ) = 0 under the assumption that the root L n of W ′ (L) is simple.
(ii) The Jacobian H associated with the n-th two-pulse solution φ(z) ≡ φ n (z) given by the decomposition (2.4) has the form: In the limit of ǫ = 0, where ǫ = O(e −κ 0 L ), the Jacobian H has a double negative eigenvalue and a double zero eigenvalue. By a linear combination of eigenfunctions, one can construct one even and one odd eigenfunctions for each of the double eigenvalues. By continuity of eigenvalues of self-adjoint operators, the double negative eigenvalue splits but the two negative eigenvalues remain negative if ǫ is small. By reversibility of the system, eigenfunctions for simple eigenvalues are either even or odd and by continuity of eigenfunctions, there is exactly one even and one odd eigenfunctions for the two negative eigenvalues. By the translation invariance, the double zero eigenvalue splits into a simple zero eigenvalue which corresponds to the odd eigenfunction φ ′ n (z) and to a small non-zero eigenvalue of the order of O(ǫ) which corresponds to an even eigenfunction. The splitting of the double zero eigenvalue in the problem Hv = µv is considered by the perturbation theory, where (α 1 , α 2 ) are coordinates of the projections to the kernel of H at ǫ = 0 and V (z) is the remainder term of the order of O(ǫ). By projecting the eigenvalue problem Hv = µv to the kernel of H and neglecting the terms of the order of O(ǫ 2 ), one can obtain a reduced eigenvalue problem: Since one eigenvalue must be zero with the odd eigenfunction φ ′ n (z), the zero eigenvalue corresponds to the eigenfunction (2.13) with α 1 = α 2 up to the leading order. By looking at the linear system, we find that the zero eigenvalue corresponding to α 1 = α 2 exists only ifW = W ′′ (L n ). The other eigenvalue at the leading order is µ = −2W ′′ (L n )/Q(c). The non-zero eigenvalue of the reduced eigenvalue problem corresponds to the even eigenfunction (2.13) with α 1 = −α 2 . The eigenvalue µ has the order of O(ǫ). By continuity of isolated eigenvalues H with respect to perturbation terms of the Lyapunov-Schmidt reductions, we obtain the result (2.7).
(iii) In the limit of ǫ = 0, where ǫ = O(e −κ 0 L ), the linearized operator ∂ z H for the n-th two-pulse solution φ n (z) has a four-dimensional algebraic kernel according to the two-dimensional kernel of the one-pulse solution (2.2). By the translation invariance, the two-dimensional algebraic kernel survives for small ǫ with the eigenfunctions {φ ′ n (z), ∂ c φ n (z)}. Two eigenvalues λ of the operator ∂ z H may bifurcate from the zero eigenvalue at the order of O( √ ǫ). The splitting of the zero eigenvalue in the problem where (α 1 , α 2 , β 1 , β 2 ) are coordinates of the projections to the algebraic kernel of ∂ z H at ǫ = 0 and V (z) is the remainder term of the order of O(ǫ). By projecting the eigenvalue problem ∂ z Hv = λv to the algebraic kernel of the adjoint operator −H∂ z and neglecting the terms of the order of O(λǫ, ǫ 2 ), we immediately find at the order of O(λ) that β j = λα j , j = 1, 2, while (α 1 , α 2 ) satisfy a reduced eigenvalue problem at the order of O(λ 2 , ǫ): where P (c) = Φ 2 L 2 andW = W ′′ (L n ). The non-zero eigenvalue λ 2 at the leading order is Isolated eigenvalues ∂ z H are continuous with respect to perturbation terms of the Lyapunov-Schmidt reductions, so that we immediately obtain the result (2.8) for λ ∈ R when W ′′ (L n ) < 0. In order to prove (2.8) for λ ∈ iR when W ′′ (L n ) > 0, we compute asymptotically the energy quadratic form where v(z) is given by the eigenfunction (2.14) with α 1 = −α 2 = 1 and β j = λα j , j = 1, 2. When λ ∈ iR and W ′′ (L n ) > 0, we have (Hv, v) < 0 up to the leading order, such that λ is an eigenvalue of negative Krein signature. Continuity of the eigenvalues of negative Krein signature (even although the eigenvalues λ ∈ iR are embedded into the continuous spectrum of ∂ z H) follows from Pontryagin's Invariant Subspace Theorem (Theorem 1 in [CP06]). By using the exponentially weighted spaces [PW94], continuity of eigenvalues λ ∈ iR holds and the bound (2.8) is obtained for an eigenvalue λ ∈ C in a neighborhood of the point λ ∈ iR. Note that the continuity of one pair of eigenvalues λ ∈ C satisfying (2.8) does not exclude the possibility of existence of another pair of eigenvalues λ ∈ C (in the exponentially weighted space with a negative weight) which also satisfies (2.8).
Remark 2.4 Theorem 2.3 is a modification of Theorems 1 and 2 in [S98] (see also [L90]) in applications to the ODE (1.2) and the spectral problem (1.5). The proof is lighten compared to the original proofs in [L90, S98]. We note that the persistence of purely imaginary eigenvalues (2.8) in the full problem (1.5) cannot be proved with the Lyapunov-Schmidt reduction method since the essential spectrum of ∂ z H occurs on the imaginary axis (contrary to the standard assumption of Theorem 2 in [S98] that the essential spectrum occurs in the left half-plane.) Remark 2.5 The following claim from the Gorshkov-Ostrovsky perturbative procedure [GO81,GOP84] illustrates the role of W (L) as the effective interaction potential for the slow dynamics of a two-pulse solution: Claim: Let C 1 , C 2 be some positive constants. For the initial time interval 0 ≤ t ≤ C 1 e κ 0 L/2 and up to the leading order O(e −κ 0 L ), the two-pulse solutions of the fifth-order KdV equation (1.1) can be written as the decomposition where U L ∞ ≤ C 2 e −κ 0 L and the slow dynamics of L(t) = 2s(t) is represented by the Newton's particle law: Although rigorous bounds on the time interval and the truncation error of the Newton's particle law were recently found in the context of NLS solitons in external potentials (see [JFGS06] and references therein), the above claim was not proved yet in the context of two-pulse solutions of the fifth-order KdV equation (1.1). We note that perturbation analysis that leads to the Newton's particle law (2.15) cannot be used to claim persistence and topological equivalence of dynamics of the Newton's particle to the full dynamics of two-pulse solutions.
According to Theorem 2.3, an infinite set of extrema of W (L) generates a sequence of equilibrium configurations for the two-pulse solutions in Theorem 2.2. Since P ′ (c) > 0 by Theorem 2.1(iii), the maxima points of W (L) correspond to a pair of real eigenvalues λ of the spectral problem (1.5), while the minima points of W (L) correspond to a pair of purely imaginary eigenvalues λ. The two-pulse solutions at the maxima points are thus expected to be linearly and nonlinearly unstable. The twopulse solutions at the minima points are stable within the leading-order approximation (2.8) and within the Newton's particle law (2.15) (a particle with the coordinate L(t) performs a periodic oscillation in the potential well). Correspondence of these predictions to the original PDE (1.1) is a subject of the present article. We will compute the interaction potential W (L) and the sequence of its extrema points {L n } n∈N , as well as the numerical approximations of the two-pulse solutions of the ODE (1.2) and of the eigenvalues of the operator ∂ z H in (1.5).

Iterations of the Petviashvili's method for two-pulse solutions
We address the Petviashvili's iteration method for numerical approximations of solutions of the fourthorder ODE (1.2) with c > 0. See review of literature on the Petviashvili's method in [PS04]. By using the standard Fourier transform in L 1 (R) ∩ L 2 (R) we reformulate the ODE (1.2) as a fixed-point problem in the Sobolev space H 2 (R): where φ 2 (k) can be represented by the convolution integral ofφ(k) to itself. An even real-valued solution φ(−z) = φ(z) of the ODE (1.2) in H 2 (R) is equivalent to the even real-valued solutionφ(−k) =φ(k) of the fixed-point problem (3.1). Let us denote the space of all even functions in H 2 (R) by H 2 ev (R) and consider solutions of the fixed-point problem (3.1) in H 2 ev (R). Let {û n (k)} ∞ n=0 be a sequence of Fourier transforms in H 2 ev (R) defined recursively bŷ whereû 0 (k) ∈ H 2 ev (R) is a starting approximation and M n ≡ M [û n ] is the Petviashvili's stabilizing factor defined by (3. 3) It follows from the fixed-point problem (3.1) that M [φ] = 1 for any solutionφ ∈ H 2 ev (R). The following theorem was proved in [PS04] and reviewed in [DS06]. Proof. We review the basic steps of the proof, which is based on the contraction mapping principle in a local neighborhood ofφ in H 2 ev (R). The linearization of the iteration map (3.2) at the solution φ is rewritten in the physical space z ∈ R as follows: where α n is a projection of v n onto φ 2 in L 2 (R): such that u n = φ + v n and M n = 1 − α n to the linear order. The operator T = (c − ∂ 2 z + ∂ 4 z ) −1 H is a self-adjoint operator in Pontryagin space Π 0 defined by the inner product See [CP06] for review of Pontryagin spaces and the Pontryagin Invariant Subspace Theorem. Since c > 0, the Pontryagin space Π 0 has zero index and, by the Pontryagin Theorem, the operator T in Π 0 has exactly one negative eigenvalue, a simple kernel and infinitely many positive eigenvalues. The eigenfunctions for the negative and zero eigenvalues are known exactly as Due to orthogonality of the eigenfunctions in the Pontryagin space Π 0 and the relation we observe that α n is a projection of v n to φ in Π 0 , which satisfies the trivial iteration map: Projection of v n to φ ′ in Π 0 is zero since v n ∈ H 2 ev (R). As a result, the linearized iteration map (3.5) defines a contraction map if the maximal positive eigenvalue of T in L 2 (R) is smaller than 2. However, If φ(z) ≥ 0 on z ∈ R, the right-hand-side of (3.6) is zero. Otherwise, the right-hand-side of (3.6) is bounded from above by 2 c |inf z∈R φ(z)|, which leads to the condition (3.4). (ii) N -the number of terms in the partial sum for the truncated Fourier series such that the grid size h of the discretization is h = 2d/N ; (iii) ε -the small tolerance distance that measures deviation of M n from 1 and the distance between two successive approximations, such that the method can be terminated at the iteration n if andφ = u n (z) can be taken as the numerical approximation of the solution φ(z).
The numerical approximation depends weakly of the three numerical parameters, provided (i) d is much larger than the half-width of the one-pulse solution, (ii) N is sufficiently large for convergence of the Fourier series, and (iii) ε is sufficiently small above the level of the round-off error. Indeed, the constraint (i) ensures that the truncation error is exponentially small when the one-pulse solution is replaced by the periodic sequence of one-pulse solutions in the trigonometric approximation [SS00]. The constraint (ii) ensures that the remainder of the Fourier partial sum is smaller than any inverse power of N (by Theorem 2.1(i), all derivatives of the function φ(z) are continuous) [T00]. The constraint (iii) specifies the level of accuracy achieved when the iterations of the method (3.2)-(3.3) are terminated. While we do not proceed with formal analysis of the three numerical factors (see [DS06] for an example of this analysis), we illustrate the weak dependence of three numerical factors on the example of the numerical approximationφ(z) of the exact one-pulse solution (1.3), which exists for c = 36 169 . Numerical    Figure 2 (right) displays convergence of the errors E M = |M n − 1| and E ∞ = u n+1 − u n L ∞ computed dynamically at each n as n increases. We can see that the error E M converges to zero much faster than the error E ∞ , in agreement with the decomposition of the linearized iterative map (3.5) into the one-dimensional projection α n and the infinite-dimensional orthogonal compliment (see the proof of Theorem 3.1). In all further approximations, we will use the error E ∞ for termination of iterations and detecting its minimal values since E ∞ is more sensitive compared to E M . Since the numerical approximationsφ(z) of one-pulse solutions can be computed for any value of c > 0, one can useφ(z) for a given c and compute the effective interaction potential (2.5), which defines the extremal values {L n } n∈N . Theorem 2.3 guarantees that the two-pulse solution φ n (z) consists of two copies of the one-pulse solutions separated by the distance L near the point L n where W ′ (L n ) = 0 and W ′′ (L n ) = 0. Table 1 shows the first four values of the sequence {L n } ∞ n=1 for c = 1 (where s n = L n /2 is the half-distance between the pulses). It also shows the corresponding values from the first four numerical approximations of two-pulse solutions φ n (z) (obtained below) and the computational error computed from the difference of the two numerical approximations. We can see that the error decreases for larger indices n in the sequence {L n } n∈N since the Lyapunov-Schmidt reductions of Theorem 2.3 become more and more accurate in this limit.  By Theorem 2.3(ii), the Jacobian operator H associated with a two-pulse solution φ(z) has one finite negative eigenvalue in the space of even functions H 2 ev (R) and one small eigenvalue which is either negative or positive depending on the sign of W ′′ (L n ). This small eigenvalue leads to either weak divergence or weak convergence of the Petviashvili method in a local neighborhood of φ in H 2 ev (R). Even if the small eigenvalue is positive and the algorithm is weakly convergent, the truncation error from the numerical discretization may push the small eigenvalue to a negative value and lead thus to weak divergence of the iterations. where U 0 (z) is a starting approximation of a sequence {u n (z)} n∈N which converges to the one-pulse solution Φ(z) and s is a parameter defined near L n /2 for the n-th two-pulse solution φ n (z). The left panel shows iterations for s near s 1 and the right panel shows iterations for s near s 2 . Since W ′′ (L 1 ) > 0 and W ′′ (L 2 ) < 0, the two-pulse solution φ 1 (z) leads to the weak divergence of the iteration method (3.2)-(3.3), while the two-pulse solution φ 2 (z) leads to the weak convergence of the method.
At the initial stage of iterations, both errors E M and E ∞ quickly drops to small values, since the starting iterations U 0 (z ∓ s) converge to the one-pulse solutions Φ(z ∓ s) while the contribution from the overlapping tails of Φ(z ∓ s) is negligible. However, at the later stage of iterations, both errors either start to grow (the left panel of Figure 3) or stop to decrease (the right panel). As it is explained above, this phenomenon is related to the presence of zero eigenvalue of H in H 2 ev (R) which bifurcates to either positive or negative values due to overlapping tails of Φ(z ∓ s) and due to the truncation error. At the final stage of iterations on the left panel of Figure 3, the numerical approximation u n (z) converges to the one-pulse solution Φ(z) centered at z = 0 and both errors quickly drop to the numerical zero, which occurs similarly to the right panel of Figure 2. No transformation of the solution shape occurs for large n on the right panel of Figure 3.
The following theorem defines an effective numerical algorithm, which enables us to compute the two-pulse solutions from the weakly divergent iterations of the Petviashvili's method (3.2)-(3.3).
where L(ǫ) has a unit eigenvalue at ǫ = 0, N (v n , ǫ) is C k in v n ∈ H 2 ev with k ≥ 2 such that N (0, 0) = D v N (0, 0) = 0, v n is a perturbation of u n to the fixed point φ n , and ǫ is a small parameter for two-pulse solutions defined in Theorem 2.3. By the Center Manifold Reduction for quasi-linear discrete systems (Theorem 1 in [J03]), there exists a one-dimensional smooth center manifold in a local neighborhood of φ n in H 2 ev (R). Let ξ be a coordinate of the center manifold such that ξ ∈ R, ξ = 0 corresponds to v = 0, and the dynamics on the center manifold is where µ(ǫ) with µ(0) = 1 is an eigenvalue of the linearized operator T (ǫ) at φ n in H 2 ev (R) defined in Theorem 3.1 and f (ξ n , ǫ) is C k in ξ ∈ R with k ≥ 2 and f (0, 0) = ∂ ξ f (0, 0) = 0. Consider the one-parameter starting approximation u 0 (z) = Φ(z − s) + Φ(z + s) in a neighborhood of φ n in H 2 ev (R), where s is close to the value s = s n defined in Theorem 2.3. By the time evolution of the hyperbolic component of v n (see Lemma 2 in [J03]), the sequence v n approaches to the center manifold with the coordinate ξ n . Iterations of ξ n are sign-definite in a neighborhood of ξ = 0. Moreover, there exists s 1 < s n and s 2 > s n , such that the sequences {ξ n (s 1 )} n∈N and {ξ n (s 2 )} n∈N are of opposite signs. By smoothness of v n and ξ n from parameter s, there exists a root s * in between s 1 < s * < s 2 such that ξ n (s * ) = 0 for all n ∈ N.
Remark 3.4 The proof of Theorem 3.3 does not require that the root s * be unique for the oneparameter starting approximation u 0 (z) = Φ(z − s) + Φ(z + s). Our numerical computations starting with a more general approximation (3.7) show that the root s * is unique in a neighborhood of s n .
In order to capture the two-pulse solutions according to Theorem 3.3, we compute the minimum of the error E ∞ for different values of s and find numerically a root s = s * of the function where n 0 is the first iterations after which the value of E infty increases (in case of the left panel of Figure  3) or remains unchanged (in case of the right panel of Figure 3). The numerical root s = s * is found by using the secant method: .   Table 1. We note that the number of iterations N h of the secant method (3.8) decreases with larger values of n, such that N h = 14 for n = 1, N h = 12 for n = 2, N h = 10 for n = 3 and N h = 9 for n = 4, while the number of iterations of the Petviashili's method for each computation does not exceed 100 iterations. Figure 5 shows numerical approximations of the two-pulse solutions for c = 1 and c = 4. We can see from the right panel that two-pulse solutions with c = 4 resemble the two copies of the one-pulse solutions from the left panel of Figure 2, separated by the small-amplitude oscillatory tails.
We note that comparison of the numerical approximations of two-pulse solutions φ n (z) (obtained with a shooting method) and the asymptotic approximations from the effective potential (2.5) is reported in [BC97]. Their numerical results (compare Table 1 in [BC97] with our Table 1) show larger deviations between the two approximations (since their ODE has a different normalization compared to our (1.2)) but the same tendency that the two approximations become identical in the limit of large separation distances s n between the two pulses.
Finally, the three-pulse and multi-pulse solutions of the fixed-point problem (1.2) cannot be approximated numerically with the use of the Petviashili's method (3.2). In the space of even functions H 2 ev (R), the three-pulse solutions have two finite negative eigenvalues and one small eigenvalue, while the stabilizing factor of Theorem 3.1 and the root finding algorithm of Theorem 3.3 can only be useful for one finite negative eigenvalue and one zero eigenvalue. The additional finite negative eigenvalue introduces a strong divergence of the iterative method (3.2) which leads to failure of numerical approximations for three-pulse solutions. This numerical problem remains open for further analysis.   By Theorem 2.3(ii), operator H has two finite negative eigenvalue, a simple kernel and one small eigenvalue, which is negative when W ′′ (L n ) > 0 and positive when W ′′ (L n ) < 0. Persistence (structural stability) of these isolated eigenvalues beyond the leading order (2.7) is a standard property of perturbation theory of self-adjoint operators in Hilbert spaces (see Section IV.3.5 in [K95]). By Theorem 2.3(iii), operator ∂ z H has a pair of small eigenvalues, which are purely imaginary when W ′′ (L n ) > 0 and real when W ′′ (L n ) < 0. We first prove that no other eigenvalues may induce instability of two-pulse solutions (i.e. no other bifurcations of eigenvalues of ∂ z H with Re(λ) > 0 may occur in H 2 (R)). We then prove persistence (structural stability) of the purely imaginary eigenvalues beyond the leading order (2.8). Combined together, these two results lead to the theorem on spectral stability of the two-pulse solution φ n (z) which corresponds to L n with W ′′ (L n ) > 0. is the corresponding eigenfunction for λ ∈ iR + . Assume that no multiple imaginary eigenvalues exist, the kernel of H is simple and P ′ (c) > 0, where P = φ 2 L 2 . Then, where n(H) is the number of negative eigenvalues of H in H 2 (R).
Corollary 4.2 Let φ(z) ≡ Φ(z) be a one-pulse solution defined by Theorem 2.1. Then, it is a spectrally stable ground state in the sense that N real = N comp = N − imag = 0.

Remark 4.3
It is shown in Lemma 4.12 and Remark 4.14 in [CP06] that multiple imaginary eigenvalues may only occur if (Hv, v) = 0 such that n(H) ≥ 2 is a necessary condition for existence of multiple eigenvalues (with P ′ (c) > 0). No multiple imaginary eigenvalues exists for the one-pulse solution Φ(z). Proof. It follows from Theorems 2.1 and 2.3 for sufficiently small O(e −κ 0 Ln ) that P ′ (c) > 0, the kernel of H is simple for W ′′ (L n ) = 0, and the only pair of imaginary eigenvalues with (Hv, v) < 0 in the case W ′′ n (L n ) > 0 is simple. Therefore, assumptions of Theorem 4.1 are satisfied for the two-pulse solutions φ n (z) with W ′′ (L n ) = 0. By the count of Theorem 2.3(ii), n(H) = 3 for W ′′ (L n ) > 0 and n(H) = 2 for W ′′ (L n ) < 0. Furthermore, persistence (structural stability) of simple real eigenvalues of the operator ∂ z H in H 2 (R) follows from the perturbation theory of isolated eigenvalues of non-self-adjoint operators (see Section VIII.2.3 in [K95]).
There exists one uncertainty in Corollary 4.4(ii) since it is not clear if the eigenvalue of negative Krein signature in Theorem 2.3(iii) remains imaginary in N − imag or bifurcates to a complex eigenvalue in N comp . This question is important for spectral stability of the corresponding two-pulse solutions since the former case implies stability while the latter case implies instability of solutions. We will remove the uncertainty and prove that N − imag = 1 and N comp = 0 for sufficiently small O(e −κ 0 L ). To do so, we rewrite the linearized problem (1.5) in the exponentially weighted space [PW94]: The linearized operator ∂ z H transforms to the form which acts on the eigenfunction v α (z) = e αz v(z) ∈ H 2 (R). The absolute continuous part of the spectrum of L α is located at λ = λ α (k), where A simple analysis shows that d dk Re(λ α (k)) = −2αk(10k 2 − 10α 2 + 3), d dk Im(λ α (k)) = c − 3α 2 + 5α 4 + 3k 2 (1 − 10α 2 ) + 5k 4 .
The following lemma gives a precise location of the dispersion relation λ = λ α (k) on λ ∈ C.
Proof. By Lemma 4.6, λ 0 is also an eigenvalue of L α in H 2 (R) for sufficiently small α. Let α be fixed in the bound (4.5). There exists a small neighborhood of λ 0 , which is isolated from the absolute continuous part of the spectrum of L α . By the perturbation theory of isolated eigenvalues of non-selfadjoint operators (see Section VIII.2.3 in [K95]), there exists a simple eigenvalue λ δ of ∂ z (H + δV (z)) in H 2 α (R) for the same value of α and sufficiently small δ in a local neighborhood of λ 0 , such that |λ δ − λ 0 | ≤ Cδ for some C > 0.
Corollary 4.9 Let φ(z) ≡ φ n (z) be a two-pulse solution defined by Theorem 2.3 that corresponds to L n with W ′′ (L n ) > 0. Then, it is spectrally stable in the sense that N real = N comp = 0 and N − imag = 1 for sufficiently small O(e −κ 0 Ln ).
Remark 4.10 Using perturbation theory in H 2 α (R) for a fixed value α > 0, one cannot a priori exclude the shift of eigenvalue λ 0 to λ δ with Re(λ δ ) > 0. Even if v 0 (z) for λ 0 contains no term e ik 0 z as z → −∞ (see Lemma 4.6), the eigenfunction v δ (z) for λ δ may contain the term e ik δ z as z → −∞ with Im(k δ ) < 0 and lim δ→0 k δ = k 0 ∈ R. However, when Theorem 4.8 holds (that is under the assumptions that v 0 ∈ H 2 (R) and (Hv 0 , v 0 ) < 0), the eigenvalue λ δ remains on iR and the eigenfunction v δ (z) must have no term e ik δ z with k δ ∈ R as z → −∞ for any sufficiently small δ. The hypothetical bifurcation above can however occur if v 0 / ∈ H 2 (R) but v 0 ∈ H 2 α (R) with α > 0. We do not know any example of such bifurcation in the case v 0 / ∈ H 2 (R) but (Hv 0 , v 0 ) ≥ 0.
We confirm results of Corollaries 4.4 and 4.9 with numerical computations of eigenvalues in the linearized problem (1.5). Throughout computations, we use the values α = 0.04 and c = 1, which satisfy the constraint (4.5). The spectra of the operators H in H 2 (R) and ∂ z H in H 2 α (R) are computed by using the Fourier spectral method. This method is an obvious choice since the solution φ(z) is obtained by using the spectral approximations in the iterative scheme (3.2)-(3.3). As in the previous section, we use numerical parameters d = 100, h = 0.01 and ε = 10 −15 for the Petviashvili's method (3.2)-(3.3).
Eigenvalues of the discretized versions of the operators H and L α are obtained with the MATLAB eigenvalue solver eig. The spectra are shown on Figure 6 for the two-pulse solution φ 1 (z) and on Figure 7 for the two-pulse solution φ 2 (z). The inserts show zoomed eigenvalues around the origin and the dotted line connects eigenvalues of the discretized operators that belong to the absolutely continuous part of the spectra. Figures 6 and 7 clearly illustrate that the small eigenvalue of H is negative for φ 1 (z) and positive for φ 2 (z), while the pair of small eigenvalues of L α is purely imaginary for φ 1 (z) and purely real for φ 2 (z). This result is in agreement with Corollaries 4.4 and 4.9. We have observed the same alternation of small eigenvalues for two-pulse solutions φ 3 (z) and φ 4 (z), as well as for other values of parameters c and α.
The numerical discretization based on the Fourier spectral method shifts eigenvalues of the operators H and L α . In order to measure the numerical error introduced by the discretization, we compute the numerical value for the "zero" eigenvalue corresponding to the simple kernel of H and the double zero eigenvalue of L α . Table II shows numerical values for the "zero" and small eigenvalues for two-pulse solutions φ n (z) with n = 1, 2, 3, 4. It is obvious from the numerical data that the small eigenvalues are still distinguished (several orders higher) than the numerical approximations for zero eigenvalues for n = 1, 2, 3 but they become comparable for higher-order two-pulse solutions n ≥ 4. This behavior is understood from Theorem 2.3 since the small eigenvalues becomes exponentially small for larger values of s (larger n) in the two-pulse solution (2.4) and the exponentially small contribution is negligible compared to the numerical error of discretization.
We have confirmed numerically the analytical predictions that all two-pulse solutions corresponding to the points L n with W ′′ (L n ) < 0 (which are maxima of the effective interaction potential) are unstable with a simple real positive eigenvalue, while all two-pulse solutions corresponding to the points L n with W ′′ (L n ) > 0 (which are minima of the effective interaction potential) are spectrally stable. The stable two-pulse solutions are not however ground states since the corresponding linearized problem has a pair of purely imaginary eigenvalues of negative Krein signature.

Nonlinear dynamics of two-pulse solution
The Newton's particle law (2.15) is a useful qualitative tool to understand the main results of our article. Existence of an infinite countable sequence of two-pulse solution {φ n (z)} n∈Z is related to existence of extremal points {L n } n∈Z of the effective potential function W (L), while alternation of stability and instability of the two-pulse solutions is related to the alternation of minima and maxima points of W (L). It is natural to ask if the Newton's law (2.15) extends beyond the existence and spectral stability analysis of two-pulse solutions in the fifth-order KdV equation (1.1). In particular, one can ask if the  purely imaginary (embedded) eigenvalues of the linearized problem (1.5) lead to nonlinear asymptotic stability of two-pulse solutions or at least to their nonlinear stability in the sense of Lyapunov. From a more technical point of view, one can ask whether the Newton's law (2.15) serves as the center manifold reduction for slow nonlinear dynamics of two-pulse solutions in the PDE (1.1) and whether solutions of the full problem are topologically equivalent to solutions of the Newton's law. While we do not attempt to develop mathematical analysis of these questions, we illustrate nonlinear dynamics of two-pulse solutions with explicit numerical simulations.
The numerical pseudo-spectral method for solutions of the fifth-order KdV equation (1.1) is described in details in [MV99]. The main idea of this method is to compute analytically the linear part of the PDE (1.1) by sing the Fourier transform and to compute numerically its nonlinear part by using an ODE solver. Letû(k, t) denote the Fourier transform of u(x, t) and rewrite the PDE (1.1) in the Fourier domain:û t = i(k 3 + k 5 )û − ik u 2 . (5.1) In order to compute u 2 (k, t) we evaluate u 2 (x, t) on x ∈ R and apply the discrete Fourier transform. Substitutionû = s(k, t)e i(k 3 +k 5 )t transforms the evolution equation (5.1) to the form: s t = −ike −i(k 3 +k 5 )t u 2 (k, t). (5. 2) The fourth-order Runge-Kutta method is used to integrate the evolution equation (5.2) in time with time step △t. To avoid large variations of the exponent for large values of k and t, the substitution above is updated after m time steps as follows: u = s m (k, t)e i(k 3 +k 5 )(t−m△t) , m△t ≤ t ≤ (m + 1)△t. (5. 3) The greatest advantage of this numerical method is that no stability restriction arising from the linear part of (5.1) is posed on the timestep of numerical integration. On contrast, the standard explicit method for the fifth-order KdV equation (1.1) has a serious limitation on the timestep of the numerical integration since the fifth-order derivative term brings stiffness to the evolution problem. The small timestep would be an obstacle for the long time integration of the evolution problem due to accumulation of computational errors.
Numerical simulations of the PDE (5.1) are started with the initial condition: u(x, 0) = Φ(x − s) + Φ(x + s), (5.4) where Φ(x) is the one-pulse solution and 2s is the initial separation between the two pulses. The onepulse solution Φ(x) is constructed with the iteration method (3.2)-(3.3) for c = 4. The numerical factors of the spectral approximation are L = 100, N = 2 12 , ε = 10 −15 , while the timestep is set to △t = 10 −4 . Figure 8 shows six individual simulations of the initial-value problem (5.1) and (5.4) with s = 2.3, s = 2.8, s = 3.6, s = 4.2, s = 4.5 and s = 4.7. Figure 9 brings these six individual simulations on the effective phase plane (L,L) computed from the distance L(t) between two local maxima (humps) of the two-pulse solutions.
When the initial distance (s = 2.3) is taken far to the left from the stable equilibrium point (which corresponds to the two-pulse solution φ 1 (x)), the two pulses repel and diverge from each other (trajectory 1). When the initial distance (s = 2.8) is taken close to the left from the stable equilibrium point, we observe small-amplitude oscillations of two pulses relative to each other (trajectory 2). When the initial distances (s = 3.6) and (s = 4.2) are taken to the right from the stable equilibrium point, we continue observing stable oscillations of larger amplitudes and larger period (trajectories 3 and 4). The oscillations are destroyed when the initial distances are taken close to the unstable equilibrium point (which corresponds to the two-pulse solution φ 2 (x)) from either left (s = 4.5) or right (s = 4.7). In either case, the two pulses repel and diverge from each other (trajectories 5 and 6). Ripples in the pictures are due to radiation effect and the numerical integration does not make sense after t ≈ 500 because the ripples reach the left end of the computational interval and appear from the right end due to periodic boundary conditions. The numerical simulations of the full PDE problem (1.1) indicate the validity of the Newton's particle law (2.15). Due to the energy conservation, all equilibrium points in the Newton's law are either centers or saddle points and the center points are surrounded by closed periodic orbits in the interior of homoclinic loops from the stable and unstable manifolds of the saddle points. Trajectories 2,3, and 4 are taken inside the homoclinic orbit from the saddle point corresponding to φ 2 (x) and these trajectories represent periodic oscillations of two-pulse solutions near the center point corresponding to φ 1 (x). Trajectories 1 and 6 are taken outside the homoclinic orbit and correspond to unbounded dynamics of two-pulse solutions. The only exception from the Newton's law (2.15) is trajectory 5, which is supposed to occur inside the homoclinic loop but turns out to occur outside the loop. This discrepancy can be explained by the fact that the Newton's law (2.15) does not exactly represent the dynamics of the PDE (5.1) generated by the initial condition (5.4) but it corresponds to an asymptotic solution after the full solution is projected into the discrete and continuous parts and the projection equations are truncated (see details in [JFGS06] in the context of the NLS equations).
Summarizing, we have studied existence, spectral stability and nonlinear dynamics of two-pulse solutions of the fifth-order KdV equation. We have proved that the two-pulse solutions can be numerically approximated by the Petviashili's iterative method supplemented with a root finding algorithm. We have also proved structural stability of embedded eigenvalues with negative Krein signature and this result completes the proof of spectral stability of two-pulse solutions related to the minima points of the effective interaction potential. The validity of the Newton's particle law is illustrated by the full numerical simulations of the fifth-order KdV equation (1.1) which show agreement of slow nonlinear dynamics of two-pulse solutions with predictions of the Newton's particle law.   Figure 8, where L is the distance between two pulses. The six trajectories correspond to the six individual simulations in the order described in the text. The black dots denote stable and unstable equilibrium points which correspond to the two-pulse solutions φ 1 (x) and φ 2 (x).