# Numerical inverse scattering transform for the derivative nonlinear $\rm
Schr\ddot{o}dinger$ equation

Shikun Cui  Zhen Wang<EMAIL_ADDRESS>School of Mathematical Sciences,
Dalian University of Technology, Dalian, 116024, China  School of Mathematical
Sciences, Beihang University, Beijing, 100191, China

###### Abstract

In this paper, we develop the numerical inverse scattering transform (NIST)
for solving the derivative nonlinear $\rm Schr\ddot{o}dinger$ (DNLS) equation.
The key technique involves formulating a Riemann-Hilbert problem (RHP) that is
associated with the initial value problem and solving it numerically. Before
solving the RHP, two essential operations need to be carried out. Firstly,
high-precision numerical calculations are performed on the scattering data.
Secondly, the RHP is deformed using the Deift-Zhou nonlinear steepest descent
method. The DNLS equation has a continuous spectrum consisting of the real and
imaginary axes and features three saddle points, which introduces complexity
not encountered in previous NIST approaches. In our numerical inverse
scattering method, we divide the $(x,t)$-plane into three regions and propose
specific deformations for each region. These strategies not only help reduce
computational costs but also minimize errors in the calculations. Unlike
traditional numerical methods, the NIST does not rely on time-stepping to
compute the solution. Instead, it directly solves the associated Riemann-
Hilbert problem. This unique characteristic of the NIST eliminates convergence
issues typically encountered in other numerical approaches and proves to be
more effective, especially for long-time simulations.

###### keywords:

Inverse scattering transform, Riemann-Hilbert problem, Derivative nonlinear
$\rm Schr\ddot{o}dinger$ equation, Numerical method

††journal: XXX

2.5cm2.5cm1cm2cm

###### Contents

  1. 1 Introduction
  2. 2 Construction of the Riemann-Hilbert Problem
    1. 2.1 Spectral analysis
    2. 2.2 A Riemann-Hilbert problem associated with initial value problem
  3. 3 Numerical direct scattering
    1. 3.1 Calculating the discrete eigenvalues
    2. 3.2 Calculating the reflection coefficient
    3. 3.3 Calculating the normalized constant
  4. 4 Numerical inverse scattering
    1. 4.1 Deformations for the RHP (the non-soliton case)
      1. 4.1.1 Deformation in Region 1(-$c_{1}\sqrt{t}\leq x\leq c_{1}\sqrt{t}$)
      2. 4.1.2 Deformation in Region 2 ($x<-c_{1}\sqrt{t}$)
      3. 4.1.3 Deformation in Region 3($x>c_{1}\sqrt{t}$)
    2. 4.2 Removing the jump matrix $B$ (the non-soliton case)
      1. 4.2.1 Removing the jump matrix $B$ in Region 1
      2. 4.2.2 Removing the jump matrix $B$ in Region 2
      3. 4.2.3 Removing the jump matrix $B$ in Region 3
    3. 4.3 Deformations for the RHP (the soliton case)
      1. 4.3.1 Deformation in Region 1
      2. 4.3.2 Deformation in Region 2
      3. 4.3.3 Deformation in Region 3
    4. 4.4 Extra deformations for the RHP (soliton cases)
    5. 4.5 Numerical solution for the RHP
  5. 5 Numerical Results and Comparison
  6. 6 Conclusions and Discussions
  7. Acknowledgment
  8. References

##  1 Introduction

The nonlinear Schr$\rm\ddot{o}$dinger equation (NLS) is a significant
integrable equation that finds applications in describing the behavior of
optical solitons, Bose-Einstein condensation, and other physical phenomena[1,
2, 3, 4]. The derivative nonlinear Schr$\rm\ddot{o}$dinger (DNLS) equation is
an important variant of the NLS equation. Within the class of DNLS equations,
three well-known equations are the Kaup-Newell (KN) equation[5]

$iq_{t}+q_{xx}+i(|q|^{2}q)_{x}=0,$ (1)

the Chen-Lee-Liu (CLL) equation[6]

$iq_{t}+q_{xx}+i|q|^{2}q_{x}=0,$ (2)

and the Gerdzhikov-Ivanov (GI) equation[7]

$iq_{t}+q_{xx}-iq^{2}\bar{q}_{x}+\frac{1}{2}|q|^{4}q=0.$ (3)

. These equations have their own unique properties and are extensively studied
in various fields.

We refer to the KN equation, the CLL equation, and the GI equation as DNLSE
$\rm I$, DNLSE $\rm II$, and DNLSE $\rm III$, respectively. These equations
are interconnected through gauge transformations[8, 9, 10], allowing for
conversion between them. Therefore, our main focus lies in investigating the
specific properties and characteristics of DNLSE $\rm III$.

The inverse scattering transform (IST), originally proposed by GGKM, is a
significant approach for solving integrable equations with initial value
problems[11]. In the case of the DNLS equation, the inverse scattering
transform was introduced by Kaup and Newell[5]. A more contemporary method
that has gained widespread popularity is the Riemann-Hilbert method, which
serves as a modern version of the inverse scattering transform. By employing
the inverse scattering transform, Zhou and Huang were able to obtain
$N$-soliton solutions for the DNLS equation[12]. Additionally, Chen and Lam
extended the application of the inverse scattering transform to the DNLS
equation with non-zero boundary conditions[13]. Utilizing the Riemann-Hilbert
method, Guo and Ling derived the $N$-soliton solution for the coupled
derivative Schr$\rm\ddot{o}$dinger equation[14]. Zhang and Yan further
developed the inverse scattering transform for the DNLS equation, considering
both zero and non-zero boundary conditions[15]. Ma and Kuang expanded the
scope of the inverse scattering transform by applying it to the nonlocal DNLS
equation[16]. Lastly, Liu et al. successfully obtained triple-pole soliton
solutions of the DNLS equation using the inverse scattering transform[17].

The inverse scattering transform is not only applicable in solving the DNLS
equation but also in exploring the existence of solutions[22, 23, 24, 25].
Furthermore, it plays a crucial role in studying the long-term behavior of the
solution. Deift and Zhou introduced a nonlinear steepest descent method for
determining the long-term behavior of integrable equations[26]. Miller and
Bilman proposed the robust inverse scattering transform to analyze the
asymptotic properties of rogue waves[27]. Xu et al. investigated the long-term
asymptotics of the DNLS equation with step-like initial values[28]. Arruda et
al. studied the long-term asymptotics of the DNLS equation on the half-
line[29]. Jenkis et al. presented the soliton resolution for the derivative
nonlinear Schr$\rm\ddot{o}$dinger equation[30]. Tian examined the long-term
asymptotic behavior of the generalized derivative nonlinear
Schr$\rm\ddot{o}$dinger equation[31]. Yang and Fan analyzed the long-term
asymptotic behavior of the DNLS equation with initial data of finite density
type[32].

Methods for solving the DNLS equation can be classified into two categories:
analytical methods[5, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] and numerical
methods[34, 35, 36]. Analytical methods, such as the inverse scattering
transform, are primarily useful for calculating pure soliton solutions of the
DNLS equation. However, they fail to provide accurate results when the
reflection coefficient is non-zero, limiting their applicability to analyzing
the long-term asymptotic behavior of the solution[26, 29, 28, 30, 31, 32]. In
contrast, numerical methods offer a viable framework for solving general
initial value problems in the DNLS equation. Traditional numerical methods
rely on time-stepping to compute the solution, which accumulates calculation
errors over time and often leads to convergence issues. Hence, it becomes
crucial to develop a reliable numerical method capable of accurately
calculating the long-term evolution of the solution.

Recently, there has been a notable development in the field of solving
integrable equations with the introduction of the numerical inverse scattering
transform (NIST) by Deconinck, Olver, and Trogdon[38]. This powerful numerical
method has been successfully applied to various equations such as the focusing
and defocusing nonlinear Schr$\rm\ddot{o}$dinger equations by Trogdon and
Olver[39], the Toda lattice by Bilman and Trogdon[40], the sine-Gordon
equation by Deconinck, Trogdon, and Yang[41], and the Kundu-Eckhaus equation
by Cui and Wang [42]. The NIST can be divided into two main components:
numerical direct scattering and numerical inverse scattering. In the numerical
direct scattering step, numerical methods are employed to compute the
scattering data. Subsequently, in the numerical inverse scattering step, the
calculated scattering data is utilized to recover the solutions for the DNLS
equation. The NIST offers distinctive characteristics, which can be summarized
as follows:

  * 1.

Unlike traditional numerical methods, the NIST does not require time-stepping,
thus making the computational cost independent of time and space.

  * 2.

The error generated by the NIST remains unaffected by increasing time, thereby
avoiding numerical convergence issues.

  * 3.

While traditional numerical methods are more effective for short-time
simulations, the NIST proves to be more effective for long-time simulations.

  * 4.

It is important to note that currently, the NIST is primarily applicable to
integrable equations, and its implementation is complex.

The NIST represents a significant advancement in numerical methods for
integrable equations, demonstrating its advantages such as independence from
time-stepping, error stability, and increased effectiveness in long-time
simulations. However, it is worth mentioning that its application is confined
to integrable equations, and its implementation can be challenging.

In this paper, the numerical inverse scattering transform of the DNLS equation
is developed. In numerical direct scattering, we need to solve a quadratic
eigenvalue problem rather than a standard eigenvalue problem, which is
different from the AKNS system. Introducing $\widehat{\psi}=k\psi$, we
standardize the quadratic eigenvalue problem, then the high-precision
numerical method is used to solve the standard eigenvalue problem. The
continuous spectrum of the DNLS equation is $i\mathbb{R}\cup\mathbb{R}$ rather
than $\mathbb{R}$, which has never appeared in the previous NIST. Considering
that the original Riemann-Hilbert problem is oscillatory, we must deform it.
There are three saddle points for the DNLS equation, so deforming the original
Riemann-Hilbert problem is complex and difficult. Based on the relationship
between $x$ and $t$, we divide the $(x,t)$-plane into three regions and
propose three schemes for deforming the original Riemann-Hilbert problem.
These deformations are proposed after considering the computational cost and
accuracy comprehensively. At the same time, the rationality of the
deformations are demonstrated.

Our paper is structured as follows:

$\bullet$ Section 2: Construction of the Riemann-Hilbert Problem

We review the key findings related to the construction process of the Riemann-
Hilbert problem.

$\bullet$ Section 3: Numerical Direct Scattering

We implement numerical direct scattering methods. Through these techniques, we
calculate essential values such as the reflection coefficient, discrete
eigenvalues, and normalized constant for the DNLS equation.

$\bullet$ Section 4: Numerical Inverse Scattering

We employ numerical inverse scattering methods. In section 4.1 and section
4.2, we focus on deforming the Riemann-Hilbert problem, assuming no soliton
evolution from the initial value of the DNLS equation. In section 4.3 and
section 4.4, we examine the deformation of the Riemann-Hilbert problem,
considering scenarios where the initial value of the DNLS equation can evolve
into solitons. In section 4.5, we utilize the Chebyshev collocation method to
solve the deformed Riemann-Hilbert problem.

$\bullet$ Section 5: Numerical Results and Comparison

We present the numerical results obtained using the Numerical Inverse
Scattering Transform (NIST). Additionally, we compare the performance of the
NIST with the Fourier spectral method (FSM).

$\bullet$ Section 6: Conclusions and Discussions

We summarize the conclusions derived from our study. Furthermore, we engage in
discussions regarding the performance and effectiveness of the NIST in
comparison to the FSM. We also propose potential avenues for future research,
and emphasize the contributions our study makes.

##  2 Construction of the Riemann-Hilbert Problem

###  2.1 Spectral analysis

Considering the initial value problem for the DNLSE $\rm III$,

$\left\\{\begin{array}[]{lr}iq_{t}+q_{xx}-iq^{2}\bar{q}_{x}+\frac{1}{2}|q|^{4}q=0,\\\
q(x,0)=q_{0}(x),\end{array}\right.$ (4)

where $q_{0}(x)$ is a given function defined in Schwartz space
$\mathcal{S}(\mathbb{R})$, subscripts $x$ and $t$ represent the partial
derivatives with respect to space and time respectively, and $\bar{q}$ is the
complex conjugate of $q$.

The DNLS $\rm III$ is completely integrable and admits the following Lax
pair[24],

$\left\\{\begin{array}[]{lr}\psi_{x}+ik^{2}\sigma_{3}\psi=U\psi,\vspace{0.15cm}\\\
\psi_{t}+ik^{4}\sigma_{3}\psi=V\psi,\end{array}\right.$ (5)

where $k\in\mathbb{C}$ is the spectral parameter, $\psi$ is the eigenfunction,

$U=kQ+\frac{i}{2}|q|^{2}\sigma_{3},\
V=2k^{3}Q+ik^{2}|q|^{2}\sigma_{3}+\frac{i}{4}|u|^{4}\sigma_{3}-ikQ_{x}\sigma_{3}+\frac{1}{2}(q\bar{q}_{x}-q_{x}\bar{q})\sigma_{3},\
Q=\left(\begin{array}[]{cc}0&q\\\ -\bar{q}&0\\\ \end{array}\right).$

Here $\sigma_{3}$ is one of the Pauli matrices

$\sigma_{1}=\left(\begin{array}[]{cc}0&1\\\ 1&0\\\
\end{array}\right),\sigma_{2}=\left(\begin{array}[]{cc}0&i\\\ -i&0\\\
\end{array}\right)\ {\rm and}\ \sigma_{3}=\left(\begin{array}[]{cc}1&0\\\
0&-1\\\ \end{array}\right).$

The eigenfunction $\psi$ admits the following asymptotics

$\psi\rightarrow{\rm e}^{-{i}\theta\sigma_{3}}\ {\rm as}\
|x|\rightarrow\infty,{\rm where}\ \theta=k^{2}x+2k^{4}t.$

Doing transformation $\mu=\psi{\rm e}^{i\theta\sigma_{3}}$, we get a new Lax
pair

$\left\\{\begin{array}[]{lr}\mu_{x}+{i}k^{2}[\sigma_{3},\mu]=U\mu,\vspace{0.15cm}\\\
\mu_{t}+2{i}k^{4}[\sigma_{\rm 3},\mu]=V\mu,\end{array}\right.$ (6)

where $[\sigma_{3},\mu]=\sigma_{3}\mu-\mu\sigma_{3}$.

The eigenfunction $\mu$ admits the following asymptotics

$\mu\rightarrow I\ {\rm as}\ |x|\rightarrow\infty,$

where $I$ is the identity matrix.

######  Proposition 2.1.

The solution of the DNLSE $\rm III$ is given by

$q=2{i}\mu^{(1)}_{(1,2)},$ (7)

where $\mu^{(1)}_{(1,2)}$ represents the first row and second column element
of matrix $\mu^{(1)}$, $\mu^{(1)}=\lim_{k\rightarrow\infty}(k\mu).$

###### Proof.

Expand $\mu$ in $k=\infty$, and we get

$\mu=\mu^{(0)}+\frac{\mu^{(1)}}{k}+\frac{\mu^{(2)}}{k^{2}}+\cdots$ (8)

Substitute Eq.(8) into Eq.(6), we can prove Eq.(7) by matting the $O(k^{0})$,
$O(k^{1})$ and $O(k^{2})$ terms.

∎

The Lax pair (6) can be written into a total differential form,

${\rm d}({\rm e}^{i\theta\hat{\sigma}_{3}}\mu)={\rm
e}^{i\theta\hat{\sigma}_{3}}(U{\rm d}x+V{\rm d}t)\mu,$

where ${\rm e}^{i\theta\hat{\sigma}_{3}}\mu={\rm e}^{i\theta\sigma_{3}}\mu{\rm
e}^{-{i}\theta\sigma_{3}}.$

The eigenfunctions of the Lax pair (6) are obtained, which satisfy the
Volterra integral equation,

$\begin{array}[]{lr}\mu_{1}(x,t,k)=I+\int_{-\infty}^{x}{\rm
e}^{{-i}k^{2}(x-y)\hat{\sigma}_{3}}U(y,t,k)\mu_{1}(y,t,k)dy,\\\
\mu_{2}(x,t,k)=I+\int_{+\infty}^{x}{\rm
e}^{{-i}k^{2}(x-y)\hat{\sigma}_{3}}U(y,t,k)\mu_{2}(y,t,k)dy.\end{array}$ (9)

The eigenfunctions of Eq. (9) exhibit specific symmetries and analyticity
properties. These characteristics play a crucial role in constructing the
Riemann-Hilbert problem associated with the initial value problem described by
equation (4). These symmetries and analyticity properties allow for the
formulation of a well-posed mathematical problem known as the Riemann-Hilbert
problem. The Riemann-Hilbert problem involves finding a function that
satisfies certain conditions on the boundary of a given domain. By leveraging
the symmetries and analyticity of the eigenfunctions, we can successfully
construct the Riemann-Hilbert problem and subsequently analyze the associated
initial value problem (4). This enables a deeper understanding of the dynamics
and behavior of the system under investigation.

######  Proposition 2.2 (Analyticity).

The eigenfunction $\mu_{j}(x,t,k)$ $(j=1,2)$ satisfies the following
analyticity,

  * 1.

$[\mu_{1}]_{2}$ and $[\mu_{2}]_{1}$ are bounded and analytic in area
$D_{1}$=$\\{k\in\mathbb{C}|\ {\rm Im}(k^{2})<0\\},$

  * 2.

$[\mu_{1}]_{1}$ and $[\mu_{2}]_{2}$ are bounded and analytic in area
$D_{2}$=$\\{k\in\mathbb{C}|\ {\rm Im}(k^{2})>0\\},$

  * 3.

${\rm det}[\mu_{1}(x,t,k)]={\rm det}[\mu_{2}(x,t,k)]=1,$

where $[\mu_{j}]_{k}$ is the $k$-th column of $\mu_{j}(x,t,k)$. The area
$D_{1}$ and the area $D_{2}$ are shown in Figure 1.

Figure 1: The area $D_{1}$ and area $D_{2}$ on complex-$k$ plane.

######  Proposition 2.3 (Symmetry).

The eigenfunction $\mu_{j}(x,t,k)$ $(j=1,2)$ satisfies the following
symmetries,

  * 1.

$-\sigma\overline{\mu_{j}(x,t,\bar{k})}\sigma=\mu_{j}$,

  * 2.

$\sigma_{3}\mu_{j}(x,t,-k)\sigma_{3}=\mu_{j}$,

where $\sigma$= $\left(\begin{array}[]{cc}0&1\\\ -1&0\\\ \end{array}\right)$,
$\sigma_{\rm 3}$= $\left(\begin{array}[]{cc}1&0\\\ 0&-1\\\
\end{array}\right)$.

###### Proof.

Considering the $x$-part of the Lax pair (6)

$(\ \mu_{j}(k)\ )_{x}+ik^{2}[\sigma_{3},\mu_{j}(k)]=(kQ+q_{1})\mu_{j}(k),$
(10)

where $q_{1}=\begin{cases}0\sigma_{3},&{\rm for\ the\ DNLS\ I},\\\
-\frac{i}{4}|q|^{2}\sigma_{3},&{\rm for\ the\ DNLS\ II},\\\
\frac{i}{2}|q|^{2}\sigma_{3},&{\rm for\ the\ DNLS\ III}.\end{cases}$

When we replace $k$ with $\bar{k}$ and perform a conjugate transformation of
Eq. (10), the new equation becomes:

$(\ \overline{\mu_{j}(\bar{k})}\
)_{x}-{i}k^{2}[\sigma_{3},\overline{\mu_{j}(\bar{k})}]=(k\overline{Q}-q_{1})\overline{\mu_{j}(\bar{k})}.$
(11)

Multiply $\sigma$ on both sides of Eq. (11), and we get

$(\ \sigma\overline{\mu_{j}(\bar{k})}\sigma\
)_{x}+{i}k^{2}[\sigma_{3},\sigma\overline{\mu_{j}(\bar{k})}\sigma]=(kQ+q_{1})\sigma\overline{\mu_{j}(\bar{k})}\sigma.$
(12)

Both $-\sigma\overline{\mu_{j}(\bar{k})}\sigma$ and $\mu_{j}(k)$ are the
solution of Eq. (10), and they have the same asymptotics,

${\rm so}\quad-\sigma\overline{\mu_{j}(x,t,\bar{k})}\sigma=\mu_{j}(x,t,k),\
(j=1,2).$

Replace $k$ with $-k$ in Eq. (10), and we get

$(\ \mu_{j}(-k)\ )_{x}+ik^{2}[\sigma_{3},\mu_{j}(-k)]=(-kQ+q_{1})\mu_{j}(-k).$
(13)

Multiplying $\sigma_{3}$ on both sides of Eq. (13), we get

$(\sigma_{3}\ \mu_{j}(-k)\sigma_{3}\
)_{x}+ik^{2}[\sigma_{3},\sigma_{3}\mu_{j}(-k)\sigma_{3}]=(kQ+q_{1})\
\sigma_{3}\mu_{j}(-k)\sigma_{3}.$

Both $\sigma_{3}\mu_{j}(-k)\sigma_{3}$ and $\mu_{j}(k)$ are the solution of
Eq. (10), and they have the same asymptotics

$\sigma_{3}\mu_{j}(-k)\sigma_{3}=\mu_{j}(k),\ (j=1,2).$

∎

Both $\mu_{1}(x,t,k)$ and $\mu_{2}(x,t,k)$ are the solution of the Lax pair
(6), so there exists a linear relation between $\mu_{1}(x,t,k)$ and
$\mu_{2}(x,t,k)$,

$\mu_{2}(x,t,k)=\mu_{1}(x,t,k){\rm e}^{-i\theta\hat{\sigma}_{3}}S(k),$ (14)

where $S(k)$ is independent of $t$ and $x$.

######  Proposition 2.4.

$S(k)$ satisfies the following symmetries:
$S(k)=-\sigma\overline{S(\bar{k})}\sigma$, $S(k)=\sigma_{\rm
3}S(-k)\sigma_{3},$

$S(k)=\left(\begin{array}[]{cc}\tilde{a}(k)&b(k)\\\ \tilde{b}(k)&a(k)\\\
\end{array}\right),$ $\tilde{a}(k)=\overline{a(\bar{k})},$
$\tilde{b}(k)=-\overline{b(\bar{k})}.$

###### Proof.

With the help of $\mathbf{Proposition\ 2.3}$, the symmetries of $S(k)$ can be
proved easily.

∎

###  2.2 A Riemann-Hilbert problem associated with initial value problem

Rewriting Eq. (14) into

$([\mu_{2}]_{1},[\mu_{2}]_{2})=([\mu_{1}]_{1},[\mu_{1}]_{2}){\rm
e}^{-i\theta\hat{\sigma}_{3}}\left(\begin{array}[]{cc}\hat{a}(k)&b(k)\\\
\hat{b}(k)&a(k)\\\ \end{array}\right),$ (15)

where $[\mu_{1}]_{2}$ and $[\mu_{2}]_{1}$ are bounded and analytic in $D_{1}$,
$[\mu_{1}]_{1}$ and $[\mu_{2}]_{2}$ are bounded and analytic in $D_{2}$. The
area $D_{1}$ and $D_{2}$ are shown in Figure 1.

We define $M(x,t,k)$ by

$M(x,t,k)=\left\\{\begin{array}[]{lr}\Big{(}\frac{[\mu_{2}]_{1}}{\tilde{a}(k)},[\mu_{1}]_{2},\Big{)},\quad
k\in D_{1},\\\ \Big{(}[\mu_{1}]_{1},\frac{[\mu_{2}]_{2}}{a(k)}\Big{)},\quad
k\in D_{2},\end{array}\right.$

then we get a Riemann-Hilbert problem (RHP):

Riemann-Hilbert problem

  * 1.

$M(x,t,k)$ is meromorphic in $\mathbb{C}\setminus\Sigma_{1}$,

  * 2.

$M_{+}(x,t,k)=M_{-}(x,t,k)T(x,t,k)$, $k\in\Sigma_{1}$,

  * 3.

$M(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where $\Sigma_{1}=i\mathbb{R}\cup\mathbb{R}$ , the jump matrix

$T(x,t,k)=\left(\begin{array}[]{cc}1&{r}(k){\rm e}^{-2i\theta(k)}\\\
-\tilde{r}(k){\rm e}^{2i\theta(k)}&1-r(k)\tilde{r}(k)\\\ \end{array}\right),$

where $r(k)=\frac{b(k)}{a(k)}$ and
$\tilde{r}(k)=\frac{\tilde{b}(k)}{\tilde{a}(k)}$.

######  Remark 1.

If $a(k)$ exists the one-order zero point $\kappa_{j}$ in the first quadrant
of $\mathbb{C}$, Riemann-Hilbert problem will satisfy the pole conditions,

  * 1.

${\rm
Res}\\{M(x,t,k),k=\kappa_{j}\\}=\lim_{k\rightarrow\kappa_{j}}M(x,t,k)\left(\begin{array}[]{cc}0&C_{j}{\rm
e}^{-2i\theta(\kappa_{j})}\\\ 0&0\\\ \end{array}\right),$

  * 2.

Res{M(x,t,k), $k=\bar{\kappa}_{j}$} =
$\lim_{k\rightarrow\bar{\kappa}_{j}}M(x,t,k)\left(\begin{array}[]{cc}0&0\\\
\bar{C}_{j}{\rm e}^{2{i}\theta(\bar{\kappa}_{j})}&0\\\ \end{array}\right),$

  * 3.

Res{M(x,t,k), $k=-\bar{\kappa}_{j}$} =
$\lim_{k\rightarrow-\bar{\kappa}_{j}}M(x,t,k)\left(\begin{array}[]{cc}0&0\\\
\bar{C}_{j}{\rm e}^{2i\theta(-\bar{\kappa}_{j})}&0\\\ \end{array}\right),$

  * 4.

Res{M(x,t,k), $k=-\kappa_{j}$} =
$\lim_{k\rightarrow-\kappa_{j}}M(x,t,k)\left(\begin{array}[]{cc}0&C_{j}{\rm
e}^{-2i\theta(-\kappa_{j})}\\\ 0&0\\\ \end{array}\right).$

######  Proposition 2.5.

The solution of the $\rm DNLSE\ I$ $\rm(\ref{initial_problem})$ is given by
the Riemann-Hilbert problem

$q(x,t)=2im_{(12)},$

where $m_{(12)}$ represents the first row and second column element of matrix
$m$, $m=\lim_{k\rightarrow\infty}\Big{(}kM(x,t,k)\Big{)}$.

##  3 Numerical direct scattering

In this section, the numerical methods are used to calculate the scattering
data. The scattering data consists of the reflection coefficient $r(k)$, the
discrete eigenvalue $\kappa_{j}(j=1,2,...,n)$ and the normalized constant
$C_{j}$. It is important to note that the condition $a(\kappa_{j})=0$ is
equivalent to $\kappa_{j}$ being a discrete eigenvalue of the DNLS equation.
If $k=\kappa_{j}$ is the one-order pole of $a(k)$, we define the normalized
constant by

$C_{j}=\frac{b_{j}}{a^{\prime}(\kappa_{j})},$

where $a^{\prime}(k)$ is the derivative of $a(k)$ and $b_{j}$ is a constant.
Set $k=\kappa_{j}$, Eq. (15) gives

$[\mu_{2}]_{2}(x,t,\kappa_{j})=b_{j}{\rm
e}^{-2i\theta(\kappa_{j})}[\mu_{1}]_{1}(x,t,\kappa_{j}).$

###  3.1 Calculating the discrete eigenvalues

There are two main methods available for calculating the discrete eigenvalues
of the DNLS equation[33]. The first method involves an iterative approach to
find the zero points of $a(k)$, while the second method revolves around
solving a matrix eigenvalue problem. When $\kappa_{j}$ is identified as a
discrete eigenvalue for the DNLS equation, $\mathbf{Proposition\ 2.3}$ asserts
that $\bar{\kappa}_{j}$, $-\bar{\kappa}_{j}$, and $-\kappa_{j}$ also qualify
as discrete eigenvalues. This property simplifies the search for zero points
of a(k) to the first quadrant of the complex plane. However, even limiting the
search to the first quadrant of the complex-$k$ plane can still be complex and
computationally expensive when using the iterative method. Therefore, the more
efficient approach of solving the matrix eigenvalue problem is typically
employed to compute the eigenvalues for the DNLS equation.

The eigenvalues for the DNLS equation are calculated from the following
eigenvalue problems

$\displaystyle\psi_{x}+ik^{2}\sigma_{3}\psi=kQ\psi,\quad$ for the DNLSE I,
(16)
$\displaystyle\psi_{x}+ik^{2}\sigma_{3}\psi=kQ\psi+\frac{i}{4}Q^{2}\sigma_{3}\psi,\quad$
for the DNLSE II, (17)
$\displaystyle\psi_{x}+ik^{2}\sigma_{3}\psi=kQ\psi-\frac{i}{2}Q^{2}\sigma_{3}\psi,\quad$
for the DNLSE III . (18)

Eqs. (16), (17) and (18) can be classified as quadratic eigenvalue problems,
which we refer to as DNLSE eigenvalue problems. In order to solve these
problems more effectively, we must convert the DNLSE eigenvalue problems into
standard eigenvalue problems.

######  Proposition 3.1.

The DNLSE eigenvalue problems (16), (17) and (18) can be written as the
standard eigenvalue problems by introducing transformation
$\widehat{\psi}=k\psi$.

###### Proof.

We only take the eigenvalue problem (16) as an example, the eigenvalue
problems (17) and (18) are similar to (16). Set $\widehat{\psi}=k\psi$, we get

$\displaystyle\psi_{x}+{i}k\sigma_{3}\widehat{\psi}=Q\widehat{\psi},$ (19)
$\displaystyle\widehat{\psi}=k\psi,$

where $\psi$ and $\widehat{\psi}$ are two vectors of size $2\times 1$.

Eq. (LABEL:linear_1) is rewritten into

$\left(\begin{array}[]{cccc}0&-iq&i\partial_{x}&0\\\
-i\bar{q}&0&0&-i\partial_{x}\\\ 1&0&0&0\\\ 0&1&0&0\\\
\end{array}\right)\left(\begin{array}[]{c}\widehat{\psi}\\\ \psi\\\
\end{array}\right)=k\left(\begin{array}[]{c}\widehat{\psi}\\\ \psi\\\
\end{array}\right),$ (20)

where $q$ is the potential function, $\partial_{x}$ represents the
differential operation, $\widehat{\psi}$ and $\psi$ are the eigenfunctions. It
is obvious that Eq. (20) is a standard eigenvalue problem.

∎

Proposition 3.1 gives a method to standardize the DNLSE eigenvalue problems.
The $\rm DNLSE\ {I}$ eigenvalue problem is transformed into Eq. (20). The $\rm
DNLSE\ {II}$ eigenvalue problem is transformed into

$\left(\begin{array}[]{cccc}0&-{i}q&i\partial_{x}-\frac{1}{4}|q|^{2}&0\\\
-i\bar{q}&0&0&-i\partial_{x}-\frac{1}{4}|q|^{2}\\\ 1&0&0&0\\\ 0&1&0&0\\\
\end{array}\right)\left(\begin{array}[]{c}\widehat{\psi}\\\ \psi\\\
\end{array}\right)=k\left(\begin{array}[]{c}\widehat{\psi}\\\ \psi\\\
\end{array}\right).$ (21)

The $\rm DNLSE\ {III}$ eigenvalue problem is transformed into

$\left(\begin{array}[]{cccc}0&-iq&i\partial_{x}+\frac{1}{2}|q|^{2}&0\\\
-i\bar{q}&0&0&-i\partial_{x}+\frac{1}{2}|q|^{2}\\\ 1&0&0&0\\\ 0&1&0&0\\\
\end{array}\right)\left(\begin{array}[]{c}\widehat{\psi}\\\ \psi\\\
\end{array}\right)=k\left(\begin{array}[]{c}\widehat{\psi}\\\ \psi\\\
\end{array}\right).$ (22)

By transforming the quadratic eigenvalue problems into standard eigenvalue
problems, we simplify the solution process and can leverage existing numerical
methods for solving standard eigenvalue problems to find the eigenvalues and
eigenvectors of the DNLSE equations. There are some methods can be used to
solve these eigenvalue problems[34, 35, 37, 43], we use the Chebyshev
collocation method[43] to solve the DNLSE eigenvalue problems. Before
introducing the numerical method, some marks are introduced in Table 1.

Table 1: Annotations of some marks used in this paper. Mark | Formula of the mark | Note
---|---|---
$\vec{\chi}$ | $\vec{\chi}=\bigg{(}-1,{\rm cos}(\frac{n_{1}-2}{n_{1}-1}\pi),\ \cdots\,\ {\rm cos}(\frac{1}{n_{1}-1}\pi),1\bigg{)}^{\top}$ | Chebyshev nodes
$n_{1}$ |  | The number of Chebyshev nodes
$T_{k}(x)$ |  $T_{k}={\rm cos}(k\cdot{\rm arccos}x)$, $k=0,\ldots,n_{1}-1$ | Chebyshev polynomials
$T(x)$ | $T(x)=[T_{0}(x),T_{1}(x),\ldots,T_{n_{1}-1}(x)]$ |
$H(x)$ |  $H(x)={\rm tanh}(\hat{a}x)$, $0<\hat{a}<1$ | Hyperbolic tangent map
$H^{-1}(x)$ |  | The inverse mapping of $H(x)$
$\mathcal{D}$ |  | Chebyshev derivative matrix
$f(\vec{\chi})$ | $f(\vec{\chi})=\left(\begin{array}[]{c}f(x)\mid_{x=-1}\\\ f(x)\mid_{x={\rm cos}(\frac{n_{1}-2}{n_{1}-1}\pi)}\\\ \vdots\\\ f(x)\mid_{x=1}\\\ \end{array}\right)$ |
$\mathcal{F}$ | The discrete cosine transform of $f(\vec{\chi})$ | Chebyshev coefficient matrix
$\mathrm{D}\Big{[}f(\vec{\chi})\Big{]}$ | ${\rm D}\Big{[}f(\vec{\chi})\Big{]}=\left(\begin{array}[]{ccc}f(x)\mid_{x=-1}&&\\\ &\ddots&\\\ &&f(x)\mid_{x=1}\\\ \end{array}\right)$ | Diagonal matrix

We need to transform the DNLSE eigenvalue problems into matrix eigenvalue
problems. For a given function $g(x)$ defined in real field $\mathbb{R}$, we
can approximate $g(x)$ and $\frac{\partial g}{\partial x}$ by

$\displaystyle g(x)=T^{\mathbb{R}}(x)\mathcal{F}g(\hat{\chi}),x\in\mathbb{R},$
(23) $\displaystyle\frac{\partial g(x)}{\partial x}=\frac{\partial
H(x)}{\partial x}T^{\mathbb{R}}(x)\mathcal{D}\mathcal{F}g(\hat{\chi}),{\rm
for}\ x\in\mathbb{R},$ (24)

where $T^{\mathbb{R}}(x)=T(H(x))=[T_{0}(H(x)),\cdots,T_{n_{1}-1}(H(x))]$ and
$\hat{\chi}=H^{-1}(\vec{\chi})$.

Using Eq. (23) and Eq. (24), $\psi=\left(\begin{array}[]{c}\psi_{1}\\\
\psi_{2}\\\ \end{array}\right)$,
$\widehat{\psi}=\left(\begin{array}[]{c}\widehat{\psi_{1}}\\\
\widehat{\psi_{2}}\\\ \end{array}\right)$, $\psi_{x}$ and the potential
function $q(x)$ are approximated by the Chebyshev polynomials with $n_{1}$
nodes,

$\displaystyle\psi_{j}=T^{\mathbb{R}}(x)\mathcal{F}\psi_{j}(\hat{\chi}),$ (25)
$\displaystyle\widehat{\psi_{j}}=T^{\mathbb{R}}(x)\mathcal{F}\widehat{\psi_{j}}(\hat{\chi}),$
$\displaystyle\frac{\partial\psi_{j}}{\partial x}=\frac{\partial
H(x)}{\partial x}T^{\mathbb{R}}(x)\mathcal{D}\mathcal{F}\psi_{j}(\hat{\chi}),$
$\displaystyle q(x)=T^{\mathbb{R}}(x)\mathcal{F}q(\hat{\chi}).$

Substitute Eq. (25) into Eqs. (20), (21) and (22), and set
$x=\hat{\chi}=H^{-1}(\vec{\chi})$, then the DNLSE eigenvalue problem is
transformed into a matrix eigenvalue problem

$\left[\begin{array}[]{cccc}0&i{\rm D}\Big{[}q(\hat{\chi})\Big{]}&\quad{\rm
D}\Big{[}q_{1}(\hat{\chi})\Big{]}+i{\rm D}\Big{[}\frac{\partial
H(\vec{\chi})}{\partial x}\Big{]}\mathcal{F}^{-1}\mathcal{D}\mathcal{F}&0\\\
i{\rm D}[\bar{q}(\hat{\chi})]&0&0&{\rm D}\Big{[}q_{1}(\hat{\chi})\Big{]}-i{\rm
D}\Big{[}\frac{\partial H(\vec{\chi})}{\partial
x}\Big{]}\mathcal{F}^{-1}\mathcal{D}\mathcal{F}\\\ I&0&0&0\\\
0&I&0&0\end{array}\right]\left[\begin{array}[]{c}\widehat{\psi_{1}}(\hat{\chi})\\\
\widehat{\psi_{2}}(\hat{\chi})\\\ \psi_{1}(\hat{\chi})\\\
\psi_{2}(\hat{\chi})\\\
\end{array}\right]=k\left[\begin{array}[]{c}\widehat{\psi_{1}}(\hat{\chi})\\\
\widehat{\psi_{2}}(\hat{\chi})\\\ \psi_{1}(\hat{\chi})\\\
\psi_{2}(\hat{\chi})\\\ \end{array}\right],$ (26)

where $q_{1}=\begin{cases}{}0&{\rm for\ the\ DNLSE\ {I}},\\\
-\frac{1}{4}|q|^{2}&{\rm for\ the\ DNLSE\ {II}},\\\ \frac{1}{2}|q|^{2}&{\rm
for\ the\ DNLSE\ {III}}.\end{cases}$

The shape of the matrix eigenvalue problems (26) is $4n_{1}\times 4n_{1}$, and
we can solve it by the ‘eig’ function in MATLAB.

For the DNLSE $I$, Kaup and Newell got the soliton solutions by the inverse
scattering transform[5]. The function

$q_{KN}(x,0)=4\frac{{\exp}(2x)\Big{(}{\exp}(4x)+{\exp}(-\frac{\pi
i}{2})\Big{)}}{\Big{(}{\exp}(4x)+{\exp}(\frac{\pi i}{2})\Big{)}^{2}}$ (27)

is the initial value for a pure soliton solution. There are four discrete
eigenvalues $k_{1}=\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}{i}$,
$\bar{k}_{1}=\frac{\sqrt{2}}{2}-\frac{\sqrt{2}}{2}{i}$,
$-\bar{k}_{1}=-\frac{\sqrt{2}}{2}+\frac{\sqrt{2}}{2}{i}$ and
$-k_{1}=-\frac{\sqrt{2}}{2}-\frac{\sqrt{2}}{2}{i}$ for (27). We calculate the
eigenvalues of (27) for the DNLSE $\rm I$ under $\hat{a}=0.1$ and $n_{1}=600$,
the calculated results are shown in Figure 2. The absolute error between the
calculated result and the exact result is $4.80\times 10^{-15}$ at $k=k_{1}$,
the absolute error is $1.94\times 10^{-15}$ at $k=\bar{k}_{1}$, the absolute
error is $9.89\times 10^{-15}$ at $k=-\bar{k}_{1}$, and the absolute error is
$6.34\times 10^{-15}$ at $k=-k_{1}$.

Figure 2: The calculated eigenvalues of the DNLSE $\rm I$ with (27) potential.

For the DNLSE $\rm III$, Fan got the soliton solutions by the Darboux
transformation[18]. The DNLSE $\rm III$ exists the follow one-solution
solution

$q_{GI}(x,t)=-8\frac{\xi_{1}\eta_{1}}{(\xi_{1}+i\eta_{1}){\rm
e}^{(X_{1}+iY_{1})}+(\xi_{1}-i\eta_{1}){\rm e}^{(-X_{1}+iY_{1})}},$ (28)

where $\xi_{1}$ is the real part of the discrete eigenvalue $k_{2}$ and
$\eta_{1}$ is the imaginary part of the discrete eigenvalue $k_{2}$,

$X_{1}=-4\xi_{1}\eta_{1}x-16\xi_{1}\mu_{1}(\xi_{1}^{2}-\eta_{1}^{2})t,$
$Y_{1}=2(\xi_{1}^{2}-\eta_{1}^{2})x+4(\xi_{1}^{4}+\eta_{1}^{4}-6\xi_{1}^{2}\eta_{1}^{2})t.$

There are four discrete eigenvalues $k_{2}=1+0.5i$, $\bar{k}_{2}=1-0.5i$,
$-\bar{k}_{2}=-1+0.5i$ and $-k_{2}=-1-0.5i$ for

$q_{GI}(x,0)=\frac{-4}{(1+0.5i){\rm e}^{-2x+1.5ix}+(1-0.5i){\rm
e}^{2x+1.5ix}}.$ (29)

We calculate the eigenvalues of $q_{GI}(x,0)$ potential with $\hat{a}=0.1$ and
$n_{1}=400$, and the calculated results are shown in Figure 3. The absolute
error between the calculated result and the exact result is $1.09\times
10^{-15}$ at $k=k_{2}$, the absolute error is $4.66\times 10^{-15}$ at
$k=\bar{k}_{2}$, the absolute error is $1.01\times 10^{-14}$ at
$k=-\bar{k}_{2}$, and the absolute error is $2.50\times 10^{-15}$ at
$k=-k_{2}$.

Figure 3: The eigenvalues for the DNLSE $\rm III$ with (29). Blue part: the
exact result; red part: the calculated result.

Although the numerical method performs well for the discrete eigenvalues and
performs poorly for the continuous spectrum, it is still feasible because the
continuous spectrum of the DNLSE is known. For the DNLSE $\rm III$, we explore
the spectral convergence for the discrete eigenvalue $k_{2}$, the calculated
result is shown in Figure 4. We learn that this method is effective and has
machine accuracy.

Figure 4: Demonstration of spectral convergence for the discrete eigenvalue
$k_{2}$.

######  Remark 2.

In this numerical method, $\hat{a}$ should be a relatively small number and is
related to the potential function $q_{0}(x)$. We can initially set $\hat{a}$
to 0.1 and use the convergency to determine whether it is necessary to reduce
$\hat{a}$.

###  3.2 Calculating the reflection coefficient

$S(k)$ is determined by $\mu_{1}$ and $\mu_{2}$

$\mu_{2}(x,t,k)=\mu_{1}(x,t,k){\rm e}^{-{i}\theta\hat{\sigma}_{3}}S(k).$

Because $S(k)$ is independent of $x$ and $t$, we can set $x=0$ and $t=0$. We
mark $\mu_{1}(0,0,k)$ and $\mu_{2}(0,0,k)$ as $\mu_{1}(k)$ and $\mu_{2}(k)$
respectively. The eigenfunctions $\mu_{1}(k)$ and $\mu_{2}(k)$ are calculated
from the $x$-part of the Lax pair (6),

$(\ \mu_{j}(k)\ )_{x}+{i}k^{2}[\sigma_{3},\mu_{j}(k)]=(kQ+Q_{1})\mu_{j}(k),\
j=1,2,$ (30)

where $Q_{1}=\begin{cases}0\sigma_{3}\quad&{\rm for\ the\ DNLS\ I},\\\
-\frac{i}{4}|q(x,0)|^{2}\sigma_{3}\quad&{\rm for\ the\ DNLS\ II},\\\
\frac{i}{2}|q(x,0)|^{2}\sigma_{3}\quad&{\rm for\ the\ DNLS\ III}.\end{cases}$
The eigenfunction $\mu_{1}(k)\rightarrow I$ when $x\rightarrow-\infty$, and
the eigenfunction $\mu_{2}(k)\rightarrow I$ when $x\rightarrow+\infty$.

Set $\Phi_{j}(k)$=$\mu_{j}(k)-I$, ($j$=1, 2), then $\Phi_{1}(k)\rightarrow 0$
when $x\rightarrow-\infty$, $\Phi_{2}(k)\rightarrow 0$ when
$x\rightarrow+\infty$. $\Phi_{j}(k)$ satisfies

$(\ \Phi_{j}(k)\
)_{x}+{i}k^{2}[\sigma_{3},\Phi_{j}(k)]-(kQ+Q_{1})\Phi_{j}(k)=(kQ+Q_{1}),\
j=1,2.$ (31)

Set

$\Phi_{j}=\left(\begin{array}[]{cc}a&b\\\ c&d\\\ \end{array}\right),$

then

$\left(\begin{array}[]{cc}a_{x}&b_{x}\\\ c_{x}&d_{x}\\\
\end{array}\right)+ik^{2}\left(\begin{array}[]{cc}0&2b\\\ -2c&0\\\
\end{array}\right)-\left(\begin{array}[]{cc}q_{2}&kq\\\ -k\bar{q}&-q_{2}\\\
\end{array}\right)\left(\begin{array}[]{cc}a&b\\\ c&d\\\
\end{array}\right)=\left(\begin{array}[]{cc}q_{2}&kq\\\ -k\bar{q}&-q_{2}\\\
\end{array}\right),$

where $q_{2}=\frac{i}{2}|q(x,0)|^{2}$.

The following two equations can be obtained

$\left(\begin{array}[]{cc}a_{x}\\\ c_{x}\\\
\end{array}\right)-\left(\begin{array}[]{cc}q_{2}&kq\\\
-k\bar{q}&-q_{2}+2{i}k^{2}\\\ \end{array}\right)\left(\begin{array}[]{cc}a\\\
c\\\ \end{array}\right)=\left(\begin{array}[]{cc}q_{2}\\\ -k\bar{q}\\\
\end{array}\right),$ (32)

$\left(\begin{array}[]{cc}b_{x}\\\ d_{x}\\\
\end{array}\right)-\left(\begin{array}[]{cc}q_{2}-2{i}k^{2}&kq\\\
-k\bar{q}&-q_{2}\\\ \end{array}\right)\left(\begin{array}[]{cc}b\\\ d\\\
\end{array}\right)=\left(\begin{array}[]{cc}kq\\\ -q_{2}\\\
\end{array}\right).$ (33)

When $j$=1, $\left(\begin{array}[]{cc}a&b\\\ c&d\\\
\end{array}\right)\rightarrow\left(\begin{array}[]{cc}0&0\\\ 0&0\\\
\end{array}\right)$ as $x\rightarrow-\infty$, when $j$=2,
$\left(\begin{array}[]{cc}a&b\\\ c&d\\\
\end{array}\right)\rightarrow\left(\begin{array}[]{cc}0&0\\\ 0&0\\\
\end{array}\right)$ as $x\rightarrow+\infty$.

The Chebyshev collocation method can be used to solve Eq. (32) and Eq. (33).
Using the knowledge in section 3.1, Eq. (32) becomes

$\begin{split}\left[\begin{array}[]{cc}{\rm D}[\frac{{\rm
d}H(\vec{\chi})}{{\rm d}x}]\mathcal{F}^{-1}\mathcal{D}\mathcal{F}-{\rm
D}[q_{2}(\hat{\chi})]&-k{\rm D}[q(\hat{\chi})]\\\ k{\rm
D}[\bar{q}(\hat{\chi})]&{\rm D}[\frac{{\rm d}H(\vec{\chi})}{{\rm
d}x}]\mathcal{F}^{-1}\mathcal{D}\mathcal{F}+{\rm
D}[q_{2}(\hat{\chi})]-2ik^{2}I\\\
\end{array}\right]\left[\begin{array}[]{c}a(\hat{\chi})\\\ c(\hat{\chi})\\\
\end{array}\right]=\left[\begin{array}[]{c}q_{2}(\hat{\chi})\\\
-k\bar{q}(\hat{\chi})\\\ \end{array}\right].\end{split}$ (34)

Eq. (33) becomes

$\begin{split}\left[\begin{array}[]{cc}2ik^{2}I+{\rm D}[\frac{{\rm
d}H(\vec{\chi})}{{\rm d}x}]\mathcal{F}^{-1}\mathcal{D}\mathcal{F}-{\rm
D}[q_{2}(\hat{\chi})]&-k{\rm D}[q(\hat{\chi})]\\\ k{\rm
D}[\bar{q}(\hat{\chi})]&{\rm D}[\frac{{\rm d}H(\vec{\chi})}{{\rm
d}x}]\mathcal{F}^{-1}\mathcal{D}\mathcal{F}+{\rm D}[q_{2}(\hat{\chi})]\\\
\end{array}\right]\left[\begin{array}[]{c}b(\hat{\chi})\\\ d(\hat{\chi})\\\
\end{array}\right]=\left[\begin{array}[]{c}kq(\hat{\chi})\\\
-q_{2}(\hat{\chi})\\\ \end{array}\right].\end{split}$ (35)

######  Remark 3.

The mapping $H(x)$ used in this subsection is different from subsection 3.1.
When $j=1$, $H(x)=2{\rm tanh}(\hat{a}x)+1$, when $j=2$, $H(x)=2{\rm
tanh}(\hat{a}x)-1$.

For the DNLSE $\rm III$, the numerical method ($n_{1}=201$, $\hat{a}=0.05$) is
used to calculate the reflection coefficient $r(k)$ of $q(x,0)=1.5{\rm
exp}(-x^{2})$ potential. The calculated results are shown in Figure 5. It is
worth noting that the continuous spectrum of the DNLSE is $\mathbb{R}\cup
i\mathbb{R}$ rather than $\mathbb{R}$. In Figure 5, the blue line is the
density of $r(k)$, the red line is the real part of $r(k)$, the green line is
the imaginary part of $r(k)$. We also explore the convergency of $r(k)$ at
$k=0.5$, $k=-1$, $k=0.5i$ and $k=-i$. The calculated results are shown in
Figure 6. The ‘Error’ of Figure 6 is defined by

${\rm Error}=\Big{|}\ r(n_{1},k)-r(500,k)\ \Big{|},$

where $r(n_{1},k)$ is the reflection coefficient $r(k)$ calculated by $n_{1}$
Chebyshev nodes. From Figure 6, we learn that the calculated results converge
quickly, so the numerical method is effective.

(a) The density of $r(k)$ in $k\in\mathbb{R}$ (b) The density of $r(k)$ in
$k\in i\mathbb{R}$

Figure 5: The calculated reflection coefficient $r(k)$ for $q(x,0)=1.5{\rm
exp}(-x^{2})$. Blue line: the density of $r(k)$; red line: the real part of
$r(k)$; green line: the imaginary part of $r(k)$

(a) $k=0.5$ (b) $k=-1$ (c) $k=0.5i$ (d) $k=-i$

Figure 6: Demonstration of spectral convergence for the reflection coefficient
$r(k)$.

######  Remark 4.

Considering that we calculate the reflection coefficient in Schwartz space
$\mathcal{S}(\mathbb{R})$, it is also a feasible framework to truncate the
calculation interval into a finite interval[39, 42].

###  3.3 Calculating the normalized constant

When $k=\kappa_{j}$, we need to calculate the normalized constant $C_{j}$. If
$k=\kappa_{j}$ is the one-order pole of $a(k)$, the normalized constant is
defined by

$C_{j}=\frac{b_{j}}{a^{\prime}(\kappa_{j})},$

where $\kappa_{j}$ is the discrete eigenvalue located in the first quadrant of
the complex field $\mathbb{C}$. Set $k=\kappa_{j}$, the $b_{j}$ can be
obtained by Eq. (14). Considering that $a(k)$ is analytic in $D_{2}$, we give
three frames for calculating $a^{\prime}(\kappa_{j})$:

  * 1.

(Central Difference Formula):

$a^{\prime}(\kappa_{j})=\lim_{|\Delta|\rightarrow
0}\frac{a(\kappa_{j}+\vec{\Delta})-a(\kappa_{j}-\vec{\Delta})}{\vec{\Delta}}.$

  * 2.

(Fourier Transform):

$a^{\prime}(z)=\mathbb{F}^{-1}[{i}k\mathbb{F}[a(z)]],$

where $\mathbb{F}$ represents the Fourier transform and $\mathbb{F}^{-1}$
represents the inverse Fourier transform. In the practical calculation, we use
the fast Fourier transform (FFT) to calculate $a^{\prime}(\kappa_{j})$.

  * 3.

(Chebyshev Derivative Matrix)

${\rm If}\ a(z)=T^{\mathbb{R}}(z)\mathcal{F}a(\hat{\chi}),\ {\rm then}\
a^{\prime}(z)=\frac{\partial H(z)}{\partial
z}T^{\mathbb{R}}(z)\mathcal{D}\mathcal{F}a(\hat{\chi}).$

In this paper, the third method is used to calculate $a^{\prime}(\kappa_{j})$.

##  4 Numerical inverse scattering

Numerical inverse scattering consists of two steps, one is deforming the RHP
in the spirit of Deift-Zhou nonlinear steepest descent[26], another one is
solving the deformed Riemann-Hilbert problem numerically by the Chebyshev
collocation method. Considering the original Riemann-Hilbert problem is
oscillatory, the deformation is necessary when time is large.

This section is organized as follows. In section 4.1, we deform the RHP under
not large time and the non-soliton case. In section 4.2, we deform the RHP
under large time and the non-soliton case. In section 4.3, we deform the RHP
under not large time and the soliton case. In section 4.4, some extra
deformations are given. In section 4.5, the numerical method for solving the
RHP is introduced.

######  Remark 5.

The green line is used to represent the jump contours of the Riemann-Hilbert
problem in the complex-$k$ plane, and the red line is used to represent
$f(k)=0$ in the complex-$k$ plane, which applies to all figures presented in
this section. The formula for $f(k)$ will be given in section 4.1.

###  4.1 Deformations for the RHP (the non-soliton case)

Assuming that the solution of the DNLS equation will not evolve into the
soliton. In this condition, the DNLS equation does not exist discrete
eigenvalues, and we do not need to calculate the normalized constant.

In $\mathbb{R}\cup i\mathbb{R}$, ${\rm e}^{2i\theta(k)}$ and ${\rm
e}^{2i\theta(k)}$ are oscillatory, especially when $t$ is large. The jump
matrix $T(k)$ has the oscillation term ${\rm e}^{2i\theta(k)}$ and ${\rm
e}^{-2i\theta(k)}$ at the same time and the jump contours of Riemann-Hilbert
problem is $\mathbb{R}\cup{i}\mathbb{R}$. Due to the existence of oscillation
terms, it is difficulty to numerically solve Riemann-Hilbert problem.
Therefore, we need to deform Riemann-Hilbert problem to keep the jump contours
away from $\mathbb{R}\cup{i}\mathbb{R}$. The method for deforming the
oscillatory Riemann-Hilbert problem is derived from Deift-Zhou nonlinear
steepest descent method[26].

######  Remark 6.

When time is small, we can use the NIST to calculate the solution of the DNLS
equation without any deformation. When time is big, the deformations of the
Riemann-Hilbert problem are necessary.

We need to find the saddle points of
$\varphi(k)={i}\theta={i}(k^{2}x+2k^{4}t).$ Define $\lambda=k^{2}$, then

$\varphi(\lambda)={i}(\lambda x+2\lambda^{2}t).$

Set $\varphi^{\prime}(\lambda)=0$, and we get

$\lambda_{0}=-\frac{x}{4t}$ (36)

which satisfies $\varphi^{\prime}(\lambda_{0})=0$ and
$\varphi^{\prime\prime}(\lambda_{0})\neq 0$ . We get two saddle points

$k_{1}=\sqrt{-\frac{x}{4t}}\ {\rm and}\ k_{2}=-\sqrt{-\frac{x}{4t}}.$

It is obvious that

$k_{3}=0$

is also the saddle point of $\varphi(k)={i}\theta={i}(k^{2}x+2k^{4}t).$

Setting $\lambda={\rm Re}\lambda+i{\rm Im}\lambda$, we get

${\rm Re}({i}\theta)=-4{\rm Im}\lambda\cdot({\rm Re}\lambda-\lambda_{0}).$
(37)

Substitute $\lambda=k^{2}$ into Eq. (37), we get

${\rm Re}({i}\theta)=-4{\rm Im}\ k^{2}\cdot({\rm Re}{\ k^{2}}-\lambda_{0}).$
(38)

Eq. (38) shows that the sign of $\rm Re({i}\theta)$ is determined by ${\rm
Im}\ {k^{2}}$ and $({\rm Re}\ {k^{2}}-\lambda_{0})$.

Given that

${\rm Im}\ (k^{2})={\rm Im}\big{[}({\rm Re}\ k+i{\rm Im}\ k)^{2}\big{]}=2{\rm
Re}\ k\cdot{\rm Im}\ k,$

the sign of ${\rm Im}\ k^{2}$ is shown in Figure 7.

Figure 7: The sign of ${\rm Im}\ k^{2}$ on complex-$k$ plane. Gray area: $\rm
Im\ {\it k^{2}}>0$. Orange area: $\rm Im\ {\it k^{2}}<0$.

After determining the sign of ${\rm Im}\ k^{2}$, we need to determine the sign
of ${\rm Re}\ k^{2}-\lambda_{0}$. We define $f(k)$ by

$f(k)={\rm Re}\ k^{2}-\lambda_{0}={\rm Re}\Big{[}({\rm Re}\ k+i{\rm Im}\
k)^{2}\Big{]}-\lambda_{0}=({\rm Re}\ k)^{2}-({\rm Im}\ k)^{2}-\lambda_{0}.$
(39)

The shape of $f(k)=0$ and the sign of ${\rm Re}\ k^{2}-\lambda_{0}$ are
related to $x$.

  * 1.

If $x<0$, $f(k)=0$ is a hyperbola whose focus is located in $\mathbb{R}$,
whose asymptotes are Re $k$= $\pm$ Im $k$.

  * 2.

If $x=0$, $f(k)=0$ are two straight lines Re $k$= $\pm$ Im $k$.

  * 3.

If $x>0$, $f(k)=0$ is a hyperbola whose focus is located in $i\mathbb{R}$,
whose asymptotes are Re $k$= $\pm$ Im $k$.

The sign of ${\rm Re}\ {k^{2}}-\lambda_{0}$ is shown in Figure 8 where the red
line is $f(k)=0$.

(a) $x<0$ (b) $x=0$ (c) $x>0$

Figure 8: The sign of ${\rm Re}\ {k^{2}}-\lambda_{0}$ in the complex-$k$
plane. Gray area: ${\rm Re}\ {k^{2}}>\lambda_{0}$. Orange area: ${\rm
Re}{k^{2}}<\lambda_{0}$.

The sign of ${\rm Re}(i\theta)$ in the complex-$k$ plane is shown in Figure 9.
The gray area of the Figure 9 is where ${\rm Re}(i\theta)>0$, and we record it
as $\mathrm{F}_{1}$. The orange area of the Figure 9 is where ${\rm
Re}(i\theta)<0$, and we record it as $\mathrm{F}_{2}$.

(a) $x<0$ (b) $x=0$ (c) $x>0$

Figure 9: The sign of ${\rm Re}({i}\theta)$ in the complex-$k$ plane. Gray
area $\mathrm{F}_{1}$: ${\rm Re}(i\theta)>0$. Orange area $\mathrm{F}_{2}$:
${\rm Re}(i\theta)<0$

After determining the sign of ${\rm Re}(i\theta)$, we start to deform Riemann-
Hilbert problem. Two decompositions of the jump matrix $T(x,t,k)$ are
considered, where

$T=\left(\begin{array}[]{cc}1&r(k){\rm e}^{-2i\theta(k)}\\\ -\tilde{r}(k){\rm
e}^{2i\theta(k)}&1-\tilde{r}(k)r(k)\\\ \end{array}\right).$

The first decomposition is $T=LU$, where

$L=\left(\begin{array}[]{cc}1&0\\\ -\tilde{r}(k){\rm e}^{2i\theta(k)}&1\\\
\end{array}\right),\ U=\left(\begin{array}[]{cc}1&r(k){\rm
e}^{-2i\theta(k)}\\\ 0&1\\\ \end{array}\right).$

The second decomposition is $T=ABC$, where

$\ A=\left(\begin{array}[]{cc}1&\frac{r(k){\rm
e}^{-2i\theta(k)}}{1-\tilde{r}(k)r(k)}\\\ 0&1\\\ \end{array}\right),\
B=\left(\begin{array}[]{cc}\frac{1}{1-\tilde{r}(k)r(k)}&0\\\
0&1-\tilde{r}(k)r(k)\\\ \end{array}\right),C=\left(\begin{array}[]{cc}1&0\\\
\frac{-\tilde{r}(k){\rm e}^{2i\theta(k)}}{1-\tilde{r}(k)r(k)}&1\\\
\end{array}\right).$

In $\mathrm{F}_{2}$ , ${\rm e}^{2i\theta(k)}$ decays exponentially; in
$\mathrm{F}_{1}$, ${\rm e}^{-2i\theta(k)}$ is decays exponentially. Matrices
$A$ and $U$ can be extended to the orange region $\mathrm{F}_{1}$, matrices
$C$ and $L$ can be extended to the grey region $\mathrm{F}_{2}$.

There are three saddle points, and the position of the saddle points are
determined by $x$ and $t$. It is very complex to deform the original Riemann-
Hilbert problem, and we will give specific expressions in different
situations. For the DNLS equation, we divide the $(x,t)$-plane into three
regions($c_{1}>0$):

  * 1.

Region 1: $-c_{1}\sqrt{t}\leq x\leq c_{1}\sqrt{t}$.

  * 2.

Region 2: $x<-c_{1}\sqrt{t}$.

  * 3.

Region 3: $x>c_{1}\sqrt{t}$.

In different regions for the $(x,t)$-plane, the deformations are different.
After the proper deformations of the jump contours, the oscillatory Riemann-
Hilbert problem is transformed into a standard Riemann-Hilbert problem or a
weak-oscillatory Riemann-Hilbert problem, the jump matrices are exponentially
decaying away from the saddle points.

####  4.1.1 Deformation in Region 1(-$c_{1}\sqrt{t}\leq x\leq c_{1}\sqrt{t}$)

In this condition, we direct deform the jump contours at $k=0$, not at the
saddle points $k_{1}$ and $k_{2}$.

We define $M^{(1)}(x,t,k)$ by

$M^{(1)}(x,t,k)=M(x,t,k)\begin{cases}I,&{\rm if}\
k\in\mathrm{E}_{3}\cup\mathrm{E}_{6}\cup\mathrm{E}_{9}\cup\mathrm{E}_{12},\\\
C^{-1},&{\rm if}\ k\in\mathrm{E}_{5}\cup\mathrm{E}_{11},\\\ A,&{\rm if}\
k\in\mathrm{E}_{4}\cup\mathrm{E}_{10},\\\ L,&{\rm if}\
k\in\mathrm{E}_{2}\cup\mathrm{E}_{8},\\\ U^{-1},&{\rm if}\
k\in\mathrm{E}_{1}\cup\mathrm{E}_{7},\end{cases}$ (40)

where $\mathrm{E}_{j}\ (j=1,2,...,12)$ is shown in Figure 10.

Then, we get Riemann-Hilbert problem I:

  * 1.

$M^{(1)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus\Sigma_{1}$,
$\Sigma_{1}$=$\Sigma^{(1)}_{1}\cup\cdots\cup\Sigma^{(8)}_{1}\cup\mathbb{R}$,

  * 2.

$M^{(1)}_{+}(x,t,k)=M^{(1)}_{-}(x,t,k)T^{(1)}(x,t,k)$, if $k\in\Sigma_{1}$,

  * 3.

$M^{(1)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where $\Sigma_{1}$ is shown in Figure 10,

$T^{(1)}(x,t,k)=\begin{cases}L,&{\rm if}\
k\in\Sigma^{(3)}_{1}\cup\Sigma^{(7)}_{1},\\\ A,&{\rm if}\
k\in\Sigma^{(4)}_{1}\cup\Sigma^{(8)}_{1},\\\ C,&{\rm if}\
k\in\Sigma^{(1)}_{1}\cup\Sigma^{(5)}_{1},\\\ U,&{\rm if}\
k\in\Sigma^{(2)}_{1}\cup\Sigma^{(6)}_{1},\\\ B,&{\rm if}\
k\in\mathbb{R}.\end{cases}$ (41)

Figure 10: The jump contours for Riemann-Hilbert problem I.

To ensure that this decomposition is reasonable, we give the following
proposition:

######  Proposition 4.1 (Uniform boundedness).

If $-c_{1}\sqrt{t}\leq x\leq c_{1}\sqrt{t}$, there is a constant $C_{0}$ s.t.

  * 1.

$|{\rm e}^{2i\theta(k)}|\leq C_{0}$, for any
$k\in\Sigma^{(1)}_{1}\cup\Sigma^{(3)}_{1}\cup\Sigma^{(5)}_{1}\cup\Sigma^{(7)}_{1}$,

  * 2.

$|{\rm e}^{-2i\theta(k)}|\leq C_{0}$, for any
$k\in\Sigma^{(2)}_{1}\cup\Sigma^{(4)}_{1}\cup\Sigma^{(6)}_{1}\cup\Sigma^{(8)}_{1}$.

###### Proof.

Assume that the argument of $\Sigma^{(1)}_{1}$ is ${\rm arctan}(\frac{1}{10})$
and $c_{1}=4$.

First, we consider $x=0$.

In this condition, all of $\Sigma^{(j)}_{1}(j=2,4,6,8)$ are located in the
orange area $\mathrm{F}_{2}$. In the orange area $\mathrm{F}_{2}$, ${{\rm
Re}(i\theta)}\leq 0$, so $|{\rm e}^{2i\theta(k)}|\leq 1$. Similarly, all of
$\Sigma^{(j)}_{1}(j=1,3,5,7)$ are located in the grey area $\mathrm{F}_{1}$.
In the grey area $\mathrm{F}_{1}$, ${{\rm Re}(-i\theta)}\leq 0$, then we get
$|{\rm e}^{-2i\theta(k)}|\leq 1$.

Then, we consider $-4\sqrt{t}\leq x<0$.

In this condition, $|{\rm e}^{2i\theta(k)}|\leq 1$ for $\forall
k\in\Sigma^{(j)}_{1}(j=3,7)$ and $|{\rm e}^{-2i\theta(k)}|\leq 1$ for $\forall
k\in\Sigma^{(j)}_{1}(j=2,6)$. We need to prove that $|{\rm
e}^{2i\theta(k)}|\leq C_{0}$ for $\forall k\in\Sigma^{(j)}_{1}(j=1,5)$ and
$|{\rm e}^{-2i\theta(k)}|\leq C_{0}$ for $\forall
k\in\Sigma^{(j)}_{1}(j=4,8)$, where $C_{0}$ is a constant. Considering
$\Sigma^{(1)}_{1}$ and $\Sigma^{(5)}_{1}$ are symmetric about $k=0$,
$\Sigma^{(8)}_{1}$ and $\Sigma^{(4)}_{1}$ are symmetric about $k=0$, we only
need to prove that $|{\rm e}^{2i\theta(k)}|\leq C_{0}$ for $\forall
k\in\Sigma^{(1)}_{1}$ and $|{\rm e}^{-2i\theta(k)}|\leq C_{0}$ for $\forall
k\in\Sigma^{(8)}_{1}$.

Both $\Sigma^{(1)}_{1}$ and $\Sigma^{(7)}_{1}$ consist of two parts, one is
located in the gray area $\mathrm{F}_{1}$ and another one is located in the
orange area $\mathrm{F}_{2}$. So $\Sigma^{(1)}_{1}$ and $\Sigma^{(8)}_{1}$ can
be written in the following form

$\Sigma^{(1)}_{1}=\Big{(}\Sigma^{(1)}_{1}\cap\mathrm{F}_{2}\Big{)}\cup\Big{(}\Sigma^{(1)}_{1}\cap\mathrm{F}_{1}\Big{)},\
\Sigma^{(8)}_{1}=\Big{(}\Sigma^{(8)}_{1}\cap\mathrm{F}_{2}\Big{)}\cup\Big{(}\Sigma^{(8)}_{1}\cap\mathrm{F}_{1}\Big{)}.$

Then

$\begin{cases}\Sigma^{(1)}_{1}=\\{\ k\ |\ k=10b+bi\\},&b\in(0,+\infty),\\\
\Sigma^{(1)}_{1}\cap\mathrm{F}_{1}=\\{\ k\ |\ k=10b+b{i}\\},&b\in(0,b_{0}),\\\
\Sigma^{(1)}_{1}\cap\mathrm{F}_{2}=\\{\ k\ |\
k=10b+b{i}\\},&b\in(b_{0},+\infty),\end{cases}$ (42)

where $b_{0}={\rm Im}(k_{0})$ and $k_{0}=a_{0}+b_{0}{i}$ is the intersection
of $\Sigma^{(1)}_{1}$ and $f(k)$. By calculating

$\begin{cases}a_{0}=10b_{0},\\\ a_{0}^{2}-b_{0}^{2}=\lambda_{0},\\\
b_{0}>0,\end{cases}$ (43)

we get $b_{0}=\frac{1}{\sqrt{99}}\sqrt{-\frac{x}{4t}}$ and
$k_{0}=\frac{10}{\sqrt{99}}\sqrt{-\frac{x}{4t}}+\frac{1}{\sqrt{99}}\sqrt{-\frac{x}{4t}}{i}.$
So $k\in\Sigma^{(1)}_{1}\cap\mathrm{F}_{2}$ is equivalent to
$b\in(0,\frac{1}{\sqrt{99}}\sqrt{-\frac{x}{4t}})$.

In the orange area $\mathrm{F}_{2}$, the density of ${\rm e}^{{i}\theta(k)}$
is smaller than one due to ${\rm Re}(\ i\theta(k)\ )<0$. So there is $|{\rm
e}^{2i\theta(k)}|\leq 1$ for $\forall
k\in\mathrm{\Sigma^{(1)}_{1}}\cap\mathrm{F}_{2}$. We only need to prove that
$|{\rm e}^{2i\theta(k)}|\leq C$ for $\forall
k\in\Sigma^{(1)}_{1}\cap\mathrm{F}_{1}$. Substitute $k=10b+bi$ into
$i\theta(k)$, we get

${\rm Re}{(i\theta)}=-7920tb^{4}-20xb^{2}.$ (44)

To prove that ${\rm Re}{(i\theta)}$ is bounded in
$\Sigma^{(1)}_{1}\cap\mathrm{F}_{1}$, we only need to prove that
$f_{1}(b)$=${\rm Re}{(i\theta)}=-7920tb^{4}-20xb^{2}$ is bounded for $\forall
b\in[0,b_{0}]$. For $b\in[0,b_{0}]$,
$max\\{f_{1}(b)\\}=f_{1}(\frac{1}{\sqrt{198}}\sqrt{-\frac{x}{4t}})=\frac{5x^{2}}{396t}.$
If $-4\sqrt{t}\leq x\leq 4\sqrt{t}$, then $f_{1}(b)\leq\frac{20}{99}$, there
is $|{\rm e}^{2i\theta(k)}|\leq{\rm e}^{\frac{40}{99}}$ for $\forall
k\in\Sigma^{(1)}_{1}\cap\mathrm{F}_{2}$.

Similarly, we can prove that $|{\rm e}^{\rm-2i\it\theta(k)}|\leq 1$ for
$\forall k\in\Sigma^{(8)}_{1}\cap\mathrm{F}_{1}$ and $|{\rm
e}^{-2i\theta(k)}|\leq{\rm e}^{\frac{40}{99}}$ for $\forall
k\in\Sigma^{(8)}_{1}\cap\mathrm{F}_{2}$.

Last, we consider $0<x\leq 4\sqrt{t}$.

In this condition, $|{\rm e}^{2i\theta(k)}|\leq 1$ for $\forall
k\in\Sigma^{(j)}_{1}(j=1,5)$ and $|{\rm e}^{-2i\theta(k)}|\leq 1$ for $\forall
k\in\Sigma^{(j)}_{1}(j=4,8)$. We need to prove that $|{\rm
e}^{2i\theta(k)}|\leq C_{0}$ for $\forall k\in\Sigma^{(j)}_{1}(j=3,7)$ and
$|{\rm e}^{-2i\theta(k)}|\leq C_{0}$ for $\forall
k\in\Sigma^{(j)}_{1}(j=2,6)$. Considering $\Sigma^{(3)}_{1}$ and
$\Sigma^{(7)}_{1}$ are symmetric about $k=0$, $\Sigma^{(2)}_{1}$ and
$\Sigma^{(6)}_{1}$ are symmetric about $k=0$, we only need to prove that
$|{\rm e}^{2i\theta(k)}|\leq C_{0}$ for $\forall k\in\Sigma^{(3)}_{1}$ and
$|{\rm e}^{-2i\theta(k)}|\leq C_{0}$ for $\forall k\in\Sigma^{(2)}_{1}$.

Both $\Sigma^{(2)}_{1}$ and $\Sigma^{(3)}_{1}$ consist of two parts, one in
the gray area $\mathrm{F}_{1}$ and one in the orange area $\mathrm{F}_{2}$. So
$\Sigma^{(2)}_{1}$ and $\Sigma^{(3)}_{1}$ can be written in the following form

$\Sigma^{(2)}_{1}=\Big{(}\Sigma^{(2)}_{1}\cap\mathrm{F}_{1}\Big{)}\cup\Big{(}\Sigma^{(2)}_{1}\cap\mathrm{F}_{2}\Big{)},\
\Sigma^{(3)}_{1}=\Big{(}\Sigma^{(3)}_{1}\cap\mathrm{F}_{1}\Big{)}\cup\Big{(}\Sigma^{(3)}_{1}\cap\mathrm{F}_{2}\Big{)}.$

Assume that the argument of $\Sigma^{(3)}_{1}$ is $\frac{\pi}{2}+{\rm
arctan}(\frac{1}{10})$. Then

$\begin{cases}\mathrm{\Sigma^{(3)}_{1}}=\\{\ k\ |\
k=-0.1b+bi\\},&b\in(0,+\infty),\\\
\mathrm{\Sigma^{(3)}_{1}}\cap\mathrm{F}_{1}=\\{\ k\ |\
k=-0.1b+bi\\},&b\in(0,b_{0}),\\\
\mathrm{\Sigma^{(3)}_{1}}\cap\mathrm{F}_{2}=\\{\ k\ |\
k=-0.1b+bi\\},&b\in(b_{0},+\infty),\end{cases}$ (45)

where $k_{0}=a_{0}+b_{0}i$ is the intersection of $\Sigma^{(3)}_{1}$ and
$f(k)$, and $b_{0}={\rm Im}(k_{0})$. By calculating

$\begin{cases}a_{0}=-0.1b_{0},\\\ a_{0}^{2}-b_{0}^{2}=\lambda_{0},\\\
b_{0}>0,\end{cases}$ (46)

we get $b_{0}=\frac{10}{\sqrt{99}}\sqrt{\frac{x}{4t}}$ and
$k_{0}=-\frac{1}{\sqrt{99}}\sqrt{\frac{x}{4t}}+\frac{10}{\sqrt{99}}\sqrt{\frac{x}{4t}}i.$

In the orange area $\mathrm{F}_{2}$, the density of ${\rm exp}(\ i\theta(k)\
)$ is smaller than one due to ${\rm Re}(\ i\theta(k)\ )<0$. So there is $|{\rm
e}^{\rm 2i\it\theta(k)}|\leq 1$ for $\forall
k\in\Sigma^{(3)}_{1}\cap\mathrm{F}_{2}$. We only need to prove that $|{\rm
e}^{2i\theta(k)}|\leq C_{0}$ for $\forall
k\in\Sigma^{(3)}_{1}\cap\mathrm{F}_{1}$. Substitute $k=-0.1b+bi$ into
$i\theta(k)$, and we get

${\rm Re}{(i\theta)}=-0.792tb^{4}+0.2xb^{2}.$ (47)

To prove that ${\rm Re}{(i\theta)}$ is bounded in
$\Sigma^{(3)}_{1}\cap\mathrm{F}_{1}$, we only need to prove that
$f_{2}(b)$=${\rm Re}{(i\theta)}=-0.792tb^{4}+0.2xb^{2}$ is bounded for
$\forall b\in(0,b_{0})$. For $b\in(0,b_{0})$,
$max\\{f_{2}(b)\\}=f_{2}(\frac{10}{\sqrt{198}}\sqrt{-\frac{x}{4t}})=\frac{5x^{2}}{396t}.$
If $-4\sqrt{t}\leq x\leq 4\sqrt{t}$, then $f_{2}(b)\leq\frac{20}{99}$ and
$|{\rm e}^{2i\theta(k)}|\leq{\rm e}^{\frac{40}{99}}$ for $\forall
k\in\Sigma^{(3)}_{1}\cap\mathrm{F}_{2}$. Similarly, we can prove that $|{\rm
e}^{-2i\theta(k)}|\leq 1$ for $\forall k\in\Sigma^{(2)}_{1}\cap\mathrm{F}_{1}$
and $|{\rm e}^{-2i\theta(k)}|\leq{\rm e}^{\frac{40}{99}}$ for $\forall
k\in\Sigma^{(2)}_{1}\cap\mathrm{F}_{2}$.

Our proof is complete, the constant $C_{0}$ is ${\rm e}^{\frac{40}{99}}$.

∎

Our decomposition transforms the original Riemann-Hilbert problem into a week-
oscillatory Riemann-Hilbert problem. Proposition 4.1 gives the uniform
boundedness of the jump matrices. This decomposition can both ensure numerical
accuracy and minimize computational costs.

####  4.1.2 Deformation in Region 2 ($x<-c_{1}\sqrt{t}$)

If we directly deform the jump contours at $k=0$, the jump matrices are not
uniformly bounded, which will lead to numerical accuracy issues. When
$x<-c_{1}\sqrt{t}$, the saddle points $k=k_{1}$ and $k=k_{2}$ are located in
$\mathbb{R}$, we need to decompose the jump contours at the saddle points
$k=0$, $k=k_{1}$ and $k=k_{2}$.

We define $M^{(2)}(x,t,k)$ by

$M^{(2)}(x,t,k)=M(x,t,k)\begin{cases}I,&{\rm if}\
k\in\mathrm{E}^{2}_{3}\cup\mathrm{E}^{2}_{6}\cup\mathrm{E}^{2}_{9}\cup\mathrm{E}^{2}_{12},\\\
U^{-1},&{\rm if}\
k\in\mathrm{E}^{2}_{1}\cup\mathrm{E}^{2}_{7}\cup\mathrm{E}^{2}_{14}\cup\mathrm{E}^{2}_{16},\\\
L,&{\rm if}\
k\in\mathrm{E}^{2}_{2}\cup\mathrm{E}^{2}_{8}\cup\mathrm{E}^{2}_{13}\cup\mathrm{E}^{2}_{15},\\\
A,&{\rm if}\ k\in\mathrm{E}^{2}_{4}\cup\mathrm{E}^{2}_{10},\\\ C^{-1},&{\rm
if}\ k\in\mathrm{E}^{2}_{5}\cup\mathrm{E}^{2}_{11},\end{cases}$ (48)

where the area $\mathrm{E}^{2}_{j}\ (j=1,2,...,16)$ is shown in Figure 11.

Figure 11: The area $\mathrm{E}^{(2)}_{j}\ (j=1,2,...,16)$ in the complex-$k$
plane.

We get Riemann-Hilbert problem II:

  * 1.

$M^{(2)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus\Sigma_{2}$,
$\mathrm{\Sigma_{2}}$=$\Sigma^{(1)}_{2}\cup\cdots\cup\Sigma^{(16)}_{2}\cup(-\infty,k_{2})\cup(k_{1},+\infty)$,

  * 2.

$M^{(2)}_{+}(x,t,k)=M^{(2)}_{-}(x,t,k)T^{(2)}(x,t,k)$, if
$k\in\mathrm{\Sigma_{2}}$,

  * 3.

$M^{(2)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where the jump contour $\Sigma^{(j)}_{2}(j=1,2,...,16)$ is shown in Figure 12,

$T^{(2)}(x,t,k)=\begin{cases}A,&{\rm if}\
k\in\Sigma^{(4)}_{2}\cup\Sigma^{(8)}_{2},\\\ C,&{\rm if}\
k\in\Sigma^{(1)}_{2}\cup\Sigma^{(5)}_{2},\\\ U,&{\rm if}\
k\in\Sigma^{(2)}_{2}\cup\Sigma^{(6)}_{2}\cup\Sigma^{(13)}_{2}\cup\Sigma^{(14)}_{2}\cup\Sigma^{(15)}_{2}\cup\Sigma^{(16)}_{2},\\\
L,&{\rm if}\
k\in\Sigma^{(3)}_{2}\cup\Sigma^{(7)}_{2}\cup\Sigma^{(9)}_{2}\cup\Sigma^{(10)}_{2}\cup\Sigma^{(11)}_{2}\cup\Sigma^{(12)}_{2},\\\
B,&{\rm if}\ k\in(-\infty,k_{2})\cup(k_{1},+\infty).\end{cases}$ (49)

Figure 12: The jump contours for Riemann-Hilbert problem II in the complex-$k$
plane.

There are 18 jump contours for Riemann-Hilbert problem II. Compare to Riemann-
Hilbert problem I, Riemann-Hilbert problem II has 8 extra jump contours. The
extra jump contours are $\Sigma^{(j)}_{2}(j=9,10,...,16)$.

####  4.1.3 Deformation in Region 3($x>c_{1}\sqrt{t}$)

When $x>c_{1}\sqrt{t}$, the saddle points are located in ${i}\mathbb{R}$, we
need to deform the jump contours at the saddle points $k=0$, $k=k_{1}$ and
$k=k_{2}$.

We define $M^{(3)}(x,t,k)$ by

$M^{(3)}(x,t,k)=M(x,t,k)\begin{cases}I,&{\rm if}\
k\in\mathrm{E}^{3}_{3}\cup\mathrm{E}^{3}_{6}\cup\mathrm{E}^{3}_{9}\cup\mathrm{E}^{3}_{12},\\\
U^{-1},&{\rm if}\
k\in\mathrm{E}^{3}_{1}\cup\mathrm{E}^{3}_{7}\cup\mathrm{E}^{3}_{14}\cup\mathrm{E}^{3}_{16},\\\
L,&{\rm if}\
k\in\mathrm{E}^{3}_{2}\cup\mathrm{E}^{3}_{8}\cup\mathrm{E}^{3}_{13}\cup\mathrm{E}^{3}_{15},\\\
A,&{\rm if}\ k\in\mathrm{E}^{3}_{4}\cup\mathrm{E}^{3}_{10},\\\ C^{-1},&{\rm
if}\ k\in\mathrm{E}^{3}_{5}\cup\mathrm{E}^{3}_{11},\end{cases}$ (50)

where the area $\mathrm{E}^{3}_{j}\ (j=1,2,...,16)$ is shown in Figure 13.

Figure 13: The area $\mathrm{E}^{3}_{j}\ (j=1,2,...,16)$ in the complex-$k$
plane.

We get Riemann-Hilbert problem III:

  * 1.

$M^{(3)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus\Sigma_{3}$,
$\Sigma_{3}=\Sigma^{(1)}_{3}\cup\cdots\cup\Sigma^{(16)}_{3}\cup(0,k_{2})\cup(0,k_{1})\cup\mathbb{R}$,

  * 2.

$M^{(3)}_{+}(x,t,k)=M^{(3)}_{-}(x,t,k)T^{(3)}(x,t,k)$, if $k\in\Sigma_{3}$,

  * 3.

$M^{(3)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where the jump contour $\Sigma^{(j)}_{3}(j=1,2,...,16)$ is shown in Figure 14,

$T^{(3)}(x,t,k)=\begin{cases}A,&{\rm if}\
k\in\Sigma^{(4)}_{3}\cup\Sigma^{(8)}_{3}\cup\Sigma^{(9)}_{3}\cup\Sigma^{(10)}_{3}\cup\Sigma^{(11)}_{3}\cup\Sigma^{(12)}_{3},\\\
C,&{\rm if}\
k\in\Sigma^{(1)}_{3}\cup\Sigma^{(5)}_{3}\cup\Sigma^{(13)}_{3}\cup\Sigma^{(14)}_{3}\cup\Sigma^{(15)}_{3}\cup\Sigma^{(16)}_{3},\\\
L,&{\rm if}\ k\in\Sigma^{(3)}_{3}\cup\Sigma^{(7)}_{3},\\\ U,&{\rm if}\
k\in\Sigma^{(2)}_{3}\cup\Sigma^{(6)}_{3},\\\ B,&{\rm if}\
k\in(0,k_{2})\cup(0,k_{1})\cup\mathbb{R}.\end{cases}$ (51)

Figure 14: The jump contours for Riemann-Hilbert problem III in the
complex-$k$ plane.

###  4.2 Removing the jump matrix $B$ (the non-soliton case)

Both Riemann-Hilbert problem I, Riemann-Hilbert problem II and Riemann-Hilbert
problem III are not oscillatory. When $t$ is not too large, we can solve them
by the Chebyshev collocation method. When $t$ is large, we need to remove the
jump matrix $B$.

####  4.2.1 Removing the jump matrix $B$ in Region 1

For Riemann-Hilbert problem I, we need to remove the jump on $(0,-\infty)$ and
$(0,+\infty)$. We consider a scalar Riemann-Hilbert problem:

  * 1.

$\delta_{j}(k)(j=1,2)$ is meromorphic in $\mathbb{C}\setminus{i}\mathbb{R}$,

  * 2.

$\begin{cases}\delta_{j,+}(k)=\delta_{j,-}(k)(1+r(k)\overline{r(\bar{k})}),\begin{cases}k\in(0,-\infty)\
{\rm for}\ j=1,\\\ k\in(0,+\infty)\ {\rm for}\ j=2,\end{cases}\\\
\delta_{j,+}(k)=\delta_{j,-}(k),\begin{cases}k\in(0,+\infty)\ {\rm for}\
j=1,\\\ k\in(0,-\infty)\ {\rm for}\ j=2,\end{cases}\end{cases}$

  * 3.

$\delta_{j}(k)\rightarrow 1$, $k\rightarrow\infty$.

When $j=1$, the Riemann-Hilbert problem exists the solution

$\delta_{1}(k)={\exp}\Big{[}{\frac{1}{2\pi i}}\int_{0}^{-\infty}\frac{{\rm
ln}(1+r(z)\overline{r(\bar{z})})}{z-k}{\rm d}z\Big{]},$ (52)

when $j=2$, the Riemann-Hilbert problem exists the solution

$\delta_{2}(k)={\exp}\Big{[}{\frac{1}{2\pi i}}\int_{0}^{+\infty}\frac{{\rm
ln}(1+r(z)\overline{r(\bar{z})})}{z-k}{\rm d}z\Big{]}.$ (53)

Set $\Delta_{(j)}=\left(\begin{array}[]{cc}\frac{1}{\delta_{j}(k)}&0\\\
0&\delta_{j}(k)\\\ \end{array}\right)(j=1,2)$, $\Delta_{(1)}$ can be used to
remove the jump on $(0,-\infty)$, $\Delta_{(2)}$ can be used to remove the
jump on $(0,+\infty)$.

For the RHP

$M^{(1)}_{+}(x,t,k)=M^{(1)}_{-}(x,t,k)B,k\in(-\infty,0),$ (54)

we multiply the solution of the RHP (54) by $\Delta_{(1)}^{-1}$ to remove the
jump of $(-\infty,0)$. In $(-\infty,0)$,

$M^{(1)}_{+}\Delta_{(1),+}^{-1}=M^{(1)}_{-}B\Delta_{(1),+}^{-1}=M^{(1)}_{-}\Delta_{(1),-}^{-1}\Delta_{(1),-}B\Delta_{(1),+}^{-1}=M^{(1)}_{-}\Delta_{(1),-}^{-1}.$

Indeed, we see that there is no jump for $M^{(1)}\Delta_{(1)}^{-1}$. We can
prove that $\delta_{1}(k)=\delta_{2}(-k)$ by the symmetries for $r(k)$.

Doing the transformation

$M^{(4)}(x,t,k)=M^{(1)}(x,t,k)\Delta_{1}^{-1},$ (55)

where
$\Delta_{1}^{-1}=\Delta_{(1)}^{-1}\Delta_{(2)}^{-1}=\left(\begin{array}[]{cc}\frac{1}{\delta_{1}(k)\delta_{1}(-k)}&0\\\
0&\delta_{1}(k)\delta_{1}(-k)\\\ \end{array}\right)$, we get Riemann-Hilbert
problem IV:

  * 1.

$M^{(4)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus\Sigma_{4}$,
$\Sigma_{4}=\Sigma^{(1)}_{1}\cup\cdots\cup\Sigma^{(8)}_{1}$,

  * 2.

$M^{(4)}_{+}(x,t,k)=M^{(4)}_{-}(x,t,k)\begin{cases}\Delta_{1}L\Delta_{1}^{-1},&{\rm
if}\ k\in\Sigma^{(3)}_{1}\cup\Sigma^{(7)}_{1},\\\
\Delta_{1}A\Delta_{1}^{-1},&{\rm if}\
k\in\Sigma^{(4)}_{1}\cup\Sigma^{(8)}_{1},\\\ \Delta_{1}C\Delta_{1}^{-1},&{\rm
if}\ k\in\Sigma^{(1)}_{1}\cup\Sigma^{(5)}_{1},\\\
\Delta_{1}U\Delta_{1}^{-1},&{\rm if}\
k\in\Sigma^{(2)}_{1}\cup\Sigma^{(6)}_{1}.\end{cases}$

  * 3.

$M^{(4)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where $\Sigma_{1}^{(j)}(j=1,2,...,8)$ is shown in Figure 10.

######  Remark 7.

For the DNLSE I, $k=0$ is not the singularity of $\delta_{1}(k)$, so the
deformation (55) is sufficient.

###### Proof.

We divide $\delta_{1}(0)$ into two parts,

$\delta_{1}(0)={\rm
exp}\Big{[}{\frac{1}{2\pi{i}}}\int^{\varepsilon}_{-\infty}\frac{{\rm
log}(1+r(z)\overline{r(\bar{z})})}{z}{\rm d}z\Big{]}+{\rm
exp}\Big{[}{\frac{1}{2\pi{i}}}\int^{0}_{\varepsilon}\frac{{\rm
log}(1+r(z)\overline{r(\bar{z})})}{z}{\rm d}z\Big{]},$

where $\varepsilon$ is a small number.

Using the asymptotic of the eigenfunctions, we learn that
$r(z)=O(\frac{1}{z})\ {\rm as}\ z\rightarrow\infty,$ so

${\rm exp}\Big{[}{\frac{1}{2\pi{i}}}\int^{\varepsilon}_{-\infty}\frac{{\rm
log}(1+r(z)\overline{r(\bar{z})})}{z}{\rm d}z\Big{]}={\rm
exp}\Big{[}{\frac{1}{2\pi{i}}}\int^{\varepsilon}_{-\infty}\frac{{\rm
log}(1+|r(z)|^{2})}{z}{\rm d}z\Big{]}$

converges.

Considering that $r(0)=0$, so

$\frac{{\rm d}[{\rm log}(1+|r(z)|^{2})]}{{\rm
d}z}\mid_{z=0}=\lim_{z\rightarrow 0}\frac{{\rm log}(1+|r(z)|^{2})-0}{z-0},$

so $\int^{0}_{\varepsilon}\frac{{\rm log}(1+|r(z)|^{2})}{z}{\rm d}z$
convergences when $\frac{{\rm d}[{\rm log}(1+|r(z)|^{2})]}{{\rm d}z}|_{z=0}$
is bounded.

∎

####  4.2.2 Removing the jump matrix $B$ in Region 2

In this condition, we need to remove the jump on $(k_{2},-\infty)$ and
$(k_{1},+\infty)$. We consider a scalar Riemann-Hilbert problem:

  * 1.

$\delta_{j}(k)(j=3,4)$ is meromorphic in $\mathbb{C}\setminus\mathbb{R}$,

  * 2.

$\begin{cases}\delta_{j,+}(k)=\delta_{j,-}(k)(1+r(k)\overline{r(\bar{k})}),\begin{cases}k\in(k_{2},-\infty)\
{\rm for}\ j=3,\\\ k\in(k_{1},+\infty)\ {\rm for}\ j=4,\end{cases}\\\
\delta_{j,+}(k)=\delta_{j,-}(k),\begin{cases}k\in(k_{2},+\infty)\ {\rm for}\
j=3,\\\ k\in(k_{1},-\infty)\ {\rm for}\ j=4,\end{cases}\end{cases}$

  * 3.

$\delta_{j}(k)\rightarrow 1$, $k\rightarrow\infty$.

The Riemann-Hilbert problem exists the solution

$\delta_{j}(k)=\begin{cases}{\rm
exp}\Big{[}{\frac{1}{2\pi{i}}}\int_{k_{2}}^{-\infty}\frac{{\rm
ln}(1+r(z)\overline{r(\bar{z})})}{z-k}{\rm d}z\Big{]},&{\rm for}j=3,\\\ {\rm
exp}\Big{[}{\frac{1}{2\pi{i}}}\int_{k_{1}}^{+\infty}\frac{{\rm
ln}(1+r(z)\overline{r(\bar{z})})}{z-k}{\rm d}z\Big{]},&{\rm
for}j=4,\end{cases}$ (56)

and $\delta_{3}(k)=\delta_{4}(-k)$. We can use
$\left(\begin{array}[]{cc}\frac{1}{\delta_{3}(k)}&0\\\ 0&\delta_{3}(k)\\\
\end{array}\right)$ to remove the jump on $(k_{2},-\infty)$, and we can use
$\left(\begin{array}[]{cc}\frac{1}{\delta_{4}(k)}&0\\\ 0&\delta_{4}(k)\\\
\end{array}\right)$ to remove the $\Delta_{(2)}$ can be used to remove the
jump on $(k_{1},+\infty)$. It is worth noting that
$\delta(k)=\delta_{3}(k)\delta_{3}(-k)$ has singular points $k=k_{1}$ and
$k=k_{2}$, so the jump contour should be away from singularities.

We define $M^{(5)}(x,t,k)$ by

$M^{(5)}(x,t,k)=M^{(2)}(x,t,k)\begin{cases}I,&\quad{\rm if}\ k\in D_{21},\\\
A,&\quad{\rm if}\ k\in D_{22},\\\ L^{-1}A,&\quad{\rm if}\ k\in D_{23},\\\
C^{-1}B^{-1},&\quad{\rm if}\ k\in D_{24},\\\ B^{-1},&\quad{\rm if}\ k\in
D_{25}.\\\ \Delta_{2}^{-1},&\quad{\rm if}\ k\in{\rm others},\end{cases}$ (57)

where
$\Delta_{2}=\left(\begin{array}[]{cc}\frac{1}{\delta_{3}(k)\delta_{3}(-k)}&0\\\
0&\delta_{3}(k)\delta_{3}(-k)\\\ \end{array}\right)$, the area
$D_{2i}(i=1,2,...,5)$ is shown in Figure 15.

We get Riemann-Hilbert problem V:

  * 1.

$M^{(5)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus\Sigma_{5}$,
$\Sigma_{5}$=$\Sigma_{2}^{(1)}\cup\cdots\cup\Sigma_{2}^{(21)}$,

  * 2.

$M^{(5)}_{+}(x,t,k)=M^{(5)}_{-}(x,t,k)T^{(5)}(x,t,k)$, if $k\in\Sigma_{5}$,

  * 3.

$M^{(5)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where $\Sigma_{2}^{(1)},\cdots,\Sigma_{2}^{(21)}$ are shown in Figure 12 and
Figure 15,

$T^{(5)}(x,t,k)=\begin{cases}\Delta_{2}C\Delta_{2}^{-1},&{\rm if}\
k\in\Sigma^{(1)}_{2}\cup\Sigma^{(5)}_{2},\\\ \Delta_{2}U\Delta_{2}^{-1},&{\rm
if}\
k\in\Sigma^{(2)}_{2}\cup\Sigma^{(6)}_{2}\cup\Sigma^{(13)}_{2}\cup\Sigma^{(14)}_{2}\cup\Sigma^{(15)}_{2}\cup\Sigma^{(16)}_{2},\\\
\Delta_{2}L\Delta_{2}^{-1},&{\rm if}\
k\in\Sigma^{(3)}_{2}\cup\Sigma^{(7)}_{2}\cup\Sigma^{(9)}_{2}\cup\Sigma^{(10)}_{2}\cup\Sigma^{(11)}_{2}\cup\Sigma^{(12)}_{2},\\\
\Delta_{2}A\Delta_{2}^{-1},&{\rm if}\
k\in\Sigma^{(4)}_{2}\cup\Sigma^{(8)}_{2},\\\ \Delta_{2},\quad&{\rm if}\
k\in\Sigma_{2}^{(17)},\\\ \Delta_{2}A,\quad&{\rm if}\
k\in\Sigma_{2}^{(18)},\\\ \Delta_{2}L^{-1}A,\quad&{\rm if}\
k\in\Sigma_{2}^{(19)},\\\ \Delta_{2}C^{-1}B^{-1},\quad&{\rm if}\
k\in\Sigma_{2}^{(20)},\\\ \Delta_{2}B^{-1},\quad&{\rm if}\
k\in\Sigma_{2}^{(21)}.\end{cases}$ (58)

Figure 15: The jump contours for the Riemann-Hilbert problem VI in the
complex-$k$ plane.

####  4.2.3 Removing the jump matrix $B$ in Region 3

In this condition, we need to remove the jump on $(-\infty,0)$, $(+\infty,0)$,
$(0,k_{1})$ and $(0,k_{2})$. First, we remove the jump on $(0,k_{1})$ and
$(0,k_{2})$. Then, we remove the jump on $(-\infty,0)$ and $(+\infty,0)$, the
method for removing the jump on $(-\infty,0)$ and $(+\infty,0)$ is introduced
in section 4.2.1.

In order to remove the jump on $(0,k_{1})$ and $(0,k_{2})$, we need two
functions

$\delta_{5}(k)={\rm exp}\Big{[}{\frac{1}{2\pi{i}}}\int^{k_{2}}_{0}\frac{{\rm
log}(1+r(z)\overline{r(\bar{z})})}{z-k}{\rm d}z\Big{]},\delta_{6}(k)={\rm
exp}\Big{[}{\frac{1}{2\pi{i}}}\int^{k_{1}}_{0}\frac{{\rm
log}(1+r(z)\overline{r(\bar{z})})}{z-k}{\rm d}z\Big{]}.$

First, we define $M^{(6)}(x,t,k)$ by

$M^{(6)}(x,t,k)=M^{(3)}(x,t,k)\begin{cases}I,&\quad{\rm if}\ k\in D_{11},\\\
A,&\quad{\rm if}\ k\in D_{12},\\\ L^{-1}A,&\quad{\rm if}\ k\in D_{13},\\\
C^{-1}B^{-1},&\quad{\rm if}\ k\in D_{14},\\\ B^{-1},&\quad{\rm if}\ k\in
D_{15}.\\\ \Delta_{3}^{-1},&\quad{\rm if}\ k\in{\rm others},\end{cases}$ (59)

where
$\Delta_{3}=\left(\begin{array}[]{cc}\frac{1}{\delta_{5}(k)\delta_{5}(-k)}&0\\\
0&\delta_{5}(k)\delta_{5}(-k)\\\ \end{array}\right)$, and the area
$D_{1i}(i=1,2,...,5)$ is shown in Figure 16.

Then, we do a transformation

$M^{(6)}(x,t,k)=M^{(6)}(x,t,k)\Delta_{1}^{-1}.$

We get Riemann-Hilbert problem VI:

  * 1.

$M^{(6)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus\Sigma_{6}$,
$\Sigma_{6}$=$\Sigma_{3}^{(1)}\cup\cdots\cup\Sigma_{3}^{(21)}$,

  * 2.

$M^{(6)}_{+}(x,t,k)=M^{(6)}_{-}(x,t,k)T^{(6)}(x,t,k)$, if $k\in\Sigma_{6}$,

  * 3.

$M^{(6)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where $\Sigma_{3}^{(1)},\cdots,\Sigma_{3}^{(21)}$ are shown in Figure 14 and
Figure 16,

$T^{(6)}(x,t,k)=\begin{cases}\Delta_{1}\Delta_{3}U\Delta_{3}^{-1}\Delta_{1}^{-1},&{\rm
if}\ k\in\Sigma^{(2)}_{3}\cup\Sigma^{(6)}_{3},\\\
\Delta_{1}\Delta_{3}C\Delta_{3}^{-1}\Delta_{1}^{-1},&{\rm if}\
k\in\Sigma^{(1)}_{3}\cup\Sigma^{(5)}_{3}\cup\Sigma^{(13)}_{3}\cup\Sigma^{(14)}_{3}\cup\Sigma^{(15)}_{3}\cup\Sigma^{(16)}_{3},\\\
\Delta_{1}\Delta_{3}A\Delta_{3}^{-1}\Delta_{1}^{-1},&{\rm if}\
k\in\Sigma^{(4)}_{3}\cup\Sigma^{(8)}_{3}\cup\Sigma^{(9)}_{3}\cup\Sigma^{(10)}_{3}\cup\Sigma^{(11)}_{3}\cup\Sigma^{(12)}_{3},\\\
\Delta_{1}\Delta_{3}L\Delta_{3}^{-1}\Delta_{1}^{-1},&{\rm if}\
k\in\Sigma^{(3)}_{3}\cup\Sigma^{(7)}_{3},\\\
\Delta_{1}\Delta_{3}\Delta_{1}^{-1},\quad&{\rm if}\ k\in\Sigma_{3}^{(17)},\\\
\Delta_{1}\Delta_{3}A\Delta_{1}^{-1},\quad&{\rm if}\ k\in\Sigma_{3}^{(18)},\\\
\Delta_{1}\Delta_{3}L^{-1}A\Delta_{1}^{-1},\quad&{\rm if}\
k\in\Sigma_{3}^{(19)},\\\
\Delta_{1}\Delta_{3}C^{-1}B^{-1}\Delta_{1}^{-1},\quad&{\rm if}\
k\in\Sigma_{3}^{(20)},\\\ \Delta_{1}\Delta_{3}B^{-1}\Delta_{1}^{-1},\quad&{\rm
if}\ k\in\Sigma_{3}^{(21)}.\end{cases}$ (60)

Figure 16: The jump contours for the Riemann-Hilbert problem VI in the
complex-$k$ plane.

###  4.3 Deformations for the RHP (the soliton case)

Lastly, we deal with the addition of solitons for the case of the DNLS
equation. We need to transform the jump conditions into the jump conditions.
We assume that the set $\\{\kappa_{j}\\}(j=1,2,...)$ is bounded away from
${i}\mathbb{R}\cup\mathbb{R}$.

Fix $\varepsilon>0$, so that the circles
$A^{(1)}_{j}=\\{k\in\mathbb{C}:|k-\kappa_{j}|=\varepsilon\\}$,
$A^{(2)}_{j}=\\{k\in\mathbb{C}:|k-\bar{\kappa}_{j}|=\varepsilon\\}$,
$A^{(3)}_{j}=\\{k\in\mathbb{C}:|k+\bar{\kappa}_{j}|=\varepsilon\\}$ and
$A^{(4)}_{j}=\\{k\in\mathbb{C}:|k+\kappa_{j}|=\varepsilon\\}$ do not intersect
each other and $i\mathbb{R}\cup\mathbb{R}$.

We define $N(x,t,k)$ by

$N(x,t,k)=M(x,t,k)\begin{cases}\left(\begin{array}[]{cc}1&\frac{C_{j}{\rm
e}^{-2i\theta{(\kappa_{j})}}}{k-\kappa_{j}}\\\ 0&1\\\ \end{array}\right),&{\rm
if}\ |k-\kappa_{j}|<\varepsilon,\\\ \left(\begin{array}[]{cc}1&0\\\
\frac{\bar{C}_{j}{\rm
e}^{{2i}\theta{(\bar{\kappa}_{j})}}}{k-\bar{\kappa}_{j}}&1\\\
\end{array}\right),&{\rm if}\ |k-\bar{\kappa}_{j}|<\varepsilon,\\\
\left(\begin{array}[]{cc}1&0\\\ \frac{\bar{C}_{j}{\rm
e}^{{2i}\theta{(-\bar{\kappa}_{j})}}}{k+\bar{\kappa}_{j}}&1\\\
\end{array}\right),&{\rm if}\ |k+\bar{\kappa}_{j}|<\varepsilon,\\\
\left(\begin{array}[]{cc}1&\frac{C_{j}{\rm
e}^{-{2i}\theta{(-\kappa_{j})}}}{k+\kappa_{j}}\\\ 0&1\\\
\end{array}\right),&{\rm if}\ |k+\kappa_{j}|<\varepsilon,\\\ I,&{\rm if}\
k\in{\rm others},\end{cases}$ (61)

and the pole conditions are replaced by the jump conditions

$N_{+}(x,t,k)=N_{-}(x,t,k)\begin{cases}\left(\begin{array}[]{cc}1&\frac{C_{j}{\rm
e}^{-{2i}\theta{(\kappa_{j})}}}{k-\kappa_{j}}\\\ 0&1\\\
\end{array}\right),&k\in A^{(1)}_{j},\\\ \left(\begin{array}[]{cc}1&0\\\
\frac{\bar{C}_{j}{\rm
e}^{{2i}\theta{(\bar{\kappa}_{j})}}}{k-\bar{\kappa}_{j}}&1\\\
\end{array}\right),&k\in A^{(2)}_{j},\\\ \left(\begin{array}[]{cc}1&0\\\
\frac{\bar{C}_{j}{\rm
e}^{{2i}\theta{(-\bar{\kappa}_{j})}}}{k+\bar{\kappa}_{j}}&1\\\
\end{array}\right),&k\in A^{(3)}_{j},\\\
\left(\begin{array}[]{cc}1&\frac{C_{j}{\rm
e}^{-{2i}\theta{(-\kappa_{j})}}}{k+\kappa_{j}}\\\ 0&1\\\
\end{array}\right),&k\in A^{(4)}_{j},\end{cases}$ (62)

where $A^{(1)}_{j}$ and $A^{(4)}_{j}$ have counter-clockwise orientations,
$A^{(2)}_{j}$ and $A^{(3)}_{j}$ have clockwise orientations.

We get Riemann-Hilbert problem N0

  * 1.

$N(x,t,k)$ is meromorphic in
$\mathbb{C}\setminus(i\mathbb{R}\cup\mathbb{R}\cup A^{(1)}_{j}\cup
A^{(2)}_{j}\cup A^{(3)}_{j}\cup A^{(4)}_{j})$,

  * 2.

$N_{+}(x,t,k)=N_{-}(x,t,k)\begin{cases}T,z\in i\mathbb{R}\cup\mathbb{R},\\\
U^{1}_{j},k\in A^{(1)}_{j},\\\ U^{2}_{j},k\in A^{(2)}_{j},\\\ U^{3}_{j},k\in
A^{(3)}_{j},\\\ U^{4}_{j},k\in A^{(4)}_{j},\end{cases}$

  * 3.

$N(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where $U^{1}_{j}=\left(\begin{array}[]{cc}1&\frac{C_{j}{\rm
e}^{-{2i}\theta{(\kappa_{j})}}}{k-\kappa_{j}}\\\ 0&1\\\ \end{array}\right),$
$U^{2}_{j}=\left(\begin{array}[]{cc}1&0\\\ \frac{\bar{C}_{j}{\rm
e}^{{2i}\theta{(\bar{\kappa}_{j})}}}{k-\bar{\kappa}_{j}}&1\\\
\end{array}\right)$, $U^{3}_{j}=\left(\begin{array}[]{cc}1&0\\\
\frac{\bar{C}_{j}{\rm
e}^{{2i}\theta{(-\bar{\kappa}_{j})}}}{k+\bar{\kappa}_{j}}&1\\\
\end{array}\right)$ and $U^{4}_{j}=\left(\begin{array}[]{cc}1&\frac{C_{j}{\rm
e}^{-{2i}\theta{(-\kappa_{j})}}}{k+\kappa_{j}}\\\ 0&1\\\ \end{array}\right).$

The jump matrix $T$ is oscillatory in ${i}\mathbb{R}\cup\mathbb{R}$, so we
need to deform the jump contours of Riemann-Hilbert problem N0. In soliton
cases, the method for deforming the jump contours is the same as the non-
soliton cases. We need to ensure that the circles $A^{(m)}_{j}(m=1,2,3,4)$ do
not intersect with any jump contours. In soliton cases, we also divide the
$(x,t)$-plane into three regions:

  * 1.

Region 1: $-c_{1}\sqrt{t}\leq x\leq c_{1}\sqrt{t}$.

  * 2.

Region 2: $x<-c_{1}\sqrt{t}$.

  * 3.

Region 3: $x>c_{1}\sqrt{t}$.

####  4.3.1 Deformation in Region 1

In Region 1, we directly deform the jump contours at $k=0$. We define
$N^{(1)}(x,t,k)$ by

$N^{(1)}(x,t,k)=N(x,t,k)\begin{cases}I,&{\rm if}\
k\in\mathrm{E}_{3}\cup\mathrm{E}_{6}\cup\mathrm{E}_{9}\cup\mathrm{E}_{12},\\\
C^{-1},&{\rm if}\ k\in\mathrm{E}_{1}\cup\mathrm{E}_{7},\\\ A,&{\rm if}\
k\in\mathrm{E}_{2}\cup\mathrm{E}_{8},\\\ L,&{\rm if}\
k\in\mathrm{E}_{4}\cup\mathrm{E}_{10},\\\ U^{-1},&{\rm if}\
k\in\mathrm{E}_{5}\cup\mathrm{E}_{11},\end{cases}$ (63)

where $\mathrm{E}_{j}\ (j=1,2,...,12)$ is shown in Figure 10.

We get Riemann-Hilbert problem N1:

  * 1.

$N^{(1)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus(\Sigma_{1}\cup
A^{(0)}_{j})$,

  * 2.

$N^{(1)}_{+}(x,t,k)=N^{(1)}_{-}(x,t,k)\begin{cases}L,\quad&{\rm if}\
k\in\Sigma^{(3)}_{1}\cup\Sigma^{(7)}_{1},\\\ A,\quad&{\rm if}\
k\in\Sigma^{(4)}_{1}\cup\Sigma^{(8)}_{1},\\\ C,\quad&{\rm if}\
k\in\Sigma^{(1)}_{1}\cup\Sigma^{(5)}_{1},\\\ U,\quad&{\rm if}\
k\in\Sigma^{(2)}_{1}\cup\Sigma^{(6)}_{1},\\\ B,\quad&{\rm if}\
k\in\mathbb{R},\\\ U^{m}_{j},\quad&{\rm if}\ k\in
A^{(m)}_{j},(m=1,2,3,4),\end{cases}$

  * 3.

$N^{(1)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where $A^{(0)}_{j}=A^{(1)}_{j}\cup A^{(2)}_{j}\cup A^{(3)}_{j}\cup
A^{(4)}_{j}$, $\Sigma^{(j)}_{1}(j=1,2,...,8)$ is shown in Figure 12.

####  4.3.2 Deformation in Region 2

In Region 2, we need to deform the jump contours at the saddle points 0,
$k_{1}$ and $k_{2}$.

We define $N^{(2)}(x,t,k)$ by

$N^{(2)}(x,t,k)=N(x,t,k)\begin{cases}I,&{\rm if}\
k\in\mathrm{E}^{2}_{3}\cup\mathrm{E}^{2}_{6}\cup\mathrm{E}^{2}_{9}\cup\mathrm{E}^{2}_{12},\\\
C^{-1},&{\rm if}\
k\in\mathrm{E}^{2}_{1}\cup\mathrm{E}^{2}_{7}\cup\mathrm{E}^{2}_{14}\cup\mathrm{E}^{2}_{16},\\\
A,&{\rm if}\
k\in\mathrm{E}^{2}_{2}\cup\mathrm{E}^{2}_{8}\cup\mathrm{E}^{2}_{13}\cup\mathrm{E}^{2}_{15},\\\
L,&{\rm if}\ k\in\mathrm{E}^{2}_{4}\cup\mathrm{E}^{2}_{10},\\\ U^{-1},&{\rm
if}\ k\in\mathrm{E}^{2}_{5}\cup\mathrm{E}^{2}_{11},\end{cases}$ (64)

where the area $\mathrm{E}^{2}_{j}\ (j=1,2,...,16)$ is shown in Figure 11.

We get Riemann-Hilbert problem N2:

  * 1.

$N^{(2)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus(\Sigma_{2}\cup
A^{(0)}_{j})$,

  * 2.

$N^{(2)}_{+}(x,t,k)=N^{(2)}_{-}(x,t,k)\begin{cases}A,\quad&{\rm if}\
k\in\Sigma^{(1)}_{2}\cup\Sigma^{(5)}_{2},\\\ C,\quad&{\rm if}\
k\in\Sigma^{(4)}_{2}\cup\Sigma^{(8)}_{2},\\\ U,\quad&{\rm if}\
k\in\Sigma^{(2)}_{2}\cup\Sigma^{(6)}_{2}\cup\Sigma^{(13)}_{1}\cup\Sigma^{(14)}_{2}\cup\Sigma^{(15)}_{2}\cup\Sigma^{(16)}_{2},\\\
L,\quad&{\rm if}\
k\in\Sigma^{(3)}_{2}\cup\Sigma^{(7)}_{2}\cup\Sigma^{(9)}_{1}\cup\Sigma^{(10)}_{2}\cup\Sigma^{(11)}_{2}\cup\Sigma^{(12)}_{2},\\\
B,\quad&{\rm if}\ k\in\cup(-\infty,k_{2})\cup(k_{1},+\infty),\\\
U^{m}_{j},\quad&{\rm if}\ k\in A^{(m)}_{j},(m=1,2,3,4),\end{cases}$

  * 3.

$N^{(2)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where the jump contour $\Sigma^{(j)}_{2}(j=1,2,...,16)$ is shown in Figure 12.

####  4.3.3 Deformation in Region 3

In Region 3, we need to deform the jump contours at the saddle points $k=0$,
$k=k_{1}$ and $k=k_{2}$.

We define $N^{(3)}(x,t,k)$ by

$N^{(3)}(x,t,k)=N(x,t,k)\begin{cases}I,&{\rm if}\
k\in\mathrm{E}^{3}_{3}\cup\mathrm{E}^{3}_{6}\cup\mathrm{E}^{3}_{9}\cup\mathrm{E}^{3}_{12},\\\
C^{-1},&{\rm if}\
k\in\mathrm{E}^{3}_{1}\cup\mathrm{E}^{3}_{7}\cup\mathrm{E}^{3}_{14}\cup\mathrm{E}^{3}_{16},\\\
A,&{\rm if}\
k\in\mathrm{E}^{3}_{2}\cup\mathrm{E}^{3}_{8}\cup\mathrm{E}^{3}_{13}\cup\mathrm{E}^{3}_{15},\\\
L,&{\rm if}\ k\in\mathrm{E}^{3}_{4}\cup\mathrm{E}^{3}_{10},\\\ U^{-1},&{\rm
if}\ k\in\mathrm{E}^{3}_{5}\cup\mathrm{E}^{3}_{11},\\\ \end{cases}$ (65)

where the area $\mathrm{E}^{3}_{j}\ (j=1,2,...,16)$ is shown in Figure 13.

We get Riemann-Hilbert problem N3:

  * 1.

$N^{(3)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus(\Sigma_{3}\cup
A^{(0)}_{j})$,
$\Sigma_{3}=\Sigma^{(1)}_{3}\cup\cdots\cup\Sigma^{(16)}_{3}\cup(-\infty{i},k_{2})\cup(\infty{i},k_{1})$,

  * 2.

$N^{(3)}_{+}(x,t,k)=N^{(3)}_{-}(x,t,k)\begin{cases}L,&{\rm if}\
k\in\Sigma^{(3)}_{3}\cup\Sigma^{(7)}_{3}\cup\Sigma^{(9)}_{3}\cup\Sigma^{(10)}_{3}\cup\Sigma^{(11)}_{3}\cup\Sigma^{(12)}_{3},\\\
U,&{\rm if}\
k\in\Sigma^{(4)}_{3}\cup\Sigma^{(8)}_{3}\cup\Sigma^{(13)}_{3}\cup\Sigma^{(14)}_{3}\cup\Sigma^{(15)}_{3}\cup\Sigma^{(16)}_{3},\\\
C,&{\rm if}\ k\in\Sigma^{(1)}_{3}\cup\Sigma^{(5)}_{3},\\\ A,&{\rm if}\
k\in\Sigma^{(2)}_{3}\cup\Sigma^{(6)}_{3},\\\ B,&{\rm if}\
k\in(-\infty{i},k_{2})\cup(\infty{i},k_{1}),\\\ U^{m}_{j},\quad&{\rm if}\ k\in
A^{(m)}_{j},(m=1,2,3,4),\end{cases}$

  * 3.

$N^{(3)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where the jump contour $\Sigma^{(j)}_{3}(j=1,2,...,16)$ is shown in Figure 14.

###  4.4 Extra deformations for the RHP (soliton cases)

When the initial value for the DNLS equation can evolve into the soliton, some
extra deformations are needed.

######  Remark 8.

For some $x$ and $t$, $|{\rm e}^{2i\theta{(\kappa_{j})}}|$ and $|{\rm
e}^{-2i\theta{(\bar{\kappa}_{j})}}|$ may become large, an extra deformation is
used to deal with $U^{m}_{j}$.

We define $N^{(4)}(x,t,k)$ by

$N^{(4)}(x,t,k)=N^{(m_{1})}(x,t,k)\begin{cases}\left(\begin{array}[]{cc}1&\frac{k-\kappa_{j}}{C_{j}{\rm
e}^{-{2i}\theta{(\kappa_{j})}}}\\\ -\frac{C_{j}{\rm
e}^{-{2i}\theta{(\kappa_{j})}}}{k-\kappa_{j}}&0\\\
\end{array}\right)K,&{|k-\kappa_{j}|<\varepsilon},\\\
\left(\begin{array}[]{cc}1&\frac{k+\kappa_{j}}{C_{j}{\rm
e}^{2i\theta{(-\kappa_{j})}}}\\\ -\frac{C_{j}{\rm
e}^{2i\theta{(-\kappa_{j})}}}{k+\kappa_{j}}&0\\\
\end{array}\right)K,&{|k+\kappa_{j}|<\varepsilon},\\\
\left(\begin{array}[]{cc}0&\frac{\bar{C}_{j}{\rm
e}^{2i\theta{(\bar{\kappa}_{j})}}}{k-\bar{\kappa}_{j}}\\\
-\frac{k-\bar{\kappa}_{j}}{\bar{C}_{j}{\rm
e}^{2i\theta{(\bar{\kappa}_{j})}}}&1\\\
\end{array}\right)K,&{|k-\bar{\kappa}_{j}|<\varepsilon},\\\
\left(\begin{array}[]{cc}0&\frac{\bar{C}_{j}{\rm
e}^{-2i\theta{(-\bar{\kappa}_{j})}}}{k+\bar{\kappa}_{j}}\\\
-\frac{k+\bar{\kappa}_{j}}{\bar{C}_{j}{\rm
e}^{-2i\theta{(-\bar{\kappa}_{j})}}}&1\\\
\end{array}\right)K,&{|k+\bar{\kappa}_{j}|<\varepsilon},\\\ K,&k\in{\rm
others},\end{cases}$ (66)

where
$K=\left(\begin{array}[]{cc}\frac{(z-\kappa_{j})(z+\kappa_{j})}{(z-\bar{\kappa}_{j})(z+\bar{\kappa}_{j})}&0\\\
0&\frac{(z-\bar{\kappa}_{j})(z+\bar{\kappa}_{j})}{(z-\kappa_{j})(z+\kappa_{j})}\\\
\end{array}\right)$ and $m_{1}$ is equal to 1, 2, or 3.

In Region 1, we get Riemann-Hilbert problem N4:

  * 1.

$N^{(4)}(x,t,k)$ is meromorphic in $\mathbb{C}\setminus(\Sigma_{1}\cup
A^{(0)}_{j})$,

  * 2.

$N^{(4)}_{+}(x,t,k)=N^{(4)}_{-}(x,t,k)\begin{cases}K^{-1}LK,\quad&{\rm if}\
k\in\Sigma^{(3)}_{1}\cup\Sigma^{(7)}_{1},\\\ K^{-1}AK,\quad&{\rm if}\
k\in\Sigma^{(4)}_{1}\cup\Sigma^{(8)}_{1},\\\ K^{-1}CK,\quad&{\rm if}\
k\in\Sigma^{(1)}_{1}\cup\Sigma^{(5)}_{1},\\\ K^{-1}UK,\quad&{\rm if}\
k\in\Sigma^{(2)}_{1}\cup\Sigma^{(6)}_{1},\\\ K^{-1}BK,\quad&{\rm if}\ k\in
i\mathbb{R},\\\ K^{-1}V^{m}_{j}K,\quad&{\rm if}\ k\in
A^{(m)}_{j},(m=1,2,3,4),\end{cases}$

  * 3.

$N^{(4)}(x,t,k)\rightarrow I$, $k\rightarrow\infty$,

where $V^{1}_{j}=\left(\begin{array}[]{cc}1&0\\\ \frac{k-\kappa_{j}}{C_{j}{\rm
e}^{2i\theta{(\kappa_{j})}}}&1\\\ \end{array}\right),$
$V^{2}_{j}=\left(\begin{array}[]{cc}1&\frac{k-\bar{\kappa}_{j}}{\bar{C}_{j}{\rm
e}^{-2i\theta{(\bar{\kappa}_{j})}}}\\\ 0&1\\\ \end{array}\right)$,
$V^{3}_{j}=\left(\begin{array}[]{cc}1&\frac{k+\bar{\kappa}_{j}}{\bar{C}_{j}{\rm
e}^{-2i\theta{(-\bar{\kappa}_{j})}}}\\\ 0&1\\\ \end{array}\right)$ and
$V^{4}_{j}=\left(\begin{array}[]{cc}1&0\\\ \frac{k+\kappa_{j}}{C_{j}{\rm
e}^{2i\theta{(-\kappa_{j})}}}&1\\\ \end{array}\right)$.

The Riemann-Hilbert problem of Region 2 or Region 3 is similar to Region 1,
and we do not give its specific expressions.

######  Remark 9.

When $t$ is large, we also need to remove the jump matrix $B$, which is the
same as methods of section 4.2.

######  Remark 10.

When choosing the jump contours, we need to consider the analytical property
of the scattering data, which have been widely applied in previous NIST[39,
43].

###  4.5 Numerical solution for the RHP

After deforming the Riemann-Hilbert problem, it is necessary to solve the
deformed problem using a numerical method. The Chebyshev collocation method,
developed by Olver, is commonly employed for solving the Riemann-Hilbert
problem numerically[44, 45, 46].

Give $x$ and $t$, and we get the following RHP,

$\begin{array}[]{lr}\Phi^{+}(k)=\Phi^{-}(k)G(k),k\in\Gamma,\\\
\Phi(\infty)=I,\end{array}$ (67)

where $\Gamma$=$\Gamma_{1}\cup\Gamma_{2}\cup\cdots\cup\Gamma_{n_{1}}$, and
$\Gamma_{1},\ldots,\Gamma_{n_{1}}$ are non-self-intersecting arcs.

$\Phi(k)$ can be written into the following form

$\Phi(k)=I+\mathcal{C}_{\Gamma}V(k),$ (68)

where $\mathcal{C}_{\Gamma}$ represents the Cauchy transformation

$\mathcal{C}_{\Gamma}V(z)=\frac{1}{2i\pi}\int_{\Gamma}\frac{V(\xi)}{\xi-z}{\rm
d}\xi.$

According to the Plemelj formula, we get

$\displaystyle\Phi^{+}(k)-\Phi^{-}(k)=V(k),$ (69)
$\displaystyle\Phi^{+}(k)+\Phi^{-}(k)=-i\mathcal{H}V(k),$ (70)

where $\mathcal{H}$ is the Hilbert transform.

We need to search the suitable $\rm M\ddot{o}bius$ transform
$M_{j}(j=1,2,\ldots,n_{1})$ which mapping $\Gamma_{j}$ into unit interval
$\mathbb{I}$.

Eq. (68) is written into

$\Phi(k)=I+\mathcal{C}_{\Gamma_{1}}V_{1}(k)+\mathcal{C}_{\Gamma_{2}}V_{2}(k)+\ldots+\mathcal{C}_{\Gamma_{n_{1}}}V_{n_{1}}(k).$
(71)

For $V_{j}(k)=\left(\begin{array}[]{cc}V^{(11)}_{j}(k)&V^{(12)}_{j}(k)\\\
V^{(21)}_{j}(k)&V^{(22)}_{j}(k)\\\ \end{array}\right),$ we appropriate it by
the Chebyshev polynomials

$\begin{cases}V^{(1,1)}_{j}(k)=T(M_{j}^{-1}(k))\mathcal{F}v^{(11)}_{j},\\\
V^{(1,2)}_{j}(k)=T(M_{j}^{-1}(k))\mathcal{F}v^{(12)}_{j},\\\
V^{(2,1)}_{j}(k)=T(M_{j}^{-1}(k))\mathcal{F}v^{(21)}_{j},\\\
V^{(2,2)}_{j}(k)=T(M_{j}^{-1}(k))\mathcal{F}v^{(22)}_{j},\end{cases}$ (72)

where $v^{(\ )}_{j}$ is the unknown coefficient.

The unknown coefficient $v^{(\ )}_{j}$ can be calculated by Eq. (67).
According to Eq. (68), we get

$\displaystyle\lim_{k\rightarrow\infty}k\Phi(k)_{(1,2)}=\lim_{k\rightarrow\infty}\frac{1}{2i\pi}\int_{\Gamma}\frac{k}{\xi-k}V^{(12)}(\xi){\rm
d}\xi$ (73)
$\displaystyle=-\frac{1}{2i\pi}\Big{[}\int_{\Gamma_{1}}V_{1}^{(12)}(\xi){\rm
d}\xi+\int_{\Gamma_{2}}V_{2}^{(12)}(\xi){\rm
d}\xi+\cdots+\int_{\Gamma_{k_{1}}}V_{k_{1}}^{(12)}(\xi){\rm d}\xi\Big{]}.$

We only summarize the key steps of the Chebyshev collocation method and do not
give specific details of the method. The specific details of the Chebyshev
collocation method can be found in Olver’s paper[44, 45]. We give the
following remarks to help readers to implement this method better.

######  Remark 11.

The keys for the Chebyshev collocation method are calculating the Cauchy
transformation of the Chebyshev polynomials and constructing Cauchy matrices,
which allows us to turn (67) into a linear system with unknown coefficient
$v^{(\ )}_{j}$.

######  Remark 12.

There are some relationships between certain Cauchy matrices, perhaps they are
completely equal, or their imaginary parts are equal and their real parts are
opposite. These relationships are determined by the value of $k$ and the
jumping contour $\Gamma_{j}(j=1,2,...,k_{1})$. Using these relationships can
reduce computational costs.

##  5 Numerical Results and Comparison

In this section, the NIST is used to solve the DNLS equation, and we compare
the NIST with the traditional numerical method. It is necessary to introduce
traditional numerical methods before presenting numerical results. One of the
most effective methods for solving the DNLS equation is the Fourier spectral
method (FSM). In the FSM, the $x$-space is truncated into $[-L,L]$ and the
periodic boundary $q(L,t)=q(-L,t)$ is applied. Doing the Fourier
transformation for the the DNLSE $\rm III$, we get

$\mathbb{F}[q]_{t}=-ik^{2}\mathbb{F}[q]+\mathbb{F}\Big{[}q^{2}\mathbb{F}^{-1}\big{[}ik\mathbb{F}[\bar{q}]\big{]}\Big{]}+\frac{i}{2}\mathbb{F}[|q|^{4}q],$
(74)

where $\mathbb{F}$ represents the Fourier transformation and $\mathbb{F}^{-1}$
represents the inverse Fourier transformation. We rewrite Eq. (74) into

$\frac{{\rm d}{\rm e}^{ik^{2}t}\mathbb{F}[q]}{{\rm d}t}={\rm
e}^{ik^{2}t}\mathbb{F}\Big{[}q^{2}\mathbb{F}^{-1}\big{[}ik\mathbb{F}[\bar{q}]\big{]}\Big{]}+{\rm
e}^{ik^{2}t}\frac{i}{2}\mathbb{F}[|q|^{4}q].$ (75)

Eq. (75) can be solved by the finite difference method (FDM) or the Runge-
Kutta method (RKM). The solution for the DNLS equation is given by the inverse
Fourier transformation of $\mathbb{F}[q]$.

There is no doubt that the FSM is efficient in the short time, but the FSM
will lose accuracy once the wave arrives the calculated boundary. The error
generated by the FSM will increase over time which will lead to numerical
convergence issues. The FSM is used to calculate the DNLS equation with
$\exp(-x^{2}){\rm sech}(x)$ initial profile, and the calculated results are
shown in Figure 17. In Figure 17(a)-(c), the calculated results are normal.
When the time reaches 15, there is a significant oscillation in the numerical
results, and the oscillation becomes more severe as time increases. When the
time reaches 30, the FSM cannot effectively calculate the solution of the DNLS
equation. These phenomena can be clearly observed from Figure 17 (d)-(h).

(a) $t=0$ (b) $t=5$ (c) $t=10$ (d) $t=15$ (e) $t=20$ (f) $t=25$ (g) $t=30$ (h)
$t=35$

Figure 17: The numerical solution for the $\exp(-x^{2}){\rm sech}(x)$ initial
profile calculated by the FSM.

There is no discrete eigenvalue for the $\exp(-x^{2}){\rm sech}(x)$ initial
profile, so we don’t need to calculate the discrete eigenvalues and the
normalized constant. The NIST is used to calculate the evolution for
$\exp(-x^{2}){\rm sech}(x)$ initial profile in the short time, and the
calculated results are marked as $q_{\rm nist1}(x,t)$.

In the short time, we don’t need to deform the Riemann-Hilbert problem. The
evolution of $\exp(-x^{2}){\rm sech}(x)$ initial profile calculated by the
NIST is shown in Figure 18. Figure 18(a) shows the density for the calculated
result; Figure 18(b) shows the real part for the calculated result; Figure
18(c) shows the imaginary part for the calculated result. From Figure 18, we
can clearly see the spatiotemporal evolution of the initial value in the short
time.

(a) Density (b) Real part (c) Imaginary part

Figure 18: The spatiotemporal evolution of the $\exp(-x^{2}){\rm sech}(x)$
initial profile calculated by the NIST in the short time.

Figure 19 shows the numerical result $q_{\rm nist1}(x,t)$ for the
$\exp(-x^{2}){\rm sech}(x)$ initial profile at $t=0$, $t=1$ and $t=2$. The
blue line is the density for $q_{\rm nist1}(x,t)$, the red line is the real
part for $q_{{\rm nist1}}(x,t)$ and the green line is the imaginary part for
$q_{\rm nist1}(x,t)$.

(a) $t=0$ (b) $t=1$ (c) $t=2$

Figure 19: The numerical solution $q_{\rm nist1}(x,t)$ for $\exp(-x^{2}){\rm
sech}(x)$ initial profile calculated by the NIST. Blue line: the density for
$q_{\rm nist1}(x,t)$. Red line: the real part for $q_{\rm nist1}(x,t)$. Green
line: the imaginary part for $q_{\rm nist1}(x,t)$.

When time is small, the calculation results of these two methods match
theoretically. To verify this conclusion, we used two methods to solve the
numerical solution for $\exp(-x^{2}){\rm sech}(x)$ initial profile in the
short time, and the density for the numerical solutions are shown in Figure
20, Figure 21 and Figure 22. Figure 20 shows the density of the numerical
solutions for $\exp(-x^{2}){\rm sech}(x)$ initial profile. Figure 21 shows the
real part of the numerical solutions for $\exp(-x^{2}){\rm sech}(x)$ initial
profile. Figure 22 shows the imaginary part of the numerical solutions for
$\exp(-x^{2}){\rm sech}(x)$ initial profile. The absolute errors calculated by
the two methods are shown in Figure 23, and ‘Error’ is defined by

${\rm Error}=\big{|}q_{\rm fsm1}(x,t)-q_{\rm nist1}(x,t)\big{|},$

where $q_{\rm fsm1}(x,t)$ is the numerical solution for $\exp(-x^{2}){\rm
sech}(x)$ initial profile calculated by the FSM. It is obvious that the two
method are highly consistent, so both of them are effective.

(a) $t=0.5$ (b) $t=1$ (c) $t=1.5$

Figure 20: The density of the numerical solution for $\exp(-x^{2}){\rm
sech}(x)$ initial profile. Red line: the numerical result calculated by the
FSM. Green line: the numerical result calculated by the NIST.

(a) $t=0.5$ (b) $t=1$ (c) $t=1.5$

Figure 21: The real part of the numerical solution for $\exp(-x^{2}){\rm
sech}(x)$ initial profile. Red line: the numerical result calculated by the
FSM. Green line: the numerical result calculated by the NIST.

(a) $t=0.5$ (b) $t=1$ (c) $t=1.5$

Figure 22: The imaginary part of the numerical solution for $\exp(-x^{2}){\rm
sech}(x)$ initial profile. Red line: the numerical result calculated by the
FSM. Green line: the numerical result calculated by the NIST.

(a) $t=0.5$ (b) $t=1$ (c) $t=1.5$

Figure 23: The absolute errors between the NIST and the FSM at different time.

The NIST is used to calculate $q_{\rm nist1}(0,t)$ in $t\in[0,7]$, the
calculated result is shown in Figure 24(a). When $t>2$, $q_{\rm nist1}(0,t)$
changes from flat to oscillating. As time increases, the oscillation for
$q_{\rm nist1}(0,t)$ becomes more severe. We have mentioned the reasons for
the above phenomenon in section 4. As time increases, the oscillation of the
jump matrix becomes more severe. For highly oscillatory functions, it is
difficult to achieve their uniform approximation, and we can not solve the
original Riemann-Hilbert problem with high accuracy. So when $t\geq 2$, we
will deform the original Riemann-Hilbert problem to ensure numerical accuracy.
After deformation, the NIST is used to calculate $q_{\rm nist1}(0,t)$, the
calculated result is shown in Figure 24(b). After deformation, the oscillation
disappears and the NIST becomes effective.

Figure 24: Subfigure(a): the numerical results of $q_{\rm nist1}(0,t)$ initial
profile without deformation. Subfigure(b): the numerical results for $q_{\rm
nist1}(0,t)$ initial profile after deformation. Blue line: the density for
$q_{\rm nist1}(0,t)$. Red line: the real part for $q_{\rm nist1}(0,t)$. Green
line: the imaginary part for $q_{\rm nist1}(0,t)$. Black line: $t=2$.

The NIST (Nonlinear Inverse Scattering Transform) method is utilized to
calculate the long-term evolution of the DNLS equation. The calculated results
are visualized in Figure 25. The reliability of these calculations highlights
the effectiveness of the NIST approach. In comparison, traditional numerical
methods like the finite difference method and the FSM struggle to achieve
high-precision numerical calculations over extended periods.

(a) $t=100$ (b) $t=200$ (c) $t=300$

Figure 25: The long-time evolution $q_{\rm nist1}(0,t)$ for $\exp(-x^{2}){\rm
sech}(x)$ initial profile calculated by the NIST. Blue line: the density for
$q_{\rm nist1}(0,t)$. Red line: the real part for $q_{\rm nist1}(0,t)$. Green
line: the imaginary part for $q_{\rm nist1}(0,t)$.

##  6 Conclusions and Discussions

In summary, we have successfully demonstrated the effectiveness of the
numerical inverse scattering transform (NIST) for the DNLS equation with zero-
boundary condition. Through numerical examples, we have validated the method
and compared it with the FSM. The DNLS equation exhibits a continuous spectrum
of $i\mathbb{R}\cup\mathbb{R}$, with three saddle points excluding $x=0$.
These characteristics posed challenges in applying the numerical inverse
scattering transform, but we have successfully addressed them. By dividing the
$(x,t)$-plane into three regions, we devised specific approaches to deform the
jump contours accordingly. To ensure the accuracy of our calculations, we have
mathematically established that the deformed jump matrices remain uniformly
bounded. This provides a guarantee for the reliability of our results.
Moreover, we have established a feasible framework for computing the
scattering data of the DNLS system.

During our calculations, we observed that the NIST method imposes strict
limitations on the space of initial values. Specifically, the initial value
must belong to the Schwartz space $\mathcal{S}(\mathbb{R})$ or even spaces
with faster decay. This restriction narrows down the range of valid initial
values for the NIST. A significant focus of future research will involve
expanding the space of initial values to encompass a more general class within
the NIST framework. By reducing the NIST’s dependence on the decay rate of the
initial value, we can enhance its applicability and accommodate a wider range
of initial conditions. Addressing and overcoming these limitations will be a
key area of focus in our future work. Such advancements will help advance the
capabilities of the numerical inverse scattering transform for the DNLS
equation and broaden its range of applications.

## Acknowledgment

This work is supported by National Science Foundation of China (52171251,
U2106225, 52231011) and Dalian Science and Technology Innovation Fund
(2022JJ12GX036).

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Data availability

Data will be made available on request.

## References

## References

  * [1] Han S.H., and Park Q.H. (2011) Effect of self-steepening on optical solitons in a continuous wave background. Phys. Rev. E 83(6) 066601
  * [2] He J.S., Charalampidis E.G., Kevrekidis P.G. and Frantzeskakis D.J. (2014) Rogue waves in nonlinear Schr$\rm\ddot{o}$dinger models with variable coefficients: application to Bose-Einstein condensates. Phys. Lett. A 378(5-6) 577-583
  * [3] Kodama Y. (1985) Optical solitons in a monomode fiber. J. Stat. Phys. 39 597-614
  * [4] de Oliveira J.R. and de Moura M.A. (1998) Analytical solution for the modified nonlinear Schr$\rm\ddot{o}$dinger equation describing optical shock formation. Phys. Rev. E 57(4) 4751
  * [5] Kaup D.J. and Newell A.C. (1978) An exact solution for a derivative nonlinear Schr$\rm\ddot{o}$dinger equation. J. Math. Phys. 19(4) 798-801
  * [6] Chen H.H., Lee Y.C. and Liu C.S. (1979) Integrability of nonlinear Hamiltonian systems by inverse scattering method. Phys. Scr. 20(3-4) 490
  * [7] Gerdzhikov V.S., Ivanov M.I. and Kulish P.P. (1980) Quadratic bundle and nonlinear equations. Theor. Math. Phys. 44 784-795
  * [8] Kundu A. (1987) Exact solutions to higher-order nonlinear equations through gauge transformation. Physica D 25(1-3) 399-406
  * [9] Wadati M. and Sogo K. (1983) Gauge transformations in soliton theory. J. Phys. Soc. Japan 52(2) 394-398
  * [10] Clarkson P.A. and Cosgrove C.M. (1987) Painleve analysis of the non-linear Schr$\rm\ddot{o}$dinger family of equations. J. Phys. A Math. Theor. 20(8) 2003
  * [11] Gardner C.S., Greene J.M., Kruskal M.D. and Miura R.M. (1967) Method for solving Korteweg-de Vries equation. Phys. Rev. Lett. 19 1095-1097
  * [12] Zhou G.Q. and Huang N.N. (2007) An $N$-soliton solution to the DNLS equation based on revised inverse scattering transform. J. Phys. A Math. Theor. 40(45) 13607
  * [13] Chen X.J. and Lam W.K. (2004) Inverse scattering transform for the derivative nonlinear schr$\rm\ddot{o}$dinger equation with nonvanishing boundary conditions. Phys. Rev. E 69 066604
  * [14] Guo B.L. and Ling. L.M. (2012) Riemann-Hilbert approach and $N$-soliton formula for coupled derivative Schr$\rm\ddot{o}$dinger equation. J. Math. Phys. 53(7) 073506
  * [15] Zhang G. and Yan Z. (2020) The derivative nonlinear Schr$\rm\ddot{o}$dinger equation with zero/nonzero boundary conditions: inverse scattering transforms and $N$-double-pole solutions. J. Nonlinear Sci. 30(6) 3089-3127
  * [16] Ma X.X. and Kuang Y.H. (2022) Inverse scattering transform for a nonlocal derivative nonlinear Schr$\rm\ddot{o}$dinger equation. Theor. Math. Phys. 210(1) 31-45
  * [17] Liu N., Xuan Z. and Sun J. (2022) Triple-pole soliton solutions of the derivative nonlinear Schr$\rm\ddot{o}$dinger equation via inverse scattering transform. Appl. Math. Lett. 125 107741
  * [18] Fan E. (2001) Explicit $N$-fold Darboux transformations and soliton solutions for nonlinear derivative Schr$\rm\ddot{o}$dinger equations. Commun. Theor. Phys. 35(6) 651
  * [19] Guo B., Ling L. and Liu Q.P. (2013) High-order solutions and generalized Darboux transformations of derivative nonlinear Schr$\rm\ddot{o}$dinger equations. Stud. Appl. Math. 130(4) 317-344
  * [20] Chen J. and Pelinovsky D.E. (2021) Rogue waves on the background of periodic standing waves in the derivative nonlinear Schr$\rm\ddot{o}$dinger equation. Phys. Rev. E 103(6) 062206
  * [21] Yang B., Chen J. and Yang J. (2020) Rogue waves in the generalized derivative nonlinear Schr$\rm\ddot{o}$dinger equations. J. Nonlinear Sci. 30(6) 3027-3056
  * [22] Liu J.Q., Perry P.A. and Sulem C. (2016) Global existence for the derivative nonlinear Schr$\rm\ddot{o}$dinger equation by the method of inverse scattering. Commun. Partial Differ. Equ. 41(11) 1692-1760
  * [23] Pelinovsky D.E., Saalmann A. and Shimabukuro Y. (2017) The derivative NLS equation: global existence with solitons. Dyn. Partial Differ. Equ. 14 (3) 271-294
  * [24] Pelinovsky D.E. and Shimabukuro Y. (2018) Existence of global solutions to the derivative nls equation with the inverse scattering transform method. Int. Math. Res. Not. 2018(18) 5663-5728
  * [25] Jenkins R., Liu J., Perry P. and Sulem C. (2020) Global existence for the derivative nonlinear Schr$\rm\ddot{o}$dinger equation with arbitrary spectral singularities. Anal. PDE 13(5) 1539-1578
  * [26] Deift P. and Zhou X. (1993). A steepest descent method for oscillatory Riemann–Hilbert problems. Asymptotics for the MKdV equation. Ann. Math. 137(2) 295-368
  * [27] Bilman D. and Miller P.D. (2019). A robust inverse scattering transform for the focusing nonlinear Schr$\rm\ddot{o}$dinger equation. Commun. Pure. Appl. Math. 72(8) 1722-1805
  * [28] Xu J., Fan E. and Chen Y. (2013) Long-time asymptotic for the derivative nonlinear Schr$\rm\ddot{o}$dinger equation with step-like initial value. Math. Phys. Anal. Geom. 16 253-288
  * [29] Arruda L.K. and Lenells, J. (2017) Long-time asymptotics for the derivative nonlinear Schr$\rm\ddot{o}$dinger equation on the half-line. Nonlinearity 30(11) 4141
  * [30] Jenkins R., Liu J., Perry P. and Sulem C. (2018) Soliton resolution for the derivative nonlinear schr$\rm\ddot{o}$dinger equation. Commun. Math. Phys. 363 1003-1049
  * [31] Tian S.F. (2021) Riemann-Hilbert problem to a generalized derivative nonlinear Schr$\rm\ddot{o}$dinger equation: Long-time asymptotic behavior. Sci. Sin. Math. 51 1-38
  * [32] Yang Y. and Fan E. (2022) Long-time Asymptotic Behavior for the Derivative Schr$\rm\ddot{o}$dinger Equation with Finite Density Type Initial Data. Chin. Ann. Math. Ser. B 43(6) 893-948
  * [33] Yousefi M.I. and Kschischang F.R. (2014) Information transmission using the nonlinear Fourier transform, Part II: Numerical methods. IEEE Trans. Inf. Theory 60(7) 4329-4345
  * [34] Boyd J.P. (2001) Chebyshev and Fourier spectral methods. Courier Corporation.
  * [35] Yang J.K. (2010) Nonlinear waves in integrable and nonintegrable systems. Society for Industrial and Applied Mathematics
  * [36] Jie S., Tao T. and Wang L.L (2011) Spectral Methods: Algorithms, Analysis and Applications.
  * [37] Trogdon T. (2021) Scattering and inverse scattering for the AKNS system: A rational function approach. Stud. Appl. Math. 147(4) 1443-1480.
  * [38] Trogdon T., Olver S. and Deconinck B. (2012) Numerical inverse scattering for the Korteweg-de Vries and modified Korteweg-de Vries equations. Physica D 241(11) 1003-1025
  * [39] Trogdon T. and Olver S. (2013) Numerical inverse scattering for the focusing and defocusing nonlinear Schr$\rm\ddot{o}$dinger equations. Proc. Math. Phys. Eng. Sci. 469(2149) 20120330
  * [40] Bilman D. and Trogdon T. (2017) Numerical inverse scattering for the Toda lattice. Commun. Math. Phys. 352 805-879
  * [41] Deconinck B., Trogdon T. and Yang X. (2019) Numerical inverse scattering for the sine-Gordon equation. Physica D 399 159-172
  * [42] Cui S. and Wang Z. (2023) Numerical inverse scattering transform for the focusing and defocusing Kundu-Eckhaus equations. Physica D 454 133838
  * [43] Cui S. and Wang Z. (2023) Efficient method for calculating the eigenvalues of the Zakharov-Shabat system. Chin. Phys. B (In Press) DOI: 10.1088/1674-1056/acd686
  * [44] Olver S. (2011) Numerical solution of Riemann-Hilbert problems: $\rm Painlev\acute{e}$ II. Found. Comut. Math. 11(2) 153-179
  * [45] Olver S. (2012) A general framework for solving Riemann-Hilbert problems numerically. Numer. Math. 122(2) 305-340
  * [46] Olver S. (2011) Computing the Hilbert transform and its inverse. Math. Comput. 80(275) 1745-1767