Fast inverse nonlinear Fourier transformation using exponential one-step methods: Darboux transformation

This paper considers the non-Hermitian Zakharov-Shabat (ZS) scattering problem which forms the basis for defining the $\mathrm{SU}\left(2\right)$ nonlinear Fourier transformation (NFT). The theoretical underpinnings of this generalization of the conventional Fourier transformation are quite well established in the Ablowitz-Kaup-Newell-Segur 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 $O\left(KN+N{log}^{2}N\right)$ such that the error in the computed $N$-samples of the $K$-soliton vanishes as $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 classical DT (CDT) algorithm which has a complexity of $O\left(K^{2}N\right)$. We further present an 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 and 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.

Figure 1

The figure shows evolution of the nonlinear Fourier spectrum along the length of the fiber. Here, the sequence $({\zeta }_k,b_k)$ denotes the discrete spectrum and $\rho \left(\xi \right),\xi \in \mathbb{R}$ is the continuous spectrum also known as the reflection coefficient.

Figure 2

The figure shows the schematic of the classical Darboux transformation for a given discrete spectrum, ${\mathfrak{S}}_K$, where the discrete spectrum of the seed solution is empty, i.e., ${\mathfrak{S}}_0=\varnothing$, and the given grid point, $x_n$. The part enclosed within broken lines is referred to as the complete DT block (labeled as $\mathrm{DT}\left[{\mathfrak{S}}_K\right]$). The sole input to this block is the seed Jost solution, $v_0(x_n,t;\zeta )$. The output of the DT block consists of the augmented Jost solution, $v_K(x_n,t;\zeta )$, and the difference between the augmented and the seed potential, $q_K(x_n,t)-q_0(x_n,t)$. Here, $\mathrm{\Delta }{q}_j(x_n,t)=q_j(x_n,t)-q_{j-1}(x_n,t)$ and ${\mathfrak{S}}_j={\mathfrak{S}}_{j-1}\cup \left\{({\zeta }_j,b_j)\right\}$ with ${\mathfrak{S}}_1=\left\{({\zeta }_1,b_1)\right\}\cup {\mathfrak{S}}_0$ where $({\zeta }_j,b_j),j=1,2,...,K$, are the distinct elements of ${\mathfrak{S}}_K$ (see Sec. 2b).

Figure 3

The figure depicts the sequential discrete forward and inverse scattering algorithm in (a) and (b), respectively. The forward scattering algorithm is identical to the transfer matrix approach used to solve wave-propagation problems in dielectric layered media. The inverse scattering algorithm shown here is also known as the layer-peeling algorithm. It consists of using ${\mathit{P}}_{n+1}\left(z^2\right)$ to determine the transfer matrix ${\tilde{M}}_{n+1}\left(z^2\right)$ so that the entire step can be repeated with ${\mathit{P}}_n\left(z^2\right)$ as depicted in (b).

Figure 4

The figure shows the binary-tree structure obtained as a result of applying a divide-and-conquer strategy to the conventional layer-peeling method. The node label depicts the range of indices of the segments or layers ordered from left to right in the computational domain.

Figure 5

The figure depicts the CDT and the FDT algorithm in (a) and (b), respectively, for computing the augmented potential, $q_n^{\left(K\right)}$ [which is the discrete approximation to $q_K\left(x_n\right)]$, by adding a given discrete spectrum, ${\mathfrak{S}}_K$, to the seed potential, $q_n^{\left(0\right)}=q_0\left(x_n\right)$. The DT block, labeled as $\mathrm{DT}\left[{\mathfrak{S}}_K\right]$, is described in Fig. 2. The quantities $v_n^{\left(j\right)}\left(z^2\right)$ refer to the discrete approximation of $v_j(x_n,0;\zeta )$ in Fig. 2.

Figure 6

The figure shows the discrete spectrum, ${\mathfrak{S}}_{16}$ [defined by (141)], where the eigenvalues and the norming constants are shown in (a) and (b), respectively.

Figure 7

The figure depicts the convergence analysis for the norming constants (top row) and the discrete eigenvalues (bottom row). The numerical test (as described in Sec. 7b1) is carried out for a fixed number of eigenvalues ($K\in \{12,16,20\}$).

Figure 8

The figure depicts the convergence analysis (top row) and the run-time behavior (bottom row) for multisolitons. The numerical test (as described in Sec. 7b2) is carried out with a fixed number of eigenvalues ($K\in \{12,16,20\}$).

Figure 9

The figure depicts the relative numerical error (top row) and the run time (bottom row) for multisolitons as a function of number of eigenvalues. The numerical test (as described in Sec. 7b2) is carried out with a fixed number of samples ($N\in \{2^{12},2^{14},2^{16}\}$).

Figure 10

The figure compares the multisoliton potential corresponding to the discrete spectrum, ${\mathfrak{S}}_{16}$ [defined by (141)], generated using the CDT and the FDT algorithm after scaling as discussed in Sec. 7a2.

Figure 11

The figure shows the absolute value of the Lubich coefficients for the Darboux matrix elements corresponding to the discrete spectrum, ${\mathfrak{S}}_{16}$ [defined by (141)]. The underlying one-step method for the Lubich method here is the trapezoidal rule (TR).

Figure 12

The figure depicts convergence analysis (top row) and run-time behavior for the general Darboux transformation. The numerical test is carried out with a soliton-free sech potential as seed and fixed number of eigenvalues ($K\in \{12,16,20\}$) to be added to the profile (see Sec. 7b3).

Figure 13

The figure depicts the relative numerical error (top row) and the run time (bottom row) for the general Darboux transform as a function of number of eigenvalues. The numerical test (as described in Sec. 7b3) is carried out with a fixed number of samples ($N\in \{2^{12},2^{14},2^{16}\}$).