An $hr$-Adaptive Method for the Cubic Nonlinear Schr\"{o}dinger Equation

TL;DR

This paper introduces an $hr$-adaptive numerical method combining mesh refinement and moving mesh strategies to accurately solve the nonlinear Schrödinger equation with steep gradients, demonstrating improved accuracy and convergence.

Contribution

The paper presents a novel $hr$-adaptive approach for solving the NLSE, effectively combining $r$-adaptive and $h$-adaptive methods for enhanced accuracy and error control.

Findings

Achieves excellent solution accuracy compared to other methods

Controls spatial error based on user-defined tolerance

Demonstrates second-order spatial convergence

Abstract

The nonlinear Schr\"{o}dinger equation (NLSE) is one of the most important equations in quantum mechanics, and appears in a wide range of applications including optical fibre communications, plasma physics and biomolecule dynamics. It is a notoriously difficult problem to solve numerically as solutions have very steep temporal and spatial gradients. Adaptive moving mesh methods ($r$-adaptive) attempt to optimise the accuracy obtained using a fixed number of nodes by moving them to regions of steep solution features. This approach on its own is however limited if the solution becomes more or less difficult to resolve over the period of interest. Mesh refinement methods ($h$-adaptive), where the mesh is locally coarsened or refined, is an alternative adaptive strategy which is popular for time-independent problems. In this paper, we consider the effectiveness of a combined method ($hr$-adaptive) to solve the NLSE in one space dimension. Simulations are presented indicating excellent solution accuracy compared to other moving mesh approaches. The method is also shown to control the spatial error based on the user's input error tolerance. Evidence is also presented indicating second-order spatial convergence using a novel monitor function to generate the adaptive moving mesh.

Tables

Table 1. Table 1 : Butcher array for SDIRK2 method.

$\gamma$	$\gamma$	$0$
$1$	$1-\gamma$	$\gamma$
	$1-\gamma$	$\gamma$

Table 2. Table 2 : Acronyms used to characterise simulation performance.

NHR	Number of $h$-refinements performed
NMAX	Largest number of nodes used
NMIN	Least number of nodes used
NSTP	Number of time steps used
JACS	Number of Jacobian computations throughout the simulation
	(when computing ${𝒌}_1$ and ${𝒌}_2$)
BS	Number of back solves needed when an $LU$
	factorisation is used in the quasi-Newton iterations.
ETF	Number of times the time step is halved due to fail in error test i.e.
	ERR $>$ ETOL or mesherr $>$ MESHTOL (§2.3)
CTF	Number of convergence test fails in Newton iteration i.e.
	${𝒌}_1$ or ${𝒌}_2$ fail to converge

Table 3. Table 3 : Performance of the h ​ r ℎ 𝑟 hr -algorithm for the propagation of a single soliton.

NHR	NMAX	NMIN	NSTP	BS	CTF	ETF	JACS
$0$	$78$	$78$	$434$	$2538$	$0$	$0$	$2538$

Table 4. Table 4 : Tolerance proportionality of h ​ r ℎ 𝑟 hr -algorithm for single soliton problem.

RTOL	$N^0$	${\Vert e\Vert }_{L_2}$
$1.5\times {10}^{-2}$	$67$	$2.95\times {10}^{-3}$
$3.75\times {10}^{-3}$	$175$	$8.87\times {10}^{-4}$
$9.37\times {10}^{-4}$	$352$	$2.43\times {10}^{-4}$

Table 5. Table 5 : Comparison of conserved quantities using h ​ r ℎ 𝑟 hr -adaptivity and a fixed uniform mesh.

Approach	$N$	$\Delta t$	$\left|Q_h^0-Q\right|$	$\left|{\overline{Q}}_h-Q\right|$	$\left|E_h^0-E\right|$	$\left|{\overline{E}}_h-E\right|$
$hr$-adaptivity	$78$	$\approx {10}^{-1}$	$7.9\times {10}^{-3}$	$7.0\times {10}^{-3}$	$1.4\times {10}^{-2}$	$2.8\times {10}^{-2}$
Uniform grid	$78$	$5\times {10}^{-3}$	$4.5\times {10}^{-2}$	$4.5\times {10}^{-2}$	$3.2\times {10}^{-1}$	$3.2\times {10}^{-1}$
MCN	$78$	$5\times {10}^{-3}$	$4.5\times {10}^{-2}$	$4.5\times {10}^{-2}$	$3.2\times {10}^{-1}$	$3.2\times {10}^{-1}$
Uniform grid	$400$	$5\times {10}^{-3}$	$5.3\times {10}^{-15}$	$4.6\times {10}^{-7}$	$2.1\times {10}^{-2}$	$2.1\times {10}^{-2}$

Table 6. Table 6 : L 2 subscript 𝐿 2 L_{2} error at T = 1 𝑇 1 T=1 using r 𝑟 r -method for single soliton problem.

	${\Vert e\Vert }_{L_2}$
$N$	$r$-method	Revilla  [34]	[36]	Uniform
$50$	$1.8\times {10}^{-3}$	$1.8\times {10}^{-2}$	$6.8\times {10}^{-2}$	$9.4\times {10}^{-2}$
$100$	$4.6\times {10}^{-4}$	$3.3\times {10}^{-3}$	$1.6\times {10}^{-2}$	$6.1\times {10}^{-2}$
$200$	$1.2\times {10}^{-4}$	$2.1\times {10}^{-3}$	$3.9\times {10}^{-3}$	$1.4\times {10}^{-2}$

Table 7. Table 7 : Performance of h ​ r ℎ 𝑟 hr -algorithm for the interaction of two solitons.

NHR	NMAX	NMIN	NSTP	JACS	BS	ETF	CTF
$7$	$197$	$134$	$909$	$4524$	$4524$	$0$	$0$

Table 8. Table 8 : Computational statistics for bound state of three solitons over 5 periods.

NHR	NMAX	NMIN	NSTP	JACS	BS	ETF	CTF
$8$	$332$	$100$	$856$	$5157$	$5157$	$7$	$0$

Equations

$$\left.\begin{array}[]{ll}i\psi_{t}+\psi_{xx}+q|\psi|^{2}\psi=0,&\quad-\infty<x<\infty,\>0<t\leq T,\\ \qquad\psi(x,0)=\psi_{0}(x),&\quad-\infty<x<\infty,\end{array}\right\}$$

$$i{\varphi }_t+{\varphi }_{xx}=0,$$

$$\varphi (x,t)=\exp [i(kx-W(k)t)],$$

$$Q$$

$$E$$

$$\left.\begin{array}[]{l}u_{t}+v_{xx}+q(u^{2}+v^{2})v=0,\\ v_{t}-u_{xx}-q(u^{2}+v^{2})u=0,\end{array}\right\}\>x\in(x_{l},x_{r}),\;t\in(0,T],$$

$$\left.\begin{array}[]{l}u(x,0)=\psi_{0R}(x),\quad v(x,0)=\psi_{0I}(x),\quad x\in(x_{l},x_{r}),\\ u(x_{l},t)=u(x_{r},t)=0,\\ v(x_{l},t)=v(x_{r},t)=0,\quad t\in(0,T],\end{array}\right\}$$

$$x=x(\xi ,t),\xi \in {\Omega }_c\equiv [0,1],t\in (0,T],$$

$${\xi }_i=\frac{i}{N},i=0,1,\dots ,N,$$

$$x_l=x_0<x_1(t)<\cdots <x_{N-1}(t)<x_N=x_r,$$

$$x_i(t)=x({\xi }_i,t),i=0,1,\dots ,N.$$

$$\dot{u}$$

$$\dot{v}$$

$${\dot{U}}_i$$

$${\dot{V}}_i$$

$$U_0=0,U_N=0,V_0=0\text{ and }{V}_N=0.$$

$$\frac{\partial x}{\partial t}=\frac{1}{\tau }{\left(M\frac{\partial x}{\partial \xi }\right)}^{-2}\frac{\partial }{\partial \xi }\left(M\frac{\partial x}{\partial \xi }\right),\xi \in {\Omega }_c,x(0,t)=x_l,x(1,t)=x_r,$$

$$M=\frac{1}{2}(M_u+M_v),$$

$$M_u={\phi }_u(t)+\sqrt{u_{xx}(x,t)},{\phi }_u(t)=\frac{1}{|{\Omega }_p|}{\int }_{{\Omega }_p}\sqrt{u_{xx}^2(x,t)}\mathrm{d}x$$

$$M_v={\phi }_v(t)+\sqrt{v_{xx}(x,t)},{\phi }_v(t)=\frac{1}{|{\Omega }_p|}{\int }_{{\Omega }_p}\sqrt{v_{xx}^2(x,t)}\mathrm{d}x.$$

$$g_{i+\frac{1}{2}}=\frac{U_{i+1}-U_i}{h_{i+1}},$$

$$w_i=\sqrt{2\left|\left(\frac{g_{i+\frac{1}{2}}-g_{i-\frac{1}{2}}}{h_{i+1}+h_i}\right)\right|}.$$

$${\phi }_u=\frac{1}{x_R-x_L}\sum _{i=0}^{N-1}\frac{h_{i+1}}{2}(w_{i+1}+w_i).$$

$$(M_u{)}_{i+\frac{1}{2}}={\phi }_u+\frac{1}{2}(w_{i+1}+w_i),i=0,1,\dots ,N-1.$$

$$M_{i+\frac{1}{2}}=\frac{1}{2}((M_u{)}_{i+\frac{1}{2}}+(M_v{)}_{i+\frac{1}{2}}).$$

$${\tilde{M}}_{i+\frac{1}{2}}=\frac{\sum _{k=i-p}^{i+p}{M}_{k+\frac{1}{2}}{\left(\frac{\gamma }{\gamma +1}\right)}^{|k-i|}}{\sum _{k=i-p}^{i+p}{\left(\frac{\gamma }{\gamma +1}\right)}^{|k-i|}},$$

$${\dot{x}}_i=\frac{4}{\tau }({\tilde{M}}_i(h_{i+1}+h_i){)}^{-2}({\tilde{M}}_{i+\frac{1}{2}}{h}_{i+1}-{\tilde{M}}_{i-\frac{1}{2}}{h}_i),$$

$${\tilde{M}}_i=\frac{{\tilde{M}}_{i-\frac{1}{2}}(x_{i+\frac{1}{2}}-x_i)+{\tilde{M}}_{i+\frac{1}{2}}(x_i-x_{i-\frac{1}{2}})}{x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}},$$

$$x_{i+\frac{1}{2}}=\frac{x_{i+1}+x_i}{2}.$$

$$\eta (t)={\left(\frac{1}{N}{\int }_{{\Omega }_p}Mdx\right)}^2,$$

$$\beta \text{ RTOL}\le \eta (t)\le \alpha \text{ RTOL}.$$

$${\tilde{N}}^{n+1}=N^n\times \min \left(\text{maxfac},\max \left[\text{minfac},\kappa {\left(\frac{\eta (t)}{RTOL}\right)}^{\frac{1}{2}}\right]\right),$$

Full text

An $hr$-Adaptive Method for the Cubic Nonlinear Schrödinger Equation333© 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/

J.A. Mackenzie111Corresponding author: Department of Mathematics and Statistics, University of Strathclyde, Livingstone Tower, 26 Richmond Street, Glasgow G1 1XH, Scotland

  

W.R. Mekwi222School of Computing, Engineering and Physical Sciences, University of the West of Scotland, Paisley Campus, Paisley PA1 2BE, Scotland

Abstract

The nonlinear Schrödinger equation (NLSE) is one of the most important equations in quantum mechanics, and appears in a wide range of applications including optical fibre communications, plasma physics and biomolecule dynamics. It is a notoriously difficult problem to solve numerically as solutions have very steep temporal and spatial gradients. Adaptive moving mesh methods ($r$-adaptive) attempt to optimise the accuracy obtained using a fixed number of nodes by moving them to regions of steep solution features. This approach on its own is however limited if the solution becomes more or less difficult to resolve over the period of interest. Mesh refinement methods ($h$-adaptive), where the mesh is locally coarsened or refined, is an alternative adaptive strategy which is popular for time-independent problems. In this paper, we consider the effectiveness of a combined method ($hr$-adaptive) to solve the NLSE in one space dimension. Simulations are presented indicating excellent solution accuracy compared to other moving mesh approaches. The method is also shown to control the spatial error based on the user’s input error tolerance. Evidence is also presented indicating second-order spatial convergence using a novel monitor function to generate the adaptive moving mesh.

Keywords

Adaptivity, moving mesh method, $hr$-adaptivity, cubic nonlinear Schrödinger equation

1 Introduction

The nonlinear Schrödinger equation (NLSE) is an important model of mathematical physics and has numerous applications in the physical sciences such as in optical fibre communications, plasma physics, biological and atomic physics, biomolecule dynamics, hydrodynamics etc. (see e.g. [40, 44] and the references therein). Its properties are also of great mathematical interest and solving it numerically has always proved to be challenging due to the solitons that appear (and disappear) in the solution.

When written in a reference frame moving at the group velocity of the carrying wave, the cubic NLSE takes the form

[TABLE]

where $\psi:\mathbb{R}\times\mathbb{R}^{+}\rightarrow\mathbb{C}$, $i^{2}=-1$ and $q$ is a given real constant. This work only considers the focusing case for which $q$ is positive.

The analytical properties of (1.1) are well known (see [18, 38]). In particular, the linear Schrödinger equation

[TABLE]

provides a model for the propagation of dispersive waves given by

[TABLE]

where $W(k)=k^{2}$. The phase speed is defined by $W(k)/k$ and clearly depends on $k$. The wave is thus dispersive [18, 38]. It can be shown that the solutions of (1.2) have an amplitude which decays like $t^{-1/2}$ for $t,x\to\infty$ with $x/t$ fixed (see [47, §11.3]).

The cubic term in (1.1) opposes dispersion and hence makes it possible for the NLSE to possess solutions where the competing forces of nonlinearity and dispersion balance each other exactly. It has been shown that, using the initial condition $\psi(x,0)=$ sech $x$, this is achieved when $q=2$ [18]. This balance happens for $q=8,18,$ or in general when $q=2N^{2}$ for integer $N$. These states correspond to bound states of $N$ solitons or solitary waves. Solitons are formed when a certain balance between nonlinearity and dispersion is reached.

The pure initial-value problem (1.1) has been shown to possess an infinite number of conservation laws (or so-called invariants of motion) [35]. The most common of these are the charge $Q$, and the energy $E$ [14], given by

[TABLE]

These conservation laws play an important role in the analysis and dynamics of the NLSE and many numerical schemes have been built [13, 18, 37, 40, 3] that attempt to conserve either or both of these quantities during simulations. Striving for methods to conserve energy is very important as conservation of energy implies the $L^{2}$-boundedness of the solution, thus preventing blow-up of the computed solution [38]. However, although lack of exact (energy) conservation may lead to nonlinear blow-up, time integrators which are energy-conserving may still not perform desirably for this problem [38].

Adaptive mesh methods for the solution of partial differential equations (PDEs) generally attempt to optimise the number and/or placement of mesh cells. These techniques have been proposed for a few decades now and can broadly be categorised as $h$-, $p$- and $r$- adaptive methods. For $p$-adaptive strategies, the order of polynomials which represent the solution locally, is varied in order to achieve better solution accuracy. Mesh refinement or $h$-adaptive methods are the most widely used adaptive methods. These methods start with a fixed number of mesh elements and then, based on some local measure of problem difficulty or a posteriori error estimate, elements are added or removed as necessary during the time integration of the problem. This approach is particularly useful when the solution is required to attain a specified accuracy. It is common to combine $h$- and $p$-adaptive strategies and this approach is well understood and is widely used. However, the method is complex and the underlying a posteriori estimates used to drive them can be rather difficult to obtain for strongly nonlinear problems [8, 30].

In contrast, moving mesh or $r$-adaptive methods generally aim to use a given number of mesh cells efficiently by prioritising areas of large solution variation and placing fewer mesh cells where the solution is smooth. This approach has been used with much success in a wide range of applications (see e.g. [17, 6, 10, 28, 8] and references therein). When compared with a fixed uniform mesh approach, $r$-adaptive methods are able to maintain similar or achieve better levels of accuracy using fewer mesh points [17]. The main shortcoming of an $r$-adaptive approach is that the solution accuracy is limited due to the fixed number of nodes.

Moving mesh methods for time-dependent problems may be further divided into static and dynamic methods. With static methods the PDE is solved between time steps assuming that the mesh is fixed. The mesh is then modified and the physical solution is transferred between meshes using interpolation or an appropriate projection technique. For dynamic methods, the number and connectivity of mesh elements is kept fixed for most of the simulation and the mesh is moved at every time step. A mesh equation is usually used to compute node speeds in order to move the mesh. This has the advantage that one is not required to transfer the solution between meshes because the PDE is reformulated to take into account the fact that the nodes are moving [49, 42, 41, 43, 26, 15]. Furthermore, these methods are thought to do a good job of reducing “dispersive errors”, a property that is useful for this problem. The experiments presented here employ a dynamic technique, using a moving mesh PDE (MMPDE) to drive the movement of the mesh.

Even though numerous studies of the NLSE have been undertaken [36, 14, 7, 11, 48, 46, 3] these methods usually focus on the conservative properties of the NLSE or only apply $r$-refinement in order to track or resolve the emerging solitons. Some authors have employed a combined $hr$-refinement approach to solve problems from applications including combustion and fluid mechanics [1, 27, 32, 19], elastostatics and fracture mechanics [2], heat transfer problems [25] and ocean modelling [33]. In previous work [28] we proposed an $hr$-algorithm and applied it to a number of one-dimensional time-dependent PDEs. To our knowledge, few authors have employed an $hr$-refinement strategy for the NLSE. In combining these methods, we aim to efficiently track the travelling soliton(s) (mesh movement) while controlling solution accuracy appropriately using mesh refinement.

The layout of this paper is as follows: section 2 describes the components of the $hr$-adaptive method for the NLSE and presents the algorithm. Numerical experiments are reported and discussed in section 3. Conclusions and future work are given in section 4.

2 Numerical method

We employ an adaptive mesh technique for the NLSE which involves adjusting both the number and position of nodes of the mesh, a so-called $hr$-adaptive approach. In addition, we also vary the time-step size during the computation. This section gives details of how the algorithm has been implemented.

2.1 Spatial discretisation

The first step in approximating the solution is re-defining the pure initial value problem (1.1) as an initial-boundary value problem in $x_{l}\leq x\leq x_{r}$, since over the time interval $[0,T]$ under consideration, the solutions of (1.1) are assumed to be negligibly small outside $[x_{l},x_{r}]$. At these boundaries, it is customary to pose homogeneous Dirichlet or Neumann boundary conditions [18, 34, 36]; our experiments use the former.

The complex function $\psi$ is first decomposed into its real and imaginary parts $u$ and $v$ respectively, resulting in the coupled system of PDEs

[TABLE]

[TABLE]

where $\psi_{0R}$ and $\psi_{0I}$ are, respectively, the real and imaginary parts of $\psi_{0}$.

To include node movement, (2.5) is rewritten with respect to a moving reference frame and the problem is recast in terms of the independent variables $\xi$ and $t$, using the transformation

[TABLE]

from computational space $\Omega_{c}\times(0,T]$ to physical space $\Omega_{p}\times(0,T]$. A uniform mesh covering $\Omega_{c}$ is given by

[TABLE]

and the corresponding (nonuniform) mesh on $\Omega_{p}$ is

[TABLE]

where

[TABLE]

In what follows, we define \mbox{\boldmathx}(t)=\left\{x_{i}(t)\right\}_{i=0}^{N}.

To incorporate mesh movement, it is convenient to express the Eulerian time derivatives in (2.5) in terms of the moving reference frame. The resulting equations are given by

[TABLE]

where $\dot{u}$, $\dot{v}$ and $\dot{x}$ denote derivatives with respect to $t$ where $\xi$ is held constant.

After applying central differencing for the spatial derivatives, the semi-discrete system on a moving grid is given by

[TABLE]

for $i=1,\ldots,N$, where $h_{i}=x_{i}-x_{i-1}$, and $U_{i}(t)$, $V_{i}(t)$ are the approximations to $u(x_{i},t)$ and $v(x_{i},t)$. Further, we define $\mathbf{U}=(U_{0},U_{1},\ldots,U_{N-1},U_{N})^{T}$ and $\mathbf{V}=(V_{0},V_{1},\ldots,V_{N-1},V_{N})^{T}$.

Application of the zero Dirichlet boundary conditions yields

[TABLE]

The discretisation of $\dot{x}$ is described in section 2.3.

2.2 Mesh adaptation

The mesh adaptation process is a twofold process involving node movement and mesh refinement.

2.2.1 Mesh movement

The aim of node movement is to place more nodes where the solution has features which are difficult to resolve. We achieve this using a monitor function. This is typically a function of the solution gradient or curvature, and its choice is crucial to the success of the mesh movement process [20]. The placing of nodes is based on an equidistribution principle, where nodes are located to ensure that some measure of error is equally distributed over the mesh [12].

The mesh velocity $\dot{x}$ is computed using a mesh generating equation. We use a moving mesh PDE (MMPDE) to achieve this. The impact of MMPDEs on the moving mesh method has been studied in some detail by Huang and his collaborators ([21, 9, 23, 24]). Our experiments here use a one-dimensional version of one of the two-dimensional MMPDEs proposed by Huang and Russell [23] given by:

[TABLE]

where $\tau$ is a temporal smoothing parameter and $M({\nu}(x,t))$ is a positive monitor function. The initial condition $x(\xi,0)$ is obtained by equidistribution of the a monitor function based on the initial condition at $t=0$.

Recognising that the NLSE is a system of equations for the real and imaginary components $u$ and $v$, and motivated by the work in [4, 45], our experiments utilise the second derivative monitor function

[TABLE]

where

[TABLE]

and

[TABLE]

The role of $\varphi_{u,v}(t)$ is to ensure that some mesh points are placed in regions that do not have steep solution features e.g. outside boundary and interior layers. Some analysis on the effect of using monitor functions in general has been carried out in [4, 5, 22, 24].

Discretisation of (2.15) requires the evaluation of the monitor function (2.16). If we let

[TABLE]

then an approximation $w_{i}$ to $\sqrt{|u_{xx}|}$ is given by

[TABLE]

The floor, $\varphi_{u}$, on the monitor function, $M_{u}$, is given by a quadrature approximation

[TABLE]

The discrete approximation of the monitor function (2.17) is then given by

[TABLE]

The approximation of $(M_{v})_{i+{\textstyle\frac{1}{2}}}$ is calculated similarly and finally

[TABLE]

To improve the robustness of the moving mesh method a smoothed monitor function is obtained as in [31] by setting

[TABLE]

with $\gamma=2$ and $p=3$ in all our experiments. In addition, the summations only use terms that are well defined $(0\leq k\leq N-1)$. Using second-order central differences for the spatial derivatives in (2.15), we obtain a semi-discrete system of moving mesh equations defined by

[TABLE]

for $i=1,2,\ldots,N-1$, with $x_{0}=x_{l}$ and $x_{N}=x_{r}$. The term $\tilde{M}_{i}$ in (2.24) is given by

[TABLE]

where

[TABLE]

2.2.2 Mesh refinement

In choosing the number of nodes to use at each time step, the objective is to keep some measure of the spatial solution error within predetermined bounds. To implement this strategy, we define an error monitor, which broadly measures global difficulty of the problem under consideration. Let us define

[TABLE]

where $M$ is given by (2.16). If RTOL is a user-prescribed tolerance and $\alpha$ and $\beta$ are given such that $\alpha>1$ and $0<\beta<1$, then we would like to ensure that

[TABLE]

Note that $\eta(t)$ is simply a cheap, heuristic means that enables us to determine when to change the number of nodes. Also, the number of remeshings performed during a calculation is inversely proportional to the difference between $\alpha$ and $\beta$.

To calculate the number of nodes required at time $t^{n+1}$ we first evaluate

[TABLE]

where maxfac, minfac and $\kappa$ are user-defined parameters. The exponent of $\frac{1}{2}$ in the formula is related to the assumption that the scheme is second-order accurate in space. In our experiments we have set these values to be: maxfac $=2.0$, minfac $=1.2$ or $0.3$ for mesh enrichment or mesh coarsening respectively, and $\kappa\geq 1$. The new number of mesh nodes is then

[TABLE]

where $\lfloor\cdot\rfloor$ denotes “integer part of”. The prediction of the number of nodes (2.28) has been motivated using the same idea as has been used to control the time step in (2.36).

2.3 Temporal adaptation

Integration from $t=t^{n}$ to $t=t^{n+1}$ is accomplished following the approach of Beckett *et al. *[5]. The semi-discretised systems (2.13), (2.14) and (2.24) for the solution and the mesh are clearly coupled. For efficiency, we solve them using a decoupled iterative procedure which we outline here.

We discretise (2.24) in time using the first-order backward Euler method. This equation is the governing equation for \mbox{\boldmathx}\equiv\{x_{i}\}_{i=0}^{N}, for which $\tilde{M}_{i\pm\frac{1}{2}}$ and $(\tilde{M}_{i}(h_{i+1}+h_{i}))^{-2}$ are known for $i=1,2,\ldots,N-1$. We use an iterative approach where, initially, $\tilde{M}_{i\pm\frac{1}{2}}$ and $(\tilde{M}_{i}(h_{i+1}+h_{i}))^{-2}$ are evaluated at $t=t^{n}$. Subsequent iterations use values from the preceding cycle. The linear system which arises at each iteration is solved directly. The final mesh that is used is obtained by combining this mesh and the mesh from the previous iteration using under-relaxation so that the $(\nu+1)$th iterate is given by

[TABLE]

where \mbox{\boldmathx}^{[n+1,\nu+1]}_{*} is the mesh obtained from solving (2.24) and $\omega=0.8$ in all our computations.

When integrating from $t=t^{n}$ to $t=t^{n+1}$, the system (2.13) and (2.14) is considered a system of equations for the approximation of $\mathbf{U}$ and $\mathbf{V}$, where the node locations are available at $t=t^{n}$ and $t=t^{n+1}$. The vector \dot{\mbox{\boldmathx}} is replaced by (\mbox{\boldmathx}^{n+1}-\mbox{\boldmathx}^{n})/\Delta t^{n} and $x$ is evaluated in $[t^{n},t^{n+1}]$ using the linear interpolant

[TABLE]

where \mbox{\boldmathx}^{n} is the approximation to $x$ at $t=t^{n}$. We also note that, while some authors (see e.g. [11]) use an upwinding approach for $\dot{x}$ in (2.11) and (2.12), such treatment is not used here.

The equations (2.13) and (2.14) can be written as a coupled set of ODEs

[TABLE]

where \mbox{\boldmathf}:\mathbb{R}\times\mathbb{R}^{2(N+1)}\to\mathbb{R}^{2(N+1)} and

[TABLE]

The solution of (2.31) at $t=t^{n+1}$ is computed using a second-order singly diagonally implicit Runge-Kutta method (SDIRK2) [5] which possesses excellent stability properties. The Butcher array is shown in Table 1, where $\gamma=(2-\sqrt{2})/2$.

Integration from $t=t^{n}$ to $t=t^{n+1}=t^{n}+\Delta t^{n}$ is given by the solution of

[TABLE]

where \mbox{\boldmathw}^{n} denotes the value of $w$ at $t=t^{n}$.

The values of \mbox{\boldmathk}_{1} and \mbox{\boldmathk}_{2} are obtained using a Newton iteration, terminating when convergence is achieve in the $l_{\infty}$ norm to within a user-prescribed tolerance, KTOL.

Our algorithm also employs adaptive time stepping. To implement this, we use an error indicator based on the computed solution and a first order solution that can be cheaply obtained as a by product of the SDIRK computation used. Indeed, if \hat{\mbox{\boldmathw}}^{n+1} is a first order approximation to $w$ at time $t=t^{n+1}$, we obtain \hat{\mbox{\boldmathw}}^{n+1} using

[TABLE]

and the error indicator ERR is calculated using the mesh-dependent $L_{2}$ error measure

[TABLE]

where

[TABLE]

In the event that ${\rm ERR}>{\rm ETOL}$, where ETOL is a user-supplied tolerance, the time step is rejected and repeated with a smaller time step. On the other hand if ${\rm ERR}<{\rm ETOL}$, then the predicted suitable time step, based on the solution, for the next time step is given by the formula (see [16])

[TABLE]

where typically $1.5\leq$ maxfac $\leq 3$, $\eta\sim 0.6,$ and minfac$=0.1$.

A similar computation based on the mesh is carried out to predict the time step. Here, we use an indicator of grid accuracy, mesherr, given by

[TABLE]

where \mbox{\boldmathx}^{[n+1,\upsilon]} is the final mesh at $t=t^{n+1}$ obtained from (2.30), and \mbox{\boldmathx}^{[n+1,\upsilon-1]} denotes the mesh at the previous pass in the iterative procedure. If MESHTOL is a user-supplied mesh tolerance and ${\rm mesherr}>{\rm MESHTOL}$, then the time step is rejected and repeated with a smaller time step. If ${\rm mesherr}<{\rm MESHTOL}$, then the predicted mesh time step is calculated using:

[TABLE]

where MESHBAL is a user-chosen parameter and we require that MESHBAL $<$ MESHTOL. We finally take the minimum of the solution time step (2.36) and the predicted mesh time step (2.38) so that

[TABLE]

2.4 The complete $hr$-algorithm

The algorithm which we have used to solve the NLSE is shown in Algorithm 1 below, where $\psi^{n}$ represents the numerical solution at time $t^{n}$. Note that, in the absence of steps 1-1 and the computation of $N^{0}$ in step 1, the method performs $r$-adaptation with the initial number of nodes provided.

The initial mesh $N^{0}$ is computed automatically so that $\eta(0)$ satisfies (2.27). This is achieved by equidistribution of the initial condition using the de Boor algorithm [12] which iterates until the difference between successive meshes is less than some prescribed tolerance, GTOL.

In the event that $h$-refinement is performed, the new grid \displaystyle\mbox{\boldmathx}^{n+1}_{N^{n+1}} is obtained via equidistribution of the monitor function defined on the mesh \displaystyle\mbox{\boldmathx}^{n+1}_{N^{n}}. Using cubic interpolation, the solution is then transferred from \displaystyle\mbox{\boldmathx}^{n+1}_{N^{n}} onto the newly computed grid, \mbox{\boldmathx}_{N^{n+1}}^{n+1}.

3 Numerical experiments

We now consider the performance of the $hr$-algorithm on three increasingly difficult test problems.

3.1 Propagation of a single soliton

If the initial condition is given by

[TABLE]

then the solution is given by the single soliton

[TABLE]

where $c$ is the speed of the soliton and $a$ is a real parameter which determines the amplitude. The modulus $\rho=|\psi(x,t)|$ is characterised by a single soliton of amplitude $\sqrt{\frac{2a}{q}}$, which is located initially at $x=x_{0}$, and travels to the right with speed $c$.

The computations shown here used the values $c=1$, $a=1$, $q=1$ and $x_{0}=0$. For the spatial domain we set $x_{l}=-30$ and $x_{r}=70$ and integrated up to time $T=30$. The simulations were performed using the following parameters and tolerances: $\alpha=1.4$, $\beta=0.8$, RTOL$=1.5\times 10^{-2}$, ETOL$=5\times 10^{-3}$ and $\tau=10^{-3}$. Based on equidistribution of the initial condition, the algorithm determined that $N^{0}=78$ grid nodes were needed to resolve the initial condition to the required tolerance. Plots of the numerical and exact solutions at $t=0,10,15,20,25$ and $30$ are shown in Figure 1. We can see that the algorithm does an excellent job of resolving the travelling soliton.

The grid point trajectories shown in Figure 2 (showing alternate grid points) indicate that the $r$-adaptive component of the adaptive algorithm performs well tracking the movement of the front throughout the time integration period.

The performance of the $hr$-algorithm is shown in Table 3, where the acronyms used are described in Table 2. We can see that no $h$-refinement steps were needed throughout the integration process as the travelling wave was well resolved using the $r$-refinement component of the adaptive strategy.

To compare our results with previous work, we computed the $L_{2}$ error at time $t^{n}$ as follows:

[TABLE]

where $e_{i}^{n}=\tilde{\rho}(x_{i},t^{n})-\rho(x_{i},t^{n})$ and $\tilde{\rho}$ is the exact solution. The first column of Figure 3 shows the evolution of the spatial error indicator $\eta$, the time step history, and the $L_{2}$ norm of the error when a value ${\rm RTOL}=6\times 10^{-2}$ is used. One can see that the algorithm ensures that $\beta\>\rm{RTOL}\leq\eta(t)\leq\alpha\>\rm{RTOL}$, as expected. In addition, the time step is chosen to be reasonably large and is relatively constant throughout the simulation. One can also see that the error remains approximately constant throughout the simulation.

The second and third columns of Figure 3 show how $\eta$, the time step, and the $L_{2}$ norm of the solution error at $t=T$, vary as RTOL is decreased by a factor of four. We can see from Table 4 that $N^{0}$ is approximately doubled as ${\rm RTOL}$ is quartered and that this results in a reduction in the $L_{2}$-norm of the error and the rate of decrease is approaching a factor of four. This behaviour is consistent with the expectation that the overall algorithm is spatially second-order accurate.

For this problem the exact values of the conserved quantities are $Q=4$ and $E=-1/3$. The conserved quantities $Q$ and $E$ are approximated at time $t^{n}$ by

[TABLE]

and

[TABLE]

We define the mean values $\bar{Q}_{h}$ and $\bar{E}_{h}$ obtained from the numerical solution over the time integration period as

[TABLE]

where $NT$ is the total number of time steps taken. Table 5 shows the initial values $Q_{h}^{0}$ and $E_{h}^{0}$, and the mean values $\bar{Q}_{h}$ and $\bar{E}_{h}$ using the $hr$-adaptive algorithm. The equivalent quantities are also shown when two uniform meshes are used: one using the SDIRK2 time integration scheme and the other using the energy-conserving modified Crank-Nicolson scheme (MCN) [13]. We have also shown the error in $Q$ and $E$ when a finer uniform mesh is used. We can see that the initial approximations $Q_{h}^{0}$ and $E_{h}^{0}$ are much more accurate using the adaptive mesh compared to the schemes using a uniform mesh with the same number of mesh points; this is especially true for the approximation of $E$ which relies on the accurate approximation of solution derivatives. We can see that there is no noticeable difference in the approximation to $Q$ and $E$ using the conservative scheme; the error being equal to its initial error as expected. We see that the value of $Q^{n}_{h}$ is almost constant throughout the computation using the $hr$-algorithm. The conservation of $E$ is less well maintained but $\bar{E}_{h}$ is far more accurately predicted with the $hr$-algorithm compared to a uniform mesh with the same number of mesh nodes. It requires at least $N=400$ nodes and a time step of size $\Delta t=5\times 10^{-3}$ to achieve a comparable level of solution accuracy (in the $L_{2}$ norm) to an adaptive grid method, which utilises $N=78$ nodes and $\Delta t\approx 10^{-1}$.

In Figure 4(a), we plot the evolution of $L_{2}$-norm of the error. For comparison we have also included the error calculated using the MCN scheme on a uniform mesh with $N=78$. It is clear that the error using the moving mesh method is considerably smaller than the error obtained using the MCN scheme and this is despite the fact that the energy is not conserved with the moving mesh method as can be seen in Figure 4(b). This example therefore suggests that it is important to accurately resolve spatial gradients and that it is not sufficient to use a numerical method which conserves an approximation to the energy if the spatial mesh is not dense enough.

Table 6 shows the $L_{2}$ error at $T=1$ for various values of $N$ using the moving mesh method when mesh enrichment and coarsening are switched off. These results are comparable with those obtained by previous authors [34, 36, 39]. The results obtained by Revilla [34] for a similar experiment are also displayed. Whereas the $hr$-algorithm we use utilises interpolation only when the mesh is coarsened/enriched, the method used by Revilla obtains the numerical solution on the new grid via an interpolation based method. It is possible that frequent interpolation is a source of spatial error in the Revilla approach which could account for the discrepancy between the results. The approach in [36] uses the arc-length monitor function based on first-derivatives of the numerical solution and interpolation to transfer the solution between meshes. It is clear that our method outperforms both of these moving mesh approaches and is far more accurate than the use of a uniform mesh.

3.2 Interaction of two solitons

We next consider the initial condition

[TABLE]

where

[TABLE]

In the simulations below we have chosen the parameters $a_{1}=0.2$, $c_{1}=1$, $x_{01}=0$ and $a_{2}=0.5$, $c_{2}=-0.2$, $x_{02}=25$ for $\psi_{1}$ and $\psi_{2}$, respectively with $q=1$. This initial condition consists of two solitons of different amplitudes, located initially at $x_{01}=0$ and $x_{02}=25$ which then move to the right and left with speeds $c_{1}$ and $c_{2}$, respectively. The solitons interact as if they were particle-like entities, i.e., they exhibit elastic collisions from which they emerge with the same shape [44].

We integrated the problem up to the final time $T=45$ and set the artificial boundaries $x_{l}=-20$ and $x_{r}=80$. For this computation we used ETOL$=5\times 10^{-4}$, $\tau=10^{-2}$, RTOL$=1\times 10^{-2}$, $\alpha=1.2$ and $\beta=0.8$. The solutions at $t=0,20,30,45$ are shown in Figure 5. The reference solution was obtained using a fixed uniform grid with $2000$ nodes. To plotting accuracy, the $hr$-algorithm has done an excellent job of capturing the interacting solitons. Details of the performance of the $hr$-algorithm are given in Table 7.

In Figure 6(a) we can see from the grid point trajectories that between $t=15$ and $t=25$ the two solitons meet and cross. It is within this time interval that refinements/derefinements of the grid are performed as can be seen from the grid point trajectories and evolution of $N$ in Figure 6(b). Note that the behaviour of $N$ in Figure 6 is what one would anticipate. The amplitudes of the solitons at $t=0$ and $t=45$ are about the same and so we would expect roughly the same number of nodes to be used. We can see that the algorithm utilised a very similar value for $N$ at these times. Furthermore, more nodes were introduced just as the problem got more challenging and nodes were removed when they were not needed.

It is clear from Figure 7(a) that the objective of keeping $\eta$ within the given bounds has been successfully achieved. Also, we see from Figure 7(b) that the algorithm has chosen the time step to be fairly constant throughout the simulation - it reduces appropriately when the solitons interact and recovers afterwards. A global plot of the solution profiles is given in Figure 8 as well as a plot of the error in the energy. We can see that there is a slight change in the energy as the solitons interact but this error does not increase dramatically through the simulation and is not far off its initial value.

For comparison in Figure 9 we show the computed solution obtained using the MCN scheme with a time step $\Delta t=10^{-3}$ and a uniform mesh using the average number of points used with the adaptive $hr$ scheme. For this example this resulted in a value of $N=177$. We can see clearly that the computed solutions are poor in comparison with the $hr$-adaptive method. The inaccuracy again is due to the lack of spatial resolution of the rapidly changing behaviour of the interacting solitions and this is not compensated for by the discrete conservation of energy.

3.3 Bound state of three solitons

The bound state of multiple solitons is a class of problems with initial condition

[TABLE]

and $q=2N_{s}^{2}$ in (1.1), where $N_{s}$ is a positive integer. For this case, $\psi(x,t)$ is a bound state of $N_{s}$ solitons [29].

As in [18, 34, 39] we study the bound state of three solitons, i.e., $q=18$, and use artificial boundaries located at $x_{l}=-20$ and $x_{r}=20$. The solution of this problem is periodic in time with the period being approximately $T=0.8$. The solution develops extremely large spatial and temporal gradients thus posing a stringent test to any numerical scheme.

For the purpose of comparison with previous work, we integrated over five periods and took $T=4$. For this simulation, we used the values RTOL$=10^{-3}$, $\alpha=3$, $\beta=0.4$ and ETOL$=5\times 10^{-3}$. For this problem we set $\varphi_{u,v}=10^{-3}$ throughout the computation. Note that the choice (2.17) and (2.18) places approximately half of the available nodes outside the region where we have high solution activity. For this problem, this leads to insufficient resolution in the soliton region. The value of $\varphi_{u,v}=10^{-3}$ is much smaller than that predicted by the algorithm thus placing more nodes in the centre of the domain.

The solutions at the times $t=0.2,0.4,0.6$ and $0.8$ are shown in Figure 10. The reference solution was obtained using a fixed uniform grid with $N=2000$. The mesh trajectories and the evolution of $\eta$ are displayed in Figure 11. The trajectories only cover the region $-5\leq x\leq 5$ where nodes experience significantly more activity. As for the previous cases, the algorithm tracks accurately the solitons as they appear and disappear during the integration. One can also see that the monitor captures regions of high curvature as expected. The error indicator $\eta$ is also kept within the required bounds. Computational statistics are shown in Table 8.

The evolution of $N$ and the time step history are shown in Figure 12. As expected, the algorithm adds nodes when they are needed (in this case when there are two solitons) and removes them when they are not needed. Furthermore, the time step is kept relatively large over the integration, another benefit of using an adaptive approach.

A global view of the solution trajectories is shown in Figure 13 as well as the evolution of the energy error. We can see that the solution is well approximated over the first few periods but that eventually the solution accuracy deteriorates. The energy error also increases more rapidly for this much more demanding test case as can be seen also in Figure 13. However, the computed solutions are far more accurate to those obtained using the MCN scheme on a uniform mesh as seen in Figure 14, where again we have used the average number of points used with the $hr$ method and a time step $\Delta t=10^{-4}$. As with the two previous examples, it’s clear that a lack of spatial resolution results in poor solution accuracy even though the scheme preserves the inital approximation to the energy.

4 Conclusions

We have developed an $hr$-adaptive algorithm which combines adaptive mesh movement and mesh enrichment for the cubic nonlinear Schödinger equation. Numerical experiments demonstrate that the method works well for problems with large spatial and temporal variations in comparison with other moving mesh methods as well as energy-conserving uniform grid methods with the equivalent number of grid nodes.

There is scope to improve the algorithm further and we hope to do this by considering more closely the effect and the choice of MMPDEs and monitor functions. It would also be useful to develop conservative versions of the moving mesh method which could potentially improve solution accuracy over longer time integration periods. For simplicity, in this work we have used linear interpolation to transfer solutions between meshes when $h$-adaptation is required. In future we will also consider the use of $L_{2}$ projection for transferring solutions between meshes. While we have focussed in this paper on the application of method to the cubic NLSE, it also has applicability to other non-linear dispersive wave equations. Another area of development would be to investigate the merits of using an $hr$-adaptive scheme in a multidimensional context.