# Direct nonlinear Fourier transform algorithms for the computation of
solitonic spectra in focusing nonlinear Schrödinger equation

A. Vasylchenkova∗ [ J. E. Prilepsky  D. Shepelsky  A. Chattopadhyay  Aston
Institute of Photonic Technologies, Aston University, B4 7ET, Birmingham, UK
System Analytics Research Institute, Aston University, B4 7ET, Birmingham, UK
B. Verkin Institute for Low Temperature Physics and Engineering, Kharkiv
61103, Ukraine  V. N. Karazin Kharkiv National University, Kharkiv 61022,
Ukraine  Mathematics, Aston University, B4 7ET, Birmingham, UK

###### Abstract

Starting from a comparison of some established numerical algorithms for the
computation of the eigenvalues (discrete or solitonic spectrum) of the non-
Hermitian version of the Zakharov-Shabat spectral problem, this article
delivers new algorithms that combine the best features of the existing ones
and thereby allays their relative weaknesses. Our algorithm is modeled within
the remit of the so-called direct nonlinear Fourier transform (NFT) associated
with the focusing nonlinear Schrödinger equation. First, we present the data
for the calibration of existing methods comparing the relative errors
associated with the computation of the continuous NF spectrum. Then each
method is paired with different numerical algorithms for finding zeros of a
complex-valued function to obtain the eigenvalues. Next we describe a new
class of methods based on the contour integrals evaluation for the efficient
search of eigenvalues. After that we introduce a new hybrid method, one of our
main results: the method combines the advances of contour integral approach
and makes use of the iterative algorithms at its second stage for the refined
eigenvalues search. The veracity of our new hybrid algorithm is established by
estimating the convergence speed and accuracy across three independent test
profiles. Along with the development of a new approach for the computation of
the eigenvalues, our study also addresses the problem of computation of the
so-called norming constants associated with the eigenvalues. We show that our
formalism effectively amounts to accurate and fast enough computation of
residues of the reflection coefficient in the upper complex half-plane of the
spectral parameter.

###### keywords:

nonlinear Schrödinger equation , inverse scattering method , numerical
algorithms , signal processing

###### MSC:

[2010] 37K15, 65M12, 35C08

††journal: Commun. Nonlinear Sci. Numer. Simulat.

<EMAIL_ADDRESS>
##  1 Introduction

Nonlinear Schrödinger equation (NLSE) has traditionally played the de facto
subservient modelling role over a wide range of topics that relate to the
dynamical evolution of states and associated variables [1, 2]. Such
applications have been rampant in nonlinear physics, notably in photonics. In
this article, our interest is in arriving at a generic master equation based
formalism, structured around the NLSE, that governs the propagation of the
slow-varying complex optical field envelope $q(z,t)$ along a single-mode
lossless optical fibre [2, 3]:

$i\frac{\partial q(t,z)}{\partial
z}+\frac{1}{2}\,\frac{\partial^{2}q(t,z)}{\partial
t^{2}}+|q(t,z)|^{2}q(t,z)=0.$ (1)

For optical transmission-related problems, $z$ plays the role of the distance
along the fibre while $t$ is the ubiquitous time variable, nomenclature that
we will adhere to in the remainder of this article. It may be worth noting
that similarly named quantities could have entirely different interpretations
in other physical systems. In our study, we will only use suitably normalised
dimensionless variables (both independent and dependent) to enable easy
mapping of our results to other models within the same levels of description.
The NLSE (1) has been explicitly written for the so-called focusing case, a
term that is associated with “anomalous dispersion” in optical fibre studies
[2].

The celebrated work of Zakharov and Shabat [4] revealed that Eq. (1) belongs
to a class of closed form integrable systems that can be completely solved by
the inverse scattering transform method, subject to imposition of appropriate
additional constraints on $q(t,z)$. In optical transmission literature, this
method is popularly referred to as the Nonlinear Fourier Transform (NFT)
method [3], an allusion to the similarity with conventional Fourier transform
that applies to linear PDEs. Such models have been shown to have robust
applications in initial-value problems associated with nonlinear PDEs [5]. In
our study, the NFT operation is landscaped to arrive at a complete set of
nonlinear spectral data (“NF spectrum”) at a given spatial point $z=z_{0}$ by
decomposing the known profile $q(t,z_{0})$; here $q(t,z_{0})$ effectively acts
as the initial condition. The $z$-evolution of the individual NF spectral
components is then decoupled and turns to be linear [1, 3, 4]. To find the
space-time profile at a desired point $z=z_{1}$, we will need to solve the set
of such linearised decoupled equations governing the NF spectrum evolution and
then recover the spectrum distribution at $z_{1}$. $q(t,z_{1})$ as a solution
is uniquely evaluated by the reciprocal (inverse) NFT operation [1, 3, 4, 5,
6, 7].

In this work, our primary aim is to compute the forward NFT assuming a bounded
evolving initial condition, represented as follows:

$\int_{-\infty}^{\infty}\\!|q(t,z_{0})|\,dt<\infty.$ (2)

This condition is automatically satisfied for each example profile studied in
the remainder of this article.

The spectral characteristics of the initial data consists of three parts:
continuous data, defined in terms of the spectral (scattering) functions of
real spectral variable, eigenvalues (distinguished discrete values of spectral
parameter), and respective norming constants associated with the eigenvalues.
Accurate estimation of both parts of discrete NFT data is critically important
as they contain information about the bound states (solitons). Integrability
of Eq. (1) implies that if $q(t,z)$ is a solution of Eq. (1) with prescribed
initial data $q(t,z_{0})$ satisfying Eq. (2), then $q(t,z)$ satisfies Eq. (2)
for all $z$ [8]. Moreover, the spectral data associated with $q(t,z)$ are
independent of $z$, so hereafter, we remove the explicit reference of $q$ in
$z$.

Inverse NFT is conventionally formulated over a subset of potentials $q(t)$,
whose discrete spectrum satisfy the following additional constraints [8]:

  1. (i)

complete lack of any real eigenvalue, and

  2. (ii)

structurally simple eigenvalues.

The second conditions could be technically regarded as a “soft” constraint but
serves as a helpful guide during the numerical evaluation of these numbers.
Also potentials satisfying these assumptions are generic. They form an open
dense set in the appropriate phase space, that relates to a subspace of
(linear) operators originating from this functional space that could be
directly applied to analyze stochastic communication problems, for example, in
studying transmission over a quenched network.

We note here that all three aforementioned components of spectral data
(continuous spectrum, eigenvalues, and norming constants) are interesting from
the perspective of their usage in essentially nonlinearity-free optical
transmission [3, 6]. This is fundamentally important as nonlinearity is known
to play the critical negative role of a “dampener” in high spectral efficiency
optical communication systems, necessitating adequate measures to mitigate its
impact on signal propagation, a key job for communication engineers [9, 10,
11]. This is where NFT-based methods could be highly beneficial, since the
data transmission within the NFT-based framework occurs inside the NFT domain
and thus is not hampered by nonlinearity-related pitfalls (recall that the NFT
modes evolve linearly). While more complex technically, for distortionless
communication chores, NFT-based signal processing could become a highly
efficient alternative for the existing methods [3, 6].

A number of various communication systems based on the modulation of the
different parts of NF spectrum have been proposed and studied recently. In the
eigenvalue communications [12] (name coined after the celebrated work by
Hasegawa and Nyu [13]), the complex soliton eigenvalues are used for
modulation and transmission. The procedure is a natural generalization of the
soliton-based methods [2], the progenitor of NFT-based transmission. The
further extension of this direction involves the utilization of norming
constants together with eigenvalues [6, 14, 15, 16, 17]. Our approach conforms
closely with this methodology; the algorithms and signal processing NFT
methods developed in our article target fast efficient computation of both the
eigenvalues and norming constants. In addition to the eigenvalue based NFT
transmission methods, we also refer to a group of methods that deal with the
continuous part of the NF spectrum and non-solitonic NFT modes, either direct
modulation based [18, 19] or within the so-called digital back-propagation
framework [20, 21]. The computation of the continuous spectral functions is
also addressed in our study. We use them to calibrate our methods. Of late,
there has been a significant progress in the experiments related to the NFT-
based optical transmission [15, 19, 22] that could potentially relate to the
theoretical framework that we propose to develop here. Finally, we note that
Wahls and Poor recently proposed a faster algorithm for the computation of the
NF spectrum [23, 24], where they used a structure similar to the famous fast
Fourier transformation (FFT) architecture [25]. The computational time of this
method grows as $n(\log n)^{2}$, where $n$ is the number of discretization
points (samples). Our current study focuses exclusively on “conventional”
methods, addressing the issue of their improvement, excluding for now the
possibility of adapting our mechanism to undertake similar fast realizations
as in [25].

The paper is organised as follows. Section 2 introduces the mathematical
formalism of the NFT procedure and describes the quality metrics used in our
further analysis. Section 3 introduces several NFT computation methods
together with some of their modifications. The accuracy of the spectral data
are presented as follows:

  1. (i)

as a point-by-point variation of the spectral parameter $\xi$,

  2. (ii)

by using mean squared relative error, and

  3. (iii)

by using the energy associated with the continuous data of the reflection
coefficient $r(\xi)$ (defined later).

Section 4 contains the analysis of iterative and contour integration
algorithms applied for the eigenvalue search, together with analysis of
convergence and stability. A new hybrid approach for stable and accurate
computation of the eigenvalues is detailed in this section. Section 5 deals
with the norming constants encompassing two possible approaches leading to the
calculation of the residue of the associated reflection coefficient and
several computational methods based on these approaches. We show that our new
approach is more accurate and stable (that is, converges faster) compared to
the existing options. Conclusions and future directions are summarized in
Section 6.

##  2 Mathematical formulation, model signals, and performance metrics

###  2.1 Formulation of the direct NFT operation and definition of NF spectrum
quantities

The NFT decomposition of a given pulse $q(t)$ is defined in terms of dedicated
solutions of the Zakharov-Shabat system (ZSS) of ODEs [4]

$\frac{d}{dt}\left(\begin{matrix}u_{1}(t,\xi)\\\
u_{2}(t,\xi)\end{matrix}\right)=\left(\begin{matrix}-i\xi&q(t)\\\
-q^{*}(t)&i\xi\end{matrix}\right)\left(\begin{matrix}u_{1}(t,\xi)\\\
u_{2}(t,\xi)\end{matrix}\right),$ (3)

where $q(t)$, for communication problems, is the signal to process (assuming
the role of the ZSS potential decaying at the asymptotic limits
$t\to\mp\infty$, in accordance with (2)). $\xi$ is the so-called spectral
parameter which can be understood as a nonlinear analogue of frequency.
Asterisk in Eq. (3) and below denotes complex conjugates of corresponding
quantities. Apart from the communication and other NFT applications mentioned
in Sec. 1, the ZSS (3) is interesting by itself as it appears as a master
equation in a coupled mode theory, describing, e.g. the scattering of waves in
the fiber Bragg gratings [26], in electrical circuit related problems [27],
etc.

To retrieve the spectral data associated with the given profile $q(t)$, we fix
the so-called Jost solutions $\Phi$ and $\Psi$ of Eq. (3) imposing the
asymptotic conditions at $t\to\pm\infty$:

$\Phi(t,\xi)\equiv\left(\begin{matrix}\phi_{1}\\\
\phi_{2}\end{matrix}\right)\xrightarrow[t\rightarrow-\infty]{}\left(\begin{matrix}e^{-i\xi
t}\\\
0\end{matrix}\right),\quad\Psi(t,\xi)\equiv\left(\begin{matrix}\psi_{1}\\\
\psi_{2}\end{matrix}\right)\xrightarrow[t\rightarrow\infty]{}\left(\begin{matrix}0\\\
e^{i\xi t}\end{matrix}\right).$ (4)

Using Eq. (4), $\Phi(t,\xi)$ and $\Psi(t,\xi)$ can be equivalently represented
as solutions of Volterra integral equations (see, e.g., [8]), from which it is
seen that they are determined for all $\xi\in\overline{{\mathbb{C}}_{+}}$ (the
upper complex half-plane closed by the real axis), are analytic in
$\mathbb{C}_{+}$ and continuous in $\overline{{\mathbb{C}}_{+}}$. For all
$\xi\in\overline{{\mathbb{C}}_{-}}$, we additionally define
$\hat{\Psi}(t,\xi)=\big{(}\psi_{2}^{*}(t,\xi^{*}),-\psi_{1}^{*}(t,\xi^{*})\big{)}^{T}$.

The goal of the NFT pulse decomposition is to find the continuous and discrete
spectral quantities associated with $q(t)$. Since $\Psi(t,\xi)$ and
$\hat{\Psi}(t,\xi)$, considered for $\xi\in\mathbb{R}$, constitute a
fundamental system of solutions of (3), the representation of $\Phi(t,\xi)$ as
a linear combination of these solutions introduces the spectral (scattering)
coefficients $a(\xi)$ and $b(\xi)$ for $\xi\in\mathbb{R}$ by

$\Phi(t,\xi)=\hat{\Psi}(t,\xi)\>a(\xi)+\Psi(t,\xi)\>b(\xi),\qquad\xi\in{\mathbb{R}}$
(5)

with

$a(\xi)=\det(\Phi(t,\xi),\Psi(t,\xi))\quad\text{and}\quad
b(\xi)=\det(\hat{\Psi}(t,\xi),\Phi(t,\xi)).$ (6)

In view of (4), they can also be expressed by

$a(\xi)=\lim_{t\rightarrow+\infty}\phi_{1}(t,\xi)e^{i\xi t},\qquad
b(\xi)=\lim_{t\rightarrow+\infty}\phi_{2}(t,\xi)e^{-i\xi t}.$ (7)

They can also be characterized in terms of a single function, the reflection
coefficient [8]:

$r(\xi)=b(\xi)/a(\xi).$ (8)

The solitonic degrees of freedom are associated with the discrete spectral
data consisting of the set of complex-valued eigenvalues $\\{\xi_{j}\\}$ of
Eq. (3) that have positive imaginary parts (the set is finite due to condition
(i) given at the end of Sec. 1), together with complex-valued norming
constants $\\{c_{j}\\}$. Since $\Phi(t\to-\infty,\xi)\to 0$ and
$\Psi(t\to\infty,\xi)\to 0$, for any $\xi$ with $\Im\xi>0$, it follows that an
eigenvalue $\xi_{j}$ is characterized by the linear dependence of
$\Phi(t,\xi_{j})$ and $\Psi(t,\xi_{j})$, i. e. by the existence of a non-zero
constant $b_{j}\in\mathbb{C}$ such that

$\Phi(t,\xi_{j})=\Psi(t,\xi_{j})\>b_{j}.$ (9)

Hence the eigenfunction $\bigl{(}u^{(j)}_{1}(t),u^{(j)}_{2}(t)\bigr{)}^{T}$ of
Eq. (3), associated with $\xi_{j}$, is given by
$\bigl{(}u^{(j)}_{1}(t),u^{(j)}_{2}(t)\bigr{)}^{T}=\Phi(t,\xi_{j})=\Psi(t,\xi_{j})\>b_{j}$.

On the other hand, as shown in Eq. (6), the eigenvalues can be equivalently
characterized as zeros of $a(\xi)$ in the upper half-plane: $a(\xi_{j})=0$.
Generally, for potentials satisfying only Eq. (2), the constants ${b_{j}}$ are
independent of the continuous scattering functions, but if $\hat{\Psi}(t,\xi)$
(and thus $b(\xi)$) also admits analytic continuation into the upper half-
plane (or at least in the domain $0<\text{Im}(\xi)<d$, for $d>0$, that is
greater than the imaginary part of all eigenvalues), the parameters
$\\{b_{j}\\}$ represent the values of $b(\xi)$ evaluated at $\xi=\xi_{j}$ (cf.
(5)):

$b_{j}=b(\xi_{j}).$ (10)

In this case, the norming constants can be defined as the residues of the
reflection coefficient $r(\xi)$ at its poles $\\{\xi_{j}\\}$:

$c_{j}=\mathrm{Res}[r(\xi)]|_{\xi=\xi_{j}}=\frac{b(\xi_{j})}{a^{\prime}(\xi_{j})}$
(11)

(where we have assumed that the zeros of $a(\xi)$ are simple. A sufficient
condition ensuring analytic continuation of the above relates to a decay rate
estimate for our profile $q(t)$:

$|q(t)|<De^{-d|t|},\quad\text{for all}~{}t\in(-\infty,\infty),$

for $D,d>0$. Particularly, for finitely supported $q(t)$ (which is the case of
the computational statement of the problem), $b(\xi)$ is analytic in the whole
plane and thus definition (10) holds.

Individual soliton parameters can be directly extracted from the discrete
spectral data [4]: the soliton amplitude is given by $2\text{Im}(\xi_{j})$,
and the soliton frequency is $-2\text{Re}(\xi_{j})$. The norming constant
defines the remaining two soliton parameters: the center position of the
individual soliton,

$j{\text{-th soliton centre
position}}=-\frac{1}{2\text{Im}(\xi_{j})}\log\frac{|c_{j}|}{2\text{Im}(\xi_{j})},$

and the solitonic phase that is proportional to the phase of the norming
constant: $\varphi=-\text{arg}[i\,c_{j}]$. The aforementioned four real
parameters completely characterise each solitonic degree of freedom. Further
details on the ZSS properties and soliton solutions can be found in Refs. [1,
4, 6] and [8].

For our purpose, it will be useful to rewrite ZSS (3) for the wave envelope
functions $\chi_{1,2}$ defined through the relations

$\phi_{1}=\chi_{1}e^{-i\xi t},\qquad\phi_{2}=\chi_{2}e^{i\xi t}.$ (12)

Then the ZSS for the envelope vector
$X(t,\xi)=\big{(}\chi_{1}(t,\xi),\chi_{2}(t,\xi)\big{)}^{T}$ becomes

$\frac{d}{dt}X(t,\xi)=\left(\begin{matrix}0&qe^{2i\xi t}\\\ -q^{*}e^{-2i\xi
t}&0\end{matrix}\right)X(t,\xi).$ (13)

In terms of $X$, the spectral coefficients are given by

$a(\xi)=\lim_{t\to+\infty}\chi_{1}(t,\xi),\qquad
b(\xi)=\lim_{t\to+\infty}\chi_{2}(t,\xi).$ (14)

Since the initial conditions for $X$ in Eq. (13) do not involve exponentials:

$(\chi_{1}(t,\xi),\chi_{2}(t,\xi))^{T}\to(1,0)^{T}\qquad\text{ as}\ \
t\to-\infty,$

the definition of spectral coefficients via $X$ turns to be convenient for
some numerical methods described below.

###  2.2 Model signals for test purposes

Our algorithms will be tested against three independent test profiles. To
ensure maximum possible variation, we have resorted to wide-spread model
signals, where the analytical expressions for the spectral data can be written
explicitly.

  1. (i)

The over-soliton potential [28] (or Satsuma-Yajima pulses) is given by

$q_{\text{over}}(t)=A\operatorname{sech}t.$ (15)

It is characterized by a single real amplitude parameter $A>0$. The associated
spectral functions are as follows (all quantities with suffix “over” represent
oversolitons):

$a_{\text{over}}(\xi)=\frac{\Gamma^{2}\left(\frac{1}{2}-i\xi\right)}{\Gamma\left(-A-i\xi+\frac{1}{2}\right)\Gamma\left(A-i\xi+\frac{1}{2}\right)},$
(16) $b_{\text{over}}(\xi)=-\sin(\pi A)\operatorname{sech}(\pi\xi),$ (17)

and

$r_{\text{over}}(\xi)=-\frac{\sin(\pi
A)\,\text{sech}(\pi\xi)\,\Gamma\left(-A-i\xi+\frac{1}{2}\right)\,\Gamma\left(A-i\xi+\frac{1}{2}\right)}{\Gamma^{2}\left(\frac{1}{2}-i\xi\right)},$
(18)

where $\Gamma(\ldots)$ is the Euler Gamma function.

Depending on the value of $A$, the discrete spectrum attributed to the
oversoliton (15) consists of simple eigenvalues

$\xi_{k}=(A-1/2-k)i,~{}k=0\ldots\left[A-\frac{1}{2}\right],$ (19)

where $[\ldots]$ denotes the integer part. If $A$ is exactly half-integer,
then $r(\xi)=0$ and the total energy is completely concentrated in the
solitonic modes. The norming constant corresponding to the highest eigenvalue
$\xi_{0}=(A-1/2)i$ is

$c_{\text{over}}=i\Gamma(2A)/\Gamma^{2}(A).$ (20)

More details on the NFT properties of the profiles (15) can be found in [28].

  2. (ii)

The ZSS for the rectangular potential

$q(t)=\begin{cases}A,~{}-L\leq t\leq L\\\ 0,~{}\text{otherwise}\end{cases}$
(21)

can also be solved analytically [30, 31]. The associated scattering
coefficients are given by

$a_{\text{rec}}(\xi)=e^{2i\xi
L}\left(\cos\Big{[}2\sqrt{\xi^{2}+A^{2}}\,L\Big{]}-\frac{i\xi}{\sqrt{\xi^{2}+a^{2}}}\sin\Big{[}2\sqrt{\xi^{2}+A^{2}}\,L\Big{]}\right),$
(22)
$b_{\text{rec}}(\xi)=\frac{A}{\sqrt{\xi^{2}+A^{2}}}\sin\Big{[}2\sqrt{\xi^{2}+A^{2}}\,L\Big{]},$
(23)

and

$r_{\text{rec}}(\xi)=\frac{A\exp(-2i\xi
L)}{i\xi-\sqrt{\xi^{2}+A^{2}}\cot\Big{[}2\sqrt{\xi^{2}+A^{2}}\,L\Big{]}}.$
(24)

The discrete eigenvalues $\\{\xi_{rec}\\}$ for the rectangle profile are given
by the roots of the following transcendental equation

$\tan\Big{[}2\sqrt{A^{2}+\xi_{rec}^{2}}\,L\Big{]}=\frac{\sqrt{A^{2}+\xi_{rec}^{2}}}{i\xi_{rec}}$
(25)

for $\xi_{rec}$ in the upper half-plane. The norming constant for $\xi_{rec}$
is given by the expression

$c_{\text{rec}}(\xi)=-\frac{i\left(A^{2}+\xi^{2}\right)e^{-2iL\xi}}{A\left(2L\sqrt{A^{2}+\xi^{2}}\cot\Big{[}2\sqrt{A^{2}+\xi^{2}}\,L\Big{]}-1\right)}.$
(26)

  3. (iii)

For our tests, we have also used the solitonic potential ($r(\xi)=0$,
$\xi\in{\mathbb{R}}$) with unit amplitude and phase [31]:

$q_{\text{sol}}=\exp(-it)\operatorname{sech}(t).$ (27)

It has a single eigenvalue $\xi_{\text{sol}}=0.5+0.5i$ with the associated
norming constant $c_{\text{sol}}=i$. This potential allows us to check the
behavior of our methods in the case of eigenvalues having a non-zero real
part.

###  2.3 Accuracy and performance metrics

Basically, any quantity that can be found analytically for ZSS solution can be
used for the numerical NFT methods’ accuracy assessment. For discrete $y$, we
use relative error as an accuracy descriptor:

$\epsilon=\frac{|y^{(\text{computed})}-y^{(\text{analytical})}|}{|y^{(\text{analytical})}|}.$
(28)

For a continuous spectral function $\varphi(\xi)$ ($\varphi$ can be $a$, $b$
or $r$), we compare the analytical and computed values using the mean squared
relative error (MSRE):

$\epsilon_{\varphi}=\frac{1}{N}\sum_{k=1}^{N}\frac{|\varphi_{k}^{\text{(computed)}}-\varphi_{k}^{\text{(analytical)}}|^{2}}{|\varphi_{k}^{\text{(analytical)}}|^{2}},$
(29)

where $\varphi_{k}=\varphi(\xi_{k})$. For
$\varphi^{\text{(analytical)}}(\xi)=0$, we use the ordinary squared difference
${\big{[}|\varphi_{k}^{\text{(computed)}}-\varphi_{k}^{\text{(analytical)}}|\big{]}}^{2}$.

Following [23, 24], we assess the performance of our algorithms using
normalized runtime $\tau/n$, where $\tau$ measures the computer run time and
$n$ is the number of discretization points (samples) defined over the finite
processing interval $t\in[-L\ldots L]$.

For numerical evaluations, we used four nodes cluster with Intel(R) Xeon(R)
CPU 2.50 GHz.

##  3 Methods and algorithm calibration using continuous spectral functions

In this section, we present several NFT methods and compare their quality and
performance using the aforementioned model potentials over a range of test
parameters. All of our methods use the potential truncation and its
discretization over a finite time interval.

For the truncation of model potentials and discretization, $q(t)$ is
analytically represented to define samples within the interval $t\in[-L,L]$.
The interval is divided into $n$ equal subintervals of length $\Delta t=2L/n$,
where the ${m}$-th subinterval is $t\in[t_{m}-\Delta t/2,t_{m}+\Delta t/2]$.
Outside of the interval $[-L,L]$, the potential is assumed to be exactly zero.
For all methods considered below, our signal (the ZSS potential) is
approximated by a constant value along a single step: $q_{m}=q(t_{m})$.

The vector of wave functions $\Phi(t)=(\phi_{1}(t),\phi_{2}(t))^{T}$ is fixed
by imposing the initial conditions on the left edge of the truncation interval
($t=-L$): according to (4), $\Phi(-L)=(e^{i\xi L},0)^{T}$. The vector of
envelopes $X(t)=(\chi_{1}(t),\chi_{2}(t))^{T}$ is fixed correspondingly as
$X(-L)=(1,0)^{T}$. The spectral functions $a(\xi)$ and $b(\xi)$ are defined on
the right edge of the processing interval ($t=L$): according to (7) and (14),
$a(\xi)=\phi_{1}(L)e^{i\xi L}$ and $b(\xi)=\phi_{2}(L)e^{-i\xi L}$, or,
alternatively, $a(\xi)=\chi_{1}(L)$ and $b(\xi)=\chi_{2}(L)$.

###  3.1 Transfer matrix methods for NF spectrum computation

To propagate the incident wave towards the end of the processing interval,
many approaches use the transfer matrix $T_{m}$ for propagating the ZSS
solution over a single ${m}$-th discretization step, inside which the
potential is considered as a constant, $q_{m}=\text{const}$:

$\Phi_{m+1}=T_{m}\Phi_{m}.$ (30)

Performing the consequent iterations of Eq. (30) from $m=1$ to $m=n$, we
eventually find the desired values of the Jost functions at the end of the
processing interval and compute the NFT parameters $a(\xi)$ and $b(\xi)$.

  1. (i)

For the Bofotta-Osborn (BO) method [3, 30, 31], we define $T_{m}$ evaluating
the matrix exponential of the matrix in ZSS with constant potential $q_{m}$:

$T_{m}^{\text{(BO)}}=\exp\left[\Delta t\left(\begin{matrix}-i\xi&q_{m}\\\
-q_{m}^{*}&i\xi\end{matrix}\right)\right]\\!,$ (31)

which can be evaluated explicitly:

$T_{m}^{\text{(BO)}}=\left(\begin{matrix}\cosh\kappa\Delta
t-i\xi/\kappa\sinh\kappa\Delta t&q_{m}/\kappa\sinh\kappa\Delta t\\\
-q_{m}^{*}/\kappa\sinh\kappa\Delta t&\cosh\kappa\Delta
t+i\xi/\kappa\sinh\kappa\Delta t\end{matrix}\right),$ (32)

with $\kappa=\sqrt{-|q_{m}|^{2}-\xi^{2}}$.

  2. (ii)

For the Ablowitz-Ladik (AL) method, we use the normalized discretization of
ZSS [23, 24, 29] and apply Euler method, then substituting $1\pm i\xi\Delta t$
with $\exp(\pm i\xi\Delta t)$, we get the transfer matrix:

$T_{m}^{\text{(AL)}}=\frac{1}{\sqrt{1+\Delta
t^{2}|q_{m}|^{2}}}\left(\begin{matrix}e^{-i\xi\Delta t}&\Delta tq_{m}\\\
-\Delta tq_{m}^{*}&e^{i\xi\Delta t}\end{matrix}\right).$ (33)

The AL method with the norming factor $1/\sqrt{1+\Delta t^{2}|q_{m}|^{2}}$
provides higher stability and accuracy compared to that without the
normalisation [23].

In this article, we also introduce and study two novel modified versions of BO
and AL algorithms, for which the similar approaches are applied to the ZSS for
the envelope functions (13). Evolution over each step $\Delta t$ is again
performed using the transfer matrices:

$X_{m+1}=T_{m}X_{m}.$ (34)

  1. (iii)

For the modified BO method, applying the matrix exponential to the matrix of
ZSS for the envelope functions defined in Eq. (13), one can get the transfer
matrix in the form:

$T_{m}^{\text{(BOmod)}}=\left(\begin{matrix}\cos|q_{m}\Delta
t|&\sin|q_{m}\Delta t|e^{i(\theta_{q_{m}}+2\xi t)}\\\ -\sin|q_{m}\Delta
t|e^{-i(\theta_{q_{m}}+2\xi t)}&\cos|q_{m}\Delta t|\end{matrix}\right),$ (35)

where $\theta_{q_{m}}$ is $\textrm{arg}[q_{m}]$.

  2. (iv)

In the case of the modified AL method, we have

$T_{m}^{\text{(ALmod)}}=\frac{1}{\sqrt{1+\Delta
t^{2}|q_{m}|^{2}}}\left(\begin{matrix}1&q_{m}\Delta te^{2i\xi t}\\\
-q_{m}^{*}\Delta te^{-2i\xi t}&1\end{matrix}\right).$ (36)

  3. (v)

Finally, for the (non-modified) Crank-Nicolson (CN) method [29, 23], the
transfer matrix entering Eq. (30) is given by

$T_{m}^{\text{(CN)}}=(I-\frac{\Delta t}{2}P_{m+1})^{-1}(I+\frac{\Delta
t}{2}P_{m}),$ (37)

where $I$ is the $2\times 2$ identity matrix and

$P_{m}=\left(\begin{matrix}-i\xi&q_{m}\\\ -q_{m}^{*}&i\xi\end{matrix}\right).$
(38)

Notice that in the NFT related works (see e.g. [29, 31]) the algorithms for
the solutions of ordinary ODEs (the Runge-Kutta scheme) were also studied in
the application of the ZSS analysis. The Runge-Kutta fourth-order algorithm
(RK) for the solution of ZSS is better applicable to the envelope system Eq.
(13), as in this case the rapid oscillations of the Jost functions (for the
region of $\xi$ with a large real part) are included into the effective
potential functions. We do not describe the RK scheme here as it is quite
standard (see e.g. [31]) but below we present a comparative analysis of the RK
algorithm with the algorithms mentioned above.

###  3.2 Results for the continuous spectral functions

Now we perform a comparative accuracy analysis of the described methods using
the three descriptors: $a(\xi)$, $b(\xi)$, and $r(\xi)$, Eqs. (7), (8). For
the real values of $\xi$ the MSRE (29) is used for the accuracy assessment. We
also address the behavior of the NFT methods in analyzing the dependence of
the method’s accuracy on the variation of amplitude: in all the following
figures, the amplitude changes are depicted using the transparency scale of
the corresponding curves, see Fig. 1, i.e. the curves for different amplitudes
are plotted by the colored areas changing from an almost transparent part (the
lowest amplitudes) to an almost solid line (the highest amplitude). Captions
to the plots provide information about the range and step of amplitude
variation.

Figure 1: Transparency scale for amplitude variation in grey tones;
progressively darker shades indicate higher amplitudes

Our analysis confirms that the BO method gives the best accuracy among all
methods studied, see Fig. 2. The AL and the modified BO algorithms display
similar behavior with the change of amplitude $A$ and of the number of points
$n$, whereas the CN algorithm came up with the worst accuracy and convergence
rate. Both AL and BO methods have the same convergence (inclination of the
curves on logarithmic scale plots), implying that they all have the same order
of accuracy. This conclusion complies with the results summarized in [3] and
earlier studies. As expected, the fourth order RK method converges faster.
This method can be better than BO for big enough values of $n$. At the same
time, for a smaller number of dicretization points, the RK method’s MSRE is
excessively large. The BO method shows the weakest dependence on the amplitude
variation, whereas the RK is the most sensitive to it (the error increment can
reach several orders of magnitude in the range of amplitudes that we used for
our plots). We also note that the BO method in application to the rectangular
potential gives the solution, which coincides with explicit analytical
expression, so that the main source of errors here is the computational error
in evaluation of $\cosh$ and $\sinh$ from (32). This offsets the surprising
increase in error with increase in the number of points for this particular
potential and method.

|
---|---
a) $a(\xi)$ error for $q_{\text{over}}$ | b) $b(\xi)$ error for $q_{\text{over}}$
|
c) $a(\xi)$ error for $q_{\text{rec}}$ | d) $b(\xi)$ error for $q_{\text{rec}}$
Figure 2: MSRE for the computation of the NF spectrum parameters (the specific
spectral parameters are given in the captions) as a function of the number of
dicretization points $n$, evaluated for different NFT methods, and shown in
a)–b) panels: for the over-soliton potential $q_{\text{over}}$, where
amplitude changes in the range $[2.25,\,\ldots,\,5.25]$ with the increment
step $\Delta A=0.5$ and $L=30$; c)–d) panels: for the rectangular potential
$q_{\text{rec}}$, where amplitude changes in the range $[2,\,\ldots,\,5]$ with
the increment step $\Delta A=0.5$ and $L=1$.

The energy embedded into the continuous NF spectrum can also be a convenient
quantity for methods’ accuracy assessment: in particular, it was used as a
qualifying metric in [31]. This energy is defined through the spectral
functions $a(\xi)$ or $r(\xi)$ as follows

$\varepsilon=-\frac{1}{\pi}\int_{-\infty}^{\infty}\log|a(\xi)|^{2}d\xi=\frac{1}{\pi}\int_{-\infty}^{\infty}\log(1+|r(\xi)|^{2}).$
(39)

(Recall that $|a(\xi)|^{2}+|b(\xi)|^{2}=1$ for $\xi\in{\mathbb{R}}$).

We analyze the relative error in the calculation of energy (28) versus the
number of points and amplitude, see Fig. 3. The energy analysis shows
qualitatively similar result as the MSRE analysis does, which confirms its
correctness for the NFT methods’ accuracy assessment. The RK algorithm
converges more rapidly than all other methods, especially for the pure
solitonic potential, where the double precision numbers are not enough to find
out the difference between the analytical and computed energy values. At the
same time, this test reveals yet another disadvantage of the RK method: it is
extremely slow in comparison with all other methods (see Table 1). The second
worst in terms of time consumption is the CN method, whereas both the AL and
the modified BO are similar in terms of computational time consumption and are
the fastest among the all discussed methods. These methods are approximately
two times faster than the ordinary BO algorithm, but the latter has an
impressively high accuracy. The results for the amplitude dependence when
using the energy as a metric, Fig. 3, are similar to our findings when the
MSRE for the NFT continuous data was used, Fig. 2.

|  |
---
b) $q_{\text{rec}}$, $L=1$, $A=\pi/2$
a) $q_{\text{over}}$, $L=30$ | c) $q_{\text{sol}}$, $L=30$
Figure 3: Relative error $\epsilon$ (28) in the continuous spectrum energy
$\varepsilon$ (39) versus the number of discretization points $n$ for the
different test potentials. The variation with the potential amplitude change
is shown for over-soliton signal only (left pane), where $A$ changes in the
range $[2.25,\,\ldots,\,5.25]$ with the increment step $\Delta A=0.5$.  Table
1: Runtimes (in seconds) of continuous energy evaluation for different NFT
algorithms for $n=2^{13}$ for different test profiles NFT method
---
BO | BOmod | AL | ALmod | CN | RK
$q_{\text{over}}$, $L=30$, $A=5.25$
125.28 | 71.8 | 70.11 | 71.33 | 232.31 | 1552.99
$q_{\text{rec}}$, $L=1$, $A=\pi/2$
177.43 | 104.21 | 121.03 | 124.99 | 307.55 | 2799.38
$q_{\text{sol}}$, $L=30$
89.99 | 49.52 | 56.14 | 58.32 | 145.67 | 1308.59

In order to analyze the NFT methods stability in dependence on the nonlinear
frequency bandwidth, we investigate the accuracy of our methods along the
nonlinear frequency $\xi$ axis point-by-point, see Fig. 4. As concluded
before, BO and RK methods typically show higher accuracy. In estimating the
accuracy dependence of the frequency, we observe that the numerical error
deviates dramatically for the above-mentioned methods, while for AL, ALmod,
BOmod and CN methods the numerical error is a lot more stable. This effect is
more pronounced for the over-soliton potential. This tendency does not change
significantly when we are tuning the amplitude of the potential. The CN method
demonstrates the worth accuracy among the all methods studied.

|
---|---
a) $a(\xi)$ for $q_{\text{rec}}$ | b) $b(\xi)$ for $q_{\text{rec}}$
|
d) $a(\xi)$ for $q_{\text{over}}$ | e) $b(\xi)$ for $q_{\text{over}}$
Figure 4: Relative error versus $\xi$ (with step size $\Delta\xi=0.1$) and
amplitude, panels a)–b): for $q_{\text{rec}}$ with $L=1$ and amplitude changes
in the range $[2,\,\ldots,\,5]$ with the increment step $\Delta A=0.5$; panels
c)–d): for $q_{\text{over}}$ with $L=20$ and amplitude changes in the range
$[2.25,\,\ldots,\,5.25]$ with the increment step $\Delta A=0.5$. Profiles are
sampled with $n=2^{10}$.

Now let us compare our findings with the previous results. The accuracy
assessment of computing the continuous spectral data was done in [31]. The BO
and RK methods were compared there with regard to the continuous spectrum
energy computation: the convergence of the methods was studied and their
runtimes were analyzed. It was shown that for smooth solitonic potentials, the
RK method was better than the BO method, but the authors [31] attributed this
finding to the properties of the CPU architecture used for their computations.
For the rectangular potential, which has sharp edges, the RK method was shown
to perform slower than the BO one. We have generally observed that in terms of
runtime and resulting accuracy, the BO method typically outperforms the other
options; we have also noticed the excessively large runtime required by the RK
methods in comparison to any transfer-matrix based approach of the type shown
in Eq. (30).. The authors of [23] and [29] presented the detailed description
of various NFT algorithms, but they were mainly focused on the discrete
eigenvalues computation accuracy. The authors of [23] also proposed the fast
implementation of the AL algorithms based on the FFT-type ideology for the
matrix product computation and multipoint fast polynomial evaluation, and this
study was continued in [24] for the periodic NFT potential. The comparison of
the “conventional” non-fast NFT methods accuracy for the case of periodic
potentials and discrete eigenvalues were also presented in [32]. Typically,
the qualitative behavior of the accuracy as a function of $n$ followed the
scenario described above in this section.

##  4 Computation of eigenvalues

The solitonic eigenvalues of ZSS under the assumptions listed above are
defined as zeros of the analytical extension of $a(\xi)$, see (7), into the
upper half-plane $\mathds{C}^{+}$ of the complex $\xi$ plane. In this section
we present two principally different approaches to the computation of location
of solitonic eigenvalues and then compare their accuracy, stability and
performance.

###  4.1 Iterative methods for eigenvalue search

For the eigenvalues computation, the most popular option among the NFT related
works is to apply some iterative scheme to identify the complex zero(s) of a
function (namely, $a(\xi)$ from Eq. (7)) in the case considered [29]. In [24,
29, 31], the dependence of the NFT method performance on a particular
iterative scheme usage was somewhat overlooked and an arbitrarily chosen
method was usually employed without a particular motivation for the choice. In
our paper, we fill this gap and compare the existent iterative schemes: we
analyze the convergence and rapidity of the computational methods in
application to the spectral data computation using the test profiles from
Subs. 2.2.

Traditionally, the most common iterative approaches are the secant and Newton-
Raphson (NR) methods. Assuming an initial guess for the location of the zero
$x_{0}$ of some function, say $f(x)$, the consecutive iteration scheme for the
NR method is given by

$x_{k+1}^{\text{(NR)}}=x_{k}-\frac{f_{k}}{f^{\prime}_{k}},$ (40)

where $f_{k}:=f(x_{k})$. This method has a quadratic convergence rate [39]
(numerical method is said to have convergence rate $p$, if $|x_{k+1}-x|\leq
C|x_{k}-x|^{p}$, here and below the orders of convergence are given under the
assumption that all roots are simple). The main disadvantage of the NR method
is the necessity to know the value of the function derivative at each
iteration step. For the purpose of brevity, we introduce a shorthand notation
for the so-called divided differences:

$f[x_{1},x_{2}]=\frac{f(x_{2})-f(x_{1})}{x_{2}-x_{1}},$ (41)

which can be recursively generalized for an arbitrary number of arguments:

$f[x_{k},\ldots,x_{k+m}]=\frac{f[x_{k+m},\ldots,x_{k+1}]-f[x_{k+m-1},\ldots,x_{k}]}{x_{k+m}-x_{k}}.$
(42)

In the secant method, the expression for the derivative in each iteration is
swapped over to the divided difference, leading to:

$x_{k+1}^{\text{(secant)}}=x_{k}-\frac{f_{k}}{f[x_{k-1},x_{k}]}=x_{k}-f_{k}\frac{x_{k}-x_{k-1}}{f_{k}-f_{k-1}}.$
(43)

The convergence rate of the secant method is approximately 1.618, so it is
worse than that for the NR method, but a single step computation using the
secant method can be faster since it does not require computing the
derivatives. Sidi [33] generalised the idea of the derivative approximation:
the function derivative $f^{\prime}_{k}$ can be replaced by the derivative of
a fitting polynomial $p(x)$ of degree $j$:

$p^{\prime}_{j}(x)=f[x_{k},x_{k-1}]+\sum_{i=2}^{j}f[x_{k},...x_{k-i}]\prod_{l=1}^{i-1}(x_{k}-x_{k-l}).$
(44)

The next iteration is given by

$x_{k+1}^{\text{(Sidi)}}=x_{k}-\frac{f_{k}}{p^{\prime}_{j}}.$ (45)

In our study, we use a cubic polynomial approximation in (44), i.e. $j=3$. For
this particular case, the convergence rate of the method is $\approx 1.93$
[33].

Steffensen’s method [34] uses the following iterative formula:

$x_{k+1}^{\text{(Steffensen)}}=x_{k}-\frac{f_{k}^{2}}{f(x_{k}+f_{k})-f_{k}}.$
(46)

It allows us to reach the convergence rate $2$, same as that for the NR
method.

The Muller method [35] has an advantage in that it allows us to find complex
roots from a real initial guess. Defining the auxiliary quantities

$\displaystyle w=f[x_{k-1},x_{k}]+f[x_{k-2},x_{k}]-f[x_{k-2},x_{k-1}],\qquad
g=f[x_{k-2},x_{k-1},x_{k}],$ (47) $\displaystyle
d=\text{max}\left[w-\sqrt{w^{2}-4f_{k}g},~{}w+\sqrt{w^{2}-4f_{k}g}\right],$

(the maximum is determined by comparing the absolute values), the iteration
step of the Muller method is given by

$x_{k+1}^{\text{(Muller)}}=x_{k}-2\cdot f_{k}/d.$ (48)

The order of convergence for this method is approximately $1.84$ [39], which
is better than that for the secant method.

All iterative algorithms applied for the eigenvalues computation require good
initial guess. In order to understand how the choice of the initial value
influences the result of the eigenvalues search, we investigate the
convergence of all iterative schemes in dependence on the initial guess point
value: Fig. 5 contains the results referring to the solitonic potential with
phase $q_{\text{sol}}$, and in Fig. 6 we depict the results for the
rectangular potential $q_{\text{rec}}$ with $A=\pi/2$ where only a single
eigenvalue is present. On these plots we show the border of the regions in the
complex plane of spectral parameter $\xi$ (marked with the closed lines of
different colour), where the relative error of zero location (estimated by
using (28)) is less than $0.01$. We run the iterative algorithm until it
reaches the pre-set precision in the difference between the function values
for consequent iterations or until it exceeds the pre-set number of iterations
(these pre-sets were correspondingly $10^{-10}$ and $10^{3}$). We also add to
each line on the plots the corresponding average runtime that the computation
of the eigenvalue takes when the initial guess point is positioned inside the
respective regions.

|  |
---|---|---
a) Newton-Raphson | b) secant | c) Sidi

|
---|---
d) Muller | f) Steffensen

Figure 5: The borders of the regions of initial assumptions for the zero
approximation, from which the corresponding iterative algorithm reaches its
correct value (marked as black point) with less than 1% relative error.
Computation was performed for the soliton with phase factor potential
$q_{\text{sol}}$ given by Eq. (27) for the values $n=2^{10},~{}L=30$. The
digits in the legends near each curve identify the average runtime for the
computation of an eigenvalue, when the initial guess point was taken inside
the respective regions. Perfect vertical and horizontal edges of some basins
mean that regions of convergence extend out of ranges where we make these
probes.

|  |
---|---|---
a) Newton-Raphson | b) secant | c) Sidi

|
---|---
d) Muller | f) Steffensen

Figure 6: The borders of the regions of initial guesses for zero
approximation, from which the corresponding iterative algorithm reaches its
best value (marked as black point) with less than 1% relative error.
Computation was performed for rectangular potential $q_{\text{rec}}$
($n=2^{10},~{}L=1,~{}A=\pi/2$). The digits in the legends near each particular
line identify the average runtime for the computation of an eigenvalue, when
the initial guess point was taken inside the respective regions. Perfect
vertical and horisontal edges of some basins mean that region of convergence
extends out of ranges where we make these probes.

The result of our analysis of the aforementioned five iterative methods
combined with different NFT algorithms can be summarized as follows: the NR
and secant methods, used in [24, 29], have smaller regions of convergence than
those of the Sidi and Muller methods (the Muller method shows the largest
convergence basin for all example profiles used), whereas the runtimes for all
four approaches are similar. In the case of the Sidi algorithm, we observe a
better convergence, but with significantly higher runtime, especially when the
iterative method is coupled with the AL and CN NFT algorithms. The Steffensen
iterative method shows even smaller region of convergence, which is worse than
that for the NR and secant methods.

The drawback of iterative methods in the application to real transmission
problems (when the position of solitonic eigenvalues is usually not known
apriori) is that we cannot predict the computational runtime since we cannot
estimate the number of iterations required to find the zero point with a
satisfactory accuracy. Moreover, the methods can fail to converge at all, so
that some additional precautions have to be taken. The previous works devoted
to the eigenvalue search algorithms assumed iterative schemes for estimating
discrete eigenvalues. The main results referred to particular features of the
NFT method but not necessarily to a particular iterative scheme. Authors of
[29] studied the one-soliton and multisoliton cases for the AL, Euler and RK
methods (also using the CN method for some cases, but evidently revealing the
weakness of the CN method). They also found out that all studied root-search
techniques resulted in similar accuracy regimes for the NFT data computed.
According to [31], the RK algorithm can converge faster than the BO, but the
authors there used the grid search for the location of eigenvalues, and this
resulted in a relatively high overall runtime. It was found that in the case
of rectangular potential, the RK method failed in the correct computation for
the total number of zeros (the authors used the total phase increment along
the $\text{Re}(\xi)$ axis for this purpose). Our analysis of the iterative
algorithms reveals that they are not sufficiently stable and manageable for
the eigenvalue computations in realistic applications. This fact motivates us
to seek for principally different methods and options applicable for the
location of eigenvalues.

What we have not studied here, though, are the class of methods that involve
matrix diagonalization, e.g. the Fourier collocation method [29]. The idea of
the Fourier collocation is the decomposition of the ZSS problem in the Fourier
series and the subsequent diagonalization of a specific complex-valued non-
Hermitian block matrix, recasting the zero search as an eigenvalue problem.
According to [29], this method can provide a good accuracy (its accuracy is
spectral, in contrast to the methods studied in our current paper), but it has
so far been adapted to find the eigenvalues only. Usually it also requires a
considerable runtime ($\sim n^{3}$) to reach the result.

###  4.2 Contour integration methods for eigenvalues search

In this subsection, we present a new alternative technique for eigenvalue
computation. Delver and Lyness in [37] presented an approach for the location
of the complex zeros of nonlinear functions based on the contour integral
evaluation, see also [38] for more resent results on this approach. This
method will be identified as DL (or ordinary DL) in the rest of this text.
Within the method of [37], we start from the relation:

$\frac{1}{2\pi i}\int_{C}z^{p}\frac{f^{\prime}(z)}{f(z)}dz=\sum z_{i}^{p},$
(49)

that allows us to make up the root-search scheme that would be applicable for
finding all zeros $z_{i}$ of a function $f(z)$ inside the closed contour $C$
drawn in the complex plane. Within this approach, the zeros $z_{i}$ will
emerge as the roots of some specially constructed polynomial. First, setting
$p=0$ in (49) gives us the total number of zeros, $N$, located inside the
contour $C$. Next, evaluating (49) for $p=1\ldots N$ we can readily find the
sums of $z_{i}^{p}$ up to $p=N$:

$s_{p}=\sum z_{i}^{p},\qquad p=1\ldots N.$ (50)

Having found the set of $s_{p}$, we can write down the equation system to
evaluate the so-called Newton’s identities, $\sigma_{p}$:

$\displaystyle\sigma_{1}$ $\displaystyle=-\sum z_{i}:~{}s_{1}+\sigma_{1}=0$
(51) $\displaystyle\sigma_{2}$
$\displaystyle=z_{1}z_{2}+z_{2}z_{3}+...+z_{n-1}z_{N}:~{}s_{2}+s_{1}\sigma_{1}+2\sigma_{2}=0$
$\displaystyle...$ $\displaystyle\sigma_{N}$
$\displaystyle=(-1)^{N}z_{1}z_{2}...z_{N}:~{}s_{N}+s_{N-1}\sigma_{1}+...+s_{1}\sigma_{N-1}+N\sigma_{N}=0.$

This system of equations can be solved recurrently using the values of the set
{$\sigma_{i}$} obtained in the previous iteration round:

$\sigma_{p}=\frac{1}{p}\left(\sum_{j=1}^{p-1}s_{j}\sigma_{p-j}+s_{p}\right).$
(52)

Next, using the Newton’s identities we can construct the following polynomial:

$P(z)=z^{N}+\sigma_{1}z^{N-1}+\sigma_{2}z^{N-2}+\ldots+\sigma_{N-1}z+\sigma_{N}.$
(53)

The polynomial $P(z)$ has exactly the same roots as the initial function
$f(z)$. Therefore, using any polynomial root-finding technique, the desired
set of roots of function $f(z)$, $z_{i}$, can be estimated. If the initial
function has multiple roots, then they will be presented repeatedly along with
the set of polynomial roots.

We remark that the integrand in (49) can be approximated using the discrete
difference in place of the derivative term:

$z^{p}\frac{f^{\prime}(z)}{f(z)}dz\rightarrow
z_{k}^{p}\frac{f^{\prime}(z_{k})\Delta z}{f(z_{k})}\approx
z_{k}^{p}\frac{f(z_{k})-f(z_{k-1})}{f(z_{k})}=z_{k}^{p}\left(1-\frac{f(z_{k-1})}{f(z_{k})}\right).$
(54)

We will refer to the contour integration method used for the search of
solitonic eigenvalues with integrand approximated as in (54) (i.e. without
explicit derivatives) as aDL.

More recently, Kravanja et al in [38] presented an improved version of the DL
approach. Their algorithm relied on a recursive construction of the so-called
formal orthogonal polynomials, which have the roots that coincide with the
zeros of our function $f(z)$. We have analyzed the accuracy of this improved
method for eigenvalues computation too, but we do not present these results
here on separate plots because our study did not reveal any noticeable
difference between results obtained from this newer method compared to those
of the ordinary DL method (at least, for the set of our test profiles).
However, we note that for some real-world applications, where signals are not
smooth and often significantly corrupted by noise [3], the approach proposed
in [38] might demonstrate a better performance.

Since all zeros inside the given contour can be located simultaneously
together with their multiplicity, the contour integration method can guarantee
more stability of the overall algorithm, especially when multiple eigenvalues
are to be found. To ensure that we successfully localize all desired zeros, a
large enough contour needs to be defined in the $\xi$ half-plane. In addition,
the contour integration method runtime depends only insignificantly on the
number of eigenvalues insofar as the values of logarithmic derivatives along
the contour from (49) can be computed just once and then saved for further
processing. On the contrary, iterative algorithms, Subs. 4.1, can
uncontrollably scan the complex plane and arrive at a zero value in the lower
half-plane of $\xi$. In the case of several solitons, iterative algorithms
evaluate all zeros independently, and that fact can increase the computational
time. However, we note that for some well-defined cases, the contour
integrals’ computation can turn out to be more time-consuming in comparison
with the iterative algorithms, because the former requires performing the
computation of the value of $a(\xi)$ (possibly $a^{\prime}(\xi)$ as well) for
the whole set of discrete points along the contour $C$, while for a good
initial guess the iterations may only require several repetitive evaluations
of the function and its derivative.

We analyzed the dependence of contour integration algorithms’ accuracy on the
number of points along the contour and on the contour shape, Fig. 7. We study
the behavior of the aDL method with the approximation scheme represented in
Eq. (54) and compare it with the ordinary DL method’s behavior. As it was
expected, the latter works more accurately due to a more accurate calculation
of the derivative, see the blue and orange lines and compare them with red and
green ones in Fig. 7. We also observed that the particular contour shape
chosen in Eq. (49) also influences the resulting accuracy of the eigenvalues
found. In particular, we checked the behavior of the methods using the
rectangular contour in the upper half-plane of $\xi$, fixing the contour
borders along $\text{Re}(\xi)$-axis and $\text{Im}(\xi)$-axis. Another option
that we tested was to define the ring sector in the $\xi$-plane, fixing the
borders for absolute value $\rho$ and for argument $\theta$ of $\xi$ written
in polar representation as $\xi=\rho e^{i\theta}$.

Figure 7: The dependence of the relative error on the number of contour
discretization points for the rectangular potential ($L=1$, $A=\pi/2$ implies
a single eigenvalue) and for the soliton potential (27) with phase factor
($L=20$), for the rectangular integration contour: $\text{Re}(\xi)$ changes in
the range $[-1\ldots 1]$, $\text{Im}(\xi)$ changes in the range $[0.1\ldots
2]$; and ring sector shape of the contour: $\xi=\rho e^{i\theta}$, $\rho$
changes in the range $[0.1\ldots 2]$, $\theta$ changes in the range
$[\pi/12\ldots 11\pi/12]$.

The integrals in all of our methods have been evaluated using the trapezoidal
rule. Since the runtime of this algorithm depends linearly on the number of
points along the contour, we found sufficiently enough to compare the runtimes
of all approaches for the largest number of points (see Table 2). We found
that results from the ordinary DL method and its approximated aDL version
differed: the runtime of the ordinary DL method is typically 1.5-2 times
larger than that of the aDL one, but the ordinary DL method allows us to reach
a smaller relative error. From Fig. 7, we can also readily see that the
rectangular contour gives a better accuracy, at least for the test profiles
that we studied. In our tests the computation time of the contour approach was
typically higher than that for iterative algorithms from Subs. 4.1 (cf. the
runtimes for different methods given in Figs. 6, 5 and 7). This result is well
explainable, since we were taking the guess points that were close enough to
the zero point, such that the iterative algorithm was able to reach a zero in
just several iterative steps, ensuring a lower computational time.

Table 2: Runtimes (in seconds) of discrete eigenvalues’ evaluation for different test profiles, contour integral approaches and shapes of contours: the rectangular integration contour: $\text{Re}(\xi)$ changes in the range $[-1\ldots 1]$, $\text{Im}(\xi)$ changes in the range $[0.1\ldots 2]$; and ring sector shape of the contour: $\xi=\rho e^{i\theta}$, $\rho$ changes in the range $[0.1\ldots 2]$, $\theta$ changes in the range $[\pi/12\ldots 11\pi/12]$ for 1600 points along the integration contour. rectangle DL | rectangle aDL | sector DL | sector aDL
---|---|---|---
$q_{\text{sol}}$
88.27 | 51.49 | 88.27 | 51.37
$q_{\text{rec}}$, $L=1$, $A=\pi/2$
109.65 | 69.75 | 109.65 | 69.64

###  4.3 Derivative computation

Both the contour integration and the iterative algorithms (in particular, the
NR method) require us to find the value of $a^{\prime}(\xi)$ together with the
value of $a(\xi)$ at the same point $\xi$. Since $a(\xi)$ is holomorphic in
the upper half-plane of $\xi$ (see, for example, [1]), the divided difference
can be used to approximate the derivative. Hoverer, it is also possible to
find the value of the derivative more accurately in the same programming loop
together with the function computation itself. A similar approach is described
in [17, 29, 30, 31]. The idea here is to evolve the derivative of
$\Phi(t,\xi)$ along $t$ together with $\Phi(t,\xi)$ itself. In the case of
ordinary ZSS (3), the evolution starts with the “initial condition” for the
derivative defined at $t\to-\infty$:
$\Phi^{\prime}(t\to-\infty,\xi)\to(-ite^{-i\xi t},0)^{T}$. After the
truncation of the $t$-interval, the initial condition takes the form
$\Phi^{\prime}(-L,\xi)=(iLe^{i\xi L},0)^{T}$. For the ZSS written for the
envelope function (13), we have $X^{\prime}(t\to-\infty,\xi)=(0,0)^{T}$, or,
after the truncation, $X^{\prime}(-L,\xi)=(0,0)^{T}$. Evolution of the wave
function derivative over a single step is performed by applying the relation:

$\Phi^{\prime}_{m+1}=T_{m}^{\prime}\Phi_{m}+T_{m}\Phi^{\prime}_{m},\quad\text{or}\quad
X^{\prime}_{m+1}=T_{m}^{\prime}X_{m}+T_{m}X^{\prime}_{m}.$ (55)

Thus, for all transfer matrix methods from Subs. 3.1, the matrix
$T^{\prime}_{m}$ can be easily determined. At the same time, the RK method
does not have that advantage: here we have to numerically solve the ZSS for
the derivatives separately, equipped with different initial conditions.
However, we recall that the RK method is the most time-consuming according to
our study, see Subs. 3.2, and so we do not deal with it further.

Our algorithm requires the use of the derivatives of the transfer matrices of
each NFT algorithm described in Subs. 3.1. In the case of the BO matrix, the
calculation of the derivative yields

$T_{m}^{{}^{\prime}\text{(BO)}}=\left(\begin{matrix}A_{+}&B\\\
C&A_{-}\end{matrix}\right),$ (56)

where

$\displaystyle A_{\pm}=\pm\frac{i\xi^{2}\Delta t}{\kappa^{2}}\cosh\kappa\Delta
t\mp\sinh\kappa\Delta t\left(\pm\frac{\xi\Delta
t}{\kappa}+\frac{i}{\kappa}+\frac{i\xi^{2}}{\kappa^{3}}\right),$ (57)
$\displaystyle B=\frac{q_{n}\xi}{\kappa^{3}}\sinh\kappa\Delta
t-\frac{q_{n}\xi\Delta t}{\kappa^{2}}\cosh\kappa\Delta t,$ $\displaystyle
C=-\frac{q_{n}^{*}\xi}{\kappa^{3}}\sinh\kappa\Delta t+\frac{q_{n}^{*}\xi\Delta
t}{\kappa^{2}}\cosh\kappa\Delta t,$

where we keep the notations from Eq. (32):
$\kappa=\sqrt{-|q_{m}|^{2}-\xi^{2}}$. For the AL method, the derivative of the
transfer matrix acquires the following form:

$T_{m}^{{}^{\prime}\text{(AL)}}=\frac{1}{\sqrt{1+\Delta
t^{2}q_{m}^{2}}}\left(\begin{matrix}-i\Delta te^{-i\xi\Delta t}&0\\\ 0&i\Delta
te^{i\xi\Delta t}\end{matrix}\right).$ (58)

For the modified BO method we have

$T_{m}^{{}^{\prime}\text{(BOmod)}}=2it\sin|q_{m}\Delta
t|\left(\begin{matrix}0&e^{i(\theta_{q_{m}}+2\xi t)}\\\
e^{-i(\theta_{q_{m}}+2\xi t)}&0\end{matrix}\right),$ (59)

where $\theta_{q_{m}}$ is an argument of $q_{m}$. For the modified AL matrix
we have:

$T_{m}^{{}^{\prime}\text{(ALmod)}}=\frac{2it\Delta t}{\sqrt{1+\Delta
t^{2}|q_{m}|^{2}}}\left(\begin{matrix}0&q_{m}e^{2i\xi t}\\\ q_{m}^{*}e^{-2i\xi
t}&0\end{matrix}\right).$ (60)

The CN-method transfer matrix derivative has the form [29]:

$T_{m}^{{}^{\prime}\text{(CN)}}=\frac{1}{2}P_{m+1}(I-\frac{\Delta
t}{2}P_{m+1})^{-2}(I+\frac{\Delta t}{2}P_{m})+\frac{1}{2}(I-\frac{\Delta
t}{2}P_{m+1})^{-1}P_{m},$ (61)

where we again keep the notations from (38). In order to find
$a^{\prime}(\xi)$, an additional step must be taken at the end of the
algorithm. For the ordinary ZSS (3), we have the following expression for the
derivative of $a(\xi)$ that involves the elements of the Jost solution and its
derivatives:
$a^{\prime}(\xi)=(\phi^{\prime}_{1}(L,\xi)+iL\phi_{1}(L,\xi))e^{i\xi L}$; in
the case of the envelope ZSS (13), we arrive at:
$a^{\prime}(\xi)=\chi^{\prime}_{1}(L,\xi)$.

To check the accuracy of the derivatives evaluated via different algorithms,
we resort to similar mechanisms as described above. We primarily use the
scattering coefficient $a(\xi)$ for the rectangular (22) and over-soliton (16)
potentials. In Figs. 8, we present the dependence of the relative error for
the derivative $a^{\prime}(\xi)$ on $\xi$ and on the number of subintervals
$n$ for the different methods described above. Our computations confirms that
the BO method is again the most accurate one, while the accuracy of the CN is
the lowest among all methods studied.

|
---|---
a) | b)
|
c) | d)
Figure 8: Relative error for $a^{\prime}(\xi)$ computation versus $\xi$ (with
step size $\Delta\xi=0.1$, $n=2^{10}$) and $n$ using different NFT algorithms
from Subs. 4.3, panels a)–b): for the rectangular potential with $L=1$,
amplitude changes in the range $[2,\,\ldots,\,5]$ with the step $\Delta
A=0.5$; panels c)–d): for the over-soliton potential with $L=20$, amplitude
changes in the range $[2.25,\,\ldots,\,5.25]$ with the step of the increment
$\Delta A=0.5$.

###  4.4 Computation of the total number of eigenvalues

A serious difficulty in dealing with iterative algorithms is that they require
an additional adaptation for the multisolitonic case. By definition, the
iterative algorithms seek for only one zero of $f(x)$, which is usually the
closest one with respect to the initial guess point. If the input signal
comprises of more than one eigenvalue (e.g. when we have sufficiently large
amplitudes for the rectangle or over-soliton potentials), it leads to the
following requirements for iterative algorithms:

  * 1.

The total number of eigenvalues must be known before the search routine is
executed. The number of eigenvalues can either be known from the system
properties or, alternatively, must be computed.

  * 2.

When a solution (zero value of the function) is found, it must be eliminated
from the next search runs, to exclude the possibility of repetition.

In order to find the total number of zeros, we can use the following relation
(logarithmic derivative):

$N=\frac{1}{2\pi
i}\int_{C}\frac{f^{\prime}(x)}{f(x)}dx=\frac{1}{2\pi}\Delta_{C}\text{arg}f(x).$
(62)

In Eq. (62) $\Delta_{C}\text{arg}f(x)$ represents the incremental change in
the value of the argument of a complex-valued function $f(x)$, when $x$
traverses the closed contour $C$ in the complex plane. We can manually choose
a sufficiently large $C$ in order to keep all possible zeros inside the
contour. Alternatively, we can calculate phase increment along the real axis,
assuming that the contour contains a part of the real axis and that the other
part of the contour located in the upper half-plane makes negligible
contribution to the total increment value. We tested the applicability of all
described ways for the computation of the number $N$ of zeros and found out
(see Fig. 9) that the precise logarithmic derivative integration was the
slowest option for determining $N$.

Figure 9:  Effective “number of zeros”, obtained via the computation of phase
increment, versus runtime diagram for different number of zeros evaluation
techniques (the BO NFT computation method was used for the illustration) at
the following profiles: the soliton with phase factor $q_{\text{sol}}$, over-
soliton $q_{\text{over}}$ and rectangular potentials, both with $A=2$, $2.5$
and $3$.

Then, the computation of the phase increment and the approximate integration
using Eq. (54) as integrand in Eq. (62), both work faster and produce a
smaller deviation from an integer number; the latter are displayed in Fig. 9
as thick vertical lines. To plot this diagram, we check the number of zeros
for all presented model signals from Subs. 2.2 with randomly chosen
amplitudes. We notice that for the computation of the total number of
eigenvalues, $N$, it is not necessary to compute the expression (62) very
accurately; the numerical error just needs to be small enough to distinguish
between two successive integer values of $N$.

In order to eliminate an estimated zero, $x_{i}$, from the next search run in
the iterative schemes, the function $f(x)$ has to be redefined as

$f(x)\rightarrow\frac{f(x)}{x-x_{i}}.$ (63)

We note that after the redefinition of $f(x)$, there can arise some
convergence problems, if the search path goes in the vicinity of already
located zeros. In addition, the approach cannot help us to eliminate multiple
zeros, they will be found according to their multiplicity. The latter could be
useful only if we need to know the eigenvalues multiplicity, otherwise the
algorithms will waste time to find an already found zero. We remark that for
the real optical applications, it is unlikely to have multiple zeros of the
scattering function $a(\xi)$ due to the presence of noise that should
typically split the multiple eigenvalue. Because of that issue we do not pay
much attention to the multiplicity in the computation of eigenvalues in our
present study.

At this point, it is pertinent to make a remark on the region, in which the
eigenvalues are sought for. The values of real and imaginary parts of an
eigenvalue can be limited by physical and mathematical reasons. By definition,
the eigenvalues lie in upper half-plane, so $\text{Im}(\xi)>0$ gives one
geometrical bound for the search region. At the same time, the imaginary part
corresponds to the amplitude of the soliton, which is associated with its
energy. Hence, a knowledge of the energy of the signal allows us to impose an
upper limit on the boundary of $\text{Im}(\xi)$. We note that for the contour
integration algorithms we seek the eigenvalues inside the manually defined
region of interest (ROI) in the complex $\xi$-plane. For iterative algorithms,
the iteration paths are not usually controlled (unless some constraints inside
the search routine are additionally imposed) and thus the iteration algorithm
scan incidentally go out of the ROI; this can be an additional drawback for
that group of methods.

###  4.5 Multisolitonic test

Multisolitonic potentials are patently interesting from the perspective of
practical applications, with subsets of multisolitonic eigenvalues being
specifically proposed for optical communication purposes [12, 13, 14, 17]. In
this subsection, we investigate the performance of both iterative and contour
integration methods for the case when the NFT pulse decomposition involves
several eigenvalues, and we need to retrieve their values. For our tests, we
choose the over-soliton potential as in Eq. (15) for the fixed amplitude
$A=5$. This signal has five eigenvalues in its discrete NF spectrum:

$\xi_{k}=(4.5-k)i,\quad k=0\ldots 4.$ (64)

The full decomposition also contains the non-zero continuous spectral data, so
that the situation considered in this section is quite general. To locate the
eigenvalues, it is convenient to use the rectangular ROI for the contour
integration methods because, as it was found in Subs. 4.2, such a contour
ensures a greater accuracy in results. We choose the contour having the shape
of a rectangle in the complex $\xi$-plane with the dimensions: $[-1,1]$ along
the real axis and $[0.1,5]$ along the imaginary axis. Such a rectangle
encompasses all eigenvalues defined by Eq. (64) for our over-soliton profile.

Now we note that our iterative algorithms need adaptation for the
multisolitonic search task, as described in the previous subsection. The
adapted multisoliton-search iterative algorithm scheme is given in Alg. 1.

Algorithm 1 Adaptation of an iterative algorithm to the multisolitonic case
Input: Define function $f(x)$, make initial guess $x_{0}$,
---
expected number of zeros $N$, ROI.
Step 0: Initialize an empty output array $x_{out}$.
Step 1: Define a current function $f_{c}(x)=f(x)$
and a current guess $x_{c}=x_{0}$.
Step 2: Launch the iterative algorithm for $f_{c}(x)$, $x_{c}$
and limited number of iteration steps.
Step 3: Check if zero $x_{i}$ was successfully located in the previous step.
If yes, go to step 4.
If no, go to step 7.
Step 4: Check if the located zero is inside the ROI.
If yes, concatenate the located zero to the output array
$x_{out}=[x_{out},x_{i}]$
and go to step 5.
If no, go to step 7.
Step 5: Check if all expected zeros are located using size($x_{out}$)==$N$.
If yes, go to Output.
If no, go to step 6.
Step 6: Redefine the function to eliminate located zero
$f_{c}(x)=f_{c}(x)/(x-x_{i})$.
Step 7: Choose an initial guess number $x_{c}$ randomly from inside the ROI
and go to step 2.
Output: Estimate the array of zeros $x_{out}$.

In our specific case, we started with the initial guess $\xi_{0}=i$, followed
by up to three attempts to locate the same zero over $10^{5}$ iteration steps
to reach each particular zero point. The resultant error versus runtime
diagram is represented as a bubble chart in Fig. 10 for the different
combinations of the transfer-matrix NFT algorithms from Subs. 3.1 and the
root-finding iterative methods from Subs. 4.1 and 4.2. Here bubble sizes are
inversely proportional to the runtime with the numbers therein indicating the
number of zeros; respective colors identify the numerical algorithm chosen.

Figure 10: The bubble chart shows the relative error of the eigenvalue
location computed via the different root-finding iterative (NR, se, St, Mu,
Si) and contour integral (DL, aDL) methods for over-soliton profile with five
embedded solitonic eigenvalues as in Eq. (64). The digit in each bubble means
the number of ultimately located zeros (the maximum is 5), the error was
assessed as the mean relative error via Eq. (28), the formula was applied for
located zeros only. The size of each particular bubble relatively shows the
inverse runtime, a smaller bubble means a longer run and vice versa; the logic
applies to concentric and overlapping circles. Runtime changes from 0.025 s
(the largest bubble) to 4192 s (the smallest one). The color of each bubble
identifies to which NFT method the root-finding algorithm was coupled.

We see that, in accordance with the remarks made in the previous subsection,
the iterative algorithms equipped with the elimination procedure as in Eq.
(63) may fail to estimate the entire cluster of zeros even when coupled with
the most accurate BO NFT method. At the same time, the contour integration
methods, especially aDL, give a higher error margin.

###  4.6 A hybrid method

As we noted in the previous parts of this section, both types of eigenvalue-
finding approaches have some inherent disadvantages. The iterative algorithms
from Subs. 4.1 can be unstable: we cannot be confident that the iterations
will eventually lead to the correct eigenvalue points (see Figs. 5 and 6).
They also require some additional adaptation to the multisolitonic case to
incorporate the elimination of previously found zeros, and the restriction on
the search region has to be generally imposed. The contour integration
algorithms from Subs. 4.2 do not allow us to reach a high accuracy and take a
comparatively long time.

In this subsection, we present the new hybrid method, which allows us to take
advantage of the best from both approaches whilst simultaneously getting rid
of their respective drawbacks. The main idea of our hybrid method is that we
can use the result of contour integration as the initial guess that is then
supplied to the consequential iterative algorithm (the same strategy was
mentioned in [37] as a way of root-finding refining). This combination allows
us to reach almost any accuracy up to the limitations imposed by the NFT
computation method itself, see Subs. 3.1. The hybrid method presented here
guarantees locations for all eigenvalues as opposed to iterative methods that
often found most eigenvalues but not necessarily all. The description of the
consecutive steps for the hybrid algorithm is as follows.

  * 1.

First, find the location of the approximate zeros’ using one of the contour
integration method drawing a large enough contour. Two key remarks on this
step here:

    * (a)

the integration result does not require to be really accurate, so the
computation time can be reduced;

    * (b)

it allows us to find a good approximation for all zeros that lie inside our
ROI;

  * 2.

Second, we apply a particular iterative method to find more precise location
for each eigenvalue, using the results of the previous step as guess points
(and, eventually, also employing some other data that can be obtained by the
contour integration, i.e. the multiplicity).

To test our new hybrid algorithm, we again employ the over-soliton potential
from the previous subsection 4.5 with five solitonic modes, and plot the error
of the eigenvalues computed in Fig. 11. We test the different combinations of
the NFT algorithms from Subs. 3.1 (we do not use the RK method in this section
as it is too time consuming) combined with contour integration methods from
Subs. 4.2 and iterative methods from Subs. 4.1.

|
---|---
a) DL | b) aDL

Figure 11:  The bubble chart shows the relative error of solitonic eigenvalue
location by using the hybrid method; the results are given for the over-
soliton profile with five eigenvalues. Panel a) corresponds to the DL contour
integration to find the guess values; panel b) does the same by using the aDL
method. As before, the bubble size shows inverse runtimes, a smaller bubble
indicates a longer run and vice versa, the logic applies to both concentric
circles. Runtime changes from 3.27 s (the largest bubble) to 286 s (the
smallest one). The color of each bubble identifies to which NFT method the
root-finding algorithm was coupled.

We see that, in general, runtime reduces comparably with contour integrals
approach (see table 2), but it is still larger than runtime of some iterative
algorithms (see runtime labels in Figs. 5 and 6). We see that the largest
runtime is demonstrated by the Sidi and Steffensen methods, but the latter
show a surprisingly high accuracy when coupled with the BO NFT algorithm. The
NR, Muller and secant algorithms perform fairly similar in terms of runtime,
but the Muller method displays a worse accuracy. If we combine the NR and
secant methods with any of the integral methods, then the utilization of both
the DL and aDL approaches gives similar results. It means, that both the DL
and aDL methods supply a sufficiently good initial approximation to reach our
desired zero.

##  5 Computation of norming constants

The third and last component of the NF spectrum is the norming constants (11)
attributed to each eigenvalue. These parameters define the phase and the
center position of each solitonic degree of freedom [1, 4, 6]. Under the
assumptions listed in Secs. 1, 2, the norming constants are expressed as the
residues of the reflection coefficient $r(\xi)$ from (8) calculated at its
simple poles, Eq. (11), i.e. at the solitonic eigenvalues addressed in the
previous section. So at this point we assume that some appropriate method from
Sec. 4 has been executed and the plausible values for all eigenvalues are now
known with high enough precision.

###  5.1 Ways of the residues evaluation

For the computation of the norming constants, we can equally use each of the
two expressions for the residue. First, the norming constant can be found via
the contour integration:

$c_{j}=\frac{1}{2\pi i}\int_{\gamma_{j}}\\!\\!r(\xi)\,d\xi,$ (65)

where $\gamma_{j}$ is a sufficiently small contour in the $\xi$-plane
encircling the single pole $\xi_{j}$ (considered to have been already
located). Alternatively [3, 6], the norming constant is given by fraction

$c_{j}=\frac{b(\xi_{j})}{a^{\prime}(\xi_{j})},$ (66)

where $a^{\prime}(\xi_{j})$ is derivative of the scattering coefficient
$a(\xi)$ from (7) with respect to $\xi$ evaluated in its simple zero
$\xi_{j}$. The computation of $a(\xi)$, $b(\xi)$ or $r(\xi)$ in any point of
$\overline{{\mathbb{C}}^{+}}$ can be performed using the methods from Subs.
3.1. The computation of $a^{\prime}(\xi)$ can be carried out in the same
program cycle in parallel with the other NFT quantities that do not involve
derivatives as it is described in Subs. 4.3. However, as we shall see, the
computation of the norming constants brings about some additional problems
[17, 40].

###  5.2 Test of the norming constants straightforward computation

Explicit analytical expressions for the norming constants corresponding to the
eigenvalues of three model potentials that we use in this paper are presented
in Subs. 2.2. For the forthcoming tests we use the over-soliton profile (15),
for which the expressions for both eigenvalues and norming constants are
explicit and relatively simple. As for the rectangular potential (21), the
eigenvalues are defined through the solutions of a transcendental equation,
which, in turn, can be found only with a finite accuracy, and this can
restrict our comparative analysis. Therefore we do not use the rectangular
potential for testing the norming constant methods accuracy.

In our test, we compare the performance of the different NFT transfer-matrix
algorithms from Subs. 3.1 and 4.3; these are combined with the two residue
expressions given in the previous subsection. From Fig. 12a) we can see that
the fraction formula (66) employed for the residue computation, improves the
result accuracy as compared to the integral definition (65) (both combined
with different NFT methods). The main feature that we can extract from the
analysis of Fig. 12 is that beyond a critical value of the amplitude
($A\lesssim 1$ for this particular case), the computational error starts to
rise rapidly. Manipulating the number of discretization points for the
computation of the integral from (65) does not change this error signature,
but it evidently influences the runtime significantly. The integral formula
(65) gives indistinguishably similar results independent of the NFT algorithm
chosen. The fraction formula (66) gives the best accuracy for the BO method
and the worst accuracy attributed to the CN one. In general, our current
results reveal that both formulae of the residue combined with different
available NFT methods are relatively inaccurate, and this can be an important
degrading factor in practical applications.

|
---|---
a) | b)

Figure 12: Residue relative error on over-soliton signal’s amplitude for
different NFT algorithms and residues computation methods ($n=2^{12}$, $L=20$)
for pane a) fraction and contour integral formulae (number of points along
integration contour rise from transparent to solid lines in the range
$n_{p}=[20\dots 100]$ with the increment step $\Delta n_{p}=20$); panel b)
fraction formula with application of improved scheme of $b(\xi)$ evaluation.

Authors of [29] also showed that the AL, CN, Euler and central differences
methods applied for the ZSS solution gave the similar unimpressively-small
accuracy of norming constant computation. In [31], it was also found that here
the RK method fares worse than the BO approach. In [40], the author used the
BO method and compared its performance with AL, CN for the norming constants
computation.

###  5.3 Discussion and the improvement of norming constants computation
accuracy

The problem of the numerical computation of residues lies in the properties of
the spectral functions $a(\xi)$, $b(\xi)$ and $r(\xi)$. When performing the
analytical continuation of these functions into the complex $\xi$-plane, we
should carefully check when this operation is indeed legitimate. As stated
above, the analytic continuation of $b(\xi)$ is limited, in general, by the
rate of decay of the ZSS potential. On the other hand, when dealing with the
function $a(\xi)$, its analytic continuation can be performed over the whole
$\mathds{C}^{+}$ [1], where it is bounded and, moreover, $a(\xi)\to 1$ as
$\xi\to\infty$, and so this property adds on our motivation to use the
function $a(\xi)$ for the solitonic eigenvalues.

In the case when $b(\xi)$ admits the analytic continuation into the whole
complex plane or a part of it covering the location of the eigenvalues
$\\{\xi_{j}\\}$, the norming constants $b_{j}$ associated to $\xi_{j}$ can be
calculated as (see (9)) $b_{j}=b(\xi_{j})$. However, $b(\xi)$ is likely to
grow exponentially when the imaginary part of $\xi$ increases. Numerical
computations are also limited by the speed of the processor and by the
available memory size (e.g. double precision numbers are bounded by
approximately 1.8e308). When we reach this limiting value, the accuracy of the
computation would be naturally affected. In order to understand the scale of
problem arising in the accurate norming constant computation, we compared the
limitations caused by the analytic continuation problems and by the
computational reasons listed above.

For the over-soliton potential with amplitude $A=5.25$, we calculate the
values of $b(\xi)$ for a purely imaginary $\xi$. The curves for the
numerically computed $b(\xi)$ (see Fig. 13) grow extremely fast, and so it is
not possible to distinguish individual lines for the different NFT algorithms.
It is also seen that the rise of $b(\xi)$ starts well before the decay rate
$d$ coming to play (it is responsible for the band in the $\mathbb{C}$, where
$b(\xi)$ is certainly analytical and shown by the green line). Therefore, our
NFT algorithms are unable to evaluate the correct value of $b(\xi)$ (the gray
line on the plot) even for sufficiently small $\text{Im}\xi$. On the other
hand, the evaluation of the norming constants in terms of $b(\xi)$ requires
the values of $b(\xi)$ at the eigenvalues (marked by the black lines), where,
as we see, the spectral function $b(\xi)$ is computed incorrectly.

Figure 13: The computed values of $b(\xi)$ (red line) along imaginary axis of
$\xi$, the value of the explicit expression (gray line), eigenvalues (black
lines), and the potential decay rate (green line), calculated for the over-
soliton profile (15) with $A=5.25$.

However, we can adapt the NFT algorithms to provide a more accurate
computation of $\\{b_{j}\\}$ associated with the eigenvalues $\\{\xi_{j}\\}$.
Indeed, the detailed control of each computational step shows that if one uses
Eq. (14) for the calculation of $b$, the aforementioned divergence happens at
the truncation edges of the potential (exponential factor grows fast for large
$\text{Im}(\xi)$ and $L$). Thus it is possible to simplify the evaluation of
$b_{j}$ (i) by estimating it at the centre of the evaluation interval,
avoiding the edges, and (ii) by getting rid of the exponentials in the initial
conditions as in Eq. (4). This idea was effectively harnessed in [17] for the
solution of ZZS employing a more straightforward Euler method and in [40] for
BOmod transfer-matrix approach. It can be realised by considering the
evolution of the so-called left wave (defined by its asymptotic as
$t\to-\infty$) from the left edge towards the centre of the interval, and the
evolution of the right wave (defined by its asymptotic as $t\to+\infty$) from
the right edge to the centre towards the left wave, and then by using Eq. (9)
at $t=0$ to calculate $b_{j}$.

Using the wave function envelopes for both $\Phi(t,\xi)$ and $\Psi(t,\xi)$,
defined as

$G(t,\xi)=\Phi(t,\xi)e^{i\xi t},\qquad H(t,\xi)=\Psi(t,\xi)e^{-i\xi t},$ (67)

we can evaluate the coupled systems for both vectors $G(t)$ and $H(t)$:

$\frac{d}{dt}G(t,\xi)=\left(\begin{matrix}0&q(t)\\\
-q(t)^{*}&2i\xi\end{matrix}\right)G(t,\xi),$ (68)
$\frac{d}{dt}H(t)=\left(\begin{matrix}-2i\xi&q(t)\\\
-q(t)^{*}&0\end{matrix}\right)H(t,\xi).$ (69)

After the potential truncation, the initial value of the vector
$H(L,\xi)=(0,1)^{T}$ evolves from $t=L$ towards $t=0$, whilst the vector
$G(-L,\xi)=(1,0)^{T}$ evolves from $t=-L$ to $t=0$. At the point $t=0$, the
desired quantity $b_{j}$ can be obtained from the following relation (cf.
(9)):

$G(t=0,\xi_{j})=H(t=0,\xi_{j})b_{j}.$ (70)

The evolution of $G(t,\xi)$ and $H(t,\xi)$ can be performed similarly to Subs.
3.1, i.e. using the transfer matrices. The implementation of the BO approach
(we used matrix exponential of ODE matrix, see Eqs. (31), (32)) leads to the
following transfer matrix for the left wave

$T_{m}^{\text{(BOleft)}}=e^{i\xi\Delta t}\left(\begin{matrix}\cos\kappa\Delta
t-i\xi/\kappa\sin\kappa\Delta t&q_{n}/\kappa\sin\kappa\Delta t\\\
-q_{n}^{*}/\kappa\sin\kappa\Delta t&\cos\kappa\Delta
t+i\xi/\kappa\sin\kappa\Delta t\end{matrix}\right)$ (71)

and

$T_{m}^{\text{(BOright)}}=e^{-i\xi\Delta
t}\left(\begin{matrix}\cos\kappa\Delta t-i\xi/\kappa\sin\kappa\Delta
t&q_{n}/\kappa\sin\kappa\Delta t\\\ -q_{n}^{*}/\kappa\sin\kappa\Delta
t&\cos\kappa\Delta t+i\xi/\kappa\sin\kappa\Delta t\end{matrix}\right)$ (72)

for the right wave (here $\kappa=\sqrt{|q_{m}|^{2}+\xi^{2}}$). The AL-type
approach ((using Euler method and exponent first order decomposition, see
(33)) leads to the following transfer matrices for the left and right
envelopes, correspondingly:

$T_{m}^{\text{(ALleft)}}=\frac{1}{\sqrt{1+\Delta
t^{2}q_{m}^{2}}}\left(\begin{matrix}1&q_{m}\Delta t\\\ -q_{m}^{*}\Delta
t&e^{2i\xi\Delta t}\end{matrix}\right),$ (73)
$T_{m}^{\text{(ALright)}}=\frac{1}{\sqrt{1+\Delta
t^{2}q_{m}^{2}}}\left(\begin{matrix}e^{-2i\xi\Delta t}&q_{m}\Delta t\\\
-q_{m}^{*}\Delta t&1\end{matrix}\right).$ (74)

We emphasize that the presented way of definition and computation of $b$ is
appropriate only for the discrete ZSS eigenvalues.

We tested our new methods of calculation of $\\{b_{j}\\}$ to find the residues
through the fraction formula presented in Eq. (66) (with $b(\xi_{j})$ replaced
by $b_{j}$); our results are summarised in Fig. 12b). We can readily observe
that by employing the method described above we get more accurate norming
constants for the wider range of the amplitude of the potential. In contrast
to our previous results, the CN method shows better accuracy than both AL and
ALmod methods, but the CN method takes more time than all other methods. The
BO method gives the best accuracy and the weakest dependence of the outcome on
the amplitude.

##  6 Conclusion

In this study, we have formulated a comprehensive list of options for NFT
computation, and then have made categorical comparisons of their relative
advantages and disadvantages, principally focusing on runtime and accuracy
optimization. The focus here has been to obtain best possible numerical
approximation algorithm of the given spectral data.

First, the calibration of different transfer-matrix algorithms’ performance
was done using the continuous part of the spectral data as a measure for the
method’s accuracy analysis. In particular, it was shown that the BO method is
usually superior to other alternatives in terms of the runtime and the
accuracy of the result obtained, a conclusion that complies with some earlier
existing studies on the NFT methods performance.

However, a major incentive of this work was the efficient computation of the
eigenvalues and norming constants associated with each solitonic degree of
freedom (i.e. in the full discrete spectral data associated with a given
profile). For the computation of the eigenvalues, we first applied different
iterative algorithms (involving the derivative computation or avoiding it at
each step) combined with different transfer-matrix methods for the ZSS
solution, using three different profiles with the known NF spectrum. At the
beginning, we analyze the regions of convergence of each NFT method combined
with a particular iterative scheme that relied on a reasonably accurate guess
point that is then used for a search run. Surprisingly enough, we find that
the Muller method typically provides the largest convergence region; the
relatively large convergence basin is observed for the derivative-free Sidi
method as well, while the Newton-Raphson, secant and Steffensen’s iterative
approaches show progressively poorer convergence.

This is followed by a description of a new class of methods for the search of
eigenvalues (the zeros of the spectral function $a(\xi)$) based on the
evaluation of contour integrals in the complex plane of the spectral parameter
$\xi$. Although we found that these methods are often slower and less accurate
compared to “more traditional” iterative search algorithms, they generally
allow us the freedom of a lax choice for the initial guess point during the
eigenvalue search. This task is often difficult to fulfil in realistic
applications, where the properties of the priories are not known a priori. In
this section, we present a new hybrid algorithm that allows us to combine
merits of both the iterative and contour integral method classes. The new
method relies on the approximate initial evaluation of the NFT data using the
contour integrals, and then using the obtained data as a good guess for an
iterative algorithm following the initial search. This scheme leads to faster
evaluation of the eigenvalues with high enough accuracy while largely being
able to avoid convergence and initial guess issues. This new method for
eigenvalues search is one of the main results of our present paper.

Finally, we address the problem of the computation of the norming constants
using different algorithms. As same as the earlier studies suggest [17, 29],
the straightforward application of the single-directional NFT algorithms
usually results in lower accuracy of the estimated norming constants. This
adverse finding can be attributed to the properties of the analytic
continuation of the function $b(\xi)$ into the upper half-plane of $\xi$,
which is used in the norming constant evaluation. To rectify this problem, we
utilise a numerical scheme first suggested in [17] that employs simultaneous
right and left scattering. We then use this bi-directional approach combined
with different transfer-matrix methods. Our results prove that this approach
results in an improved accuracy of the norming constants computation.

We believe that this new regime of modified NFT algorithms has led to a
framework that will concomitantly optimize speed against accuracy, thereby to
drive an adjustable numerical routine for the evaluation of NF spectral data
tailored specifically to the problem in hand. This subjectivity aspect can be
highly valuable in view of the active current progress of optical
communication methods based on the NF spectrum modulation [3].

## Acknowledgement

JEP acknowledges the support from the UK EPSRC Programme Grant UNLOC
EP/J017582/1. JEP and DS are thankful to the Erasmus+ ERC mobility programme
between the Aston University and Kharkiv National University that helped us to
launch the collaborative activity. AC acknowledges the RISE-FRAMED grant.
Discussions with Sergei Turytsyn and Sotos Generalis are gratefully
acknowledged.

Data used in the figures in this paper are available on
http://doi.org/10.17036/researchdata.aston.ac.uk.00000321.

## References

## References

  * [1] M. J. Ablowitz and H. Segur, Solitons and the Inverse Scattering Transform (SIAM, 1981).
  * [2] L. F. Mollenauer and J. P. Gordon, Solitons in Optical Fibers: Fundamentals and Applications (Academic Press, 2006).
  * [3] S. K. Turitsyn, J. E. Prilepsky, S. T. Le, S. Wahls, L. L. Frumin, M. Kamalian, and S. A. Derevyanko, “Nonlinear Fourier transform for optical data processing and transmission: advances and perspectives,” Optica 4, 307-322 (2017).
  * [4] V. E. Zakharov and A. B. Shabat, “Exact theory of 2-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media,” Sov. Phys.-JETP 34, 62–69 (1972).
  * [5] R. Beals and R. R. Coifman, ”Inverse scattering and evolution equations”, Comm. Pure Appl. Math. 38(1), 29-42 (1985).
  * [6] M. I. Yousefi and F. R. Kschischang, “Information transmission using the nonlinear Fourier transform, Part I: Mathematical Tools.” IEEE Trans. Inform. Theory 60, 4312–4328 (2014).
  * [7] E. Stein, ”Harmonic Analysis, Real variable methods, orthogonality, and oscillatory integrals”, Princeton Mathematical Series 43, Princeton University Press, NJ (1993).
  * [8] L. D. Faddeev and L. Takhtajan, Hamiltonian Methods in the Theory of Solitons (Springer-Verlag, 1987).
  * [9] R.-J. Essiambre, G. Kramer, P.J. Winzer, G.J. Foschini, and B. Goebel, “Capacity limits of optical fiber networks,” J. Lightwave Technol. 28, 662–701 (2010).
  * [10] P. Bayvel, R. Maher, T. Xu, G. Liga, N. A. Shevchenko, D. Lavery, A. Alvarado, and R. I. Killey, “Maximizing the optical network capacity,” Phil. Trans. R. Soc. A 374, 20140440 (2016).
  * [11] J. C. Cartledge, F. P. Guiomar, F. R. Kschischang, G. Liga, and M. P. Yankov, “Digital signal processing for fiber nonlinearities [Invited],” Optics Express 25, 1916-1936 (2017).
  * [12] H. Terauchi and A. Maruta, “Eigenvalue Modulated Optical Transmission System Based on Digital Coherent Technology,” in The 10th Conference on Lasers and Electro-Optics Pacific Rim, and the 18th OptoElectronics and Communications Conference/Photonics in Switching (CLEO-PR & OECC/PS), Kyoto, Japan, Paper WR2-5, 2013.
  * [13] A. Hasegawa and T. Nyu, “Eigenvalue communication,” J. Lightwave Technol. 11, 395–399 (1993).
  * [14] S. Hari, F. Kschischang, and M. Yousefi, “Multi-eigenvalue communication via the nonlinear Fourier transform,” in 27th Biennial Symposium on Communications (QBSC), Kingston, ON, Canada, pp. 92–95, 2014.
  * [15] Z. Dong, S. Hari, T. Gui, K. Zhong, M. I. Yousefi, C. Lu, P.-K. A. Wai, F. R. Kschischang, and A. P. T. Lau, “Nonlinear Frequency Division Multiplexed Transmissions based on NFT,” IEEE Photon. Tech. Lett. 27, 1621–1623 (2015).
  * [16] S. Hari, M. I. Yousefi, and F. R. Kschischang, “Multieigenvalue Communication”, J. Lightwave Technol. 34, 3110–3117 (2016).
  * [17] S. Hari and F. R. Kschischang, “Bi-Directional Algorithm for Computing Discrete Spectral Amplitudes in the NFT,” J. Lightwave Technol. 34, 3529–3537 (2016).
  * [18] J. E. Prilepsky, S. A. Derevyanko, K. J. Blow, I. Gabitov, and S. K. Turitsyn, “Nonlinear Inverse Synthesis and Eigenvalue Division Multiplexing in Optical Fiber Channels.” Phys. Rev. Lett. 113, 013901 (2014).
  * [19] S. T. Le, I. D. Philips, J. E. Prilepsky, P. Harper, A. D. Ellis, and S. K. Turitsyn, “Demonstration of Nonlinear Inverse Synthesis Transmission over Transoceanic Distances,” J. Lightwave Technol. 34, 2459–2466 (2016).
  * [20] E. G. Turitsyna and S. K. Turitsyn, “Digital signal processing based on inverse scattering transform.” Opt. Lett. 38, 4186–4188 (2013).
  * [21] S. Wahls, S. T. Le, J. E. Prilepsky, H. V. Poor, and S. K. Turitsyn, “Digital Backpropagation in the Nonlinear Fourier Domain,” in Proceedings of IEEE 16th International Workshop in Signal Processing Advances in Wireless Communications (SPAWC), Stockholm, Sweden, pp. 445–449, 2015.
  * [22] H. Bülow, “Experimental Demonstration of Optical Signal Detection Using Nonlinear Fourier Transform,” J. Lightwave Technol. 33, 1433–1439 (2015).
  * [23] S. Wahls and H. V. Poor, “Introducing the fast nonlinear Fourier transform,” in Proceedings of International Conference on Acoustics, Speech, and Signal Processing (ICASSP), Vancouver, Canada, pp. 5780–5784, 2013.
  * [24] S. Wahls and H. V. Poor, “Fast Numerical Nonlinear Fourier Transforms,” IEEE Trans. Inform. Theory 61, 6957–6974 (2015).
  * [25] J. W.Cooley and J. W. Tukey.“An Algorithm for the Machine Calculation of Complex Fourier Series”, Math. Comp. 19, 297–301 (1965).
  * [26] R. Kashyap, Fiber Bragg Gratings, 2nd ed. (Elsvier, 2010).
  * [27] M. Jaulent, “The inverse scattering problem for LCRG transmission lines,” J. Math. Phys. 23, 2286 (1982).
  * [28] J. Satsuma and N. Yajima, “Initial Value Problems of One-Dimensional Self-Modulation of Nonlinear Waves in Dispersive Media,” Prog. Theor. Phys. Suppl. 55, 284–306 (1974).
  * [29] M. I. Yousefi and F. R. Kschischang, “Information transmission using the nonlinear Fourier transform, Part II: Numerical Methods,” IEEE Trans. Inform. Theory 60, 4329–4345 (2014).
  * [30] G. Boffetta and A. Osborne, “Computation of the direct scattering transform for the nonlinear schroedinger equation,” J. Comput. Phys. 102, 252–264 (1992).
  * [31] S. Burtsev, R. Camassa, and I. Timofeyev, “Numerical algorithms for the direct spectral transform with applications to nonlinear Schrödinger type systems,” J. Comput. Phys. 147, 166–186 (1998).
  * [32] M. Kamalian, J. E. Prilepsky, S. T. Le, and S. K. Turitsyn, “Periodic nonlinear Fourier transform for fiber-optic communications, Part I: theory and numerical methods,” Opt. Express 24, 18353–18369 (2016).
  * [33] A. Sidi, “Generalization Of The Secant Method For Nonlinear equations,” Applied Mathematics E-Notes 8, 115–123 (2008).
  * [34] L. W. Johnson and D. R. Scholz, “On Steffensen’s Method,” SIAM J. Numer. Anal. 5, 296–302 (1968).
  * [35] D. E. Muller, “A Method for Solving Algebraic equations Using an Automatic Computer,” Mathematical Tables and Other Aids to Computation, 10, 208–215 (1956).
  * [36] M. I. Yousefi and F. R. Kschischang, “Information Transmission Using the Nonlinear Fourier Transform, Part III: Spectrum Modulation” IEEE Trans. Inform. Theory 60, 4346–4369 (2014).
  * [37] L. M. Delves and J. N. Lyness, “A numerical method for locating the zeros of an analytic function,” Math. Comp. 21, 543–560 (1967).
  * [38] P. Kravanja, T. Sakurai, and M. van Barel, “On Location Clusters of Zeros of Analytic Functions” BIT 39, 646–682 (1999).
  * [39] M.K. Jain “Numerical Methods For Scientific And Engineering Computation” (New Age International, 2003).
  * [40] V. Aref “Control and Detection of Discrete Spectral Amplitudes in Nonlinear Fourier Spectrum”, preprint, arXiv:1605.06328v1 (2016).