Novel Numerical Algorithm with Fourth-Order Accuracy for the Direct Zakharov-Shabat Problem

TL;DR

This paper introduces a new fourth-order numerical algorithm for solving the Zakharov-Shabat spectral problem, improving accuracy and efficiency over existing methods for continuous and discrete spectra.

Contribution

The paper presents a novel high-precision, fourth-order algorithm that generalizes the Boffetta-Osborne scheme for the Zakharov-Shabat problem.

Findings

Achieves fourth-order accuracy in solving the Zakharov-Shabat system

Enhances efficiency in spectral problem computations

Applicable to both continuous and discrete spectra

Abstract

We propose a new high-precision algorithm for solving the initial problem for the Zakharov-Shabat system. This method has the fourth order of accuracy and is a generalization of the second order Boffetta-Osborne scheme. It is allowed by our method to solve more effectively the Zakharov-Shabat spectral problem for continuous and discrete spectra.

Equations

$$\frac{d}{dt}\Psi (t)=Q(t)\Psi (t),$$

$${\Psi}(t)=\left(\begin{array}[]{c}\psi_{1}(t)\\ \psi_{2}(t)\end{array}\right),\quad Q(t)=\left(\begin{array}[]{cc}-i\zeta&q(t,z_{0})\\ -\sigma q^{*}(t,z_{0})&i\zeta\end{array}\right),$$

$$\left(\begin{array}[]{c}\psi_{1}\\ \psi_{2}\end{array}\right)=\left(\begin{array}[]{c}e^{-i\zeta t}\\ 0\end{array}\right)[1+o(1)],\quad t\to-\infty,$$

$$a(\xi )=\underset{t\to \infty }{\lim }{\psi }_1(t,\xi )e^{i\xi t},b(\xi )=\underset{t\to \infty }{\lim }{\psi }_2(t,\xi )e^{-i\xi t}.$$

$$r_k={\frac{b(\zeta )}{a^{\prime }(\zeta )}|}_{\zeta ={\zeta }_k},\text{mbox}where{a}^{\prime }(\zeta )=\frac{da(\zeta )}{d\zeta };$$

$$r(\xi )=\frac{b(\xi )}{a(\xi )},\xi \in \mathbb{R}.$$

$$\frac{d}{dt}\left(|{\psi }_1(t){|}^2+|{\psi }_2(t){|}^2\right)=0.$$

$$|{\psi }_1(t){|}^2+|{\psi }_2(t){|}^2=1.$$

$$C_n=-\frac{1}{\pi }\underset{-\infty }{\overset{\infty }{\int }}(2i\xi {)}^n\ln |a(\xi ){|}^{2}d\xi ++\sum _{k=1}^K\frac{1}{(n+1)}\left[(2i{\zeta }_k^{\ast }{)}^{n+1}-(2i{\zeta }_k{)}^{n+1}\right],$$

$$C_0=\underset{-\infty }{\overset{\infty }{\int }}|q{|}^{2}dt,C_1=\underset{-\infty }{\overset{\infty }{\int }}q{q}_t^{\ast }dt,C_2=\underset{-\infty }{\overset{\infty }{\int }}(q{q}_{tt}^{\ast }+|q{|}^4)dt,$$

$$C_3=\underset{-\infty }{\overset{\infty }{\int }}(q{q}_{ttt}^{\ast }+4|q{|}^{2}q{q}_t^{\ast }+|q{|}^2{q}^{\ast }{q}_t)dt.$$

$$\frac{d}{dt}\Psi (t)=Q(t_n)\Psi (t),Q_n=Q(t_n),$$

$${\Psi }_{n+\frac{1}{2}}=e^{\tau Q_n}{\Psi }_{n-\frac{1}{2}}.$$

$$\Psi (t)=e^{t{Q}_n}Y(t),$$

$$\frac{d}{dt}Y(t)=L(t)Y(t),L(t)=e^{-t{Q}_n}(Q(t)-Q_n)e^{t{Q}_n}.$$

$$\frac{d}{dt}Y(t)=0,$$

$$Y_{n+\frac{1}{2}}=Y_{n-\frac{1}{2}}.$$

$${\Psi }_{n-\frac{1}{2}}=e^{\left(t_n-\frac{\tau }{2}\right)Q_n}{Y}_{n-\frac{1}{2}},{\Psi }_{n+\frac{1}{2}}=e^{\left(t_n+\frac{\tau }{2}\right)Q_n}{Y}_{n+\frac{1}{2}},$$

$$\frac{Y_{n+\frac{1}{2}}-Y_{n-\frac{1}{2}}}{\tau }=L_n\frac{Y_{n+\frac{1}{2}}+Y_{n-\frac{1}{2}}}{2},$$

$$\frac{d}{dt}Y\left(t_n\right)-L(t_n)Y(t_n)\approx \frac{{\tau }^2}{24}\frac{d^3}{d{t}^3}Y(t_n).$$

$$\frac{Y_{n+\frac{1}{2}}-Y_{n-\frac{1}{2}}}{\tau }=\left(\alpha L_{n+1}+\beta L_{n-1}\right)Y_{n+\frac{1}{2}}+\left(\gamma L_{n+1}+\delta L_{n-1}\right)Y_{n-\frac{1}{2}}.$$

$$\frac{d}{dt}Y(t_n)-L(t_n)Y(t_n)\approx \frac{{\tau }^4}{24}(\beta -\alpha )\frac{d{L}_n}{dt}\frac{d^2{L}_n}{d{t}^2}{Y}_n-\frac{{\tau }^4}{5760}\left(17\frac{d^4{L}_n}{d{t}^4}+12\frac{d^2{L}_n}{d{t}^2}\frac{d{L}_n}{dt}\right)Y_n$$

$$\gamma =\frac{1}{24}-\alpha ,\delta =\frac{1}{24}-\beta$$

$$\left[I-\tau \alpha L_{n+1}-\tau \beta L_{n-1}\right]Y_{n+\frac{1}{2}}=\left[I+\tau \left(\frac{1}{24}-\alpha \right)L_{n+1}+\tau \left(\frac{1}{24}-\beta \right)L_{n-1}\right]Y_{n-\frac{1}{2}},$$

$$L_{n+1}=e^{-(t_n+\tau )Q_n}\left(Q_{n+1}-Q_n\right)e^{(t_n+\tau )Q_n},$$

$$L_{n-1}=e^{-(t_n-\tau )Q_n}\left(Q_{n-1}-Q_n\right)e^{(t_n-\tau )Q_n}.$$

$$\left[I-\tau \alpha M_{n+1}-\tau \beta M_{n-1}\right]e^{-\frac{\tau }{2}{Q}_n}{\Psi }_{n+\frac{1}{2}}=\left[I+\tau \left(\frac{1}{24}-\alpha \right)M_{n+1}+\tau \left(\frac{1}{24}-\beta \right)M_{n-1}\right]e^{\frac{\tau }{2}{Q}_n}{\Psi }_{n-\frac{1}{2}},$$

$$M_{n+1}=e^{-\tau Q_n}\left(Q_{n+1}-Q_n\right)e^{\tau Q_n},$$

$$M_{n-1}=e^{\tau Q_n}\left(Q_{n-1}-Q_n\right)e^{-\tau Q_n}.$$

$$T=e^{\frac{\tau }{2}{Q}_n}{\left[I-\frac{\tau }{48}\left(M_{n+1}+M_{n-1}\right)\right]}^{-1}\left[I+\frac{\tau }{48}\left(M_{n+1}+M_{n-1}\right)\right]e^{\frac{\tau }{2}{Q}_n}$$

Full text

K****eywords Zakharov-Shabat problem, inverse scattering transform, nonlinear Schrödinger equation, numerical methods

The solution of the direct problem for the Zakharov-Shabat problem (ZSP) is the first step in the inverse scattering transform (IST) for solving the nonlinear Schrödinger equation (NLSE) [1]. The numerical implementation of the IST has gained great importance and attracted special attention since Hasegawa and Tappert [2] proposed to use soliton solutions as a bit of information for fiber optic data transmission.

The direct scattering problem is solved by spectral data. To calculate them it is necessary to solve the initial problem for the Zakharov-Shabat system. Therefore, a lot of effort was made to find effective numerical methods for solving this problem. An overview of the methods used can be found in [3, 4, 5]. Currently, one of the most effective methods for solving the ZSP is the Boffetta-Osborne (BO) method [6], which has the second order of approximation. Comparisons of this method with other methods were carried out in [5, 7].

Besides the approximation accuracy, it is necessary to have an algorithm requiring a minimum computational time to get a discrete set of spectral parameters with sufficient accuracy. This direction is implemented in the fast algorithm (FNFT) for solving the direct ZSP using the modified Ablowitz-Ladik method [8, 4]. The BO method does not allow a direct application of the fast algorithm. But the fast method can be applied to the exponential approximation of the transition matrix of the BO method [9]. In this Letter, we will focus on building the method of the fourth order of accuracy on a uniform grid. For a non-uniform grid, a fourth order scheme [10] was applied in [11]. In perspective, the exponential approximation can be applied to our scheme so that we can use the fast algorithm.

We write the Zakharov-Shabat system in a matrix form

[TABLE]

where $\Psi(t)$ is a complex vector function of the real argument $t$, $Q(t)$ is a complex matrix

[TABLE]

where $\sigma=\pm 1$ for anomalous and normal dispersion, $z_{0}$ plays the role of a parameter and will not be used further. The asterisk means the complex conjugation.

Consider the Jost initial conditions

[TABLE]

which define the Jost solutions for real $\zeta=\xi$. The coefficients of the scattering matrix $a(\xi)$ and $b(\xi)$ are obtained as limits

[TABLE]

The function $a(\xi)$ can be extended to the upper half-plane $\xi\to\zeta$, where $\zeta$ is a complex number with the positive imaginary part $\eta=\mbox{Im}\,\zeta>0$. The spectral data are determined by $a(\zeta)$ and $b(\zeta)$ in the following way:

(1) the zeros of $a(\zeta)=0$ define the discrete spectrum $\{\zeta_{k}\}$, $k=1,...,K$ of ZSP (1) and phase coefficients

[TABLE]

(2) the continuous spectrum is determined by the reflection coefficient

[TABLE]

The matrix $Q(t)$ in the system (1) becomes the skew-Hermitian ($Q^{*}=-Q^{T}$) when the spectral parameter $\zeta=\xi$ is real and $\sigma=1$. Therefore, the system (1) preserves the integral

[TABLE]

Taking into account the boundary conditions (2), we have

[TABLE]

In addition, the trace formula is valid [12]

[TABLE]

which connects the NLSE integrals $C_{n}$ with the coefficient $a(\xi)$ and the discrete spectrum $\zeta_{k}$. The first integrals have the form

[TABLE]

[TABLE]

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.

We solve the system (1). The matrix $Q(t)$ linearly depends on the complex function $q(t)$, which is given in the whole nodes of the uniform grid $t_{n}=-L+\tau n$ with a step $\tau$ on the interval $[-L,L]$. If the total number of points is $2M+1$, then the grid step is $\tau=L/M$. Since the matrix $Q(t)$ is specified only on the grid, Boffetta and Osborne suggested replacing the original system on the interval $[t_{n}-\frac{\tau}{2},t_{n}+\frac{\tau}{2}]$ with an approximate system with constant coefficients [6]

[TABLE]

which is easily solved on the selected interval and gives the transition matrix from the layer $n-\frac{1}{2}$ to the layer $n+\frac{1}{2}$:

[TABLE]

This method has proven itself well, but nonetheless one would like to get a more accurate solution. So we formulate our task: it is required on the interval $[t_{n}-\frac{\tau}{2},t_{n}+\frac{\tau}{2}]$ to build the transition matrix from $\Psi_{n-\frac{1}{2}}$ to $\Psi_{n+\frac{1}{2}}$ with maximum accuracy and minimum computational cost.

The first step is a change of variables

[TABLE]

so the initial system takes the form

[TABLE]

In this form, the linear matrix $L(t)$ becomes zero at $t=t_{n}$, and the derivative of $Y(t)$ is zero at this point. This means that the solution is almost constant in the neighborhood of $t_{n}$. If the original system is replaced by an approximate

[TABLE]

then the transition from $Y_{n-\frac{1}{2}}$ to $Y_{n+\frac{1}{2}}$ becomes trivial:

[TABLE]

Returning to the original values of the variable $\Psi$ at the points $t_{n}-\frac{\tau}{2}$ and $t_{n}+\frac{\tau}{2}$, we get

[TABLE]

which, taking into account (10), exactly gives a transition in the BO scheme (7). Since we are interested in the values of $Y$ only in the grid nodes, the solution (10) can be interpreted as a solution to the difference equation

[TABLE]

which is an approximation of the continuous equation (8) given that $L_{n}=L(t_{n})=0$. By decomposing (12) into a Taylor series at the point $t=t_{n}$, we get the second order of approximation

[TABLE]

Thus, we have shown that the BO scheme corresponds to the simplest finite-difference approximation (12).

There are two possibilities to construct more complex approximations for the equation (8) on the interval $[t_{n}-\frac{\tau}{2},t_{n}+\frac{\tau}{2}]$. The first is to build a finite difference analog for this equation. The second possibility is to construct an approximation of the operator $L(t)$ on the entire interval $[t_{n}-\frac{\tau}{2},t_{n}+\frac{\tau}{2}]$ according to the existing values of $Q_{n}$ on a regular grid and the subsequent solution of such a system by any analytical method.

Consider the first approach. Since we want to refine the BO scheme, we take the function $Y$ only in two nodes of the grid $Y_{n-\frac{1}{2}}$ and $Y_{n+\frac{1}{2}}$. For the matrix $L$, we take the three nearest values $L_{n-1}$, $L_{n}$ and $L_{n+1}$. Using these values, we will look for a scheme using the method of uncertain coefficients

[TABLE]

Here we used the condition $L_{n}=0$, to drop the terms with $L_{n}$ in the right-hand side. By decomposing (Novel Numerical Algorithm with Fourth-Order Accuracy for the Direct Zakharov-Shabat Problem) into a Taylor series at the point $t=t_{n}$ and using the equation (8) at this point and its time derivatives, we get that the expression (Novel Numerical Algorithm with Fourth-Order Accuracy for the Direct Zakharov-Shabat Problem) has at least the fourth order approximation in $\tau$:

[TABLE]

for

[TABLE]

and arbitrary $\alpha$ and $\beta$. The resulting scheme can be rewritten as

[TABLE]

where

[TABLE]

[TABLE]

In the original variables, this scheme will take the form

[TABLE]

where

[TABLE]

[TABLE]

For real values $\zeta=\xi$, energy conservation (5) is important, so if $\alpha=\beta=1/48$, then the transition operator

[TABLE]

becomes a unitary matrix that conserves quadratic energy (5). Since the spectrum of matrices $Q_{n}$ is purely imaginary, the expression with square brackets is the Cayley formula. The transmission matrix (Novel Numerical Algorithm with Fourth-Order Accuracy for the Direct Zakharov-Shabat Problem) was obtained using a transformation of variables, and it conserves the energy for the real spectral parameters; therefore the corresponding scheme will be called as the fourth-order conservative transformed scheme (CT4).

Remark 1. The spectral parameter $\xi$ is included only through the exponent exponents, as in the BO method; therefore, the use of the fast algorithm (FNFT) is difficult [8], but it is possible after exponential approximation [9].

Remark 2. An open question is how to use the free parameters $a$ and $b$ for computation with complex spectral parameters. Although the preservation of high-frequency oscillations is important, for the eigenvalues near the imaginary axis, another criterion for the scheme may be needed.

The following formula was used to calculate the approximation order $m$:

[TABLE]

where $\tau_{i}$, $i=1$, $2$ are the steps of computational grids for two calculations with one spectral parameter $\zeta$ and $\tau_{1}>\tau_{2}$, $\tilde{\Psi}_{i}(L)$ is a deviation of the calculated value $\Psi_{i}(L)$ from the exact analytical value $\bar{\Psi}_{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.

The scheme (Novel Numerical Algorithm with Fourth-Order Accuracy for the Direct Zakharov-Shabat Problem) was tested for different model signals, where the analytical expressions for spectral data were known. In particular, there were calculations for the oversoliton from [13] for a small number of discrete eigenvalues. However, to present our scheme (Novel Numerical Algorithm with Fourth-Order Accuracy for the Direct Zakharov-Shabat Problem), we chose calculations for one soliton, because this solution is smooth and not only spectral data are known for it, but also eigenfunctions [14].

Here we present numerical results for the best known potential $q(t)=\mbox{sech}(t)$. It has a single eigenvalue $\zeta_{1}=0.5i$, $b(\zeta_{1})=-1$. Since this potential is purely solitonic $b(\xi)=0$, the continuous spectrum energy $E_{c}=-\frac{1}{\pi}\int\limits_{-\infty}^{\infty}\,\ln|a(\xi)|^{2}\,d\xi=0$, while $a(\xi)=(\xi-0.5i)/(\xi+0.5i)$.

Figure 1 demonstrates the approximation order $m$ of both schemes with respect to a spectral parameter $\xi\in[-20,20]$. Each line was calculated by the formula (13) using two embedded grids with a doubled grid step $\tau=L/M$, $L=40$. For the blue line, coarse and fine grids were defined by $M=2^{10}$ and $2^{11}$. For the red one, the grid was refined one more time, namely $M=2^{11}$ and $2^{12}$. Let us remind that the total number of points in the whole domain $[-L,L]$ is $2M+1$.

Figure 2 shows the continuous spectrum errors for the fixed value of the spectral parameter $\xi=20$. Black dashed line in Fig. 2 marks the minimum number of grid nodes $M_{\min}$ that guarantee a good approximation. Actually, when calculating the continuous spectrum, it is necessary to choose a time step $\tau=L/M$ to describe correctly the fastest oscillations. For a fixed value of $\xi$, the local frequency $\omega(t;\xi)=\sqrt{\xi^{2}+|q(t)|^{2}}$ of the system (1) varies from $\omega_{\min}=|\xi|$ to $\omega_{\max}=\sqrt{\xi^{2}+q_{\max}^{2}}$, where $q_{\max}=\max\limits_{t}|q(t)|$ is the maximum absolute value of the potential $q(t)$. Therefore, step $\tau$ 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:

[TABLE]

Therefore, any difference schemes will approximate the solutions of the original continuous system (1) if the inequality is fulfilled for the number of points $M\geq M_{\min}=2\,L\,\omega_{\max}/\pi$.

The calculation errors for the continuous spectrum energy are compared in Fig. 3. It is important to define the size of the spectral domain $L_{\xi}$ and the corresponding grid step $d\xi$ for the calculation of the continuous spectrum energy. According to the conventional discrete Fourier transform, we take the same number of points $N_{\xi}=N$ in the spectral domain and define a spectral step as $d\xi=\pi/(2L)$. So the size of the spectral interval is $L_{\xi}=\pi/(2\tau)$. The energy integral was computed by the trapezoid rule.

The discrete spectrum errors are presented in Fig. 4. The parameters $a(\zeta)$ and $b(\zeta)$ were computed for the analytically known eigenvalue $\zeta_{1}=0.5i$. In this test, we did not use any numerical algorithm to find the eigenvalue but compute $a(\zeta)$ and $b(\zeta)$ at the exact point $\zeta=\zeta_{1}$ 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.

All the errors in Figs.  2–4 are calculated using the Euclidean 2-norm.

Figures 2, 4 also demonstrate a comparison of the computational time. One can see that CT4 scheme allows getting a better accuracy faster than the BO scheme.

In this Letter, we proposed the family of fourth-order finite-difference one-step schemes for solving the direct Zakharov-Shabat problem on a uniform grid. Among this family, a quadratic integral preserving scheme for the continuous spectrum was distinguished. Numerical experiments for the soliton potential confirmed the theoretical order of approximation and demonstrated a significant advantage of our conservative scheme over the Boffetta-Osborne scheme. The proposed scheme works for uniform grids that can be useful when processing optical signals recorded at the receiver at regular time intervals.

Funding. This work was supported by the Russian Science Foundation (grant No. 17-72-30006).