PUBLICATIONS DE L’INSTITUT MATHÉMATIQUE Nouvelle série, tome 96 (110) (2014), 125–141 DOI: 10.2298/PIM1410125F

EMERGING PROBLEMS IN APPROXIMATION THEORY FOR THE NUMERICAL SOLUTION

OF THE NONLINEAR SCHRÖDINGER EQUATION

L. Fermo, C. Van der Mee, and S. Seatzu

Abstract. We present some open problems pertaining to the approximation theory involved in the solution of the Nonlinear Schrödinger (NLS) equation. For this important equation, any Initial Value Problem (IVP) can be theoretically solved by the Inverse Scattering Transform (IST) technique whose main steps involve the solution of Volterra equations with structured kernels on unbounded domains, the solution of Fredholm integral equations and the identification of coefficients and parameters of monomial-exponential sums. The aim of the paper is twofold: propose a method for solving the above mentioned problems under particular hypothesis; arise interest in the issues illustrated to achieve an effective method for solving the problem under more general assumptions

1. Introduction

The class of Nonlinear Partial Differential Equations (NPDEs) of integrable type is important in mathematics as in several applicative areas of physics, biology and engineering [1–3, 9, 12]. For this special class of NPDEs, the nonlinear Schrödinger (NLS) equation, which arises in modeling electromagnetic waves in optical fibers as well as waves on the surface of deep water, has a special role in mathematics [4, 11, 13, 15]. Firstly, we recall that the NLS equation is expressed as

(1.1)  iu_{t} + u_{xx} ± 2|u|^{2}u = 0, u = u(x, t), x ∈ ℝ, t > 0

where i denotes the imaginary unit, the subscripts x and t denote the partial derivatives with respect to position x and time t and the sign ± depends on the symmetry properties of the model we are addressing. In particular, the plus sign appears in the focusing case and the minus sign in the defocusing case, which represent the two most important situations.

2010 Mathematics Subject Classification: 41A46, 65R20, 35P25. The research was partially supported by the National Group for Scientific Computing and

the National Group for Mathematical Physics of the National Institute for Advanced Mathematics (INdAM-GNCS, INdAM-GNFM) and by Autonomous Region of Sardinia under grant L.R.7/2007 “Promozione della Regione Scientifica e della Innovazione Tecnologica in Sardegna”.

We are interested in the initial value problem (IVP) for the NLS, that is in considering (1.1), given the initial solution

u_{0}(x) = u(0, x), x ∈ ℝ, u_{0} ∈ L^{1}(ℝ).

Following the path of the IST [2], its solution can be obtained by solving, in order, the following three problems:

(a) determine the initial scattering data, given its initial solution;

(b) propagate the initial scattering data in time;

(c) solve two systems of integral equations whose kernels codify the initial scattering data evolved in time.

From the numerical point of view, the problem of most interest is the first one, as effective methods to solve the two other problems have been developed recently, under the assumption that the initial scattering data are known [5]. Some attempts to approximate the solution of (1.1), following the path of the IST, have been made in [6] under assumptions too restrictive in many applicative areas.

Let us now illustrate the organization of the paper. In Section 2 we discuss the Zakharov–Shabat (ZS) system, which gives a complete characterization of the scattering data associated to the NLS we want to compute, that is the transmission coefficient, the reflection coefficients (from the left and from the right), the bound states and the norming constants. Section 3 is devoted to the introduction and characterization of the auxiliary functions whose approximation is basic to evaluating all of the scattering data mentioned above. In Section 4 we introduce the initial Marchenko kernels, which codify the scattering data and that can be computed by solving Volterra integral equations. In Section 5 we propose a numerical method to compute the scattering data in the reflectionless case. In Section 6 we present the numerical results which confirm the effectiveness of the method in this particular case. Section 7 is devoted to conclusions and perspectives.

2. Initial scattering data

The characterization of the initial scattering data is based on the spectral analysis of the ZS system associated to the NLS equation, which in turn is represented by an ordinary differential equation of first order [2].

In fact, assuming that u_{0} ∈ L^{1}(ℝ), it can be expressed in the following way:

(2.1)  iJ\frac{∂X}{∂x} (λ, x) − V(x)X(λ, x) = λX(λ, x), x ∈ ℝ

where λ ∈ ℂ is a spectral parameter,

J = ( \begin{matrix}1 & 0 \\ 0 & −1 \end{matrix} ) , V(x) = ( \begin{matrix}0 & iu_{0}(x) \\ ±iū_{0}(x) & 0 \end{matrix}) .

Here the bar is used to denote complex conjugation.

The initial scattering data are the entries of the so-called scattering matrix and the coefficients and parameters of two spectral sums. Denoting by

S(λ) = ( \begin{matrix}T (λ) & L(λ) \\ R(λ) & T (λ) \end{matrix}) ,

the scattering matrix, T (λ) represents the (initial) transmission coefficient, while L(λ) and R(λ) stand for the initial reflection coefficients from the left and from the right, respectively. If T (λ) has no poles in the complex upper half plane ℂ^{+} , there are no spectral sums to identify.

Otherwise, denoting by λ_{1}, …, λ_{n} the so-called bound states that is the finitely many poles of T (λ) in ℂ^{+} and by m_{1}, …, m_{n} the corresponding multiplicities, we have to identify the coefficients {(Γ_{ℓ})_{js} , (Γ_{r})_{js} } as well the parameters {n, m_{j} , λ_{j} } of the initial spectral sums from the left and from the right

(2.2)  S_{ℓ}(α) = ∑_{j=1}^{n} e^{−λ_{j} α} ∑_{s=0}^{m_{j} −1} (Γ_{ℓ})_{js} \frac{α^{s}}{s!} , α ⩾ 0

S_{r}(α) = ∑_{j=1}^{n} e^{λ_{j} α} ∑_{s=0}^{m_{j} −1} (Γ_{r})_{js} \frac{α^{s}}{s!} , α ⩽ 0

where the coefficients (Γ_{ℓ})_{js} and (Γ_{r})_{js} are the so-called norming constants from the left and from the right, respectively, and 0^{0} = 1.

In the IST technique, a crucial role is played by the initial Marchenko kernels from the left Ω_{ℓ}(α) and from the right Ω_{r}(α), which are connected to the above spectral coefficients and spectral sums as follows

Ω_{ℓ}(α) = ρ(α) + S_{ℓ}(α), for α ⩾ 0
Ω_{r}(α) = ℓ(α) + S_{r}(α), for α ⩽ 0

where

ρ(α) = \frac{1}{2π} ∫_{−∞}^{+∞} R(λ)e^{iλα} dλ = F^{−1} \{R(λ)\}

is the inverse Fourier transform of the reflection coefficient from the right R(λ) and

ℓ(α) = \frac{1}{2π} ∫_{−∞}^{+∞} L(λ)e^{−iλα} dλ = \frac{1}{2π} F\{L(λ)\},

apart from the factor 1/2π, is the Fourier transform of the reflection coefficient from the left L(λ).

We note that Ω_{ℓ}(α) and Ω_{r}(α), respectively, reduce to

(a) S_{ℓ}(α) and S_{r}(α) if the reflection coefficients vanish (reflectionless case);

(b) ρ(α) and ℓ(α) if there are no bound states.

3. Auxiliary functions

Let us now introduce, for y ⩾ x, the two pairs of unknown auxiliary functions

K̄(x, y) ≡ ( \begin{matrix}K̄^{up} (x, y) \\ K̄^{dn} (x, y) \end{matrix}) , K(x, y) ≡ ( \begin{matrix}K^{up}(x, y) \\ K^{dn}(x, y) \end{matrix}) ,

and, for y ⩽ x, the two other pairs of unknown auxiliary functions

M̄(x, y) ≡ ( \begin{matrix}M̄^{up} (x, y) \\ M̄^{dn} (x, y) \end{matrix}) , M(x, y) ≡ ( \begin{matrix}M^{up}(x, y) \\ M^{dn}(x, y) \end{matrix}) .

Each of these pair of functions, given the initial solution, is the solution of a system of two structured Volterra integral equations [7, 8, 14].

More precisely, in the focusing case, which is the case we are addressing in this paper, for y ⩾ x, the unknown pair (K̄^{up} , K̄^{dn} ) is the solution of the system

(3.1)

K̄^{up} (x, y) = − ∫_{x}^{∞} u_{0}(z)K̄^{dn} (z, z + y − x) dz,
K̄^{dn} (x, y) = \frac{1}{2} ū_{0}(\frac{1}{2}(x + y)) + ∫_{x}^{(x+y)/2} ū_{0}(z)K̄^{up} (z, x + y − z) dz,

as well as the pair (K^{up}, K^{dn}) it is of the system

(3.2)

K^{up}(x, y) = −\frac{1}{2}u_{0}(\frac{1}{2}(x + y)) − ∫_{x}^{(x+y)/2} u_{0}(z)K^{dn}(z, x + y − z) dz,
K^{dn}(x, y) = ∫_{x}^{∞} ū_{0}(z)K^{up}(z, z + y − x) dz.

Similarly, for y ⩽ x, the unknown pair (M^{up}, M^{dn}) is the solution of the system

(3.3)

M^{up}(x, y) = ∫_{−∞}^{x} u_{0}(z)M^{dn}(z, z + y − x) dz
M^{dn}(x, y) = −\frac{1}{2} ū_{0}(\frac{1}{2}(x + y)) − ∫_{(x+y)/2}^{x} ū_{0}(z)M^{up}(z, x + y − z) dz

as well as the pair (M̄^{up} , M̄^{dn} ) it is of the system

(3.4)

M̄^{up} (x, y) = \frac{1}{2}u_{0}(\frac{1}{2}(x + y)) + ∫_{(x+y)/2}^{x} u_{0}(z)M̄^{dn} (z, x + y − z) dz
M̄^{dn} (x, y) = − ∫_{−∞}^{x} ū_{0}(z)M̄^{up} (z, z + y − x) dz.

From the computational point of view, it is important to note that each auxiliary function is uniquely determined on the bisector y = x, by the initial solution or its partial integral energy.

In fact, setting y = x in each of the four Volterra systems, we immediately obtain

(3.5)  K̄^{dn} (x, x) = \frac{1}{2}ū_{0}(x)  K̄^{up} (x, x) = −\frac{1}{2} ∫_{x}^{∞} |u_{0}(z)|^{2} dz,
(3.6)  K^{up}(x, x) = −\frac{1}{2}u_{0}(x), K^{dn}(x, x) = −\frac{1}{2} ∫_{x}^{∞} |u_{0}(z)|^{2} dz,
(3.7)  M^{dn}(x, x) = −\frac{1}{2}ū_{0}(x), M^{up}(x, x) = −\frac{1}{2} ∫_{−∞}^{x} |u_{0}(z)|^{2} dz,
(3.8)  M̄^{up} (x, x) = \frac{1}{2}u_{0}(x),  M̄^{dn} (x, x) = −\frac{1}{2} ∫_{−∞}^{x} |u_{0}(z)|^{2} dz.

Remark 3.1. If u_{0}(x) is real we need to solve uniquely (3.1) and (3.3) as K^{up}(x, y)= −K̄^{dn} (x, y) and K^{dn}(x, y)= K̄^{up} (x, y) as well as M̄^{up} (x, y)= M^{dn}(x, y) and M̄^{dn} (x, y) = −M^{up}(x, y).

4. Initial Marchenko kernels and scattering matrix

Once the auxiliary functions have been computed, both Ω_{ℓ} and Ω_{r} can be approximated by solving a Volterra integral equation. In fact, using the strong connection between Ω_{ℓ} and the function pairs (K̄^{up} , K̄^{dn} ) and (K^{up}, K^{dn}) as well as that between Ω_{r} and the function pairs (M̄^{up} , M̄^{dn} ) and (M^{up}, M^{dn}), we have [14]

(4.1) K̄^{dn} (x, y) + Ω_{ℓ}(x + y) + ∫_{x}^{∞} K^{dn}(x, z)Ω_{ℓ}(z + y) dz = 0, for y ⩾ x ⩾ 0,

(4.2) M^{dn}(x, y) − Ω_{r}(x + y) − ∫_{−∞}^{x} M̄^{dn} (x, z)Ω_{r}(z + y) dz = 0, for y ⩽ x ⩽ 0.

Given the auxiliary vectors K̄, K, M̄ and M, relations (4.1)–(4.2) can be interpreted as Volterra integral equations having Ω_{ℓ} and Ω_{r} as unknowns.

Remark 4.1. We point out that, from the computational point of view, each initial Marchenko kernel can be treated as a function of only one variable, as we only have to deal with the sum of the two variables.

Following the procedure proposed in [14], the entries of the scattering matrix S(λ) can be computed as follows

(4.3)  T (λ) = \frac{1}{a_{ℓ4}(λ)} = \frac{1}{a_{r1}(λ)} ,
(4.4)  L(λ) = \frac{a_{ℓ2}(λ)}{a_{ℓ4}(λ)} = −\frac{a_{r2}(λ)}{a_{r1}(λ)} ,
(4.5)  R(λ) = \frac{a_{r3}(λ)}{a_{r1}(λ)} = −\frac{a_{ℓ3}(λ)}{a_{ℓ4}(λ)} ,

where

a_{ℓ1}(λ) = 1 − ∫_{ℝ^{+}} e^{−iλz} ( ∫_{ℝ} u_{0}(y)K̄^{dn} (y, y + z) dy) dz,
a_{ℓ2}(λ) = − ∫_{ℝ} e^{2iλy} u_{0}(y) dy − ∫_{ℝ} e^{iλz} ( ∫−∞^{\frac{z}{2}} u_{0}(y)K^{dn}(y, z − y) dy) dz,
a_{ℓ3}(λ) = ∫_{ℝ} e^{−2iλy} ū_{0}(y) dy + ∫_{ℝ} e^{−iλz} ( ∫−∞^{\frac{z}{2}} ū_{0}(y)K̄^{up} (y, z − y) dy) dz,
a_{ℓ4}(λ) = 1 + ∫_{ℝ^{+}} e^{iλz} ( ∫_{ℝ} ū_{0}(y)K^{up}(y, y + z) dy) dz,
a_{r1}(λ) = 1 + ∫_{ℝ^{+}} e^{iλz} ( ∫_{ℝ} u_{0}(y)M^{dn}(y, y − z) dy) dz,
a_{r2}(λ) = ∫_{ℝ} e^{2iλy} u_{0}(y) dy + ∫_{ℝ} e^{iλz} ( ∫ +∞_{\frac{z}{2}} u_{0}(y)M̄^{dn} (y, z − y) dy) dz,

a_{r3}(λ) = − ∫_{ℝ} e^{−2iλy} ū_{0}(y) dy − ∫_{ℝ} e^{−iλz} ( ∫ +∞_{\frac{z}{2}} ū_{0}(y)M^{up}(y, z − y) dy) dz,

a_{r4}(λ) = 1 − ∫_{ℝ^{+}} e^{−iλz} ( ∫_{ℝ} ū_{0}(y)M̄^{up} (y, y − z) dy) dz.

Remark 4.2. If the solution is a soliton or a multisoliton, the only spectral data to be computed are the coefficients and the parameters of the spectral sums S_{ℓ}(α) and S_{r}(α). The coefficients and parameters of S_{ℓ}(α) can be computed by applying the matrix-pencil method recently proposed in [10] to a sufficiently large set of equispaced data of Ω_{ℓ}(α).

5. Computational strategy

In this section we propose a numerical procedure to evaluate the initial scattering data in the case

(5.1)  u_{0}(x) = 0, for |x| > L.

We note that (5.1) can be considered acceptable whenever u_{0}(x) → 0 for |x| → ∞, provided that L be large enough.

Assuming for computational simplicity that the reflection coefficients are zeros, we must solve

1. the systems of Volterra integral equations (3.1)–(3.4);

2. the Volterra integral equations (4.1)–(4.2);

3. a nonlinear approximation problem [10].

5.1. Auxiliary functions computation. Hereafter we assume u_{0} to be real, though the algorithms remain essentially the same in the complex case. In this case (Remark 3.1) we need to solve only systems (3.1) and (3.3), instead of systems (3.1)–(3.4).

Let us first consider system (3.1) to identify the supports of K̄^{up} (x, y) and K̄^{dn} (x, y).

Lemma 5.1. Under hypothesis (5.1), the following properties hold true

(1) For x > L, K̄^{up} (x, y) and K̄^{dn} (x, y) are both zero;

(2) If x ⩽ L and x + y ⩾ 2L then K̄^{up} (x, y) = K̄^{dn} (x, y) = 0;

(3) If x < −L and y − x > 4L then K̄^{up} (x, y) = 0;

(4) If x < −L and x + y < 2L then K̄^{dn} (x, y) = 0.

Proof. To prove (1) it is sufficient to note that, as y ⩾ x and then \frac{1}{2}(x+y) ⩾ 0, u_{0}(\frac{1}{2}(x + y)) = 0 as well as u_{0}(z) = 0 for z ∈ (x, \frac{1}{2}(x + y)). About (2), we consider the sequence x_{i} = ih, i = 0, ±1, ±2, ... where h = \frac{L}{n}, n ∈ ℕ, and then collocate system (3.1) in the node (x_{n}, x_{n+2} ). Noting that both u_{0}(\frac{1}{2}(x + y)) = 0 and u_{0}(z) = 0 for z ⩾ x_{n+1} , we can write:

K̄^{up} (x_{n}, x_{n+2} ) = − ∫_{x_{n}}^{x_{n+1}} u_{0}(z)K̄^{dn} (z, z + 2h) dz,
K̄^{dn} (x_{n}, x_{n+2} ) = ∫_{x_{n}}^{x_{n+1}} u_{0}(z)K̄^{up} (z, x_{2n+2} − z) dz.

(5.2)

Figure 1. Supports of the auxiliary functions K̄^{up} (to the left) and K̄^{dn} (to the right)

Applying then the composite trapezoidal quadrature formula to the computation of the two integrals, and denoting by K̄_{n,n+r}^{up} and K̄_{n,n+r}^{dn} , the approximate values of K̄^{up} (x_{n}, x_{n+r} ) and K̄^{dn} (x_{n}, x_{n+r} ), respectively, we obtain the nonsingular homogeneous system

K̄_{n,n+2}^{up} + \frac{h}{2} u_{0}(x_{n})K̄_{n,n+2}^{dn} = 0

−\frac{h}{2} u_{0}(x_{n})K̄_{n,n+2}^{up} + K̄_{n,n+2}^{dn} = 0 ,

whose solution is K̄_{n,n+2}^{up} = K̄_{n,n+2}^{dn} = 0, for any fixed h value. Applying recursively the same procedure to the nodal points (x_{j} , x_{2n+2ℓ−j} ), with j = n−1, n−2, … and ℓ = 0, 1, 2, … and considering that h is arbitrary, the result follows immediately.

For (3) we note that (2) implies that K̄^{dn} (z, z + y − x) = 0, as z ∈ [−L, L] and y − x > 4L, that is 2z + y − x > 2L.

Result (4) is immediate, as the integration domain appearing in the second equation of (3.1) is null and u_{0}(\frac{1}{2}(x + y)) = 0. □

Taking into account the above properties, we can say that the supports of the auxiliary functions are those represented in Figure 1.

For the numerical solution of system (3.1), the following properties are also important

1. If x ⩽ −L, whatever h, K̄^{up} (x, y) is constant on the line y = x + h. This can be seen taking into account that supp(u_{0}) = [−L, L] and noting that

K̄^{up} (x, x + h) = − ∫_{−L}^{L} u_{0}(z)K̄^{dn} (z, z + h) dz

only depends on h. For this reason we put K̄^{up} (x, x + h) = C_{K̄ ,h}^{up} for each given value h.

Figure 2. Additional properties of K̄^{up} (to the left) and K̄^{dn} (to the right)

2. If x < −L and x + y > −2L, K̄^{dn} (x, y) is constant on each line x + y = −2(L + h) for each 0 < h < 2L. In fact, by the second equation of (3.1) it results that

K̄^{dn} (x, y) = \frac{1}{2}u_{0}(−L + h) + ∫_{−L}^{L+h} u_{0}(z)K̄^{up} (z, −2(L + h) − z) dz = C_{K̄ ,h}^{dn}.

These two results are graphically represented in Figure 2.

A visualization of the area where we need to compute K̄^{up} and K̄^{dn} is given by the gray triangle represented in Figure 3. In the remaining areas of the respective supports their values are immediately obtained by using those of the gray triangle. The gray line shows, in particular, the values of the gray triangle we use to know K̄^{up} and K̄^{dn} in the gray point of the gray area.

Fixed n ∈ ℕ, let x_{i} = ih, i = −n, −n + 1, …, n and y_{j} = jh with h = L/n, j = −n, −n + 1, …, 3n. The algorithm first requires the computation of K̄^{up} and K̄^{dn} in the nodal points of the bisector (x_{i}, y_{j}) of the gray triangle. Recalling (3.5) and denoting by K̄_{r,s}^{up} , K̄_{r,s}^{dn} and u_{0,r} the approximation of K̄^{up} (x, y), K̄^{dn} (x, y) and u_{0}(x) in the nodal points (x_{r}, y_{s}), we can write:

K̄_{n−i+1,n−i+1}^{dn}  = \frac{1}{2}u_{0,n−i+1} ,
K̄_{n−i+1,n−i+1}^{up}  = −\frac{1}{2} ∫_{x_{n−i+1}}^{L} |u_{0}(z)|^{2} dz, i = 1, 2, …, 2n + 1.

The approximation of K̄_{n−i+1,n−i+1}^{up} is then obtained recursively by applying the composite Simpson’s quadrature formula. Fixed j = 1, 2, …, 2n, system (3.1) is collocated in the nodal points

D_{j} = {(x_{n−j−i+1} , x_{n+j−i+1} ) ∈ ℝ^{2} | i = 1, …, 2n + 1}.

We note that in each collocation point we have to compute K̄^{up} and K̄^{dn} , given their values in the bisector y = x. For this reason, recalling that K̄^{up} (x_{i}, y_{j}) and

Figure 3. Geometrical visualization of the computational area

K̄^{dn} (x_{i}, y_{j}) are zero for x_{i} + y_{j} > 2L, we collocate system (3.1) following the order depicted in Figure 4.

Our numerical algorithm is based on the approximation of integrals (3.1) by means of the composite trapezoidal formula. In this way, following the ordering depicted in Figure 4, for each collocation point we have only to solve recursively a sequence of nonsingular 2 × 2 linear systems. In fact, for each collocation point, we must solve the following system:

K̄_{n−j−i+1,n+j−i+1}^{up} + \frac{h}{2}u_{0,n−j−i+1} K̄_{n−j−i+1,n+j−i+1}^{dn}

= −h(1 − δ_{i−1,0} ) ∑_{ℓ=1}^{i−1} u_{0,n−j−i+1+ℓ} K̄_{n−j−i+ℓ+1,n+j−i+ℓ+1}^{dn}

− \frac{h}{2}u_{0,n−j−i+1} K̄_{n−j−i+1,n+j−i+1}^{up} + K̄_{n−j−i+1,n+j−i+1}^{dn}

= \frac{u_{0,n−i+1}}{2} + \frac{h}{2} ∑_{ℓ=1}^{j} (2 − δ_{ℓj})u_{0,n−j−i+ℓ+1} K̄_{n−j−i+ℓ+1,n+j−i−ℓ+1}^{up}

where δ_{ij} is the Kronecker symbol.

Let us now consider system (3.3) and identify the supports of the auxiliary functions M^{up} and M^{dn}. Taking into account hypothesis (5.1), they are characterized by the following

Lemma 5.2. Under hypothesis (5.1), the following statements hold true:

(1) For x < −L both M^{up}(x, y) and M^{dn}(x, y) are zero;

(2) For x ⩾ −L and x + y ⩽ −2L, M^{up}(x, y) = M^{dn}(x, y) = 0;

(3) If x > L and y − x < −4L then M^{up}(x, y) = 0;

(4) If x > L and y + x > 2L, then M^{dn}(x, y) = 0.

Figure 4. Sorting visualization of collocation points in the triangle

Proof. (1) It is enough to note that \frac{1}{2}(x + y) < −L and then M^{up}(x, y) = M^{dn}(x, y) = 0.

In order to prove (2), we observe that collocating system (3.3) in (x_{−n} , x_{−n−2} ) and taking into account that u_{0}(\frac{1}{2}(x + y)) = 0 and u_{0}(z) = 0 for z < x_{−n} , one can write

M_{−n,−n−2}^{up} − \frac{h}{2} u_{0}(x_{n})M_{−n,−n−2}^{dn} = 0
\frac{h}{2} u_{0}(x_{−n} )M_{−n,−n−2}^{up} + M_{−n,−n−2}^{dn} = 0 ,

whose solution is M_{−n,−n−2}^{up} = M_{−n,−n−2}^{dn} for any given h value. Iterating the procedure in the nodal points (x_{−j} , x_{−j−2ℓ} ) with j = n−1, n−2, … and ℓ = 1, 2, … and considering that h is arbitrarily, we obtain that M_{−j,−j−2ℓ}^{up} = M_{−j,−j−2ℓ}^{dn} = 0, for the above fixed value of j and ℓ.

Concerning (3), it is sufficient to note that, as the domain of integration is [−L, L] and 2z + y − x < −2L, M^{up}(x, y) = 0 as a consequence of (2).

Finally (4) is immediate considering that M^{dn}(x, y) = 0 for x + y > 2L as well as u_{0}(\frac{1}{2}(x + y)). □

A geometrical representation of these supports is given in Figure 5.

Moreover, in order to find a numerical solution of (3.3), it is important to note that the unknown functions M^{up} and M^{dn} have also the following properties, represented in Figure 6:

1. If x ⩾ L, whatever h, M^{up} is constant on the line y = x+h. This property can immediately be verified, noting that

M^{up}(x, y) = ∫_{−L}^{L} u_{0}(t)M^{dn}(z, z + h) dz

only depends on h. Hence, to make evident this property, for x ⩾ L, we write M^{up}(x, y) = C_{M,h}^{up} .

Figure 5. Supports of the auxiliary functions M^{up} (to the left) and M^{dn} (to the right)

Figure 6. Additional properties of M^{up} (to the left) and M^{dn} (to the right)

2. If x ⩾ L and x + y < 2L, M^{dn} is constant on each line x + y = 2(L − h), for 0 < h < 2L. In fact using the second equation of (3.3), we have

M^{dn}(x, y) = −\frac{1}{2}u_{0}(L − h) − ∫_{L−h}^{L} u_{0}(z)M^{up}(z, 2(L − h) − z) dz = C_{M,h}^{dn} .

As a result, as for system (3.1), we need only to compute the unknowns in the gray triangle depicted in Figure 7, since in the remainder of the support we can apply the properties discussed above.

The computational strategy developed for the numerical solution of this system is essentially the same adopted for system (3.1). Hence, at first by using (3.7) we compute

M_{n−i+1,n−i+1}^{up} = −\frac{1}{2} ∫_{−∞}^{x_{n−i+1}} |u_{0}(z)|^{2}dz, M_{n−i+1,n−i+1}^{dn} = −\frac{1}{2}u_{0,n−i+1} .

Figure 7. Geometrical visualization of the computational area

After that, fixed j = 1, 2, …, 2L, system (3.3) is collocated in the nodal points

D_{j} = {(x_{−n+j+i−1} , x_{−n−j+i−1} ) ∈ ℝ^{2} | i = 1, …, 2n − j + 1},

and the integrals are approximated by using the composite trapezoidal rule. Operating in this way, we obtain the sequence of 2 × 2 systems

M_{i+j−n−1,i−j−n−1}^{up} − \frac{h}{2}u_{0,i+j−n−1} M_{i+j−n−1,i−j−n−1}^{dn}

= h(1 − δ_{i−1,0} ) ∑_{ℓ=1}^{i−1} u_{0,i+j−n−1−ℓ} M_{i+j−n−ℓ−1,i−j−n−ℓ−1}^{dn}

\frac{h}{2}u_{0,i+j−n−1} M_{i+j−n−1,i−j−n−1}^{up} + M_{i+j−n−1,i−j−n−1}^{dn}

= −\frac{u_{0,i−n−1}}{2} − \frac{h}{2} ∑_{ℓ=1}^{j} (2 − δ_{ℓj})u_{0,i+j−n−ℓ−1} M_{i+j−n−ℓ−1,i+ℓ−j−n−1}^{up}

that we solve recursively, by following the ordering depicted in Figure 8.

5.2. Numerical solution of Marchenko equations. Once K̄^{up} and K̄^{dn} have been computed, we have to solve the integral equations (4.1) and (4.2) that, by Remark 3.1, can be written as

K̄^{dn} (x, y) + Ω_{ℓ}(x + y) + ∫_{x}^{∞} K̄^{up} (x, z)Ω_{ℓ}(z + y) dz = 0 (5.3)

(5.4)  M^{dn}(x, y) − Ω_{r}(x + y) + ∫_{−∞}^{x} M^{up}(x, z)Ω_{r}(z + y) dz = 0 ,

for y ⩾ x ⩾ 0 and y ⩽ x ⩽ 0, respectively. Let us first consider (5.3). Fixed a steplenght h and taken a set of nodal points (x_{i}, x_{j}) with x_{i} = ih, i = 0, 1, … and

Figure 8. Sorting visualization of collocation points in the triangle

j = i, i + 1, we collocate (5.3) in (x_{i}, x_{j}), taking a great advantage by the support of K̄^{up} and K̄^{dn} as explained in the following

Lemma 5.3. If u_{0}(x) = 0 for |x| > L, then Ω_{ℓ}(z) = 0 for z > 2L.

Proof. Given the steplenght h and assuming x_{n} = nh = L, we set K̄^{dn} (x_{i}, x_{j}) = K̄_{ij}^{dn} and K̄^{up} (x_{i}, x_{j}) = K̄_{ij}^{up} as well as Ω_{ℓ}(x_{i} + y_{j}) = (Ω_{ℓ})_{i+j} . Recalling that, as proved in Lemma 5.1, K̄_{n,n+m}^{dn} = K̄_{n,n+m}^{up} = 0, m = 1, 2, … and using (5.3) it is immediate to state that (Ω_{ℓ})_{2n+m} = 0, m = 1, 2, … Moreover, considering that K̄^{up} (x, z) = 0 for x ⩾ x_{n} and recalling (3.5), equation (5.3) implies that

(5.5)  (Ω_{ℓ})_{2n} = (Ω_{ℓ})(x_{n}, x_{n}) = −K̄_{n,n}^{dn} = −\frac{1}{2} ū_{0}(x_{n}).

As a result, if the support of u_{0} is [−L, L] then the support of Ω_{ℓ}(α) is [0, 2L]. We can similarly argue that supp(Ω_{r}(z)) ⊂ [−2L, 0] if supp(u_{0}) ⊂ [−L, L]. □

To compute Ω_{ℓ} in [0, 2L], first we collocate (5.3) in (x_{n−1} , x_{n}) and approximate the integral by the trapezoidal rule. Proceeding in this way we obtain that (Ω_{ℓ})_{2n−1} is the solution of the equation

(1 + \frac{h}{2} K̄_{n−1,n−1}^{up} )(Ω_{ℓ})_{2n−1} = −K̄_{n−1,n}^{dn} + \frac{h}{4} K̄_{n−1,n}^{up} ū_{0}(x_{n}).

Collocating now in (x_{n−1} , x_{n−1} ) and adopting the same procedure, we obtain that for h small enough (Ω_{ℓ})_{2n−2} is the solution of the equation

(1 + \frac{h}{2} K̄_{n−1,n−1}^{up} )(Ω_{ℓ})_{2n−2} = −K̄_{n−1,n−1}^{dn} − \frac{h}{2} K̄_{n−1,n}^{up} (Ω_{ℓ})_{2n−1} .

Iterating the procedure we collocate (5.3) in (x_{n−j} , x_{n−j+k} ) with j = 1, 2, …, n and k = 0, 1 (if j = n, then k = 0). So we obtain that (Ω_{ℓ})_{2(n−j)+k} is a solution of

(5.6) ( 1 + \frac{h}{2} K̄_{n−j,n−j}^{up} ) (Ω_{ℓ})_{2(n−j)+k}

= −K̄_{n−j,n−j+k}^{dn} − \frac{h}{2} ∑_{ℓ=1}^{j−1} K̄_{n−j,n−j+ℓ}^{up} (Ω_{ℓ})_{2(n−j)+k+ℓ} .

As |K̄^{dn} (x, x)| is decreasing, h can be considered sufficiently small if 1− \frac{h}{2} K̄_{0,0}^{up} > 0.

The same approach allows us to compute Ω_{r}, that is to solve (5.4) in a sequence of discretization points of its support x_{i} = ih, i = 0, −1, …, −2n, x_{−n} = −L. To this end, recalling that M^{up}(x, z) = 0 if x + z < −2L and using (5.4), first we note that in the bisector y = x

(Ω_{r})_{−2n} = Ω_{r}(x_{−n} , x_{−n} ) = M_{−n,−n}^{dn} = −\frac{1}{2} ū_{0}(x_{−n} ).

Collocating (5.4) in (x_{−n+1} , x_{−n} ) and approximating the integral value by the trapezoidal rule we obtain the equation

(1 − \frac{h}{2}M_{−n+1,−n+1}^{up} )(Ω_{r})_{−2n+1} = M_{−n+1,−n}^{dn} + \frac{h}{4}M_{−n+1,−n}^{up} u_{0}(x_{−n} ).

Iterating the procedure, we collocate (5.4) in (x_{−n+j} , x_{−n+j−k} ) with j = 1, 2, …, n and k = 0, 1 (if j = n, then k = 0), we obtain that (Ω_{r})_{2(j−n)−k} is a solution of

(1 − \frac{h}{2}M_{−n+j,−n+j}^{up} )(Ω_{r})_{2(j−n)−k}

= M_{−n+j,−n+j−k}^{dn} + \frac{h}{2} ∑_{ℓ=1}^{j−1} M_{j−n,j−n−ℓ}^{up} (Ω_{r})_{2(j−n)−k−ℓ} ,

for any h small enough, that is for 1 + \frac{h}{2}M_{0,0}^{up} > 0.

5.3. Bound states and norming constants. Let us give a brief description of the matrix-pencil method that we have recently developed for the identification of the bound states and the norming constants [10]. Setting z_{j} = e^{−λ_{j}} , the spectral function sum S_{ℓ}(α) introduced in (2.2), can be represented as the monomial-power sum

S_{ℓ}(α) = ∑_{j=1}^{n} ∑_{s=0}^{m_{j} −1} c_{js} α^{s}z_{j}^{α} , 0^{0} ≡ 1.

Setting M = m_{1} + ⋯+ m_{n}, the method allows to compute the parameters {n, m_{j} , z_{j} } and the coefficients {c_{js} } given S_{ℓ}(α) in 2N integer values (N > M)

α = α_{0}, α_{0} + 1, …, α_{0} + 2N − 1, with α_{0} ∈ ℕ^{+} = {0, 1, 2, …}.

The basic idea of the method is the interpretation of S_{ℓ}(α) as the general solution of a homogeneous linear difference equation ∑ _{k=0}^{M} p_{k}S_{k+α_{0}} = 0 of order M whose characteristic polynomial (Prony polynomial)

P (z) = ∏_{j=1}^{n} (z − z_{j})^{m_{j}} = ∑_{k=0}^{M} p_{k}z^{k}, p_{M} ≡ 1

is uniquely characterized by the z_{j} values we are looking for. The identification of the zeros {z_{j} } allows to compute the coefficients c_{js} by solving a linear system.

For the computation of {z_{j} } and then of the bound states λ_{j}, the given data are arranged in the two Hankel matrices of order N

(S_{ℓ}^{0})_{ij} = S_{ℓ}(i + j − 2), (S_{ℓ}^{1})_{ij} = S_{ℓ}(i + j − 1), i, j = 1, 2, …, N.

To these matrices we associate the M × M matrix-pencil

S_{MM} (z) = (S_{NM}^{0} )^{∗}(S_{NM}^{1} − zS_{NM}^{0} )

where the asterisk denotes the conjugate transpose. As proved in [10], the zeros z_{j} of the Prony polynomial, with their multiplicities, are exactly the generalized eigenvalues of the matrix-pencil S_{MM} (z). Then, by applying the generalized singular value decomposition to the matrices S_{NM}^{0} and S_{NM}^{1} , the algorithm developed allows to compute the zeros z_{j} and then the bound states λ_{j}, as λ_{j} = − logz_{j} .

The vector of coefficients c = [c_{1 0}, …, c_{1 m_{1} −1} , …, c_{n 0} , …, c_{n m_{n} −1} ]^{T} is then computed by solving (in the least square sense) the overdetermined linear system K_{NM}^{0} c = S_{ℓ}^{0} where S_{ℓ}^{0} = [S_{ℓ}(0), S_{ℓ}(1), …, S_{ℓ}(N − 1)]^{T} and K_{NM}^{0} is the Casorati matrix associated to the monomial powers {k^{s}z_{j}^{k}} for k = 1, …, N − 1.

If m_{j} ≡ 1, the Casorati matrix K_{Nn}^{0} reduces to the Vandermonde matrix (V)_{ij} = z_{j}^{i+1} of order N × n associated to the zeros z_{1}, …, z_{n}. The solution of the Casorati system allows us to immediately compute the norming constants as (Γ_{ℓ})_{js} = s!c_{js} .

The coefficients {(Γ_{r})_{js} } are then obtained by solving, in the least square sense, a linear system whose vector of known data is given by Ω_{r}(α) evaluated in a set of N points, with a sufficiently large N > M.

6. Numerical results

In order to access the effectiveness of our method in the approximation of Ω_{ℓ}, we adopted the following estimator for the relative error

E_{n} = \frac{∥Ω̃_{ℓ,n} − Ω_{ℓ}∥_{∞}}{∥Ω_{ℓ}∥_{∞}}

where Ω̃_{ℓ,n} denotes the array of the approximated Marchenko kernel computed by means of equations (5.6) in 2n equispaced points belonging to the interval [0, 2L], Ω_{ℓ} is the array of the exact Marchenko kernel evaluated at the same points and ∥ · ∥_{∞} is the usual infinity norm. Let us now show our results in two cases in which Ω_{ℓ} is analytically known.

Test 1. Let us consider as initial potential the soliton given by

(6.1)  u_{0}(x) = \frac{−2ce^{−2ax}}{1 + \frac{|c|^{2}}{4p^{2}} e^{−4px}} , x ∈ ℝ,

where c, a and p are real parameters [15]. In this case the Marchenko kernel on the left is Ω_{ℓ}(x) = e^{−ax} . Considering that u_{0}(x) ⩽ 10^{−13} if |x| > 15, in our computation we assumed L = 15. Hence taken n and a steplenght h such that nh = 2L, we computed Ω_{ℓ}(α) in the n points α_{i} = ih, i = 0, 1, …, n and we reported in Table 1 the relative error E_{n} for different values of n.

Test 2. Let us now take as initial potential a multisoliton represented by four solitons which interact each other nonlinearly [14], namely

(6.2)  u_{0}(x) = \{ \begin{matrix}−2b^{∗}[e^{2xA^{∗}} + Qe^{−2xA} N]^{−1} c^{∗} & x ⩾ 0 \\ −2c [e^{−2xA} + Ne^{2xA^{∗}} Q]^{−1} b & x < 0 \end{matrix}

where b and c are column and row vectors, respectively, A is a matrix with eigenvalues α_{j} having positive real parts and N and Q are two matrices obtained by Table 1. Relative errors E_{n} for Ω_{ℓ} in the soliton case (left) and in the multisoliton case (right)

\begin{matrix}n & E_{n} \\ 300 & 1.03e-03 \\ 600 & 2.62e-04 \\ 900 & 1.17e-04 \\ 1200 & 6.60e-05 \end{matrix}  \begin{matrix}n & E_{n} \\ 300 & 8.42e-03 \\ 600 & 2.08e-03 \\ 900 & 9.27e-04 \\ 1200 & 5.21e-04 \end{matrix}

solving the respective Lyapunov equations:

QA + A^{∗}Q = c^{∗}c, AN + NA^{∗} = bb^{∗}.

In this case, it is possible to prove [14] that the initial Marchenko kernel Ω_{ℓ} is given by Ω_{ℓ}(α) = ce^{−αA} b, the reflection coefficients R(λ) = L(λ) = 0, the bound state terms are λ_{j} = ia_{j} and the norming constants Γ_{ℓ,j} = b_{j} c_{j}.

As in the previous example we assumed L = 15 as u_{0}(x) ⩽ 10^{−13} for |x| > 15. Moreover, for simplicity we considered A = diag([1, 2, 3, 4]), b = [1, 2, −2, −1]^{T}, c = [2, 1, 1, 2], which implies that Ω_{ℓ}(α) = ∑ _{j=1}^{4} b_{j} c_{j} e^{−αa_{j}} . The results reported in Table 1 show that, as expected, the relative error E_{n} decreases with respect to n.

7. Conclusions and perspectives

The effectiveness of the numerical solution of the direct scattering problem in the NLS, as probably in various other NPDEs of integrable type, basically depends on the effectiveness of the numerical solution of Volterra’s systems of integral equations with structured kernels on unbounded domains and then on the identification of parameters in monomial-exponential sums. Our experiments show that in the reflectionless case our matrix-pencil method for the identification of spectral parameters is fully reliable whenever the relative error coming from the solution of systems of Volterra is small enough.

A challenging mathematical problem is to develop effective algorithms to approximate quite well the reflection coefficients, that is to compute ratios (4.4)–(4.5) or to generate alternative formula for their evaluation. This challenging task is devoted to another paper, as well the generation of a new family of numerical methods for computing more efficiently the auxiliary functions and the Marchenko kernels.

References

1. M. J. Ablowitz, P. A. Clarkson, Solitons, Nonlinear Evolution Equations and Inverse Scattering, Cambridge University Press, Cambridge, 1991.

2. M. J. Ablowitz, B. Prinari, A. D. Trubatch, Discrete and Continuous Nonlinear Schrödinger Systems, Cambridge University Press, Cambridge, 2004.

3. M. J. Ablowitz, H. Segur, Solitons and the Inverse Scattering Transform, SIAM, Philadelphia, 1981.

4. G. P. Agrawal, Nonlinear Fiber Optics, Academic Press, New York, 2001.

5. A. Aricò, G. Rodriguez, S. Seatzu, Numerical solution of the nonlinear Schrödinger equation, starting from the scattering data, Calcolo 48(1) (2011), 75–88.

6. A. Aricò, C. van der Mee, S. Seatzu, Structured matrix numerical solution of the nonlinear Schrödinger equation by the inverse scattering transform, Electron. J. Differ. Equ. 15 (2009), 1–21.

7. F. Demontis, C. Van der Mee, Explicit solutions of the cubic matrix nonlinear Schrödinger equation, Inverse Probl. 24 (2008), 02520, 16 pp.

8. , Marchenko equations and norming constants of the matrix Zakharov–Shabat system, Oper. Matrices 2 (2008), 79–113.

9. L. D. Faddeev, L. A. Takhtajan, Hamiltonian Methods in the Theory of Solitons, Classics in Mathematics, Springer, New York, 1987.

10. L. Fermo, C. Van der Mee, S. Seatzu, Parameter estimation of monomial-exponential sums, ETNA, Electron. Trans. Numer. Anal. 41 (2014), 249–261.

11. S. V. Manakov, On the theory of two-dimensional stationary self-focusing of electromagnetic waves, Sov. Phys. JETP 38 (1974), 248–253; translated from Zh. Eksp. Teor. Fiz. 65 (1973), 505–516 [Russian].

12. S. P. Novikov, S. V. Manakov, L. B. Pitaevskii, V. E. Zakharov, Theory of Solitons. The Inverse Scattering Method, Plenum Press, New York, 1984.

13. J. K. Shaw, Mathematical principles of optical fiber communications, CBMS-NSF Regional Conference Series 76 SIAM Philadelphia, 2004.

14. C. Van der Mee, Nonlinear Evolution Models of Integrable Type, SIMAI e-Lect. Notes 11, 2013.

15. V. E. Zakharov, A. B. Shabat, Exact theory of two-dimensional self-focusing and one dimensional self-modulation of waves in nonlinear media, Sov. Phys. JETP 34 (1972), 62–69; translated from Zh. Eksp. Teor. Fiz. 61 (1971), 118–134 [Russian].

Department of Mathematics and Computer Science University of Cagliari Cagliari Italy

 fermo@unica.it
 cornelis@krein.unica.it
 seatzu@unica.it