{'fraction_non_alphanumeric': 0.06291119331900849, 'fraction_numerical': 0.020556896971845026, 'mean_word_length': 3.798469387755102, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computation of Large-Genus Solutions of the Korteweg-de Vries Equation

We consider the numerical computation of finite-genus solutions of the Korteweg-de Vries equation when the genus is large. Our method applies both to the initial-value problem when spectral data can be computed and to dressing scenarios when spectral data is specified arbitrarily. In order to compute large genus solutions, we employ a weighted Chebyshev basis to solve an associated singular integral equation. We also extend previous work to compute period matrices and the Abel map when the genus is large, maintaining numerical stability. We demonstrate our method on four different classes of solutions. Specifically, we demonstrate dispersive quantization for"box"initial data and demonstrate how a large genus limit can be taken to produce a new class of potentials.


INTRODUCTION
Consider the Korteweg-de Vries (KdV) equation, written in the form subject to periodic boundary conditions. A main outcome of this work is an efficient numerical method for the computation of the inverse scattering transform associated with (1), that is, a numerical inverse scattering transform for the Schrödinger operator with a periodic, piecewise smooth, potential.
Following [21,31,34], we formulate a Riemann-Hilbert problem for the so-called Bloch eigenfunctions of the Schrödinger operator. In the finite-gap case, two eigenfunctions are classically used to construct the associated Baker-Akhiezer function on a hyperelliptic Riemann surface. A key improvement we make here over the numerical approach in [31] is that through a transformation z 2 = λ we pose Riemann-Hilbert problem with jumps supported on the gaps. This idea was used with limited scope in [33].
The key numerical innovations come from the use of the Chebyshev-U and Chebyshev-V (third and fourth kind) polynomials and their weighted Cauchy transforms. These weighted Cauchy transforms encode the singularity structure of the solution of the Riemann-Hilbert problem we pose and allow for an extremely sparse representation of the solution. This is demonstrated in Figure 13 below.
The developed numerical method can handle high-genus potentials -Riemann-Hilbert problems with jump matrices supported on thousands of intervals. This ability stems from the choice of Chebyshev-U and V basis, which is one of the new ideas in this work. But there are additional developments that are required to even pose that Riemann-Hilbert problem in the inverse scattering context. These developments are related to computing the period matrix for a basis of holomorphic differentials when the genus is high. Through a judicious choice of the basis of 1 u(x, 1.03π) x FIGURE 1. u(x, 1.03π) in (170) with (171) solving KdV equation in the form u t − uu x + u xxx = 0 (as in [7]). This plot shows dispersive quantization. The solution appears to be piecewise smooth at rational-times-π times and fractal otherwise. For the KdV equation, this was first observed by Chen and Olver, see [7], for example. These plots are produced using a genus g = 300 approximation. More details concerning the computation of this solution can be found in Section 5.4. holomorphic differentials, we develop an approach that appears stable as g → ∞. Furthermore, our use of Chebyshev-U and Chebyshev-V polynomials and their weights is predicated on having a potential that produces a Baker-Akhiezer function with poles at band ends. As this is not the generic setting, following [31], we construct a parameterix Baker-Akhiezer function to move the poles to the band ends, without loss of generality. A simplification explained in Section 2.3 allows for g to be large with this approach.
In this paper we do not discuss, in detail, the computation of the direct scattering transform for the Schrödinger operator with a periodic potential -the computation of the periodic, antiperiodic and Dirichlet spectra. We do accomplish this using existing standard techniques but consider any improvement on these approaches as important future research topics. In the case of a "box" potential (see (171) below) we can compute the Bloch eigenfunctions explicitly and apply simple root-finding routines to compute the requisite spectra. In Figure 1 we plot the evolution of this infinite-genus box potential to time t = 1.03π using a genus g = 300 approximation. Realization of dispersive quantization in a nonlinear setting is clear -the solution appears to be piecewise smooth whenever t is a rational multiple of π [3,7,23,28].
1.1. Relation to other work. We emphasize that the computation of finite-genus solutions is a nontrivial matter. Lax's foundational paper [20] includes an appendix by Hyman, where solutions of genus 2 were obtained through a variational principle. The classical approach to their computation goes through their algebro-geometric description in terms of Riemann surfaces, see [8] or [17], for instance. While very effective, this approach has only been applied to relatively small genus Riemann surfaces. Yet another approach is by the numerical solution of the so-called Dubrovin equations [2,10]. And the finite-genus solution is easily recovered from the solution of the Dubrovin equations [22,27]. We do not take this approach again because 1) the dimensionality involved may pose possible stability issues and 2) one has to time-step the solution to get to large times. The Riemann-Hilbert problem we pose has x and t as explicit parameters, and therefore the complexity associated with computing the solution at any given (x, t) value is independent of (x, t).
As mentioned above, a numerical Riemann-Hilbert approach was introduced in [31] (see also [34]). While the approach in [31] should be seen as the precursor to the current work, it was only successful for small genus solutions and was too inefficient when the genus is larger than, say, 10.

1.2.
Outline of the paper. The paper is laid out as follows. In Section 2 we review the inverse spectral theory for the Schrödinger operation with a periodic or finite-gap potential, connecting it to an underlying Riemann surface (in the finite-gap case) and the associated Baker-Akhiezer function. In Section 2.3 we discuss the parametrix Baker-Akhiezer function that allows the movement of poles and in Section 2.4 we begin formulating a Riemann-Hilbert problem satisfied by the planar representation of the Baker-Akhiezer function. In Section 3 we convert the Riemann-Hilbert problem to a singular integral equation on a collection of intervals. We look for solutions in a weighted L 2 space. In Section 4 we discuss the numerical solution of the singular integral equation from the previous section, discussing both preconditioning and adaptivity of grid points. In Section 5 we discuss the comptuation of various solutions of the KdV equation. Specifically, in Section 5.1 and 5.2 we compute solutions with prescribed spectral data. In Section 5.2 we give a formal universality result that demonstrates how primitive potentials can be obtained in a largegenus limit. In Section 5.3 we solve the initial-value problem for the KdV equation with smooth initial data. In Section 5.4, we give an extensive treatment of the numerical solution of the KdV equation with "box" initial data.
This work gives rise to many interesting questions. The work here, while empirically valid, comes with no rigorous error bounds and the full numerical analysis of the method is an open problem. Similarly, we provide no error bounds for the approximation of an infinite genus potential by one of finite genus. The reconstruction formula (53) appears to imply that the errors will be small if one removes gaps such that α j+1 − β j is small. But this removal has a non-trivial impact on γ k (x) for k = j and that error needs to be estimated. This leads to the question of understanding both the large g limit of the period matrix of our basis of holomorphic differentials and the large g limit of the singular integral equation we formulate. These issues will be addressed in future work. There is also some room for improvement in the complexity of the numerical method. A significant improvement on the complexity would be to put it inside a matrix-free framework using some incarnation of the fast multiple method [6]. Code used to generate the plots in the current paper can be found at [30].
Before we proceed, we give a remark that details our notational conventions.
Remark 1.1 (Notational conventions). We use capital boldface letters, e.g., M, to denote 1 × 2 rowvectors and to denote matrices, with the exception of the Pauli matrices, and denote the identity matrix or identity operator by I. We use lowercase boldface letters, e.g., u to denote column-vectors that are of arbitrary dimension. We use the capitalized Greek characters, e.g., Ψ, to denote functions defined on a hyperelliptic Riemann surface. Given such Ψ, we denote by ψ ± the (scalar-valued) components of its planar representation in the form of a rowvector which is denoted by the boldface version Ψ of Ψ. We use superscripts f ± (z) to denote the boundary values of f at a point z on an oriented contour taken from the left (+) and the right (−) side of the contour with respect to the orientation. We use fraktur a and b to denote the cycles on a Riemann surface. Lastly, for a function f : C → C we use f (u) to denote f applied entrywise to the vector u.

INVERSE SCATTERING TRANSFORM FOR PERIODIC SOLUTIONS
In this section we give a summary of the well-known inverse scattering transform associated with (1) and define the quantities relevant to the method we develop in this work, along with particular choices we make. The KdV equation in the form (1) is the λ-independent compatibility condition for the linear problems, i.e., the Lax pair, where is L is the Schrödinger operator with the time-dependent potential −q( , t), and P is the skew-symmetric operator The compatibility condition for the system of linear problems (3)-(4) yields the operator equation, referred to as the Lax equation [19], in the form which is equivalent to the KdV equation (1) in the sense that the left-hand side defines an operator of multiplication by the function −(q t + 6qq x + q xxx ), where [L, P ] := LP − P L is the operator commutator. As q evolves in time according to the KdV equation (1), (7) defines an isospectral deformation of the Schrödinger operator L.
2.1. The spectrum and the Riemann surface. For fixed t ≥ 0, consider the problem (3) for the Schrödinger operator with the time-independent potential −q( , t) = −q( ): for real periodic q with minimal period L > 0: q(x + L) = q(x). The Bloch spectrum σ B (q) associated with the periodic potential −q for the Schrödinger operator (5) is For real-valued smooth (and periodic) q, the Bloch spectrum is a countable union of real intervals We refer to the intervals [α k , β k ] ⊂ σ B (q) as bands and (β k , α k+1 ) as gaps. If the number of intervals g + 1 is finite, β N := +∞ and the last interval is [α g+1 , +∞), in which case the associated −q is called a finite-gap potential. The endpoints α j and β j , j = 1, 2, . . . , N, remain invariant as q( , t) evolves according to the KdV equation (1), and hence σ B (q 0 ) = σ B (q) for q( , t) solving (1) with q( , 0) = q 0 . The following well-known symmetry transformations associated with the KdV equation play a role in various choices we will make in this work.
Remark 2.1 (Two symmetry groups of KdV). Suppose that q(x, t) is a solution of (1).
• Galilean transformation: The functioñ is also a solution of (1) for any constant c.
• Scaling transformation: The functioñ is also a solution of (1) for any constant c.
Given q = q( , t) and α 1 = min(σ B (q)), using the Galilean symmetry transformation (12) with c = α 1 lets one map q( , t) toq( − 6ct, t) + c for which min(σ B (q)) = 0. Doing so becomes useful in the formulation of a Riemann-Hilbert problem (and of the associated singular integral equation). This transformation is employed in the numerical implementation of our method: once α 1 ∈ σ B (q 0 ) is computed for given q 0 at t = 0, we perform the spectral shift described above and then invert it to obtain q( , t) fromq( − 6ct, t) + c at a later time t > 0. Accordingly, we take α 1 = 0 without loss of generality in the remainder of this paper.
For our (computational) purposes, we restrict the theory to the finite-gap case. For q 0 giving rise to g + 1 bands (14) σ and g gaps, g ≥ 2, consider the monic polynomial P(λ) of degree 2g + 1 given by and define Σ to be the hyperelliptic (elliptic, if g = 1) nonsingular Riemann surface associated with the zero locus of F(λ, w) := w 2 − P(λ). The points (α j , 0), (β j , 0), j = 1, 2, . . . , g, and (α g+1 , 0) on Σ are branch points for the projection (z, w) → z and there is a single point at ∞ on Σ. For P 0 = (λ 0 , w 0 ) ∈ Σ we have the following choices of a local coordinate ζ: • If P 0 is not a branch point and not ∞ ∈ Σ, then for (λ, w) near P 0 we may take essentially λ to be a local coordinate, so we write for |ζ| sufficiently small • If P 0 = E k = (λ k , 0) for some k, then for (λ, w) near P 0 we may write • Finally, if P 0 = ∞ ∈ Σ, then for (λ, w) near P 0 we may write In all three cases λ(ζ) and w(ζ) are locally holomorphic functions of ζ in a neighborhood of ζ = 0 with non-zero derivatives at ζ = 0, making them locally injective, and ζ(P 0 ) = 0. Define the branch of square root R(λ) of P(λ) to be the (unique) single-valued function that is analytic in C \ σ B (q 0 ) and that satisfies R(λ) 2 = P(λ) along with the asymptotic behavior Importantly, the boundary values R(λ) taken on the bands σ B (q) are real-valued and R(λ) is purely imaginary on the gaps, namely for λ ∈ (β j , α j+1 ), j = 1, 2, . . . , g. Assuming that the bands constituting σ B (q 0 ) are oriented in the increasing direction of the real line, we define the boundary values R ± (λ) of R taken on σ B (q 0 ) as: We denote by +R(λ) the function that's defined for all λ ∈ C which coincides with R(λ) for λ ∈ C \ σ B (q 0 ) and also satisfies +R(λ) = R + (λ) for λ ∈ σ B (q 0 ); and define the sheets Σ ± of Σ by The Riemann surface Σ is topologically equivalent to a sphere with g handles, which is obtained by gluing the two sheets Σ ± along the edges of the cuts (the bands) [α j , β j ], j = 1, 2, . . . , g, and [α g+1 , ∞]. We define the cycles {a j , b j } g j=1 , which constitute a homology basis for Σ, as depicted in Figure 2; and we denote by {ν 1 , ν 2 , . . . , ν g } the basis of normalized holomorphic differentials on Σ that satisfy The g × g matrix B defined by is the Riemann matrix for Σ. B is symmetric with Re(B) < 0. Although finite-gap solutions of (1) have, in principle, representations given in terms Riemann theta functions Θ which are based on the Riemann matrix B, our method does not require at all any explicit knowledge of the Riemann matrix B.
An illustration of the hyperelliptic Riemann surface Σ and the choices for the a-and the b-cycles.
We will use a particular basis of differentials that is ideal for our computational purposes, so some observations are in order. Let {ω 1 , ω 2 , . . . , ω g } be an arbitrary basis of holomorphic differentials on Σ and define the g × g matrices A andB by (25) A jk := It is well-known and easy to see that A is nonsingular since otherwise one can find a nontrivial linear combination of {ω j } g j=1 that has vanishing a-cycles, implying that the resulting holomorphic differential is identically zero, and hence contradicting the independence of the differentials ω j . Note that there exists scalars C lk , 1 ≤ l, k ≤ g, such that ν l = ∑ g j=1 C lj ω j . Let C be the matrix of these scalars. Then we have (26) 2πiδ lk = a k implying that CA = 2πiI, and hence C = 2πiA −1 , which yields the relation A classical choice for the basis {ω 1 , . . . , ω g } is But from a computational point-of-view, when g is large, this basis is ill-conditioned. It is better to chose 2.2. The Baker-Akhiezer function. Our approach to compute solutions of (1) is based on numerical solution of a RH problem satisfied by a suitable renormalization of a Baker-Akhiezer function. We now give a series of definitions and then give construction of the relevant Baker-Akhiezer function built from certain solutions of the spectral problem Definition 2.2. For the hyperelliptic Riemann surface Σ defined by w 2 = P(λ), a divisor D on Σ is a formal sum D = n 1 P 1 + n 2 P 2 + · · · + n m P m , (30) where n j ∈ Z and P j ∈ Σ for j = 1, 2, . . . , m. A divisor is called positive if n j > 0 for all j, and the degree of a divisor is the number ∑ m j=1 n j .

Definition 2.3 (Abel map)
. Fix an arbitrary point P 0 on the Riemann surface Σ defined by w 2 = P(λ) and let D = n 1 P 1 + n 2 P 2 + · · · + n m P m be a divisor on Σ. The Abel map A(D) is given by where the path of integration from P 0 to P j is chosen to be the same for each = 1, 2, . . . , g.
(2) There exists a neighborhood of every point Q j , j = 1, . . . , n, such that the product Ψ(P) exp −q j ζ j (P) −1 is analytic in this neighborhood.
Such a function Ψ(P) is called a Baker-Akhiezer function.
The following theorem from [2, Theorem 2.24] is quite useful in this work. We do not define all the quantities that arise in its statement but only highlight the components that are crucial for us to proceed.
Theorem 2.5. The space V (D; Q 1 , . . . , Q n ; q 1 , . . . , q n ) is one-dimensional for a non-special divisor 1 D and polynomials q j with sufficiently small coefficients. Its basis is described explicitly by where Ω(P) is a normalized Abelian integral of the second kind 2 with poles at the points Q 1 , . . . , Q n , the principle parts of which coincide with the polynomials q j (z j ), j = 1, . . . , n, Θ is Riemann's theta function, and v is a vector of the b-periods of the integrals of Ω(P): is the Abel map and k is a vector of Riemann constants, and the integration path for the integrals Ω(P) = P P 0 dΩ and A(P) (34) is chosen to be the same.
We now focus our attention to solutions of (29). First, fix x 0 ∈ R and define a set of fundamental solutions {c(x; λ), s(x; λ)} of (29) that are determined by It is easy to verify that these solutions solve the following Volterra integral equations This demonstrates that these two solutions are entire functions of λ for given x because cosine and sine are even and odd functions, respectively, and the paths of integration are finite. We omit the parametric dependence of these solutions on x 0 in our notation. Next, we define the monodromy 1 The divisors D that we encounter in this work will always be of the form D = P 1 + P 2 + . . . + P g for distinct P j . Such divisors are non-special [13]. 2 A normalized Abelian differential is a meromorphic differential that is normalized to integrate to zero over all the a-cycles. operator T for (29) by (T ψ)(x) := ψ(x + L) and represent the action of T on the set of fundamental solutions constructed above. Consider the first-order system equivalent to (29) along with its fundamental solution matrix which is unimodular (see [13, §1.6]) and satisfies F(x 0 ; λ) = I. As x → c(x + L; λ) and x → s(x + L; λ) also define solutions of (29) thanks to the periodicity of q 0 , for some x-independent 2 × 2 matrix T(λ). Evaluating both sides at x = x 0 yields which called the monodromy matrix. It is an entire matrix-valued function of λ and det(T(λ, x 0 , L)) = 1. Then, for a solution ψ( ; λ) of (29), we have for some a(λ) ∈ C 2×1 , and hence for any positive integer n. This and the unimodularity of F(x; λ) imply that for given λ, µ(λ) is an eigenvalue of T with an eigenfunction ψ( ; λ) in the (two dimensional) solution space of (29) if and only if µ(λ) is an eigenvalue of T(λ) with an eigenvector a(λ). It can be shown that λ ∈ σ B (q 0 ) if and only if T has an eigenvalue µ(λ) (in the solution space of (29)) with |µ(λ)| = 1, and that |µ(λ)| = 1 implies λ ∈ R (see [13, Lemma 1.6.4, Lemma 1.6.7]). To find µ(λ), set ∆(λ) := 1 2 tr(T(λ)) and note that ∆(λ) is entire and does not depend on x 0 . Then the eigenvalues of T(λ) where the square root is taken as the principal branch (positive for −1 < ∆(λ) < 1), and the Bloch spectrum σ B (q 0 ) consists of λ ∈ R such that |∆(λ)| ≤ 1. The Bloch eigenfunctions ψ ± ( ; λ) are bounded solutions of (29) that are eigenfunctions of T for λ ∈ σ B (q 0 ) with the normalization ψ ± (x 0 ; λ) = 1, hence they are obtained by choosing the first row of a(λ) in (43) to be equal to 1 in solving Setting , it follows that χ ± (x; λ) are periodic functions of x with period L and are independent of the choice of x 0 (see [12,Lemma 1.1]). Using the independence of the Wronskian Wron(ψ ± (x; λ), c(x; λ)) from x yields the formula For the finite-gap case (i.e., g < ∞) we are considering, ∆(λ) 2 − 1 = 0 has finite (and odd) number of simple roots, which are the band endpoints [12]. From the representation [12, Eqn.
Remark 2.6. The zeros of Ψ(x; P) are at the points where λ = γ j (x), and they lie also in the gaps. It's well-known that the potential q 0 (x) can be recovered via the formula see, for example, [14]. Our method for obtaining q 0 from Ψ makes no reference to this formula, and hence avoids root-finding. (3) can be constructed at a given fixed time t as q( , t) evolves according to the KdV equation. Let ψ [t] ± (x; λ) denote these solutions and we have ψ (48) and the normalization ψ ± (x 0 ; λ) = 1, they do not provide a set of simultaneous solutions of (3)-(4) as they do not satisfy (4). A calculation identical to [21,Proposition 6.2] shows that for x-independent coefficients d ± (t; λ) that are given by These are obtained by evaluating (55) at x = x 0 . Again from [21, Proposition 6.2] (see also [21,Proposition 3.3]) we have the asymptotic behavior Following [21], one uses the solution φ ± (t; λ) of ± (x; λ)φ ± (t; λ) define a set of simultaneous solutions of (3)- (4). As proved in [21, Proposition 6.3], φ ± (t; λ) satisfy Moreover, the product in (59) fixes the poles of ψ ± in time, see [21,Proposition 6.3]. Thus, with ψ ± (x, t; λ) we introduce a Baker-Akhiezer function Ψ(x, t; P) on the Riemann surface with all the same properties as (51) with the exception of the replacement of the asymptotics with 2.3. Moving poles to the band endpoints. The procedure described here resembles what was employed in the earlier works [31,32] (see also [34,Chapter 11]). However, we make an observation that enables the treatment of the case when the genus g is large.
is then the Abel map restricted to the second sheet. The following properties of the theta function are now needed where e j is the jth column of the g × g identity matrix and B is the Riemann matrix. Then note that for a vector n of ones and zeros. Then we compute Note that from (51), for a non-special divisor D = ∑ where, as before, k is the vector of Riemann constants with base point α 1 . So, for two non-special divisors D = ∑ g j=1 P j and D = ∑ For P = (λ, w) ∈ Σ, define π(λ, w) = λ. Then it follows that if P j is on the first (second) sheet of Σ then ∆(λ; d, v) has a zero at π(P j ) in its first (second) column. Similarly, if P j is on the first (second) sheet of Σ then ∆(λ; d, v) has a pole at π(P j ) in its first (second) column. Now suppose that λ is not a pole of either column of ∆. Then Choose the divisor now has poles at the right endpoints of the gaps, namely the points where λ = α 2 , α 3 , . . . , α g+1 . We arrive at the following proposition.
An illustration of the choices for the b-and theb-cycles on the hyperelliptic Riemann surface Σ.

The Riemann-Hilbert problem.
Let Ψ( ; x, t) : C 1×2 \ σ B (q) → C denote the row vector planar representation of the Baker-Akhiezer function associated with q(x, t): which satisfies the "twist" jump condition and has the asymptotic behavior Here the power function λ → λ 1 2 is defined to be analytic on C \ [0, +∞), satisfying λ 2 (x + 4λt) and observe that θ(λ; x, t) has a jump discontinuity across the half-line (0, +∞), which we orient from λ = 0 to λ = +∞. Based on the considerations in the previous section, we assume that the poles of Ψ occur at the points in the divisor (65).
Remark 2.8. An important calculation to make here is to define Then apply [33, Theorem 2.1] to see that K (and therefore M, and hence Ψ(x, t; P)) is a simultaneous solution of an appropriate version of the Lax pair for the KdV equation.
To control the oscillatory factors in the jump matrices above, we seek a function G(λ; x, t) that is analytic for λ ∈ C \ [0, +∞) satisfying for some constants Ω j , and normalized to satisfy G(λ) → 0 as λ → ∞. It is easy to see that is analytic for λ ∈ C \ [0, +∞), admits continuous boundary values on [0, +∞) which satisfy the jump conditions (78)-(79). Observe that Thus, in order to have G(λ) = o(1) as λ → ∞, we need to have m k ≡ 0 for k = 1, 2, . . . , g, which yields the conditions This is a linear system of g equations for the constants Ω j = Ω j (x, t), j = 1, 2, . . . , g. Taking a linear combination of these equations we can instead consider for any basis {p 0 , . . . , p g−1 } for polynomials of degree at most g − 1. Taking into account the orientation of the a-cycles depicted in Figure 2 and the sign change that occurs from passing from one sheet to the other, we have from (25), choosing p j−1 so that p j (ζ)/R(ζ) dζ = ω j , and hence (83) reads Thus, the coefficient matrix for the linear system (83) is nothing but a constant multiple of the matrix A of a-cycles of the basis of differentials {ω k } g k=1 , which is nonsingular. Therefore, the system (83) is uniquely solvable. Moreover, because θ + (ζ; x, t) is real-valued and R(ζ) is purely imaginary for ζ ∈ (β j , α j+1 ), j = 1, 2, . . . , g, it follows that Ω j , j = 1, 2, . . . , g, are all real valued. Using the basis of differentials in (28) results in a linear system which can be solved in a numerically stable fashion as g becomes large.
Remark 2.9. The computation of the integrals that appear in this section and the computation of the Abel map is discussed in [31]. There is a numerical subtlety here. If one computes the Abel map A(λ) for λ near a branch point and λ is known to within an error , that error may be amplified to be on the order of √ . So, if q 0 is such that γ j (x) is near a branch point, the computation of the of v j in Proposition 2.7 may suffer increased errors. In practice, one can choose x 0 in the initial scattering theory to move this away from the branch point. Various schemes can be employed to find a good choice of x 0 . Choosing x 0 randomly is often sufficient.

A SINGULAR INTEGRAL EQUATION FOR THE SOLUTION OF THE A RIEMANN-HILBERT PROBLEM ON CUTS
In this section we describe a numerical method to approximate (87). But before we do that, we need to make one more transformation to remove the non-trivial jump that has infinite extent. Recall that from the discussion following Remark 2.1 that we take α 1 = 0 without loss of generality. Define For convenience we write α j =: a 2 j , a j ≥ 0, j = 1, 2, . . . , g + 1, (91) Note that for Re(z) > 0 if Im(z) > 0 then Im(z 2 ) > 0. And similarly, if Re(z) > 0, Im(z) < 0 implies Im(z 2 ) < 0. And if Re(z) < 0, then these implications are flipped. From this we find that S(z; x, t) only has jumps on the (symmetric) collection of intervals where it satisfies: where we have reoriented the intervals (−a j+1 , −b j ) so that all of the intervals in (93) are oriented from their left endpoint to the right endpoint. Moreover, S(z; x, t) is normalized so that S(z; x, t) = 1 1 + O(z −1 ) as |z| → ∞.
Next, we will want a formula to recover q(x, t), the solution of the KdV equation (1), directly from a representation of S(λ; x, t) as function of λ. As |z| → ∞, write Supposing this limit is taken in the upper-half of the z-plane, this then implies that Then we recall that where the entries of the 1 × 2 vector Ψ(λ, x, t) are solutions of L(t)ψ = λψ, see (71). So, consider the function m(λ; x, t) = ψ + (λ; x, t)e −iθ(λ;x,t) : Adding these so as to eliminate the ∂ x ψ + term, we find It follows from (32) that both ∂ xx m(λ; x, t) and ∂ x m(λ; x, t) decay at infinity, giving the recovery formula In other words, we have as λ → ∞ where m g+1 (x, t) denotes the coefficient of the term proportional to λ −g−1 in the expansion (81). Thus, we arrive at 3.1. Weighted spaces. We now formulate a singular integral equation on a direct sum of weighted L 2 spaces. Define Then set, for y ∈ I j (106) w j (y) := where each weight w j (y) is understood to vanish outside its domain of definition. Then define It is convenient to order the component functions (each of which is 2 × 1 row-vector-valued) for i.e., right multiplication by the jump matrix and division by the weight on I j , respectively.
We define the boundary values of (111) whenever the following limits exist: When the domain of the weight w is clear from context we write C w , C ± w . When w ≡ 1 we write C Γ , C ± Γ . These operators are understood to apply to vectors component-wise.

Definition 3.1. A function S(z; x, t) is a solution of the Riemann-Hilbert problem
for U ∈ L 2 w j I j and the jump condition (113) is satisfied for a.e. z ∈ I j for each j. Further, for j = k we use the notation Theorem 3.2. Suppose S(z; x, t) satisfies the following 4 (2) The jump condition (113) is satisfied for a.e. z ∈ I j for each j. Proof. The first condition imposes that S − C is an element of the Hardy space of the upper-half and lower-half planes [15]. This implies that the boundary values from above and below exist a.e. Furthermore, it also implies that S(z; x, t) is given by the Cauchy integral of its boundary values: Imposing the jump condition S + (z; x, t) = S − (z; x, t)J j , z ∈ I j for each j results in the following system of singular integral equations that are satisfied by U k , k ∈ {±1, ±2, . . . , ±g}, It is important to note that C w j U j (z) = C − w j U j (z) = C + w j U j (z) if z ∈ I j . Then we consider the following block operator on L 2 w j I j 4 Here · is any norm on C m×2 .
We now state some observations that motivate the preconditioning we employ in solving (118) numerically, which is described in Section 4.3. The linear system obtained from discretization of (118) upon preconditioning ends up being extremely well-conditioned; see Figure 12. First, one can prove the following lemma. The following is then immediate and is the heuristic that motivates the use of the aforementioned preconditioner in the numerical procedure. We will present the proof of Lemma 3.3 along with an analytical justification of the preconditioning and the convergence of the numerical method proposed in this work to solve (118) in a forthcoming paper.

NUMERICAL INVERSE SCATTERING
In this section we develop a numerical method to solve the Riemann-Hilbert problem in Definition 3.1. We consider the Chebyshev-V and Chebyshev-W polynomials which are also known as the Chebyshev polynomials of the third and fourth kind, respectively. The polynomials (V n (y)) n≥0 satisfy 0 < lim y→∞ y −n V n (y) < ∞ as well as for the Kronecker delta, δ n,m . Similarly, the polynomials (W n (y)) n≥0 satisfy 0 < lim y→∞ y −n W n (y) < ∞ as well as The same identity also holds for the case b > a > 0. In other words, Cauchy integrals over general intervals with these weights can be computed by first mapping a function to the interval [−1, 1], computing the Cauchy integral for the mapped function and then mapping back. We do note that b a P n (y; 4.1. Computing Cauchy integrals. As is well-known, real orthonormal polynomials (p n ) n≥0 (with positive leading coefficients), on the real axis, with respect to a probability measure µ satisfy a three-term recurrence relation yp n (y) = A n p n (y) + B n p n+1 (y) + B n−1 p n−1 , (126) for recurrence coefficients (A n ) n≥0 , (B n ) n≥0 . What is maybe less well-known is the weighted Cauchy transforms c n (z) = C µ p n (z) := 1 2πi R p n (y) y − z µ(dy), n ≥ 0, (128) satisfy the same recurrence with different initial conditions, and in particular Remark 4.1. For Chebyshev-V and Chebyshev-W polynomials we have, respectively, There are some subtleties in solving the recurrence for the Cauchy transforms. For z in the complex plane, away from the support of µ, (p n (z)) n≥0 represents an exponentially growing solution of the three-term recurrence while (c n (z)) n≥0 is an exponentially decreasing solution. Thus, evaluating c n (z) by forward recurrence is inherently unstable. Consider the case where µ(dy) has its support on [−1, 1]. In practice, the following is effective [25]: (1) For z inside a Bernstein ellipse with minor axis O(1/n), solve for c n (z) by forward recurrence allowing one to easily compute the boundary values of c n on [−1, 1] from above and below. (2) For z ouside a Bernstein ellipse with minor axis O(1/n), solve the boundary value problem with the adaptive QR algorithm [26].
We also define projection operator of evaluation of a function at points where to be truly precise, S should be an ordered set. Now suppose f : [a, b] → C can be written as where the choice of m for our purposes is discussed in Section 4.4. The discretized versions of C ± w , w = w a,b are given by (136) Note that each row of these matrices can be constructed either by forward recurrence or via the back substitution step of the adaptive QR algorithm, depending on where the evaluation points are located in the complex plane. We now demonstrate the discretization of (118) in the case g = 1.
Example 4.3. We seek vector-valued functions U ±1 : I ±1 → C 1×2 . So, write We write out the full system of equations for scalar-valued functions explicitly: For z ∈ I 1 and for z ∈ I −1 In block-operator form: The discretized version is then where each entry in the right-hand side vector in (144) is a constant vector for given (x, t).
Lastly, we need to consider the computation of ∂ x s 1 (x, t) from (96) in order to use (104). First, we note that if u j,1 in (118) is given by then by the orthogonality of the polynomials implying that we need to solve for the x-derivative of the coefficients in the expansion (147). To do this, the linear system Aγ = ω can be differentiated to find 4.3. Preconditioning. The discretization described in Example 4.3 is easily extended to find a discretization of the operators S and diag(S ) where we expect the discretization of diag(S ) to become a good preconditioner for S in light of Lemma 3.4. In practice, we find it works well to use a discretization of diag(S (I 1 ; Ω 1 ), S(I 2 ; Ω 2 ), . . . , S(I g ; Ω g )) (150) as a block-diagonal preconditioner. We find that with this preconditioner, for a fixed tolerance, the GMRES algorithm requires a bounded number of iterations, independent of x and t. We explore this more in Figure 12 It can be shown that the solution S(z; x, t) can, in an appropriate sense, be analytically continued off the interval [−1, 1] [33]. For example, one expects the solution u j on I j of (118) to have an analytic continuation to any ellipse with foci at the endpoints of I j provided that ellipse does not intersect any other I for = j. And, it is well known that rate of exponential convergence of a Chebyshev interpolant can be estimated based on this ellipse [29]: Then for one has and consequently and similarly for [a , b ]. So, set So, we expect u j • T a,b to have an analytic extension to B ρ for any ρ such that B ρ ∩ R ⊂ [−1 − δ, 1 + δ]. To be conservative in our estimates, we use δ/2 instead and find that ρ should be chosen to be: So, given an estimate for M, and a tolerance > 0, we can choose m so that and this provides an a priori guide as to how to choose m in the discretizations C ±  compute the associated potential and its evolution under the KdV flow. To demonstrate this we make first make the following choice for the gaps in the λ-plane.
While, as we demonstrate in the next section, we can compute the corresponding solution q(x, t) for g large, the solution oscillates wildly and is difficult to visualize. For this reason we plot the solution for smaller values of g over short time ranges. See Figure 4.
Proof. Expanding log(z) at z = 1 gives By the mean value theorem, Then because is uniformly continuous, for any > 0 there exists a g 0 > 0 such that | (ξ j ) − (β j )| < for all j if g > g 0 , so that and the claim follows. Now, if in the above lemma, w j = v(β j ) for some continuous function v : [β 1 , α] → (0, 1) we find that Thus, the distribution of individual locations α j , β j do not influence the limiting behavior of B(λ) as g becomes large. But rather, the distribution of the lengths of the bands is the most important quantity. To see how B(λ) will arise in a Riemann-Hilbert problem consider the above choice for α j and β j , for given functions σ and v. Previously, we have moved poles in the gap [β j , α j+1 ] on one sheet of the Riemann surface to the point λ = α j . This was for numerical convenience. Here, for analytical convenience, we put the poles at λ = β j . We diagonalize the twist jump matrix for Ψ(λ; x, t): can be used 6 to remove the jumps of Ψ. Defině where Ω is the disk 7 centered at λ = α 1 = 0 with radius α + c, c > 0, and we orient the circle ∂Ω counter-clockwise. Then one finds thatΨ satisfies the following jump conditioň For λ ∈ ∂Ω, supposing that β 1 → β, 0 ≤ β < α, one has Bringing the jump from ∂Ω back to the real axis, we find the following Riemann-Hilbert problem for a limiting Ψ ∞ and has the asymptotic behavior This construction can be immediately generalized to where v j : [β j , α j+1 ] → (0, 1) is continuous and Ψ ∞ and has the asymptotic behavior (164). While full exploration of such Riemann-Hilbert problems is beyond the scope of the current paper, potentials for v(λ) ≡ 1/2 are given in Figure 5 in the case where β 1 → 0 as g → ∞ and the evolution is plotted in Figure 6.
q(x, 0) x FIGURE 5. Numerical solutions of the KdV equation as g increases with v in Section 5.2 being constant, equal to 1/2 and α g+1 = 2 + 1/g, β = 3/g. Here we also, for simplicity, take to be uniform on [β 1 , α]. The number of collocation points on each interval is obtained adaptively using the methodology in Section 4.4. It is clear from these panels that the solution converges as g → ∞. The limit may be related to the so-called primitive potentials [16] but the connection is not immediate.. Error g FIGURE 7. The error in the computation of u(x, 1) where u is the solution of (165) as measured by comparing the numerical solution with that obtained by direct integration of the KdV equation using an exponential time integration method, see [18] for a general reference and [4] for the precise method used. The error is on the order of 10 −14 for sufficiently large genus and we are unable to distinguish which numerical solution is giving the dominant contribution to the error. The number of collocation points per contour is chosen using the methodology in Section 4.4.

5.3.
Initial-value problem with smooth data. We consider the classical problem of Zabusky-Kruskal [35] u t + uu x + δ 2 u xxx = 0, (165) for x ∈ (−2, 2]. Based on Remark 2.1, since we are set to solve q x + 6qq x + q xxx = 0, we choose and then u(x, t) = aq(bx, ct). We choose δ = 0.08 and use an error tolerance of 10 −13 (see Section 4.4) to choose the number of collocation points on each interval I j . We then plot the error in computing u(x, 1) as g increases. To estimate the true error we use the exponential integrator method discussed in [4] motivated from the work in [18] to compute the "true" solution. Exponential convergence is seen in Figure 7. The evolution of the corresponding solution is given in Figure 8.
Remark 5.2. To be able to compute this solution, one needs to be able to compute the spectrum. We use the Fourier-Floquet-Hill method [9] to compute the periodic/anti-periodic eigenvalues and use a Chebyshev method to compute the Dirichlet eigenvalues. This latter method can be found implemented in both Chebfun [1] and ApproxFun [24].
x FIGURE 8. The numerically compute evolution of (165) using a genus g = 12 approximation. The number of collocation point on each contour is chosen adaptively using the methodology in Section 4.4. 5.4. Initial-value problem with "box" data. To be able to solve the initial value problem for the box initial condition, q 0 (x) = q 0 (x + L), we need to compute the forward spectral theory for (5). In this case a basis of solutions c(x; λ) and s(x; λ) to (3) normalized as in (35) and (36) can be computed explicitly, and these solutions contain all information needed to compute the forward spectral theory. In particular, we have The Dirichlet eigenvalues γ 1 (x 0 = 0) < γ 2 (x 0 = 0) < γ 3 (x 0 = 0) < . . . of (5) are then the zeros of (168) c(L; λ) The band ends α 0 < β 0 < · · · α g < β g < · · · are the zeros of (169) ∆(λ) 2 − 1.
Error: u(x, 0; g) ∞ g FIGURE 9. A test of the convergence of the numerical approximation of u(x, 0; g) of u(x, 0) as g increases. Because of expected non-uniform convergence due to a Gibbs-like phenomenon, we take a uniform grid x on [0.1, π − .1] and use u(x, 0; g) ∞ as a proxy for the error (u(x, 0) = 0 on this interval). Due to increased oscillations as g increases, this can only be thought of as an estimate for the true error. The best fit line is given by (g) ≈ 0.89 × g −0.93 which appears to be consistent with O(g −1 ) convergence.
These are easily computed using standard root-finding techniques. In practice, we find it convenient to use high-precision arithmetic here so that one is sure where future errors are incurred. Since this is an infinite genus potential, we specify a finite g to truncate the spectrum, setting β g+1 = ∞, resulting in an approximate solution u(x, t; g). The convergence of this approximation is slow, but reliable, and this is investigated in Figure 9.
We investigate various aspects of computing q(x, t) with initial data q 0 as above. We focus on the case discussed by Chen and Olver [7]. Specifically, if one chooses w = π/ √ 6, L = 2π/ √ 6, h = −1/2, then u(x, t) := −q(6 −1/2 x, 6 −3/2 t), is the solution of u t + u xxx = uu x with initial data extended periodically. This allows us to reproduce much of the phenomenon in [7]. In Figure 10 we demonstrate dispersive quantization. When t is a rational multiple of π, u(x, t) as a function of x appears to be piecewise smooth and slowly varying. When t is an irrational multiple of π the solution appears to have a fractal nature. One interesting observation we make here is that while Chen and Olver remark in [7] that it is not clear from their numerical method if there are truly oscillations between jumps at rational-times-π times, our numerics in Figure 11 indicate that these oscillations will disappear as the genus increases.
We also use this problem to illustrate some important aspects of the numerical method we have developed. First, recall the matrix A in (145). We plot the eigenvalues of A(0, 0) in the left panel of Figure 12. Here we use 10 collocation points for I j if |j| ≤ 4 and 3 collocation points otherwise. Then in the right panel of Figure 12 we show the preconditioned matrixÃ(0, 0) −1 A(0, 0) whereÃ is obtained from a discretization of (150). It becomes clear that the eigenvalues become localized near λ = 2. This problem is extremely well conditioned and GMRES will converge in just a few iterations.
Lastly, in Figure 13 we display how the magnitude of the computed Chebyshev coefficients γ i,j in (147) depends on I j . In the top-left panel of Figure 13 we see that γ 0,j = O(1). This is not unexpected because the recovery formula (148) weights these coefficients by the gap lengths. What is rather surprising is how, for large |j|, the second coefficient in (147) decays rapidly, see the top-right panel of Figure 13. The third coefficient decays even more rapidly, see the bottom panel of Figure 13.
u(x, 0) x u(x, 0.1π) x u(x, 0.1) x u(x, 0.5π) x u(x, 0.01π) x FIGURE 10. The evolution of u in (170) with (171). These plots show dispersive quantization where the solution appears to be piecewise smooth at rational-timesπ times and fractal otherwise. This was first observed by Chen and Olver, see [7], for example. These plots are produced using a genus g = 300 approximation and using 10 collocation points on I j if |j| ≤ 4 and two collocation points otherwise. This choice is justified by Figure 13.
u(x, 0.1π) x g = 500 u(x, 0.1π) x g = 900 u(x, 0.1π) x FIGURE 11. A zoomed view of u(x, 0.1π) as g increases. These plots indicate that the amplitude of the oscillations decrease as the genus increases. This leads to the conjecture that the limiting solution profile is piecewise smooth and slowly varying.

Re λ
Im λ Re λ FIGURE 12. Left panel: The eigenvalues of A(0, 0) from (145) for the potential (171). We use 10 collocation points for I j if |j| ≤ 4 and 3 collocation points otherwise. Right panel: The preconditioned matrixÃ(0, 0) −1 A(0, 0) whereÃ is obtained from a discretization of (150). The eigenvalues become localized near λ = 2. This problem is extremely well conditioned and GMRES will converge in just a few iterations.