Recent Advances in the Numerical Solution of the Nonlinear Schrödinger Equation

1. Introduction

In this paper we review the main facts about a recent approach for efficiently solving the following general form of (non-dimensional) Nonlinear Schrödinger Equation (NLSE),

$$\mathrm{i}{\psi }_t(x,t)+{\psi }_{xx}(x,t)+f^{\prime }{\left(\right|\psi (x,t)|}^2)\psi (x,t)=0,(x,t)\in [a,b]\times [0,T],$$

(1)

coupled with initial data and periodic boundary conditions, with $\mathrm{i}$ the imaginary unit, $f:\mathbb{R}\to \mathbb{R}$ a differentiable function, and $f^{\prime }$ its derivative.

Generally speaking, there are basically two ways in which NLSE enters physical and applied sciences. The first one is the appearance of the NLSE as the first-order approximation of a nonlinear dispersive system, such as the wave equation in a nonlinear medium. To this category belong the NLSE used in Fiber Optics [1], Plasma Physics [2,3,4], Geophysics [5,6] and Biology [7]. The other way NLSE appears in physics is as the single-particle, mean-field approximation of a many-particle linear system. This is typically the case of quantum mechanics, where the NLSE arises in several contexts, the most important being the description of the dynamics of a Bose-Einstein condensate. In this case the suitable NLSE is also known as Gross-Pitaevskii equation [8]. In most such applications, the NLSE is used in its more common form, i.e., with cubic nonlinearity, corresponding to $f\left(x\right)=x^2$. However, higher-degree nonlinearities appear in important applications, among which it is worth mentioning the cubic-quintic NLSE, corresponding to $f\left(x\right)=\alpha x^2-\beta x^3$, which is used to describe the propagation of ultrashort pulses in optical fibers [1,9,10,11] or high-density effects in Bose-Einstein condensates [12,13,14]. The NLSE is of course of great interest also from a purely mathematical viewpoint, especially in view of its profound connections with the theory of integrable PDEs, solitons and inverse scattering [15,16,17,18,19].

For all these reasons, the numerical solution of the NLSE has been the subject of investigation since a long time (see, e.g., [20,21,22,23,24] for early approaches). More recently, geometric numerical methods have been considered: as an example, we mention multi-symplectic methods (see, e.g., [25,26,27,28,29,30,31,32,33,34,35,36]), symplectic methods (see, e.g., [37,38,39,40,41,42,43,44]), and so forth (see, e.g., [45,46,47,48,49,50,51]). More recently, methods able to conserve the energy or other invariants (see, e.g., [52,53,54,55,56,57,58,59,60,61]) have been investigated: the methods we shall deal with, are exactly placed in this latter setting. The basic approach follows what suggested in [62, page 187], namely by suitably using the method of lines: if the PDEs are of Hamiltonian type, [...] the space discretization should be carried out in such a way that the resulting system of ODEs is Hamiltonian (for a suitable Poisson bracket). The latter system is then solved by using a geometric method for the time integration. Usually, symplectic methods [62,63,64,65,66] have been used for this purpose. Instead, in the present case, we shall consider energy-conserving methods, namely methods able to conserve the energy. Energy-conserving methods have been studied by many authors (see, e.g., [67,68,69,70]): we shall here consider the class of Hamiltonian Boundary Value Methods (HBVMs), which is a family of Runge-Kutta (RK) methods recently devised for the efficient numerical solution of Hamiltonian ODEs. We refer to [68,71] (and, in particular, to the monograph [67] 1) for the derivation and analysis of such methods, which have been devised within the framework of the so called line integral methods (see, e.g., the review paper [72]). We mention that generalizations and extensions of such approach have been also considered in [73,74,75,76,77,78,79,80,81,82,83,84,85]. HBVMs have been also considered for the efficient numerical solution of a number of Hamiltonian PDEs (see, e.g., [86,87,88,89,90]), among which the NLSE [91] and Manakov systems [92].

With this premise, the structure of the paper is as follows: in Section 2 we recall the basic facts about the NLSE along with its space semi-discretization; in Section 3 we review HBVMs and their properties; in Section 4 we give some details about the efficient implementation of such methods for numerically solving the NLSE (1); Section 5 contains a few numerical tests; at last some concluding remarks are given in Section 6.

2. Basic Facts about the NLSE and Its Space Semi-Discretization

We shall hereafter consider the real form of (1) by setting

$$\mathsf{\Psi }(x,t)=u(x,t)+\mathrm{i}v(x,t).$$

(2)

Omitting, for the sake of brevity, the arguments $(x,t)$ when unnecessary, one then obtains:

$$u_t=-v_{xx}-f^{\prime }(u^2+v^2)v,v_t=u_{xx}+f^{\prime }(u^2+v^2)u,(x,t)\in [a,b]\times [0,T],$$

(3)

with the initial conditions

$$u(x,0)=u_0\left(x\right),v(x,0)=v_0\left(x\right),x\in [a,b],$$

(4)

and periodic boundary conditions (consequently, also $u_0\left(x\right)$ and $v_0\left(x\right)$ are assumed to be periodic). By defining the Hamiltonian functional

$$\mathcal{H}[u,v]\left(t\right)=\frac{1}{2}{\int }_a^b\left[u_x{(x,t)}^2+v_x{(x,t)}^2-f(u{(x,t)}^2+v{(x,t)}^2)\right]\mathrm{d}x,$$

(5)

and considering that for any variations $\eta ,\nu \in L^2[a,b]$, its functional derivatives are given by:

$$\begin{array}{c}\hfill {\displaystyle \underset{\epsilon \to 0}{lim}\frac{\mathcal{H}[u+\epsilon \eta ,v]\left(t\right)-\mathcal{H}[u,v]\left(t\right)}{\epsilon }=:{\int }_a^b\frac{\delta }{\delta u\left(x\right)}\mathcal{H}[u,v]\left(t\right)\eta \left(x\right)\mathrm{d}x}\\ \equiv -{\int }_a^b\left[u_{xx}(x,t)+f^{\prime }(u{(x,t)}^2+v{(x,t)}^2)u(x,t)\right]\eta \left(x\right)\mathrm{d}x,\\ \hfill {\displaystyle \underset{\epsilon \to 0}{lim}\frac{\mathcal{H}[u,v+\epsilon \nu ]\left(t\right)-\mathcal{H}[u,v]\left(t\right)}{\epsilon }=:{\int }_a^b\frac{\delta }{\delta v\left(x\right)}\mathcal{H}[u,v]\left(t\right)\nu \left(x\right)\mathrm{d}x}\\ \equiv -{\int }_a^b\left[v_{xx}(x,t)+f^{\prime }(u{(x,t)}^2+v{(x,t)}^2)v(x,t)\right]\nu \left(x\right)\mathrm{d}x,\end{array}$$

one deduces that (3) has a Hamiltonian structure. In more detail, by defining

$$y=\left(\begin{array}{c}u\\ v\end{array}\right),J=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)=-J^{\top }=-J^{-1},\nabla \mathcal{H}\left[y\right]=\left(\begin{array}{c}\frac{\delta \mathcal{H}}{\delta u}\\ \frac{\delta \mathcal{H}}{\delta v}\end{array}\right),$$

(6)

it is straightforward deriving that (3) can be rewritten as

$$y_t=J\nabla \mathcal{H}\left[y\right],$$

(7)

which is formally similar to a usual system of Hamiltonian ODEs.

Because of the periodic boundary conditions, it can be proved that $\mathcal{H}[u,v]\left(t\right)$ is a conserved quantity for all $t\ge 0$, as well as the mass

$${\mathcal{M}}_1[u,v]\left(t\right)={\int }_a^b\left[u{(x,t)}^2+v{(x,t)}^2\right]\mathrm{d}x,$$

(8)

and the momentum

$${\mathcal{M}}_2[u,v]\left(t\right)={\int }_a^b\left[u(x,t)v_x(x,t)-v(x,t)u_x(x,t)\right]\mathrm{d}x.$$

(9)

Let us sketch the conservation of (8): the conservation of (5) and (9) can be similarly obtained (see, e.g., [91]):

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{M}}_1[u,v]\left(t\right)& =& 2{\int }_a^b\left[u(x,t)u_t(x,t)+v(x,t)v_t(x,t)\right]\mathrm{d}x\hfill \\ & =& 2{\int }_a^b\left\{u(x,t)\left[-v_{xx}(x,t)-f^{\prime }(u{(x,t)}^2+v{(x,t)}^2)v(x,t)\right]\right.\hfill \\ & & \left.+v(x,t)\left[u_{xx}(x,t)+f^{\prime }(u{(x,t)}^2+v{(x,t)}^2)u(x,t)\right]\right\}\mathrm{d}x\hfill \\ & =& 2{\int }_a^b\frac{\mathrm{d}}{\mathrm{d}x}\left[v(x,t)u_x(x,t)-u(x,t)v_x(x,t)\right]\mathrm{d}x=0,\hfill \end{array}$$

because of the periodic boundary conditions.

2.1. Space semi-discretization

Since $u(x,t)$ and $v(x,t)$ in (2) are assumed to be periodic functions of $x\in [a,b]$, we consider their expansion in space along the following orthonormal basis on $[a,b]$:2

$$\begin{array}{ccc}\hfill {\omega }_{2j}\left(x\right)& =& \sqrt{\frac{2-{\delta }_{j0}}{b-a}}cos\left(j\frac{x-a}{b-a}2\pi \right),\hfill \\ \hfill {\omega }_{2j+1}\left(x\right)& =& \sqrt{\frac{2}{b-a}}sin\left((j+1)\frac{x-a}{b-a}2\pi \right),j=0,1,2,\dots ,\hfill \end{array}$$

so that

$${\int }_a^b{\omega }_i\left(x\right){\omega }_j\left(x\right)={\delta }_{ij},i,j=0,1,2,\dots ,$$

with ${\delta }_{ij}$ denoting, as usual, the Kronecker delta. One then obtains that, for suitable time-dependent coefficients $q_j\left(t\right),p_j\left(t\right)$, $j\ge 0$,

$$u(x,t)=\sum _{j\ge 0}{q}_j\left(t\right){\omega }_j\left(x\right),v(x,t)=\sum _{j\ge 0}{p}_j\left(t\right){\omega }_j\left(x\right),$$

(10)

which can be more compactly rewritten, by introducing the infinite vectors

$$\mathsf{\omega }\left(x\right)=\left(\begin{array}{c}{\omega }_0\left(x\right)\\ {\omega }_1\left(x\right)\\ {\omega }_2\left(x\right)\\ ⋮\end{array}\right),\mathit{q}\left(t\right)=\left(\begin{array}{c}{q}_0\left(t\right)\\ q_1\left(t\right)\\ q_2\left(t\right)\\ ⋮\end{array}\right),\mathit{p}\left(t\right)=\left(\begin{array}{c}{p}_0\left(t\right)\\ p_1\left(t\right)\\ p_2\left(t\right)\\ ⋮\end{array}\right),$$

(11)

as

$$u(x,t)=\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\left(t\right),v(x,t)=\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\left(t\right).$$

(12)

Further, by considering that

$${\mathsf{\omega }}^{\prime }\left(x\right)=D_1\mathsf{\omega }\left(x\right),{\mathsf{\omega }}^{\prime \prime }\left(x\right)=D_1^2\mathsf{\omega }\left(x\right)\equiv -D\mathsf{\omega }\left(x\right),$$

(13)

where, by setting hereafter $I_r$ the identity matrix of dimension r, and with reference to the matrix J defined in (6),

$$\begin{array}{ccc}\hfill D_1& =& \frac{2\pi }{b-a}\left(\begin{array}{c}0\\ & 1·J\\ & & 2·J\\ & & & \ddots \end{array}\right),D=\frac{4{\pi }^2}{{(b-a)}^2}\left(\begin{array}{c}0\\ & 1^2·I_2\\ & & 2^2·I_2\\ & & & \ddots \end{array}\right),\hfill \\ \hfill D_1& =& -D_1^{\top },D=D^{\top },\hfill \end{array}$$

(14)

one deduces that it is possible to rewrite (3) as:

$$\begin{array}{ccc}\hfill \mathsf{\omega }{\left(x\right)}^{\top }\dot{\mathit{q}}\left(t\right)& =& \mathsf{\omega }{\left(x\right)}^{\top }D\mathit{p}\left(t\right)-f^{\prime }\left({\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\left(t\right)\right)}^2+{\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\left(t\right)\right)}^2\right)\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\left(t\right),\hfill \\ \hfill \mathsf{\omega }{\left(x\right)}^{\top }\dot{\mathit{p}}\left(t\right)& =& -\mathsf{\omega }{\left(x\right)}^{\top }D\mathit{q}\left(t\right)+f^{\prime }\left({\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\left(t\right)\right)}^2+{\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\left(t\right)\right)}^2\right)\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\left(t\right),\hfill \end{array}$$

with

$$\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\left(0\right)=u_0\left(x\right),\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\left(0\right)=v_0\left(x\right),x\in [a,b].$$

Considering that

$${\int }_a^b\mathsf{\omega }\left(x\right)\mathsf{\omega }{\left(x\right)}^{\top }\mathrm{d}x=I,$$

(15)

the infinite dimensional identity matrix, left multiplication of the above equations by $\mathsf{\omega }\left(x\right)$ then gives the following infinite system of ODEs (hereafter, we shall omit the argument t for $\mathit{q}$ and $\mathit{p}$):

$$\begin{array}{ccc}\hfill \dot{\mathit{q}}& =& D\mathit{p}-{\int }_a^b\mathsf{\omega }\left(x\right)f^{\prime }\left({\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\right)}^2+{\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\right)}^2\right)\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\mathrm{d}x,\hfill \\ \hfill \dot{\mathit{p}}& =& -D\mathit{q}+{\int }_a^b\mathsf{\omega }\left(x\right)f^{\prime }\left({\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\right)}^2+{\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\right)}^2\right)\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\mathrm{d}x,\hfill \\ & & t\in [0,T].\hfill \end{array}$$

(16)

with the initial conditions:

$$\mathit{q}\left(0\right)={\int }_a^b\mathsf{\omega }\left(x\right)u_0\left(x\right)\mathrm{d}x,\mathit{p}\left(0\right)={\int }_a^b\mathsf{\omega }\left(x\right)v_0\left(x\right)\mathrm{d}x.$$

(17)

The following result holds true.

Theorem 1. 

The ODE system (16) is Hamiltonian w.r.t. the Hamiltonian

$$H(\mathit{q},\mathit{p})=\frac{1}{2}\left[{\mathit{q}}^{\top }D\mathit{q}+{\mathit{p}}^{\top }D\mathit{p}-{\int }_a^{b}f\left({\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\right)}^2+{\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\right)}^2\right)\mathrm{d}x\right].$$

(18)

This latter, in turn, is equivalent to the Hamiltonian functional (5).

Proof. 

Since

$$\dot{\mathit{q}}=\frac{\partial }{\partial \mathit{p}}H(\mathit{q},\mathit{p}),\dot{\mathit{p}}=-\frac{\partial }{\partial \mathit{q}}H(\mathit{q},\mathit{p}),$$

the first part of the statement clearly holds true. The second part follows from (12)–(15), considering that $D=D_1^{\top }{D}_1$, so that

$$\begin{array}{ccc}\hfill {\mathit{q}}^{\top }D\mathit{q}& =& {\mathit{q}}^{\top }{D}_1^{\top }{D}_1\mathit{q}={\int }_a^b{\mathit{q}}^{\top }{D}_1^{\top }\mathsf{\omega }\left(x\right)\mathsf{\omega }{\left(x\right)}^{\top }{D}_1\mathit{q}\mathrm{d}x\hfill \\ & =& {\int }_a^b{\left(\mathsf{\omega }{\left(x\right)}^{\top }{D}_1\mathit{q}\right)}^2\mathrm{d}x={\int }_a^b{u}_x{(x,t)}^2\mathrm{d}x,\hfill \end{array}$$

and, similarly,

$${\mathit{p}}^{\top }D\mathit{p}={\int }_a^b{v}_x{(x,t)}^2\mathrm{d}x.$$

  □

By using similar arguments, it is possible to prove the following result.

Theorem 2. 

The two functionals (8) and (9) can be rewritten, respectively, as:

$$M_1(\mathit{q},\mathit{p})={\int }_a^b{\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{q}\right)}^2+{\left(\mathsf{\omega }{\left(x\right)}^{\top }\mathit{p}\right)}^2\mathrm{d}x,M_2(\mathit{q},\mathit{p})=-2{\mathit{q}}^{\top }{D}_1\mathit{p},$$

(19)

with $D_1$ the infinite matrix defined in (14).

2.2. Truncating the infinite expansion

Clearly, we cannot actually solve the infinite system of ODEs (16): in order to obtain a computational procedure, we need to truncate the infinite expansions (10) to finite sums. I.e., for a conveniently large N,

$$u(x,t)=\sum _{j=0}^{2N}{q}_j\left(t\right){\omega }_j\left(x\right),v(x,t)=\sum _{j=0}^{2N}{p}_j\left(t\right){\omega }_j\left(x\right),$$

(20)

where, for sake of brevity, hereafter we continue to use the same symbols u and v to denote the truncated approximations. Consequently, the infinite vectors in (11) are replaced with the following truncated ones,

$$\mathsf{\omega }\left(x\right)=\left(\begin{array}{c}{\omega }_0\left(x\right)\\ {\omega }_1\left(x\right)\\ ⋮\\ {\omega }_{2N}\left(x\right)\end{array}\right),\mathit{q}\left(t\right)=\left(\begin{array}{c}{q}_0\left(t\right)\\ q_1\left(t\right)\\ ⋮\\ q_{2N}\left(t\right)\end{array}\right),\mathit{p}\left(t\right)=\left(\begin{array}{c}{p}_0\left(t\right)\\ p_1\left(t\right)\\ ⋮\\ p_{2N}\left(t\right)\end{array}\right).$$

(21)

In so doing (12) and (13) continue formally to hold, upon replacing the infinite matrices (14) with the following truncated ones of dimension $(2N+1)\times (2N+1):$

$$D_1=\frac{2\pi }{b-a}\left(\begin{array}{c}0\\ & 1·J\\ & & \ddots \\ & & & N·J\end{array}\right),D=\frac{4{\pi }^2}{{(b-a)}^2}\left(\begin{array}{c}0\\ & 1^2·I_2\\ & & \ddots \\ & & & N^2·I_2\end{array}\right).$$

(22)

As is clear, the functions $u(x,t)$ and $v(x,t)$ in (20) belong, for each fixed t, to the finite-dimensional functional subspace  ${\mathcal{V}}_N=\mathrm{span}\left\{{\omega }_0\left(x\right),{\omega }_1\left(x\right),\dots ,{\omega }_{2N}\left(x\right)\right\}$,  and will no longer, in general, satisfy the original continuous equations (3). However, by requiring the corresponding residual be orthogonal to ${\mathcal{V}}_N$, one formally retrieves again the equations (16)-(17) which, now, amount to a Hamiltonian system of $(4N+2)$ ODEs with the Hamiltonian formally still given by (18) and mass and momentum given by (19), respectively.

We shall discuss the efficient numerical solution of the new finite-dimensional ODE-IVP (16)-(17) in the next sections. Preliminarily, we mention that:

• the value of N is generally very large, in view of obtaining a spectrally accurate space semi-discretization;
• having fixed N, the composite trapezoidal rule over the equally spaced points

  $$x_i=a+i\frac{b-a}{m},i=0,1,\dots ,m,m:=2N+1,$$

  (23)

  is conveniently used to retrieve the Fourier coefficients in (20). In fact, as an example,

  $$q_j\left(t\right)={\int }_a^b{\omega }_j\left(x\right)u(x,t)\mathrm{d}x=\sum _{k=0}^{2N}{q}_k\left(t\right){\int }_a^b{\omega }_j\left(x\right){\omega }_k\left(x\right)\mathrm{d}x,$$

  with the latter integrands being, by virtue of the prosthaphaeresis formulae, trigonometric polynomials of degree at most $2N$. Consequently, the composite trapezoidal rule over $2N+2\equiv m+1$ equally spaced points is exact for computing them [93, page 155].3 Consequently, with reference to (23), and taking into account of the periodicity of the functions, one has:

  $${\int }_a^b{\omega }_j\left(x\right){\omega }_k\left(x\right)\mathrm{d}x=\frac{1}{m}\sum _{i=0}^{m-1}{\omega }_j\left(x_i\right){\omega }_k\left(x_i\right)={\delta }_{jk},\forall j,k=0,\dots ,2N\equiv m-1.$$
  Similar arguments apply, of course, for retrieving $p_j\left(t\right)$.

For sake of completeness, let us give some more operative detail about the latter point. Because of the periodic boundary conditions, according to what exposed above, with reference to (20) one obtains (see (23)):

$$q_j\left(t\right)=\frac{1}{m}\sum _{i=0}^{m-1}{\omega }_j\left(x_i\right)u(x_i,t),p_j\left(t\right)=\frac{1}{m}\sum _{i=0}^{m-1}{\omega }_j\left(x_i\right)v(x_i,t),j=0,\dots ,2N.$$

(24)

These quantities can be efficiently computed via the inverse DFT (e.g., the Matlab^© function ifft) of the vectors

$$\mathit{u}=(u(x_0,t),\dots ,u(x_{m-1},t))\mathrm{and}\mathit{v}=(v(x_0,t),\dots ,v(x_{m-1},t)),$$

(25)

respectively. Specifically, considering that $m-1=2N$, if

$$({\varphi }_0,\dots ,{\varphi }_{2N}):=\mathtt{ifft}\left(\mathit{u}\right)\mathrm{and}({\xi }_0,\dots ,{\xi }_{2N}):=\mathtt{ifft}\left(\mathit{v}\right),$$

then, by using a Matlab^© -like formalism, and avoiding for sake of brevity the argument t,

$$q_0={\varphi }_0,q_{2j-1}=\sqrt{2}·\mathtt{imag}\left({\varphi }_j\right),q_{2j}=\sqrt{2}·\mathtt{real}\left({\varphi }_j\right),j=1,\dots ,N,$$

and, similarly,

$$p_0={\xi }_0,p_{2j-1}=\sqrt{2}·\mathtt{imag}\left({\xi }_j\right),p_{2j}=\sqrt{2}·\mathtt{real}\left({\xi }_j\right),j=1,\dots ,N.$$

Conversely, if (see (23))

$$\mathsf{\mu }=({\mu }_0,\dots ,{\mu }_{m-1}):=\mathtt{real}\left(\mathtt{fft}([0,q_2+\mathrm{i}{q}_1,\dots ,q_{2N}+\mathrm{i}{q}_{2N-1}],m)\right)$$

and

$$\mathsf{\chi }=({\chi }_0,\dots ,{\chi }_{m-1}):=\mathtt{real}\left(\mathtt{fft}([0,p_2+\mathrm{i}{p}_1,\dots ,p_{2N}+\mathrm{i}{p}_{2N-1}],m)\right),$$

then (see (25)),

$$\mathit{u}=q_0+\sqrt{2}·\mathsf{\mu },\mathrm{and}\mathit{v}=p_0+\sqrt{2}·\mathsf{\chi }.$$

One then concludes that the evaluation of the Fourier coefficients, as well as the reconstruction of the approximated solution, can be done with a complexity of $O(mlogm)$ floating-point operations: this, in turn, allows using large values of N.

3. Hamiltonian Boundary Value Methods

In order to obtain a fully discrete method, we now consider the numerical solution of problem (16)-(17). In view of a more compact representation, let us introduce the block vector (see (21) and (23))

$$\mathit{y}=\left(\begin{array}{c}\mathit{q}\\ \mathit{p}\end{array}\right)\in {\mathbb{R}}^{2m},$$

(26)

and denote by $H\left(\mathit{y}\right):=H(\mathit{q},\mathit{p})$ the Hamiltonian (18). Consequently, we can rewrite (16)-(17) as

$$\dot{\mathit{y}}=J_m\nabla H\left(\mathit{y}\right),t\in [0,T],\mathit{y}\left(0\right)={\mathit{y}}_0:=\left(\begin{array}{c}\mathit{q}\left(0\right)\\ \mathit{p}\left(0\right)\end{array}\right),$$

(27)

having set (see (6))

$$J_m:=J\otimes I_m\equiv \left(\begin{array}{cc}& I_m\\ -I_m\end{array}\right)=-J_m^{\top }.$$

Since

$$\frac{\mathrm{d}}{\mathrm{d}t}H\left(\mathit{y}\right)=\nabla H{\left(\mathit{y}\right)}^{\top }\dot{\mathit{y}}=\nabla H{\left(\mathit{y}\right)}^{\top }{J}_m\nabla H\left(\mathit{y}\right)=0,$$

(28)

due to the skew-symmetry of $J_m$, we again retrieve the conservation of H and, hence, our aim is to derive energy-conserving methods, i.e., able to conserve the energy along the numerical trajectory. Such methods are derived within the so called framework of line integral methods (see, e.g., the monograph [67] or the review paper [72] ). In more details, if h is the considered timestep, we shall look for a suitable path $\sigma$ such that:

$$\sigma \left(0\right)={\mathit{y}}_0,\sigma \left(h\right)=:{\mathit{y}}_1\approx \mathit{y}\left(h\right),{\int }_0^1\nabla H{\left(\sigma \left(ch\right)\right)}^{\top }\dot{\sigma }\left(ch\right)\mathrm{d}c=0.$$

(29)

These properties imply energy-conservation, since

$$H\left({\mathit{y}}_1\right)-H\left({\mathit{y}}_0\right)=H\left(\sigma \left(h\right)\right)-H\left(\sigma \left(0\right)\right)=h{\int }_0^1\nabla H{\left(\sigma \left(ch\right)\right)}^{\top }\dot{\sigma }\left(ch\right)\mathrm{d}c=0.$$

Remark 1. 

We observe that for the continuous solution (28) the last property in (29) derives from the fact that the integrand identically vanishes. This is no more the case for the path σ characterizing the given method (29).

A straightforward way to derive paths satisfying (29) relies on the expansion of the vector field (27) along a suitable orthonormal basis, which we choose as the Legendre polynomial basis ${\left\{P_j\right\}}_{j\ge 0}$ [68,71]:

$$P_i\in {\mathsf{\Pi }}_i,{\int }_0^1{P}_i\left(x\right)P_j\left(x\right)\mathrm{d}x={\delta }_{ij},i,j=0,1,2,\dots ,$$

(30)

where, as usual, ${\mathsf{\Pi }}_i$ denotes the vector space of polynomials of degree i. In so doing, one obtains:

$$\dot{\mathit{y}}\left(ch\right)=\sum _{j\ge 0}{P}_j\left(c\right){\gamma }_j\left(\mathit{y}\right),c\in [0,1],\mathit{y}\left(0\right)={\mathit{y}}_0,$$

(31)

with the Fourier coefficients given by:

$${\gamma }_j\left(\mathit{y}\right)=J_m{\int }_0^1{P}_j\left(\tau \right)\nabla H\left(\mathit{y}\left(\tau h\right)\right)\mathrm{d}\tau ,j=0,1,2,\dots .$$

(32)

Integrating (31) side by side, and imposing the initial condition, one then obtains:

$$\mathit{y}\left(ch\right)={\mathit{y}}_0+h\sum _{j\ge 0}{\int }_0^c{P}_j\left(x\right)\mathrm{d}x{\gamma }_j\left(\mathit{y}\right),c\in [0,1].$$

(33)

In particular, by taking into account (30) and (27),

$$\mathit{y}\left(h\right)={\mathit{y}}_0+h{\gamma }_0\left(\mathit{y}\right)\equiv {\mathit{y}}_0+h{J}_m{\int }_0^1\nabla H\left(\mathit{y}\left(\tau h\right)\right)\mathrm{d}\tau \equiv {\mathit{y}}_0+{\int }_0^h\dot{\mathit{y}}\left(t\right)\mathrm{d}t,$$

i.e., one retrieves the usual Fundamental Theorem of Calculus. Concerning the Fourier coefficients (32), the following general result holds true.

Lemma 1. 

Let $G:[0,h]\to V$, with V a vector space, admit a Taylor expansion at 0. Then,

$${\int }_0^1{P}_j\left(\tau \right)G\left(\tau h\right)\mathrm{d}\tau =O\left(h^j\right),j=0,1,2,\dots .$$

Proof. 

See [71, Lemma 1]. □

More in general, it can be proved the following, more general, result whose proof is omitted.

Lemma 2. 

Let $G:[0,h]\to V$, with V a vector space, admit a Taylor expansion at 0 with remainder at the $O\left(h^k\right)$ term. Then,

$${\int }_0^1{P}_j\left(\tau \right)G\left(\tau h\right)\mathrm{d}\tau =O\left(h^{min\{j,k\}}\right),j=0,1,2,\dots .$$

In order to obtain a polynomial approximation $\sigma \in {\mathsf{\Pi }}_s$ to $\mathit{y}\left(ch\right)$, it is enough truncating the infinite series in (31) and (33) to finite sums with s terms:

$$\dot{\sigma }\left(ch\right)=\sum _{j=0}^{s-1}{P}_j\left(c\right){\gamma }_j\left(\sigma \right),c\in [0,1],\sigma \left(0\right)={\mathit{y}}_0,$$

(34)

and

$$\sigma \left(ch\right)={\mathit{y}}_0+h\sum _{j=0}^{s-1}{\int }_0^c{P}_j\left(x\right)\mathrm{d}x{\gamma }_j\left(\sigma \right),c\in [0,1].$$

(35)

with ${\gamma }_j\left(\sigma \right)$ defined, similarly as in (32), as

$${\gamma }_j\left(\sigma \right)=J_m{\int }_0^1{P}_j\left(\tau \right)\nabla H\left(\sigma \left(\tau h\right)\right)\mathrm{d}\tau ,j=0,\dots ,s-1.$$

(36)

In such a case, according to (29), the new approximation is given by

$${\mathit{y}}_1:=\sigma \left(h\right)={\mathit{y}}_0+h{\gamma }_0\left(\sigma \right).$$

(37)

As stated in the next theorem, the approximation procedure is energy-conserving.

Theorem 3. 

With reference to (18), (26)-(27), and (29), one has: $H\left({\mathit{y}}_1\right)=H\left({\mathit{y}}_0\right)$.

Proof. 

By virtue of (34), and taking into account (36), one obtains:

$$\begin{array}{c}{\displaystyle H\left({\mathit{y}}_1\right)-H\left({\mathit{y}}_0\right)=H\left(\sigma \left(h\right)\right)-H\left(\sigma \left(0\right)\right)=h{\int }_0^1\nabla H{\left(\sigma \left(ch\right)\right)}^{\top }\dot{\sigma }\left(ch\right)\mathrm{d}c}\\ =h{\int }_0^1\nabla H{\left(\sigma \left(ch\right)\right)}^{\top }\sum _{j=0}^{s-1}{P}_j\left(c\right){\gamma }_j\left(\sigma \right)\mathrm{d}c\hfill \\ =h\sum _{j=0}^{s-1}{\left[{\int }_0^1{P}_j\left(c\right)\nabla H\left(\sigma \left(ch\right)\right)\mathrm{d}c\right]}^{\top }{J}_m\left[{\int }_0^1{P}_j\left(c\right)\nabla H\left(\sigma \left(ch\right)\right)\mathrm{d}c\right]=0,\hfill \end{array}$$

due to the skew-symmetry of $J_m$. □

Further, the following result can be proved, stating that we have derived an order $2s$ approximation procedure.

Theorem 4. 

With reference to (26) and (29), one has: ${\mathit{y}}_1=\mathit{y}\left(h\right)+O\left(h^{2s+1}\right)$.

Proof. 

See [71, Theorem 1]. □

Interestingly enough, the problem of determining $\sigma$ can be cast into the problem of finding the Fourier coefficients ${\gamma }_j\left(\sigma \right)$. In fact, from (35) and (36), one deduces the following system of (generally nonlinear) equations:

$${\gamma }_j\left(\sigma \right)=J_m{\int }_0^1{P}_j\left(c\right)\nabla H\left({\mathit{y}}_0+h\sum _{\ell =0}^{s-1}{\int }_0^c{P}_{\ell }\left(x\right)\mathrm{d}x{\gamma }_{\ell }\left(\sigma \right)\right)\mathrm{d}c,j=0,\dots ,s-1.$$

(38)

Once this system is solved, the new approximation is formally given by (37), as above explained.

3.1. Discretization of the Fourier coefficients

Quoting Dahlquist and Bijörk [94, page 521], as is well known, even many relatively simple integrals cannot be expressed in finite terms of elementary functions, and thus must be evaluated by numerical methods. Within our framework, this quite obvious statement means that we cannot in general directly compute the Fourier coefficients (36): instead, we can approximate them by using a suitable quadrature rule, which we choose as the Gauss-Legendre formula of order $2k$, for a convenient $k\ge s$, with abscissae $\left\{c_i\right\}$ given by the zeros of $P_k$, i.e., such that $P_k\left(c_i\right)=0$, $i=1,\dots ,k$, and weights

$$b_i={\int }_0^1\prod _{j\ne i}\frac{c-c_j}{c_i-c_j}\mathrm{d}c,i=1,\dots ,k.$$

In so doing, we obtain:

$${\widehat{\gamma }}_j:=J_m\sum _{i=1}^k{b}_i{P}_j\left(c_i\right)\nabla H\left(\sigma \left(c_{i}h\right)\right)\equiv {\gamma }_j\left(\sigma \right)+{\Delta }_j\left(h\right),j=0,\dots ,s-1,$$

(39)

with ${\Delta }_j\left(h\right)$ the quadrature error. For sake of brevity, we shall continue to use the same symbol $\sigma$ for the new polynomial approximation obtained by formally substituting ${\widehat{\gamma }}_j$ to ${\gamma }_j\left(\sigma \right)$, $j=0,\dots ,s-1$, in (34) and (35), with the new approximation now given by

$${\mathit{y}}_1:=\sigma \left(h\right)={\mathit{y}}_0+h{\widehat{\gamma }}_0,$$

(40)

in place of (37), and the discrete problem

$${\widehat{\gamma }}_j=J_m\sum _{i=1}^k{b}_i{P}_j\left(c_i\right)\nabla H\left({\mathit{y}}_0+h\sum _{\ell =0}^{s-1}{\int }_0^{c_i}{P}_{\ell }\left(x\right)\mathrm{d}x{\widehat{\gamma }}_{\ell }\right),j=0,\dots ,s-1,$$

(41)

in place of (38).

Definition 1. 

The method defined by (39)–(41) is named Hamiltonian Boundary Value Method with parameters k and s: in short,HBVM$(k,s)$.

Concerning the quadrature error in (39), it is quite straightforward to prove that, assuming a suitably regular Hamiltonian function H, i.e., $f\left(x\right)$ in (18), then:

• in the relevant case where $f\in {\mathsf{\Pi }}_{\nu }$, with $\nu \ge 1$,4 then the integrand in (36) is a polynomial of degree at most:

  $$(\nu -1)2s+2s-1=2\nu s-1.$$
  Consequently, the quadrature is exact, provided that

  $$\nu \le k/s;$$

  (42)
• conversely, since the quadrature error is proportional to the $2k$-th derivative of the integrand, then

  $${\Delta }_j\left(h\right)=O\left(h^{2k-j}\right),j=0,\dots ,s-1.$$

  (43)

The above arguments, allows us to state the following conservation result.

Theorem 5. 

With reference to (18) and (26), if $f\in {\mathsf{\Pi }}_{\nu }$ with ν satisfying (42), then $H\left({\mathit{y}}_1\right)=H\left({\mathit{y}}_0\right)$. Conversely, $H\left({\mathit{y}}_1\right)=H\left({\mathit{y}}_0\right)+O\left(h^{2k+1}\right)$

Proof. 

The first part of the statement follows from Theorem 3, due to the fact that the quadrature is exact. The second part follows similarly, by considering that, by virtue of Lemma 1,

$$\left[{\int }_0^1{P}_j\left(c\right)\nabla H\left(\sigma \left(ch\right)\right)\mathrm{d}c\right]=O\left(h^j\right),j=0,\dots ,s-1.$$

Consequently, from (39) and (43) one has:

$$\begin{array}{c}{\displaystyle H\left({\mathit{y}}_1\right)-H\left({\mathit{y}}_0\right)=H\left(\sigma \left(h\right)\right)-H\left(\sigma \left(0\right)\right)=h{\int }_0^1\nabla H{\left(\sigma \left(ch\right)\right)}^{\top }\dot{\sigma }\left(ch\right)\mathrm{d}c}\hfill \\ =h{\int }_0^1\nabla H{\left(\sigma \left(ch\right)\right)}^{\top }\sum _{j=0}^{s-1}{P}_j\left(c\right){\widehat{\gamma }}_j\mathrm{d}c\hfill \\ =h{\int }_0^1\nabla H{\left(\sigma \left(ch\right)\right)}^{\top }\sum _{j=0}^{s-1}{P}_j\left(c\right)\left[{\gamma }_j\left(\sigma \right)+{\Delta }_j\left(h\right)\right]\mathrm{d}c\hfill \\ =h\sum _{j=0}^{s-1}{\left[{\int }_0^1{P}_j\left(c\right)\nabla H\left(\sigma \left(ch\right)\right)\mathrm{d}c\right]}^{\top }{J}_m\left[{\int }_0^1{P}_j\left(c\right)\nabla H\left(\sigma \left(ch\right)\right)\mathrm{d}c\right]\hfill \\ +h\sum _{j=0}^{s-1}{\left[{\int }_0^1{P}_j\left(c\right)\nabla H\left(\sigma \left(ch\right)\right)\mathrm{d}c\right]}^{\top }{\Delta }_j\left(h\right)=h\sum _{j=0}^{s-1}O\left(h^j\right)O\left(h^{2k-j}\right)=O\left(h^{2k+1}\right),\hfill \end{array}$$

due to the skew-symmetry of $J_m$. □

Remark 2. 

From the result of Theorem 5 it follows that:

● either an exact

● or apractical

energy conservation can always be gained, by choosing k large enough. In fact, in the latter case, it is enough that the energy error be within the roundoff error level.

3.2. Runge-Kutta form

The HBVM$(k,s)$ method (39)–(41) can be seen to be a k-stage Runge-Kutta (RK) method. As matter of fact,

$$Y_i:=\sigma \left(c_{i}h\right)={\mathit{y}}_0+h\sum _{\ell =0}^{s-1}{\int }_0^{c_i}{P}_{\ell }\left(x\right)\mathrm{d}x{\widehat{\gamma }}_{\ell },i=1,\dots ,k,$$

(44)

can be regarded as its i-th stage. In fact, by plugging it at the r.h.s. that of (41), and considering (44), one obtains (by slightly adapting the indices) :

$$\begin{array}{ccc}\hfill Y_i& =& {\mathit{y}}_0+h\sum _{\ell =0}^{s-1}{\int }_0^{c_i}{P}_{\ell }\left(x\right)\mathrm{d}x\sum _{j=1}^k{b}_j{P}_{\ell }\left(c_j\right)J_m\nabla H\left(Y_j\right)\hfill \\ & =& {\mathit{y}}_0+h\sum _{j=1}^k\underset{=:a_{ij}}{\underbrace{b_j\left[\sum _{\ell =0}^{s-1}{\int }_0^{c_i}{P}_{\ell }\left(x\right)\mathrm{d}x{P}_{\ell }\left(c_j\right)\right]}}{J}_m\nabla H\left(Y_j\right),i=1,\dots ,k,\hfill \end{array}$$

(45)

with the new approximation given, considering that $P_0\left(c\right)\equiv 1$, by:

$${\mathit{y}}_1={\mathit{y}}_0+h\sum _{i=1}^k{b}_i{J}_m\nabla H\left(Y_i\right).$$

(46)

Consequently, by setting $\mathbf{b}={(b_1,\dots ,b_k)}^{\top }$ and $\mathbf{c}={(c_1,\dots ,c_k)}^{\top }$ the vectors of the weights and the abscissae of the quadrature, and $A=\left(a_{ij}\right)\in {\mathbb{R}}^{k\times k}$ (see (45)) the Butcher matrix, we have derived the k-stage RK method

(47)

Moreover, concerning the Butcher matrix, the following result holds true.

Theorem 6. 

Setting the matrices

$$\begin{array}{ccc}\hfill {\mathcal{P}}_s& =& \left(\begin{array}{ccc}{P}_0\left(c_1\right)& \dots & P_{s-1}\left(c_1\right)\\ ⋮& & ⋮\\ P_0\left(c_k\right)& \dots & P_{s-1}\left(c_k\right)\end{array}\right),\Omega =\left(\begin{array}{c}{b}_1\\ & \ddots \\ & & b_k\end{array}\right),\hfill \\ \hfill {\mathcal{I}}_s& =& \left(\begin{array}{ccc}{\int }_0^{c_1}{P}_0\left(x\right)\mathrm{d}x& \dots & {\int }_0^{c_1}{P}_{s-1}\left(x\right)\mathrm{d}x\\ ⋮& & ⋮\\ {\int }_0^{c_k}{P}_0\left(x\right)\mathrm{d}x& \dots & {\int }_0^{c_k}{P}_{s-1}\left(x\right)\mathrm{d}x\end{array}\right),\hfill \end{array}$$

(48)

one has that the Butcher matrix in (47) is given by:

$$A={\mathcal{I}}_s{\mathcal{P}}_s^{\top }\Omega .$$

(49)

Proof. 

The statement immediately follows by considering the $(i,j)$-th entry of the matrices at each member of (49). □

By defining the block vectors

$$Y=\left(\begin{array}{c}{Y}_1\\ ⋮\\ Y_k\end{array}\right),\nabla H\left(Y\right)=\left(\begin{array}{c}\nabla H\left(Y_1\right)\\ ⋮\\ \nabla H\left(Y_k\right)\end{array}\right),$$

and $e={(1,\dots ,1)}^{\top }\in {\mathbb{R}}^k$, the stage equations (45) can be cast in vector form, by virtue of (49), as

$$Y=e\otimes {\mathit{y}}_0+h{\mathcal{I}}_s{\mathcal{P}}_s^{\top }\Omega \otimes J_m\nabla H\left(Y\right).$$

(50)

This latter equation can be, in turn, rewritten as

$$Y=e\otimes {\mathit{y}}_0+h{\mathcal{I}}_s\otimes I_{2m}\widehat{\mathbf{\gamma }},$$

having set (see (41) and (44))

$$\widehat{\mathbf{\gamma }}\equiv \left(\begin{array}{c}{\widehat{\gamma }}_0\\ ⋮\\ {\widehat{\gamma }}_{s-1}\end{array}\right)={\mathcal{P}}_s^{\top }\Omega \otimes J_m\nabla H\left(Y\right),$$

the block vector with the (approximate) Fourier coefficients (41). By combining the last two equations, one then obtains

$$\widehat{\mathbf{\gamma }}={\mathcal{P}}_s^{\top }\Omega \otimes J_m\nabla H\left(e\otimes {\mathit{y}}_0+h{\mathcal{I}}_s\otimes I_{2m}\widehat{\mathbf{\gamma }}\right).$$

(51)

Remark 3. 

In other words, even though a HBVM$(k,s)$ method is a k-stage RK method and, therefore, its stage equation (50) has (block)-dimension k, the actual discrete problem to be solved at each time-step can be cast as (51), thus having (block)-dimension s,independentlyof k.

What exposed above, in turn, allows using relatively large values of k, in view of getting an (either exact of practical) conservation of the Hamiltonian, according to Theorem 5. The efficient solution of the discrete problem (51) will be discussed in Section 4.

We conclude this section by reporting the following additional result, concerning a HBVM$(k,s)$ method. For the corresponding proofs we refer, .e.g., to the monograph [67].

Theorem 7. 

For all $k\ge s$, a HBVM$(k,s)$ method is symmetric and has order $2s$. When $k=s$ it reduces to the s-stage Gauss-collocation method.

3.3. HBVMs as spectral methods in time

An interesting interpretation of HBVMs is that they can be regarded as spectral methods in time. In more detail, let us recall the expansion (31)-(32). Assuming that $\dot{\mathit{y}}\in L^2[0,h]$, it then follows that

$$\parallel \dot{\mathit{y}}{\parallel }_{L^2}^2=\sum _{j\ge 0}\parallel {\gamma }_j{\left(\mathit{y}\right)\parallel }^2<\infty \Rightarrow \parallel {\gamma }_j\left(\mathit{y}\right)\parallel \to 0,j\to \infty .$$

Consequently, if we are using a finite precision arithmetic with machine epsilon $\epsilon$, and we choose the degree s of the polynomial approximation $\sigma$ in (34)–(36) such that

$$\forall j\ge s:\parallel {\gamma }_j\left(\mathit{y}\right)\parallel <\epsilon ·\underset{i=0,\dots ,s-1}{max}\parallel {\gamma }_i\left(\mathit{y}\right)\parallel ,$$

(52)

it then follows that

$$\sigma \left(ch\right)\doteq \mathit{y}\left(ch\right),c\in [0,1],$$

(53)

where ≐ means equal within roundoff errors.

Definition 2. 

We shall refer to a HBVM$(k,s)$ method such that k and s, with $k>s$, are large enough such that (see (36) and (39)) both

$${\widehat{\gamma }}_j\doteq {\gamma }_j\left(\sigma \right),j=0,\dots ,s-1,$$

and (52) hold true, so that (53) follows, as a Spectral HBVMor, in short, SHBVM.5

The use of HBVMs as spectral methods in time has been studied in [95]. In particular, the choice of s has been investigated in [96], whereas, hereafter we shall consider the following choice for k [89]:

$$k=max\{20,s+2\}.$$

(54)

4. Efficient application of HBVMs to the NLSE

Solving the discrete problem (51) could be, in general, a severe computational issue, because of many reasons:

• the high dimensionality, $2m=4N+2$ (see (23), (24), and (26)), of the state space, if a spectrally accurate space semi-discretization is targeted;
• the possible large (block)-size s of $\widehat{\mathbf{\gamma }}$, in view of the use of HBVMs as spectral methods in time;
• the unpractical use of (51) to derive a straightforward fixed-point iteration, which would require the use of a very small stepsize h.

Consequently, we need an efficient Newton-type iteration for this purpose. To begin with, let us start considering the simplified Newton iteration for solving

$$G\left(\widehat{\mathbf{\gamma }}\right):=\widehat{\mathbf{\gamma }}-{\mathcal{P}}_s^{\top }\Omega \otimes J_m\nabla H\left(e\otimes {\mathit{y}}_0+h{\mathcal{I}}_s\otimes I_{2m}\widehat{\mathbf{\gamma }}\right)=\mathbf{0},$$

(55)

which, in view of the fact that (see (48))

$${\mathcal{P}}_s^{\top }\Omega {\mathcal{I}}_s=\left(\begin{array}{cc}{\xi }_0& -{\xi }_1\\ {\xi }_1& 0& \ddots \\ & \ddots & \ddots & -{\xi }_{s-1}\\ & & {\xi }_{s-1}& 0\end{array}\right)=:X_s,{\xi }_i=\frac{1}{2\sqrt{|4i^2-1|}},$$

(56)

reads

$$\left[I_{2sm}-h{X}_s\otimes \left(J_m{\nabla }^{2}H\left({\mathit{y}}_0\right)\right)\right]\Delta {\widehat{\mathbf{\gamma }}}^r=-G\left({\widehat{\mathbf{\gamma }}}^r\right),r=0,1,2,\dots ,$$

(57)

with ${\widehat{\mathbf{\gamma }}}^0=\mathbf{0}$ a convenient initial guess. In this form, the iteration (57) would require the factorization of a $2sm\times 2sm$ matrix at each time-step, which would be too expensive.

In this respect, a first improvement can be obtained by observing that the problem (16) has a leading linear part, according to (22), when N is large (as in the case when dealing with a spectrally accurate space semi-discretization). Consequently, we could consider the approximation

$${\nabla }^{2}H\left({\mathit{y}}_0\right)\approx J\otimes D,$$

with matrix J defined in (6). In so doing (57) simplifies to

$$\left[I_{2sm}-h{X}_s\otimes J\otimes D\right]\Delta {\widehat{\mathbf{\gamma }}}^r=-G\left({\widehat{\mathbf{\gamma }}}^r\right),r=0,1,\dots ,$$

(58)

which has the advantage of having the same coefficient matrix for all time-steps. As is clear, however, this matrix has still dimension $2sm\times 2sm$, which can be again unpractical, when using HBVMs as SHBVMs, and considering that its factors are in general full.

The second, and more decisive, improvement can be obtained via a blended iteration for solving (58). The basic idea, which has been initially developed in a series of papers (see [97,98,99,100]) focused on block implicit methods for ODEs, and adapted in [101] (see also [102]) for HBVMs, relies on the combination of the equation (58) with an equivalent formulation, i.e.,

$${\rho }_s\left[X_s^{-1}\otimes I_{2m}-h{I}_s\otimes J\otimes D\right]\Delta {\widehat{\mathbf{\gamma }}}^r=-\left({\rho }_s{X}_s^{-1}\otimes I_{2m}\right)G\left({\widehat{\mathbf{\gamma }}}^r\right),r=0,1,\dots ,$$

(59)

with the scalar parameter ${\rho }_s$ defined as 6

$${\rho }_s=\underset{\lambda \in \sigma \left(X_s\right)}{min}\left|\lambda \right|,$$

(60)

combined (i.e., blended) with weights

$$I_s\otimes \Gamma \mathrm{and}{I}_s\otimes (I_{2m}-\Gamma ),$$

respectively, where

$$\Gamma ={\left(I_{2m}-h{\rho }_{s}J\otimes D\right)}^{-1}.$$

The good news is now twofold:

• in so doing, the inverse of the coefficient matrix, required by the iterative procedure, can be seen to be given by $I_s\otimes \Gamma$, i.e., the same matrix defined above;
• moreover, $\Gamma$ is a very sparse matrix. As matter of fact, one directly verifies that:

  $$\Gamma =I_2\otimes {(I_m+B^2)}^{-1}+J\otimes \left[B{(I_m+B^2)}^{-1}\right],B=h{\rho }_{s}D,$$

  (61)

  which can be stored in two vectors, containing the main diagonals of ${(I_m+B^2)}^{-1}$ and $B{(I_m+B^2)}^{-1}$, respectively.

For sake of completeness, the resulting final algorithm is listed in Table 1, where we assume that $\Gamma$ has been computed (once for all time-steps) through (60)-(61).

Remark 4. 

It must be stressed that, in order to obtain energy-conservation and spectral accuracy in time, the blendediteration in Table 1 has to be carried out until full machine accuracy is gained. This is, indeed, the implementation considered in the numerical tests.

5. Numerical examples

We now report a few numerical results, aimed at assessing the theoretical findings, and proving the effectiveness of HBVMs, especially when used as SHBVMs. All numerical tests have been carried out in Matlab^© (Rel. 2023b) using a Silicon M2-based computer, with 16 GB of shared memory.

5.1. Example 1

The first test problem is given by the so called focusing case of the classical NLSE, for which $f\left(x\right)=x^2$ in (1). In such a case, the initial condition at $t=0$ is taken from the (known) solution,

$$\psi (x,t)=\mathrm{sech}(x-x_0-2\kappa t)exp\left(\mathrm{i}((1-{\kappa }^2)t+\kappa x)\right),(x,t)\in [-160,160]\times [0,20],$$

(62)

where we have used the parameters $x_0=-100$ and $\kappa =5$. Such solution, named bright soliton [18,103], is depicted in Figure 1.

We observe that, since $f\left(x\right)$ is a polynomial of degree 2, then any HBVM$(2s,s)$ is exactly energy-conserving and of order $2s$, according to Theorems 5 and 7. In particular, in such a case the HBVM(2,1) method coincides with the so called AVF method[70], used for solving the semi-discrete Hamiltonian problem (16)-(17). We shall solve the problem over the given time interval $[0,20]$ by using the following methods (all implemented in the same code):

• HBVM(2,1) (i.e., the AVF method), with stepsizes $h=2^{-n}{10}^{-1}$, $n=3,4\dots ,12$;
• HBVM(4,2), with stepsizes $h=2^{-n}{10}^{-1}$, $n=1,2\dots ,8$;
• HBVM(6,3), with step sizes $h=2^{-n}{10}^{-1}$, $n=0,1\dots ,6$;
• HBVM(20,18), which is spectrally accurate in time, when using the stepsize $h=0.1$ (i.e., it provides a SHBVM).

For all such methods we use a value of $N=1200$ in (20)–(23). Consequently, we obtain a system of 4802 Hamiltonian ODEs (16). By using this value of N, the space semi-discretization error falls within the double precision IEEE roundoff error level. Further, the constant values of Hamiltonian, mass, and momentum in (18) and (19) are respectively given by:

$$H=\frac{74}{3},M_1=2,M_2=-\frac{1}{32}.$$

(63)

We also compare the previous methods with the two explicit symplectic composition methods described in [37], coupling the standard second-order difference of the second space-derivative, with parameter $\Delta x$, with a composition method of order 1 (S1) or 2 (S2) in time, so that, when using a stepsize $\Delta t$, the error is either $O(\Delta x^2+\Delta t)$ or $O(\Delta x^2+\Delta t^2)$, respectively. However, because of stability constraints, $\Delta t$ is not independent of $\Delta x$. In particular, we use $\Delta t=\frac{1}{2}\Delta x^2$. Consequently, both methods are actually only first order accurate, since the error is $O(\Delta x^2)=O(\Delta t)$ for both of them. In light of this, we consider the methods S1 and S2 with discretization parameters:

$$\Delta x=\frac{0.32}{n},\Delta t=\frac{1}{2}{\left(\frac{0.32}{n}\right)}^2,n=5,10,20,40,80,$$

so that we have a space semi-discretization with ${10}^{3}n+1$ points, for a total of $[20/\Delta t]$ timesteps, where $[·]$ denotes the rounding of the argument. For all methods, we compare the errors in the last timestep.

The two plots in Figure 2 show the obtained results:

• in the left-plot, we have the error versus the used stepsize, showing the predicted order:

  -

  1 for S1 and S2,

  -

  2 for HBVM(2,1) (i.e., the AVF method),

  -

  4 for HBVM(4,2),

  -

  6 for HBVM(6,3).

  Clearly, the SHBVM reaches full machine accuracy for the considered stepsize;
• in the right-plot, we have the so called work-precision diagram, where the errors are plotted against the execution times (in sec). As one may see, the higher the order of the method, the better its performance, with the SHBVM reaching full accuracy in a moderate time.

From the obtained results, it clearly follows that the SHBVM is the method with the best performance. This conclusion is further confirmed by the plots in Figure 3, where there are:

• the relative errors in the Hamiltonian, mass, and momentum for the AVF and HBVM(4,2) methods used with stepsize $h={10}^{-2}$ (upper plots);
• the relative errors in the invariants for the SHBVM used with stepsize $h=0.1$ (lower left-plot), and the modulus of the absolute error (lower right-plot).

From the pictures, one can see that all methods are energy-conserving, as expected. However, the AVF and HBVM(4,2) methods only approximately conserve the other invariants (though with no drift), unlike the SHBVM, for which they are conserved up to roundoff. For this latter method, also the solution error can be seen to be within the roundoff error level.

Figure 3. Example 1: errors in the invariants for the AVF and HBVM(4,2) methods used with stepsize $h={10}^{-2}$ (upper plots); errors in the invariants (lower left-plot) and in the solution (lower right-plot) for SHBVM used with the stepsize $h=0.1$.

Figure 3. Example 1: errors in the invariants for the AVF and HBVM(4,2) methods used with stepsize $h={10}^{-2}$ (upper plots); errors in the invariants (lower left-plot) and in the solution (lower right-plot) for SHBVM used with the stepsize $h=0.1$.

5.2. Example 2

We now consider the same equation used in the first example, but with a different initial condition. Namely (compare with (62)),

$$\psi (x,0)=\mathrm{sech}(x-x_0)exp\left(\mathrm{i}\kappa x\right)+\mathrm{sech}(x+x_0)exp(-\mathrm{i}\kappa x),x\in [-160,160],$$

(64)

with the same parameters $x_0=-100$ and $\kappa =5$ considered for (62). In this case, we have two bright solitons that collide at $t=10$ and emerge after the crossing. In such a case, the values of the invariants are (compare with (63)):

$$H=\frac{148}{3},M_1=4,M_2=0.$$

(65)

As in the previous example, we solve the problem by using a trigonometric polynomial of degree $N=1200$ for the space semi-discretization, and covering the time interval $[0,20]$ by using a stepsize:

• $h={10}^{-2}$, for the AVF (i.e., HBVM(2,1)) and HBVM(4,2) methods;
• $h=0.1$, for the HBVM(20,18) method, thus obtaining a spectrally accurate space-time method (SHBVM).

The obtained results are depicted in Figure 4, where one may find:7

• the errors in the invariants for the AVF and HBVM(4,2) methods (upper plots), which exactly conserve the energy H and the momentum $M_2$, whereas there is a numerical “peak” in the mass $M_1$, when the two solitons collide;
• the errors in the invariants for the SHBVM method, along with the modulus of the computed solution, where one may see that the two solitons emerge after the collision at $t=10$ (lower plots).

As one may see, the first two methods exhibit an artificial peak in the numerical mass, where the two solitons collide, whereas the SHBVM method has all the errors in the invariants within the roundoff error level and, remarkably enough, there are no peaks at $t=10$, where the two solitons cross.

Figure 4. Example 2: errors in the invariants for the AVF and HBVM(4,2) methods used with stepsize $h={10}^{-2}$ (upper plots, with the Hamiltonian and momentum errors, both within roundoff, overlapping); errors in the invariants and modulus of the computed solution for the SHBVM method used with stepsize $h=0.1$ (lower plots).

Figure 4. Example 2: errors in the invariants for the AVF and HBVM(4,2) methods used with stepsize $h={10}^{-2}$ (upper plots, with the Hamiltonian and momentum errors, both within roundoff, overlapping); errors in the invariants and modulus of the computed solution for the SHBVM method used with stepsize $h=0.1$ (lower plots).

5.3. Example 3

The next test problem is given by the so called defocusing case of the classical NLSE, for which $f\left(x\right)=-x^2$ in (1). It is known that this problem is much more challenging than the focusing case. The initial condition is taken as:

$$\psi (x,0)=(1-\mathrm{sech}\left(x\right))(1-exp(\lambda x^2+\mathrm{i}\theta x)),x\in [-120,120].$$

(66)

We choose the parameters $\lambda =-2$ and $\theta =\frac{1}{2}$. The values of the invariants are given by:8

$$H\approx 118,M_1\approx 236,M_2\approx 0.$$

(67)

We solve this problem by using $N=1200$ for the space semi-discretization and, for the time integration, a HBVM(22,20) with stepsize $h=0.1$ for 200 timesteps (thus, covering the interval $[0,20]$). The method turns out to be a SHBVM, for the considered stepsize. In Figure 5 we plot the obtained result: in the upper-left, the modulus of $\psi (x,t)$; in the upper-right the errors in the invariants (relative, for the Hamiltonian and the mass, absolute for the momentum, according to (67)); in the middle plots are the real and imaginary parts of the solution, u and v, respectively; for comparison, we also plot the errors in the invariants for the AVF (i.e., HBVM(2,1)) and the HBVM(4,2) methods (lower plots), when using a stepsize $h={10}^{-2}$. As is clear, all methods conserve the energy, but the SHBVM exhibits a much better conservation of mass and momentum (always within the roundoff error level).

5.4. Example 4

As a last example, we consider the cubic-quintic nonlinear Schrödinger equation, introduced in the introduction (see also [104]), with $f\left(x\right)=x^2-\frac{1}{3}{x}^3$. We use the initial condition $\psi (x,0)$ from (62), and the same space-time domain $[-160,160]\times [0,20]$. In such a case, any HBVM$(3s,s)$ method turns out to be energy conserving, according to Theorem 5. The values of the invariants are now given by:9

$$H\approx 25,M_1=2,M_2\approx 0.$$

(68)

In Figure 6 we plot the errors in the invariants for:10

• the AVF (i.e., HBVM(3,1)) and the HBVM(6,2) methods used with a stepsize $h=5·{10}^{-3}$ (upper plots);
• the HBVM(20,18) method, which is spectrally accurate in time (i.e., it is a SHBVM), when using a stepsize $h=5·{10}^{-2}$. We also plot the modulus of the computed solution (lower plots).

In all cases, we have used a trigonometric polynomial approximation of degree $N=1200$ which, as in the case of the previous examples, turns out to provide a spectrally accurate space semi-discretization. From the obtained results, one infers that all methods are energy-conserving. However, only the SHBVM has all the errors in the invariants within the roundoff error level, whereas this is not the case for the other methods, though no numerical drift is observed.

Figure 6. Example 4: errors in the invariants for the AVF and HBVM(6,2) methods used with stepsize $h=5·{10}^{-3}$ (upper plots); errors in the invariants (lower left-plot) and modulus of the computed solution (lower right-plot) for SHBVM used with the stepsize $h=5·{10}^{-2}$.

Figure 6. Example 4: errors in the invariants for the AVF and HBVM(6,2) methods used with stepsize $h=5·{10}^{-3}$ (upper plots); errors in the invariants (lower left-plot) and modulus of the computed solution (lower right-plot) for SHBVM used with the stepsize $h=5·{10}^{-2}$.

6. Conclusions

In this paper, we have reviewed in major detail an effective approach for the numerical solution of the NLSE, when equipped with initial and periodic boundary conditions. After a space semi-discretization, obtained by using a classical Fourier expansion (because of the periodic boundary conditions), the time integration is carried out by using energy-conserving methods in the HBVMs class. The actual implementation of the methods, for the problem at hand, has been explained in full details. Also their use as spectral methods in time (SHBVM) has been recalled. In particular, the numerical tests here reported clearly show the advantage of SHBVMs over other existing methods. As a possible development of this approach, we mention the possibility of a parallel-in-time solution of the problem, following the approach in [105,106], as well as the extension to more space dimensions.