In a series of papers [1, 2, 3], it was shown recently that exponential linear multistep methods (LMMs) provide a natural setting for higher-order convergent fast nonlinear Fourier transform (NFT). This followed from a simple observation that the transfer matrices obtained are amenable to FFT-based fast polynomial arithmetic. Note that in the earliest works on fast NFTs , the Ablowitz-Ladik method can be interpreted as the exponential Euler method. In this paper, we use the exponential Runge-Kutta methods to obtain a family of fast NFTs (provided that the nodes are equispaced). In particular, we present two fast NFTs with fourth order accuracy based on fourth order Runge-Kutta methods. The structure of the transfer matrix reveals that such methods are superior to those based on LMMs in terms of complexity while the numerical tests reveal that they also have a superior accuracy-complexity trade-off. The first algorithm is based on the classical fourth order (explicit) Runge-Kutta method which has been studied by several authors in the context of NFTs [5, 6]. The second method uses the three-stage Lobatto IIIA (implicit) Runge-Kutta method.
For moderately small step-sizes, most fast methods yield poor accuracy specially corresponding to the large values of the spectral parameter. An error analysis of such integrating factor based methods  shows that the error terms contain positive powers of the spectral parameter which necessitates the use of smaller step-sizes in order to keep the error low. On the contrary, integrators based on Magnus series expansion are known to have error terms that contain negative powers of the spectral parameter. In this letter, we follow the recipe presented by Blanes et al.  to develop a fourth order Magnus method and a fourth order commutator-free (CF) Magnus method , both of which take samples of the potential on an equispaced grid in order to compute the NFT. Let us emphasize that our CF method is different from those considered in 
where the samples of the potential are needed on the Gauss-Legendre nodes (the authors generate the samples by interpolation on an equispaced grid, locally). The fourth-order (standard) Magnus method happens to be faster than the corresponding CF method on account of the fact that there is an additional matrix exponential introduced in the CF method in order to avoid the use of commutators. The accuracy-complexity trade-off, however, is similar for the two methods. Finally, we also present a fast variant of the CF method of order four (formally) by employing the fourth order splitting on the lines of[10, 6]. Despite the well-known limitation imposed on the order and stability of such techniques as demonstrated by Sheng , we do not find any reduction of order within the double precision arithmetic. It is noteworthy that the authors in  found the aforementioned splitting worsen in accuracy after a certain step-size. It is also not clear from their analysis if this scheme is convergent in their setting.
We begin our discussion with a brief review of the scattering theory closely following the formalism presented in . The nonlinear Fourier transform of any signal is defined via the Zakharov-Shabat (ZS) scattering problem which can be stated as follows: For and ,
where . The potential is defined by and with (). Here, is known as the spectral parameter and is the complex-valued signal. The solution of the scattering problem (1), henceforth referred to as the ZS problem, consists in finding the so called scattering coefficients which are defined through special solutions of (1) known as the Jost solutions which are linearly independent solutions of (1) such that they have a plane-wave like behavior at or . The Jost solution of the second kind, denoted by , has the asymptotic behavior as . The asymptotic behavior as determines the scattering coefficients and for . In this letter, we primarily focus on the continuous spectrum, also referred to as the reflection coefficient, which is defined by for .
Ii The Numerical Scheme
Ii-a Runge-Kutta Method
In this section, we will develop the integrating factor based exponential Runge-Kutta (RK) method for the numerical solution of the ZS problem. Following [1, 3], we begin with the transformation so that (1) becomes with whose entries are and . Let the step size be and the quantities be ordered so that the nodes within the step can be stated as . For the potential sampled at these nodes, we use the convention , and . In order for the resulting discrete system to be amenable to FFT-based fast polynomial arithmetic, it is sufficient to have each of the
’s belong to the set of uniformly distributed nodes in. A -stage RK method is characterized by the nodes , and, the weights and . Introducing the intermediate stage quantities for , we have
This system of equations can be solved in any computer algebra system to obtained the transfer matrix connecting the vectorsto .
Setting , the Lobatto IIIA method (labelled as IRK) of order  simplifies to
The fourth order classical RK method (labelled as ERK) simplifies to the form (3) with and the entries of the transfer matrix are given by
Ii-A1 Scattering coefficients
Let the computational domain be and set . Let the number of steps be and . Also, let and . The grid is defined by with where is the number of samples. Let the potential be supported in and we assume for convenience. In order to represent the Jost solutions, we introduce the polynomial vector
and the polynomial . Consider the Jost solution . For the Lobatto IIIA method, the Jost solution can be stated as , with
where is the constant part of and
From , it follows that the discrete scattering coefficients are given by
where . For the classical fourth order RK method, ; therefore, and . The principal branch for the discrete scattering coefficients here works out to be . This follows from the principle branch of the individual transfer matrices. The nodes lead to which is not in the standard form for FFT algorithms to be used. Therefore, we would like to work with nodes so that
and pad the input vector with zeros. Following as in, the complexity of computing the scattering coefficient works out to be .
Ii-B Standard and commutator-free Magnus methods
Ii-B1 Magnus method
Let us assume that the solution of the ODE (1) can be written as for , then has a series representation known as the Magnus series . Truncating this series to achieve the desired order of accuracy yields a family of numerical schemes known as Magnus method. For the ZS problem, Magnus integrators preserve the Lie group structure of the Jost solution and its accuracy does not worsen with increasing . In designing a fourth order Magnus integrator (labelled as M), we follow the method due to Blanes et al. [7, 8]: Defining
the method proceeds by expanding in terms of the quantities using the Magnus series. For the fourth order method, setting , we have [7, 8] . Evaluating and upto fourth order accuracy using the three-point Gauss quadrature involving Legendre-Gauss-Lobatto (LGL) nodes, , the numerical scheme for the ZS problem can be stated as where,
where , and . For the purpose of comparison, we would also like to consider the Magnus method with one point Gauss quadrature (labelled as M) which is of order .
Ii-B2 Commutator-free Magnus method
Blanes and Moan  have constructed fourth-order commutator-free (CF) methods that are based on Magnus method. Using the quantities defined above, the CF method (labelled as CF) can be stated as
Evaluating and upto fourth-order accuracy using the three-point Gauss quadrature involving LGL nodes, we have where ,
and . Given that there are two matrix exponentials involved, the CF Magnus method has higher complexity than that of the standard Magnus method.
Ii-B3 A fast variant
The CF Magnus method further allows us to obtain a fast NFT algorithm (labelled as SCF) via splitting of the matrix exponential. Consider a formally fourth order splitting :
where the entries of the matrix are
where . The discrete scattering coefficients can be written as and . The principal branch for the discrete scattering coefficients here works out to be . This again follows from the principle branch of the individual transfer matrices. As before, the nodes lead to which is not in the standard form for FFT algorithms to be used. Therefore, we would like to work with nodes so that and pad the input vector with zeros.
Iii Numerical Tests and Conclusion
For the numerical experiments, we employ the well-known secant-hyperbolic potential given by , () for which the scattering coefficients are given in . We set the computational domain to be and let . Let be the principal branch; then, the error in computing is quantified by
where the integrals are computed using the trapezoidal rule. Similar consideration applies to . For the purpose of testing, we include the implicit Adams method presented in  which are labelled as IA with . The method IA is identical to the trapezoidal rule, therefore, we use the label TR. The convergence analysis is carried out in Fig. 2 and the trade-off between accuracy and complexity is presented in Fig. 2. In terms of accuracy, the CF outperforms every other method with M being a close second. However, the accuracy-complexity trade-off is similar for the two methods. The ‘fast’ methods evidently lower complexity at the cost of accuracy. The RK methods (ERK and IRK) outperform all the other ’fast’ methods in terms of accuracy-complexity trade-off (see Fig. 2); however, with increasing signal strength, the ‘slow’ methods becoming equally competitive. In fact, at moderately small step-sizes, the ‘slow’ methods far outperform the ‘fast’ methods (see Fig. 3) with increasing signal strength.
-  V. Vaibhav, “Fast inverse nonlinear Fourier transformation using exponential one-step methods: Darboux transformation,” Phys. Rev. E, vol. 96, p. 063302, 2017.
-  ——, “Fast inverse nonlinear Fourier transform,” Phys. Rev. E, vol. 98, p. 013304, 2018.
-  ——, “Higher order convergent fast nonlinear Fourier transform,” IEEE Photonics Technol. Lett., vol. 30, no. 8, pp. 700–703, 2018.
-  S. Wahls and H. V. Poor, “Fast numerical nonlinear Fourier transforms,” IEEE Trans. Inf. Theory, vol. 61, no. 12, pp. 6957–6974, 2015.
-  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., vol. 147, no. 1, pp. 166–186, 1998.
-  S. Chimmalgi, P. J. Prins, and S. Wahls, “Fast nonlinear Fourier transform algorithms using higher order exponential integrators,” 2018, arXiv:1812.00703[eess.SP]. [Online]. Available: https://arxiv.org/abs/1812.00703
-  S. Blanes, F. Casas, and J. Ros, “Improved high order integrators based on the magnus expansion,” BIT Numer. Math., vol. 40, no. 3, pp. 434–450, 2000.
-  S. Blanes and P. C. Moan, “Fourth- and sixth-order commutator-free Magnus integrators for linear and non-linear dynamical systems,” Appl. Numer. Math., vol. 56, no. 12, pp. 1519–1537, 2006.
-  S. Chimmalgi, P. J. Prins, and S. Wahls, “Nonlinear Fourier transform algorithm using a higher order exponential integrator,” in Signal Processing in Photonic Communications, 2018, p. SpM4G.5.
-  S. Descombes, “Convergence of a splitting method of high order for reaction-diffusion systems,” Mathematics of Computation, vol. 70, no. 236, pp. 1481–1501, 2001. [Online]. Available: http://www.jstor.org/stable/2698737
Q. Sheng, “Solving linear partial differential equations by exponential splitting,”IMA J. Numer. Anal., vol. 9, no. 2, pp. 199–212, 1989.
-  M. J. Ablowitz, D. J. Kaup, A. C. Newell, and H. Segur, “The inverse scattering transform - Fourier analysis for nonlinear problems,” Stud. Appl. Math., vol. 53, no. 4, pp. 249–315, 1974.
J. C. Butcher,
Numerical Methods for Ordinary Differential Equations, 2nd ed. San Francisco: John Wiley & Sons, Ltd, 2003.