Fast Inverse Nonlinear Fourier Transformation using Exponential One-Step Methods, Part I: Darboux Transformation

TL;DR

This paper introduces a fast, FFT-based numerical framework for the nonlinear Fourier transform and Darboux transformation, significantly improving computational efficiency and accuracy for soliton-based scattering problems.

Contribution

It develops a unified exponential one-step method framework for the nonlinear Fourier transform and proposes a fast Darboux transformation algorithm with superior complexity.

Findings

The FDT algorithm has complexity $ ext{O}(KN + N ext{log}^2 N)$.

The error in the $N$-sample $K$-soliton computation decreases as $ ext{O}(N^{-p})$.

The proposed methods outperform classical algorithms in efficiency and accuracy.

Abstract

This paper considers the non-Hermitian Zakharov-Shabat (ZS) scattering problem which forms the basis for defining the SU$(2)$-nonlinear Fourier transformation (NFT). The theoretical underpinnings of this generalization of the conventional Fourier transformation is quite well established in the Ablowitz-Kaup-Newell-Segur (AKNS) formalism; however, efficient numerical algorithms that could be employed in practical applications are still unavailable. In this paper, we present a unified framework for the forward and inverse NFT using exponential one-step methods which are amenable to FFT-based fast polynomial arithmetic. Within this discrete framework, we propose a fast Darboux transformation (FDT) algorithm having an operational complexity of $\mathscr{O}\left(KN+N\log^2N\right)$ such that the error in the computed $N$-samples of the $K$-soliton vanishes as $\mathscr{O}\left(N^{-p}\right)$ where $p$ is the order of convergence of the underlying one-step method. For fixed $N$, this algorithm outperforms the the classical DT (CDT) algorithm which has a complexity of $\mathscr{O}\left(K^2N\right)$. We further present extension of these algorithms to the general version of DT which allows one to add solitons to arbitrary profiles that are admissible as scattering potentials in the ZS-problem. The general CDT/FDT algorithms have the same operational complexity as that of the $K$-soliton case and the order of convergence matches that of the underlying one-step method. A comparative study of these algorithms is presented through exhaustive numerical tests.

Figures

Equations

$${\sigma }_1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\sigma }_2=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right),{\sigma }_3=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right),$$

$$\mathbf{f}(\tau )=\frac{1}{2\pi }{\int }_{\Gamma }F(\zeta )e^{-i\zeta \tau }d\zeta ,$$

$$i{\partial }_{Z}q={\partial }_T^{2}q+2|q{|}^{2}q,(T,Z)\in \mathbb{R}\times {\mathbb{R}}_{+},$$

$$|\rho (\xi )|\le \frac{C}{1+|\xi |},$$

$$i{q}_t=q_{xx}+2|q{|}^{2}q,(x,t)\in \mathbb{R}\times {\mathbb{R}}_{+},$$

$${\mathbf{v}}_x=-i\zeta {\sigma }_3\mathbf{v}+U\mathbf{v},$$

$${\mathbf{v}}_t=2i{\zeta }^2{\sigma }_3\mathbf{v}+[-2\zeta U+i{\sigma }_3(U^2-U_x)]\mathbf{v},$$

$$U=\left(\begin{array}{cc}0& q(x,t)\\ r(x,t)& 0\end{array}\right),r(x,t)=-q^{\ast }(x,t),$$

$$\begin{array}{rl}\mathbf{\varphi }(x;\zeta )& =a(\zeta )\overline{\mathbf{\psi }}(x;\zeta )+b(\zeta )\mathbf{\psi }(x;\zeta ),\\ \overline{\mathbf{\varphi }}(x;\zeta )& =-\overline{a}(\zeta )\mathbf{\psi }(x;\zeta )+\overline{b}(\zeta )\overline{\mathbf{\psi }}(x;\zeta ).\end{array}$$

$$\begin{array}{rl}\mathbf{\psi }(x;\zeta )& =-a(\zeta )\overline{\mathbf{\varphi }}(x;\zeta )+\overline{b}(\zeta )\mathbf{\varphi }(x;\zeta ),\\ \overline{\mathbf{\psi }}(x;\zeta )& =\overline{a}(\zeta )\mathbf{\varphi }(x;\zeta )+b(\zeta )\overline{\mathbf{\varphi }}(x;\zeta ).\end{array}$$

$$W\left(\mathbf{u},\mathbf{v}\right)=\left(\mathbf{u},\mathbf{v}\right)=u_1{v}_2-v_1{u}_2,$$

$$\begin{array}{rl}& a(\zeta )=W\left(\mathbf{\varphi },\mathbf{\psi }\right),b(\zeta )=W\left(\overline{\mathbf{\psi }},\mathbf{\varphi }\right),\\ & \overline{a}(\zeta )=W\left(\overline{\mathbf{\varphi }},\overline{\mathbf{\psi }}\right),\overline{b}(\zeta )=W\left(\overline{\mathbf{\varphi }},\mathbf{\psi }\right).\end{array}$$

$$\begin{array}{r}\overline{\mathbf{\psi }}(x;\zeta )=i{\sigma }_2{\mathbf{\psi }}^{\ast }(x;{\zeta }^{\ast })=\left(\begin{array}{c}{\psi }_2^{\ast }(x;{\zeta }^{\ast })\\ -{\psi }_1^{\ast }(x;{\zeta }^{\ast })\end{array}\right),\\ \overline{\mathbf{\varphi }}(x;\zeta )=i{\sigma }_2{\mathbf{\varphi }}^{\ast }(x;{\zeta }^{\ast })=\left(\begin{array}{c}{\varphi }_2^{\ast }(x;{\zeta }^{\ast })\\ -{\varphi }_1^{\ast }(x;{\zeta }^{\ast })\end{array}\right),\end{array}$$

$${\mathfrak{S}}_K=\{({\zeta }_k,b_k)\in {\mathbb{C}}^2|\mathrm{Im}{\zeta }_k>0,k=1,2,\dots ,K\}.$$

$$v(x,t;\zeta )=(\mathbf{\varphi },\mathbf{\psi })=\left(\begin{array}{cc}{\varphi }_1& {\psi }_1\\ {\varphi }_2& {\psi }_2\end{array}\right).$$

$$v_K(x,t;\zeta )={\mu }_K(\zeta )D_K(x,t;\zeta ,{\mathfrak{S}}_K)v_0(x,t;\zeta ),\zeta \in {\overline{\mathbb{C}}}_{+},$$

$$D_K(x,t;\zeta ,{\mathfrak{S}}_K)=\sum _{k=0}^K{D}_k^{(K)}(x,t;{\mathfrak{S}}_K){\zeta }^k,$$

$$D_k^{(K)}=\left(\begin{array}{cc}{d}_0^{(k,K)}& d_1^{(k,K)}\\ -d_1^{(k,K)\ast }& d_0^{(k,K)\ast }\end{array}\right),k=0,1,\dots ,K-1.$$

$$\begin{array}{rl}{a}_K(\zeta )& =\det \left[v_K(x,t;\zeta )\right]\\ & =[{\mu }_K(\zeta ){]}^2\det \left[D_K(x,t;\zeta ,{\mathfrak{S}}_K)\right]a_0(\zeta ).\end{array}$$

$$\det \left[D_K(x,t;\zeta ,{\mathfrak{S}}_K)\right]=\prod _{k=1}^K(\zeta -{\zeta }_k)(\zeta -{\zeta }_k^{\ast }),$$

$$a_K(\zeta )=a_0(\zeta )\prod _{k=1}^K\left(\frac{\zeta -{\zeta }_k}{\zeta -{\zeta }_k^{\ast }}\right),$$

$${\mu }_K(\zeta )=\prod _{k=1}^K\frac{1}{(\zeta -{\zeta }_k^{\ast })}.$$

$$D_K(x,t;{\zeta }_k,{\mathfrak{S}}_K)[{\mathbf{\varphi }}_0(x,t;{\zeta }_k)-b_k(t){\mathbf{\psi }}_0(x,t;{\zeta }_k)]=0.$$

$$[D_K{v}_0{]}_x-(-i\zeta {\sigma }_3+U_K)D_K{v}_0=0,$$

$$[{\partial }_x{D}_K-(-i\zeta {\sigma }_3+U_K)D_K+D_K(-i\zeta {\sigma }_3+U_0)]v_0=0.$$

$$[{\partial }_x{D}_K-(-i\zeta {\sigma }_3+U_K)D_K+D_K(-i\zeta {\sigma }_3+U_0)]=0.$$

$$\begin{array}{rl}{U}_K& =U_0+i[{\sigma }_3,D_{K-1}^{(K)}]\\ & =U_0+\left(\begin{array}{cc}0& 2i{d}_1^{(K-1,K)}\\ 2i{d}_1^{(K-1,K)\ast }& 0\end{array}\right).\end{array}$$

$${\beta }_0(x,t;{\zeta }_1,b_1)=\frac{{\varphi }_1^{(0)}(x,t;{\zeta }_1)-b_1(t){\psi }_1^{(0)}(x,t;{\zeta }_1)}{{\varphi }_2^{(0)}(x,t;{\zeta }_1)-b_1(t){\psi }_2^{(0)}(x,t;{\zeta }_1)},$$

$$\begin{array}{rl}& D_1(x,t;\zeta ,{\mathfrak{S}}_1|{\mathfrak{S}}_0)\\ & =\zeta {\sigma }_0-\left(\begin{array}{cc}{\beta }_0& 1\\ 1& -{\beta }_0^{\ast }\end{array}\right)\left(\begin{array}{cc}{\zeta }_1& 0\\ 0& {\zeta }_1^{\ast }\end{array}\right){\left(\begin{array}{cc}{\beta }_0& 1\\ 1& -{\beta }_0^{\ast }\end{array}\right)}^{-1}\\ & =\zeta {\sigma }_0-\frac{1}{1+|{\beta }_0{|}^2}\left(\begin{array}{cc}|{\beta }_0{|}^2{\zeta }_1+{\zeta }_1^{\ast }& ({\zeta }_1-{\zeta }_1^{\ast }){\beta }_0\\ ({\zeta }_1-{\zeta }_1^{\ast }){\beta }_0^{\ast }& {\zeta }_1+{\zeta }_1^{\ast }|{\beta }_0{|}^2\end{array}\right).\end{array}$$

$$q_1(x,t)=q_0(x,t)-2i\frac{({\zeta }_1-{\zeta }_1^{\ast }){\beta }_0}{1+|{\beta }_0{|}^2}.$$

$$\begin{gathered} D_K(\zeta ,{\mathfrak{S}}_K|{\mathfrak{S}}_0)=D_1(\zeta ,{\mathfrak{S}}_K|{\mathfrak{S}}_{K-1}) \\ \times D_1(\zeta ,{\mathfrak{S}}_{K-1}|{\mathfrak{S}}_{K-2})\times \dots \times D_1(\zeta ,{\mathfrak{S}}_1|{\mathfrak{S}}_0), \end{gathered}$$

Full text

In this paper, we present a unified framework for the forward and inverse NFT using exponential one-step methods which are amenable to FFT-based fast polynomial arithmetic. Within this discrete framework, we propose a fast Darboux transformation (FDT) algorithm having an operational complexity of $\mathop{\mathscr{O}}\left(KN+N\log^{2}N\right)$ such that the error in the computed $N$-samples of the $K$-soliton vanishes as $\mathop{\mathscr{O}}\left(N^{-p}\right)$ where $p$ is the order of convergence of the underlying one-step method. For fixed $N$, this algorithm outperforms the the classical DT (CDT) algorithm which has a complexity of $\mathop{\mathscr{O}}\left(K^{2}N\right)$. We further present extension of these algorithms to the general version of DT which allows one to add solitons to arbitrary profiles that are admissible as scattering potentials in the ZS-problem. The general CDT/FDT algorithms have the same operational complexity as that of the $K$-soliton case and the order of convergence matches that of the underlying one-step method. A comparative study of these algorithms is presented through exhaustive numerical tests.

Notations

The set of real numbers (integers) is denoted by $\mathbb{R}$ ($\mathbb{Z}$) and the set of non-zero positive real numbers (integers) by $\mathbb{R}_{+}$ ($\mathbb{Z}_{+}$). The set of complex numbers are denoted by $\mathbb{C}$, and, for $\zeta\in\mathbb{C}$, $\operatorname{Re}(\zeta)$ and $\operatorname{Im}(\zeta)$ refer to the real and the imaginary parts of $\zeta$, respectively. The complex conjugate of $\zeta\in\mathbb{C}$ is denoted by $\zeta^{*}$ and $\sqrt{\zeta}$ denotes its square root with a positive real part. The upper-half (lower-half) of $\mathbb{C}$ is denoted by $\mathbb{C}_{+}$ ($\mathbb{C}_{-}$) and it closure by $\overline{\mathbb{C}}_{+}$ ($\overline{\mathbb{C}}_{-}$). The set $\mathbb{D}=\{z|\,z\in\mathbb{C},\,|z|<1\}$ denotes an open unit disk in $\mathbb{C}$ and $\overline{\mathbb{D}}$ denotes its closure. The set $\mathbb{T}=\{z|\,z\in\mathbb{C},\,|z|=1\}$ denotes the unit circle in $\mathbb{C}$. The Pauli’s spin matrices are denoted by, $\sigma_{j},\,j=1,2,3$, which are defined as

[TABLE]

where $i=\sqrt{-1}$. For uniformity of notations, we denote $\sigma_{0}=\text{diag}(1,1)$. Matrix transposition is denoted by $(\cdot)^{\intercal}$ and $I$ denotes the identity matrix. For any two vectors $\bm{u},\bm{v}\in\mathbb{C}^{2}$, $\operatorname{\mathscr{W}}(\bm{u},\bm{v})\equiv(u_{1}v_{2}-u_{2}v_{1})$ denotes the Wronskian of the two vectors and $[A,B]$ stands for the commutator of two matrices $A$ and $B$. Partial derivatives with respect to $x$ are denoted by $\partial_{x}$ or $(\cdot)_{x}$ while repeated derivatives by $\partial^{2}_{x}$. The support of a function $f:\Omega\rightarrow\mathbb{R}$ in $\Omega$ is defined as $\operatorname*{supp}f=\overline{\{x\in\Omega|\,f(x)\neq 0\}}$. The Lebesgue spaces of complex-valued functions defined in $\mathbb{R}$ are denoted by $\mathsf{L}^{p}$ for $1\leq p\leq\infty$ with their corresponding norm denoted by $\|\cdot\|_{\mathsf{L}^{p}}$ or $\|\cdot\|_{p}$.

The inverse Fourier-Laplace transform of a function $F(\zeta)$ analytic in $\overline{\mathbb{C}}_{+}$ is defined as

[TABLE]

where $\Gamma$ is any contour parallel to the real line.

I Introduction

This paper considers the two-component non-Hermitian scattering problem first studied by Zakharov and Shabat (ZS) Zakharov and Shabat (1972), which forms the basis for defining the SU$(2)$-nonlinear Fourier transformation (NFT). For certain integrable nonlinear equations whose general description is provided by the AKNS-formalism Ablowitz et al. (1974); Ablowitz and Segur (1981), the NFT offers a powerful means of solving the corresponding initial-value problem (IVP). One such example is the nonlinear Schrödinger equation (NSE) that is commonly used to model channels for optical fiber communication. The propagation of optical field in a loss-less single mode fiber under Kerr-type focusing nonlinearity is governed by the NSE Kodama and Hasegawa (1987); Agrawal (2013) which can be cast into the following standard form

[TABLE]

where $q(T,Z)$ is a complex valued function associated with the slowly varying envelope of the electric field, $Z\in\mathbb{R}_{+}$ is the position along the fiber and $T$ is the retarded time. This equation also provides a satisfactory description of optical pulse propagation in the guiding-center or path-averaged formulation Hasegawa and Kodama (1990, 1991); Turitsyn et al. (2012) when more general scenarios such as presence of fiber losses, lumped or distributed periodic amplification are included in the mathematical model of the physical channel. The IVP corresponding to (1) consists in finding the evolved field $q(T,Z)$ for a given initial condition $q(T,0)$ under vanishing boundary conditions.

For a given initial condition $q(T,0)$, the nonlinear Fourier spectrum consists of (i) a continuous part $\rho(\xi),\,\,\xi\in\mathbb{R}$, and, (ii) a discrete part given by $\mathfrak{S}_{K}=\{(\zeta_{k},b_{k})\in\mathbb{C}^{2}|\,\operatorname{Im}{\zeta_{k}}>0,\,k=1,2,\ldots,K\}$ which is an ordered pair of eigenvalues $\zeta_{k}$ and the respective norming constants $b_{k}$ (see Ablowitz et al. (1974) or Sec. II for a complete introduction). The discrete spectrum is associated with the solitonic components of the potential which will be referred to as bound states in the rest of this paper. The energy in these states does not disperse away as in the case of linear waves, a phenomenon which is adequately characterized by the term “bound states”. The evolution of the nonlinear Fourier (NF) spectrum depicted in Fig. 1 is reminiscent of evolution in a linear channel–a property which is attributed to the integrablity of the nonlinear channel.

In passing, we also note that the ZS-problem appears in various other physical systems, for instance, grating-assisted co-directional couplers (GACCs), a device used to couple light between two different guided modes of an optical fiber (see Feced and Zervas (2000); Brenne and Skaar (2003) and references therein), and, NMR spectroscopy where design of frequency selective pulses requires solution of a ZS-problem Rourke and Morris (1992a, b); Rourke and Saunders (1994).

Amongst the key physical effects that affect the performance of an optical fiber communication system, namely, chromatic dispersion, Kerr-type nonlinearity and optical noise, it is the latter two that have become the principle factors limiting the spectral efficiency of wavelength–division–multiplexed (WDM) networks at high signal powers. The reason behind this is largely the transmission methodologies that assume a linear model of the channel. The NF spectrum, in contrast, offers a novel way of encoding information in optical pulses where the nonlinear effects are adequately taken into account as opposed to being treated as a source of signal distortion. The idea to use discrete eigenvalues of the NF spectrum was first proposed by Hasegawa and Nyu Hasegawa and Nyu (1993) which they termed as the eigenvalue communication. Recently, Yousefi and Kschischang Yousefi and Kschischang (2014) have proposed nonlinear signal multiplexing in multi-user channels in order to mitigate the problem of nonlinear cross-talk that occurs in WDM systems. We note that the most general modulation technique uses both the discrete as well as the continuous part of the NF spectrum which was recently demonstrated in Aref et al. (2016). We refer the reader to a comprehensive review article Turitsyn et al. (2017) and the references therein for an overview of the progress in theoretical as well as experimental aspects of various NFT-based optical communication methodologies. It must be noted that practical implementation of NFT-based transmission is still quite far from becoming a reality Turitsyn et al. (2017) and there are other potential ways to combat nonlinear signal distortions in optical fibers Cartledge et al. (2017); Dar and Winzer (2017)111At this point, studies indicate that there is no clear winner as far as mitigation of impairments due to nonliearity in optical fibers is concerned; therefore, we continue our efforts to improve NFT-based approach. Also noteworthy is the fact that ZS-problem appears in various other systems of physical significance where this research can find application..

In any NFT-based modulation technique, the importance of low-complexity NFT-algorithms cannot be over emphasized. In this paper, we focus on the development of fast algorithms for various modulation scenarios of a NFT-based transmission system. As noted in Turitsyn et al. (2017), many of the existing numerical approaches tend to become inaccurate as the signal power increases. While this maybe attributed to lack of numerical precision, it could also be numerical ill-conditioning or as a result of naive implementation. It is difficult to fully address these problems in this work, but let us remark that stability and convergence of the numerical algorithm plays a key role in determining its performance in realistic scenarios. We discuss these two aspects quite rigorously in this work.

Our primary goal here is to provide a theoretical foundation for the algorithms reported in Vaibhav and Wahls (2017) where we also showcased our preliminary results demonstrating the first fast inverse NFT. The specific problems for which we seek fast algorithms in this work are as follows:

Problem I.1** (Generation of multi-solitons).**

Given an arbitrary discrete spectrum $\mathfrak{S}_{K}$ ($K$ being its cardinality or, in other words, the number of bound states), compute the corresponding multi-soliton potential.

Problem I.2** (Addition of bound states).**

Given an arbitrary potential $q_{\text{seed}}(x)$ referred to as the “seed” potential (assumed to be admissible as a scattering potential in the ZS-problem) and a given discrete spectrum $\mathfrak{S}_{K}$, compute the “augmented” potential such that its discrete spectrum is given by $\mathfrak{S}_{\text{aug.}}=\mathfrak{S}_{\text{seed}}\cup\mathfrak{S}_{K}$ where $\mathfrak{S}_{K}$ is known to be disjoint with $\mathfrak{S}_{\text{seed}}$, the discrete spectrum of the seed potential.

Problem I.3** (Inversion of continuous spectrum).**

Given an arbitrary continuous spectrum $\rho(\xi),\,\,\xi\in\mathbb{R}$, such that there exists a positive constant $C>0$ for which the estimate

[TABLE]

holds, compute the potential such that its continuous spectrum is $\rho(\xi)$ and the discrete spectrum is empty.

Problem I.4** (Inverse NFT).**

Given an arbitrary continuous spectrum $\rho(\xi),\,\,\xi\in\mathbb{R}$, satisfying the estimate in Prob. I.3 and a given discrete spectrum $\mathfrak{S}_{K}$, compute the potential such that its continuous spectrum is $\rho(\xi)$ and its discrete spectrum is $\mathfrak{S}_{K}$.

The first two of these problems can be solved, at least in principle, using the Darboux transformations (DT) Lin (1990); Gu et al. (2005). The Prob. I.1 can be solved with machine precision using DT with null-potential as the seed. The resulting complexity is $\mathop{\mathscr{O}}\left(K^{2}N\right)$ where $K$ is the number of eigenvalues and $N$ is the number of samples of the potential. The scenario in Prob. I.1 corresponds to the modulation of the discrete NF spectrum which has been explored by a number of groups Hari et al. (2014); Dong et al. (2015); Hari and Kschischang (2016) and it has also been experimentally demonstrated Terauchi and Maruta (2013); Bülow (2014, 2015); Matsuda et al. (2014); Aref et al. (2015); Bülow et al. (2016).

The Prob. I.2 cannot be solved without resorting to numerical methods for the ZS-problem because the so called Jost solutions (which are required in DT) are not known in a closed form for any arbitrary seed potential. Prob. I.3 and I.4 will be treated in a sequel to this paper.

The numerical techniques for solving Prob. I.1–I.4 developed in this work are based on exponential (linear) one-step methods Cox and Matthews (2002); Gautschi (2012) for the discretization of the ZS-problem. The method yields a discrete framework for solving the ZS-problem which resembles the transfer matrix approach for solving wave-propagation problems in dielectric layered media (Born and Wolf, 1999, Chap. 1). These transfer matrices have polynomial entries–a form that is amenable to the FFT-based polynomial arithmetic Henrici (1993) and is also compatible with the layer-peeling algorithm Bruckstein and Kailath (1987). All the methods considered in this article exhibit either a first order or a second order of convergence, i.e., the numerical errors vanish as $\mathop{\mathscr{O}}\left(N^{-p}\right)$ where $p$ is the order of the one-step method222The discrete system corresponding to Ablowitz-Ladik (AL) scattering problem (Ablowitz et al., 2004, Chap. 3) is also amenable to FFT-based fast polynomial arithmetic and satisfies the layer-peeling property Wahls and Poor (2015a); however, it does not illuminate on how to obtain a general recipe that could be applied to the ZS-problem in order to obtained a similar discrete system possessing a given order of convergence..

Within this discrete framework, we develop two algorithms: (a) the classical Darboux transform (CDT) which addresses Prob. I.2, and, (b) the fast Darboux transformation (FDT) which addresses Prob. I.1 and I.2 both333It is worth noting that an alternative fast method of solving Prob. I.1 is reported in Wahls and Poor (2015b) and it can be readily adapted to the discrete framework considered in this work. However, this method offers no control over the norming constants, therefore, we do not address this algorithm here.. The CDT algorithm is a direct numerical implementation of the DT in the continuum case where the seed Jost solutions are computed by numerically solving the scattering problem resulting in an overall complexity of $\mathop{\mathscr{O}}\left(K^{2}N\right)$. The FDT algorithm is entirely new and it is based on the pioneering work of Lubich on convolution quadrature Lubich (1988a, b, 1994). In order to ensure compatibility with Lubich’s construction, we restricted ourselves to the implicit Euler method and the trapezoidal rule. This algorithm has an operational complexity of $\mathop{\mathscr{O}}\left(N(K+\log^{2}N)\right)$ and an order of convergence that matches that of the underlying one-step method, i.e., $\mathop{\mathscr{O}}\left(N^{-p}\right)$ where $p=1$ (implicit Euler), $p=2$ (trapezoidal rule). With increasing number of eigenvalues, FDT clearly outperforms CDT. The numerical tests and error analysis of the numerical scheme suggests that CDT is useful only for smaller number of eigenvalues. These tests further reveal that FDT is not only more accurate for the general case, it also has superior numerical conditioning with increasing number of eigenvalues as opposed to the CDT algorithm which becomes unstable.

I.1 Outline of the paper

This paper is organized as follows: In Sec. II, we summarize the basic scattering theory and the Darboux transformation in the continuous regime.

The discrete scattering framework for the ZS-problem is developed in Sec. III where the numerical discretization in the spectral domain is described in Sec. III.1, and, properties of the numerical Jost solutions are discussed in Sec. III.2. We formulate the layer-peeling scheme in Sec. III.3 which is based on the discrete framework developed in Sec. III.1. Algorithmic aspects are addressed in Sec. III.4 and III.5 where we describe the sequential algorithm and its fast version obtained using a divide-and-conquer strategy, respectively. The sections III.6 to III.8 contain the main contribution of this paper: The method of inversion of continuous scattering coefficients using the Lubich’s method is discussed in Sec. III.6. In Sec. III.7, we apply the Lubich’s method to obtain the FDT algorithm for $K$-soliton potentials. Finally, the general version of the CDT algorithm and the FDT algorithm is discussed in Sec. III.8.

The benchmarking methods that used for comparison are discussed in Sec. IV. The necessary and sufficient condition for discrete inverse scattering is discussed in Sec. V. The stability and convergence analysis of the numerical schemes developed in earlier section is carried out in Sec. VI. The numerical experiments and results are discussed in the Sec. VII which is followed by Sec. VIII which concludes the paper.

II The AKNS System

In order to describe the fundamental basis of the nonlinear Fourier transform (NFT), we briefly review the scattering theory for a $2\times 2$ AKNS system corresponding to the NSE. Because the NSE shows up in various disciplines, we choose to present the theory in a form that is independent of the context and conform to the way it appears in the classical texts on the scattering theory. For a complex valued field $q(x,t)$, we will work with the standard form of NSE which reads as

[TABLE]

where $t>0$ is the evolution parameter identified as a time-like variable (this turns out to be the propagation distance $Z$ for the fiber model) and $x\in\mathbb{R}$ is the domain over which the field is defined (which is the retarded time $T$ for the fiber model). Henceforth, we closely follow the formalism developed in Ablowitz et al. (1974); Ablowitz and Segur (1981) for the exposition in this article. The NFT of the complex-valued field $q(x,t)$ is introduced via the associated Zakharov-Shabat scattering problem Zakharov and Shabat (1972) which can be stated as follows: Let $\zeta\in\mathbb{R}$ and $\bm{v}=(v_{1},v_{2})^{\intercal}\in\mathbb{C}^{2}$, then

[TABLE]

where

[TABLE]

is identified as the scattering potential. The second relation above corresponds to the focusing-type of nonlinearity for the NSE. The compatibility condition ($\bm{v}_{xt}=\bm{v}_{tx}$) between (3) and (4), assuming $\zeta$ is independent of $t$, produces the NSE as stated in (2).

The solution of the scattering problem (3), henceforth referred to as the ZS-problem, consists in finding the so called scattering coefficients which are defined through special solutions of (3) known as the Jost solutions described in the next subsection. These Jost solutions also play an important role in defining the Darboux transformation (DT) which is a powerful technique for constructing more complex potentials (as well as their Jost solutions) from simpler ones–this will be discussed in the final part of this section. There, we will be primarily interested in studying the form of DT which allows one to add bound states to a given potential.

II.1 Jost solutions

The Jost solutions are linearly independent solutions of (3) such that they have a plane-wave like behavior at $+\infty$ or $-\infty$. In the following, we set $t=0$ and suppress the time-dependence of the solutions for the sake of brevity.

• •

First kind: The Jost solutions of the first kind, denoted by $\bm{\psi}(x;\zeta)$ and $\overline{\bm{\psi}}(x;\zeta)$, are the linearly independent solutions of (3) which have the following asymptotic behavior as $x\rightarrow\infty$: $\bm{\psi}(x;\zeta)e^{-i\zeta x}\rightarrow(0,1)^{\intercal}$ and $\overline{\bm{\psi}}(x;\zeta)e^{i\zeta x}\rightarrow(1,0)^{\intercal}$.

• •

Second kind: The Jost solutions of the second kind, denoted by $\bm{\phi}(x,\zeta)$ and $\overline{\bm{\phi}}(x,\zeta)$, are the linearly independent solutions of (3) which have the following asymptotic behavior as $x\rightarrow-\infty$: $\bm{\phi}(x;\zeta)e^{i\zeta x}\rightarrow(1,0)^{\intercal}$ and $\overline{\bm{\phi}}(x;\zeta)e^{-i\zeta x}\rightarrow(0,-1)^{\intercal}$.

The evolution of the Jost solutions in time is governed by the equation (4) for $t\in\mathbb{R}_{+}$ under the asymptotic boundary conditions prescribed above. On account of the linear independence of $\bm{\psi}$ and $\overline{\bm{\psi}}$, we have

[TABLE]

Similarly, using the pair $\bm{\phi}$ and $\overline{\bm{\phi}}$, we have

[TABLE]

The coefficients appearing in the equations above can be written in terms of the Jost solutions by using the Wronskian relations444For any pair of linearly independent vectors, $\bm{v},\,\bm{u}\in\mathbb{C}^{2}$, their Wronskian which is defined as

$\operatorname{\mathscr{W}}\left(\bm{u},\bm{v}\right)=\left(\bm{u},\bm{v}\right)=u_{1}v_{2}-v_{1}u_{2},$

is non-zero. If $\bm{u},\bm{v}$ also qualify as Jost solutions, then their Wronskian is independent of $x$ Ablowitz et al. (1974).:

[TABLE]

These coefficients are known as the scattering coefficients and the process of computing them is referred to as forward scattering. As it turns out, we would also be interested in studying the analytic continuation of the Jost solutions with respect to $\zeta$, which in turn also determines the analytic continuation of the scattering coefficients. The motivation behind this is threefold: First, the inversion of the scattering coefficients cannot be done in general by knowing the value of the scattering coefficients over the real line (i.e. $\zeta\in\mathbb{R}$). Second, the knowledge of analyticity and decay properties of these functions in the complex plane allows us to establish certain theoretical estimates with greater ease. Lastly, in many cases, the knowledge of the analytic form introduces a certain redundancy in the system that can be exploited by the numerical algorithms to improve its numerical conditioning and stability.

In order to discuss the analytic continuation of the Jost solution with respect to $\zeta$, let us specify the following two classes of functions for the scattering potential (at $t=0$): Let $q(\cdot,0)\in\mathsf{L}^{1}$ such that $\operatorname*{supp}q(\cdot,0)\subset\Omega=[L_{1},L_{2}]$ or $|q(x,0)|\leq C\exp[-2d|x|]$ almost everywhere in $\mathbb{R}$ for some constants $C>0$ and $d>0$. In the former case, the Jost solutions have analytic continuation in whole of the complex plane with respect to $\zeta$. Consequently, the scattering coefficients $a(\zeta)$, $b(\zeta)$, $\overline{a}(\zeta)$, $\overline{b}(\zeta)$ are analytic functions of $\zeta\in\mathbb{C}$. In the latter case, the analyticity property can be summarized as follows (Ablowitz et al., 1974, Sec. IV.A): The functions $e^{-i\zeta x}\bm{\psi}$ and $e^{i\zeta x}\bm{\phi}$ are analytic in the half-space $\{\zeta\in\mathbb{C}|\,\operatorname{Im}{\zeta}>-d\}$. The functions $e^{i\zeta x}\overline{\bm{\psi}}$ and $e^{-i\zeta x}\overline{\bm{\phi}}$ are analytic in the half-space $\{\zeta\in\mathbb{C}|\,\operatorname{Im}{\zeta}<d\}$. In this case, the coefficient $a(\zeta)$ is analytic for $\operatorname{Im}{\zeta}>-d$ while the coefficient $b(\zeta)$ is analytic in the strip defined by $-d<\operatorname{Im}\zeta<d$. More will be said about the analyticity and decay properties of the scattering coefficients in Sec. VI.1.

Furthermore, the symmetry properties

[TABLE]

yield the relations $\overline{a}(\zeta)=a^{*}(\zeta^{*})$ and $\overline{b}(\zeta)=b^{*}(\zeta^{*})$.

II.2 Scattering data and the nonlinear Fourier spectrum

The scattering coefficients introduced in the last section together with certain quantities defined below that facilitate the recovery of the scattering potential are collectively referred to as the scattering data. The nonlinear Fourier spectrum can then be defined as any of the subsets which qualify as the “primordial” scattering data (Ablowitz et al., 1974, App. 5), i.e., the minimal set of quantities sufficient to determine the scattering potential, uniquely.

In general, the nonlinear Fourier spectrum for the potential $q(x,0)$ comprises a discrete and a continuous spectrum. The discrete spectrum consists of the so called eigenvalues $\zeta_{k}\in\mathbb{C}_{+}$, such that $a(\zeta_{k})=0$, and, the norming constants $b_{k}$ such that $\bm{\phi}(x;\zeta_{k})=b_{k}\bm{\psi}(x;\zeta_{k})$. For convenience, let the discrete spectrum be denoted by the set

[TABLE]

For compactly supported potentials, $b_{k}=b(\zeta_{k})$. Note that some authors choose to define the discrete spectrum using the pair $(\zeta_{k},\rho_{k})$ where $\rho_{k}=b_{k}/\dot{a}(\zeta_{k})$ is known as the spectral amplitude corresponding to $\zeta_{k}$ ($\dot{a}$ denotes the derivative of $a$).

The continuous spectrum, also referred to as the reflection coefficient, is defined by $\rho(\xi)={b(\xi)}/{a(\xi)}$ for $\xi\in\mathbb{R}$. The coefficient $a(\zeta)$ and consequently the discrete eigenvalues do not evolve in time. The rest of the scattering data evolves according to the relations $b_{k}(t)=b_{k}e^{-4i\zeta_{k}^{2}t}$ and $\rho(\xi,t)=\rho(\xi)e^{-4i\xi^{2}t}$.

II.3 The Darboux transformation

The Darboux transformation provides a purely algebraic means of adding bound states to a seed solution Neugebauer and Meinel (1984); Lin (1990); Gu et al. (2005). In doing so the $b$-coefficient of the potential remains invariant Lin (1990) while the $a$-coefficient gets modified to reflect the addition of the bound states. In particular, starting from the “vacuum” solution (i.e. the solution for the null-potential), one can compute reflectionless potentials also referred to as the multi-soliton or, more precisely, the $K$-soliton potential with the desired discrete spectrum. The Darboux transformation is carried out by means of Darboux matrices which is described in the following paragraphs.

Let $\mathfrak{S}_{K}$ as defined by (8) be the discrete spectrum to be added to the seed potential. Define the matrix form of the Jost solutions as

[TABLE]

The augmented matrix Jost solution ${v}_{K}(x,t;\zeta)$ can be obtained from the seed solution $v_{0}(x,t;\zeta)$ using the Darboux matrix as

[TABLE]

where $\mu_{K}(\zeta)$ is to be determined. In the following, we summarize the approach proposed by Neugebauer and Meinel Neugebauer and Meinel (1984) which requires the Darboux matrix to be written as

[TABLE]

where the coefficient matrices are such that (for the special case $r=-q^{*}$) $D^{(K)}_{K}=\sigma_{0}$ and

[TABLE]

From the Wronskian relation, we know $a_{0}(\zeta)=\det[v_{0}]$; hence, it follows that

[TABLE]

It is shown in Neugebauer and Meinel (1984) that $\det[D_{K}(x,t;\zeta,\mathfrak{S}_{K})]$ is independent of $(x,t)$. Further, the symmetry imposed by the condition $r=-q^{*}$, requires

[TABLE]

which combined with the fact that Lin (1990)

[TABLE]

yields

[TABLE]

From $\bm{\phi}_{K}(x,t;\zeta_{k})=b_{k}(t)\bm{\psi}_{K}(x,t;\zeta_{k})$, we have

[TABLE]

Note that $\bm{\phi}_{0}(x,t;\zeta_{k})-b_{k}(t)\bm{\psi}_{0}(x,t;\zeta_{k})\neq 0$ on account of $a_{0}(\zeta_{k})\neq 0$, i.e., $\zeta_{k}$ is not an eigenvalue of the seed potential. The $2K$ system of equations in (10) can be used to compute the $2K$ unknown coefficients of the Darboux matrix. Let $U_{K}$ and $U_{0}$ correspond to the augmented potential $q_{K}$ and the seed potential $q_{0}$, respectively; then using the fact that $v_{K}(x,t;\zeta)$ is a Jost solution, we have

[TABLE]

which expands to

[TABLE]

Given that $v_{0}$ is invertible, we must have

[TABLE]

Equating the coefficient of $\zeta^{K}$ to zero, we have

[TABLE]

II.3.1 Darboux matrix of degree one

For the sake of simplicity, let the us consider the seed solution with empty discrete spectra. Let us define the successive discrete spectra $\emptyset=\mathfrak{S}_{0}\subset\mathfrak{S}_{1}\subset\mathfrak{S}_{2}\subset\ldots\subset\mathfrak{S}_{K}$ such that ${\mathfrak{S}}_{j}=\{(\zeta_{j},b_{j})\}\cup{\mathfrak{S}}_{j-1}$ for $j=1,2,\ldots,K$ where $(\zeta_{j},b_{j})$ are distinct elements of $\mathfrak{S}_{K}$.

For single bound state, described by $\mathfrak{S}_{1}$, putting

[TABLE]

the solution of the corresponding linear system (10) yields the Darboux matrix of degree one given by

[TABLE]

The augmented potential then works out as

[TABLE]

The Jost solutions for this new potential can be obtained via the Darboux matrix and the entire procedure can be repeated for adding another bound state to the augmented potential. Suppressing the $x$ and $t$ dependence for the sake of brevity, it follows that the Darboux matrix of degree $K>1$ can be factorized into Darboux matrices of degree one as

[TABLE]

where $D_{1}(\zeta,\mathfrak{S}_{j}|\mathfrak{S}_{j-1}),\,j=1,\ldots,K$ are the successive Darboux matrices of degree one with the convention that $(\zeta_{j},b_{j})=\mathfrak{S}_{j}\cap\mathfrak{S}_{j-1}$ is the bound state being added to the seed solution whose discrete spectra is $\mathfrak{S}_{j-1}$. Using the expression in (12), we have

[TABLE]

where

[TABLE]

for $(\zeta_{j},b_{j})\in\mathfrak{S}_{K}$ and the successive Jost solutions, ${v}_{j}=(\bm{\phi}_{j},\bm{\psi}_{j})$, needed in this ratio are computed as

[TABLE]

The successive potentials are given by

[TABLE]

See Fig. 2 for a schematic representation of the DT.

If the seed Jost solution $v_{0}(x,t;\zeta)$ corresponding to the seed potential $q_{0}(x,t)$ is known, then the Darboux transformations can be readily carried out over any set of grid points $\{x_{n}\}\subset\mathbb{R}$ in order to compute the augmented potential at these grid points. The resulting order of operational complexity, excluding the cost of evaluating the seed potential and the seed Jost solution, works out to be $\mathop{\mathscr{O}}\left(K^{2}N\right)$ where $N$ is the number of samples of the augmented potential. For the special case of $K$-solitons, the seed potential as well as the seed Jost solutions are trivially known; therefore, this method provides us with an algorithm for computing the $K$-soliton potentials with machine precision. In general, closed form solutions are rarely known for arbitrary potentials; nevertheless, this procedure can be carried out with numerically computed Jost solutions in any discrete framework. This scheme will be referred to as the classical Darboux transformation (CDT) in the rest of the article. The error analysis of this method is carried out in Sec. VI.5.

For multi-solitons, the asymptotic form of the potential as $x\rightarrow\infty$ works out to be

[TABLE]

and as $x\rightarrow-\infty$

[TABLE]

where $a_{j}(\zeta)=a(\zeta;\mathfrak{S}_{j})$ are the successive $a$-coefficients. Therefore, $q_{K}(x,t)$ exhibits exponential decay with a decay constant that is given by $d_{K}=\min_{1\leq j\leq K}\operatorname{Im}{\zeta}_{j}$. This observation allows us to conclude that round off errors in the CDT scheme can be minimized if the eigenvalues are “added” in the decreasing order of the magnitude of their imaginary parts Vaibhav and Wahls (2016). Further, the knowledge of the decay constant can be used to choose an optimal computational domain so that the numerical errors due to domain truncation is minimized (see Sec. VII.1.1).

II.3.2 Effective support of multi-soliton potentials

A multi-soliton potential has an unbounded support, therefore, in any practical application it is mandatory to introduce an effective support with desired energy content. Posed conversely, one may also be interested in choosing the discrete spectrum which leads to a prescribed effective support with desired energy content initially or over a finite duration of evolution.

In case of multi-solitons, the energy content of the side lobe which we wish to truncate is trivially available in the CDT scheme and it can be used as a truncation criteria. Let $\chi_{\Omega}$ denote the characteristic function of $\Omega$ and let $[-L,L]$ ($L>0$) be the domain that needs to be determined so that

[TABLE]

Suppressing the dependence on $t$ for the sake of brevity, the asymptotic expansion of $\phi_{K}(-L;\zeta)e^{-i\zeta L}$ with respect to $\zeta$ yields (Ablowitz et al., 1974, Sec. IV.A)

[TABLE]

and that corresponding to $\psi_{K}(L;\zeta)e^{-i\zeta L}$ yields

[TABLE]

These relationships are also known as the nonlinear Parseval’s relationships. Asymptotic estimates when $L\gg 1$ can be easily obtained from the above relations:

[TABLE]

This allows us to obtain an asymptotic formula for the effective support of a $K$-soliton potential. Define $L=L(\epsilon;\mathfrak{S}_{K})>0$ such that

[TABLE]

then

[TABLE]

under the assumption $\epsilon\sum_{j=1}^{K}\eta_{j}\ll\sum_{j=1}^{K}\omega_{j}\eta_{j}$ where

[TABLE]

Finally, let us note that a binary search algorithm (bisection method) can be devised to solve the nonlinear equation (14) for $L=L(\epsilon,\mathfrak{S}_{K})$ where $[0,W]$ can be taken as the bracketing interval for the root555Numerical tests indicates that $[-W,W]$ tends to over estimate the effective support.. The complexity of such an algorithm (for fixed $t$) works out to be $\mathop{\mathscr{O}}\left(mK^{2}\right)$ where $m$ is the number of bisection steps needed.

II.3.3 Scattering coefficients of a truncated multi-soliton

Let $x=0$ be taken as the point of truncation. Then a multi-soliton potential can be seen as comprising a left-sided profile (supported in ${\mathbb{R}}_{-}\cup\{0\}$) and a right-sided profile (supported in $\{0\}\cup{\mathbb{R}}_{+}$). The respective scattering coefficients of each of the truncated potentials turn out to be a rational function of $\zeta$. These observations were already made by several authors Lamb (1980); Rourke and Morris (1992b); Rourke and Saunders (1994); Steudel (2002); Steudel and Kaup (2008) and a number of different methods do exist for inversion of the scattering data which exploit the rational character of the truncated scattering coefficients. Our numerical scheme also exploits this property; therefore, we discuss this case in some detail below.

Let us consider the left-sided profile, denoted by $q^{(-)}(x,t)$. The Jost solution $\bm{\phi}^{(-)}(x,t;\zeta)$ at $x=0$ can be computed using the Darboux transformation as described above. The Jost solution $\bm{\psi}^{(-)}(x,t;\zeta)$ at $x=0$ corresponds to that of a null-potential, i.e., $\bm{\psi}^{(-)}(0,t;\zeta)=(0,1)^{\intercal}$. The scattering coefficients for the left-sided profile, therefore, works out to be

[TABLE]

This corresponds to the first column of the Darboux matrix $D_{K}(0,t;\zeta,\mathfrak{S}_{K})$, therefore, a purely rational function of $\zeta$ analytic in $\overline{\mathbb{C}}_{+}$. Now, let us consider the right-sided profile, denoted by $q^{(+)}(x,t)$. The Jost solution $\bm{\psi}^{(+)}(x,t;\zeta)$ at $x=0$ can be computed using the Darboux transformation as before while the Jost solution $\bm{\phi}^{(+)}(x,t;\zeta)$ at $x=0$ is given by $\bm{\phi}^{(+)}(0,t;\zeta)=(1,0)^{\intercal}$. Therefore, the relevant scattering coefficients for the right-sided profile works out to be

[TABLE]

This corresponds to the second column of the Darboux matrix $D_{K}(0,t;\zeta,\mathfrak{S}_{K})$ and, therefore, a purely rational function of $\zeta$ analytic in $\overline{\mathbb{C}}_{+}$.

Remark II.1** (Conjugation and reflection).**

The inverse scattering problem for the right-sided profile can be transformed to that of a left-sided profile in the following way: putting $y=-x$, we have

[TABLE]

where $\bm{w}(y)=\sigma_{1}\bm{v}(-y;\zeta)$. Denote the Jost solutions of the new system (i.e. with potential $U^{*}(-y)$) by $\bm{\Psi}(y;\zeta)$, $\overline{\bm{\Psi}}(y;\zeta)$ (first kind) and $\bm{\Phi}(y;\zeta)$, $\overline{\bm{\Phi}}(y;\zeta)$ (second kind), then

[TABLE]

Let $A(\zeta)$, $B(\zeta)$, $\overline{A}(\zeta)$ and $\overline{B}(\zeta)$ be the scattering coefficients for the new system, then

[TABLE]

The discrete eigenvalues do not change, however, the norming constants change as $B_{k}=1/b_{k}$. Now, the scattering coefficients for the left-sided profile obtained as result of truncating the new potential from the right at $x=0$ work out to be

[TABLE]

Therefore, an implementation for the case of left-sided profile is sufficient to solve problems of general nature encountered in forward/inverse NFT.

Remark II.2** (Translation).**

Let us note that there is no loss of generality in choosing the point of truncation to be $x=0$ on account of the translational properties of the discrete spectrum. If we wish to choose the point of truncation to be $x=x_{0}$, we can consider the transformation $x=y+x_{0}$. Define the new potential to be $\widetilde{U}(y)=U(y+x_{0})$ so that

[TABLE]

where $\bm{w}(y;\zeta)=\bm{v}(y+x_{0};\zeta)$. Denote the Jost solutions of the new system by $\bm{\Psi}(y;\zeta)$, $\overline{\bm{\Psi}}(y;\zeta)$ (first kind) and $\bm{\Phi}(y;\zeta)$, $\overline{\bm{\Phi}}(y;\zeta)$ (second kind), then

[TABLE]

Let $A(\zeta)$, $B(\zeta)$, $\overline{A}(\zeta)$ and $\overline{B}(\zeta)$ be the scattering coefficients for the new system, then

[TABLE]

The discrete eigenvalues do not change, however, the norming constants change as $B_{k}=b_{k}e^{-2i\zeta_{k}x_{0}}$.

III Discrete Forward and Inverse Scattering

In this section, we discuss certain discretization schemes for the scattering problem in (3) such that they are amenable to FFT-based fast polynomial arithmetic Henrici (1993). This method of obtaining a discrete scattering problem is referred to as the spectral-domain approach666See Bruckstein and Kailath (1987) for alternative approaches.. We begin with the transformation $\tilde{\bm{v}}=e^{i\sigma_{3}\zeta x}\bm{v}$ so that (3) becomes

[TABLE]

or,

[TABLE]

The next step is to apply linear one-step method Gautschi (2012) to (18) in order to setup a recurrence relation initialized by the given initial condition. Let us note that the method of numerical integration just described above is identified as the exponential integrator based on linear one-step methods, in particular, the integrating factor (IF) method Cox and Matthews (2002). One of the advantages of the transformation carried out above in arriving at (18) is that the “vacuum” solution obtained from the discrete problem is exact.

Remark III.1**.**

In the literature, the usage of the terms “forward scattering” and “inverse scattering” is not made precise; for instance, “forward scattering” could refer to computation of the scattering coefficients $a$ and $b$ or the nonlinear Fourier spectrum. In order to avoid any confusion arising in the usage of these terms, we follow the convention that the term “forward scattering” refers to the computation of the Jost solutions while the term “inverse scattering” refers to the process of recovering the samples of the scattering potential from (the polynomial form of) the Jost solutions. Note that in almost all cases, knowledge of the Jost solutions trivially allows one to compute the truncated discrete scattering coefficients and vice versa, therefore, no confusion should arise in what constitutes as input to the inverse scattering process.

III.1 Discretization in the spectral-domain

In order to discuss various discretization schemes, we take an equispaced grid defined by $x_{n}=L_{1}+nh,\,\,n=0,1,\ldots,N,$ with $x_{N}=L_{2}$ where $h$ is the grid spacing. Define $\ell_{-},\ell_{+}\in\mathbb{R}$ such that $h\ell_{-}=-L_{1}$, $h\ell_{+}=L_{2}$. Further, let us define $z=e^{i\zeta h}$ and treat $\zeta$ as a fixed parameter. For the potential functions sampled on the grid, we set $q_{n}=q(x_{n},t)$, $r_{n}=r(x_{n},t)$ where the time-dependence is suppressed. Using the same convention, $U_{n}=U(x_{n},t)$ and $\widetilde{U}_{n}=\widetilde{U}(x_{n},t)$.

III.1.1 Forward Euler method

The forward Euler (FE) method is the simplest of the finite-difference schemes. It can be stated as

[TABLE]

Setting $Q_{n}=hq_{n}$, $R_{n}=hr_{n}$ and $\Theta_{n}=(1-Q_{n}R_{n})$, we have

[TABLE]

or, equivalently,

[TABLE]

Let us note that the transfer matrix can be transformed to a form that resembles that of the implicit Euler method described in the next section: Putting $\bm{w}_{n}=e^{i\sigma_{3}\zeta h}{\bm{v}}_{n}$, we have

[TABLE]

III.1.2 Implicit Euler method

The backward differentiation formula of order one (BDF1) is also known as the implicit Euler method. The discretization of (18) using this method reads as

[TABLE]

Setting $Q_{n}=hq_{n}$, $R_{n}=hr_{n}$ and $\Theta_{n}=(1-Q_{n}R_{n})$, this scheme can be stated as follows:

[TABLE]

or, equivalently,

[TABLE]

III.1.3 Trapezoidal rule

The trapezoidal rule (TR) happens to be one of the most popular methods of integrating ODEs numerically. The discretization of (18) using this method reads as

[TABLE]

Setting $2Q_{n}=hq_{n}$, $2R_{n}=hr_{n}$ and $\Theta_{n}=1-Q_{n}R_{n}$, this scheme can be stated as follows:

[TABLE]

or, equivalently,

[TABLE]

III.2 Jost solutions and scattering coefficients

In order to express the discrete approximation to the Jost solutions, let us define the vector-valued polynomial

[TABLE]

The Jost solutions $\bm{\psi}$ and $\bm{\phi}$, for the forward/implicit Euler method and the trapezoidal rule, can be written in the form

[TABLE]

where $m+n=N$. Note that the expressions above correspond to the boundary conditions $\bm{\psi}_{N}=z^{\ell_{+}}(0,1)^{\intercal}$ and $\bm{\phi}_{0}=z^{\ell_{-}}(1,0)^{\intercal}$ which translate to $\bm{S}_{0}=(0,1)^{\intercal}$ and $\bm{P}_{0}=(1,0)^{\intercal}$, respectively. The other Jost solutions, $\overline{\bm{\psi}}_{n}$ and $\overline{\bm{\phi}}_{n}$, can be written as

[TABLE]

The recurrence relation for the polynomial functions defined in (25) take the form

[TABLE]

where $M_{n+1}(z^{2})$ with its inverse $z^{-2}\widetilde{M}_{n+1}(z^{2})$ is determined by the respective discretization scheme. The discrete approximation to the scattering coefficients is obtained from the scattered field: $\bm{\phi}_{N}=(a_{N}z^{-\ell_{+}},b_{N}z^{\ell_{+}})^{\intercal}$ yields

[TABLE]

and $\bm{\psi}_{0}=(\overline{b}_{N}z^{\ell_{-}},a_{N}z^{-\ell_{-}})^{\intercal}$ yields

[TABLE]

The quantities $a_{N}$, $b_{N}$ and $\overline{b}_{N}$ above are referred to as the discrete scattering coefficients. Note that these coefficients can only be defined for $\operatorname{Re}\zeta\in[-{\pi}/{2h},\,{\pi}/{2h}]$.

Remark III.2**.**

For the sake of brevity, we may occasionally refer to the polynomials $\bm{S}_{m}(z^{2})$ and $\bm{P}_{n}(z^{2})$ (as opposed to $\bm{\psi}_{n}$ and $\bm{\phi}_{n}$) as the (discrete) Jost solutions.

III.2.1 Discrete spectrum

The eigenvalues are computed by forming $a_{N}(z^{2})$ and employing a suitable root-finding algorithm (see Wahls and Poor (2013) and the references therein for more details). It turns out that the computation of the norming constants by evaluating $b_{N}$ is ill-conditioned on account of the vanishingly small contribution from the solitonic components of the potential. Note that addition of bound states leaves $b$-coefficients invariant; therefore, recovery of the norming constant from $b(\zeta)$ cannot be expected to succeed in all cases. In order to remedy this problem, we use the general definition of the norming constants777Similar approach is reported in Hari and Kschischang (2016) and Aref (2016), however, it is not emphasized in these papers that the norming constants are never defined to be a value of $b(\zeta)$ unless it is guaranteed to be analytic in $\mathbb{C}_{+}$. Note that the study of the errors introduced by the numerical discretization also provides significant insight into why the evaluation of $b_{N}(z^{2})$ at complex values of $\zeta$ is ill-conditioned (see Sec. VI.3).: To this end, we proceed by computing the truncated scattering coefficients. Consider the case of potentials truncated from the right, i.e., $q^{(-)}(x)=\theta(x_{1}-x)q(x)$ where $x_{1}$ is the point of truncation and $\theta(x)$ is the Heaviside step function. The new potential now supported in $(-\infty,x_{1}]$ is interpreted as left-sided with respect to $x_{1}$. The scattering coefficient can be stated in terms of the Jost solutions of the original potential as Lamb (1980)

[TABLE]

Similarly, for potentials truncated from the left, we have

[TABLE]

Denoting the corresponding discrete scattering coefficients by $a^{(-)}_{n}$, $b^{(-)}_{n}$, $a^{(+)}_{m}$ and $\overline{b}^{(+)}_{m}$, where $m+n=\ell_{-}+\ell_{+}$, we have

[TABLE]

where $m=N-n$. Here $n$ can be chosen to be $N/2$. Once an admissible root, $z_{k}$, of $a_{N}(z^{2})$ that corresponds to a soliton is determined888Given that $z_{k}=\exp(i\zeta_{k}h)$ and $\operatorname{Im}\zeta_{k}>0$, we must have $|z_{k}|<1$., the corresponding norming constant is obtained via the proportionality of $\bm{\phi}_{n}$ and $\bm{\psi}_{n}$ which translates to

[TABLE]

The truncated potential does not share discrete eigenvalues with the original potential; therefore, $a^{(+)}_{m}(z_{k}^{2})\neq 0$ and $a^{(-)}_{n}(z_{k}^{2})\neq 0$. The computation of the truncated scattering coefficients can be accomplished by direct evaluation of transfer matrices and subsequently forming the cumulative product leading to an operational complexity of $\mathop{\mathscr{O}}\left(N\right)$ for each eigenvalue (see Sec. III.4.1).

It must be noted that our fast algorithm for forward scattering as discussed Sec. III.5.1 is entirely compatible with the approach suggested here. The scattering coefficients are easily obtainable from the truncated scattering coefficients using the Wronskian relations given in Sec. II.1 as

[TABLE]

Every polynomial multiplication involved above can be carried out efficiently using the FFT algorithm (see Sec. III.5.1).

III.3 Inversion of discrete scattering coefficients

In this section, we consider the problem of recovering the discrete samples of the scattering potential from the discrete scattering coefficients known in the polynomial form. This step is referred to as the discrete inverse scattering step. Starting from the recurrence relation (26), we develop a layer-peeling algorithm similar to that reported by Brenne and Skaar Brenne and Skaar (2003). The common aspect of the layer-peeling step for all kinds of discretization schemes is that using nothing but the knowledge of $\bm{P}_{n+1}(z^{2})$, one should be able to retrieve the samples of the potential needed to compute the transfer matrix $\widetilde{M}_{n+1}(z^{2})$ so that the entire step can be repeated with $\bm{P}_{n}(z^{2})$ until all the samples of the potential are recovered (as illustrated in Fig. 3b). In the following, we summarize the main results which facilitate the layer-peeling step corresponding to the each of the discretization schemes introduced so far. A detailed study of the recurrence relation and the proof of the necessary and sufficient conditions for discrete inverse scattering is provided in Sec. V.

III.3.1 Forward Euler method

The recurrence relation for the forward Euler method yields

[TABLE]

The layer-peeling algorithm based on the forward Euler method uses the relation

[TABLE]

where $P^{(n+1)}_{1,0}\neq 0$ on account of (33). As evident from (19), the transfer matrix, $M_{n+1}(z^{2})$, connecting $\bm{P}_{n}(z^{2})$ and $\bm{P}_{n+1}(z^{2})$ is therefore completely determined by $R_{n}$ (with $Q_{n}=-R^{*}_{n}$).

III.3.2 Implicit Euler method

The recurrence relation for the implicit Euler method yields

[TABLE]

The layer-peeling algorithm based on the implicit Euler method uses the relation

[TABLE]

where $P^{(n+1)}_{1,0}\neq 0$ on account of (35). As evident from (22), the transfer matrix, $\widetilde{M}_{n+1}(z^{2})$, connecting $\bm{P}_{n}(z^{2})$ and $\bm{P}_{n+1}(z^{2})$ is therefore completely determined by $R_{n+1}$ (with $Q_{n+1}=-R^{*}_{n+1}$).

III.3.3 Trapezoidal rule

Let us assume $Q_{0}=0$. The recurrence relation for the trapezoidal rule yields

[TABLE]

where the last relationship follows from the assumption $Q_{0}=0$. For sufficiently small $h$, it is reasonable to assume that $1+Q_{n}R_{n}=2-\Theta_{n}>0$ so that $P^{(n)}_{1,0}>0$ (it also implies that $|Q_{n}|=|R_{n}|<1$). The layer-peeling algorithm based on the trapezoidal scheme uses the relations

[TABLE]

where

[TABLE]

Note that $P^{(n+1)}_{1,0}\neq 0$ and ${P^{(n+1)}_{1,0}-Q_{n+1}P^{(n+1)}_{2,0}}\neq 0$. As evident from (23), the transfer matrix, $\widetilde{M}_{n+1}(z^{2})$, connecting $\bm{P}_{n}(z^{2})$ and $\bm{P}_{n+1}(z^{2})$ is completely determined by the samples $R_{n+1}$ and $R_{n}$ (with $Q_{n+1}=-R^{*}_{n+1}$ and $Q_{n}=-R^{*}_{n}$).

III.4 Sequential algorithm

III.4.1 Forward scattering

The computation of the Jost solution for a given value of the spectral parameter, $\zeta\in\mathbb{C}$ is considered here as the forward scattering step. The direct use of the recurrence relations obtained in Sec. VI.2 gives us a sequential algorithm (see the illustration in Fig. 3a). If $\varpi(n),\,n\in\mathbb{Z}_{+}$, denotes the complexity of computing the Jost solution $\bm{P}_{n}(z^{2})$ for a given $\zeta$, then $\varpi(n+1)=4+\varpi(n)$, counting only the multiplications involved. This recurrence relation yields $\varpi(N)=4N$. It must be noted that the sequential algorithms can be useful for computing norming constants as discussed in Sec. III.2.1 if the eigenvalues are known beforehand. If good initial guesses are known for the eigenvalues, search based methods such as Newton’s method of finding the eigenvalues can also benefit from sequential algorithms Wahls and Poor (2013).

The sequential algorithm for computing the polynomial coefficients of $\bm{P}_{N}(z^{2})$ can also be obtained in the same manner where transfer matrices are now treated as polynomial matrices. If $\varpi(n)$ denotes the complexity of computing the polynomial coefficients for the Jost solution $\bm{P}_{n}(z^{2})$, then $\varpi(n+1)=4(n+1)+\varpi(n)$, counting only the multiplications involved. This yields $\varpi(N)=2(N+1)(N+2)=\mathop{\mathscr{O}}\left(N^{2}\right)$ which is extremely prohibitive for large number of samples. This task can be accomplished much more efficiently using a divide-and-conquer strategy together with FFT-based fast polynomial arithmetic as described in Sec. III.5.1.

III.4.2 Inverse scattering

The inverse scattering step here refers to the retrieval of the samples of the scattering potential from the known polynomial form of the discrete scattering coefficients. This can be accomplished by a sequential layer-peeling algorithm as described in Sec. III.3 (see the illustration in Fig. 3b). If $\varpi(n),\,n\in\mathbb{Z}_{+}$, denotes the complexity of inversion of $\bm{P}_{n}(z^{2})$, then $\varpi(n)=4(n+1)+\varpi(n-1)$ counting only the multiplications. This again yields a complexity of $\mathop{\mathscr{O}}\left(N^{2}\right)$ for inverting $\bm{P}_{N}(z^{2})$. This task can also be accomplished much more efficiently using a divide-and-conquer strategy together with FFT-based fast polynomial arithmetic as described in Sec. III.5.2.

III.5 Fast algorithm: A divide-and-conquer strategy

III.5.1 Forward scattering

The scattering algorithm consists in forming cumulative product of, say $N$, transfer matrices. Given that the transfer matrices have polynomial entries (of maximum degree one), one can use FFT-based polynomial multiplication Henrici (1993) to obtain a fast forward scattering algorithm. In this article we restrict ourselves to the case where $N$ is a power of $2$. Most efficient use of the FFT-based multiplication can be made if we use a divide-and-conquer strategy as in Wahls and Poor (2015a) where products are formed pair-wise culminating in the full transfer matrix. The complexity of obtaining the cumulative transfer matrix from $n$ transfer matrices, denoted by $\varpi(n)$, then satisfies the recurrence relation

[TABLE]

where $\nu(n)=n(3\log_{2}2n+2)$ is the complexity of multiplying two polynomials of degree $n-1$ using the FFT algorithm. The number of pairs is given by $l=\log_{2}N$ so that the recurrence relation yields

[TABLE]

which simplifies to

[TABLE]

Therefore, the complexity of the forward scattering algorithm is $\mathop{\mathscr{O}}\left(N\log^{2}N\right)$. Note that $\varpi(1)$ denotes the cost of obtaining each of the transfer matrices.

Evaluation of $\bm{P}_{N}(z^{2})$ at an arbitrary complex point can be done using Horner’s method (Henrici, 1964, Chap. 3) which has the complexity of $\mathop{\mathscr{O}}\left(N\right)$. However, multipoint evaluation at $M\,\,(\geq N)$ Fourier nodes can be carried out with complexity $\mathop{\mathscr{O}}\left(M\log M\right)$ where $M$ is a power of $2$.

III.5.2 Inverse scattering

In this section, we describe how to obtain a fast layer-peeling algorithm by adapting McClary’s approach McClary (1983) for our discrete inverse scattering problem. Consider the grid $(x_{n})_{0\leq n\leq N}$ and let us label the segment $[x_{n},x_{n+1}]$ by $n+1$ for $n<N$. Recall that the inverse of the transfer matrix $M_{n}(z^{2})$ is $z^{-2}\widetilde{M}_{n}(z^{2})$. The cumulative transfer matrix from the $n$-th segment to the $(n-m+1)$-th segment is given by

[TABLE]

Note that in order to determine the transfer matrices for last $l$ segments starting from the $n$-th segment, it is sufficient to have a partial knowledge of the Jost solution, more specifically999We discuss the case where the underlying one-step method is the trapezoidal rule on account of the fact that the corresponding transfer matrix is the most general among the methods considered in this article., $\{\bm{P}_{n}\}_{l+1}$, where $\{\cdot\}_{l}$ denotes truncation after first $l$ coefficients. Let the complimentary polynomial vector be defined as

[TABLE]

and consider the inverse propagation relation in terms of the inverse of the transfer matrices:

[TABLE]

For every $m>0$, the first two coefficients of the polynomial $\bm{P}_{n-m}(z^{2})$ are required in order to determine the transfer matrix for the segment $n-m$; therefore, $2(l+1-m)>0$ ensures that no contribution comes from the complimentary polynomial in computing these first two coefficients. It then follows that the transfer matrices

[TABLE]

can be determined without needing the complimentary polynomial $\{\bm{P}_{n}(z^{2})\}^{c}_{l+1}$. Once the matrices are determined, the Jost solution needed to determine the transfer matrices for $n-l$ segments works out to be

[TABLE]

All polynomial multiplications can be carried out using the FFT-algorithm. The observations made above makes it clear that a divide-and-conquer strategy can be easily devised in order to speed up the layer-peeling algorithm. For the inversion of the discrete scattering coefficients, we start with the associated Jost solution $\bm{P}_{N}(z^{2})$ where $N$ is a power of $2$, we devise a divide-and-conquer strategy that reduces the original problem into two equal size (in terms of number of segments) subproblems101010Note that the analysis in Sec. III.3 reveals that the number of coefficients associated with $\bm{P}_{N}(z^{2})$ is exactly $N$.. The algorithm can be described as follows:

1. i.

Define a binary tree with the number of levels given by $l=\log_{2}N$ (see Fig. 4). Every parent node forks into two child nodes eventually terminating the tree at the leaf nodes. 2. ii.

Associate $N$ segments with the root node which is assumed to be at the level zero. Number of segments associated with every child node is half of that of the parent node. If $\mathcal{S}(k)$ denotes the number of segments associated with nodes at the $k$-th level, then $\mathcal{S}(k)=N2^{-k}$ for $k=0,1,\ldots,l-1$. 3. iii.

Every node in the binary tree is labeled by the index-coordinates $(j,k)$ where $k$ is the level and $j$ being the horizontal position of the node from the left in any particular level, say $k$, so that $0\leq j\leq k$. If the index of the last segment associated with a given node $(j,k)$ is denoted by $N_{jk}$, then $N_{jk}=2^{j}\mathcal{S}(k)$. 4. iv.

All polynomial products to be formed at any node at the $k$-th level requires executing an FFT-algorithm for vectors of length no more than $2\mathcal{S}(k)$. 5. v.

The segments associated with a node dictate the associated cumulative transfer matrix and the Jost solution (with the required number of coefficients) needed in order to determine the entries of constituting transfer matrices. For the node $(j,k)$, the associated cumulative transfer matrix is

[TABLE]

and the associated Jost solution is $\{\bm{P}_{N_{jk}}(z^{2})\}_{n+1}$. 6. vi.

Our algorithm requires exactly two types of operations to be carried out at every node except for the leaf nodes. The first is the computation of the cumulative transfer matrix once the constituting matrices are known at the child nodes. The second is computing the Jost solution needed by any of the child nodes. Both of these operations boil down to polynomial multiplications, therefore, it can be carried out efficiently using the FFT-algorithm. The samples of the potential are determined at the leaf nodes.

Denoting the complexity of multiplying two polynomials of degree $n-1$ (via the FFT-algorithm) by $\nu(n)$, the recurrence relation for the complexity of the fast layer-peeling procedure, denoted by $\varpi(n)$ (where $n=\mathcal{S}(k)$, the number of segments at level $k$), can be stated as

[TABLE]

The first term on the RHS corresponds to the determination of Jost solution for the second child node assuming that the Jost solution is known at the parent node and the cumulative transfer matrix is known at the first child node. The second term corresponds to the determination of the cumulative transfer matrix at the corresponding parent node using the transfer matrices of the child nodes. Observing

[TABLE]

where the last term on RHS is a correction for the root node since the determination of the cumulative transfer matrix at the root level is unnecessary. Using $\nu(n)=n(3\log_{2}2n+2)$, we have

[TABLE]

valid for $N\geq 4$ where $\varpi(2)$ refers to the cost of executing the leaf node. Therefore, the fast layer-peeling algorithm has the complexity of $\mathop{\mathscr{O}}\left(N\log^{2}N\right)$.

III.6 Inversion of scattering coefficients

Let us assume that the scattering coefficients $a(\zeta)$ and $b(\zeta)$ are analytic in $\overline{\mathbb{C}}_{+}$ such that for $\zeta\in\overline{\mathbb{C}}_{+}$ and some $C>0$, we have

[TABLE]

where $\breve{b}(\zeta)=b(\zeta)e^{2i\zeta L_{2}}$. The precise conditions under which such a situation may arise is discussed in theorems VI.2 and VI.3. We further assume that the potential is supported in a domain of the form $(-\infty,L_{2}]$ or $[L_{1},L_{2}]$. In this section, we would like to develop a method to compute the discrete scattering coefficients from the analytic form of the scattering coefficients so that the corresponding inverse problem can be solved numerically using the layer-peeling algorithm discussed in Sec. III.3. It turns out that this task can be efficiently accomplished using the method developed by Lubich Lubich (1988a) which is used in computing the quadrature weights for convolution-type integrals 111111The method based on the trapezoidal rule also appears in control literature where it is known as the Tustin’s method Tustin (1947)..

Introduce the function $\delta(z)$ as in Lubich (1988a) which corresponds to the A-stable one-step methods, namely, BDF1 and TR:

[TABLE]

Putting $z=e^{i\zeta h}$, let us define the coefficients $a_{k}$ and $\breve{b}_{k}$ as

[TABLE]

The coefficients can be obtained using the Cauchy integrals

[TABLE]

which can be easily computed using FFT. Note that the zeroth coefficient can be computed exactly as

[TABLE]

On account of the decay property of the scattering coefficients with respect to $\zeta$, $a_{0}=\mathop{\mathscr{O}}\left(h\right)$ and $\breve{b}_{0}=\mathop{\mathscr{O}}\left(h\right)$.

Let $f_{k}(h)$ denote either $a_{k}$ or $\breve{b}_{k}$ and let $F(z^{2})$ represent the corresponding integrand in (50). Following Lubich (1988b), we obtain the approximation $f_{k}(h;M)$ for $f_{k}(h)$ as

[TABLE]

where $F_{j}=F(\varrho e^{-i\frac{2\pi jk}{M}})$. Choosing $\varrho\leq 1$ ensures that $\operatorname{Im}\zeta\geq 0$. In order to achieve an accuracy of $\mathop{\mathscr{O}}\left(\epsilon\right)$ for computing $f_{k}(h;M)$ for $k=0,1,\ldots,N$ choose $\log\varrho=(1/N)\log\epsilon$ and $M=N\log(1/\epsilon)$. The Lubich’s method, therefore, delivers discrete scattering coefficients with $\mathop{\mathscr{O}}\left(M\log M\right)$ complexity excluding the cost of function evaluations.

Remark III.3**.**

If it is known that the scattering coefficients are also analytic in $\mathbb{C}_{-}$, say, in the strip $\mathbb{S}_{-}(\mu)=\{\zeta\in\mathbb{C}_{-}|\operatorname{Im}{\zeta}\geq-\mu\}$, then Cauchy’s estimate can be used to show that the Lubich coefficients decay exponentially with $k$. Let $\Gamma=\{z\in\mathbb{C}|\,|z|=\varrho,\,\varrho>1\}$ be such that $[i\delta(z)/2h]\in\overline{\mathbb{C}}_{+}\cup\mathbb{S}(\mu)$ for all $z\in\Gamma$. Then, Cauchy’s estimate gives

[TABLE]

where $f(\zeta)$ stands for $a(\zeta)$ or $\breve{b}(\zeta)$ and $f_{k}(h)$ denotes the $k$-th Lubich coefficients.

III.6.1 Relationship with inverse Fourier-Laplace transform

In case of rational scattering coefficients, the Lubich coefficients $a_{k}$ and $\breve{b}_{k}$ can be computed using the inverse Fourier-Laplace transform of the scattering coefficients. For rational functions121212It suffices for our purpose to consider rational functions with simple poles (See III.7)., resolution into partial fractions offers a straightforward means of computing inverse Fourier-Laplace transform. This property can be exploited to lower the cost of computing the discrete scattering coefficients as follows: Define the functions $\alpha(\tau)$ and $\breve{\beta}(\tau)$ as

[TABLE]

Note that for $\tau<0$, the contour can be closed in $\mathbb{C}_{+}$ and the integrals would evaluate to zero, therefore $\alpha(\tau)$ and $\breve{\beta}(\tau)$ are causal. According to (Lubich, 1988a, Theorem 4.1), the coefficients $a_{k}$ and $\breve{b}_{k}$ approximate the quantities $(2h)\alpha(2hk)$ and $(2h)\breve{\beta}(2hk)$ up to $\mathop{\mathscr{O}}\left(h^{p+1}\right)$, respectively, for $k>0$ (note that the zeroth coefficient is given by (51) which merely requires function evaluation). For the trapezoidal rule, this property is proven in Appendix A. It is observed that agreement between true Lubich coefficients and those computed as stated above improves with increasing $k$. Therefore, one should choose $k>N_{\text{th}}$ where $N_{\text{th}}>0$ is a suitably chosen threshold in order to switch to the partial-fraction variant of computing Lubich coefficients.

III.7 Inversion of rational scattering coefficients: Truncated multi-solitons

In order to obtain a fast version of the Darboux transformations (DT) for generating multi-solitons (Problem I.1), we would like to employ the scattering coefficients obtained as a result of truncation of a $K$-soliton potential at $x=0$. As shown in Sec. II.3.3, the scattering coefficients are rational functions of $\zeta$ with no poles in $\overline{\mathbb{C}}_{+}$. Therefore, the Lubich’s method of obtaining discrete scattering coefficients as described in Sec. III.6 is also applicable here. It must be noted that in order to obtain the complete $K$-soliton potential at a given time $t$, the truncation must be done after computing the time-evolved Darboux matrix.

Discrete inverse scattering proceeds by computing the polynomial vector $\bm{P}_{N}(z^{2})$ associated with the discrete scattering coefficients. Without the loss of generality, we assume that the truncation is done at $x=0$ (see Remark II.2). Let the discrete spectrum of the $K$-soliton be $\mathfrak{S}_{K}$ as defined in Sec. II.2. Using the notations introduced in Sec. II.3.1 (we drop the dependence of the Darboux matrices on $\mathfrak{S}_{K}$ for the sake of brevity) and setting $N_{1}=N/2\in\mathbb{Z}_{+}$, for the left-sided profile, we have

[TABLE]

where truncation after $N_{1}$ terms is implied by the notation $\{\cdot\}_{N_{1}}$. This determines $U(x)$ for $x<0$. The right-sided profile can be generated using the transformation described in Remark II.1 so that

[TABLE]

This would determine $U^{*}(-x)$ for $x<0$. Combining the two parts determines the complete multi-soliton potential. Note that the foregoing description also applies to any set of rational functions which qualify as scattering coefficients of a left-sided or a right-sided profile, respectively.

The operational complexity of this algorithm can be computed by taking into account the complexity of DT at $x=0$, which is $\mathop{\mathscr{O}}\left(K^{2}\right)$, and the complexity of computation of Lubich coefficients which is $\mathop{\mathscr{O}}\left(KM\right)+\mathop{\mathscr{O}}\left(M\log M\right)$ where $M$ is the number of nodes used in evaluating the Cauchy integral. Given that $K\ll M$ and $M=\mathop{\mathscr{O}}\left(N\right)$, the overall complexity of generating the multi-soliton including the layer-peeling step works out to be $\mathop{\mathscr{O}}\left(N(K+\log^{2}N)\right)$. The algorithm presented in this section is referred to as the fast Darboux transformation (FDT) algorithm. As pointed out in Sec. II.3.1, the CDT algorithm offers machine precision for computing $K$-soliton potentials with an operational complexity of $\mathop{\mathscr{O}}\left(K^{2}N\right)$. The fundamental difference between the CDT and the FDT algorithm is depicted in Fig. 5 where it is evident that by avoiding DT-iterations at each of the grid points (except at $x=0$) and using the fast LP algorithm, a lower complexity order algorithm can be obtained.

For any rational function, if the poles and residues are known then resolution into partial fractions offers a straightforward means of computing the inverse Fourier-Laplace transform. Let us apply this idea to the problem of generating multi-solitons as discussed in the last paragraph: Poles of the Jost solutions are known to be $\zeta_{k}^{*}$ (where $\zeta_{k}$ are the discrete eigenvalues), therefore, the resolution of the Darboux matrix into partial fractions reads as

[TABLE]

The inversion of $(\zeta-\zeta_{k}^{*})^{-1}$ leads to terms of the form $-ie^{-i\zeta_{k}^{*}\tau}$, therefore, the quantities $e^{-2ih\zeta_{k}^{*}}$ must be computed beforehand. Excluding the cost of computing the $K$ exponentials, the complexity of this algorithm is $\mathop{\mathscr{O}}\left(KN\right)$ where $N$ is the number of samples in the $\tau$-domain. In practice, replacing Lubich coefficients with that obtained by resolution into partial-fractions leads to increase in error and even failure to converge; however, for larger values of the index, the agreement between the two improves allowing us to reduce the overall complexity of computing the discrete coefficients $\bm{P}^{(N_{1})}_{k}$ by switching to the faster algorithm for $k>N_{\text{th}}$ where $N_{\text{th}}>0$ is a suitably chosen threshold131313A recipe to choose $N_{\text{th}}$ based on the number of samples $N$, the size of the computational domain $(L_{2}-L_{1})$ and the eigenvalue with the smallest imaginary part is provided in the Appendix A..

Before we conclude this section, it is worth mentioning that the case treated by Rourke et al. Rourke and Morris (1992b); Rourke and Saunders (1994) of rational reflection coefficient $\rho(\zeta)$ proceeds by reducing the problem to an equivalent problem of generating multi-solitons on a given half-space. Therefore, such cases are amenable to the method discussed in this article.

III.8 General Darboux transformation: Addition of bound sates

In this section, we address Problem I.2 introduced in the beginning of this article. To this end, let us note that the general Darboux transformation consists in adding a given discrete spectrum $\mathfrak{S}_{K}$ (as defined in Sec. II.2) to a given seed potential, $q_{\text{seed}}=q_{0}(x)$, which is assumed to be admissible as a scattering potential in the ZS-problem. The two algorithms developed for this purpose, namely, the classical Darboux transformation (CDT) and the fast Darboux transformation (FDT) meant to carry out the general Darboux transformation are described in the following subsections. For the sake of brevity of presentation, we restrict ourselves to the case $t=0$.

III.8.1 The CDT algorithm

The basic idea behind the CDT algorithm is described in Sec. II.3.1 and also depicted in Fig. 5a. In the discrete framework developed in Sec. VI.2, the seed Jost solutions (which need to be evaluated at the eigenvalues $\zeta_{j}$ to be added) can be computed via the sequential algorithm discussed in Sec. III.4.1. Using the notations introduced in Sec. II.3.1 and Sec. III.2.1, and, introducing $\beta^{(j-1)}_{n}(z_{j})$ as the discrete approximation to $\beta_{j-1}(x_{n},0;\zeta_{j})$, we have

[TABLE]

where $(\zeta_{j},b_{j})\in\mathfrak{S}_{K}$, $m+n=N$ and $z_{j}=e^{2i\zeta_{j}h}$. Noting that $v^{(j)}_{n}=(\bm{P}^{(j)}_{n},\bm{S}^{(j)}_{m})$, the rest of the steps involved are similar to that discussed in Sec. II.3.1.

The operational complexity of computing the seed Jost solutions at $K$ eigenvalues using the sequential algorithm is $\mathop{\mathscr{O}}\left(KN\right)$ so that the overall complexity of the CDT algorithm is $\mathop{\mathscr{O}}\left(K^{2}N\right)$. A final remark that we would like to make with regard to the CDT algorithm is that numerical computation of the Jost solutions for complex values of the spectral parameter $\zeta$ tends to become inaccurate on account of the $\zeta$-dependence of the truncation error coefficient as discussed in Sec. VI.2. It is therefore recommended that $\operatorname{Im}{\zeta_{k}}$ is kept below a certain threshold.

III.8.2 The FDT algorithm

The fundamental idea of the FDT algorithm is the same as that described in Sec. III.7 which is considers the problem of adding bound states, described by $\mathfrak{S}_{K}$, to a null seed potential. The difference merely lies in how we compute the seed Jost solutions required in the DT-iterations at $x=0$ for a general seed potential. Following Sec. II.3.1 and  III.2.1, note that evaluation of the Jost solutions at $\zeta=\zeta_{j}$ amounts to evaluating the approximating polynomial at $z_{j}=e^{2i\zeta_{j}h}$ (setting $x=0$ and $t=0$), so that the recursive step for computing the $\beta$-coefficients reads as

[TABLE]

where we have assumed $\ell_{-},\,\ell_{+}\in\mathbb{Z}$ for simplicity and $(\zeta_{j},b_{j})\in\mathfrak{S}_{K}$. Noting that $v^{(j)}_{\ell_{-}}=(\bm{P}^{(j)}_{\ell_{-}},\bm{S}^{(j)}_{\ell_{+}})$, other steps of the iteration are identical to that described in Sec. II.3.1. Here, our objective is not to follow the conventional Darboux transformation but merely obtain the truncated scattering coefficients (for the left-sided and the right-sided potential) at the origin so that a fast layer-peeling algorithm can be used to compute the samples of the augmented potential.

The operational complexity of this algorithm can be worked out as follows: The cost of computing the Jost solutions (as a polynomial vector) is $\mathop{\mathscr{O}}\left(N\log^{2}N\right)$ and the cost of evaluation of the Jost solutions using Horner’s scheme is $\mathop{\mathscr{O}}\left(N\right)$ for each of the eigenvalues so that the overall complexity of computing the discrete truncated scattering coefficient at $x=0$ is $\mathop{\mathscr{O}}\left(K^{2}\right)+\mathop{\mathscr{O}}\left(KM\right)+\mathop{\mathscr{O}}\left(M\log M\right)+\mathop{\mathscr{O}}\left(N\log^{2}N\right)$ where $M$ is the number of nodes used in evaluating the Cauchy integral, $N$ is the number of samples of the potential and $K$ is the number of eigenvalues to be added. Observing that $K\ll M$, $N$ and $M=\mathop{\mathscr{O}}\left(N\right)$, the overall complexity is effectively $\mathop{\mathscr{O}}\left(N(K+\log^{2}N)\right)$ including the layer-peeling step.

The convergence behavior of the FDT algorithm is studied in the Sec. VI.5 where it is shown that the Darboux matrices can be computed with the same order of accuracy as that of the underlying one-step method used in the computation of the Jost solutions of the seed potential. Further, the global order of convergence matches that of the underlying one-step method for the computation of Lubich coefficients or the layer-peeling algorithm depending on which of the two is lower.

Finally, let us conclude this section by pointing out that if a fast and sufficiently accurate means of inversion of continuous spectrum (i.e., no bound states present) is available then a fast inverse scattering algorithm can be easily obtained for the general cases using the FDT algorithm outlined in this section. The first results in this direction are reported in Vaibhav and Wahls (2017) where the trapezoidal rule is used to develop two algorithms of complexity $\mathop{\mathscr{O}}\left(N(K+\log^{2}N)\right)$ that exhibit a convergence behavior of $\mathop{\mathscr{O}}\left(N^{-2}\right)$.

IV Benchmarking methods

In this section, we discuss two of the conventional methods which are widely used for solving scattering problems. We would like to benchmark our method against these known methods. Unlike the linear one-step methods, here we employ a staggered grid configuration given by $(x_{n+1/2})_{0\leq n<N}$ such that $x_{n+1/2}=x_{n}+h/2$.

IV.1 Magnus integrator

By applying the Magnus method with one-point Gaussian quadrature (see Magnus (1954); Iserles and Nørsett (1999); Hochbruck and Lubich (2003)) to the original ZS-problem in (3), we obtain

[TABLE]

The exponential operator can be computed exactly as

[TABLE]

where $\Gamma=\sqrt{Q_{n+1/2}R_{n+1/2}-\zeta^{2}h^{2}}$ where $Q_{n+1/2}=hq(x_{n}+h/2,t)$ and $R_{n+1/2}=hr(x_{n}+h/2,t)$. We refer to this integrator as “MG1” signifying Magnus integrator with one-point Gauss quadrature. This method is also referred to as the exponential mid-point rule in the literature and it can be shown to be consistent and stable with an order $p=2$. Additionally, it also forms the part of the Lie-group methods Iserles and Nørsett (1999); Hairer et al. (2006) as it retains the SU$(2)$ structure of the Jost solution $v=(\bm{\phi},\overline{\bm{\phi}})$ for $\zeta\in\mathbb{R}$. It must be noted that this method is specially suited for highly oscillatory problems and has been employed by several authors to solve forward scattering problems Boffetta and Osborne (1992); Burtsev et al. (1998). Finally, let us mention that the method of computing the norming constants as described in Sec. III.2.1 can also be adapted to MG1.

IV.2 Split-Magnus method

A further simplification obtained by applying Strang-type splitting Strang (1968) to the exponential operator provides the right discrete framework for the layer-peeling algorithm. This simplification is achieved as follows:

[TABLE]

The order of approximation is determined by applying the Baker-Campbell-Hausdorff (BCH) formula to the exponential operators (Hairer et al., 2006, Chapter 4). Setting $\Gamma=\sqrt{Q_{n+1/2}R_{n+1/2}}$, we have

[TABLE]

Therefore, the discretization scheme works out to be

[TABLE]

where $\Theta_{n+1/2}=(1-Q_{n+1/2}R_{n+1/2})>0$. This form has been used by a number of authors in connection with the conventional layer-peeling algorithm Bruckstein et al. (1985); Feced and Zervas (2000); Brenne and Skaar (2003) as well as for the fast version of the layer-peeling algorithm Wahls and Poor (2015a, b). By employing the transformation $\bm{w}_{n}=e^{i\zeta\sigma_{3}h/2}\bm{v}_{n}$, we obtain

[TABLE]

which maybe viewed as a modification of the implicit Euler scheme. The integration scheme thus obtained is referred to as the split-Magnus (SM) method. The inverse relationship is given by

[TABLE]

The Jost solution can be put in to the form

[TABLE]

where $\bm{S}_{m}(z^{2})$ and $\bm{P}_{n}(z^{2})$ obey the same kind of transfer matrix relation as in (26) with initial condition $\bm{S}_{0}=(0,1)^{\intercal}$ and $\bm{P}_{0}=(1,0)^{\intercal}$. The scattering coefficients work out to be

[TABLE]

The layer-peeling property can be stated as

[TABLE]

with the following additional constraints:

[TABLE]

The norming constants can be computed using any of the following formulas

[TABLE]

Lastly, we note that a staggered grid configuration may prove superior for potentials with jump discontinuity at any grid point because the sampling of the potential at the points of discontinuity is avoided.

Remark IV.1**.**

It must be noted that the CDT and the FDT algorithms are incompatible with the staggered grid configuration, therefore, the SM integrator is ruled out for all DT-related algorithms.

V Discrete inverse scattering: Necessary and sufficient condition

In this section, we study the necessary and sufficient condition for the inversion of the discrete scattering coefficients within the framework of the numerical discretization introduced in Sec. III. Let $\{(\cdot)_{k}\}_{k=1}^{N}$ denote a sequence of quantities such as scalars, vectors or matrices.

Definition V.1**.**

Let $d$ be a non-negative integer. A polynomial $\bm{P}_{n}(z)$ defined as in (24) (with coefficients $\bm{P}^{(n)}_{k}\in\mathbb{C}^{2},\,k=0,1,\ldots,n$) is said to belong to the class $\mathsf{P}(d;\mathbb{C}^{2})$ if $\deg[\bm{P}_{n}(z)]\leq d$ and, for all $z\in\mathbb{T}$, we have

[TABLE]

with ${P}_{1,0}^{(n)}\in\mathbb{R}_{+}$.

For any $\bm{P}_{n}(z)\in\mathsf{P}(d;\mathbb{C}^{2})$, on equating the coefficient of the zeroth degree term on LHS and RHS of (67), we obtain

[TABLE]

Therefore, $|{P}^{(n)}_{1,k}|\leq 1$ and $|{P}^{(n)}_{2,k}|\leq 1$ for $k=0,1,\ldots,n$. Note that the condition ${P}_{1,0}^{(n)}\in\mathbb{R}_{+}$ ensures that there are no constant phase factors in $\bm{P}_{n}(z)$ because the relation (67) is insensitive to constant phase factors.

Definition V.2** (Para-conjugate).**

For any scalar valued complex function, ${f}(z)$, we define $\overline{f}(z)=f^{*}(1/z^{*})$. For any vector valued complex function, $\bm{f}(z)=(f_{1}(z),f_{2}(z))^{\intercal}$, we define

[TABLE]

For a matrix valued function, $M(z)$, we define

[TABLE]

so that the operation $\overline{(\cdot)}$ is distributive over matrix-vector and matrix-matrix products.

Based on the discrete formulation of the ZS-problem in Sec. III, we identified a discrete representation of the Jost solution which can also be stated in the form (leaving out the factors independent of $n$)

[TABLE]

such that the column vectors are linearly independent for all $z\in\mathbb{C}$. This implies, $\det[w_{n}]\neq 0$. In fact, the determinant must turn out to be independent of $z^{2}$ so that we may put $\det[w_{n}]=W_{n}$ which translates into the constraint141414Given that $\det[w_{n}]$ is a polynomial, the only way $\det[w_{n}]\neq 0$ is when it is a polynomial of degree zero.

[TABLE]

For $z\in\mathbb{T}$,

[TABLE]

This condition is necessary for $w_{n}$, defined by (68), to be a Jost solution. Further, it is easy to verify that $w_{n}$ satisfies the relation

[TABLE]

Finally, let us note that $\tilde{w}_{n}=w_{n}/{\sqrt{-W_{n}}}$ forms a $\text{SU}(2)$-valued sequence for $z\in\mathbb{T}$.

The discrete scattering problem will be assumed to be stated in the form of a recurrence relation which reads as

[TABLE]

where $M_{n}(z^{2})$ is a polynomial matrix of degree one. Note that $w_{n}$ as defined by (68) satisfies the relation $\overline{w}_{n}=-w_{n}$; therefore, in order that $w_{n+1}$ be a Jost solution, we must have $z^{-1}M_{n+1}(z^{2})=z\overline{M}_{n+1}(z^{2})$. This relationship expands to

[TABLE]

and $\det[M_{n}(z^{2})]=z^{2}C_{n}$ where $C_{n}$ is independent of $z$. Introducing the functions

[TABLE]

it follows that the general form of the transfer matrix (of degree one in $z^{2}$) can be written as

[TABLE]

with

[TABLE]

Let the inverse of $M_{n}(z^{2})$ be denoted by $z^{-2}\widetilde{M}_{n}(z^{2})$ which also satisfies a similar symmetry relation as in (72) and

[TABLE]

so that $\overline{\widetilde{M}}_{n}=z^{2}\widetilde{M}_{n}$. Further, it is straightforward to verify that, for $z\in\mathbb{T}$, the matrices $z^{-1}M_{n}/\sqrt{C_{n}}$ and $z^{-1}\widetilde{M}_{n}/\sqrt{C_{n}}$ are elements of $\text{SU}(2)$. The discrete scattering problem in its unitary form reads as

[TABLE]

Introducing $\mu_{n},A_{n},B_{n}\in\mathbb{C}$, the independent elements of transfer matrix can be put into the form

[TABLE]

so that

[TABLE]

Setting $\mu_{n}=|\mu_{n}|e^{i\theta_{n}}$, the transfer matrix admits of the following factorization

[TABLE]

For the cases considered in this article, $\theta_{n}=0$, therefore, we assume $\mu_{n}\in\mathbb{R}_{+}$ so that it does not play a role in the unitary form of the transfer matrix for $z\in\mathbb{T}$. For a given initial condition and fixed sequence of transfer matrices, the recurrence relation (71) leads to a unique polynomial associated with the Jost solution $w_{n}$. In particular, the following result is straightforward:

Lemma V.1**.**

Let $N$ be a finite positive integer. Let the vectors $\bm{A}=(A_{1},A_{2},\ldots,A_{N}),\,\bm{B}=(B_{1},B_{2},\ldots,B_{N})\in\mathbb{C}^{N}$ and $\bm{\mu}=(\mu_{1},\mu_{2},\ldots,\mu_{N})\in\mathbb{R}^{N}_{+}$ define $\{M_{n}(z^{2})\}_{n=1}^{N}$ through (79). Let $w_{0}=\sigma_{1}$, then the recurrence relation (71) determines a sequence of Jost solutions $\{w_{n}\}_{n=1}^{N}$ such that for every $n$ ($1\leq n\leq N$) there exists a unique polynomial $\bm{P}_{n}(z^{2})$ associated with $w_{n}$.

Now let us consider an arbitrary polynomial $\bm{P}_{n}(z^{2})$ satisfying (69) for $n\geq 0$. Assume $\bm{P}_{n}(z^{2})/\sqrt{-W_{n}}\in\mathsf{P}(n;\mathbb{C}^{2})$ and let $\bm{P}_{n+1}(z^{2})$ be associated with $w_{n+1}$. To understand the properties of the polynomial $\bm{P}_{n+1}(z^{2})$, we consider the recurrence relation (71). Equating the coefficients of the zeroth degree term on the RHS and the LHS of (71), we have

[TABLE]

It is straightforward to see that

[TABLE]

and

[TABLE]

Therefore, in order that $\bm{P}_{n+1}(z^{2})/\sqrt{-W_{n+1}}\in\mathsf{P}(n+1;\mathbb{C}^{2})$ where $W_{n+1}=\det[w_{n+1}]=C_{n+1}W_{n}$, we must have

[TABLE]

Lemma V.2**.**

Under the assumption of the previous lemma, setting $W_{n}=\det[w_{n}]$, the polynomial $\bm{P}_{n}(z^{2})/\sqrt{-W_{n}}\in\mathsf{P}(n;\mathbb{C}^{2})$ if and only if the sequence $\{(A_{n},B_{n})\}_{n=1}^{N}$ satisfies the constraint $(1-A_{n}B_{n+1}^{*})\in\mathbb{R}_{+}$ for $1\leq n<N$. If $B_{1}=0$, then $\bm{P}_{n}(z^{2})/\sqrt{-W_{n}}\in\mathsf{P}(n-1;\mathbb{C}^{2})$.

Proof.

Using the recurrence relation (80) and the property (81) for all $n\geq 0$, it is straightforward to see that

[TABLE]

for $n>0$ while $P^{(1)}_{1,0}=\mu_{1}$. The proof of the first part of the lemma follows from this relation.

For the second part, equating the coefficients of $(z^{2})^{n+1}$ on the RHS and the LHS of (71), we have

[TABLE]

for $n\geq 0$. These relations yield

[TABLE]

for $n>0$ while $P^{(1)}_{1,1}=-\mu_{1}A_{1}^{*}B_{1}$ and $P^{(1)}_{2,1}=\mu_{1}B_{1}$. Therefore, if $B_{1}=0$ then $\bm{P}^{(n)}_{n}=0$ for $1\leq n\leq N$. ∎

Next we would like to analyze the inverse problem described as follows: Given an arbitrary polynomial $\bm{P}_{n+1}(z^{2})$ associated with $w_{n+1}$ satisfying

[TABLE]

for $n\geq 0$ such that $\bm{P}_{n+1}(z^{2})/\sqrt{-W_{n+1}}\in\mathsf{P}(n+1;\mathbb{C}^{2})$. Find a polynomial $\bm{P}_{n}(z^{2})$ associated with $w_{n}$ and a transfer matrix $M_{n+1}(z^{2})$ of the form (79) such that $w_{n}$, defined by

[TABLE]

is a Jost solution. If such a polynomial $\bm{P}_{n}(z^{2})$ exists then it must be consistent with the recurrence relation

[TABLE]

or, equivalently,

[TABLE]

Equating the coefficient of $z^{-2}$ on the RHS of (86), we have

[TABLE]

which yields the recurrence relation (81). Equating the coefficients of $z^{0}$ on the RHS and the LHS of (86), we obtain

[TABLE]

This yields

[TABLE]

and

[TABLE]

which thanks to (78) ($\mu_{n+1}\in\mathbb{R}_{+}$) becomes identical to (80). Note that the relationship (89) can also be verified by equating the coefficients of $z^{2}$ on the RHS and the LHS of (87):

[TABLE]

This yields

[TABLE]

which is identical to (89) thanks to (78). Now, using (89) and (90), we have

[TABLE]

where

[TABLE]

So far we have found that the parameter $A_{n+1}$ of $M_{n+1}(z^{2})$ must be set according to (81) so that we may write

[TABLE]

where $\chi_{n+1}$ is known but $B_{n+1}$ is still an unknown. In order to compute $B_{n+1}$, we introduce a free parameter, $\lambda_{n}=P^{(n)}_{2,0}/P^{(n)}_{1,0}$, so that from (91), we have

[TABLE]

Now, let us observe that

[TABLE]

and

[TABLE]

so that

[TABLE]

Now, the zeroth degree coefficient given by (93) simplifies to

[TABLE]

Therefore, in order that $\bm{P}_{n}(z^{2})/\sqrt{-W_{n}}\in\mathsf{P}(n;\mathbb{C}^{2})$ where $W_{n}=W_{n+1}/C_{n+1}$, we must have

[TABLE]

The above condition can be enforced by setting $\lambda_{n}=\chi_{n+1}\omega_{n}$ where we restrict ourselves to the case $\omega_{n}\in\mathbb{R},\,\omega_{n}\geq 0$. Under this condition, the expressions for $B_{n+1}$ and $\bm{P}_{0}^{(n)}$ simplifies to

[TABLE]

and

[TABLE]

respectively. Clearly, the transfer matrix $M_{n+1}(z^{2})$ as well as the polynomial $\bm{P}_{n}(z^{2})$ is not unique as it depends on a free parameter $\omega_{n}\geq 0$. Note that the parameter $\mu_{n+1}$ turns out to be merely a scale factor which does not play a role in the unitary form of the discrete scattering problem. Finally, let us observe that in order to predict the highest degree term that is non-zero in $\bm{P}_{n}(z^{2})$, the recurrence relation for $(z^{2})^{n}\overline{\bm{P}}(z^{2})$ can be considered where the zeroth degree term is $i\sigma_{2}\bm{P}^{(n)*}_{n}$ so that

[TABLE]

Remark V.1**.**

In the discrete inverse scattering case, the two formulas (81) and (92) remain invariant under any scaling of the polynomial $\bm{P}_{n+1}(z^{2})$. Therefore, knowledge of either $\bm{P}_{n+1}(z^{2})$ or $\bm{P}_{n+1}(z^{2})/\sqrt{-W_{n+1}}$ is sufficient to determine the transfer matrix $M_{n+1}(z^{2})$.

The discussion above regarding the discrete inverse scattering step can be summarized in the following lemma:

Lemma V.3**.**

Given $\bm{P}_{n+1}(z^{2})/\sqrt{-W_{n+1}}\in\mathsf{P}(d;\mathbb{C}^{2})$ where $d\in\{n+1,n\}$ and $\omega_{n}\in\mathbb{R}_{+}$, there exists a unique unitary matrix $M_{n+1}(z^{2})/\sqrt{C_{n+1}}$ for $z\in\mathbb{T}$ and a polynomial $\bm{P}_{n}(z^{2})/\sqrt{-W_{n}}\in\mathsf{P}(d-1;\mathbb{C}^{2})$ such that

[TABLE]

Further, if $\omega_{n}\leq 1$, then

[TABLE]

Proof.

The first part of the lemma is evident from the discussion above. The second part follows from the inequality

[TABLE]

for $\omega_{n}\leq 1$. ∎

Next, we consider some of special cases where it is possible to obtain unique solution of the discrete inverse scattering problem. It is worth noting that these special cases belong to a certain choice of the values $\{\omega_{n}\}_{n\in\mathbb{Z}}$.

V.1 Case I: $A_{n}=B_{n+1}$

Let $A_{n}=B_{n+1}$ and assume $A_{n}\in\mathbb{D}$. Then the forward scattering problem described in Lemma V.1 always yields a polynomial $\bm{P}_{n}(z^{2})/\sqrt{-W_{n}}\in\mathsf{P}(n;\mathbb{C}^{2})$ on account of Lemma V.2.

For discrete inverse scattering, the condition $A_{n}=B_{n+1}$ amounts to $B_{n+1}=\chi_{n+1}\omega_{n}$. From (97), we have

[TABLE]

which yields

[TABLE]

as the admissible solution (the other root violates the positivity constraint in (96)). For this case, the expression (98) for the zeroth degree coefficient simplifies to

[TABLE]

In the Lemma V.3, we favor the case of $d=n$ so that number of (vector) coefficients associated with $\bm{P}_{n}(z^{2})$ be $n$. If the steps described in the aforementioned lemma are carried out recursively to the point $n=0$, we obtain

[TABLE]

Note that $\chi_{1}=0$ on account of $\bm{P}^{(1)}_{1}=0$; therefore,

[TABLE]

Finally, we state the main result of this section which is a now merely a consequence of the preceding lemmas applied to the case at hand:

Proposition V.4**.**

Let $\bm{A}=(A_{1},A_{2},\ldots,A_{N})\in\mathbb{D}^{N}$ be an arbitrary vector. Let the transfer matrices $\{M_{n}(z^{2})\}_{n=1}^{N}$ be determined by (79) using $\bm{A}$ together with $\bm{B}\in\mathbb{D}^{N}$ given by $B_{1}=0$ and $B_{n}=A_{n-1}$ for $1<n\leq N$. Then, corresponding to the initial condition $\bm{P}_{0}(z^{2})=(1,0)^{\intercal}$, the recurrence relation

[TABLE]

yields a unique polynomial $\bm{P}_{N}(z^{2})/\sqrt{-W_{N}}\in\mathsf{P}(N-1;\mathbb{C}^{2})$ with $(-W_{N})=\prod_{n=1}^{N}C_{n}>0$ such that

[TABLE]

Conversely, for any given polynomial $\breve{\bm{P}}_{N}(z^{2})\in\mathsf{P}(N-1;\mathbb{C}^{2})$ such that

[TABLE]

there exists a unique vector $\bm{A}=(A_{1},A_{2},\ldots,A_{N})\in\mathbb{D}^{N}$ which determines the the transfer matrices $\{\widetilde{M}_{n}(z^{2})/\sqrt{C_{n}}\}_{n=1}^{N}$ as stated above such that the recurrence relation

[TABLE]

starting from $n=N$ yields $\breve{\bm{P}}_{0}(z^{2})=(1,0)^{\intercal}$.

Putting $\breve{\bm{P}}_{N}(z^{2})=\bm{P}_{N}(z^{2})/\sqrt{-W_{N}}$, note that the condition

[TABLE]

corresponds to the fact that $A_{N}\in\mathbb{D}$ in the direct part of the last proposition. The condition above is imposed explicitly in the converse part in order to ensure $A_{N}\in\mathbb{D}$.

Corollary V.5**.**

Let $\bm{A}=(A_{1},A_{2},\ldots,A_{N})\in\mathbb{D}^{N}$ correspond to $\breve{\bm{P}}_{N}(z^{2})\in\mathsf{P}(N-1;\mathbb{C}^{2})$ as in the converse part of the last proposition. Then the following estimate holds:

[TABLE]

Proof.

The proof follows from the relation (100) of Lemma V.3. ∎

We conclude this section with the discussion of the trapezoidal rule which corresponds to the case at hand. Let $\bm{Q}=(Q_{1},Q_{2},\ldots,Q_{N})\in\mathbb{D}^{N}$. In the case of trapezoidal rule, it follows from the description in Sec. III.1.3 that the coefficients $A_{n}$ and $B_{n}$ introduced in (77) satisfy

[TABLE]

with $A_{N}=Q_{N}$ and we choose $Q_{0}=B_{1}=0$. It also follows that the quantities $\mu_{n}\in\mathbb{R}_{+}$ introduced in (77) are given by

[TABLE]

Further, we have

[TABLE]

while $C_{1}=\Theta^{-1}_{1}$.

V.2 Case II: $A_{n}\neq B_{n+1}$

First, let us assume that $B_{n}=0$. The discussion of the forward scattering problem is identical to that of the previous case. For discrete inverse scattering, this case corresponds to $\omega_{n}=1$. The expression for the zeroth degree coefficient (98) simplifies to

[TABLE]

As in the last section, we favor the case of $d=n$ in the Lemma V.3. Again, if the steps described in the aforementioned lemma are carried out recursively to the point $n=0$, it is easy to conclude that

[TABLE]

The necessary and sufficient condition for discrete inverse scattering in this case can be stated as:

Proposition V.6**.**

Let $\bm{A}=(A_{1},A_{2},\ldots,A_{N})\in\mathbb{C}^{N}$ be an arbitrary vector. Let the transfer matrices $\{M_{n}(z^{2})\}_{n=1}^{N}$ be determined by (79) using $\bm{A}$ together with $B_{n}=0$ for $1\leq n\leq N$. Then, corresponding to the initial condition $\bm{P}_{0}(z^{2})=(1,0)^{\intercal}$, the recurrence relation

[TABLE]

yields a unique polynomial $\bm{P}_{N}(z^{2})/\sqrt{-W_{N}}\in\mathsf{P}(N-1;\mathbb{C}^{2})$ with $(-W_{N})=\prod_{n=1}^{N}C_{n}>0$.

Conversely, for any given polynomial $\breve{\bm{P}}_{N}(z^{2})\in\mathsf{P}(N-1;\mathbb{C}^{2})$ there exists a unique vector $\bm{A}=(A_{1},A_{2},\ldots,A_{N})\in\mathbb{C}^{N}$ which determines the the transfer matrices $\{\widetilde{M}_{n}(z^{2})/\sqrt{C_{n}}\}_{n=1}^{N}$ as stated above such that the recurrence relation

[TABLE]

starting from $n=N$ yields $\breve{\bm{P}}_{0}(z^{2})=(1,0)^{\intercal}$.

Corollary V.7**.**

Let $\bm{A}=(A_{1},A_{2},\ldots,A_{N})\in\mathbb{C}^{N}$ correspond to $\breve{\bm{P}}_{N}(z^{2})\in\mathsf{P}(N-1;\mathbb{C}^{2})$ as in the converse part of the last proposition. Then the following estimate holds:

[TABLE]

Secondly, let us assume that $A_{n}=0$. The discussion of the forward scattering problem is identical to that of the previous case. For discrete inverse scattering, this case corresponds to $\omega_{n}=0$. The expression for the zeroth degree coefficient (98) simplifies to

[TABLE]

Here, we favor the case of $d=n+1$ in the Lemma V.3. Again, if the steps described in the aforementioned lemma are carried out recursively to the point $n=0$, it is easy to conclude that

[TABLE]

The expression for the highest degree coefficient (99), simplifies to

[TABLE]

The necessary and sufficient condition for discrete inverse scattering in this case can be stated as:

Proposition V.8**.**

Let $\bm{B}=(B_{1},B_{2},\ldots,B_{N})\in\mathbb{C}^{N}$ be an arbitrary vector. Let the transfer matrices $\{M_{n}(z^{2})\}_{n=1}^{N}$ be determined by (79) using $\bm{B}$ together with $A_{n}=0$ for $1\leq n\leq N$. Then, corresponding to the initial condition $\bm{P}_{0}(z^{2})=(1,0)^{\intercal}$, the recurrence relation

[TABLE]

yields a unique polynomial $\bm{P}_{N}(z^{2})/\sqrt{-W_{N}}\in\mathsf{P}(N;\mathbb{C}^{2})$ with $(-W_{N})=\prod_{n=1}^{N}C_{n}>0$.

Conversely, for any given polynomial $\breve{\bm{P}}_{N}(z^{2})\in\mathsf{P}(N;\mathbb{C}^{2})$ there exists a unique vector $\bm{B}=(B_{1},B_{2},\ldots,B_{N})\in\mathbb{C}^{N}$ which determines the the transfer matrices $\{\widetilde{M}_{n}(z^{2})/\sqrt{C_{n}}\}_{n=1}^{N}$ as stated above such that the recurrence relation

[TABLE]

starting from $n=N$ yields $\breve{\bm{P}}_{0}(z^{2})=(1,0)^{\intercal}$.

Corollary V.9**.**

Let $\bm{B}=(B_{1},B_{2},\ldots,B_{N})\in\mathbb{C}^{N}$ correspond to $\breve{\bm{P}}_{N}(z^{2})\in\mathsf{P}(N;\mathbb{C}^{2})$ as in the converse part of the last proposition. Then the following estimate holds:

[TABLE]

V.2.1 Implicit Euler method

Let $\bm{Q}=(Q_{1},Q_{2},\ldots,Q_{N})\in\mathbb{C}^{N}$. For the implicit Euler method, it is evident from the discussion in Sec. III.1.2 that

[TABLE]

and

[TABLE]

Further,

[TABLE]

V.2.2 Split-Magnus method

For the split-Magnus method, we consider the samples on a staggered grid so that $\bm{Q}=(Q_{1/2},Q_{3/2},\ldots,Q_{N-1/2})\in\mathbb{C}^{N}$. It is evident from the discussion in Sec. IV.2 that

[TABLE]

and

[TABLE]

Further, $C_{n}=1$.

V.2.3 Forward Euler method

Let $\bm{Q}=(Q_{0},Q_{1},\ldots,Q_{N-1})\in\mathbb{C}^{N}$. For the forward Euler method, it is evident from the discussion in Sec. III.1.1 that

[TABLE]

and $\mu_{n}=1$. Further,

[TABLE]

VI Stability and convergence analysis

The main objective of this section is to carry out an error-analysis for various steps involved in the algorithms proposed in Sec. III. We first study the analyticity properties of the scattering coefficients in order to understand the difficulties involved in transitioning from the continuous to the discrete regime. In Sec. VI.2, we study the stability and convergence of the numerical scheme for forward scattering. Note that the convergence of the layer-peeling algorithm where the input is synthesized using Lubich’s method is not discussed in this work, instead, we study it empirically. The error propagation in layer-peeling procedure has been addressed in the work of Bruckstien et al. Bruckstein et al. (1986); however, on account of the underlying assumption of piece-wise constant potential, the question of convergence beyond the first order cannot be addressed in their work. We leave these aspects for future research.

Notations

The class of $m$-times differentiable complex-valued functions is denoted by $\mathsf{C}^{m}$. A function of class $\mathsf{C}^{m}$ is said to belong to $\mathsf{C}_{0}^{m}(\Omega)$, if the function and its derivatives up to order $m$ have a compact support in $\Omega$ and if they vanish on the boundary ($\partial\Omega$). Complex-valued functions of bounded variation over over $\mathbb{R}$ is denoted by $\mathsf{BV}$ and the variation of any function $f\in\mathsf{BV}$ over $\Omega\in\mathbb{R}$ is denoted by $\mathscr{V}[f;\Omega]$. If $q\in\mathsf{BV}$, then $\partial_{x}q\in\mathsf{L}^{1}$ exists almost everywhere such that $\|\partial_{x}q\|_{\mathsf{L}^{1}}\leq\mathscr{V}[q;\Omega]$ (Jones, 2001, Chap. 16). Let $q^{(1)}$ to be equivalent to $\partial_{x}q$ so that $\|q^{(1)}\|_{\mathsf{L}^{1}}=\|\partial_{x}q\|_{\mathsf{L}^{1}}$.

Let $J=(-\infty,L]$ and $d>0$. A complex-valued function $f(x)$ is said to belong to the class $\mathsf{E}_{d}(J)$ if $\operatorname*{supp}f\subset J$ and there exists a constant $\kappa_{\infty}>0$ such that the estimate $|f(x)|{\leq}\kappa_{\infty}e^{-2d|x|}$ holds almost everywhere in $J$. Clearly, $\mathsf{E}_{d}(J)\subset\mathsf{L}^{p}(J)$ for $1\leq p\leq\infty$. Define $\mathbb{S}_{-}(\mu)=\{\zeta\in\mathbb{C}_{-}|\operatorname{Im}{\zeta}\geq-\mu\}$.

VI.1 Compactly supported and one-sided potentials

The Jost solution for compactly supported and one-sided potential are known to have analytic continuation into the upper-half of the complex plane Ablowitz et al. (1974); Ablowitz and Segur (1981). We detail some of these analyticity and decay properties of the Jost solutions required for our purpose. This discussion is motivated by the fact that our fast Darboux transformation (FDT) algorithm discussed in Sec. III.8.2 proceeds by computing the Jost solutions of a truncated potential which can be interpreted as one-sided (if it does not have a compact support). Further, the analyticity properties of the Jost solutions also determine the behavior of the Lubich coefficients as discussed in Sec. III.6.

We begin with a study of the modified Jost solutions defined by

[TABLE]

Let $\Omega=[L_{1},L_{2}]$ in the following unless stated otherwise. The system of equations (18) can be transformed into a set of Volterra integral equations of the second kind for the modified Jost solution $\widetilde{\bm{P}}(x;\zeta)$:

[TABLE]

where $\bm{\Phi}(x;\zeta)=(\Phi_{1},\Phi_{2})^{\intercal}\in\mathbb{C}^{2}$ with

[TABLE]

and the Volterra kernel $\mathcal{K}(x,y;\zeta)=\operatorname{diag}(\mathcal{K}_{1},\mathcal{K}_{2})\in\mathbb{C}^{2\times 2}$ is such that

[TABLE]

with $\mathcal{K}(x,y;\zeta)=0$ for $y>x$.

Theorem VI.1**.**

Let $q\in\mathsf{L}^{1}$ be supported in $\Omega=[L_{1},L_{2}]$ with $\kappa=\|q\|_{\mathsf{L}^{1}}$. Then the estimate

[TABLE]

holds with $C=\|\bm{D}\|\cosh\kappa$ where $\bm{D}=(\kappa^{2}/2,\kappa)^{\intercal}$.

Proof.

The proof can be obtained using the same method as in Ablowitz et al. (1974); Novikov et al. (1984). For fixed $\zeta\in\overline{\mathbb{C}}_{+}$, let $\mathscr{K}$ denote the Volterra integral operator in (110) corresponding to the kernel $\mathcal{K}(x,y;\zeta)$ such that

[TABLE]

Consider the $\mathsf{L}^{\infty}(\Omega)$-norm (Gripenberg et al., 1990, Chap. 9) of $\mathscr{K}$ given by

[TABLE]

so that $\|\mathscr{K}\|_{\mathsf{L}^{\infty}(\Omega)}\leq\kappa^{2}/2$. The resolvent $\mathscr{R}$ of this operator exists and is given by the Neumann series $\mathscr{R}=\sum_{n=1}^{\infty}\mathscr{K}_{n}$ where $\mathscr{K}_{n}=\mathscr{K}\circ\mathscr{K}_{n-1}$ with $\mathscr{K}_{1}=\mathscr{K}$. It can also be shown using the methods in Ablowitz et al. (1974); Novikov et al. (1984) that $\|\mathscr{K}_{n}\|_{\mathsf{L}^{\infty}(\Omega)}\leq{\kappa^{2n}}/{(2n)!}$, yielding the estimate $\|\mathscr{R}\|_{\mathsf{L}^{\infty}(\Omega)}\leq[\cosh(\kappa)-1]$. Therefore, for any $\bm{\Phi}(x;\zeta)\in\mathsf{L}^{\infty}(\Omega)$, the relationship $\widetilde{\bm{P}}(x;\zeta)=\bm{\Phi}(x;\zeta)+\mathscr{R}[\bm{\Phi}](x;\zeta)$ implies, for $\zeta\in\overline{\mathbb{C}}_{+}$,

[TABLE]

The result for $\overline{\mathbb{C}}_{+}$ in (113) follows from the observation that, for $\zeta\in\overline{\mathbb{C}}_{+}$, $\|\bm{\Phi}(x;\zeta)\|_{\mathsf{L}^{\infty}(\Omega)}\leq\|\bm{D}\|$ where $\bm{D}=({\kappa^{2}}/{2},\kappa)^{\intercal}$. Therefore, $C$ can be chosen to be $\|\bm{D}\|\cosh\kappa$. For the case $\mathbb{C}_{-}$ of (113), we consider $\widetilde{\bm{P}}_{-}(x;\zeta)=\widetilde{\bm{P}}(x;\zeta)e^{-2i\zeta x}$. The Volterra integral equations then reads as $\widetilde{\bm{P}}_{-}(x;\zeta)$:

[TABLE]

where $\bm{\Phi}_{-}(x;\zeta)=\bm{\Phi}(x;\zeta)e^{-2i\zeta x}\in\mathbb{C}^{2}$ and the Volterra kernel $\mathcal{K}_{-}(x,y;\zeta)=\operatorname{diag}(\mathcal{K}^{(-)}_{1},\mathcal{K}^{(-)}_{2})\in\mathbb{C}^{2\times 2}$ is such that

[TABLE]

with $\mathcal{K}_{-}(x,y;\zeta)=0$ for $y>x$. Using the approach outlined above, it is possible to show that, for $\zeta\in\mathbb{C}_{-}$, $\|\widetilde{\bm{P}}_{-}(x;\zeta)\|_{\mathsf{L}^{\infty}(\Omega)}\leq\cosh(\kappa)\|\bm{\Phi}_{-}(x;\zeta)\|_{\mathsf{L}^{\infty}(\Omega)}$. The result for the case $\zeta\in\mathbb{C}_{-}$ in (113) then follows from the observation that $\|\bm{\Phi}_{-}(x;\zeta)\|_{\mathsf{L}^{\infty}(\Omega)}\leq\|\bm{D}\|e^{2\operatorname{Im}(\zeta)L_{1}}$ for $\zeta\in\mathbb{C}_{-}$. ∎

Theorem VI.2**.**

Let $q\in\mathsf{BV}$ with support in $\Omega=[L_{1},L_{2}]$ such that $q(x)=0$ for $x\in\partial\Omega$. Then, there exists a constant $C>0$ independent of $\zeta\in\mathbb{C}$ such that the estimate

[TABLE]

holds.

Proof.

Consider the first term on the RHS of (110): Integrating by parts, we obtain

[TABLE]

so that

[TABLE]

Setting $2D_{1}=\|q\|^{2}_{2}+\|q\|^{2}_{1}+\|q\|_{1}\|q^{(1)}\|_{1}$, we have

[TABLE]

Again, integrating by parts, we have

[TABLE]

so that

[TABLE]

Putting $2D_{2}=\|q\|_{\infty}+2\|q\|_{1}+\|q^{(1)}\|_{1}$, then

[TABLE]

Now, proceeding as in the proof of Theorem VI.1, we conclude that the estimate (119) holds with $C=\|\bm{D}\|\cosh(\|q\|_{1})$ where $\bm{D}=(D_{1},D_{2})^{\intercal}$.

∎

Finally, let us extend the preceding two results to the one-sided potentials:

Theorem VI.3**.**

Let $q\in\mathsf{E}_{d}(J)$ for some $d>0$ with $J=(-\infty,L]$. Let $\kappa_{1}=\|q\|_{\mathsf{L}^{1}(J)}$ and $\kappa_{\infty}>0$ be the constant such that $|q(x)|\leq\kappa_{\infty}e^{-2d|x|}$. Then, for every $\mu\in(0,d)$, the estimate

[TABLE]

holds with constants $C_{1}$ and $C_{2}$ given by

[TABLE]

where $\bm{D}=(\kappa_{1}^{2}/2,\kappa_{1})^{\intercal}$ and $\bm{E}=(\kappa^{2}_{\infty},d\kappa_{\infty})^{\intercal}$.

In addition, if $\partial_{x}q\in\mathsf{E}_{d}(J)$, then there exists a constant $C>0$ independent of $\zeta\in\mathbb{C}$ such that the estimate

[TABLE]

holds.

VI.1.1 Error due to domain-truncation

For the purpose of numerical solution of the ZS-problem posed on a unbounded domain, it is mandatory to choose a computational domain that is bounded. This requires truncation of the original unbounded domain to a bounded domain, say, $[-L_{-},L_{+}]$ where $L_{-},L_{+}>0$. Let us observe here that the estimates obtained in Theorem VI.3 can be improved slightly in order to give us a better control of the domain truncation error. Let $\mathscr{K}_{j}$ denote the Volterra integral operator corresponding to the kernel $\mathcal{K}_{j}$ for $j=1,2$ defined in (112). Set the domain to be $J=(-\infty,-L_{-}]$ and assume the conditions stated in the first part of Theorem VI.3 to be true. Then it can be shown that, for $\zeta\in\overline{\mathbb{C}}_{+}$, we have

[TABLE]

Now in any numerical scheme, one would take $(1,0)^{\intercal}$ as the initial value for the Jost solution $\bm{\phi}(x;\zeta)e^{i\zeta x}$ at $x=-L_{-}$. This step introduces an error which is bounded by $\max(\|\widetilde{P}_{1}\|_{\mathsf{L}^{\infty}(J)},\,\|\widetilde{P}_{2}\|_{\mathsf{L}^{\infty}(J)})$. Let $L_{-}>0$ be a free parameter and assume $q\in\mathsf{E}_{d}(\mathbb{R})$. Now, if we require the maximum error to be equal to $\epsilon>0$, then it suffices to have $\sinh[\|q\chi_{(-\infty,-L_{-}]}\|_{\mathsf{L}^{1}}]=\epsilon$ which works out to be

[TABLE]

Similar result can be obtained for truncation from the right side by using the property in Remark II.1.

VI.2 Discretization in the spectral domain

Let the grid points be as defined in Sec. III.1. In this section we discuss of the stability and convergence properties of the numerical methods developed in Sec. III.1. To this end, we closely follow the terminology introduced in Gautschi (2012) adapted to the problem at hand.

The general form of a one-step method for (18) can be stated as

[TABLE]

where dependence on the spectral parameter, $\zeta$, is suppressed. We keep the spectral parameter fixed in the following discussion or allow it to vary over any compact domain of $\mathbb{C}$. The function $\Lambda(x_{n};h)$ is referred to as the update function of the one-step method. The truncation error of this method is defined as

[TABLE]

with $\tilde{\bm{v}}(x)=\bm{y}$. A method is called consistent if $\lim_{h\rightarrow 0}\bm{T}(x,\bm{y};h)=0$. The necessary and sufficient condition for consistency in this case is $\Lambda(x;0)=\widetilde{U}(x)$. A method is said to have an order $p$ if, for some vector norm $\|\cdot\|$, $\|\bm{T}(x,\bm{y};h)\|\leq Ch^{p}$ holds uniformly over $\Omega\times\Gamma$ where $\Gamma\subset\mathbb{C}^{2}$ is a compact set and $C$ is independent of $x$, $\bm{y}$ and $h$. Let $\bm{x}_{h}=(x_{n})_{0\leq n\leq N}$ represent the grid. Let us introduce a vector-valued grid-function as $\bm{\mathfrak{u}}=\{\bm{u}_{n}\}_{n=0}^{N}$ where $\bm{u}_{n}\in\mathbb{C}^{2}$ such that value of $\bm{\mathfrak{u}}$ at $x_{n}$ is $\bm{u}_{n}$. The class of such solutions is denoted by $\mathsf{G}(\bm{x}_{h})$. Define the infinity-norm of any grid-function as

[TABLE]

In order to introduce the concept of stability of the one-step method, let us define the residue operator as

[TABLE]

for any grid-function $\bm{\mathfrak{u}}\in\mathsf{G}(\bm{x}_{h})$ and $n<N$ (we set $(\mathscr{R}_{h}\bm{\mathfrak{u}})_{N}=(\mathscr{R}_{h}\bm{\mathfrak{u}})_{N-1}$). A method is said to be stable if there exists a constant $C_{0}$ for $h_{0}>0$ such that for any two arbitrary grid-functions, $\bm{\mathfrak{u}},\,\bm{\mathfrak{w}}\in\mathsf{G}(\bm{x}_{h})$, we have

[TABLE]

for all $h\leq h_{0}$.

Remark VI.1**.**

The intuition behind this definition, as explained in Gautschi (2012) is as follows: if $\bm{\mathfrak{u}}\in\mathsf{G}(\bm{x}_{h})$ denotes the grid-function obtained by the one-step method using infinite-precision arithmetic (so that $\mathscr{R}_{h}\bm{\mathfrak{u}}=0$) and if $\bm{\mathfrak{w}}\in\mathsf{G}(\bm{x}_{h})$ denotes the grid-function obtained using finite-precision arithmetic (initial conditions being the same, i.e., $\bm{u}_{0}=\bm{w}_{0}$), then any stable method must yield $\|\bm{\mathfrak{u}}-\bm{\mathfrak{w}}\|_{\infty}=\mathop{\mathscr{O}}\left(\epsilon\right)$ where $\epsilon$ is the machine precision in the latter case.

Let $q\in\mathsf{BV}(\Omega)$, then there exist a constants $C$ and $h_{0}>0$ independent of $x$ and $h$ such that $\|\Lambda(x;h)\|<C$ for all $h\in[0,h_{0}]$, (where $\|\cdot\|$ is the induced matrix norm). This shows that for any two arbitrary vectors, $\bm{u},\bm{w}\in\mathbb{C}^{2}$ and $h\in[0,h_{0}]$, the Lipschitz condition,

[TABLE]

is satisfied. Therefore, the stability of the one-step method (123) easily follows from (Gautschi, 2012, Theorem 5.3.1). Further, for any grid-function $\bm{\mathfrak{u}}\in\mathsf{G}(\bm{x}_{h})$, we have

[TABLE]

Then using the inequality $(1+Ch)^{N}<e^{ChN}$, it follows that $\|\bm{\mathfrak{u}}\|_{\infty}\leq e^{C(L_{2}-L_{1})}\|\bm{u}_{0}\|$ which also guarantees the boundedness of the numerical solution when computed using infinite precision.

Finally, consistency and stability for any given one-step method imply global convergence. Moreover, if $\tilde{\bm{\mathfrak{v}}}=\{\tilde{\bm{v}}(x_{n})\}_{n=0}^{N}$, denotes the grid-function determined by the exact solution and $\bm{\mathfrak{u}}\in\mathsf{G}(\bm{x}_{h})$ be any grid-function obtained using the one-step method (123) with initial condition $\bm{u}_{0}=\tilde{\bm{v}}(x_{0})$, then $\|\bm{\mathfrak{u}}-\tilde{\bm{\mathfrak{v}}}\|_{\infty}=\mathop{\mathscr{O}}\left(h^{p}\right)$ where $p$ is the order of the one-step method (Gautschi, 2012, Theorem 5.3.2).

VI.2.1 Implicit Euler method

Continuing from Sec. III.1.2, we have

[TABLE]

which determines the update function to be

[TABLE]

It is easy to verify that $\Lambda(x_{n};0)=\widetilde{U}_{n}$, therefore, the method is consistent. Using the Taylor’s theorem

[TABLE]

where $x\leq x^{\prime}\leq x+h$ and $\tilde{\bm{v}}(x)=\bm{y}$. Assuming that $q(x)\in\mathsf{C}^{1}_{0}(\Omega)$, we have $\partial^{2}_{x}\tilde{\bm{v}}=e^{i\sigma_{3}\zeta x}(\partial_{x}U+U^{2}+2i\zeta[\sigma_{3},U])\bm{v}$, therefore, the order of the method is $p=1$. If the Jost solution under consideration is ${\bm{v}}=\bm{\phi}$, then $\|e^{i\zeta x}\bm{\phi}\|$ is bounded for $\zeta\in\overline{\mathbb{C}}_{+}$ (see Theorem VI.1), consequently, the truncation error coefficient to the leading order in $\zeta$ is $|\zeta|he^{2\operatorname{Im}(\zeta)x}\|[\sigma_{3},U]\|$. Evidently, the method is stable which together with its consistency imply convergence (with order $p=1$).

VI.2.2 Trapezoidal rule

Continuing from Sec. III.1.3, we have

[TABLE]

so that the update function is given by

[TABLE]

Again, it is easy to verify that $\Lambda(x_{n};0)=\widetilde{U}_{n}$, therefore, the method is consistent. Using the Taylor’s theorem

[TABLE]

with $\tilde{\bm{v}}(x)=\bm{y}$. Assuming that $q(x)\in\mathsf{C}^{2}_{0}(\Omega)$, we have

[TABLE]

therefore, the order of the method is $p=2$. Again, if the Jost solution under consideration is ${\bm{v}}=\bm{\phi}$, then $\|e^{i\zeta x}\bm{\phi}\|$ is bounded for $\zeta\in\overline{\mathbb{C}}_{+}$ (see Theorem VI.1), consequently, the truncation coefficient to the leading order in $\zeta$ is $|\zeta|^{2}h^{2}e^{2\operatorname{Im}(\zeta)x}\|[\sigma_{3},[\sigma_{3},U]]\|/3$. Evidently, the method is stable which together with its consistency imply convergence (with order $p=2$).

VI.2.3 Split-Magnus method

For the convergence analysis of the split-Magnus method described in Sec. IV.2, let us observe that an equivalent form of the integrator is

[TABLE]

Using Taylor’s theorem for matrix functions, we have

[TABLE]

with $\tilde{\bm{v}}(x)=\bm{y}$. Assuming $U$ to be twice differentiable on $[x,x+h]$, we conclude that the order of the method is $p=2$. Further, this one-step method is consistent and stable, therefore, also convergent for fixed $\zeta$ (or $\zeta$ varying in a compact domain). The truncation error coefficient to the leading order in $\zeta$ is $|\zeta|^{2}h^{2}e^{2\operatorname{Im}(\zeta)x}\|[\sigma_{3},[\sigma_{3},U]]\|/6$. This value can be seen to be twice as small as that of the trapezoidal rule. Let us note that it does not seem straightforward to determine which of the two one-step methods has smaller total truncation error coefficient (for fixed $\zeta$); however, the trapezoidal rule appears to exhibit smaller total truncation error in the numerical tests.

VI.3 Computation of norming constants

In Sec. III.2.1, it was stated that the computation of norming constants from the discrete $b$-coefficients $b_{N}(z^{2})$ is ill-conditioned. This can be attributed to the nature of the truncation error coefficients in the underlying one-step method for complex values of $\zeta$. It is evidenced by the presence of a factor of the form $\exp[2\operatorname{Im}(\zeta)x]$ in the truncation error coefficient which tends to grow for $x>0$ (see Sec. VI.2). Therefore, it is better to “truncate” the scattering potential at the origin151515If the growth behavior of the potential is known before-hand, then it possible to choose an optimal point of truncation. and solve the corresponding one-sided ZS-problems as discussed in Sec. III.2.1. Finally, let us note that there are other discretization schemes such as the exponential time-differencing (ETD) scheme Cox and Matthews (2002) which may alleviate these problems; however, it may come at a cost of increased operational complexity. These ideas will be explored in a future publication.

VI.4 Lubich’s method

Starting from the functions $a(\zeta)$ and $\breve{b}(\zeta)$ analytic in the upper-half of the complex plane, Lubich’s construction as described in Sec. III.6 allows us to compute the polynomials associated with the discrete scattering coefficients $\bm{P}_{N}(z^{2})=\{\bm{P}(z^{2})\}_{N}$. Note that, in the preceding section, we discussed the necessary and sufficient condition for discrete inverse scattering with polynomials (which can be seen as a finite-support sequence of coefficients). However, Lubich’s method yields infinite series that needs to be truncated. Therefore, the compatibility of Lubich’s construction with the layer-peeling algorithm cannot be studied within the framework of finite-support sequences. However, it is possible to determine if $\bm{P}(z^{2})$ can be associated to a Jost solution prior to truncation of the series. If the coefficients of the series decay sufficiently fast, the truncation introduces a negligible error so that the layer-peeling criteria can be satisfied to a sufficient degree of accuracy.

Let us first consider the case of compactly supported potential. Define the vector $\bm{P}(z^{2})=(P_{1},P_{2})^{\intercal}$ as

[TABLE]

which can be expanded into a convergent Taylor series as in Sec. III.6 on account of the analyticity of the scattering coefficients over whole of the complex plane. Further note that

[TABLE]

Therefore, for $z\in\mathbb{T}$, we have161616For sufficiently small $h$, it can be verified that $|P_{1,0}|\neq 0$. Other conditions pertaining to the specific discretization schemes can be explicitly verified using the results in Sec. III.6.

[TABLE]

Note that here we have used the fact that $a(\zeta)\overline{a}(\zeta)+{b}(\zeta)\overline{b}(\zeta)=1$ for all $\zeta\in\mathbb{C}$, however, such a relationship would not hold if we relax the requirement of compact support of the potential.

Let $f(\zeta)$ denote either $a(\zeta)-1$ or $\breve{b}(\zeta)$. When $f(\zeta)$ is analytic in the upper-half of the complex plane, then on any compact domain $\Gamma\subset\overline{\mathbb{C}}_{+}$ the functions can be regarded as Lipschitz continuous. Observing that $\delta(e^{-h})/h=1+\mathop{\mathscr{O}}\left(h^{p}\right)$ where $p=1$ for BDF1 and $p=2$ for TR, we have

[TABLE]

on any compact domain of $\Gamma\subset\overline{\mathbb{C}}_{+}$ and $h\in(0,\bar{h}]$ ($\bar{h}>0$) where $C>0$ depends only on $\Gamma$ and $\bar{h}$. Therefore, using the estimate (129) and the Lipschitz continuity of $f(\zeta)$ one can assert that there exists a constant $C^{\prime}>0$ for a given $\Gamma$ and $h_{0}$ such that Lubich (1994)

[TABLE]

Therefore, the Wronskian relationship, $|a(\xi)|^{2}+|b(\xi)|^{2}=1$ for $\xi\in\mathbb{R}$ can only be satisfied upto $\mathop{\mathscr{O}}\left(h^{p}\right)$ on any bounded interval in $\mathbb{R}$.

Finally, as far as the truncation of the infinite series is concerned, let us note that for the kind problems considered in this article, Lubich’s method is applied to rational functions with known poles in $\mathbb{C}_{-}$ which makes it easy to determine the decay behavior of these coefficients using the method of partial-fractions (see Sec. III.7).

VI.5 Darboux transformation

In this section, we study convergence behavior of the Darboux transformation with numerically computed Jost solutions. Continuing from Sec. II.3, let $(\zeta_{k},b_{k})\in\mathfrak{S}_{K}$ denote the discrete eigenvalue and the corresponding norming constant. Define the Vandermonde matrix

[TABLE]

the diagonal matrix $\Gamma=\text{diag}(\gamma_{1},\gamma_{2},\ldots,\gamma_{K})$ and the vectors

[TABLE]

where

[TABLE]

The unknown Darboux coefficients can be put into the vector form

[TABLE]

then the linear system of equations (10) which determines the coefficients of the Darboux matrix can be written as

[TABLE]

Note that the quantities $\bm{f}$ and $F$ are known exactly while $\Gamma$ (and in turn $\bm{g}$) is determined only up to $\mathop{\mathscr{O}}\left(h^{p}\right)$, where $p$ is the order of convergence of the one-step method. Let $\|\cdot\|$ denote the Euclidean norm for vectors and the induced spectral norm for matrices. Define $\kappa(\mathcal{W})=\|\mathcal{W}^{-1}\|\cdot\|\mathcal{W}\|$ to be the condition number of $\mathcal{W}$; then, under the assumption $\|\mathcal{W}^{-1}\|\cdot\|\Delta\mathcal{W}\|<1$ (which can be satisfied for sufficiently small $h$), the standard perturbation theory (Lancaster and Tismenetsky, 1985, Chap. 11) yields the estimate

[TABLE]

Given that the perturbations are of $\mathop{\mathscr{O}}\left(h^{p}\right)$, from above equation it follows that the coefficients of the Darboux matrix can be determined up to $\mathop{\mathscr{O}}\left(h^{p}\right)$.

In order to determine the convergence behavior of the fast Darboux transformation (FDT) algorithm as described in Sec. III.8.1, we need to study the convergence of the corresponding Lubich coefficients. To this end, let us denote by $\widetilde{D}_{K}(\zeta)$ the approximation to the Darboux matrix ${D}_{K}(\zeta)$ (for the sake of brevity, we suppress the dependence on $x$, $t$ and $\mathfrak{S}_{K}$). Now, using the partial fraction expansion as in (53), we have

[TABLE]

In order to establish the relationship between the error in the Darboux matrix as stated above and the error in the coefficients of the Darboux matrix, we need the following lemma:

Lemma VI.4**.**

For a given discrete spectrum $\mathfrak{S}_{K}$ where $K$ is finite, the inequality

[TABLE]

holds for any $\zeta\in\mathbb{C}$ where

[TABLE]

Proof.

From the the definition of the Darboux matrix, we have

[TABLE]

for $\zeta\in\mathbb{C}$. Now, using the Cauchy-Schwartz inequality, we obtain

[TABLE]

Note that this inequality does not change on replacing the spectral norm ($\|\cdot\|$) with the Euclidean norm ($\|\cdot\|_{E}$) and it is easy to see

[TABLE]

which concludes the proof. ∎

Let $\mathcal{D}_{K}(\tau)$ and $\widetilde{\mathcal{D}}_{K}(\tau)$ denote the inverse Fourier-Laplace transform of $\mu_{K}(\zeta)D_{K}(\zeta)-\sigma_{0}$ and $\mu_{K}(\zeta)\widetilde{D}_{K}(\zeta)-\sigma_{0}$, respectively; then, we have the following proposition for the rate of convergence:

Proposition VI.5**.**

Consider the discrete spectrum $\mathfrak{S}_{K}$ with finite $K$. If $\|\Delta\bm{D}\|=\mathop{\mathscr{O}}\left(h^{p}\right)$ where $p$ is order of the underlying one-step method, then

[TABLE]

Proof.

Let the set of eigenvalues be $\mathfrak{E}_{K}$ corresponding to $\mathfrak{S}_{K}$ and define

[TABLE]

where $G_{K}$ is as defined in the forgoing lemma. From (135) and the forgoing lemma, we have

[TABLE]

where $\eta_{k}=\operatorname{Im}{\zeta_{k}}>0$. The result follows by setting $\|\Delta\bm{D}\|=\mathop{\mathscr{O}}\left(h^{p}\right)$. ∎

Let the matrix-valued Lubich coefficients for ${D}_{K}(\zeta)$ and $\widetilde{D}_{K}(\zeta)$ be defined as

[TABLE]

respectively. The zeroth Lubich coefficient is obtained by evaluating the Darboux matrix at $\zeta=i\delta(0)/2h$. Therefore,

[TABLE]

leading to $\|\Lambda_{0}-\tilde{\Lambda}_{0}\|=\mathop{\mathscr{O}}\left(h^{p+1}\right)$. Using the properties of the Lubich coefficients and the forgoing proposition, it follows that

[TABLE]

VII Numerical Tests

In this section, we present several numerical tests to demonstrate the performance of the numerical algorithms developed in this paper. For better numerical conditioning, we scale the scattering potential $q(x)$ of the ZS-problem by a suitable scaling parameter such that $\|q\|_{\mathsf{L}^{2}}$ is unity or close to unity. Let us briefly review the effect of this scaling to (3): For some $\kappa>0$, let $V(y)={U}(x)/\kappa$, $y=\kappa x$ and $\lambda=\zeta/\kappa$ then

[TABLE]

where $\bm{w}(y;\lambda)=\bm{v}(y/\kappa;\lambda\kappa)$.

For the sake of clarity, let us specify the acronyms used to denote the one-step methods considered in this article for testing: implicit Euler method (BDF1), trapezoidal rule (TR), Magnus method with one-point Gauss quadrature (MG1) and split-Magnus method (SM). The main focus of this section is to study the dependence of the total numerical error on the free parameters of a given algorithm together with its total run-time. In particular, we have considered the test cases that test the performance of the new methods introduced in this article against the so called benchmarking methods (MG1 and SM) wherever possible. In all of the test cases described below, $N$ represents the number of samples which is taken from the set $\mathfrak{N}=\{2^{j},\,j=10,11,\ldots,20\}$.

VII.1 Examples

Our test cases are derived from following examples for which the exact value of the quantities to be analyzed are known in a closed form or can be evaluated to the machine precision by a known method.

VII.1.1 Multi-solitons

Define a sequence of angles for $J\in\mathbb{Z}_{+}$ by choosing $\Delta\theta=(\pi-2\theta_{0})/(J-1),\,\theta_{0}>0$, and

[TABLE]

so that $\theta_{j}\in[\theta_{0},\pi-\theta_{0}]$. Then the eigenvalues for our numerical experiment are chosen as

[TABLE]

The norming constants are chosen as

[TABLE]

Here we choose, $\theta_{0}=\pi/3$ and $J=4$. Then we consider a sequence of discrete spectra defined as

[TABLE]

where $K=4,8,\ldots,32$. Let $\mathfrak{E}_{K}$ be the set of all the eigenvalues. The potential can be computed with machine precision using the CDT algorithm which is taken as the reference for error analysis in this case. For fixed $K$, the eigenvalues are scaled by the scaling parameter $\kappa=2(\sum_{k=0}^{K}\operatorname{Im}\zeta_{k})^{1/2}$. Let $\eta_{\text{min}}=\min_{\zeta\in\mathfrak{E}_{K}}\operatorname{Im}\zeta$, then the computational domain for this example is chosen as $[-L,\,L]$ where $L={11\kappa}/\eta_{\text{min}}$.

VII.1.2 Secant-hyperbolic potential

The exact solution of the ZS-problem for the secant-hyperbolic potential was first reported in Satsuma and Yajima (1974). We summarize the results required for our purpose as follows: Let the potential be written as

[TABLE]

where $A$ is referred to as the amplitude. The scattering coefficients are then given by

[TABLE]

The eigenvalues are given by

[TABLE]

where $K$ is the largest integer smaller than $\tilde{A}=(A+1/2)$. Putting $\tilde{A}_{f}=\tilde{A}-K$, the non-integer part of $\tilde{A}$, the $a$-coefficient can be written as a product of solitonic and radiative parts as follows

[TABLE]

Note that $a_{R}(\zeta)$ belongs to a secant-hyperbolic potential with amplitude $A_{R}=\tilde{A}_{f}-1/2\,\,(>0)$. The corresponding norming constants are given by $b_{k}=(-1)^{k}$.

This example allows one to test the CDT and the FDT algorithms where the seed potential can be taken as $q_{0}(x)=A_{R}\operatorname{sech}(x)$ and the sequence of discrete spectra to be added,

[TABLE]

where we set $A_{R}=0.4$. Corresponding to $\mathfrak{S}_{K}$, the amplitude of the augmented secant-hyperbolic potential is given by $A=0.4+K$ and $\tilde{A}=0.9+K$. As in the last example, for fixed $K$, the eigenvalues are scaled by the scaling parameter given by $\kappa=2(\sum_{k=0}^{K}\operatorname{Im}\zeta_{k})^{1/2}$.

In order to choose the computational domain $[-L,\,L]$ for the sech-potential (142) with the aforementioned scaling, we can use the relation (122). Choosing $\eta_{\text{min}}=\min_{\zeta\in\mathfrak{E}_{K}}\operatorname{Im}\zeta$ where $\mathfrak{E}_{K}$ is the set of all the eigenvalues and observing that

[TABLE]

we have $L\approx[\eta_{\text{min}}\log(2A/\epsilon)](\kappa/\eta_{\text{min}})$ which rounds to $L\approx 30(\kappa/\eta_{\text{min}})$ for $\epsilon=10^{-12}$.

VII.2 Test cases

VII.2.1 Discrete spectrum

For multi-soliton potentials described in Sec. VII.1.1, we test the convergence behavior with regard to the discrete spectrum for various discretization schemes, namely, BDF1, TR, SM and MG1. For the convergence behavior of the numerically computed norming constants, we assume that the eigenvalues are known exactly. The error in the norming constants is quantified by

[TABLE]

For the convergence behavior with regard to the eigenvalues, we compute $a^{(\text{num.})}(\zeta_{k})$ where the $a$-coefficient is computed numerically. The error is then quantified by

[TABLE]

For MG1, the limit is evaluated by setting $\eta=100$. For others $\lim_{\eta\rightarrow\infty}a(i\eta)=P^{(N)}_{1,0}$. Except for MG1, all other schemes are implemented using the fast forward scattering algorithm (see Sec. III.5.1). The computation of the norming constant is discussed in Sec. III.2.1 and Sec. IV.

VII.2.2 Multi-soliton potential

In this test case, we carry out the convergence analysis and a comparison of run-time (per sample) of different variants of the FDT algorithm for multi-solitons as described in Sec. VII.1.1. Note that the CDT algorithm in this case gives the exact potential which allows us to compute the total numerical error for the FDT algorithm for arbitrary discrete spectra. The error is quantified by

[TABLE]

where the integrals are evaluated numerically using the trapezoidal rule.

The different variants of the FDT algorithm are described as follows: any one-step method for the ZS-problem can be combined with any one-step method for the Lubich coefficients to obtain the FDT algorithm. In particular the relevant combinations are: BDF1-BDF1, BDF1-TR and TR-TR. We also consider the partial-fraction variant of the TR-TR combination which is labeled as TR-TR-PF. Note that the combination of a first order method for the ZS-problem with a second order method for the Lubich coefficients or vice versa should lead to a first order FDT algorithm. A second order method for the ZS-problem must be combined with a second order method or higher for the Lubich coefficients in order to obtain a second order FDT algorithm.

Parameters for the Lubich method are as follows: $M=8N$ and $N_{\text{th}}=N/8$ (for the PF-variant). For the Cauchy integral, the radius of the circular contour is $\varrho=\exp[-8/(N/2-1)]$.

VII.2.3 General Darboux transformation

In this test case, we carry out the convergence analysis and a comparison of run-time (per sample) of different variants of the CDT/FDT algorithm for the secant-hyperbolic potential as described in Sec. VII.1.2. Note that, in the case of the secant-hyperbolic potential, the soliton-free seed potential as well as the augmented potential can be stated in a closed form.

The variants of the CDT/FDT algorithm are determined by the underlying one-step method. Unlike in the last test case (Sec. VII.2.2), Lubich method uses the same one-step method as that of the ZS-problem. The total numerical error is quantified by (146). Parameters for the Lubich method are the same as in the last test case.

VII.3 Results and discussion

VII.3.1 Discrete spectrum

For a given multi-soliton, this test case was designed to assess the performance of the discretization schemes, namely, BDF1, TR, SM and MG1 in the determination of the discrete spectrum. The results are plotted in Fig. 7 where it can be easily seen that all the methods considered show convergence at a rate that is determined by the underlying one-step method. However, the rate of convergence of BDF1 with regard to discrete eigenvalues seems to be better than expected as evident from the plots in the bottom row of Fig 7. The overall accuracy of MG1 is evidently superior to that of others while TR turns out to be a close second.

VII.3.2 Multi-soliton potential

This test case was designed to study the convergence and run-time behavior of different variants of the FDT algorithm for multi-solitons. The results for fixed number of eigenvalues and varying number of samples is shown in Fig. 8. The second order of convergence of the schemes BDF1-TR, TR-TR and TR-TR-PF can be identified from the plots in the top row of Fig. 8. The scheme BDF1-TR shows better rate of convergence than expected.

The run-time behavior of CDT for fixed number of eigenvalues is clearly superior to all of the method as evident from Fig. 8 (bottom row). The scheme TR-PF however becomes very competitive with the CDT algorithm.

The relative error and the run-time (per sample) as a function of the number of eigenvalues keeping the number of samples fixed is shown in Fig. 9. Here, the FDT algorithm outperforms the CDT algorithm with TR-TR-PF variant being the fastest as evident from the plots in the bottom row of Fig. 9. It is interesting to note that the relative error as a function of number of eigenvalues as shown in plots at the top row of Fig. 9 exhibits exponentially increasing behavior. This puts an upper limit to the number of eigenvalues that can be handled with the FDT algorithm within a given precision171717Note that the CDT algorithm also suffers from this drawback. However, in order to determine the upper limit for the CDT algorithm, one requires an implementation which employs a variable precision arithmetic. This program is not followed in this article..

VII.3.3 General Darboux transformation

This test case was designed to study the convergence and run-time behavior of different variants of the CDT/FDT algorithm for a soliton-free seed potential. The results for fixed number of eigenvalues (that are meant to be added) and varying number of samples is shown in Fig. 12. The second order of convergence of the TR variant of the CDT/FDT algorithm can be identified from the plots in the top row of Fig. 12. However, the TR variant of the CDT algorithm performs not only worse as compared to that of the FDT algorithm but it also becomes unstable with increasing number of eigenvalues. Further, unlike the CDT algorithm, the BDF1 and TR variant of FDT shows convergence (at an expected rate) with increasing number of samples.

The run-time behavior of CDT for fixed number of eigenvalues is clearly superior to that of FDT as evident from Fig. 13 (bottom row). The scheme TR-PF however becomes very competitive to the CDT-TR variant. Note that CDT in this case is reliable only for small number of eigenvalues.

The relative error and the run-time (per sample) as a function of the number of eigenvalues keeping the number of samples fixed is shown in Fig. 13. Here, the FDT algorithm outperforms the CDT algorithm with TR-TR-PF variant being the fastest as evident from the plots in the bottom row of Fig. 13. As in the last test case, the relative error as a function of number of eigenvalues as shown in plots at the top row of Fig. 13 exhibits exponentially increasing behavior. Note that FDT not only outperforms CDT in terms of accuracy, it also exhibits superior numerical conditioning with increasing number of eigenvalues as evident from Fig. 13 (top row).

VIII Conclusion

To conclude, we have presented a systematic approach to discretize the non-Hermitian Zakharov-Shabat (ZS) problem which is based on exponential one-step methods. The discrete framework thus obtained is amenable to FFT-based fast polynomial arithmetic and also admits of a layer-peeling property. In this setting we have presented different variants of a fast forward/inverse SU(2)-nonlinear Fourier transformation (NFT) algorithm. As a first step towards developing a general fast inverse NFT, we have presented several ways to obtain a fast Darboux transformation (FDT) algorithm with an operational complexity of $\mathop{\mathscr{O}}\left(KN+N\log^{2}N\right)$ where $K$ is the number of eigenvalues to be added to a seed potential and $N$ is the number of samples of the potential. This algorithm exhibits an order of convergence that matches the underlying exponential one-step method. In particular, if one uses the trapezoidal rule of integration, the order of convergence is $\mathop{\mathscr{O}}\left(N^{-2}\right)$. The strength of this algorithm was demonstrated by exhaustive numerical tests where we could successfully add $32$ eigenvalues to a soliton-free seed potential. It must be noted that the FDT algorithm shows a promising route to a fast inverse NFT which is confirmed empirically in Vaibhav and Wahls (2017)–this forms the subject matter of a sequel to this paper.

Furthermore, we have also presented a second approach that naively tries to mimic the classical Darboux transformation (CDT) scheme in the discrete framework developed for the ZS-problem with an arbitrary seed potential. This algorithm affords a complexity of $\mathop{\mathscr{O}}\left(K^{2}N\right)$; however, it turns out to be less accurate and numerically unstable beyond certain number of eigenvalues.

Finally, let us emphasize that, based on the ideas presented in this paper and drawing on the pioneering work of Lubich on convolution quadrature, it seems plausible to anticipate the existence of higher-order convergent fast forward/inverse NFT algorithms using (exponential) linear multistep methods–we hope to return to this theme in the future.

Appendix A Lubich coefficients for rational functions with simple poles

Consider the simplest case of a rational function with a simple pole $E(\zeta)=(\zeta-\zeta_{0})^{-1}$ where $\operatorname{Im}\zeta_{0}<0$. It satisfies the kind of growth estimated stated in (47). The inverse Fourier-Laplace transform is given by $e(\tau)=-ie^{-i\zeta_{0}\tau}$. The Lubich coefficients corresponding to the trapezoidal rule is defined through

[TABLE]

where we note that $|{1-i\zeta_{0}h}|/|{1+i\zeta_{0}h}|<1$ on account of $\operatorname{Im}\zeta_{0}<0$. The Lubich coefficient $e_{k}$ is defined as the coefficient of $z^{2k}$ in the RHS of (147) which can be worked out explicitly: $e_{0}={-ih}/(1+i\zeta_{0}h)$ and, for $k>0$,

[TABLE]

Note that when $\operatorname{Re}\zeta_{0}=0$, then we may restrict $h\in(0,\bar{h}]$ so that $1+\eta_{0}h\neq 0$ where $\eta_{0}=\operatorname{Im}{\zeta}$. In the following, we wish to study the error involved in replacing $e_{k}$ with $-2ihe^{-2i\zeta_{0}hk}$. For $k>0$, this difference is given by

[TABLE]

Using the $[1/1]$-Padé approximant Baker and Graves-Morris (1981), we have

[TABLE]

Next, let us show that, for $h\in(0,\bar{h}]$, there exists a positive integer $n>1$ dependent only on $\bar{h}$ such that

[TABLE]

Recalling $\eta_{0}=\operatorname{Im}\zeta_{0}$, let $n$ be chosen such that

[TABLE]

From the inequality (Olver et al., 2010, Chap. 4)

[TABLE]

it follows that $n>1$. Also, observing

[TABLE]

it suffices to choose $n$ to be the smallest integer greater than $(1+|\zeta_{0}|\bar{h})^{2}$. Now, using standard inequalities for exponential function (Olver et al., 2010, Chap. 4), we have

[TABLE]

where

[TABLE]

Finally, using the last estimate and from (149), we conclude

[TABLE]

where $C^{\prime}>0$ is independent of $h$ and $k$. If in addition $(-\eta_{0}\bar{h})<1$, then for $h\in(0,\bar{h}]$ and $k>0$, we may write

[TABLE]

where $C=C^{\prime}/(1+\eta_{0}\bar{h})$ is independent of $h$ and $k$. Using this estimate, let us now show that one can make an informed choice of the parameter $N_{\text{th}}\in\mathbb{Z}_{+}$ introduced in sections III.6.1 and III.7 in connection with the partial-fraction variant of FDT. To this end, we start with

[TABLE]

where $\epsilon$ is a positive number less than unity. Putting $N_{\text{th}}=N/m$ and using $h\sim 2L/N$ where $2L=L_{2}-L_{1}$ and $N\in\mathbb{Z}_{+}$, we have

[TABLE]

Setting $\epsilon=e^{-1}$ and for $|\zeta_{0}|\bar{h}<1$ one can set $n=4$ so that $m\sim|\eta_{0}|L$. In case of multi-solitons, we would like to tune the parameter $N_{\text{th}}$ with respect to the eigenvalue with the smallest imaginary part. Here $\zeta_{0}$ must be replaced by the complex conjugate of the eigenvalue with smallest imaginary part. For example, if the smallest of all the imaginary parts of eigenvalues is unity and $L=10$, we should choose $m=10$ (or, $m=8$ so that $N_{\text{th}}$ is a power of $2$ when $N$ is a power of $2$).