## Abstract

We solve the focusing and defocusing nonlinear Schrödinger (NLS) equations numerically by implementing the inverse scattering transform. The computation of the scattering data and of the NLS solution are both spectrally convergent. Initial conditions in a suitable space are treated. Using the approach of Biondini & Bui, we numerically solve homogeneous Robin boundary-value problems on the half line. Finally, using recent theoretical developments in the numerical approximation of Riemann–Hilbert problems, we prove that, under mild assumptions, our method of approximating solutions to the NLS equations is uniformly accurate in their domain of definition.

## 1. Introduction

We consider the initial-value problem on the whole line for the nonlinear Schrödinger (NLS) equation
1.1
Here denotes the Schwartz class of rapidly decaying functions. When *λ*=1, we obtain the focusing NLS equation and for *λ*=−1, the defocusing NLS equation. The NLS equations describe physical phenomena in optics [1], Bose–Einstein condensates [2,3], as well as water waves [4]. Together with the Korteweg–de Vries (KdV) equation, these equations are the canonical examples of the 1+1-dimensional integrable partial differential equations.

For each of the focusing and defocusing NLS equations, there exists an inverse scattering transform (IST) [5,6]. The IST can be thought of as a nonlinear generalization of the Fourier transform on . The IST is, in general, equation specific in a non-trivial way. Whereas the Fourier transform solution of a linear equation will only change slightly as the dispersion relation changes, the IST changes in a much more complicated way.

The presence of an oscillatory dispersive tail is seen for both the focusing and defocusing NLS equations. Examination of the linear dispersion relationship for the PDEs indicates that small amplitude waves will travel at a speed proportional to their wavenumber. Unlike the KdV equation, solitons do not separate from dispersion asymptotically. These factors make traditional numerics inefficient to capture the solution for large time. The computational cost to compute the solution using time-stepping methods grows rapidly in time. The methods in the study of Klein [7] and Grava & Klein [8] are well suited for solving dispersive problems for small time. When a periodic approximation to the whole-line problem is used, one has to resolve the dispersive oscillations and increase the domain size as time increases. A thorough discussion of this can be found in [9] for the case of the KdV equation. In particular, it is shown that the computational costs grow at least like .

In the work of Trogdon *et al*. [9], the authors developed a numerical implementation of the IST for the KdV equation. The main benefit of this method is that the computational cost to compute the solution at a given *x* and *t* value is seen to be independent of *x* and *t*. In the work of Trogdon *et al*. [9], this claim was verified with numerical tests. In the study of Olver & Trogdon [10], the authors derived sufficient conditions for bounded computational cost.

In this paper, we evaluate the IST associated with the focusing and defocusing NLS equations numerically. In addition, we use symmetries to solve specific boundary-value problems on . To our knowledge, this is the first time the solution of a boundary-value problem has been computed through the IST. We compute unbounded solutions to the defocusing NLS equation which have poles. Lastly, we prove that our approximation, under a specific assumption on the underlying numerical scheme, approximates solutions of (1.1) uniformly in the entire (*x*,*t*) plane. It should also be noted that the scattering problem (see §2) for the focusing NLS equation is a non-self adjoint spectral problem and this complicates the asymptotic analysis of the problem [11]. The numerical method outlined here is not adversely affected by this additional complication.

All methods we employ are demonstrated to converge spectrally. Other approaches to computing the scattering data [12,13] are limited to second-order convergence. In the study of Hald [14], a method for reconstructing the solution of an integrable equation from its scattering data was developed. This method is based on the trapezoidal rule and thus does not achieve spectral convergence. Furthermore, we take advantage of contour deformation to retain accuracy for large *x* and *t* unlike the approach in [14].

## 2. Integrability and Riemann–Hilbert problems

Both the focusing and defocusing NLS equations are completely integrable [15], p. 110. This means that for each equation there exist two linear systems of ordinary differential equations (ODEs) that depend on a spectral parameter *k*,
such that *μ*_{tx}=*μ*_{xt} if and only if *q* satisfies the given equation. We refer to this system of equations as the Lax pair or the scattering problem associated with the equation. Lax pairs for both the focusing and defocusing NLS equations are found through the modified Zakharov–Shabat scattering problem given by
2.1
Here *q*,*r*,*A*,*B*,*C* and *D* are scalar functions to be determined [6] (see also [15], p. 44). Choosing,
2.2
and , we obtain Lax pairs for both the focusing and defocusing NLS equations.

The analyticity properties in *k* of the solutions of (2.1) are studied. For simplicity, we assume where
Piecewise-smooth initial data with exponential decay can also be treated but this is beyond the scope of this paper. We define two matrix solutions of (2.1) by their corresponding asymptotic behaviour:
2.3
Lioville’s formula implies that the determinants of these solutions are constant in *x*. Since the determinants are non-zero at , these matrix solutions define a linearly independent set which spans the solution space. There must exist a transition matrix
such that *μ*^{+}(*x*;*k*)=*μ*^{−}(*x*;*k*)*T*(*k*). Define *ρ*(*k*)=*b*(*k*)/*a*(*k*) to be the reflection coefficient. Some symmetry properties of *T* follow [15]:
2.4

For the focusing NLS equation (*λ*=1), there may exist values and Im *κ*_{j}>0 such that
where the subscripts refer to columns. This implies that decays at both , exponentially. Thus, is an eigenfunction of (2.1). From the symmetries (2.4), is also an eigenvalue. For these values of *k*, we define the norming constants
For the defocusing NLS equation (*λ*=−1), it is known that there are no such eigenfunctions [15]. This implies that there are no smooth soliton solutions with spatial decay for the defocusing NLS equation. We define the set
2.5
to be the scattering data, noting that the sets of eigenvalues and norming constants could be empty. We call the process of finding the scattering data *direct scattering*.

### Remark 2.1

We assume for *δ*>0. This allows *T*(*k*) to be analytically extended to a neighbourhood of the real line. But it may happen that *T*(*k*) (in particular *b*(*k*)) cannot be extended above *κ*_{j} and thus *b*(*κ*_{j}) is really an abuse of notation and should be replaced with a constant, *b*_{j}.

The solutions *μ*^{±} can be grouped and transformed in such a way that they satisfy a Riemann–Hilbert problem (RHP) [15]. Loosely speaking, an RHP is a boundary-value problem in the complex plane:

### Problem 2.1

[16] Given an oriented contour and a jump matrix , find a bounded function which is analytic everywhere in the complex plane except on *Γ*, such that
2.6
where *Φ*^{+} denotes the limit of *Φ* as *k* approaches *Γ* from the left, *Φ*^{−} denotes the limit of *Φ* as *k* approaches *Γ* from the right and . We use [*G*;*Γ*] to denote this RHP.

The RHP associated with the NLS equations is
2.7
with the residue conditions (when *λ*=1),
and

Once the solution of the RHP is known the solution to the corresponding NLS equations is given by the expression
where the subscript denotes the (1,2) component of the matrix. We call the process of solving the RHP and reconstructing the solution *inverse scattering*.

We follow the standard procedure [9] to turn the residue conditions into jump conditions. Fix 0<*ϵ*, so that the circles and do not intersect each other or the real axis. We define by
2.8
It is straightforward to show that solves the RHP:
where *A*^{−}_{j}(*A*^{+}_{j}) has clockwise(counter-clockwise) orientation.

We end this section with some final remarks on RHP. Given an oriented contour *Γ*, we define the Cauchy integral
It is well known that the operators defined by
are bounded operators from *L*^{2}(*Γ*) to itself when *Γ* is Lipschitz [17]. All contours we consider are piecewise smooth with transverse intersections and the boundedness of the operator follows for such contours. Furthermore, these operators satisfy the identity
2.9
If we assume the solution to an RHP is of the form , we can substitute this into *Φ*^{+}=*Φ*^{−}*G* and use (2.9) to obtain
2.10
This is a singular integral equation (SIE) for *u*. We use to denote the operator in (2.10). This motivates the following definition.

### Definition 2.1

*An RHP [G;Γ] is said to be well-posed if* *is boundedly invertible on L*^{2}*(Γ) and G−I∈L*^{2}*(Γ).*

In addition to , we assume that the set {*κ*_{j}} is bounded away from the real line. This is sufficient to ensure that all RHPs we address are well-posed.

## 3. Numerical direct scattering

We describe a procedure to compute the scattering data (2.5). This follows [9]. We look for solutions of the form (2.3) to (2.1). Define
and two functions
3.1
Therefore, as and as . Rewriting (2.1),
3.2
and we find that both *K* and *J* solve
For each *k*, this can be solved with a Chebyshev collocation method on for *J* and on for *K* using the appropriate boundary condition at infinity. (We conformally map the unbounded domains to the unit interval using (1+*x*)/(1−*x*) and (*x*−1)/(1+*x*).) If we use *n* collocation points, this gives two approximate solutions *J*_{n} and *K*_{n} for *J* and *K*, respectively. From *J*_{n} and *K*_{n}, we obtain *μ*^{−}_{n} and *μ*^{+}_{n}, approximations of *μ*^{−} and *μ*^{+}, respectively, by inverting (3.1). Furthermore, *μ*^{−}_{n} and *μ*^{+}_{n} share the point *x*=0 in their domain of definition. Define
This is an approximation of the transition matrix, from which we extract an approximation of the reflection coefficient.

### Remark 3.1

Since *q*_{0}(*x*) decays exponentially, we benefit from solving the ODEs on [−*L*,0] and [0,*L*] for *L* sufficiently large. In this case, solving the equations on unbounded domains provides a good check to ensure *L* is sufficiently large. It can actually be shown that the resulting approximation after this truncation is an entire function of *k*.

This procedure works well for *k* in a neighbourhood of the real line. The solutions which decay at both are all that is needed to obtain *b*(*k*) when *a*(*k*)=0. Furthermore, from the analyticity properties of *a* we have [15], p. 75
Thus, knowing *a*(*k*) on the real line and *b*(*k*) when *a*(*k*)=0 allows us to compute *C*_{j}=*b*(*κ*_{j})/*a*′(*κ*_{j}). In practice, we use the framework [18] that will be discussed in §4 to compute these Cauchy integrals. Also, *a*′(*k*) can be obtained accurately using spectral differentiation by mapping the real line to the unit circle.

The remaining problem is that of computing *κ*_{j}. We consider (3.2)
3.3
We make the change of variables , , and we obtain
We use Hill’s method [19] to compute the eigenvalues of the operator
3.4
in the space *L*^{2}([−*π*,*π*]). Following the arguments in Trogdon *et al*. [9] and using the convergence of Hill’s method, even for non-self adjoint operators [20,21], the only eigenvalues we obtain off the real line are those associated with (3.3). This allows us to compute with spectral accuracy. We can always test to make sure we have captured all eigenvalues by computing the IST at *t*=0 and checking that we recover the initial condition. The method has not been developed to deal singular limits in which lots of eigenvalues are present. There is in principle no limit on card{*κ*_{j}} but the computational cost grows with the size of this set. Future work will be in the direction of numerically investigating .

### (a) Numerical results

In this section, we present numerical results for direct scattering. First, we compare our method with a reflection coefficient that is known analytically. Next, we present numerically computed reflection coefficients, which we use later.

For the focusing NLS equation (*λ*=1), the authors Tovbis *et al*. [22] present an explicit reflection coefficient for initial conditions of the form
3.5
The components of the transition matrix take the form
where
Here *Γ* is the Gamma function [23]. The set {*κ*_{j}} is non-empty for 0≤*μ*<2. Its elements are
In figure 1, we plot the reflection coefficient for *A*=1 and *μ*=0.1. These plots demonstrate spectral convergence.

In figure 2, we show the computed reflection coefficient for for both focusing and defocusing NLS. For focusing NLS, we find *κ*_{1}=−0.5+1.11151i, see figure 2*c* for a plot of the spectrum of (3.4) found using Hill’s method.

## 4. Numerical inverse scattering

Numerical inverse scattering has two major components. The first is the use of a Chebyshev collocation method for solving RHPs [18] and the second is the deformation of contours in the spirit of the method of nonlinear steepest descent (NSD) [24]. The use of NSD is essential since the jump for the RHP (2.7) is oscillatory for large values of *x* and *t*. An interesting and desired consequence of using NSD is that the resulting numerical method is accurate for large values of *x* and *t*. More specifically, the computational cost to compute the solution at a given point, accurate to within a given tolerance, is shown to be independent of *x* and *t* (see §7). First, we describe the collocation method for RHPs and then deform the RHP (2.7) to produce an efficient numerical method.

Consider the contour where each *Γ*_{j} is a non-self-intersecting arc. We use *γ*_{0} to denote the set of self-intersections of *Γ*. Assume we have a sequence of Möbius transformations *M*_{1},…,*M*_{n} such that *M*_{k}([−1,1])=*Γ*_{k}. Define , the Chebyshev points, and let *T*_{m}(*x*) denote the *m*th Chebyshev polynomial of the first kind. Given an RHP
4.1
which satisfies compatibility conditions at the self-intersection points of *Γ*, the framework in [18] will generally return a vector *W*_{j} of values at the mapped points , so that the function defined piecewise by
satisfies

—

*I*+𝒞_{Γ}*W*is bounded, and—

*I*+𝒞_{Γ}*W*satisfies the RHP (4.1) exactly at .

All RHPs considered here satisfy the conditions needed to apply the framework in [18].

We describe briefly the approach of this framework. A closed-form expression [25] for the Cauchy transform of the basis allows the discretization of (2.10) by evaluating the Cauchy transform of the basis at the points . A modified definition for the Cauchy transform is required at the junction points *γ*_{0} (which are included in the collocation points ), at which the Cauchy transform of this basis is unbounded. By assuming that the computed *W* is in the class of functions for which *I*+𝒞_{Γ}*v* is bounded, we define the bounded contribution of the Cauchy transform of each basis element at the points *γ*_{0}. A consequence of solving the resulting linear system is that the numerically calculated *W* must be in this class of functions. Therefore, *I*+𝒞_{Γ}*W* is bounded and satisfies the RHP at *γ*_{0}, hence at all points in .

If and *W*∈*L*^{1}(*Γ*) then
4.2
by dominated convergence. The integral on the right-hand side is computed using Clenshaw–Curtis quadrature. This relationship is needed in what follows to reconstruct the solution of an NLS equation from the solution of the RHP. The framework in [18,25] gives an efficient method for computing off *Γ* as well.

The spectral convergence of this method is guaranteed provided that the norm of inverse of the discretized operator does not grow too rapidly [10,18]. We make the following assumption.

### Assumption 4.1

*If* *G* *is continuous and* [*G*;*Γ*] *is well-posed then we assume*
*where* ∥⋅∥_{2} *is the induced operator norm on* *L*^{2}(*Γ*) *and* ∥⋅∥_{n} *is an appropriate finite-dimensional operator norm* [10], §5.2.

### Remark 4.1

As previously mentioned, the jump matrix *G* needs to satisfy a compatibility condition at each point in *γ*_{0}. The description of this is beyond the scope of this paper, see [10] for details.

Furthermore, the results in [10] justify the truncation of contours when they are, to machine precision, the identity matrix. This allows us to consider all RHPs on appropriate contours of finite length.

### Remark 4.2

A Mathematica implementation of the framework in [18] is available in the electronic supplementary material [26].

### (a) Small time

When both *x* and *t* are small, the RHP needs no deformation. When *t* is small, but *x* is large, the RHP needs to be deformed. We introduce factorizations of the jump matrix. Define
4.3
Note that *G*(*k*)=*L*(*k*)*D*(*k*)*U*(*k*)=*M*(*k*)*P*(*k*). We assume the sets and are empty. Later make the proper modifications to incorporate the extra contours required otherwise.

We wish to deform contours of (2.7) off the real line so that oscillations are turned to exponential decay. The matrix *G* contains the two factors and if one decays the other must grow. This motivates separating these factors using the process of *lensing* [24]. Suppose we can factor the jump function as *ABC*. It is possible to separate the single contour into three contours, as in figure 3, assuming *A* and *C* are analytic in a sufficiently large neighbourhood of the original contour and continuous up to the new contour. If satisfies the jumps on the split contour, it is clear that we can recover *Φ* by defining between the top contour and the original contour, between the original contour and the bottom contour and everywhere else, see the bottom of figure 3. It is important to note that the limiting properties of *Φ* and at infinity are the same.

Since for some *δ*>0 it follows that *ρ* is analytic in some strip for *γ*>0. This can be seen by considering the Volterra integral equations for the eigenfunctions *μ*^{±}. The factors in (4.3) allow lensing but we need to determine where to lens. We look for stationary points of the oscillator: *θ*′(*k*)=0 when *k*=*k*_{0}=−*x*/(4*t*). We use the *LDU* factorization for *k*<*k*_{0} and *MP* for *k*>*k*_{0}. See figure 4 for this deformation and note that the contours are locally deformed in the direction of steepest descent, that is, the direction along which the jump matrix tends to the identity matrix most rapidly. We denote the solution of this lensed RHP by *Φ*_{1}.

This RHP can be solved numerically provided *t* is not too large. As *t* increases, the solution *v* on the contour is increasingly oscillatory and is not well resolved using Chebyshev polynomials.

### (b) Long time

Next, we provide a deformation that leads to a numerical method that is accurate for arbitrarily large *x* and *t*. In the light of the results in [10], we need to remove the contour *D* from the RHP so that all jumps decay to the identity matrix away from *k*_{0} as *x* and *t* become large. The matrix-valued function
satisfies

We multiply the solution of the RHP in figure 4 by *Δ*^{−1} to remove the jump. To compute the solution to the new RHP, we use
Indeed, we see that if *J*=*D* there is no jump. Define *Φ*_{2}=*Φ*_{1}*Δ*^{−1}. See figure 5 for a schematic overview of the new jumps. This deformation is not sufficient for a numerical solution as *Δ*(*k*;*k*_{0}) has a singularity at *k*=*k*_{0}. We must perform one last deformation in a neighbourhood of *k*=*k*_{0} to bound contours away from this singularity. We use a piecewise definition of a function *Φ*_{3} (figure 6*a*) and compute the jumps (figure 6*b*). This is the final RHP. It is used, after contour truncation and scaling, to compute solutions of the NLS equations for arbitrarily large time. We discuss the scaling of the RHP in more detail in §7.

We use a combination of the deformation in figure 4 for small time and the deformation in figure 6*b* to obtain an approximation to focusing or defocusing NLS when no solitons are present in the solution. Lastly, we deal with the addition of solitons for the case of focusing NLS. There are additional jumps of the form
We assume that Re(*κ*_{j})>*γ*.

For small *x* and *t*, is close to unity and the contours and jumps need to be added to one of the deformations discussed above. This will not be the case for all *x* and *t*. When |*C*_{j} e^{θ(κj)}|>1, we invert this factor through a deformation. Define the set *K*_{x,t}={*j*:|*C*_{j} *e*^{θ(κj)}|>1}. Note that the *x* and *t* dependence enters through *θ*(*κ*_{j}). Next, define the functions

Define the piecewise-analytic matrix-valued function :
Computing the jumps that satisfies we find, for *j*∈*K*_{x,t},
This turns growth of the exponential to decay to the identity matrix. To simplify notation, we define
4.4
4.5

In figure 7, we present the full small-time and long-time RHPs. We use the notation [*J*;*Σ*] to denote the RHP in figure 7*b*.

### (c) Numerical results

In figure 8, we plot the solution of the focusing NLS equation with *q*_{0} given by (3.5) with *A*=1 and *μ*=0.1. The solution is nearly reflectionless but figure 8*d* shows the important dispersive aspect of the solution. Traditional numerical methods can fail to resolve this (figure 9). In figure 10, we plot the initial condition . The solutions of the focusing and defocusing NLS equations with this initial condition are computed. See figure 11 for focusing and figure 12 for defocusing. We also note that when the initial condition is less localized the corresponding reflection coefficient is more oscillatory. This makes it more difficult to resolve the solution of the corresponding RHP. We restrict ourselves to initial data with rapid fall-off for this reason, i.e. in numerical examples we consider with *δ* large.

We show numerical results that demonstrate spectral convergence. Let *q*_{0} be given by (3.5) with *A*=1 and *μ*=0.1 so that we can assume the reflection coefficient is computed to machine precision, i.e. *n*>80 in figure 1*b*. Define *q*(*n*,*x*,*t*) to be the approximate solution such that the number of collocation points per contour is proportional to *n*. In practice, we set the number of collocation points to be *n* on shorter contours, like all contours in figure 7*b*. For larger contours, like the horizontal contours in figure 7*b*, we use 5*n* collocation points. To analyse the error, we define
4.6
Using this notation, see figure 9 for a demonstration of spectral (Cauchy) convergence. Note that we choose *x* and *t* values to demonstrate spectral convergence in both the small-time and large-time regimes.

## 5. Extension to homogeneous Robin boundary conditions on the half line

The results thus far have been for the solution of the NLS equation posed on the whole line. We switch our attention to boundary-value problems on the half line, *x*≥0. Specifically, we extend the previous method to solve the following boundary-value problem:
5.1
Here is the space of smooth functions *f* on such that
If we take *α*=0, we obtain a Neumann problem and the limit effectively produces a Dirichlet problem. A method of images approach can be used to solve this problem. The approach of Biondini & Bui [27], first introduced by Bikbaev & Tarasov [28], takes the given initial condition on and produces an extension to using a Darboux transformation. For Neumann boundary conditions, this results in an even extension and for Dirichlet boundary conditions the transformation produces an odd extension. Consider the system of ODEs
5.2
with initial conditions
and the function
It was shown in [27] that the solution of the Cauchy problem for NLS equation on with initial data , restricted to , is the unique solution of (5.1). To compute the extended initial data , we first solve the system (5.2) numerically using a combination of Runge–Kutta 4 and 5. This is implemented in NDSolve in Mathematica. The IST for the extended potential can be used in the framework of §4 to numerically solve (5.1).

### Remark 5.1

The method of Bikbaev & Tarasov [28] was used to derive asymptotics by Deift & Park [29]. Another approach would be to use the method of Fokas to compute solutions [30,31].

### (a) Numerical results

In this section, we show numerical results for a Robin problem and a Neumann problem. As noted above, we could treat the Dirichlet problem by using an odd extension of our initial condition.

—

*Robin boundary conditions.*Here we show results for the case of the focusing NLS equation (*λ*=1) with*α*=−1 and with initial condition . Note that the initial condition satisfies the boundary condition at*t*=0. In figure 13*a*, we give the extended initial condition , and in figure 13*b*, we show the corresponding reflection coefficient. For this extended initial condition, we have four poles on the imaginary axis in the RHP which corresponds to two stationary solitons:—

*Neumann boundary conditions.*To show the reflection of a soliton off the boundary at*x*=0, we solve a Neumann problem (*α*=0) with initial condition (figure 14). The extension of the initial condition can be seen in figure 15*a*. In this case, it is just the even extension. The scattering data consist of This shows that we have a pair of poles in the RHP to the right of the imaginary axis and two to left. This corresponds to one soliton moving to the left and one soliton moving to the right. The reflection coefficient is shown in figure 15*b*. See figure 16 for the evolution of the solution up to*t*=7.

## 6. Singular solutions

As mentioned above, the defocusing NLS equation does not have soliton solutions that decay at infinity. We can insert the contours (see (2.8)) into the RHP anyway. We introduce *λ* into (4.4) to obtain appropriate jump conditions for the defocusing NLS equations. Define
When *λ*=1 (focusing), this definition agrees with (4.4).

When *λ*=−1 (defocusing), we investigate how these additional contours will manifest themselves in the solution. Consider
which is the one soliton solution of the focusing NLS equation [32]. A simple calculation shows that
is a solution of defocusing NLS. We are using the term solution loosely since this function has a pole when 2*η*4*tξ*+*x*−*x*_{0}=0. We call this a singular solution or singular soliton. These solutions are also called *positons*. See [33] for a deeper discussion of these solutions with applications to rogue waves.

What we obtain when adding the above contours to the RHP associated with the defocusing NLS equation is a nonlinear combination of these singular solutions in the presence of dispersion, as in the focusing case where the soliton was non-singular. See figure 17 for plots of a solution obtained using the reflection coefficient in figure 2 along with
This corresponds to two of these singular solitons moving towards each other, until they bounce off each other (the poles never cross). They interact in the presence of dispersion. We choose large and small norming constants to have the solitons away from *x*=0 when *t*=0. It is not possible to obtain these types of solutions with traditional numerical methods. Not surprisingly, the relative accuracy of our numerical approximation breaks down near the poles. For points bounded away from the poles, we still expect uniform convergence as discussed in §7.

## 7. Accuracy uniform in *x* and *t*

In this section, we describe how to use the theory of Olver & Trogdon [10] to prove the accuracy of the numerical method for arbitrarily large time for [*J*;*Σ*] in figure 7*b*. We assume the contours of the RHP are truncated according to §4. Define and *Σ*_{2}=*Σ*∖*Σ*_{1}. Define the restrictions of *J*:
We introduce some *x*- and *t*-independent domains *Ω*_{1} and *Ω*_{2}:
We use the change of variables
7.1
and the notation . Fix the trajectory in the (*x*,*t*) plane: *x*=*c*4*t*. We wish to apply the theory of Olver & Trogdon [10]. To do this, we need to modify our solution procedure. First, we numerically solve with *n* collocation points to obtain a solution . The change of variables (7.1) is inverted, defining
Define . Then is solved numerically with *n* collocation points (for simplicity) to obtain a function *Φ*_{2,n}. Note that there is no change of variables to invert for this RHP. The function *Φ*_{n}=*Φ*_{1,n}*Φ*_{2,n} is an approximation of the solution to full RHP [*J*;*Σ*].

Since the arc length of *Σ*_{1} tends to zero for large *t* RHPs and [*J*_{2};*Ω*_{2}] decouple asymptotically. The conditions we check are [10], proposition 5.2:

— exists and is uniformly bounded in

*t*;— exists and is uniformly bounded in

*t*; and— all derivatives of and

*J*_{2}(*z*), in the*z*variable, are uniformly bounded in*t*and*z*.

We defer the details of proving each of these statements to the electronic supplementary material and state the result. Define *W*^{i} to be the exact solution of the SIE on *Ω*_{i}, see (2.10) and to be its approximation with *n* collocation points per contour. We arrive at the following theorem.

### Theorem 7.1

Fix *T*>0 and *C*>0. For every *ϵ*>0, there exists *N*>0 such that for *n*>*N* and *t*>*T*

Since the arc length of the contours is bounded we have uniform *L*^{1} convergence and (4.2) demonstrates that *q*(*x*,*t*) is approximated uniformly.

### Remark 7.1

This theorem relies heavily on assumption 4.1. If assumption 4.1 fails to be true the numerical method may not converge.

### Remark 7.2

The results of [34,35] only apply to the defocusing NLS equation. The difficulty of the focusing NLS equation is the lack of information about {*κ*_{j}} and possible accumulation points on the real line [11]. We are considering cases where we assume {*κ*_{j}} are known (again see [11]) and the analysis proceeds in a way similar to focusing NLS.

### Remark 7.3

Despite the fact that the theorem restricts to |*k*_{0}|<*C*, we still obtain uniform convergence. If , for every *ϵ*>0 there exists *C*>0 such that for |*k*_{0}|>*C*, |*q*(*x*,*t*)|<*ϵ* [34]. Thus, we approximate the solution with zero when |*k*_{0}|>*C*. The domain of accurate and non-trivial approximation {*x*:|*x*/(4*t*)|≤*C*} grows as *t* increases.

### (a) Numerical results

To demonstrate asymptotic accuracy, we use the notation in (4.6) and fix *n* and *m*. We let *x* and *t* become large along a specific trajectory. For our purposes, we use *x*=4*t*. Note that along this trajectory *q* is on the order of [36] (see also [34,35,37]). This allows us to estimate the relative error. See figure 18 for a demonstration of the accuracy of the method for large *x* and *t*. We see that the relative error is bounded and small using relatively few collocation points.

## 8. Conclusions

We are able to successfully compute the full IST on the whole line for smooth initial data with sufficiently fast decay and compute singular solutions if we use special spectral data. All spectral data can be obtained with spectral accuracy. Furthermore, numerical extensions of data on the half line can be used to solve homogeneous Dirichlet, Neumann and Robin boundary-value problems on the half line. After deforming the contours of the RHP in the spirit of Deift and Zhou and the method of NSD, we obtain a method that is accurate for arbitrarily large values of *x* and *t*. More precisely, the computational cost required to compute the solution accurate to within a given tolerance can be shown to be independent of *x* and *t*.

## Acknowledgements

We thank Bernard Deconinck for many useful conversations and Gino Biondini for suggesting the numerical solution of boundary-value problems. We also thank the referees for their input as it resulted in a more comprehensive paper. We acknowledge the National Science Foundation for its generous support through grant NSF-DMS-1008001 (T.T.). Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the funding sources.

- Received May 30, 2012.
- Accepted October 19, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.