Soliton Phase Shift Calculation for the Korteweg–De Vries Equation

Several non-linear fluid mechanical processes, such as wave propagation in shallow water, are known to generate solitons: localized waves of translation. Solitons are often hidden in a wave packet at the beginning and only reveal themselves in the far-field. With a special signal processing technique known as the non-linear Fourier transform (NFT), solitons can be detected and characterized before they emerge. In this paper, we present a new algorithm aimed at computing the phase shift of solitons in processes governed by the Korteweg–de Vries (KdV) equation. In numerical examples, the new algorithm is found to perform reliably even in cases where existing algorithms break down.


I. INTRODUCTION
Solitons are wave packets that retain their shape while propagating at a constant celerity. They appear in diverse fields of modern physics, see e.g. [1]- [9]. In physical processes governed by partial differential equations, solitons typically form when a dispersive effect, which broadens the wave packet, is compensated by a suitable non-linear steepening effect. The prototypical example is the Korteweg-de Vries (KdV) equation, which has been used to model a plethora of physical processes such as internal waves in the ocean, Rosby waves in the atmosphere, plasma waves, acoustic waves in bubbly liquids [10]; blood pressure waves in arteries [11], [12]; acoustic waves in metal [13]; and traffic flow [14].
As an illustrative example we consider the application of the KdV equation to long unidirectional surface waves in shallow 1 water [16]. A dimensional form of the KdV is where η = η(x,t) [m] is the free surface elevation as a function of the timet [s] and the position inx [m] in lab coordinates and where subscripts denote partial derivatives.
The associate editor coordinating the review of this manuscript and approving it for publication was Jing Liang. 1 Shallow: The depth is less than 0.22 times the wavelength [15].
In the case of free surface waves, the coefficients of (1) are c = gh , α = 3c 2h , β = ch 2 6 , where g [m/s 2 ] is the gravitational acceleration, h [m] the still water depth, and c [m/s] the phase celerity of the wave [15], [17]. By the changes of variables [18] we obtain a standard form of the KdV equation, to wit All variables in this standard form are unit-less and real.
The position x in (2) is expressed with respect to a frame that moves with the phase celerity c, which cancels one term compared to (1). As said, the KdV equation allows soliton solutions, which are waves of translation. In the case of the KdV equation they take the shape of a squared hyperbolic secant, sech 2 (θ ) = 1 2 e −θ + A noteworthy fact is that the KdV equation evolves any wave packet into a parade of N ≥ 0 solitons. (If there are no solitons, N = 0, the wave is fully dispersive and will vanish over time [19, §1.7.c].) That is, after a sufficiently long time t, q(x, t) ≈ N n=1 q n sech 2 (k n x − ω n t − ϕ n ) , (4) where q n = 2k 2 n and ω n = 4k 3 n [20,Eq. 17], [21,p. 83,Eq. 3.3], [22,Eq. 2.20a]. The amplitude, wavelength and celerity of a soliton are thus coupled. We call the regime where (4) holds the far field. Otherwise we are speaking of the near field. The near field can in general not be described as a linear superposition of wave components, due to the non-linear interaction between them.
Nevertheless, from near field data we can obtain the free parameters that describe the far field -the generalized wave numbers k n and the phase shifts ϕ n -long before individual solitons start to separate from a wave packet. We do that with the so-called scattering transform. It is a generalization of the common Fourier transform and therefore nowadays often called Non-linear Fourier Transform (NFT) in the literature [23]. The generalized wave numbers and phase shifts of the solitons are represented in the NFT spectrum by the so-called eigenvalues ζ 0n and norming constants b(ζ 0n ), as will be explained later, and the NFT enables us to calculate them from the normalized free surface at any fixed time. For brevity and to conform to the parlance of the NFT, we will hereafter call the normalized free surface at a fixed time, q(x, t 0 ), the potential. It has been shown that the NFT can extract features of water waves that remain hidden with linear methods of signal analysis. We remark that most of the recent research assumes periodic waves [15], [23]- [26], whereas in this paper we consider wave packets [17], [18], [27]- [29].
The main motivation of this paper is that whereas the numerical computation of the soliton amplitude, wavelength and celerity from near field data is a solved problem, the numerical computation of the phase shifts ϕ n that appear in (4) is surprisingly hard. The difficulty lies in the numerical computation of the aforementioned norming constants b(ζ 0n ). This is the problem we deal with in this paper. Recently two closely related algorithms have been proposed independently to address a similar problem, namely for the NFT that solves the Non-linear Schrödinger Equation (NSE) [30], [31]. In this paper we derive a new algorithm to compute the norming constants of the KdV-NFT that builds on these ideas. The key contribution we make is that we utilize not one but two different estimators of the norming constant. Our new approach is considerably more reliable than the existing approaches. Furthermore, to the best of our knowledge this is the first paper to address the numerical calculation of the norming constants for the KdV equation. 2 The paper is organized as follows. In Section II we recapitulate the relevant parts of the NFT for the KdV equation. In Section III we present and motivate our new algorithm to calculate the norming constants; we validate it with numerical examples in Section IV and show that it is significantly more reliable than a reapplication of the ideas of [30], [31] for the KdV equation. The main part of this paper ends with a conclusion in Section V. We have furthermore included several appendices with supporting information.

A. NOTATION
In this paper we typeset variables in italic and constants upright. The constants include j as the imaginary unit, e as Euler's number and π as the ratio of a circle's circumference to its diameter. For vectors we use lower case bold (v); for matrices upper case bold (A) and I denotes the 2 × 2 identity matrix. Scalars are typeset in lower or upper case regular (q, X 0 ), where for elements of a matrix two subscripts denote the row and column in that order (A ij ). For operators we use a sans-serif upright font (V). A hat on top of a variable (q) means an estimation or approximation of that variable. The set R is the set of real numbers; I is the set of imaginary numbers; C is the set of complex numbers; other sets are denoted with a capital in calligraphic font (X ). The symbol := means that the left-hand side is defined by the right-hand side; the symbol ∼ means 'is proportional to'. We use O as the 'big-O' Landau order symbol. We reserve log for the base 10 logarithm and use ln for the natural logarithm, with base e. The function exp denotes the exponential with base e, so exp(θ) := e θ .

B. NON-LINEAR FOURIER TRANSFORM FOR WAVE PACKETS
In this subsection, we survey the mathematical background of the NFT for wave packets. By wave packets we mean a localized real potential q(x, t 0 ) with so-called vanishing boundary conditions. Formally, 3 for some fixed time t 0 . In this paper we deal with the case that the function q(x, t) furthermore satisfies the KdV equation (2). As a magic step (refer to [32], [33] for the explanation) we use q(x, t) as a time-varying potential in the Schrödinger eigenvalue problem: Then (6) has two important types of solutions. Firstly, for all ζ ∈ R\{0} there exist solutions for which the eigenfunctions f are power signals, 4 and they comprise the so-called 3 If the KdV equation is normalized otherwise, such that q(x, t 0 ) →h as |x| → ∞ for some finite constanth, then (2) is equivalent to q t + 6q q x + q x x x = 0, where x := x − 6ht and q := q −h → 0 as |x | → ∞. If the transformed potential q satisfies (5), we can proceed with this scaled KdV equation continuous spectrum. This part of the NFT spectrum decays over time; the corresponding wave components dwindle as a dispersive wave train [22]. The continuous spectrum is thus of interest for the near-field of a wave, but not in the scope of this paper. Secondly, there is a finite number N ≥ 0 of distinct solutions for which the eigenfunctions f are energy signals. 5 These solutions comprise the so-called discrete spectrum. The values of ζ that are part of the discrete spectrum are called eigenvalues and we will denote them as ζ 0 , or when we refer to a specific eigenvalue as ζ 01 , ζ 02 , etcetera. It can be shown that all eigenvalues of (6) are non-zero imaginary numbers in the upper half plane [33, p. 251], [28,Eq. 3.21], [34, §3]. When we write hereafter ζ , we mean a non-zero number that is either real or imaginary unless explicitly indicated otherwise.
It is the discrete spectrum that has our focus in this paper, for the reason that every eigenvalue corresponds to one separated soliton in the far field. The discrete spectrum can be subdivided into two parts. The first part consists of the eigenvalues, which can be shown to remain constant while the wave evolves [35]. That is, from the near field up to and including the far field. Specifically, in (4), jk n = ζ 0n , so an eigenvalue tells us the amplitude, wavelength and celerity of the corresponding soliton in the far field and as the solitons evolve further into the far field we can observe that those indeed remain constant. However, we can calculate the eigenvalues with the NFT at any fixed t = t 0 .
The second part of the discrete spectrum contains information about the eigenfunctions. This part of the spectrum does evolve over time -in an easy to compute way -and contains the additional information we need, to calculate the phase shifts ϕ n in (4). By the virtue of (5), as |x| → ∞ the Schrödinger eigenvalue problem (6) reduces to and hence every eigenfunction of (6) reduces in this limit to a linear combination of the functions exp(±jζ x) [36, §2.8].
In particular, we can look for one set of eigenfunctions that satisfies and another set of eigenfunctions that satisfies These sets of eigenfunctions are known as the Jost solutions and it should be noted that they are also uniquely defined for all other x by (6). Each of the two sets of Jost solutions forms a linearly independent basis for the eigenfunctions of (6) and they are related as where where the scattering parameters b(ζ, t) ≡b(−ζ, t) and a(ζ ) ≡ā(−ζ ) are implicitly defined by (10); their explicit definitions can be found in Appendix A. The parameter a(ζ ) depends neither on x nor on t; the parameter b(ζ, t) does not depend on x, but evolves over time as [33,Eq. 3 Although the scattering matrix S(ζ, t) is in general complexvalued, it is real-valued when ζ is an imaginary number, so for the discrete spectrum in particular S(ζ 0 , t) ∈ R 2×2 . Another important property of the scattering matrix is [ It can be shown that the set of eigenvalues of (6) Associated to each eigenvalue ζ 0n is a norming constant that is given by b(ζ 0n , t), which is in fact not a constant, but a quantity that evolves over time according to (12). With this norming constant we can finally calculate the phase shifts ϕ n in (4) as [20,Eq. 25a where |ζ 01 | < |ζ 02 | < . . . < |ζ 0N | is required as the order of the eigenvalues and where a (ζ 0n ) := da(ζ ) dζ ζ =ζ 0n ∈ I .
Among the quantities required in (15), the norming constants are especially hard to calculate numerically. 6 We will substantiate this claim in Section III. The main goal of this paper is to still calculate these norming constants accurately, to be able to calculate the phase shifts. Not only do the norming constants contain information about the position of the solitons in the far field, but also of the evolution of the wave in the near field. The calculation of the latter is however more involved due to the non-linear interaction of the wave components. In the general case one needs to compute the inverse NFT, e.g. [34]; for a pure soliton potential (i.e. if b(ζ, t 0 ) = 0 ∀ζ ∈ R\{0}) that may contain significant non-linear interactions between the solitons, simpler methods exist, such as the one described in Appendix F and the reference therein.
In the remainder of this paper we drop the dependence on t = t 0 from our notation.
Remark 1 (Sign Inconsistency in [33]): The definitions used in this paper for bothφ(x, ζ ) andb(ζ ) differ from [33,§III] by a minus sign, in order to be consistent with [33,Appendix 3]. Appendix E, where we talk about the Nonlinear Schrödinger Equation (NSE), is an exception.

C. CALCULATION OF THE SCATTERING MATRIX
The mathematical definition of the NFT is not convenient for numerical calculations. Therefore the Schrödinger eigenvalue problem, (6), is typically rewritten as a system of first order differential equations, where and V(ζ ) is a column vector of length two of operators.
In literature different choices are made for V(ζ ). We say that each suitable choice results in an analytically equivalent calculation, expressed in a different basis for the state vector v(x, ζ ). In this paper we make the special choice that simplifies the exposition. We will refer to this choice as the S basis. In Appendix C we show how our results extend to alternative choices found in literature. For the Jost solutions with respect to the S basis we let by (10). Then we left-multiply (21) and (22) by and take the limits for x → ∓∞ respectively to obtain where we used (8) and (9). It is useful to note that If we are dealing with a potential that satisfies (5) and is furthermore zero outside some window (X − , X + ), the boundary conditions defining the Jost solutions [(8) and (9)] hold for all x ≥ X + or all x ≤ X − respectively. This allows us to replace x → −∞ in all the previous equations by x = X − and x → ∞ by x = X + . If we then multiply (25) from the left by E(−ζ X − ) and from the right by (24) we find where That is, H S (X − , X + , ζ ) is a state transition matrix of (17), a matrix that defines a bijective linear mapping from every initial state vector v S (X − , ζ ) to its corresponding final state vector v S (X + , ζ ). The calculation of the scattering matrix S(ζ ) can then be done by numerically evaluating (17) to find H S (X − , X + , ζ ), after which the scattering matrix is found according to (27).

Remark 2 (Why Not the Ubiquitous AKNS Basis?):
An alternative choice of basis results in the so-called Ablowitz-Kaup-Newell-Segur (AKNS) system [33]. The AKNS system is popular for two reasons. Firstly it is a framework that covers multiple evolution equations, among which the KdV equation and the NSE. Secondly, for most of these evolution equations, the AKNS system is the special choice of basis that simplifies the exposition, because (24) and (25) (24) and (25) respectively, in the sense that they implicitly use the change of basis matrix from the S basis to the AKNS basis, see (77) in Appendix C. It is hence by the virtue of the S basis that the exposition in this paper parallels the one for e.g. the NSE in the AKNS basis.

III. BIDIRECTIONAL ALGORITHM
Recall from the introduction that the eigenvalues ζ 0 of a potential -a normalized free surface at any fixed timesignify the amplitude, wavelength and celerity of solitons, but contain no information about their phase shifts in the far field or their location in general. Therefore we need for each eigenvalue a second parameter that encodes this information: the norming constant b(ζ 0 ).
It is notoriously hard to calculate norming constants numerically, but the topic appears to be unmentioned in literature for the KdV equation. Let us thus shortly outline the reason for this issue. The eigenfunction f (x, ζ 0 ) should by definition be an energy signal. Therefore f (x, ζ 0 ) must be bounded 7 as x → ±∞. For the KdV equation jζ 0 < 0, so the Jost solutionsφ(x, ζ 0 ) andψ(x, ζ 0 ) are by definition unbounded, see (8) and (9). Hence we can express the eigenfunction as a scalar multiple of the remaining, bounded Jost solutions: by (10) and (14). However, (28) will in general not hold exactly in a numerical calculation. Instead, an eigenfunction is typically represented likê The Jost solutionψ(x, ζ 0 ) grows exponentially as x → ∞, so whenâ(ζ 0 ) is small but not exactly zero,φ(x, ζ 0 ) grows exponentially as x → ∞. Although this is in principle an error inâ(ζ 0 ), it will make the calculation of the norming constant b(ζ 0 ) ill-conditioned, as illustrated in more detail in Appendix D.
The bidirectional algorithm is a numerical method to calculate the norming constants that lessens the influence of the aforementioned numerical errors. It was published in [30] and simultaneously discovered as the forward-backward method in [31]. The key idea is to evaluate the norming constant b(ζ 0 ) not at a boundary of the window (X − , X + ), but at a point x = X 0 in between, and enforceâ(ζ 0 ) = 0 in the numerical calculation. The algorithm was originally developed and presented for use with the Zakharov-Shabat (ZS) system (albeit with a change of variables), which is equivalent to the AKNS system for the NSE. In Section III-A we extend this method to make it usable for the KdV equation and we formulate an overdetermined equation for the calculation of norming constants b(ζ 0 ). We show in Appendix E that the equations found in [30], [31] respectively for the calculation of norming constants for the NSE are actually the two different halves of this overdetermined equation. In Section III-C we propose a new criterion for the selection of the cutting point X 0 that is observed to minimize the estimation error of the norming constantb(ζ 0 ) by utilizing the overdetermined equation as a whole. By minimizing the estimation error of the norming constant, we also minimize the estimation error of the phase shift that is calculated from the norming constant. We summarize the criteria found in literature in Section III-C as well, but (as we will show in Section IV) our criterion leads to considerably more reliable estimates of the norming constants than the existing criteria.

A. BIDIRECTIONAL ALGORITHM FOR THE KDV EQUATION
We assume again that the potential resembles a wave packet of finite length: The potential satisfies (5) and we can choose a window (X − , X + ) ⊂ R such that the potential q(x) is zero outside this window. Let us cut the potential in a point X 0 ∈ (X − , X + ) in a left and a right part: For the state transition matrices from X − to X 0 and from X 0 to X + it holds by (26) that can be calculated from the potentials q L (x) and q R (x) respectively. Then by (27), where S R (ζ ) and S L (ζ ) are the scattering matrices for the respective potentials q R (x) and q L (x). By (13) and Cramer's rule, We multiply (35) for ζ = ζ 0 from the left by S −1 R (ζ 0 ) and from the right by 1 0 T to find after substitution of a(ζ 0 ) = 0 that which is an overdetermined equation from which b(ζ 0 ) is to be solved. If S R (ζ 0 ) and S L (ζ 0 ) are the exact scattering matrices of any potentialq(x) at an exact eigenvalue thereof, (37) will be consistent. Hence, an inconsistency in (37) indicates a numerical error that is not due to e.g. discretisation of the potential. The extent to which (37) is consistent appears to depend heavily on the choice of the cutting point X 0 , an observation we exploit in Section III-C to formulate a new criterion for this choice. To facilitate the discussion thereof, we use the two rows of (37) to define two separate estimators for the norming constant: which depend implicitly on X 0 via (30) and (31). In the derivation above we have not posed any restriction on the way ζ 0 , (the first row of) S R (ζ 0 ), and (the first column of) S L (ζ 0 ) are calculated numerically. The bidirectional algorithm is hence independent on the numerical method by which (17) is solved. This may for example be an exponential integrator method (e.g. [17], [38]- [40]) or a collocation method (e.g. [41, §2.4.3], [42]). Furthermore we are free to choose the basis for this calculation. In the S basis (37) becomes, after left-multiplication by e jζ X − E(−ζ X 0 ), With the aid of Appendix C we could readily express (37) in any other basis we may wish to use. For example in the ubiquitous AKNS basis we obtain after left-multiplying (37) by e jζ X − E(−ζ X 0 ) the equivalent expression which demonstrates that the use of the AKNS basis for the KdV equation is possible at the cost of more complicated equations compared to the S basis.
In Appendix E we link the formulation of (37) to the previous work on the bidirectional algorithm, which was solely aimed at calculations for the NSE in the AKNS basis. In short, both [30] and [31] develop the bidirectional algorithm as described in this paper, but [30] finds only (38), whereas [31] finds only (39) as the estimator for the norming constant. They furthermore use different criteria to select the cutting point, as we will discuss in Section III-C. However, since we are proposing a criterion that aims to minimize the error in the phase shift, we need to discuss first -in the next subsectionhow this error is affected by an error in the norming constant.

B. PHASE SHIFT ERROR
We have seen that we can calculate the phase shifts ϕ n appearing in (4) from the norming constants b(ζ 0n ) with (15). From a numerical estimation of the norming constant,b(ζ 0n ), we can thus calculate an estimation of the phase shift,φ n . It is readily verified that the estimation error in the phase shift satisfieŝ when all other variables that appear in (15) remain the same.
Hence an additive error in the phase shift is directly related to a multiplicative error in the norming constant. We must further consider what happens when this multiplicative error is negative. (We will see in Section IV that this occurs commonly for calculations of the norming constant according to the benchmark algorithms.) In that case (42) evaluates to a complex number, which is a meaningless result. Empirically, when we reconstruct a potential from an NFT spectrum that is modified by flipping the sign of one or more norming constants, we obtain a completely different potential. Hence, a proper measure for the phase shift error is Errors in numerical algorithms for the calculation of the NFT emerge for example due to discretisation of the potential, where there is a trade-off between the error and the required number of computations. For an exponential integrator method the required number of computations depends on the number of samples D and the relative error in the result is typically of the order O(D −p ) for some positive integer p. We remark that then the error measure defined in (43) converges at the same rate, i.e.
which can be shown by Taylor expansion.

C. CHOICE OF CUTTING POINT
We have shown in Section III-A how the bidirectional algorithm can be used for the calculation of the norming constants for the KdV equation. Two questions are left to answer in this subsection: How should we choose the cutting point X 0 and how do we find the optimal estimate of a norming constant from the overdetermined equation (37)  provides the most accurate estimate of the norming constant in every case. By following the proposed criterion the error between the two estimates - (38) and (39) respectively -will generally become negligible compared to the discretisation error, in which case we can simply select any of the two as the numerical norming constant. Suppose that the potential q(x) is known on a grid x ∈ . . , D}}, where the step size ε := (X + −X − )/D. The natural candidates for the cutting point are the points exactly in between: X 0 ∈ X 0 := {X − + mε | m ∈ {1, 2, . . . , D − 1}}. We propose to select the cutting point that both minimizes the relative error between the two estimates of the norming constant (b 1 (ζ 0 ) andb 2 (ζ 0 ) as defined in (38) and (39) respectively) and minimizes the additive error of the phase shifts that can be calculated from these. That is, with E as defined in (43). The intuition behind this choice is that when the numerical error is small, the two estimates are close to each other. Conversely, when either of the two suffers from a large numerical error, the two estimates are probably far apart.
In the literature three other criteria are known to select the cutting point. These were all proposed for the calculation of norming constants with respect to the NSE, but we will use them as benchmarks for the calculation with respect to the KdV equation anyway, as there are to our best knowledge no such criteria for the KdV equation.
In case X − = −X + , two options are given in [31, §III.A]: Shift the potential such that (X − , X + ) shifts to ( X + −X − −2 , X + −X − 2 ) and correct for this space translation afterwards, or pad the potential with zero on one side to enlarge the window until X − = −X + . The choice between these two affects the outcome of (47). We will refer to these criteria as the Aref criterion with potential shift or support extension respectively.
• The 1-norm criterion, used by the Fast Non-linear Fourier Transform (FNFT) software library [43], is Apart from the bidirectional algorithm one could also estimate the norming constant with a naive calculation, directly from the scattering matrix: We treat the naive calculation hereafter in the framework of the bidirectional algorithm by choosing where we have used that the scattering matrix S R (ζ 0 ) of q R (x) ≡ 0 is the identity matrix. Note thatb 1 (ζ 0 ) is undetermined for X 0 naive .

IV. NUMERICAL EXAMPLES
In this section we compare the numerical calculation of the norming constants according to the bidirectional algorithm with the proposed criterion to the benchmark criteria listed in Section III-C. We demonstrate that the bidirectional algorithm with the proposed criterion computes the correct norming constants even in difficult examples where all other criteria fail. In each example we start with a potential q(x) of VOLUME 7, 2019 which the norming constants are analytically known. Then we calculate the norming constants numerically according to our proposed criterion as well as each of the benchmark criteria and compare the results.

A. EXAMPLE SETUP
For each example we approximate the potential q(x) as a piecewise constant functionq(x) that is 0 outside a window x ∈ (X − , X + ). The step size is a constant ε := (X + − X − )/D, where D is the number of samples. Each step has the same value as the potential at the midpoint of the step, soq(x m ) = q(x m ) for all x m = X − + mε − 1 2 ε, where m ∈ {1, 2, . . . , D}. This approximation of the potential introduces a relative error in the spectrum proportional to D −2 [37], [39], thus by (44) we expect the error according to the error measure defined in (43) to be of the order O(D −2 ). 8 For the approximated potentialq(x) we find the eigenvalues numerically according the algorithm described in [37, §4] where we take U(q, ε) = exp(εA S (x, ζ )) [see (20)] to do the calculation in the S basis for the KdV equation. Then for each eigenvalue we calculate the two norming constant estimateŝ b 1 (ζ 0 ) andb 2 (ζ 0 ) according to (38) and (39) respectively at every cutting point candidate X 0 ∈ X 0 . Finally we find the cutting point according to the proposed criterion, (45), as well as to the benchmark criteria, (46) to (48) and (50) respectively. 9 We will report both estimates of the norming constant for each of the cutting point criteria (except for the naive computation), even though their respective sources make use of only one (see Appendix E).
We display the results for each example in two different ways. Firstly, we choose a low number of samples D for which still the analytically known number of eigenvalues can be found. We vary the cutting point X 0 and plot for every eigenvalue against that the error between the two norming constants estimatesb 1 (ζ 0 ) andb 2 (ζ 0 ) as well as the error between both of these estimates and the analytically known norming constant. Secondly, we vary the number of samples D and plot for every eigenvalue against that the error between both of these estimates and the analytically known norming constant, where the cutting point is chosen according to the proposed criterion and each of the benchmark criteria respectively. We plot these cutting points as well against the number of samples.

B. EXAMPLE 1: TWO SEPARATED SOLITONS
For this example we choose a pure soliton potential with two eigenvalues: ζ 01 = 0.5j with norming constant b(ζ 01 ) = −10 10 and ζ 02 = 0.6j with b(ζ 02 ) = 10 −12 . The resulting potential consists of two well separated solitons as shown in Fig. 1 (top), with the soliton corresponding to ζ 01 at x ≈ 25 8 Since a piecewise constant interpolation leads (with a suitable x-grid) to an exact representation of a rectangular potential, i.e.q(x) ≡ q(x), such an example shows atypical results for all criteria and is therefore not included. 9 In case any of the criteria does not have a unique global minimum, we choose the leftmost (lowest) cutting point among the global minima. and the one corresponding to ζ 02 at x ≈ −25. We calculated this potential with the algorithm from [44], with an essential numerical improvement described in Appendix F.
For the numerical calculation of the NFT we approximate this potential (initially) with a coarse grid of D = 100 samples in the interval x ∈ (X − , X + ) = (−50, 50) as shown in Fig. 1 (top). Fig. 1 (middle) shows the error between the two estimatesb 1 (ζ 0 ) andb 2 (ζ 0 ) for all cutting point candidates X 0 ∈ X 0 . (Here and later, errors above a certain threshold are not shown as they can become very large.) These errors attain a minimum when X 0 is at the location of the corresponding soliton. Hence, at this location (37) is most consistent. Away from the soliton this error becomes several orders of magnitude larger. Motivated by this observation we want to know if the numerical error is minimized by choosing the cutting point for each norming constant according the proposed criterion. In Fig. 1 (bottom) we plot therefore the error between both respective estimates and the ground truth and indeed we see that the numerical error of both estimates for both solitons is minimal when we choose the cutting point X 0 according to the proposed criterion. 10 The main difference compared to Fig. 1 middle is that the error reaches an error 10 Sharp dips below the error floor, such as in Fig. 1 (bottom) for E(b 2 , b; ζ 02 ) at X 0 = −40 are cases where the error in the estimate of the norming constant coincidentally cancels the error due to the approximation of the potential. Since the norming constant is a real number this is not unlikely to happen at some cutting point candidate, but there is no way to find that point without knowledge of the true norming constant and exploit this effect.  floor in the vicinity of the soliton, thereby forming bathtub shaped curves rather than V-shaped curves. This error floor is caused by the approximation of the potential by a piecewise constant function and can be lowered by reducing the step size. This is shown in Fig. 2, where we plot (as a function of the number of samples D) X 0 according to each of the criteria as well as the error of the resulting norming constant estimatesb 1 (ζ 0 ) andb 2 (ζ 0 ) compared to the true norming constant b(ζ 0 ) for both of the eigenvalues ζ 0 . It can be seen that the bidirectional algorithm with the proposed criterion is the only one for which the error decay is consistently for all estimates O(D −2 ) -a factor 100 per decade -as expected. The other criteria find norming constants with an error that shows no convergence and that is several orders of magnitude larger for one of the eigenvalues or for one of the estimators,b 1 (ζ 0 ) orb 2 (ζ 0 ). We remark that curves that leave the graphing area vertically indicate that the neighboring data point corresponds to an estimate of the norming constant with the opposite sign, yielding infinite error by (43). The reason why the benchmark criteria perform like this, is that the solitons in the potential are separated from each other. The benchmark criteria select the cutting point for all eigenvalues either near one of the two solitons or right in the middle, which are clearly no suitable cutting points in every case if we look at Fig. 1 (bottom). Hence, even for this simple potential the proposed criterion is the only one that results in only reliable estimates for the norming constants.

C. EXAMPLE 2: SIX PARTIALLY OVERLAPPING SOLITONS
In this example we construct a pure soliton potential with six eigenvalues, ζ 0n = 0.1nj for n ∈ {1, 2, . . . , 6}, and norming constants b(ζ 0n ) = (−1) n exp(8n(−1.01) n ). The resulting potential is calculated as in Example 1 and consists of two clusters of three overlapping solitons each: The solitons for odd n cluster at x ≈ −45, those for even n cluster at x ≈ +45.
For the numerical calculation of the NFT we approximate this potential with a coarse grid of D = 316 steps in the interval x ∈ (−160, 160). The potential and its approximation are shown in Fig. 3 (first plot). The error between the two estimatesb 1 (ζ 0 ) andb 2 ,(ζ 0 ) are shown in Fig. 3 (second plot). Again we see that for all norming constants the global minimum of this error is for X 0 in the vicinity of the corresponding soliton. In Fig. 3 (last two plots) we show for each eigenvalue the error between both respective estimates and the ground truth. Similar to Example 1 we see that the numerical error of both estimates for all solitons is at a minimum when we follow the proposed criterion. Due to the larger window (X − , X + ) we can recognize the bathtub shapes of the curves better than in Example 1, and see that the ones forb 1 (ζ 0 ) are shifted to the left compared to the corresponding soliton, whereas the ones forb 2 (ζ 0 ) are shifted to the right. Furthermore we see that the bottoms of the bathtubs become narrower as the corresponding eigenvalue increases in magnitude, implying that choosing a good cutting point becomes more important as the energy of the corresponding soliton increases.
We can lower the bottoms of the bathtubs, and thereby the achievable error of the estimate of the norming constant, by increasing the number of samples, as shown in Fig. 4. Again we observe that all of the benchmark criteria return several estimates of the norming constant with an error that does not decay consistently with an increase of the number of samples, or is even high throughout. The proposed algorithm in contrast returns estimates of the norming constant of which the error decays neatly at the expected rate of O(D −2 ) for every eigenvalue and for both estimates thereof.
In Fig. 5 (middle) we show the error between the two estimates of each norming constant. We see that all three curves show a global minimum near the axis of symmetry of the potential and in Fig. 5 (bottom) we show the errors between these two estimates and the true norming constants, which have a wide global minimum near the axis of symmetry of the potential. This minimum is quite high because of the limited number of samples.
When we increase the number of samples, we obtain the results shown in Fig. 6. The Hari criterion, the 1-norm criterion, and the proposed criterion select a cutting point near the axis of symmetry of the potential and the errors of all their estimates of the norming constant decay at the expected rate of a factor 100 per decade, O(D −2 ). The flooring that is seen in all these cases from around D = 10 5 samples is because the error due to the truncation to x ∈ (0, 80) becomes dominant compared to the error due to the piecewise constant approximation itself and could hence be removed by enlarging the window (X − , X + ). For this example even the naive computation shows the expected error decay for two of the three norming constants, the third one floors at a higher level. The Aref criterion 12 selects cutting points for two of the eigenvalues at the far-right end of the potential, which results in two estimates of the norming constants with a non-decaying error. From this example we see that in the case of a potential where the solitons are all clustered together, the proposed criterion performs equally well as the best among the benchmark criteria.

V. CONCLUSION
In this paper we investigated how to calculate the phase shifts of solitons in the far field of any wave packet that evolves according to the Korteweg-de Vries (KdV) equation. For that we in particular need to compute the so-called norming constant of each soliton. The naive method to compute norming constants is however known to be unreliable. We adapted the bidirectional algorithm, which was originally designed for the computation of norming constants for the Non-linear Schrödinger Equation (NSE), to the KdV equation. Furthermore we proposed a new criterion to select the cutting point required for this algorithm. The criterion is based on the observation that the bidirectional algorithm actually provides two estimates of each norming constant. The proposed criterion is to select the cutting point that minimizes the difference between these two estimates, which is observed to minimize also their distance to the true norming constant. We demonstrated with three examples that the proposed method performs at least as good as existing algorithms, and often several orders of magnitude better.
With the reliability of our method established, we plan to implement it in the FNFT software library [43] and to apply it to real-world data in the future.

APPENDIX A DEFINITION OF THE SCATTERING PARAMETERS OF THE KDV EQUATION
By Abel's identity [45, p. 22], the Wronskian of any two solutions f 1 and f 2 of (6) for equal ζ is independent of x: Consequently, we may evaluate Wronskians of the Jost solutions, (8) and (9), at any convenient value x: Because these Wronskians are non-zero for ζ = 0, both sets of Jost solutions are linearly independent for ζ = 0. To relate these two sets, define the scattering parameters as 12 If there is a cutting point X 0 for which q(X 0 ) = 0, then it minimizes (47). Consequently the Aref criterion with support extension returns a cutting point in the zero padding of the potential, which does not give up to par results. Therefore we use the Aref criterion with potential shift in this example.

Lemma 1 (Space translation):
If the scattering matrix of a potential q(x) is S(ζ ) as in (11), then the scattering matrix of the potential q (x) ≡ q(x − x 0 ) is where E(ζ x) is defined in (23). Proof: The scattering problem for q (x) is equivalent to the scattering problem for q(x ) with x := x − x 0 . The Jost solutions in the translated coordinate are Equation (59) follows from filling these out in the definitions of the scattering parameters, (55) to (58).

Remark 3 (Change of x-coordinate frame):
Although b(ζ ) andb(ζ ) do not depend on x (see Appendix A), b(ζ ) andb(ζ ) do change when we change from x to a translated coordinate x := x − x 0 , as we see in Lemma 1. For the continuous spectrum ζ is real and then this change affects b(ζ ) and b (ζ ), which is analogous to a phase shift of the ordinary Fourier transform under a translation of the space coordinate. For the discrete spectrum ζ 0 is imaginary, and then this change affects |b(ζ 0 )| and |b(ζ 0 )| instead.
Lemma 2 (Space reversal): If the scattering matrix of a potential q(x) is S(ζ ) as in (11), then the scattering matrix of the potential q (x) ≡ q(−x) is Proof: The scattering problem for q (x) is equivalent to the scattering problem for q(x ) with x := −x. The Jost solutions in the mirrored coordinate are Equation (64) follows from filling these out in the definitions of the scattering parameters, (55) to (58). Corollary 1 (Even symmetric potential): For a potential that is even symmetric, i.e. q(x) ≡ q(−x), by Lemma 2 Furthermore, since a(ζ 0 ) = 0 for every eigenvalue ζ 0 , from (13) and (69) norming constants of even symmetric potentials satisfy b(ζ 0 ) = ±1.

APPENDIX C OTHER BASES FOR SOLVING THE SCHRÖDINGER EIGENVALUE PROBLEM
The derivations in this paper for the S basis can be translated into other bases by means of similarity transformations. Let b indicate any such basis, then the basis dependent variables are related to the S basis as where Hereafter we discuss some bases found in literature.

A. AKNS BASIS
The AKNS system for the KdV equation [33] is found with the choice where r(x) ≡ −1. The transformation matrices that relate the AKNS basis for the KdV equation to the S basis are We remark that (24) and (25)  By setting r(x) differently, the AKNS system can be used for other non-linear differential equations [33]. In particular, the choice r(x) = ±q * (x) results in a system for the NSE, which is also known as the ZS system. In this paper we mean by AKNS the AKNS system with r(x) ≡ −1 and refer to the NSE version as ZS.
A variant of the AKNS system is and is for the KdV equation related the S basis by the transformation matrices This basis leads for the KdV equation to a numerically more accurate calculation of the continuous spectrum [29, Footnote 3].

B. COMPANION BASIS
Another choice results in a companion system: The transformation matrices that relate this basis to the S basis are The advantage of this basis is that it only requires computations on real numbers for both the discrete spectrum (ζ ∈ I) and the continuous spectrum (ζ ∈ R), whereas the other bases in this paper need complex arithmetic for the continuous spectrum. This advantage is employed by e.g. [

C. OSBORNE BASIS
The following choice leads to a close relative of the S basis: The transformation matrices are given by This basis is implicitly used in [17], [28], although with all matrix equations transposed compared to this paper.

APPENDIX D NAIVE CALCULATION OF THE NORMING CONSTANT
A naive numerical calculation of the norming constant of a certain eigenvalue ζ 0 would use (10) to calculate the scattering matrix S(ζ 0 ) which contains the norming constant b(ζ 0 ). However, this calculation is ill-conditioned. As an illustration thereof we will add here a particular small perturbation to this calculation and show that this has a large effect on the result. Let us consider a potential q(x) that is zero for all x outside a window (X − , X + ) and with at least one eigenvalue ζ 0 for which we have numerically calculated the scattering matrix S(ζ 0 ), in whichâ(ζ 0 ) = 0. Now we perturb the potential near x = X + as where we assume the potential to be approximately constant for x ∈ (X + − ε, X + ). Then the scattering matrix of the perturbed potential becomeŝ where first steps back, thereby cancelling the unperturbed potential q(x), then steps forward with the perturbed potential q µ (x).
Finally after filling out the first order approximation of (89) in (86), we find The numerical error in the calculated norming constant is typically large because of the exponential factor in (91). As an example suppose that we have an even symmetric potential, then b(ζ 0 ) = ±1 for all eigenvalues (see Corollary 1 in Appendix B). Then if |δâ(ζ 0 )| ≈ 10 −15 with δ 0.1, the estimation error according to (91) is already 10% for −jζ 0 X + ≈ 16.
It may seem as if we could apply a change of variables x → x + x 0 to makeâ(ζ 0 )b −1 (ζ 0 ) exp(−2jζ 0 X + ) arbitrarily small by lowering X + . However, application of Lemma 1 in Appendix B shows that such a change of variables leaves this quantity unchanged.

APPENDIX E LINK WITH PREVIOUS WORK ON THE NON-LINEAR SCHRÖDINGER EQUATION
Reference [30] uses for the NSE a variant of the ZS basis in which any state transition matrix on an interval that covers all non-zero parts of the potential equals the scattering matrix. To see this let us write [33, Eqs The variables are changed in [30,§IIIA] such that Then for a potential that is zero for all x outside an interval (X − , X + ) it follows that [cf. (26)] H (X + , ζ ) = H H (X − , X + , ζ ) H (X − , ζ ) , where Hence propagating 1 0 T forward up till the cutting point results in a L (ζ 0 ) b L (ζ 0 ) T , whereas propagating 0 1 T backward up till the cutting point results in b R (ζ 0 ) a R (ζ 0 ) T .
Finally b(ζ 0 ) is calculated using only the first element of both results:b which isb 1 (ζ 0 ), the estimator according to the first row of (37), with the sign difference explained in Remark 1.
This isb 2 (ζ 0 ), the estimate that only makes use of is the second row of (37).

APPENDIX F GENERATION OF A MULTISOLITON POTENTIAL FOR THE KDV EQUATION
For the generation of the multisoliton potential in Section IV-B, we made use of an algorithm from [44]. However, we observed that the calculation was only well-conditioned near the centre of the solitons. In order to use it for the 'tails' of the solitons, we adapted the algorithm as described below.

A. SIMPLIFICATION OF THE DETERMINANT EQUATION
Consider from [44] the unnumbered equation between (22) and (23). For this equation to be valid, the denominator matrix must be invertible. Hence, dividing the two determinants is equivalent to taking the determinant of the product between the inverse of the denominator matrix (hereafter D) and the numerator matrix (hereafter N). Since these two matrices are equal except for the last column, the aforementioned product (hereafter Q) has a particular structure that considerably simplifies taking its determinant: where * denotes a number that is not of interest. Hence, |N|/|D| = c and this quantity c can be found by solving Dy = n for y with a suitable linear solver, where n is the last column of N and c is the last element of y.

B. SCALING OF THE MATRIX EQUATION
Consider again from [44] the unnumbered equation between (22) and (23) and the alternative calculation thereof described above. We observed that these equations become badly scaled in the 'tails' of the potential of the KdV equation, because half of the β n parameters tend to ∞. This can be solved by reformulation of the equation in terms of the α n parameters, i.e. by left-multiplying both N and D by a diagonal matrix diag(α 1 , α 2 , . . . , α 2N ) = diag(β −1 1 , β −1 2 , . . . , β −1 2N ). This results in a condition number of the denominator matrix D that tends to a constant as |x| → ∞. We remark that in some cases this scaling results in a worse condition number near the solitons centers. For such cases we compare the condition numbers between the original scaling and the one described here for each potential sample and choose the best conditioned one for the calculation of that potential sample.