Numerical Methods for Fast Nonlinear Fourier Transformation, Part I: Exponential Runge-Kutta and Linear Multistep Methods

TL;DR

This paper analyzes exponential Runge-Kutta and linear multistep methods for fast nonlinear Fourier transformation, demonstrating their stability, convergence, and computational complexity, and highlighting their suitability for different boundary conditions.

Contribution

It provides a theoretical analysis of exponential RK and LM methods for NFT, establishing their stability, convergence, and computational efficiency, and discusses their applicability to various boundary conditions.

Findings

Both methods achieve $ ext{O}(N ext{log}^2 N)$ complexity.

RK methods require equispaced nodes for FFT-based acceleration.

LM methods handle vanishing boundary conditions without starting procedures.

Abstract

The main objective of this series of papers is to explore the entire landscape of numerical methods for fast nonlinear Fourier transformation (NFT) within the class of integrators known as the exponential integrators. In this paper, we explore the theoretical aspects of exponential Runge-Kutta (RK) and linear multistep (LM) methods, in particular, the stability and convergence of these methods via the transfer matrix formulation. The analysis carried out in the paper shows that while the exponential LM methods are naturally amenable to FFT-based fast polynomial arithmetic, the RK methods require equispaced nodes to achieve that. Therefore, each these family of methods is capable of yielding a family of fast NFT algorithms such that the scattering coefficients can be computed with a complexity of $\mathscr{O}(N\log^2N)$ and a rate of convergence given by $\mathscr{O}(N^{-p})$ where $N$ is the number of samples of the signal and $p$ is order of the underlying discretization scheme. Further, while RK methods can accommodate vanishing as well as periodic boundary conditions, the LM methods can only handle the former type of boundary conditions without requiring a starting procedure. The ideas presented in this paper extend naturally to the family of integrators known as general linear methods which will be explored in a forthcoming paper.

Tables

Table 1. Table 1: Explicit Adams Method

Method	$𝜷$	Order
EA1	$(1,0)$	$1$
EA2	$(-\frac{1}{2},\frac{3}{2},0)$	$2$
EA3	$(\frac{5}{12},-\frac{16}{12},\frac{23}{12},0)$	$3$
EA4	$(-\frac{9}{24},\frac{37}{24},-\frac{59}{24},\frac{55}{24},0)$	$4$
EA5	$(\frac{251}{720},-\frac{1274}{720},-\frac{2616}{720},-\frac{2774}{720},\frac{1901}{720},0)$	$5$

Table 2. Table 2: Implicit Adams Method

Method	$𝜷$	Order
IA1	$(\frac{1}{2},\frac{1}{2})$	$2$
IA2	$(-\frac{1}{12},\frac{8}{12},\frac{5}{12})$	$3$
IA3	$(\frac{1}{24},-\frac{5}{24},\frac{19}{24},\frac{9}{24})$	$4$
IA4	$(-\frac{19}{720},\frac{106}{720},-\frac{264}{720},\frac{646}{720},\frac{251}{720})$	$5$

Table 3. Table 3: Backward differentiation formula

Method	$𝜶$	$\beta$
BDF1	$(-1,1)$	$1$
BDF2	$(\frac{1}{3},-\frac{4}{3},1)$	$\frac{2}{3}$
BDF3	$(-\frac{2}{11},\frac{9}{11},-\frac{18}{11},1)$	$\frac{6}{11}$
BDF4	$(\frac{3}{25},-\frac{16}{25},\frac{36}{25},-\frac{48}{25},1)$	$\frac{12}{25}$
BDF5	$(-\frac{12}{137},\frac{75}{137},-\frac{200}{137},\frac{300}{137},-\frac{300}{137},1)$	$\frac{60}{137}$
BDF6	$(\frac{10}{147},-\frac{72}{147},\frac{225}{147},-\frac{400}{147},\frac{450}{147},-\frac{360}{147},1)$	$\frac{60}{147}$

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 ,$$

$${\mathbf{v}}_t=-i\zeta {\sigma }_3\mathbf{v}+U(t)\mathbf{v},$$

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

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

$$\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}$$

$$|b(\zeta )|\le \sinh (\mu )\times \{\begin{array}{ll}{e}^{2T\mathrm{Im}\zeta },& \zeta \in {\overline{\mathbb{C}}}_{+},\\ e^{-2T\mathrm{Im}\zeta },& \zeta \in {\mathbb{C}}_{-},\end{array}$$

$$|\tilde{a}(\zeta )|\le [\cosh (\mu )-1]\times \{\begin{array}{ll}{e}^{2T\mathrm{Im}\zeta },& \zeta \in {\overline{\mathbb{C}}}_{+},\\ e^{-2T\mathrm{Im}\zeta },& \zeta \in {\mathbb{C}}_{-},\end{array}$$

$$\mathrm{supp}{F}^{-1}[f]\subset [-2T,2T],$$

$$\overline{\mathbf{\psi }}(t;\zeta )=\left(\begin{array}{c}{\psi }_2^{\ast }(t;{\zeta }^{\ast })\\ \kappa {\psi }_1^{\ast }(t;{\zeta }^{\ast })\end{array}\right)=\{\begin{array}{ll}{\sigma }_1{\mathbf{\psi }}^{\ast }(t;{\zeta }^{\ast }),& \kappa =+1,\\ i{\sigma }_2{\mathbf{\psi }}^{\ast }(t;{\zeta }^{\ast }),& \kappa =-1,\end{array}$$

$$\overline{\mathbf{\varphi }}(t;\zeta )=-\left(\begin{array}{c}\kappa {\varphi }_2^{\ast }(t;{\zeta }^{\ast })\\ {\varphi }_1^{\ast }(t;{\zeta }^{\ast })\end{array}\right)=\{\begin{array}{ll}-{\sigma }_1{\mathbf{\varphi }}^{\ast }(t;{\zeta }^{\ast }),& \kappa =+1,\\ i{\sigma }_2{\mathbf{\varphi }}^{\ast }(t;{\zeta }^{\ast }),& \kappa =-1;\end{array}$$

$$\overline{a}(\zeta )=a^{\ast }({\zeta }^{\ast }),\overline{b}(\zeta )=-\kappa b^{\ast }({\zeta }^{\ast }).$$

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

$$\mathbf{\varphi }(t;\zeta )=\left(\begin{array}{c}{\varphi }_1(t;\zeta )\\ {\varphi }_2(t;\zeta )\end{array}\right),\overline{\mathbf{\varphi }}(t;\zeta )=\left(\begin{array}{c}{\varphi }_2^{\ast }(t;{\zeta }^{\ast })\\ \kappa {\varphi }_1^{\ast }(t;{\zeta }^{\ast })\end{array}\right),$$

$$\mathbf{\varphi }(t_0;\zeta )=\left(\begin{array}{c}1\\ 0\end{array}\right).$$

$$\mathbf{\varphi }(t+2T)=a(\zeta ;t_0)\mathbf{\varphi }(t;\zeta )+b(\zeta ;t_0)\overline{\mathbf{\varphi }}(t;\zeta ),$$

$$\begin{array}{rl}& a(\zeta ;t_0)={\varphi }_1(t_0+2T;\zeta ),\\ & b(\zeta ;t_0)=\kappa {\varphi }_2(t_0+2T;\zeta ).\end{array}$$

$$\Phi (\zeta ;t_0)={\left(\mathbf{\varphi },\kappa \overline{\mathbf{\varphi }}\right)}_{t=t_0+2T}=\left(\begin{array}{cc}{\varphi }_1(t_0+2T;\zeta )& \kappa {\varphi }_2^{\ast }(t_0+2T;{\zeta }^{\ast })\\ {\varphi }_2(t_0+2T;\zeta )& {\varphi }_1^{\ast }(t_0+2T;{\zeta }^{\ast })\end{array}\right),$$

$$\mathbf{\psi }(t+2T;\zeta )=\lambda \mathbf{\psi }(t;\zeta ).$$

$$\mathbf{\psi }(t)=c(\zeta ;t_0)\mathbf{\varphi }(t;\zeta )+d(\zeta ;t_0)\overline{\mathbf{\varphi }}(t;\zeta ),$$

$$\Phi (\zeta ;t_0)\left(\begin{array}{c}c\\ d\end{array}\right)=\lambda \left(\begin{array}{c}c\\ d\end{array}\right).$$

$$\det \left[\Phi (\zeta ;t_0)-\lambda {\sigma }_0\right]={\lambda }^2-\mathrm{tr}[\Phi ]\lambda +1=0.$$

$$a_R(\zeta )=\frac{1}{2}\left[a(\zeta ;t_0)+\overline{a}(\zeta ;t_0)\right]=\frac{1}{2}\mathrm{tr}[\Phi ](\zeta ),$$

$$a_I(\zeta ;t_0)=\frac{1}{2i}\left[a(\zeta ;t_0)-\overline{a}(\zeta ;t_0)\right].$$

$${\tilde{\mathbf{v}}}_t=\tilde{U}\tilde{\mathbf{v}},\tilde{U}=e^{i{\sigma }_3\zeta t}U{e}^{-i{\sigma }_3\zeta t},$$

$$\begin{array}[]{c|cccc}c_{1}&a_{11}&a_{12}&\ldots&a_{1s}\\ c_{2}&a_{21}&a_{22}&\ldots&a_{2s}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ c_{s}&a_{s1}&a_{s2}&\ldots&a_{ss}\\ \hline\cr&b_{1}&b_{2}&\ldots&b_{s}\end{array}\equiv\begin{array}[]{c|c}\bm{c}&A\\ \hline\cr&\bm{b}^{\intercal}\end{array}.$$

$${\tilde{\mathbf{v}}}_{n,j}={\tilde{\mathbf{v}}}_n+h\sum _{k=1}^s{a}_{jk}{\tilde{U}}_{n+c_k}{\tilde{\mathbf{v}}}_{n,k},j=1,2,\dots ,s,$$

$${\tilde{\mathbf{v}}}_{n+1}={\tilde{\mathbf{v}}}_n+h\sum _{k=1}^s{b}_k{\tilde{U}}_{n+c_k}{\tilde{\mathbf{v}}}_{n,k}.$$

$${\widetilde{\mathbf{{\rm Y}}}}_n={\left({\tilde{\mathbf{v}}}_{n,1},{\tilde{\mathbf{v}}}_{n,2},\dots ,{\tilde{\mathbf{v}}}_{n,s}\right)}^{⊺}\in {\mathbb{C}}^{2s},$$

$${\tilde{D}}_n=\mathrm{diag}\left(h{\tilde{U}}_{n+c_1},\dots ,h{\tilde{U}}_{n+c_s}\right)\in {\mathbb{C}}^{2s\times 2s},$$

$$\left(\begin{array}{cc}{I}_{2s}-(A\otimes {\sigma }_0){\tilde{D}}_n& 0_{2s\times 2}\\ -({\mathbf{b}}^{⊺}\otimes {\sigma }_0){\tilde{D}}_n& {\sigma }_0\end{array}\right)\left(\begin{array}{c}{\widetilde{\mathbf{{\rm Y}}}}_n\\ {\tilde{\mathbf{v}}}_{n+1}\end{array}\right)=1_{s+1}\otimes {\tilde{\mathbf{v}}}_n,$$

$${\Delta }_{n+1}=\det \left(\begin{array}{cc}{I}_{2s}-(A\otimes {\sigma }_0){\tilde{D}}_n& 0_{2s\times 2}\\ -({\mathbf{b}}^{⊺}\otimes {\sigma }_0){\tilde{D}}_n& {\sigma }_0\end{array}\right)=\det \left(I_{2s}-(A\otimes {\sigma }_0){\tilde{D}}_n\right).$$

Full text

††journal: Communications in Nonlinear Science and Numerical Simulation

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. The kronecker product of two matrices $A$ and $B$ is denoted by $A\otimes B$. 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 class of continuous functions is denoted by $\mathsf{C}$ and $\mathsf{C}^{m}$ comprises function that are at least $m$-times differentiable. 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.

1 Introduction

The main objective of this series of papers is to provide a comprehensive survey of the numerical methods available for designing nonlinear Fourier transform (NFT) algorithms with desired stability properties, error behavior and (low) complexity of computation. Our primary interest is in the class of integrators known as exponential integrators which were recently shown to be promising in designing fast NFT algorithms [1, 2] that are stable as well as capable of achieving higher orders of convergence. It must be noted that the family of methods considered in the aforementioned works hardly represent the totality of the methods available for the numerical solution of ordinary differential equations (ODEs). Besides, there are several open problems still to be addressed, for instance, the work presented in [2] does not address the NFTs with periodic boundary condition, or, beyond the second order of convergence, the differential approach has not been shown to be successful in devising inverse NFT algorithms [3, 1, 4, 5]. Therefore, the purpose of this series of papers is to explore the entire landscape of available numerical methods that can yield low complexity algorithms for direct NFT (and perhaps an inverse NFT within the differential approach), and, characterize their utility under various constraints and performance requirements. In this paper, we consider two families of integrators known as the Runge-Kutta method and linear multistep methods et al. [6].

Our exposition begins by collecting some of the basic theoretical results that characterize the signal processing aspects of NFTs. The main results of this paper are contained in Sec. 2 where we present the numerical scheme. It is quite fortunate that the existing literature on the classical methods considered in this work for the solution of the Zakharov-Shabat problem (introduced later in this section) is entirely sufficient for our purpose. In this light, the primary contribution of this paper consists in unifying the approach for fast NFTs and demonstrating its connection with some of the time tested methods for ODEs. These areas of research, however, are still very alive given that the quest for designing efficient Runge-Kutta and related methods is far from complete. For instance, general linear methods which appeared in the work of J. C. Butcher [7, 8] is an evolving novel technique that essentially combines the benefits of Runge-Kutta and linear multistep methods–this would be the subject matter of a forthcoming paper in this series.

1.1 The AKNS system

We begin our discussion with a brief review of the scattering theory closely following the formalism presented in [9]. The nonlinear Fourier transform of any signal is defined via the Zakharov-Shabat (ZS) [10] scattering problem which can be stated as follows: For $\zeta\in\mathbb{R}$ and $\bm{v}=(v_{1},v_{2})^{\intercal}$,

[TABLE]

where $\sigma_{3}$ is one of the Pauli matrices defined in the beginning of this article. The potential $U(t)$ is defined by $U_{11}=U_{22}=0,\,U_{12}=q(t)$ and $U_{21}=r(t)$ with $r=\kappa q^{*}$ ($\kappa\in\{+1,-1\}$). Here, $\zeta\in\mathbb{R}$ is known as the spectral parameter and $q(t)$ is the complex-valued signal. The solution of the scattering problem (1), henceforth referred to as the ZS problem, consists in finding the so called scattering coefficients which are defined through special solutions of (1) known as the Jost solutions which are linearly independent solutions of (1) such that they have a plane-wave like behavior at $+\infty$ or $-\infty$.

• •

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

• •

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

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 relations:

[TABLE]

These coefficients are known as the scattering coefficients and the process of computing them is referred to as direct scattering.

The analytic continuation of the Jost solution with respect to $\zeta$ is possible provided the potential is decays sufficiently fast or has a compact support. If the potential has a compact support, the Jost solutions have analytic continuation in the entire complex plane. Consequently, the scattering coefficients $a(\zeta)$, $b(\zeta)$, $\overline{a}(\zeta)$, $\overline{b}(\zeta)$ are analytic functions of $\zeta\in\mathbb{C}$ [9]. Further, it was shown in [11, 5] that the scattering coefficients are entire functions of exponential type:

Theorem 1.1**.**

Let $q\in\mathsf{L}^{1}$ with support in $\Omega=[-T,T]$ ($T>0$) and set $\mu=\|q\|_{\mathsf{L}^{1}(\Omega)}$. Then the estimates

[TABLE]

where $\tilde{a}(\zeta)$ denotes $[a(\zeta)-1]e^{-2i\zeta T}$, hold.

An immediate consequence of this theorem together with the Paley-Wiener theorem [12, Chap. VI] is that

[TABLE]

where $f(\zeta)$ denotes either $b(\zeta)$ or $\tilde{a}(\zeta)$. The Fourier transformation in the above equation is understood in the sense of distributions.

Let $q(t)\in\mathsf{L}^{1}$ such that $|q(t)|\leq C\exp[-2\delta|t|]$ almost everywhere in $\mathbb{R}$ for some constants $C>0$ and $\delta>0$. Then, the functions $e^{-i\zeta t}\bm{\psi}$ and $e^{i\zeta t}\bm{\phi}$ are analytic in the half-space $\{\zeta\in\mathbb{C}|\,\operatorname{Im}{\zeta}>-\delta\}$. The functions $e^{i\zeta t}\overline{\bm{\psi}}$ and $e^{-i\zeta t}\overline{\bm{\phi}}$ are analytic in the half-space $\{\zeta\in\mathbb{C}|\,\operatorname{Im}{\zeta}<\delta\}$. In this case, the coefficient $a(\zeta)$ is analytic for $\operatorname{Im}{\zeta}>-\delta$ while the coefficient $b(\zeta)$ is analytic in the strip defined by $-\delta<\operatorname{Im}\zeta<\delta$. Further characterization of time-limited and one-sided signals can be found in [5].

Finally, assuming that the analytic continuation is guaranteed, the restriction $r=\kappa q^{*}$ yields

[TABLE]

consequently,

[TABLE]

The scattering coefficients 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 [9, 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(t)$ 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}(t;\zeta_{k})=b_{k}\bm{\psi}(t;\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}$.

In the following, we collect some of the useful results on the characterization of signals synthesized using the inverse NFT when starting from the reflection coefficient where the discrete spectrum is assumed to be empty unless otherwise stated. These results are extremely important in that their characterization remains invariant under evolution consistent with the AKNS-class of integrable nonlinear equations. An interesting result which is due to Van der Mee [13, 14] which prescribes the conditions under which $q$ is square integrable is as follows:

Theorem 1.2** (Square integrable signals [13, 14]).**

Let $\rho(\xi)\in\mathsf{L}^{2}\cap\mathsf{L}^{\infty}$. For $\kappa=+1$, it is additionally assumed that $|\rho|<1$. Then, there exists a unique scattering potential $q\in\mathsf{L}^{2}$.

Another important class of signals we would like to consider are such that their discrete spectrum is empty and the continuous spectrum is compactly supported. In analogy with bandlimited signals in conventional Fourier analysis, such signals are referred to as nonlinearly bandlimited. Such signals lend themselves well to the development of fast inverse NFT algorithm [4]. The characterization of the such signals, in part, is contained in the preceding theorem as a special case. A constructive approach based on sampling expansion presented in [15] yields somewhat stronger result:

Theorem 1.3** (Nonlinearly bandlimited signals [15]).**

Let $\rho(\xi)\in\mathsf{L}^{\infty}$ with $\operatorname*{supp}\rho\subset[-\sigma,\sigma]$. For $\kappa=+1$, it is additionally assumed that $|\rho|<1$. Then, there exists a unique scattering potential $q\in\mathsf{L}^{2}\cap\mathsf{C}^{\infty}$.

Proof.

A constructive proof of this theorem can be obtained using the approach presented in [15] where it is shown that the solution of the Gelfand-Levitan-Marchenko equation exists under the conditioned prescribed in the statement of theorem and $q\in\mathsf{L}^{2}$ by construction. The expression obtained for $q$ is manifestly differentiable to any order. ∎

The work [15] also contains several algorithms for synthesizing the aforementioned class of signals with extremely high accuracy.

1.2 Periodic Boundary Conditions

We would like to qualify our numerical methods for solution of the ZS problem when the potential is periodic; therefore, we briefly introduce the scattering theory for periodic potentials. In this section, we assume that $q(t)$ is a complex-valued periodic signal whose period is $2T$. The solution of the scattering problem proceeds by computing two linearly independent solutions

[TABLE]

such that

[TABLE]

On account of the periodicity of the potential, $\bm{\phi}(t+2T;\zeta)$ is also a solution so that it can be expanded in terms of the linearly independent solutions, $\bm{\phi}$ and $\overline{\bm{\phi}}$, as

[TABLE]

where $a(\zeta;t_{0})$ and $b(\zeta;t_{0})$ are defined to be the scattering coefficients which works to be

[TABLE]

The matrix defined by

[TABLE]

connects the solution at $t=t_{0}$ to $t=t_{0}+2T$, therefore, is referred to as the transfer matrix (or the monodromy matrix). It is easy to see that $\operatorname{\mathscr{W}}(\bm{\phi},\overline{\bm{\phi}})=\kappa\neq 0$ (a consequence of the fact that $\bm{\phi}$ and $\overline{\bm{\phi}}$ are linearly independent) which yields $\det\Phi=1$.

The next important distinction in the periodic case is the existence of the so called Bloch eigenfunctions, denoted by $\bm{\psi}(t;\zeta)$, which is defined by

[TABLE]

Consider the representation

[TABLE]

which translates into

[TABLE]

Therefore, a necessary condition for the existence of Bloch eigenfunctions is that

[TABLE]

The main spectrum of a periodic signal is defined as the set of eigenvalues $\zeta_{i}$ such that $\operatorname{tr}[\Phi](\zeta_{i})=\pm 2$ so that $\lambda=\pm 1$. These values of $\lambda$ correspond to Bloch eigenfunctions that are either periodic or anti-periodic. In the following, we adopt the terminology introduced in Ma and Ablowitz [16] for the characterization of the main spectra. To this end, define $a_{R}$ and $a_{I}(\zeta)$ to be the real and imaginary parts of $a(\zeta;t_{0})$ for $\zeta\in\mathbb{R}$, respectively, so that their analytic continuation into the complex plane reads as

[TABLE]

Note that $a_{R}(\zeta)$ is independent of $t_{0}$ and satisfies the symmetry condition $\overline{a}_{R}(\zeta)=a_{R}(\zeta)$ for all $\zeta\in\mathbb{C}$. The main spectrum can be characterized based on the zeros of $1-a_{R}^{2}(\zeta)$.

2 The Numerical Scheme

In this section, we will develop the integrating factor based exponential Runge-Kutta method for the numerical solution of the Zakharov-Shabat scattering problem. We also include a general discussion of the exponential linear multistep method which was presented in [2] and employed as a benchmarking tool in [4].

Following [1, 2], in order to develop the numerical scheme, we begin with the transformation $\tilde{\bm{v}}=e^{i\sigma_{3}\zeta t}\bm{v}$ so that (1) becomes

[TABLE]

with $\widetilde{U}_{11}=\widetilde{U}_{22}=0,\,\widetilde{U}_{12}=q(t)e^{2i\zeta t}$ and $\widetilde{U}_{21}=r(t)e^{-2i\zeta t}$. Let us remark that the general theory of the Runge-Kutta method and linear multistep methods can be found in classic texts such as Hairer et al. [6] and Butcher [7].

2.1 Runge-Kutta Method

Let the Butcher tableau of a $s$-stage Runge-Kutta method be given by

[TABLE]

Let the step size be $h>0$ and the quantities $c_{j}\in[0,1]$ be ordered so that the nodes within the step can be stated as $t_{n}\leq t_{n}+c_{1}h\leq t_{n}+c_{2}h\leq\ldots\leq t_{n}+c_{s}h\leq t_{n+1}$. For the potential sampled at these nodes, we use the convention $Q_{n+c_{k}}=hq(t_{n}+c_{k}h)$, $R_{n+c_{k}}=hr(t_{n}+c_{k}h)$ and $\widetilde{U}_{n+c_{k}}=\widetilde{U}(t_{n}+c_{k}h)$. In order for the resulting discrete system to be amenable to FFT-based fast polynomial arithmetic, it is sufficient to have each of the $c_{i}$’s to be taken from a set of uniformly distributed nodes in $[0,1]$. Introducing the intermediate stage quantities $\tilde{\bm{v}}_{n,k}$ for $k=1,2,\ldots,s$, we have

[TABLE]

with the final update given by

[TABLE]

Defining the augmented intermediate stage vector

[TABLE]

and introducing the block-diagonal matrix

[TABLE]

we have

[TABLE]

where ‘$\otimes$’ stands for the Kronecker-product of two matrices and $\bm{1}_{s+1}$ denotes a column vector of size $s+1$ with $1$ as its entries. Our goal in the following is to solve the linear system using the Cramer’s rule and find the transfer matrix relationship between $\tilde{\bm{v}}_{n}$ and $\tilde{\bm{v}}_{n+1}$. To this end, consider

[TABLE]

In order to obtain $v_{1}^{(n+1)}$ according to the Cramer’s rule, we must consider the determinant

[TABLE]

where $\bm{0}_{2s}$ denotes a column vector of size $2s$ with [math] as its entries and

[TABLE]

Now, expanding along the second last column, we have

[TABLE]

Next, in order to obtain $v_{2}^{(n+1)}$ according to the Cramer’s rule, we must consider the determinant

[TABLE]

Expanding along the second last column, we have

[TABLE]

Let us introduce the transfer matrix $\widetilde{M}_{n+1}$ so that

[TABLE]

Now, the elements of the transfer matrix can be obtained from the equations (30) and (31) which after some simplification read as follows: The diagonal elements work out to be

[TABLE]

The off diagonal elements are given by

[TABLE]

Consider the relations

[TABLE]

An interesting consequence of this is

[TABLE]

This also turns out to be true of $\det\widetilde{\Gamma}$:

[TABLE]

Consider

[TABLE]

and,

[TABLE]

These relations allow us to conclude:

Lemma 2.1**.**

The transfer matrix $\widetilde{M}_{n+1}$ seen as a function of $\widetilde{D}_{n}$ from a subspace of $\mathbb{C}^{2s\times 2s}$, which satisfies the property (35), to $\mathbb{C}^{2\times 2}$ satisfies the following relations:

[TABLE]

This lemma may prove useful in verifying the output of any computer algebra system that is used for computing the transfer matrix.

Diagonally implicit RK methods are characterized by a lower triangular matrix $A$. In this case $\Delta_{n+1}$ can be computed explicitly:

[TABLE]

For any matrix $A$, let $\|A\|$ stand for the spectral norm unless otherwise stated.

Lemma 2.2**.**

Given samples of the potential, $q_{n+c_{1}},q_{n+c_{2}},\ldots,q_{n+c_{s}}$, assumed to be bounded, there exits a constant $h_{0}>0$ such that

[TABLE]

if $0\leq h\leq h_{0}$.

Proof.

In order to prove the lemma, it suffices to show that for some $h_{0}>0$, the spectral norm of $(A\otimes\sigma_{0})\widetilde{D}_{n}$ is less than unity for all $h\leq h_{0}$ and $\zeta\in\mathbb{R}$. To this end, consider

[TABLE]

Noting that

[TABLE]

it follows that $\|\widetilde{D}_{n}\|\leq\max_{1\leq k\leq s}|Q_{n+c_{k}}|$. If $|q|<\omega$, we may write $\|\widetilde{D}_{n}\|<\omega h$ so that

[TABLE]

The proof follows if we choose $h_{0}<(\omega\|A\|)^{-1}$. ∎

In fact, it is also possible to show that by choosing $h$ sufficiently small for given $\zeta\in\mathbb{C}$, one can ensure that

[TABLE]

We defer this discussion to later part of this section. In order to study the stability of the method in the sense of [17] for fixed $\zeta\in\mathbb{R}$, we write the update relation over one-step as

[TABLE]

then, the boundedness of the matrix

[TABLE]

as $h\rightarrow 0$ is sufficient to guarantee the stability of this one-step method [17]. Observing that

[TABLE]

for all $h\leq h_{0}$ where $\omega$ and $h_{0}$ are as defined in the proof of Lemma 2.2 so that $\|\widetilde{D}_{n}\|<\omega h$. The method (41) is also consistent, i.e., the truncation error goes to [math] as $h\rightarrow 0$, follows from

[TABLE]

provided $(\bm{b}^{\intercal}\mathbf{1}_{s})=\sum_{j=1}^{s}b_{j}=1$. Consistency and stability imply convergence; therefore, the method (41) is convergent under this condition.

Proposition 2.3**.**

Let $q$ be bounded over its support which is assumed to be compact. Then, for fixed $\xi\in\mathbb{R}$, the method defined by (41) is convergent if $\sum_{j=1}^{s}b_{j}=1$.

The transfer matrix relation obtained thus far does not explicitly make it evident that its entries are rationals. To this end, we carry out further simplification so that the rational structure of the transfer matrix becomes self-evident. Let

[TABLE]

and

[TABLE]

so that

[TABLE]

The linear system in (27) can now be written as

[TABLE]

where $\bm{\Upsilon}_{n}=(I_{s}\otimes e^{-i\zeta t_{n}\sigma_{3}})\widetilde{\bm{\Upsilon}}_{n}$. Define

[TABLE]

Then the transfer matrix relation in (32) can be stated as

[TABLE]

where we note that

[TABLE]

and the diagonal elements of the transfer matrix are

[TABLE]

while the off-diagonal elements are given by

[TABLE]

Lemma 2.4**.**

Let $\kappa=-1$. Then, the transfer matrix relation between $\bm{v}_{n}(\zeta)$ and $\bm{v}_{n+1}(\zeta)$ is identical to that between $\sigma_{2}\bm{v}^{*}_{n}(\zeta^{*})$ and $\sigma_{2}\bm{v}^{*}_{n+1}(\zeta^{*})$. In other words, the transfer matrix $M_{n+1}(\zeta)$ in (48) satisfies the relation

[TABLE]

with $\Delta^{*}_{n+1}(\zeta^{*})=\Delta_{n+1}(\zeta)$. Further, for $\xi\in\mathbb{R}$,

[TABLE]

Proof.

For $\kappa=-1$, the symmetry property

[TABLE]

allows us to conclude

[TABLE]

Changing $\zeta$ to $\zeta^{*}$ in (46) and taking the complex conjugate of both sides of the equation, we have

[TABLE]

where $\breve{\bm{\Upsilon}}_{n}=(I_{s}\otimes\sigma_{2})\bm{\Upsilon}^{*}_{n}$. Further,

[TABLE]

shows that the linear system connecting $\sigma_{2}\bm{v}^{*}_{n+1}$ and $\sigma_{2}\bm{v}^{*}_{n}$ with $\zeta^{*}$ as the spectral parameter is identical to that between $\bm{v}_{n+1}$ and $\bm{v}_{n}$ with $\zeta$ as the spectral parameter.

The foregoing conclusion can also be drawn from the structure of the transfer matrix $M_{n+1}(\zeta)$. Consider

[TABLE]

Next,

[TABLE]

which implies

[TABLE]

Now writing ${M}^{(n+1)}_{12}$ in the form

[TABLE]

we have

[TABLE]

∎

Lemma 2.5**.**

Let $\kappa=+1$. Then, the transfer matrix relation between $\bm{v}_{n}(\zeta)$ and $\bm{v}_{n+1}(\zeta)$ is identical to that between $\sigma_{1}\bm{v}^{*}_{n}(\zeta^{*})$ and $\sigma_{1}\bm{v}^{*}_{n+1}(\zeta^{*})$. In other words, the transfer matrix $M_{n+1}(\zeta)$ in (48) satisfies the relation

[TABLE]

with $\Delta^{*}_{n+1}(\zeta^{*})=\Delta_{n+1}(\zeta)$. Further, for $\xi\in\mathbb{R}$,

[TABLE]

Proof.

For $\kappa=+1$, the symmetry property

[TABLE]

allows us to conclude

[TABLE]

The rest of the proof is similar to that of the preceding lemma, therefore, it is being omitted for the sake of brevity of presentation. ∎

Now, we would like to address the stability of the method (41) for fixed $\zeta=\xi+i\eta\in\mathbb{C}$. Rewrite the update relation as

[TABLE]

where

[TABLE]

Now, taking into account that $\|D_{n}\|\leq\omega he^{2|\eta|h}$ where $\omega$ is as defined in the proof of Lemma 2.2 and choosing $h_{0}<(2|\eta|+\omega\|A\|)^{-1}$, we have

[TABLE]

for all $h\leq h_{0}$. Under these conditions, we also have $\Delta_{n+1}(\zeta)\neq 0$. The proposition 2.3 can now be modified as

Proposition 2.6**.**

Let $q$ be bounded over its support which is assumed to be compact. Then, for fixed $\zeta\in\mathbb{C}$ such that $|\operatorname{Im}\zeta|<\infty$, the method defined by (55) is convergent if $\sum_{j=1}^{s}b_{j}=1$.

These observations lead to the conclusion that in order to evaluate the Jost solutions at any point $\zeta$ in the complex plane, the step size must be chosen to be much smaller than that in the case of any point $\xi$ on the real axis. The determination of the discrete spectrum is, therefore, numerically more challenging than that of the continuous spectrum.

Next we look at the recurrence relation for the Wronskian for real values of the spectral parameter. To this end, we set $\zeta=\xi\in\mathbb{R}$. The transfer matrix relation (48) can also be written for the matrix solution of the ZS problem. Let $v_{n}=(\bm{\phi}_{n},\bm{\psi}_{n})$ so that

[TABLE]

The recurrence relation for the Wronskian then works out to be

[TABLE]

Given that $\bm{\phi}_{n}$ and $\bm{\psi}_{n}$ are linearly independent solutions, their Wronskian must not be zero. For sufficiently small $h$, one can ensure that $|M^{(n+1)}_{11}(\xi)|^{2}-\kappa|M^{(n+1)}_{12}(\xi)|^{2}\neq 0$ and $\Delta^{2}_{n+1}(\xi)\neq 0$. For $\kappa=-1$, we have $|M^{(n+1)}_{11}(\xi)|^{2}-\kappa|M^{(n+1)}_{12}(\xi)|^{2}\neq 0$ for all values of $h$.

Finally, let us show that the entries of the transfer matrix are rationals. First, let us observe that

[TABLE]

Let $\nu\in\mathbb{Z}_{+}$ be such that $c_{k}=n_{k}/\nu$ where $n_{k}\in\mathbb{Z}_{+}\cup\{0\}$. Then putting $z=\exp(i\zeta h/\nu)$, we have

[TABLE]

so that $D_{n}=D_{n}(z^{2})$ with finite powers of $z^{2}$. This shows that $M_{n+1}=M_{n+1}(z^{2})$ and $\Delta_{n+1}=\Delta_{n+1}(z^{2})$ with finite powers of $z^{2}$.

2.1.1 Computation of discrete scattering coefficients

Let the computational domain be $\Omega=[T_{1},T_{2}]$ and set

[TABLE]

Let the potential $q$ be supported in $\Omega$ or it is assumed to periodic with period $2T$. Let the number of of segments be $N_{\text{seg.}}$ so that the number of samples of the potential becomes $N=\nu N_{\text{seg.}}+1$ (here $\nu+1$ is the number of distinct nodes within $[0,1]$) and $h=2T/N_{\text{seg.}}$. Consider the Jost solution $\phi(t;\zeta)$:

[TABLE]

Let the cumulative transfer matrix be given by

[TABLE]

then, the discrete scattering coefficients $a$ and $b$ are given by

[TABLE]

defined for $\operatorname{Re}\zeta\in[-\pi/2h,\pi/2h]$ where $h\ell\nu=T$, $h\overline{\ell}\nu=\overline{T}$ and $z=e^{i\zeta h/\nu}$. Also, let $h\ell_{+}=T_{2}$ and $h\ell_{-}=-T_{1}$ so that

[TABLE]

The notation above is consistent with the fact that the transfer matrices have entries that are rational functions of $z^{2}$ .

Let $N_{\text{mat.}}$ be the number of matrices and let $m$ be the number of coefficients of the polynomials involved in the rational entries of each of the matrices. In order to use FFT-based fast polynomial arithmetic [18] to compute the cumulative transfer matrix, $N_{\text{mat.}}$ must be a power of $2$ and the polynomials involved in the rational entries must have number of coefficients $m$ that is a power of $2$. Note that by introducing dummy transfer matrices (which can be the identity matrix) one can easily choose $N_{\text{mat.}}$ to be a power of $2$. The number of coefficients $m$ can also be chosen freely by introducing dummy coefficients which are identically zero. Let $\overline{N}=mN_{\text{mat.}}$. The complexity of computing the polynomial coefficients of the discrete scattering coefficient can be worked out to be $\mathop{\mathscr{O}}\left(\overline{N}\log^{2}\overline{N}\right)$. The evaluation of the scattering coefficient at $N^{\prime}\geq\overline{N}$ (where $N^{\prime}$ is a power of $2$) equispaced points amounts to an FFT operation; therefore, the complexity of computing the reflection coefficient works out to be $\mathop{\mathscr{O}}\left(N^{\prime}\log N^{\prime}\right)$.

In order to represent the Jost solutions, we introduce the polynomial vector

[TABLE]

and the polynomial

[TABLE]

For the methods considered in this paper, the form of the discrete scattering coefficients can be discussed as follows: Let $M_{n}(z^{2})$ have polynomial entries of degree $(z^{2})^{\nu}$ and

[TABLE]

where $\Theta_{n}$ is defined below.

• •

Let $\Delta_{n}(z^{2})z^{-\nu}$ be a polynomial of degree $(z^{2})^{\nu}$ whose constant term is $\Theta_{n}$ so that

[TABLE]

then,

[TABLE]

where we recall $N=\nu N_{\text{seg.}}+1$. By introducing a dummy coefficient, it is also justified to write

[TABLE]

with the domain of definition restricted by $\operatorname{Re}\zeta\in[-\pi/2h,\pi/2h]$.

• •

Let $\Delta_{n}(z^{2})=\Theta_{n}$ which is independent of $z$; then

[TABLE]

The computation of the discrete eigenvalues requires finding all the zeros of $a_{N}(z^{2})$ in $\mathbb{C}_{+}$. We do not address this problem in this paper and refer the reader to the work of Henrici [18] and, more recent, Van Barel et al. [19], Kravanja and Van Barel [20]. We will instead assume that the eigenvalues are known by design111Let us remark that the fastest non-iterative method available for computing the zeros of a polynomial requires $\mathop{\mathscr{O}}\left(N^{2}_{\text{deg.}}\right)$ operations (where $N_{\text{deg.}}$ is the degree of the polynomial) and it is based on finding the eigenvalues of a companion matrix.. The method of computation of the norming constant is adapted from [1] which introduces an additional complexity of $\mathop{\mathscr{O}}\left(KN_{\text{seg.}}\right)$ or, equivalently, $\mathop{\mathscr{O}}\left(KN\right)$ where $K$ is the number of eigenvalues.

Therefore, excluding the cost of computing the discrete eigenvalues, the total complexity of computing the nonlinear Fourier spectrum to be $\mathop{\mathscr{O}}\left(N(K+\log^{2}N)\right)$.

2.1.2 One-stage methods

The simplest example of a one-stage method which achieves second order of convergence is the implicit midpoint (IM) method:

[TABLE]

In the expanded form, we have the system of equations

[TABLE]

Solving this linear system, we have

[TABLE]

Putting $\Theta_{n}=1-\frac{1}{4}Q_{n}R_{n}$ and introducing $z=\exp(i\zeta h)$, the above expression can be written as

[TABLE]

where

[TABLE]

This method resembles the split-Magnus method [1] and it also turns out to be a simplectic method. The layer-peeling formula for this system in terms of the potentials $G$ and $H$ is identical to that of the split-Magnus method. If $H_{n+1/2}$ is known, the potential $R_{n+1/2}$ can recovered as

[TABLE]

It is worth noting that the discrete system based on (exponential) trapezoidal rule presented in [1] can also be transformed into a form similar to (75) (see [21]). This method also belongs to the class of Ablowitz-Ladik problems and admits of a discrete Darboux transformation [21].

By employing the transformation $\bm{w}_{n}=e^{i\zeta\sigma_{3}h/2}\bm{v}_{n}$, we obtain

[TABLE]

so that the scattering coefficients can be represented as

[TABLE]

where

[TABLE]

Remark 2.1** (A modified split-Magnus method).**

Let us note that a modified form of the split-Magnus method can be derived which turns out to be identical to the implicit midpoint method described above. Following [1] and applying Magnus method with one-point Gauss quadrature to the original ZS-problem in (1), we obtain

[TABLE]

Splitting the exponential operator as

[TABLE]

and setting $2\Gamma=\sqrt{Q_{n+\frac{1}{2}}R_{n+\frac{1}{2}}}$, we have

[TABLE]

From here it is evident that the resulting transfer matrix becomes identical to that of the implicit midpoint method.

2.1.3 Two-stage methods

We consider two examples of a two-stage, second order, method which are diagonally implicit. The Butcher tableau for the two-stage Lobatto IIIA and Lobatto IIIB methods are given by

[TABLE]

respectively. Setting $z=\exp(i\zeta h)$, each of these methods yield a transfer matrix relation given by

[TABLE]

where

[TABLE]

while $\Delta_{n+1}=1-\frac{1}{4}Q_{n+1}R_{n+1}$ for the Lobatto IIIA and $\Delta_{n+1}=1-\frac{1}{4}Q_{n}R_{n}$ for Lobatto IIIB method.

2.1.4 Three-stage methods

The Butcher tableau for Kutta’s third order method, which is an explicit method, reads as

[TABLE]

Setting $z=\exp(i\zeta h/2)$, this method simplifies to

[TABLE]

where entries of the transfer matrix are given by

[TABLE]

For testing purposes, we label this method as “RK3”.

The next two examples are implicit methods based on $3$-point Lobatto quadrature: The Butcher tableau for the Lobatto IIIA method of order $4$ reads as

[TABLE]

Again, setting $z=\exp(i\zeta h/2)$, the method simplifies to

[TABLE]

where

[TABLE]

and

[TABLE]

The Lobatto IIIB method of order 4 reads as

[TABLE]

This method simplifies to a form that is identical to (87) with the transfer matrix given by (88) and

[TABLE]

2.1.5 Four stage methods and above

The fourth order classical RK method is given by

[TABLE]

Setting $z=\exp(i\zeta h/2)$, this method simplifies to

[TABLE]

where the entries of the transfer matrix are given by

[TABLE]

and,

[TABLE]

where

[TABLE]

For the sake of completeness, let us discuss a general method to design implicit RK methods with arbitrary number of stages: the collocation method [6]. It is a straight-forward recipe to design implicit RK methods which leads to dense $A$ matrix. Taking the nodes to be $c_{k}=(k-1)/(s-1)$, the coefficients of the RK method are given by

[TABLE]

where

[TABLE]

Clearly $a_{1k}\equiv 0$. The order of convergence of these methods turn out to be $\mathop{\mathscr{O}}\left(h^{s}\right)$ if $s$ is even while $\mathop{\mathscr{O}}\left(h^{s+1}\right)$ if $s$ is odd. An example of a $5$-stage method of order $6$ is

[TABLE]

Note that on account of the dense nature of $A$, the computed transfer matrix (for increasing number of stages) becomes increasingly complicated and cumbersome to implement. A remedy to this situation is to design diagonally implicit RK methods on the uniform grid. It appears that such a program would be more realistically pursued in the framework of general linear methods introduced by Butcher and his school [7, 8]. This method would be taken up in a forthcoming paper of this series.

We conclude this section with the example of a $6$-stage explicit method of order $5$ due to Kutta [7]:

[TABLE]

We label this method as “RK6”. The transfer matrix for this method can be computed easily using any computer algebra system (the expressions are being omitted for the sake of brevity of presentation).

2.2 Linear multistep methods

Fast nonlinear Fourier transform based on exponential linear multistep methods (LMM) were introduced in [2] and used as benchmarking tool in [4]. This family of methods are well suited for vanishing boundary condition; however, periodic boundary conditions would require a Runge-Kutta based starting procedure. For the sake of comparison with the Runge-Kutta approach, we revisit the discretization based on a general zero-stable LMM which can be stated as

[TABLE]

where $\bm{\alpha}=(\alpha_{0},\alpha_{1},\ldots,\alpha_{m})$ and $\bm{\beta}=(\beta_{0},\beta_{1},\ldots,\beta_{m})$. Without loss of generality, we may set $\alpha_{m}\equiv 1$. Solving for $\widetilde{\bm{v}}_{n+m}$, we have

[TABLE]

or, equivalently,

[TABLE]

Let $z=e^{i\zeta h}$ and

[TABLE]

where $\Theta_{n+m}=(1-{\beta}_{m}^{2}R_{n+m}Q_{n+m})$ and

[TABLE]

On simplifying, we have

[TABLE]

which can be put into a one-step form by defining

[TABLE]

so that

[TABLE]

where $\bm{w}_{n}=z^{n}\bm{v}_{n}$ and

[TABLE]

Let us consider the Jost solution $\bm{\phi}(t;\zeta)$. Set the computational domain to be $\Omega=[T_{1},T_{2}]$ with $-h\ell_{-}=T_{1}$ and $h\ell_{+}=T_{2}$. We assume that $q(t)$ is supported in $\Omega$ and $q_{n}=0$ for $n=-m+1,-m+2,\ldots,0$ so that $\bm{\phi}_{n}=z^{\ell_{-}}z^{-n}(1,0)^{\intercal}$ for $n=-m+1,-m+2,\ldots,0$. In order to express the discrete approximation to the Jost solutions, let us define the vector-valued polynomial

[TABLE]

such that $\bm{\phi}_{n}=z^{\ell_{-}}z^{-n}\bm{P}_{n}(z^{2})$. The initial condition works out to be

[TABLE]

yielding the recurrence relation

[TABLE]

where $\boldsymbol{\mathcal{P}}_{n}(z^{2})=(\bm{P}_{n}(z^{2}),\bm{P}_{n-1}(z^{2}),\ldots,\bm{P}_{n-m+1}(z^{2}))^{\intercal}\in\mathbb{C}^{2m}$. From (62), the discrete scattering coefficients work out to be

[TABLE]

which are uniquely defined for $\operatorname{Re}\zeta\in[-{\pi}/{2h},\,{\pi}/{2h}]$. Finally, the overall computational complexity of the direct NFT excluding the cost of computing eigenvalues again works out to be $\mathop{\mathscr{O}}\left(N(K+\log^{2}N)\right)$ [2].

2.2.1 Stability and convergence analysis

The stability and convergence analysis of the multistep method can be carried out by first formulating it as a one-step method. The discussion provided below is taken from [6, Chap. 4, Part III]: Define the augmented vector

[TABLE]

and the diagonal matrix

[TABLE]

Then the method (98) can be stated as

[TABLE]

where

[TABLE]

An LMM is stable if the generating polynomial

[TABLE]

satisfies the root condition, i.e., roots of $\varrho(z)$ lie on or within the unit circle, and, the roots on the unit circle are simple. This is the well-known zero-stability criteria. This property directly implies that $\|A\otimes\sigma_{0}\|\leq 1$ where $\|\cdot\|$ is some matrix norm induced by a vector norm. Putting

[TABLE]

it is easy to conclude that the method is consistent since

[TABLE]

The method (112) now written as

[TABLE]

makes it manifestly evident that it is convergent provided $\|\widetilde{\Lambda}_{n+1}(t_{n};h)\|_{s}$ is bounded which follows from the boundedness of the potential.

2.2.2 Explicit Adams Method

The explicit Adams method (see Table 1) can be stated as

[TABLE]

where $\beta_{m}$ is identically [math] so that $\Theta_{n+m}=1$ and

[TABLE]

2.2.3 Implicit Adams Method

The implicit Adams method (see Table 2) can be stated as

[TABLE]

so that $\Theta_{n+m}=(1-{\beta}_{m}^{2}R_{n+m}Q_{n+m})$ with

[TABLE]

and, for $s<m-1$,

[TABLE]

2.2.4 Backward Differentiation Formula

Discretization using the BDF (see Table 3) can be stated as

[TABLE]

where $\beta_{m}\equiv\beta$ so that $\Theta_{n+m}=(1-{\beta}^{2}R_{n+m}Q_{n+m})$ and

[TABLE]

3 Conclusion

To conclude, we explored in this paper the theoretical aspects of the exponential Runge-Kutta and linear multistep methods for the numerical solution of the Zakharov-Shabat scattering problem. The discrete system obtained can be conveniently stated using a transfer matrix (TM) relation whose explicit form was provided for both of these family of methods. Looking at the structure of the TMs involved, it can be concluded that the Runge-Kutta methods lead to a more optimal TM relation as compared to that of the linear multistep methods; however, the TM for the linear multistep methods is far more easy to derive as well as implement. Both of these family of methods lead to a discrete system that is amenable to FFT-based fast polynomial arithmetic which yields fast algorithms for the computation of discrete scattering coefficients.

Based on the ideas presented in this paper, it is not difficult conclude that (exponential) general linear methods developed by Butcher [7, 8] can also be exploited to obtain a family of fast NFT algorithms. These ideas will be explored in forthcoming paper where we hope to present a comprehensive comparative study of these numerical algorithms.

Appendix A Some useful identities

Cosider a $2\times 2$ matrix $A$, then

[TABLE]

Further,

[TABLE]

and,

[TABLE]

If $A$, $B$, $C$ and $D$ are matrices of such size that one can form the matrix products $AC$ and $BD$, then

[TABLE]

Suppose $A$, $B$, $C$, and $D$ are matrices of dimension $n\times n$, $n\times m$, $m\times n$, and $m\times m$, respectively. If $A$ is invertible, one has

[TABLE]

while, if $D$ is invertible, we have

[TABLE]

Two interesting special cases when $m=1$ are as follows:

[TABLE]

and

[TABLE]