Exponential Fourth Order Schemes for Direct Zakharov-Shabat problem

We propose two finite-difference algorithms of fourth order of accuracy for solving the initial problem of the Zakharov-Shabat system. Both schemes have the exponential form and conserve quadratic invariant of Zakharov-Shabat system. The second scheme contains the spectral parameter in exponent only and allows to apply the fast computational algorithm.


Introduction
Recently, there has been great interest in the so-called nonlinear Fourier transform (NFT), which is a generalization of the ordinary linear Fourier transform. The NFT makes it possible to find exact solutions for nonlinear equations, such as the nonlinear Schrödinger equation (NLSE) and the Korteweg de Vries (KdV) equation, interpreting the evolution of optical fields by analogy with the linear Fourier transform as an evolution of a certain set of frequencies. For the first time, this idea for the NLSE was proposed by Zakharov and Shabat in 1971 [1]. They showed that the NLSE can be integrated by the inverse scattering problem (IST) method, previously applied to the Korteweg de Vries (KdV) equation. The NLSE describes the envelope for wave beams, therefore it is used in many areas of physics where there are wave systems. Despite a large number of articles [2,3,4] devoted to NFT, the development of the accurate and fast numerical algorithms for NFT still remains an actual mathematical problem.
In this paper, we consider the direct spectral problem which consists in the numerical solution of the Zakharov-Shabat (ZS) system. The integration of this system is the first step in the general scheme of NFT. In addition, we limit ourselves to constructing one-step finite-difference schemes for solving the ZS system and do not touch the question of methods for finding discrete spectral parameters and phase coefficients for them.
We present the general necessary conditions for the transition operator for fourth-order one-step difference schemes for linear systems of first-order differential equations. Then we give two examples of such schemes. Both schemes are the exponential fourth order ones, and we show their connection with the Magnus decomposition. The main property of such schemes is to conserve the quadratic invariant for continuous spectral parameters. The second scheme contains the spectral parameter only in the exponent. This property allows one to apply a fast algorithm for the numerical solution of the ZS system.
The final part of the article presents comparisons of numerical computations using the proposed schemes and other wellknown schemes: the Boffetta-Osborne second-order scheme [5], Runge-Kutta fourth-order scheme and the fourth-order conservative scheme [6], which is a generalization of the Boffetta-Osborn scheme.

The direct Zakharov-Shabat problem
The standard NLSE is a basic model for the pulse propagation along an ideally lossless and noiseless fiber i ∂q ∂z + σ 2 where q = q(t, z) is a slow-varying complex optical field envelope, the variable z is the distance along the optical fiber, t is a time variable; σ = −1 and σ = 1 corresponds to the normal and anomalous dispersion, respectively. Eq. (1) is written in the moving coordinate system and describes the propagation of pulses q(t, z) in optical fibers. Therefore, the initial data are almost stationary, and the Cauchy problem is solved with the initial conditions as follows: The mathematical method suggested by Zakharov and Shabat [1] allows to integrate the NLSE. The method, widely known as the Nonlinear Fourier Transform (NFT), allows transforming signal into nonlinear Fourier spectrum, which defined by the solution of the Zakharov-Shabat problem (ZSP).
The equation (1) can be written as a condition of compatibility of two linear equations where Ψ(t) is a complex vector function of a real argument t, The first equation in (3) is the eigenvalue problem for the operator L. For σ = −1 the operator L is Hermitian (L = L † ≡ (L * ) T ), therefore the complex spectral parameter ζ = ξ + iη becomes real ζ = ξ ∈ R. There is no such restriction for σ = 1. In this case the problem has the continuous and discrete spectra. The continuous spectrum lies on the real axis and the discrete spectrum is in the upper half plane Im(ζ) > 0 .
Also the first equation in (3) can be rewritten as an evolutionary system where q = q(t, z) and Here z is a parameter, that we will skip further.
For σ = 1 the matrix Q is skew-Hermitian Q = −Q † . This leads to conservation of the quadratic energy invariant for the real spectral parameters. Moreover, the system (5) can be written in a gradient form as follows: where H = |ψ 1 | 2 + σ|ψ 2 | 2 . For real spectral parameters ζ = ξ the matrix J is skew-Hermitian for any σ = ±1 and, consequently, H conserves. The value H and the matrix Q can be written using Pauli matrices σ 0 and σ 3 as follows: Assuming that q(t) decays rapidly when t → ±∞, the specific solutions (Jost functions) for ZSP (5) can be derived as: and Then we obtain the Jost scattering coefficients a(ξ) and b(ξ) as follows: The functions a(ξ) and b(ξ) can be extended to the upper half-plane ξ → ζ, where ζ is a complex number with the positive imaginary part η = Imζ > 0. The spectral data of ZSP (5) are determined by a(ζ) and b(ζ) in the following way: (1) K zeros of a(ζ) = 0 define the discrete spectrum {ζ k }, k = 0, K − 1 of ZSP (5) and phase coefficients (2) the continuous spectrum is determined by the reflection coefficient r(ξ) = b(ξ)/a(ξ), ξ ∈ R.
These spectral data were defined using the "left" boundary condition (8). Both conditions (8) and (9) can be used to calculate the coefficient b(ζ k ) of the discrete spectrum: For real values of the spectral parameter ζ = ξ we have invariant H. Taking into account the boundary conditions (8), we get H = 1.
In addition, the following trace formula is valid [7]: which connects the NLSE integrals C n with the coefficient a(ξ) and the discrete spectrum ζ k . The first integrals have the form This formula with n = 0 is called the Parseval nonlinear equality and is used to verify the numerical calculations and the consistency of the continuous and discrete spectra found. The first term on the right-hand side of Eq. (13) refers to the continuous spectrum energy: We solve the system (5). The matrix Q(t) linearly depends on the complex function q(t) which is given in the whole nodes of the uniform grid with a step τ on the interval [−L, L]. Let us note main features of the discrete problem: • Since the matrix Q is defined on a uniform grid, the unknown function Ψ must also be computed on a uniform grid with the same step. Therefore, the Runge-Kutta methods (RK) cannot be used on such grid. If, for example, we use an explicit 4th order RK scheme, then we need to take the computational grid with a double step 2τ . In this case, values of Q n will be used unequally. • For small values of the potential |q(t)| << |ζ| and Im ζ > 0, ZSP has exponentially growing and decreasing solutions, thus A-stability of finite-difference methods is required [8]. The method is called A-stable if all solutions of the equation ∂x/∂t = λx tend to zero at Re λ < 0 and fixed step h. The second barrier of Dahlquist restricts the use of multi-step methods [9]. It means that there are no explicit A-stable multishep methods for the Eq. (5), and the 2nd order of convergence is maximal for implicit multi-step methods. • The ZSP has a second order matrix, therefore, the inverse matrices and the matrix exponential can be easily calculated. This allows us to include practically any functions of the matrix Q in the difference schemes. • To calculate the spectral data, it is necessary to solve the ZSP for a large number of values ζ at a fixed potential q(t). This should be taken into account when implementing algorithms.

General theory of one-step schemes
Let us consider the problem (5) in a general case. We need to solve the equation where x = x(t) ∈ C n , using a one-step algorithm where T is the transition operator, x n = x(t n ), t n = nτ , τ is a step of the uniform grid.
We differentiate Eq. (15) and get the expressions for the derivatives D k x up to 5-th order as follows: Let us introduce the notation for the right-hand side of Eq. (17) and derivatives Using (17) and (18) we find recurrence relations for Q k as follows: Let us derive the Taylor series of x(t) at the point t, such that t n = t + sτ , t n+1 = t + sτ , s = s − 1: Then we denote the terms of (19) and (20): and write the expansion of (16) up to 5-th order After equating the terms of the same order we get Now we can derive recurrence relations for T k . From Eq. (23) we get Therefore, for first-order approximation, the expansion of transition operator T must begin with T ≈ E + τ Q and we need to know the value of Q at the point t.
From Eq. (24) we get To find T 2 we need to know the values of Q 2 and Q (1) at the point t. If we do not have the analytical expression of Q(t), we need to know the value of Q at two different points to use finite differences to calculate Q (1) . Otherwise, we can set s = 1/2 and zero the coefficient at Q 2 . Then we only need to know the value of Q 2 at the point t.
From Eq. (25) we get Equation has no real roots, therefore we can not get rid of the term with Q 3 varying s. It means that to use any scheme of order higher then 2 we must know Q (2) or values of Q at three different points.
From Eq. (26) we get Since equation has only one real root s = 1/2, this is the only way to get rid of Q 4 , that contains Q (3) .
From Eq. (27) we get Since equation has no real roots, we can not zero this coefficient varying s.

Thus we formulate
Theorem. Any one-step finite-difference scheme (16) approximates the equation (15) with a fourth order of accuracy if and only if the expansion of the transition operator T for the fixed s has a form where and the coefficients Q k are expressed through the matrix Q and its derivatives

Examples of schemes
Let us consider examples of constructing fourth-order schemes that satisfy the conditions of the theorem.

Constant matrix Q
Corollary 1. If the matrix Q is constant, then Q n = Q n and the expansion of the matrix T does not depend on s and has the form It is clear that Eq. (40) is the expansion of the matrix exponential exp(τ Q). Since it is an exact solution for the system with a constant matrix, then the one-step scheme with the exponential form has the infinity order of approximation for transition operator T = exp(τ Q).

Symmetrical case
As mentioned before, the second order scheme does not contain the derivative Q (1) if and only if s = 1/2. It means that the transition matrix T depends only on Q(t n + τ /2). Such choice of s corresponds to the center of the interval and will be called the symmetrical case. The expansion of the matrix T for the fourth order schemes in the symmetrical case gets rid of some terms and has only the dependence on Q (1) and Q (2) . Therefore, it is necessary to use Q at least at three points for the fourth order scheme.
Corollary 2. The expansion (33) of the matrix T in the symmetrical case has the form as follows: Approximating the derivatives in Eq. (41) by central finite differences of the second order we retain the fourth order of accuracy of the operator

Exponential form
The formula (41) is an expansion of the exponent where Replacing the derivatives by the finite differences (42) we get the exponential scheme with the fourth order of accuracy (ES4). The scheme conserves energy for the real spectral parameters ζ = ξ.
Exponential form (44) of the scheme allows approximating it by the Padé approximation and apply to multidimensional systems.

Magnus expansion
Exponential form of Eq. (44) follows from the Magnus expansion [10]. It provides an exponential representation of the exact evolution operator of the system (5) which is constructed as a series expansion referred to as Magnus expansion: Ω(t) = ∞ k=0 Ω k (t). First terms of this series have forms as follows: , Square brackets [A, B] = AB − BA are the matrix commutator of A and B.
If we represent a matrix Q(t) at the center of the integration interval t/2 by Taylor series and keep the main terms with respect to the small parameter t, then we get For t = τ formulas (46) coincide with (44)-(45). For the Schrödinger equation with a time-dependent operator, decomposition was obtained in a series of papers [11,12], but the authors decomposed the matrix Q in the Galerkin series and did not use difference schemes to find derivatives of Q.

Triple-exponential scheme
Equation (41) can be continued in a different way: We can continue the scheme (47) to the triple-exponential fourth-order scheme (TES4) which contains a spectral parameter ζ only at the exponential e τ Q . The exponential can be split [13], so the fast techniques (FNFT) can be applied to this scheme.
If Q is skew-Hermitian, then Q (1) and Q (2) are also skew-Hermitian. All exponentials in (48) are unitary matricies and, therefore, they preserve the quadratic invariant. The Maclaurin series of T (48) in τ gives exactly the decomposition of fourth order schemes for the symmetric case.

Numerical algorithms
We solve the system (5) on the uniform grid t n = −L + τ n with a step τ on the interval [−L, L], L = 30 unless otherwise stated. If the total number of points is 2M + 1, then the grid step is τ = L/M . We replace the original system (5) on the interval (t n − τ 2 , t n + τ 2 ) with an approximate system with constant coefficients The transition matrix T from the layer n − 1 2 to the layer n + 1 2 can be found using different numerical algorithms. Here we compared the numerical results for two new schemes presented above: the exponential scheme ES4 (44) and triple-exponential scheme TES4 (48). Then we tried the fourth-order conservative transformed scheme (CT4) with the transition operator T = e where M n+1 = e −τ Qn (Q n+1 − Q n ) e τ Qn , M n−1 = e τ Qn (Q n−1 − Q n ) e −τ Qn .
The CT4 scheme was introduced recently in [6]. Here we present new and more detailed numerical results for this scheme.
Among the well known algorithms we chose the Boffetta-Osborn second-order scheme (BO) [5] and the Runge-Kutta fourth-order algorithm (RK4). Following [14] for RK4 scheme we solve the system for the envelope χ 1 = φ 1 e iζt , χ 2 = φ 2 e −iζt . Unlike the above schemes, RK4 does not require computing the transition matrix T . Note also that the conventional Runge-Kutta algorithm uses half-steps in its description. Here we set this half-step equal τ , where τ is a grid step for the potential q(t).
The spectral data are finally defined by

Computation of the derivative of a(ζ)
To calculate the phase coefficients r k we need to find the derivatives where initial value is defined from (8) For the exponential scheme ES4 (44) the transition matrix T can be represented as T = exp(A) and calculated using Pauli matrices (72) from the appendix. Hence, we can find the derivative where , For triple-exponential scheme TES4 (48) we obtain the derivative of the transition matrix T where D = e τ Qn can be found from where σ 3 is a Pauli matrix (71).
For CT4 scheme the transition matrix T (50) can be represented as a fraction T = A −1 B, therefore .
We find the derivative of matrix exponential using (57).
For BO scheme the derivative of the transition matrix can be found in [5]. For RK4 scheme the derivative of a(ζ) can be computed using Romberg algorithm [14,15].

Model signals
We considered a model signal in the form of a chirped hyperbolic secant (59) For C = 0 it is a well-known Satsuma-Yajima signal. The detailed numerical results for this potential are presented in [4].
The analytical expressions of the spectral data of the potential (59) for anomalous dispersion are presented in [16]. But they can be obtained similarly for normal dispersion (σ = −1). Here we present general formulas using the Euler Gamma function Γ: The discrete spectrum ζ k , k = 0, K − 1 is determined by the zeros of the coefficient a(ζ) and exists only for anomalous dispersion (σ = 1): where square brackets denote the integer part of the expression.
To compute the phase coefficients r k we need to know the derivative a (ζ) only at points ζ k of the discrete spectrum. Let us write the coefficient a(ζ) in the following form: The function f (ζ) and its derivative f (ζ) have no singularities, so we have at the points of the discrete spectrum: where the function ϕ k is defined at the points of the discrete spectrum by a recurrence relation ϕ k+1 = −(k + 1)ϕ k , ϕ 0 = 1.
(63) Thus we have a formula to compute the phase coefficients To calculate energy of the discrete and continuous spectra E d , E c , we use the formula (13): Full energy of the potential (59) is easily computed as E = 2A 2 .
Let us denote K = A 2 − C 2 /4 + 1/2 as an integer part of the expression in square brackets and δ = A 2 − C 2 /4 + 1/2 as its fractional part. Then the discrete spectrum energy is and the continuous spectrum energy is (a) (c) (b) Figure 1: The approximation order with respect to the spectral parameter ξ.

Approximation order
The following formula was used to calculate the approximation order m: where τ i , i = 1, 2 are the steps of computational grids for two calculations with one spectral parameter ζ and τ 1 > τ 2 , Ψ i (L) is a deviation of the calculated value Ψ i (L) from the exact analytical valueΨ i (L) at the boundary point t = L.
The calculations were carried out for different p-norms and showed close values for the approximation orders. However, for the Euclidean 2-norm, the graphics were the smoothest.

Formulas for errors
We present the numerical errors of calculating the spectral data for continuous and discrete spectrum. To find the calculation errors of the continuous spectrum energy E c , residuals r k , and the coefficients a(ζ), b(ζ) at fixed ζ we use formula where φ can represent E c , r k , a(ζ) or b(ζ) at fixed ζ.
For the continuous spectrum we calculate the mean squared error where φ can represent a(ξ) or b(ξ). Here we suppose the spectral parameter ξ ∈ [−20, 20] with the total number of points N = 1025.   [6]. Actually, when calculating the continuous spectrum, it is necessary to choose a time step τ = L/M to describe correctly the fastest oscillations. For a fixed value of ξ, the local frequency ω(t; ξ) = ξ 2 + |q(t)| 2 of the system (5) varies from ω min = |ξ| to ω max = ξ 2 + q 2 max , where q max = max t |q(t)| is the maximum absolute value of the potential q(t). Therefore, step τ cannot be arbitrary. In order  to describe the most rapid oscillations, it is necessary to have at least 4-time steps for the oscillation period, so the inequality must be satisfied:

Numerical results for continuous spectrum
Therefore, any difference schemes will approximate the solutions of the original continuous system (5) if the inequality is fulfilled for the number of points M ≥ M min = 2 L ω max /π. Figure 2 also demonstrates a comparison of the computational time. One can see that ES4 scheme shows the best accuracy with a maximum speed. Figure 3 shows how the numerical schemes conserve energy. The numerical errors (68) for the continuous spectrum energy (14) are compared in Fig. 3 (a, c, e). To calculate the continuous spectrum energy it is important to define the size of the spectral domain L ξ and the corresponding grid step dξ [6]. According to the conventional discrete Fourier transform, we take the same number of points N ξ = N in the spectral domain and define a spectral step as dξ = π/(2L). So the size of the spectral interval is Figure 3 (b, d, f) demonstrates the deviation of the quadratic invariant H = |a| 2 + σ|b| 2 from unit with respect to the real spectral parameter ξ.
If the matrix Q is skew-Hermitian, then the matrix exp(τ Q) is unitary and the quadratic invariant conserves. For the direct ZSP this corresponds to anomalous dispersion (σ = 1) with a real spectral parameter ζ = ξ.
Let us consider a more general system with a matrix Q = KD, where K(t) is anti-Hermitian matrix depending on t, D is a constant Hermitian matrix. The system (5)  For one-step exponential methods Ψ n+ 1 2 = e RnD Ψ n− 1 2 , where R n is skew-Hermitian matrix, the quadratic invariant also conserves. It follows from the chain of equalities , DΨ n− 1 2 ). Here we used the formula De τ RnD = e τ DRn D, because for any natural p the equality is valid: From this result follows, that Boffetta-Osborn scheme (49) and the exponential scheme (44) are conservative for normal and anomalous dispersion σ = ±1. Similarly the scheme (50) also conserves the quadratic invariant, because it is the function of R n D.

Numerical results for discrete spectrum
Figures 4, 5 present the discrete spectrum errors (68). The parameters a(ζ k ), b(ζ k ) and r k were computed for the analytically known eigenvalues (61). Here we did not use any numerical algorithm to find the eigenvalue but computed spectral data at the exact point ζ = ζ k right away. It was made intentionally to estimate the error of the scheme itself and to avoid the influence of the other numerical algorithm errors.
There are well known problems with the computation of the coefficient b(ζ k ). We used the bi-directional algorithm [17] to find b(ζ k ) by the formula (11).

Conclusion
Two new forth-order exponential schemes for the numerical solution of the direct Zakharov-Shabat problem were presented and compared with the known ones. The ES4 scheme demonstrates the excellent computational speed and accuracy, but has difficulties with the direct application of the fast algorithms.
The TES4 scheme shows the same accuracy, but it requires about 2 times longer to calculate. However the main advantage of TES4 scheme is that fast algorithms can be applied to it.
The CT4 scheme also shows good accuracy comparable with two exponential schemes mentioned above, but it works about 2.5 − 3 times longer than ES4. The CT4 scheme does not allow the direct application of the fast algorithms. But it is possible after exponential approximation applied to the scheme.
All these schemes have an advantage over RK4 because they conserve the energy for the continuous spectrum parameter.