{'fraction_non_alphanumeric': 0.05568028047730348, 'fraction_numerical': 0.020374584819781032, 'mean_word_length': 4.3188026498732315, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Conjugate gradient method for finding fundamental solitary waves

The Conjugate Gradient method (CGM) is known to be the fastest generic iterative method for solving linear systems with symmetric sign definite matrices. In this paper, we modify this method so that it could find fundamental solitary waves of nonlinear Hamiltonian equations. The main obstacle that such a modified CGM overcomes is that the operator of the equation linearized about a solitary wave is not sign definite. Instead, it has a finite number of eigenvalues on the opposite side of zero than the rest of its spectrum. We present versions of the modified CGM that can find solitary waves with prescribed values of either the propagation constant or power. We also extend these methods to handle multi-component nonlinear wave equations. Convergence conditions of the proposed methods are given, and their practical implications are discussed. We demonstrate that our modified CGMs converge much faster than, say, Petviashvili's or similar methods, especially when the latter converge slowly.


Introduction
Solitary waves of most nonlinear evolution equations can be found only numerically. In one spatial dimension, a variety of numerical methods for finding solitary waves are available, with the shooting and Newton's methods being probably the most widely used ones. However, in two and three spatial dimensions, the shooting method is no longer applicable, and the complexity of programming of Newton's method increases considerably compared to the one-dimensional case. Also, the accuracy of Newton's method is only algebraic in the size of the spatial discretization step ∆x (typically, it is O (∆x) 2 ). This accuracy may be too low for some applications, since in two or three spatial dimensions, one tends to take larger discretization steps than in one dimension in order to limit the computational time and the size of stored numeric arrays. Therefore, it is desirable to have a method that would: (i) be easier to program than Newton's method; (ii) have the exponential accuracy in ∆x, as spectral methods; and (iii) converge sufficiently fast under known conditions.
A number of such methods are available. The well-known iterative method proposed by Petviashvili [1] can be used to find stationary solitary waves of equations with power-law nonlinearity: − M u + u p = 0 , u(x) → 0 as |x| → ∞, (1.1) where u is the real-valued field of the solitary wave, M is a positive definite and self-adjoint differential operator with constant coefficients, and p is a constant. For example, the solitary wave of the nonlinear Schrödinger equation in D spatial dimensions, iU t + ∇ 2 U + |U | 2 U = 0 , U (|x| → ∞) → 0 , upon the substitution U (x, t) = e iµt u(x) reduces to Eq. (1.1) with p = 3 and 3) The parameter µ > 0 is referred to as the propagation constant of the solitary wave.
Convergence conditions of the Petviashvili method for Eq. (1.1) were established in [2].
Studies of solitary waves in nonlinear photonic lattices and Bose-Einstein condensates motivated the interest in finding solitary waves of equations that have a more general form than (1.1): − M u + F (u, x) = 0 , u(x) → 0 as |x| → ∞, (1.4) where F (u, x) is a real function. For example, solitary waves in nonlinear photonic lattices with Kerr nonlinearity satisfy the following equation of the form (1.5): ∇ 2 u + V 0 (cos 2 x + cos 2 y) u + u 3 = µu , (1.5) where the second term accounts for the effect of the periodic potential provided by the lattice. Various modifications of the Petviashvili method for obtaining solutions of (1.4) were proposed; see, e.g., [3,4,5,6]. However, convergence conditions of those modifications were not studied. In Ref. [7,8], J. Yang and the present author proposed two alternative iterative methods which satisfied conditions (i)-(iii) stated before Eq. (1.1); in particular, the conditions of their convergence to fundamental (see below) [7] and both fundamental and non-fundamental solitary waves (also known as ground and excited states, respectively) [8] were given. In Refs. [9,8], a technique referred there as mode elimination was proposed, which was shown to provide considerable acceleration of iterative methods when the later converged slowly to the respective solitary wave.
In this paper, we propose yet another family of numerical methods for finding fundamental solitary waves of Hamiltonian nonlinear evolution equations. Consider a linear system Ay − b = 0, (1.6) for the unknown vector y, where A is a real symmetric matrix and b is a known vector. Then the methods proposed here are extensions of the well-known Conjugate Gradient method (CGM) for the linear system (1.6) to the nonlinear boundary-value problem (1.4), and they converge even faster than, e.g., the generalized Petviashvili method [7] accelerated by mode elimination [9].
The main part of this paper is organized as follows. In Section 2, we explain how the methods introduced in [7,8,9] for the nonlinear problem (1.4) are related to some well-known iterative methods of solving the linear system (1.6) with a real and symmetric matrix A. This comparative review will lead us to explaining the main issue that needs to be resolved in order to extend the CGM in such a way that it could find fundamental solitary waves. In Section 3, we present such an extension of the CGM for the case where the solitary wave of the single-component Eq. (1.4) has a prescribed value of the propagation constant µ. In Section 4, we present a modification of the CGM for finding fundamental waves of (1.4) with a prescribed value not of the propagation constant but of power, defined as P = u 2 (x) dx , (1.7) where the integration is over the entire spatial domain. Let us note that another method for finding solitary waves with a specified value of power was previously considered in, e.g., [10,11,12,13], with its convergence conditions for (1.4) being established in [13]. We will refer to this method as the imaginary-time evolution method. In Section 5, we generalize the modified CGMs of Sections 3 and 4 to multi-component equations. In Section 6, we compare the performances of the generalized Petviashvili and imaginary-time evolution methods accelerated by the mode elimination technique [9,8] with the performances of the respective modified CGMs from Sections 3-5. In Section 7 we summarize this work.
Appendix 1 contains proofs and discussions of some auxiliary results of Sections 3 and 4.
Appendix 2 contains an alternative method to the modified CGM of Section 3. Appendix 3 contains a sample code of a modified CGM for a two-component equation considered in Section 6. Before concluding this Introduction, we need to discuss what we will refer to as a fundamental solitary wave. In fact, we are not aware of a universal definition of this term.
In many situations, one can intuitively identify a solitary wave as being fundamental if it satisfies two easily-verifiable conditions. First, it is to have one main peak, with smaller peaks possibly existing around the main one. Second, for envelope solitary waves (e.g., those satisfying the complex-valued Eq. (1.2)), the propagation constant of the fundamental wave must lie in the semi-infinite spectral bandgap; see Fig. 1 in Section 6. (Typically, only equations with periodic potentials or systems describing more than two linearly coupled waves propagating with different group velocities have more than one spectral bandgap.) Alternatively, for carrier (e.g., Korteweg-de Vries-type) solitary waves, the velocity, rather than the propagation constant, of a fundamental wave must lie in the semi-infinite gap. More rigorous characterization of fundamental waves is possible in special cases. For instance, when the spatial operator in the wave equation is the Laplacian (as, e.g., in (1.2) or (1.5)), the fundamental wave is known to be nonzero for all x (see, e.g., [14]); however, for a more general operator as, e.g., in the Kadomtsev-Petviashvili equation, the fundamental wave may have zeros [15]. For some equations (see, e.g., [16,17]), it was rigorously shown that the fundamental solution minimizes a functional known as the action.
The property of an S-component fundamental solitary wave that we will rely on in this paper is that for many such waves, the linearized operator of the stationary equation has no more than S positive eigenvalues. To the author's knowledge, there is no general proof of this property even for a single-component Eq. (1.4), let alone for the S-component case.
(See, however, a recent review of solitary wave stability in lattices [18].) Therefore, we simply declare this property to be our working definition of a fundamental solitary wave for the purpose of this paper. Again, we do not imply that any solitary wave that has the two intuitive properties described at the beginning of the previous paragraph also satisfies the above property about the positive eigenvalues of the corresponding linearized operator. Rather, we only say that it is only for the waves that do have the latter property that the numerical methods developed in this paper will be guaranteed to converge (with possible restrictions, as described in subsequent sections). For solitary waves whose corresponding linearized operators have more positive eigenvalues than the number of the wave's component, our modified CGMs are not guaranteed to converge.
2 Comparative review of iterative methods for nonlinear problem (1.4) and for linear problem (1.6) The reader who is not interested in such a review and wants to see the description of new algorithms proposed in this paper may skip to Section 3.
For reference purposes, let us rewrite (1.4) as An important role in our analysis will be played by the linearized operator L of the nonlinear Eq. (2.1). Let u be the exact solution of (2.1), u n be the approximation of that solution obtained at the nth iteration, and u n = u +ũ n , where |ũ n | ≪ |u|. (2.2) Then For example, for the nonlinear equation (1.5), the linearized operator is In our discussion, the linear system (1.6) is a counterpart of the nonlinear equation (2.1) and matrix A is the counterpart of the linearized operator L. For Hamiltonian wave equations that give rise to (2.1), operator L is always self-adjoint, which is the counterpart of matrix A being real and symmetric. Having stated this correspondence between the nonlinear problem (2.1) and its linear counterpart (1.6), we will, from now on, refer to them interchangeably, because our statements will apply to both of them equally.
As we will explain below, the linearized forms of the iterative methods proposed in [1,7,8] for solving the nonlinear problems (1.1) and (1.4) are related to Richardson's method 1 (see, e.g., [19], p. 274) of solving the linear problem (1.6): The necessary condition for the iteration scheme (2.5) to converge to the solution, y = A −1 b, is that A be negative definite. Choosing the parameter ∆τ in (2.5) so that ∆τ < −2/λ min , where λ min is the minimum (i.e., most negative) eigenvalue of A, one ensures that Richardson's method for (1.6) with a negative definite A converges.
A naive counterpart of (2.5) for the nonlinear problem (2.1) is However, in most cases it will not converge to a solitary wave. Indeed, from (2.3), the linearized form of (2.6) isũ n+1 =ũ n + Lũ n ∆τ .
This equation is a counterpart of (2.5), and so a necessary condition for its convergence is that operator L be negative definite. However, for most fundamental single-component 1 also known as Picard's or fixed-point iterative method solitary waves, the corresponding linearized operator L has one positive eigenvalue (see the last paragraph of the Introduction), which causes divergence of small deviationsũ n in (2.7) and, thereby, divergence of the iterations (2.6).
For example, in [8], 8) so that L (mod) = −L 2 . In terms of the linear system (1.6) where A is real symmetric but not sign definite, this is equivalent to applying Richardson's method (2.5) to the so called normal equation 9) where the matrix on the left-hand side is always negative definite.
In [7], we proposed an alternative idea of constructing a negative definite L (mod) . This idea will play an important role in the development of this paper, and here we will explain it as applied to the linear system (1.6). As noted above, for most single-component fundamental solitary waves, the linearized operator L has only one positive eigenvalue. Accordingly, let λ (1) be the only positive eigenvalue of matrix A in (1.6), with the corresponding eigenvector being z (1) . In the context of the nonlinear equation (2.1), the eigenvector of L corresponding to its only positive eigenvalue is related to the solution u of (2.1) [7], as we will explain later. Hence it can be considered as approximately known once the iterative solution u n is sufficiently close to u, which we will always assume to be the case (see (2.2)). Now, instead of (1.6), consider an equivalent system where I is the identity matrix and γ is any number greater than one. Note that the term z (1) z (1) T is the projection matrix onto the unstable direction z (1) . Using the spectral decomposition of a real symmetric p × p matrix A: 11) where the eigenvectors z (i) form an orthonormal set: 12) it is straightforward to show that (2.10) and hence (1.6) is equivalent to Since λ (1) (1 − γ) < 0 and λ (i) < 0 for i ≥ 2 by our assumptions about A and γ, the matrix A (mod) in (2.13b) is negative definite, and hence Richardson's method (2.5) for it can converge.
Note that if one chooses [7] γ = 1 + 1/(λ (1) ∆τ ) (2.14) and uses the resulting A (mod) in Richardson's method (2.5), one thereby forces the projection of the errorỹ n ≡ y n − y on the direction of the unstable eigenvector z (1) to be zero at every iteration. This idea of zeroing out the projection of the iteration error on the unstable direction lies behind the original Petviashvili method [2] and its generalization for Eq. (1.4) [7]. This is also the idea that we will use later in this paper for the development of modified CGMs.
So far we have discussed how Richardson's method for (1.6) with a sign indefinite A can be made to converge. The next issue is, how fast it converges. This is quantified by the convergence factor R, such that the number of iterations required for the iteration error to decrease by a factor of e is (asymptotically) inversely proportional to logR. For the optimal value of ∆τ , the convergence factor of Richardson's method is (see, e.g., Sec. 5.2.2 in [20], or [8]): 15) where the condition number is related to the eigenvalues of A (mod) : cond(A (mod) ) = λ min /λ max (2.16) (recall that all the eigenvalues of A (mod) are negative and so cond(A (mod) ) > 1). For cond(A) ≫ 1, convergence to a prescribed accuracy iterations. That is, the greater the condition number, the slower the convergence of the iterative method. A common method to reduce the condition number is to use a preconditioner, which can drastically lower |λ min | (see, e.g., Lecture 40 in [19]). The preconditioned Richardson's method is 17) where B is the preconditioning matrix 2 . For the Richardson-type method (2.6) with L (0) replaced by L (0,mod) , a convenient preconditioning operator, N , has a form similar to that of M . For example, if M is as in (1.3), the convenient form for N is: where c is some constant.
Even if the condition number is lowered by reducing |λ min | via preconditioning, the convergence of an iterative method may still be slow due to |λ max | being close to zero.
In the context of the nonlinear problem (1.4), examplified by (1.5), this occurs when the propagation constant µ is close to the edge of the spectral bandgap of the linear potential (i.e., V 0 (cos 2 x+cos 2 y) in (1.5)); see, e.g., Fig. 4(a) in [21] or Fig. 1 in Section 6 below. In [9,8] we proposed to accelerate the iterations by eliminating the slowest-decaying eigenmode, i.e., the component of the errorũ n that corresponds to the eigenvalue λ max . This effectively removes this slowest eigenmode from the spectrum of L (mod) and thereby increases the magnitude of its maximum eigenvalue. This, in its turn, increases the convergence rate of the method via (2.16) and (2.15). The slowest mode can be eliminated in exactly the same way as the unstable mode (i.e., the mode with the positive eigenvalue λ (1) ) in (2.10). The slowest mode can be approximated by [9] z slow ∝ y n − y n−1 , (2.19) and a preconditioned Richardson's method is then applied not to (1.6) but to 20) where γ slow is computed similarly to (2.14). The explicit form of a counterpart of (2.20) for the nonlinear problem (2.1) will be given in Section 6.
In [7,8] we demonstrated that when the convergence of an iterative method is slow, the slowest mode elimination technique described above can considerably (by a factor of five to ten times) accelerate the convergence. However, it is well known (see, e.g., [19], p. 341 and [20], p. 451) that the fastest generic method for a linear system (1.6) with a real symmetric and sign definite matrix A is the Conjugate Gradient method (CGM), whose algorithm is given by Eqs. (3.1) of the next section. The convergence factor of the CGM satisfies (see, e.g., [19], p. 299): which for matrices with large condition numbers implies considerably faster convergence than does (2.15). Therefore, if one could modify the CGM so that it would be guaranteed to converge when A (or the linearized operator L) is not negative definite but has one positive eigenvalue, then such a modified CGM is expected to be faster than a Richardsontype method accelerated by mode elimination.
In this paper we present several versions of such a modified CGM for nonlinear problem (2.1) and its generalizations and compare them with the Richardson-type methods of [7,13] accelerated by mode elimination. We find that, as expected, the modified CGMs are the faster of these two groups of methods. As we noted after Eq. (2.14), the main idea behind all these versions of the modified CGM is to eliminate the component of the errorũ n that is aligned along certain unstable eigenmode(s) of the linearized iteration operator. Our modified CGMs are guaranteed to converge in almost all situations where a fundamental solitary wave with a prescribed (set of) propagation constant(s) is sought. Moreover, even for those rare situations where our guarantee is formally void, we will show a simple and practical way to still make the corresponding modified CGM converge. For the situations where a solitary wave with a prescribed (set of) power(s) is sought, our modified CGM is guaranteed to converge in the same parameter space where the Richardson-type method [13,22] converges.

Modified CGM for solitary waves with prescribed propagation constant
In this section, we will first state the algorithm of the CGM for the linear Eq. (1.6). Then, we will briefly review the generalized Petviashvili method [7] for solitary waves and re-cast it in a form more convenient for the present analysis. Finally, we will propose a modification of the CGM for finding fundamental solitary waves with a prescribed value of the propagation constant. All analysis in this section is done for the single-component case; its generalization for multi-component solitary waves is presented in Section 5.

CGM for the linear Eq. (1.6)
The original CGM algorithm for Eq. (1.6) with a real symmetric and sign definite matrix A and starting with an initial guess y 0 , is: Here a, b denotes the natural inner product between the real-valued vectors a and b.
Equations (3.1a) define the initial residual vector r 0 and the initial search direction d 0 .
Equation (3.1c) updates the solution by adding to the previous solution a vector α n d n along the search direction d n . The value of α n is set by (3.1b) to guarantee the orthogonality of the new residual, defined by (3.1d), to the search direction d n : Finally, Eq. (3.1f) updates the search direction, where β n is set by (3.1e) so that The condition that A be real symmetric is inherently required for the derivation of algorithm (3.1). In particular, it guarantees that the orthogonality relations (3.2) and (3.3) between quantities at two consecutive iterations imply more general orthogonality relations among such quantities at all iterations (see, e.g., [23] or [19], p. 295): The orthogonality relation (3.5) with a sign definite A ensures that the CGM at each iteration produces a search direction that is linearly independent from all previous search directions. This implies that the error (in a certain norm) decreases with each iteration; this fact can be restated by saying that the CGM is an optimal iterative method for Eq. (1.6) with a sign definite matrix A (see, e.g., [19], p. 296).
The condition that A be sign definite is also needed to guarantee that the denominators of α n and β n never vanish (or, for practical purposes, never become too close to zero). Thus, applying the CGM to an equation with a sign indefinite matrix A would come without the guarantee that this method would be optimal and that it would even converge. However, it should be emphasized that the CGM may converge when applied to problem (1.6) with a sign indefinite matrix A (see, e.g., [23], or [19], p. 301). This is in stark contrast to Richardson's method (2.5), which, for a generic initial condition y 0 , is guaranteed to diverge in such a case.
The CGM (3.1) is straightforwardly extended to solve a nonlinear equation (3.1d) and A with K in (3.1b) and (3.1e). In fact, for the CGM applied to a nonlinear problem, several choices of parameter β n that are equivalent in the linear approximation, will produce different versions of the method. Two most widely used such versions are known as Fletcher-Reeves and Polak-Ribière (see, e.g., [24]). In the following we will not attempt to compare these two versions because the focus of this paper is not on the extension of the linear algorithm (3.1) to the nonlinear case but on the extension of that algorithm to the case when the counterpart of matrix A has one eigenvalue of the opposite sign than the rest of the spectrum.
Let us also note that, unlike the algorithm (3.1) for linear equations, its extension to nonlinear problems described in the previous paragraph may fail if the linearized operator K is negative definite but has one or more eigenvalues close to zero. This can occur if, at some iteration, the search direction d n becomes aligned primarily along the eigenmode(s) corresponding to the small eigenvalue(s). Then the denominator in (3.1b) will be small and hence the increment α n d n of the solution in (3.1c), large. If the nonlinear equation admits more than one solution for the considered values of its parameters, which typically occurs near a bifurcation, then a large shift, α n d n , of the solution may cause the iterations to converge to a different solution than initially intended, or to diverge. Such a failure of the method can be avoided simply by changing the initial condition y 0 . We will discuss this in more detail in Section 6.

Review of the generalized Petviashvili method
In [7], we proposed the following iterative method for finding the fundamental solitary wave u of the nonlinear problem (2.1): and the integration is over the entire spatial domain. The role of the γ-term is discussed at the end of this subsection. In (3.7), a self-adjoint and positive definite differential operator N with constant coefficients is chosen so that u is an approximate eigenvalue of N −1 L corresponding to its largest eigenvalue λ (1) : The Sylvester law of inertia (see, e.g., [25]) guarantees that the signs of the eigenvalues of  [2].) The constant γ in (3.7) is given by (2.14) and In practice, N , γ, and λ (1) are computed until the iteration error reaches some predefined tolerance, after which their last-computed values are used for all subsequent iterations. Since N is a differential operator with constant coefficients, it has a simple representation in the Fourier space. Therefore, quantities like N −1 L (0) u n , N u n , etc. are easily computed using the direct and inverse Fast Fourier Transforms, which are available as built-in commands in all major computing software.
In what follows we will frequently refer to the linearized form of (3.7), which is: Note that the operator N −1 L in the leading term of (3.10) is not self-adjoint. Therefore, we introduce a commonly employed change of variables to rewrite that equation in a form involving a self-adjoint operator. Namely, we denote Note that operator K is self-adjoint, since L is self-adjoint and N is self-adjoint and positive definite. Let us stress that K and v have been introduced only for the convenience of carrying out the subsequent analysis with self-adjoint operators. Once the derivation of the algorithm in this framework is complete, we will recast it in terms of the original variable u n and operator L.
With the notations (3.11), Eqs. (3.8), (3.7), and (3.10) are rewritten as: Equations (3.13) and (3.14) are the nonlinear and linearized counterparts of the linear problem (2.10) with a symmetric matrix A. The γ-terms in these equations are needed to eliminate the component of the errorṽ n+1 along the unstable "direction" v; see (3.12) and the text after (3.8). If the equality in (3.12) (and (3.8)) is exact, which occurs only for equations with power-law nonlinearity, (1.1), this elimination is complete. However, for the majority of equations (1.4) with a more complicated nonlinear term, v is not an exact eigenmode of K and hence the γ-terms in (3.13) and (3.14) strongly suppress, but do not completely eliminate, the component along the unstable eigenmode. Nonetheless, this strong suppression is sufficient to turn the unstable eigenmode into a stable one.
In the remainder of this section, we will continue using Eq. (3.14) with the self-adjoint operator K for all derivations, and only at the final steps of those derivations will convert to the original variable u n .

A modified CGM
As we noted at the end of Sec. 3.1, the straightforward nonlinear generalization of the CGM applied to Eq. (3.6) whose linearized operator K is self-adjoint but not sign definite, may diverge. For equations that we consider in Section 6 below, this method indeed diverges. To modify the method so that it could converge, one can use the idea stated before Eq. (2.8).
Namely, replace (3.6) with an equivalent equation K (0,mod) v = 0 such that the corresponding linearized operator K (mod) is self-adjoint and sign definite. In Appendix 1 we show that the expression in (3.15a) below, which mimics that on the right-hand side of (3.13), comes close to satisfying these properties: (i) it is equivalent to K (0) v, (ii) its linearized operator K (mod) is approximately self-adjoint, and (iii) its eigenvalues that are not too close to zero are all negative. Statements (ii) and (iii) are approximate 3 rather than exact because so is relation (3.12). We will discuss the implications of the approximate character of statement (ii) after we present the algorithm of the modified CGM.
Thus, denoting we propose the following modified CGM for finding the fundamental solitary wave of (3.6).
It mimics algorithm (3.1) with the following modifications: , and vectors r and d are renamedr andd. Upon the change of variables 16) and (3.11), this modified CGM is: In a practical implementation of algorithm (3.17), iterations should be carried out with the generalized Petviashvili method (3.7) until the error reaches some predefined tolerance (usually, between 1 and 10%), so that the computed parameters of operator N and the constants λ (1) and Γ may be considered as known with the corresponding accuracy. After that, the iterations should be carried out with the modified CGM (3.17), which is expected to converge considerably faster than the generalized Petviashvili method.
Operator K (mod) in (3.15c), and hence the quadratic form d n , K (mod)d n in the counterparts of (3.17b,e), would be guaranteed to be negative definite if in relation (3.12) (or, equivalently, in (3.8)), the equality held exactly. However, this takes place only for equations with power-law nonlinearity, (1.1). Therefore, for some nonlinear wave equations, d n , K (mod)d n could be sign indefinite and, importantly, arbitrarily close to zero. In such situations algorithm (3.17) could fail, and we will now point out when this should be expected. As we explain in Appendix 1, the slight non-self-adjointness of K (mod) can cause the quadratic form to become too close to zero when some of the negative eigenvalues of the linearized operator L of (2.1) are close to zero. In single-component equations, this typically occurs when the solitary wave bifurcates from the edge of the continuous spectrum, while for multi-component equations this can also happen when, say, an asymmetric wave bifurcates from one with a symmetry (see, e.g., [26] and references therein).
Our numerical experiments confirm that algorithm (3.17) converges whenever the eigenvalues of L are sufficiently far from zero 4 , and also that is may diverge in the opposite situation, as described above. Recall, from the last paragraph of Section 3.1, that in exactly the same situation, a CGM for a nonlinear problem can fail even if the quadratic form d n , K (mod)d n were negative definite. Thus, there are two different reasons that can cause the modified CGM to fail, but, fortunately for the applications of the method, they both can do so in the same situation, namely, when L has near-zero eigenvalues. Moreover, they both require that at some iteration, the search direction become primarily "aligned" with the eigenmode of L corresponding to a small eigenvalue. Therefore, a simple and practical way to avoid such divergence is to change the initial condition in a certain manner. In Section 6 we show that this way does indeed work.
Before concluding this section, let us note that one could develop a modified version of the CGM based on a different idea than the elimination of the unstable eigenmode from the underlying operator, as in (3.15). Instead, one can employ the original operators K (0) and K but force all search directionsd n to be orthogonal to the unstable eigenmode of K approximated by the exact solitary wave v. Since (3.12) holds approximately, the search 4 Note that if the solitary wave is translationally invariant, L has a zero eigenvalue corresponding to the respective translational eigenmode. As we show in Appendix 1, such a mode presents no problem for convergence of the iterations since it would simply slightly shift the solitary wave. Therefore, in the following we ignore the possible presence of such zero eigenvalues of L.
directions can be made only approximately orthogonal to the true unstable eigenmode of K. One can show that this causes this other modified CGM to be prone to failure under the same circumstances as the modified CGM (3.17), i.e. when L has eigenvalues that are close to zero. Our numerical experiments show that the simple trick that can help algorithm (3.17) overcome that failure, is not as effective for the other modified CGM. Moreover, the coding of the latter method is slightly more complicated than that of (3.17). For these reasons, we only present that alternative modified CGM in Appendix 2, but do not discuss it further in this paper.

Modified CGM for solitary waves with prescribed power
We first review the preconditioned imaginary-time evolution method (ITEM) [13], which is a Richardson-type iterative method for finding solitary waves with a specified value of power, and then present a modified CGM that converges under the same conditions as the ITEM, but faster.

Review of the ITEM
When one seeks a solitary wave with a prescribed value of power (1.7), problem (2.1) should be redefined because the propagation constant µ, which enters L (0) via operator M (see (1.3)), is no longer known exactly. Instead, it is estimated at each iteration, as shown below. Thus, we rewrite (1.4) as and f (u) is any function of u; its choice will be specified shortly.
The preconditioned ITEM whose convergence conditions are found in [13] is: Here P is the prescribed value of the power and N is a preconditioning operator of the form (2.18) where now c is an arbitrary (unlike in Section 3) positive number.
In what follows we will need a few facts about the linearized method (4.2) [13]. First, the condition that u n , u n = P = const entails the orthogonality between the error and the exact solution: Next, the linearized form of the operator in parentheses in (4.2b) is: The ITEM (4.2) converges when its linearized operator, L, has only negative eigenvalues (see the footnote at the end of Section 3). Loosely speaking (see [13] for a more precise statement), for a wide subclass of equations (1.4), this occurs when operator L has none or one positive eigenvalues. Moreover, in the latter case, the following condition has to hold: Below we will use the change of variables (3.11a) and similar notations: where µ n is given by (4.2a). The linearization of K (0) v n is an operator When the convergence condition of the ITEM (4.2), stated near (4.5), hold, operator K is [13] self-adjoint and negative definite on the space of functions satisfying the exact form of the orthogonality relation (4.3):

A modified CGM
The idea of this modified CGM is to ensure that both the iteration errorṽ n and the search direction be orthogonal, in the sense of (4.8), to the exact solution v. Then, according to the last sentence in the previous subsection, operator K is guaranteed to be negative definite on this space of functions, and hence the following modified CGM is guaranteed to converge under the same conditions as the ITEM (4.2): The definitions ofd 0 in (4.9a) andd n+1 in (4.9f) ensure that d n , N −1 v n = 0 (4.10) at every iteration. Here we use the fact that our iterative solution has a fixed value of power: v n , N −1 v n = P (4.11) for all n. Since the errorṽ n at the nth iteration satisfies the orthogonality condition (4.8), the first of Eqs. (4.9c) and Eq. (4.10) ensure that the errorṽ n+1 at the next iteration also satisfies (4.8). The second of Eqs. (4.9c) does not change this relation; it only forces (4.11) to be satisfied exactly rather than in the linear approximation. Equations (4.9b-d) entail the counterpart of (3.2): 12) and Eqs. (4.9e,f) entail a counterpart of (3.3): We demonstrate these relations in Appendix 1.
Let us note that condition (4.10), in addition to ensuring (4.8) at each iteration, playes one other important role in this algorithm. Namely, it ensures that in the inner product on the left-hand side of (4.13), K is a self-adjoint operator: see the sentence before (4.8). This property of K entails the counterpart of (3.5), which together with negative definiteness of K implies that all search directionsd n are linearly independent and the error decreases as the iterations proceed. In Section 3.1 we noted that this means that the CGM is an optimal iterative method (see, e.g., [19], p. 296, and [24]).
Thus, the modified CGM (4.9) is guaranteed to converge under the same conditions as its Richardson-type counterpart (4.2). This, mathematically, is a more rigorous result than what we obtained for the modified CGM (3.17), which is guaranteed to converge in a slightly narrower parameter space than its Richardson-type counterpart (3.7).
Finally, we present algorithm (4.9) in the original variables: Here L, which is the linearization of is given in (4.4).

Modified CGMs for multi-component solitary waves
Here we present the multi-component versions of the modified CGMs (3.17) and (4.14) for fundamental solitary waves with prescribed values of, respectively, propagation constants and quadratic conserved quantities that are the counterparts of power (1.7). We will not give derivations of these methods because they are similar to those found in Sections 3 and

4.
We begin with stating the generalized Petviashvili method from [7] for finding a multicomponent solitary wave u = [u (1) , . . . , u (S) ] T of the equation whose linearized operator is L: Here N = diag N (1) , . . . , N (S) where each N (k) has a form similar to the linear differential part of the kth component of L (0) ; e.g., for a system of nonlinear Schrödinger-type equations, and c (k) , b (k) are computed by formulas found in [7]. Furthermore, e (k) are the approximate eigenvectors of N −1 L: 4) which are mutually orthogonal with weight N: Note that (5.4) is the multi-component counterpart of (3.8). The form of e (k) 's is related to the solution of (5.1) by e (1) n = u n , e (k) = ρ (k,1) u (1) , . . . , ρ (k,S) u (S) T , (5.6) where ρ (k,k) = 1 and the rest of ρ (k,j) are found as explained in [7]. Finally, the constants λ (k) and γ (k) in (5.2) are found similarly to (3.9) and (2.14): The multi-component version of the modified CGM of Section 3 looks as (3.17), with the scalar functions u n , r n , d n being replaced by the S-dimensional vectors u n , r n , d n , operator N being replaced by its matrix counterpart given before Eq. (5.3), and the operators L (0,mod) and L (mod) being given by expressions generalizing (3.18): A sample code of this algorithm for a two-component system considered in Section 6 is presented in Appendix 3.
Let us note that the purpose of the terms under the sum in (5.8) is to eliminate the eigenmodes e (k) of L whose eigenvalues λ (k) could be positive. By our definition of the fundamental solitary wave, stated at the end of Introduction, eliminating S such eigenmodes should be sufficient to turn the positive eigenvalues of L into negative eigenvalues of L (mod) .
As explained in Section 3, due to the approximate character of relations (5.4), this unstable mode elimination still may not always guarantee convergence of the multi-component modified CGM (3.17), (5.8) when L has small negative eigenvalues. However, a simple approach based on shifting the initial condition, as discussed in Section 6, overcomes this divergence.
We now turn to the case where a set of quadratic conserved quantities related to the solution components' powers are prescribed. We need to introduce a number of new notations. Let us denote the s-component vector of these quantities by Q, so that its kth component is q (kl) P (l) , P (l) = u (l) , u (l) , k = 1, . . . , s ≤ S , l = 1, . . . , S. (5.9) Note that the number of these conserved quantities can be less than the number of the components of the solitary wave: s ≤ S. To emphasize this fact, we use a different vector notation for Q than for u. The number of the components of Q equals the number of the propagation constants that one can prescribe independently. One example of this situation is the system of three waves interacting via quadratic nonlinearity [27] (see also [8]). Another example, which we discuss in more detail below to make our notations clear, is the system describing the evolution of pulses in a two-core nonlinear directional coupler [26]: where U (k) and U (k+2) are the orthogonal polarization components of the pulse in the kth core. Upon the substitution where u (k) can be chosen to be real-valued, system (5.10) reduces to: where µ = µ (1) , µ (2) T . Thus, in this example, S = 4, s = 2, and the vector of quadratic conserved quantities (verified from the time-dependent equations (5.10)) is Generalizing the above example, one can see that the counterpart of Eq. (5.1) when Q rather than µ is prescribed, is the following extension of Eq. (4.1a): For example, in (5.12), L (00) u is the first term (the 4 × 1 vector), and U is the first factor of the second term (the 4 × 2 matrix) on the left-hand side.
For the convenience of subsequent notations, we will assume that the matrix q (kl) in (5.9) is put into the reduced echelon form (see any textbook on undergraduate linear algebra) and, in addition, its columns are rearranged so that For example, for the 2 × 4 matrix in (5.13) this would imply switching the second and third columns. Then for the powers of the solitary wave components whose indices equal the indices of the pivot columns of q (kl) , one has from (5.9): We will use this fact in the next paragraph.
With the convention (5.15), the multi-component version of the ITEM is: In analogy to the algorithm (4.2) for single-component equations, Eq. (5.17b) does not change the linearized form of (5.17a); its role is to guarantee that the s components of vector Q equal their prescribed values exactly rather than in the linear approximation.
To our knowledge, the multi-component ITEM (5.17) has not been presented in the literature before; however, its close "relative" (related to a family of squared-operator methods) was stated in [8]. The convergence conditions for (5.17) are [22]: (i) the s × s Jacobian The corresponding multi-component version of the modified CGM (4.14) is: , k = 1, . . . , s ≤ S , This method is guaranteed to converge to a fundamental solitary wave with a prescribed vector of conserved quantities Q under the conditions stated before Eq. (5.18).

Numerical examples
Here we compare the performance of the modified versions of the CGM presented in Sections 3-5 with the performance of the corresponding Richardson-type methods accelerated by the mode elimination (ME) technique [9,8].
The model equations are (1.5) and its two-component extension: where the constants Codes seeking solitary waves with prescribed power are also much simpler than the code in Appendix 3 below, because they do not compute quantities N, e (k) , etc.) We begin by describing the simulations setup for the methods where the propagation constant(s) is (are) specified. We compare the modified CGMs (3.17) (for the single equation (1.5)) and (3.17), (5.8) (for the two-component system (6.1)) with the ME-accelerated generalized Petviashvili method. The latter method for a single equation is [9]: which for a multi-component system extends to: Here the notations are as in (5.2) above and in addition, where h is the fraction of the slowest-decaying mode that is subtracted at each iteration.
In [9] we advocated subtracting around 70% of such a mode for robust and nearly optimal performance; accordingly, here we use h = 0.7 in all examples.
For all the methods listed in the previous paragraph, we designate the following sets of parameters in Eqs. (1.5) and ( For both the modified CGMs and the ME-accelerated generalized Petviashvili methods, we start the iterations by the non-accelerated generalized Petviashvili method, and when the error, defined as In regards to the specific numeric value on the right-hand side of (6.6), we observed that the accelerated methods converge in all cases when this value is not too high (say, is less than 10 −1 ). If the error is substantially higher, the computed parameters of N etc. may be too inaccurate, so that the corresponding operator K (0,mod) would be "too far" from self-adjoint, which might then prevent the convergence of algorithms (3.17) or (5.8).
In addition, only for the modified CGM in the stiffest case, the value of the acceleration threshold should not be too low (e.g., should be higher than about 5 · 10 −3 ); otherwise, the modified CGM in this case would diverge. The explanation of this fact was sketched in the last paragraph of Section 3.1. Here we present specific details for it; for the sake of convenience, we do it for the single Eq. (1.5) and use the transformed variables (3.11), (3.16).
First, as is well known (see, e.g., [28] and references therein), when µ in (1.5) is very close to the edge of the spectral gap of the linear part of K (0) (i.e. in the "stiffest" case considered here), the shape of the exact solution of (1.5) is very similar to that of the eigenmode(s) of K with negative near-zero eigenvalue(s). In particular, the solution is broad, as illustrated in the top panel of Fig. 2. Next, from (3.17a,b) and an analog of (2.3) it follows that at the first CGM iteration, the step size α 0 along the search direction d 0 is (Here the iteration's number is referenced from the start of the CGM algorithm.) The last strong inequality follows from two facts: (i) operator K (= N −1/2 LN −1/2 ) has a band of negative eigenvalues that comes very close to zero (see, e.g., Fig. 4(a) in [21] for a related problem), and (ii) it is precisely a superposition of the corresponding eigenmodes that dominates the errorṽ 0 obtained by the generalized Petviashvili method after sufficiently many iterations (because these modes decay the slowest in Richardson-type methods; see, e.g., [8]). Having α 0 ≫ 1 and hence a large first step in the CGM adds to v 0 some superposition of near-zero eigenmodes of K and thereby makes the solution at the next iteration, v 1 , more broad and hence even more resembling the near-zero eigenmodes. Then the step size α 1 , computed at the next iteration, becomes even greater than α 0 , and the solution at the following iteration, v 2 , is made even more broad. This quickly leads to divergence of the iterations. (Practically, the solution converges to a nonlocalized two-dimensional plane wave.) On the other hand, if the acceleration by the modified CGM starts when the error v 0 is not too small and hence still contains a significant contribution from the eigenmodes whose eigenvalues are not too close to zero, the corresponding step size α 0 is not too large, and subsequent CGM iterations are able to gradually reduce the errorṽ n .
The above consideration explains why the CGM iterations should start when the iteration error is not too low. In practice, however, this is not a limitation, since in numerically stiff cases, one wants to start the acceleration by the CGM sooner rather than later since this considerably reduces the computational time, as we will demonstrate below.
The initial condition in all cases that we consider is taken as for Eq. (1.5) and as for Eqs. (6.1). The exact values of the amplitude(s) and width of the initial condition are not essential; all methods converge for a wide range of these parameters. As for the asymmetry parameters ǫ x and ǫ y , these can be also varied quite a bit as long as the shape if the initial condition resembles a pulse with one main peak. However, these parameters should not be set to zero (or too close to zero) for the modified CGMs (3.17) and (5.8) in the stiffest case. (That is, in the other two, less stiff, cases, the modified CGMs converge even for a symmetric initial condition.) The reason for this is similar to the reason why the acceleration threshold in (6.6) should not be chosen too small. Namely, iterations of the generalized Petviashvili method which start with a symmetric initial condition, i.e., (6.7) with ǫ x = ǫ y = 0, produce symmetric solutions at all subsequent iterations. These solutions increasingly resemble the near-zero eigenmode of the linearized operator L. Then when the CGM starts, the first step α 0 along the search direction d 0 becomes too large, which leads to divergence of the iterations, as explained after Eq. (6.7). On the other hand, having an asymmetric initial condition leads to the iteration error being asymmetric, and such an asymmetric error can have a considerable content of eigenmodes of L whose eigenvalues are not close to zero. This reduces the step size α 0 at the first CGM iteration, and the modified CGM is able to converge. Note that from the practical perspective, starting with an asymmetric initial condition is not a limitation of the method.
While the ME-based acceleration does not require an asymmetric initial condition for convergence, we still use the same expressions (6.8) for all methods, so as to make our performance comparison uniform.
We now comment on the choices of the "fictitious time" step size ∆τ . For each simulation, we first selected its nearly optimal value for the corresponding non-accelerated generalized Petviashvili method. This takes just a couple of trials since the optimal ∆τ is only slightly less than the maximum value of this parameter for which the generalized Petviashvili method still converges. (See, e.g., Eqs. (2.5) and (2.6) in [9] and Fig. 1(d) in [8], based on the same equations, although for another iterative method.) Once this ∆τ opt has been determined, we use a slightly smaller (specifically, by 0.1) value of the step size for the method accelerated by ME; the analysis of ME [9] predicts that this should lead to a more robust and smooth peformance than using ∆τ opt . For example, if we find that ∆τ opt = 0.9 for the non-accelerated generalized Petviashvili method, we then use ∆τ = 0.8 for the ME-accelerated version of this method. For the method accelerated by a modified

Iterative method
Case in (6.4b) non-accelerated (5.2) accel. by ME (6.2b) accel. by CGM (  CGM, we also start the iterations using ∆τ = ∆τ opt − 0.1, although this does not noticeably affect its performance (since the CGM itself does not involve ∆τ ). for each of the three methods (non-accelerated generalized Petviashvili method and its two versions accelerated by ME and CGM) in the three cases (6.4) for Eqs. (1.5) and (6.1). In accordance to the note above, the value of ∆τ is listed only for the non-accelerated method.
In Fig. 3 we plot the iteration error (6.5) versus the iteration number for the stiffest case for Eq. (1.5). The error evolutions shown there are representative of all the other cases, except that in the less stiff cases, the curves corresponding to the ME and CGM are less jagged.
We also verified that a straightforward extension of the CGM, as described at the end of Section 3.1, would diverge for all cases listed above and even for non-stiff cases (not listed here).
From these tables we see that, as expected, both ME-and CGM-based accelerations considerably improve the performance of the iterative method, with the stiffer the case, the more the improvement. Furthermore, the CGM-based acceleration performs better than the ME-based one by a factor of 2-2.5; again, the stiffer the case is, the more improvement the CGM provides over the ME.
We now present similar results when the power (1.7), or its multi-component counterpart (5.9), of the solution is specified. The non-accelerated methods are the ITEMs (4.2) and (5.17) applied to the single-and multi-component Eqs. (1.5) and (6.1), respectively. When the ME-based acceleration is applied, the lines updating the intermediate solutionsû n+1 andû n+1 of these methods become: where γ slow n and Φ slow n are defined as in (6.3) and L (0) u n is the term in parentheses in (6.9b); the last lines remain the same as in the respective methods (4.2) and (5.17). The modified CGMs for the single-and two-component equations are (4.14) and (5.20).
The three cases of increasing numerical stiffness are chosen as, For ( where the values of µ, as computed within the methods, are listed for reference only. Note that for the single-component equation, we had to choose the parameters for which the solution would satisfy the stability condition (4.5) (see Fig. 1) of the iterative methods (4.2), (6.9a), and (4.14). For the two-component equation, the parameters were chosen by trial and error.
The preconditioning operator was taken to be of the form (2.18) or, in the two-component case, as with c = 1 (which did not necessarily result in the optimal convergence rate). The "fictitious time" step ∆τ for the ME-based methods was chosen as ∆τ = ∆τ opt − 0.1, as described above, where ∆τ opt was the optimal step size for the non-accelerated ITEM. The computational domain and the initial conditions were taken as for the methods that find solutions with prescribed values of the propagation constant. Note that here, the modified CGMs (4.14) and (5.20) are guaranteed to converge to the solution for any initial conditions that resemble a two-dimensional pulse. Nonetheless we used (6.8) with ǫ x = 0.1 and ǫ y = −0.2

Iterative method
Case in (6.10a) non-accelerated (4.2) accel. by ME (6.9a) accel. by CGM (4.14) mildly stiff, ∆τ = 0.   just for uniformity of all our simulations. For the same reason, we also started the acceleration at the same threshold (6.6), even though both the ME-based methods (6.9) and the modified CGMs (4.14) and (5.20) could be employed at the very first iteration. The results are summarized in Tables 3 and 4.
Again, as expected, we see that both ME-and CGM-based accelerations considerably improve the performance of the iterative method, with the stiffer the case, the more the improvement. Furthermore, the CGM-based acceleration performs better than the MEbased one by a factor of 2-3 for the two-component equation, with the stiffer the case, the more imrovement being provided by the CGM over the ME. Interestingly, however, for the single-component equation, the modified CGM provides an improvement over the ME only for the stiffest case; for somewhat less stiff cases, the improvement (in terms of computing time) is only marginal.

Conclusions
In this work, we proposed modifications of the well-known Conjugate Gradient method The modified CGMs that we proposed in Sections 3 to 5 do not, strictly speaking, meet that goal. However, we show below that they come as closely as theoretically possible to doing so. As far as their practical applications, we demonstrated that our modified CGMs converge in the same parameter regions as earlier proposed, slower iterative methods. This convergence can be achieved by a simple shift of the initial condition, which has a theoretical explanation given in Section 6.
When finding solitary waves with prescribed values of the propagation constants, we consider an equivalent nonlinear equation whose linearized operator is modified in such a way that its only positive eigenvalue is essentially turned into a negative one; see Eqs. (3.18) and (5.8). This, however, makes the modified linearized operator "slightly" non-self-adjoint, because the generalized eigenfunction(s) of the original linearized operator, as in (3.12) and (5.4), is (are) available only approximately. This is a fundamental fact about all nonlinear equations except for their very narrow subclass, equations with power-law nonlinearity (1.1) (see also their two-component generalization in [7]), and hence it, in general, cannot be improved. Thus, this fact causes "slight" non-self-adjointness of the modified linearized operator. This, in turn, leads to sign indefiniteness of the quadratic form in (3.17b,e) and thereby prevents one from obtaining a rigoroius guarantee that the modified CGM would always converge. However, we explained in Section 6 and Appendix 1 that the method can diverge only when L has small eigenvalues, and for those cases pointed out that a simple shift of the initial condition will suffice to make the modified CGM to actually converge.
For finding solitary waves with prescribed values of the power (1.7) or, more generally, quadratic conserved quantities (5.9), the situation is different. There, the modified CGM is caused to converge not by modifying a linearized operator but by making the search directions satisfy a certain orthogonality relation; see (4.3) and (5.19). The linearized operator employed by the method can be shown [13] to be self-adjoint in the space of functions satisfying those orthogonality relations, but it is negative definite only under conditions stated before Eqs. (4.5) and (5.18). This restriction is intrinsic to a method seeking a solution with a prescribed power, and hence cannot be relaxed for any nonlinear wave equations.
Thus, we have justified why the modified CGMs proposed in this work come as closely as theoretically possible to guaranteeing that they would converge to fundamental solitary waves. In practice, however, these methods can always be forced to converge, as have been mentioned above and demonstrated in Section 6.
The modified CGMs are faster not only than Richardson-type methods, but also than those methods accelerated by the slowest-decaying mode elimination (ME) technique [9,8].
Namely, in comparison to the ME-accelerated methods, the modified CGMs provide an improvement of about a factor of two to three in terms of computing time; see Tables 1-4 in Section 6 for details. Importantly, the CGMs do so in the most numerically stiff cases, when the respective non-accelerated methods would converge extremely slowly and hence the acceleration would be most desirable. In such cases, it would pay off to use the modified CGMs instead of ME, even though the latter is a little easier to program into a code. On the other hand, in non-stiff cases, i.e. when the non-accelerated methods would converge in just a few tens of iterations, both the modified CGMs and ME would provide only modest improvement in computing time; compare (2.15) and (2.21). In such cases, it may be simpler to use ME or even the non-accelerated method. Let us point out one other, practical, advantage of the ME-based acceleration over the CGMs. Namely, the fact that we were able to construct modified CGMs is critically grounded in the existence of relations (3.18), (5.8) or (4.3), (5.19), as we explained above. On the contrary, ME can be used to accelerate any Richardson-type iterative method, e.g., methods from [3]- [6] or the ITEM with amplitude normalization [13]. Of course, while these methods do converge in many cases, their convergence conditions are not known, and hence their use with or without ME-based acceleration would come without a guarantee that they would converge.
Finally, let us briefly mention alternative iterative methods, other than Newton's method proper, which can be used to compute solitary waves, both fundamental and non-fundamental.
In a recent numerical study [29], Yang showed that a certain combination of the CGM and  (2.8) and (2.9)) and apply the original CGM to it. For linear systems, this trick has long been known (see, e.g., [19], p. 304), and some researchers used it for finding non-fundamental solitary waves [30]. Let us note, however, that "squaring" an equation leads to squaring the condition number of the corresponding linearized operator (since cond(A 2 ) = (cond(A)) 2 ), which, according to (2.21), slows down the convergence; also, more arithmetic operatons per iteration are required for a "squared" method.

Acknowledgement
This work at its early stage was supported in part by the National Science Foundation under grant DMS-0507429. The author thanks Dr. J. Yang for useful discussions and Mr.
N. Jones for assistance with computations at the early stage of this work.

Appendix 1: Technical results for Sections 3 and 4
Here we will first prove statements (i)-(iii) found in the first paragraph of Section 3.3 and also explain the choice (3.15b) for the constant Γ. Then we will prove that the fact that L may have a zero eigenvalue due to a translational eigenmode will not cause the CGM to break down. Finally, we will supply missing steps in the derivations of (4.12) and (4.13).
Statement (ii), i.e. that K (mod) is approximately self-adjoint, follows by considering the following inner product for arbitrary functions f and g: The two approximate equalities above are due to the approximate relation (3.12), and we have also used the fact that K is self-adjoint.
Statement (iii), i.e. that all eigenvalues of K (mod) that are not too close to zero are negative, follows by analogy with the statement found after Eq. (2.13b). First, due to (3.12), one eigenfunction, ψ (1) , of K (mod) is close to v, and for it, When Γ is chosen according to (3.15b), the eigenvalue corresponding to ψ (1) is approximately −1, which is near the accumulation point of the spectrum of K [13]. It is convenient to "place" this eigenvalue inside (as opposed to at the edge of) the spectrum of K (mod) because then it does not affect the condition number of the operator and hence the convergence factor of the iterative method; see (2.16), (2.15), (2.21). Thus, we have shown that the "main culprit" of sign indefiniteness of K -the eigenvalue λ (1) -has been successfully dealt with. To finish the proof of statement (iii), let us show that K (mod) can, in principle, acquire small positive eigenvalue(s). Let ψ (k) , k ≥ 2 be eigenfunctions of K with negative eigenvalues.
Since v is only an approximate eigenfunction of the self-adjoint operator K, then v, Kψ (k) = Kv, ψ (k) does not exactly equal to zero, which causes the eigenfunctions, and hence eigenvalues, of K (mod) to differ slightly from those of K. Hence small negative eigenvalues of K could, in theory, become small positive eigenvalues of K (mod) .
This, however, does not appear to be the case in practice, since if it had been, then the generalized Petviashvili method (3.7) would diverge (see the text before Eq. (3.6)), whereas our numerical experiments conducted in [7] and in this paper do not show such a divergence.
Even if K (mod) has only negative eigenvalues, the quadratic form d n , K (mod)d n may still not be negative definite because K (mod) is not self-adjoint. Since, however, K (mod) is only "slightly" non-self-adjoint, the sign indefiniteness of the above quadratic form can occur only when some eigenvalues of K (mod) are close to zero, as illustrated by the following 2 × 2 example: Note that, for ǫ ≪ 1, the vector in (A1.5) is aligned primarily with the eigenvector, (5ǫ 1) T , corresponding to the smaller eigenvalue of the matrix. The practical implication of the sign indefiniteness of the quadratic form d n , K (mod)d n is that it, and hence the denominators in (3.17b,e), can become arbitrarily close to zero. This may cause divergence of the modified CGM (3.17).
To conclude the consideration of statement (iii), note that the eigenvalues of K = N −1/2 LN −1/2 are the same as those of N −1 L. According to the Sylvester law of inertia [25], eigenvalues of N −1 L and L have the same signs. Moreover, for reasonable choices of N (i.e., for c = O(1) in (2.18)), the eigenvalues of both operators are of the same order of magnitude. Thus, small eigenvalues of K (and hence of K (mod) , see above) imply small eigenvalues of L.
Let us now show that if operator L has an eigenmode ψ trans corresponding to translational invariance of the solitary wave, this would not cause a breakdown of algorithm (3.17).
As before, we will perform the analysis in transformed variables (3.11) and (3.16) and in the linear approximation with respect toṽ n ,r n , andd n .
Finally, let us show that (4.12) and (4.13) hold. Since, as we noted in Section 4.2, the second equation in (4.9c) does not change the linearized form of the first equation, we can usev n+1 instead of v n+1 in (4.9d). Then r n+1 =r n + α n Kd n + O(ṽ 2 n ) . (A1.11) Along with (4.9b), this implies (4.12). Next, the only step that requires some explanation in the derivation of (4.13) is the omission of d n , N −1 v n+1 . By (4.10), this inner product is d n , N −1 v n+1 = d n , N −1 v n + O(ṽ 2 n ) = O(ṽ 2 n ) , (A1.12) so that its omission does not change the accuracy implied by the right-hand side of (4.13).
Step 3 implements Eqs. (4.5) of [7]. Note that α k of [7] is denoted as λ (k) in this paper. Also note that in the code, quantities like E_NE(1) denote e (1) , N e (1) , i.e. here, the superscript refers to the particular eigenvector e (k) of N −1 L.
At the first iteration of the CGM, in addition to implementing Eqs. (3.17a) of this paper, one also computes Γ (k) and Ne (k) , k = 1, 2, so that these quantities are not computed again at subsequent iterations of the CGM. Notations LL0u and LLd stand for L (0,mod) u n and L (mod) d n . The remainder of the code implements Eqs. (3.17) and (5.8) of this paper.
Note that the CGM iterations start when N, e (k) , and Γ (k) are found rather imprecisely: accuracy (i.e., the value of normDu_accel) of 5% is used in the above code and in the examples reported in Section 6, and we verified that the code still worked when this accuracy was lowered to as much as 10%. Moreover, even if N etc. were computed up to a higher accuracy, the functions e (k) could still satisfy the eigenrelations (5.4) only approximately.
( u_SumjLkjEj1*sum(u_DELu) -DELu_SumjLkjEj1*sum(u_u) ); u_Nu=c*u_u-u_DELu; % <u, Nu> componentwise Note that the reason why the inner products u_u(k) and u_DELu(k) are computed for each component of u n rather than for the entire vector function, as DELu_DELu etc, is that they are needed to compute the individual components of u_Nu. The latter are needed to compute ρ (2,1) in Step 2.
We did not test the performance of such a simplified code because the goal of this paper is the proposal of modified CGMs for solitary waves and not optimization of those methods.  Figure 3: Evolutions of the iteration error, defined by (6.5) with S = 1, for the stiffest case (6.4a) for Eq. (1.5). The curve for the non-accelerated generalized Petviashvili method (3.7) is not shown in full so as not to obscure the details of the other two curves.