Approximation of rogue waves using Malmquist-Takenaka functions

Abstract

Rogue waves are fascinating large amplitude coherent structures that abruptly appear and then disappear soon after. In certain partial differential equations these waves are modeled by rational solutions. In this work we discuss approximating rogue wave solutions in a basis of orthogonal functions known as the Malmquist-Takenaka (MT) functions. This family of rational functions can be directly mapped to a modified Fourier series, allowing the fast Fourier transform computation of the spectral MT coefficients. Spectral differentiation matrices are derived. The approximation of the various rogue wave solutions in the nonlinear Schrödinger (NLS) equation is explored. The unstable nature of the NLS equation on a constant background and its effect on destabilizing and generating rogue waves is studied. Perturbing the constant solution with certain localized functions is found to generate rogue wave-like structures.

[label1]organization=University of Colorado, Colorado Springs,addressline=1420 Austin Bluffs Parkway, city=Colorado Springs, postcode=80918, state=CO, country=USA

Approximation of rational solutions to partial differential equations in terms of a family of orthogonal rational functions called the Malmquist-Takenaka functions. For certain rational functions, these approximations converge spectrally (exponentially) fast. Spectral differentiation matrices are derived.

Using this technique, different models for rogue waves in the nonlinear Schrödinger equation are studied. The well-known Peregrine soliton is observed to be exponentially unstable to localized perturbations. A perturbation of the constant background generates rogue wave-type structures which appear to be rational in structure.

1 Introduction

Rogue, or freak, waves are brief, but intense wave formations that have the capability to inflict serious harm with little notice. A commonly accepted description is that of a coherent structure that is localized in both space and time. Since their initial detection on the Draupner platform in 1995 [1], several observations have been reported around the world [2, 3, 4]. As a result, the understanding and detection of rogue waves is important for the safety of maritime vessels and coastal infrastructure.

A canonical model for nonlinear dispersive waves is the nonlinear Schrödinger (NLS) partial differential equation (PDE). The NLS equation has been derived to describe a wide variety of physical systems, such as water waves [5], fiber optics [6], and plasmas [7], to name a few. The derivation of the model typically relies on a slowly-varying enveloped approximation, where a monochromatic linear wave is modulated by a slowly-varying amplitude function that is governed by an NLS-type equation [8]. While these models typically assume weak nonlinearity, the NLS equation often captures leading-order effects and is a reasonable place to begin a discussion on rogue waves in deep water.

Within the NLS equation, there are two well-known paradigms for describing the physical generation of rogue waves [2]. The first is rational solutions [9, 10] that are localized in space and time; the most well-known is the Peregrine soliton [11]. The second is modulational instability (MI), where perturbations with a small sideband wavenumber are linearly unstable [12, 13]. Recently, several works have highlighted the unstable nature of the Peregrine soliton [14, 15, 16], throwing into question how this state could be realized in the tempestuous ocean. The MI route is worth pursuing, but it should be performed on the real line, and not on a periodic domain, to truly capture it’s nature in large bodies of water. The numerical simulation of these two rogue wave models is the motivation of this paper.

The industry standard for numerically approximating wave-related (that is, periodic) phenomena is spectral Fourier methods. For relatively smooth solutions, Fourier spectral methods are an attractive option due to their rapid convergence rates [17, 18]. Fourier approximations can also be useful when approximating exponentially decaying solutions on the real line. In the latter case, one typically truncates the infinite line to a large but finite domain, on which an exponentially localized solution can be treated as a periodic function [19, 20]. However, this approach requires the function reach its (constant) boundary conditions well within the computation window. As such, secant hyperbolic-type soliton solutions are an ideal candidate for this approach.

However, the Peregrine soliton, a rogue wave model, is a rational function that does not decay exponentially fast. The Fourier method described above struggles with rational functions due to their (slow) algebraic decay rate. On most feasible computational domains, these rational solutions are nowhere near zero, or machine precision, in practice. Hence, the periodic extension of these functions on finite domains are not sufficiently smooth and so, in general, their Fourier coefficients decay algebraically (slow).

In order to better approximate rational functions several ideas have been explored. One idea is to use the Hermite functions [21, 22] as a basis. This is a natural approach since the Hermite functions form an orthogonal basis on the real line. However, there is no known Hermite fast transform and the truncated series only converges spectrally fast if the function decays exponentially fast at $\pm\infty$ [18]. Another approach is that of coordinate mapping [23, 24]. The idea here is to map the (infinite) real line to a finite interval by some coordinate transformation. One such mapping is the so-called algebraic type, e.g. $x=\frac{y}{\sqrt{1-y^{2}}}$ for $y\in(-1,1)$ [25]. A method that uses said mapping and approximates rational solutions of the nonlinear Schrödinger equation is the Chebyshev pseudo-spectral method given in [26]. Often scaling parameters are included in the transformation and chosen to optimize the accuracy of the method.

In this work we explore the idea of using a rational basis to approximate rational functions. A family of orthogonal rational functions was considered in [27] that consists of the algebraic mapping mentioned above substituted into the argument of the Chebyshev polynomials, i.e. $\psi_{n}(x)=T_{n}\left(\frac{x}{\sqrt{1+x^{2}}}\right),n=0,1,\dots$ for $x\in(-\infty,\infty)$ where $T_{n}$ are the Chebyshev polynomials. This orthogonal basis yields exponential convergence for functions that are analytic on the real line as well as banded differentiation matrices. A related candidate for rational basis is called the Malmquist-Takenaka (MT) functions [28, 29, 30]. It was shown in [31] that the MT functions are actually equivalent to the Chebyshev polynomials through a linear combination. As such, they inherit the many desirable properties mentioned above. Moreover, this family of functions are mutually orthogonal and the coefficients can be computed via Fourier transform.

Our method highlights the unstable nature of the Peregrine soliton. This instability is an existential threat to the formation of the Pergrine soliton solution over long time scales. We then explore on to perturbing constant background solutions. Localized perturbations are found to generate rational Peregrine-like structures that are localized in both space and time. The latter finding suggests that localized disturbances to periodic wave trains e.g. a strong wind gust, are capable of forming rogue waves.

The outline of the work is as follows. In Section 2 we introduce the MT functions and some of their properties. We discuss an expansion in terms of these functions and the rapid decay of their coefficients. In Section 3 spectral Galerkin differentiation matrices are constructed to approximate derivatives. In Section 4 we discuss a spectral Galerkin method for approximating PDEs. In particular, we establish the split-step method for approximating the dynamics of the NLS equation. In Section 5 we explore instability in the Peregrine soliton and mechanisms for generating rogue waves on a constant background.

2 The Malmquist-Takenaka functions

In this section, some properties of the MT functions are reviewed. Define the family of Malmquist–Takenaka rational functions

$$\phi_{n}(x)=\sqrt{\frac{2}{\pi}}i^{n}\left(\frac{1+2ix}{1-2ix}\right)^{n}\frac{1}{1-2ix},$$ (1)

for $n\in{\mathbb{Z}}$. First, observe that $\phi_{n}(x)$ is well-defined for $x\in\mathbb{R}$ with poles at $x=-i/2$ ($x=i/2$) when $n$ is even (odd). Several MT functions are plotted in Fig. 1. Below is a list of some properties of these functions.

Figure 1: The real and imaginary parts of $\phi_{n}(x)$ for $n\in\{-3,\dots,2\}$ indicated by blue (solid) and red (dashed) curves, respectively.

• Boundedness and decay

  The MT functions are uniformly bounded on the real line, that is

  $$|\phi_{n}(x)|=\sqrt{\frac{2}{\pi}}\frac{1}{\sqrt{1+4x^{2}}}\leq\sqrt{\frac{2}{\pi}}.$$ (2)

  The algebraic decay rate of the MT functions as $x\rightarrow\pm\infty$ is linear since

  $$|\phi_{n}(x)|=\sqrt{\frac{2}{\pi}}\frac{1}{\sqrt{1+4x^{2}}}\rightarrow\frac{1}{\sqrt{2\pi}}\frac{1}{|x|}.$$ (3)

• Oscillations

  The frequency of oscillations is directly related to the magnitude of the modal value $|n|$. As $|n|\rightarrow\infty~{}(\rightarrow 0)$, the functions become more (less) oscillatory near the origin. Express the rational function that lies on the unit circle by $e^{i\theta}=\frac{1+2ix}{1-2ix}$. Then $\phi_{n}(x)$ contains a factor of the form $e^{in\theta}=\cos(n\theta)+i\sin(n\theta)$ whose angular frequency increases with $|n|$.

• Symmetry

  When $n$ is even, the real part of $\phi_{n}(x)$ is even and the imaginary part is odd. On the other hand, if $n$ is an odd integer, then the real part of $\phi_{n}(x)$ is odd and the imaginary part is even. That is,

  $$n~{}\text{even:} ~{}~{}~{}~{}~{}~{}~{}\Re(\phi_{n}(-x))=\Re(\phi_{n}(x)),~{}~{}~{}~{}\Im(\phi_{n}(-x))=-\Im(\phi_{n}(x))$$

  $$n~{}\text{odd:} ~{}~{}~{}~{}~{}~{}~{}\Re(\phi_{n}(-x))=-\Re(\phi_{n}(x)),~{}~{}~{}~{}\Im(\phi_{n}(-x))=\Im(\phi_{n}(x))$$

  This symmetry is evident in Fig. 1.

• Orthogonality

  Define the complex $L^{2}(\mathbb{R})$ inner product

  $$\langle f,g\rangle=\int_{-\infty}^{\infty}f^{*}(x)g(x)dx,$$ (4)

  where ^∗ denotes complex conjugation. Then for distinct modes, the MT functions are mutually orthogonal in $L^{2}({\mathbb{R}})$, that is

  $$\langle\phi_{m},\phi_{n}\rangle=\int_{-\infty}^{\infty}\phi^{*}_{m}(x)\phi_{n}(x)dx=\delta_{mn},$$ (5)

  where $\delta_{mn}$ denotes the Kronecker delta function. The coefficients in (1) are chosen to also ensure unit norm, i.e. $||\phi_{n}||=1$, for each $n$.

• Derivative recurrence relation

  The MT functions satisfy the following skew-symmetric recurrence relation [30, 34]

  $$\frac{d}{dx}\phi_{n}(x)=-n\phi_{n-1}(x)+i(2n+1)\phi_{n}(x)+(n+1)\phi_{n+1}(x).$$ (6)

  This formula naturally leads to tridiagonal differentiation matrices of the first derivative. Contrast this with Fourier methods which are diagonal, i.e. $\frac{d}{dx}\varphi_{n}(x)=ikn\varphi_{n}(x),$ where $\varphi_{n}(x)=e^{iknx}$. Higher-order derivatives can be obtained by differentiating this formula.

• Hilbert transform eigenfunction

  The MT functions are eigenfunctions of the Hilbert transform. The Hilbert transform is defined by

  $$H[u](x)=\frac{1}{\pi}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\!\int_{-\infty}^{\infty}\frac{u(y)}{x-y}dy,$$

  where the integral is defined in the Cauchy principal values sense

  $$\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\!\int_{-\infty}^{\infty}f(x)dx=\underset{\epsilon\to 0^{+}}{\lim}\left[\int_{-\infty}^{b-\epsilon}f(x)dx\hskip 3.61371pt+\hskip 3.61371pt\int_{b+\epsilon}^{\infty}f(x)dx\right].$$

  The MT functions satisfy the eigenvalue problem

  $$H[\phi_{n}](x)=-i\hskip 3.61371pt\text{sgn}(n)\phi_{n}(x)~{}~{}{\rm where}~{}~{}\text{sgn}(n)=\begin{cases}1&\text{if}\hskip 21.68121ptn\geq 0\\ -1&\text{if}\hskip 21.68121ptn<0\\ \end{cases}.$$ (7)

  A natural application of the MT functions has been in the approximation of solutions to the Benjamin-Ono equation which includes the Hilbert transform [37, 41].

Relationship to Fourier series

Throughout this work we consider expanding a function in terms of the MT functions, that is

$$f(x)=\sum_{n=-\infty}^{\infty}\check{f}_{n}\phi_{n}(x),\hskip 7.22743pt\text{where}\hskip 7.22743pt\check{f}_{n}=\langle\phi_{n},f\rangle=\int_{-\infty}^{\infty}\phi^{*}_{n}(x)f(x)~{}dx,$$ (8)

for the MT functions given in (1). As can be seen, the coefficients are obtained by projecting on an arbitrary mode and exploiting the orthogonality (5). Through a natural change of variable, this MT expansion can be related to the Fourier series of a periodic function. This is useful because it allows us to use the Fast Fourier transform (FFT) to compute the MT coefficients.

Start from coordinate transformation

$$e^{i\theta}=\frac{1+2ix}{1-2ix}~{}~{}~{}\Leftrightarrow~{}~{}~{}x=\frac{1}{2}\tan\left(\frac{\theta}{2}\right),$$ (9)

which defines a map from the real line $x\in(-\infty,\infty)$ to the finite interval $\theta\in(-\pi,\pi)$. Applying this change of variable to the MT functions (1) yields

$$\phi_{n}\left(\frac{1}{2}\tan\left(\frac{\theta}{2}\right)\right)=\sqrt{\frac{2}{\pi}}i^{n}e^{in\theta}\frac{1}{1-i\tan\left(\frac{\theta}{2}\right)}.$$

Next, observe that

$$\frac{1}{1-i\tan\left(\frac{\theta}{2}\right)}=\frac{\cos(\theta/2)}{\cos(\theta/2)-i\sin(\theta/2)}=e^{i\frac{\theta}{2}}\cos\left(\frac{\theta}{2}\right),$$

so that

$$\phi_{n}(\theta)=\sqrt{\frac{2}{\pi}}i^{n}e^{i\theta\left(n+\frac{1}{2}\right)}\cos\left(\frac{\theta}{2}\right).$$ (10)

Notice that the MT function in (10) is $2\pi$-periodic in the $\theta$. Through the change of variable (9), the projection integral in (8) becomes

$$\check{f}_{n}=\int_{-\infty}^{\infty}\phi^{*}_{n}(x)f(x)dx=\frac{(-i)^{n}}{2\sqrt{2\pi}}\int_{-\pi}^{\pi}f\left(\frac{1}{2}\tan\left(\frac{\theta}{2}\right)\right)\left(1-i\tan\left(\frac{\theta}{2}\right)\right)e^{-in\theta}d\theta.$$ (11)

Notice that the MT coefficients are precisely the Fourier coefficients of the (periodic) function $f\left(\frac{1}{2}\tan\left(\frac{\theta}{2}\right)\right)\left(1-i\tan\left(\frac{\theta}{2}\right)\right)$.

The MT coefficients of the function $f(x)$ are computed via

$$\check{f}_{n}=\mathcal{MF}[f(x)]=(-i)^{n}\sqrt{\frac{\pi}{2}}\mathcal{F}\left[f\left(\frac{1}{2}\tan\left(\frac{\theta}{2}\right)\right)\left(1-i\tan\left(\frac{\theta}{2}\right)\right)\right],~{}~{}~{}~{}~{}~{}n\in\mathbb{Z}$$ (12)

and the function is summed (see Eqs. (8) and (10)) by

$$f(x)=\mathcal{MF}^{-1}[\check{f}_{n}]=\sqrt{\frac{2}{\pi}}\frac{\mathcal{F}^{-1}\left[(i)^{n}\check{f}_{n}\right]}{1-i\tan\left(\frac{\theta}{2}\right)},$$ (13)

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the Fourier and inverse Fourier (semi-discrete) transforms, respectively. We refer to $\mathcal{MF}$ and $\mathcal{MF}^{-1}$ as the modified Fourier and inverse modified Fourier transforms, respectively, which we will utilize to compute the MT coefficients of a given function. To be clear, we compute the MT coefficients directly through Fourier transforms. Note that the divisor, $1-i\tan\left(\frac{\theta}{2}\right)$, in (13), is included because (12) is computing the Fourier transform of the function $f\left(\frac{1}{2}\tan\left(\frac{\theta}{2}\right)\right)\left(1-i\tan\left(\frac{\theta}{2}\right)\right)=f\left(x\right)\left(1-2ix\right)$, so one recovers the function $f(x)$ by dividing out the term, $1-i\tan\left(\frac{\theta}{2}\right)$. The discrete version of the Fourier transform are naturally computed using the fast Fourier transform (FFT) algorithm (see B).

Decacy rate of MT coefficients

In this work we focus on rational functions that are square integrable on the real line, i.e. $L^{2}(\mathbb{R})$. Equivalently, we consider rational functions of the form $p(x)/q(x)$, where $p(x)$ and $q(x)$ are polynomials of degree $r$ and $s$, respectively, where $r<s$. We assume $q(x)$ has no zeros in $\mathbb{R}$. We can also consider the case $r=s$ when $\frac{p(x)}{q(x)}$ approaches the finite value $\frac{p_{\infty}}{q_{\infty}}$ as $x\rightarrow\pm\infty$. This latter class of rational of functions is clearly not in $L^{2}(\mathbb{R})$, however the shifted function $\frac{p(x)}{q(x)}-\frac{p_{\infty}}{q_{\infty}}$ decays to zero and is in $L^{2}(\mathbb{R})$. Moreover, the functions $\frac{p}{q}$ and $\frac{p}{q}$ - $\frac{p_{\infty}}{q_{\infty}}$ have the same derivatives. Hence, we can indirectly consider rational functions even when the degree of the number and denominator are equal.

One essential ingredient for this spectral method, or any method for that matter, to be effective is the coefficients (see (8)) decay rapidly. In the case of Fourier approximations on the real line, the coefficient decay rate is well-known to be related to the distance between the real line and poles [35]. We can apply this approach to this method too by locating the poles of the periodic function in (11).

Figure 2: (a) Numerically computed MT coefficients for the function $f(x)=(1+x^{2})^{-1}$ in blue, the exact MT coefficients (14) are in red. (b) Numerically computed MT coefficients for the function $f(x)=(1+x^{4})^{-1}$ in blue, the decay rate of MT coefficients in red. Notice for the slowly decaying rational function, the coefficients decay exponentially fast.

As an example consider the function $f(x)=(1+x^{2})^{-1}$ on the real line. Under transformation (9), the poles of the function are located at $x=\pm i$ or $\theta=\pm 2i\tanh^{-1}(2)\approx\pm(\pi+1.09861i)$. Since the distance from the real line to the pole is approximately 1.09861, the Fourier-MT coefficients defined in (11) are predicted to decay exponentially with the form $|\check{f}_{n}|\sim e^{-1.09861|n|}$ as $n\rightarrow\pm\infty$. Actually, $e^{-1.09861|n|}=3^{-|n|}$ which agrees with the analytical formula for this series

$$\frac{1}{1+x^{2}}=-\sqrt{2\pi}\sum_{n=-\infty}^{-1}\frac{i^{n}}{3^{-n}}\phi_{n}(x)+\sqrt{2\pi}\sum_{n=0}^{\infty}\frac{i^{n}}{3^{n+1}}\phi_{n}(x).$$ (14)

Another example is the function $f(x)=(1+x^{4})^{-1}$ which has poles at $x=e^{i\frac{\pi}{4}(1+2m)}$ for $m=0,1,2,3$ or $\theta=2\tan^{-1}\left(2e^{i\frac{\pi}{4}(1+2m)}\right)$. Regardless of which root one takes, the distance from the real line to the imaginary part of the poles is approximately 0.641155. Similar to the previous case, the MT-Fourier coefficients decay like $|\check{f}_{n}|\sim e^{-0.641155|n|}$ as $n\rightarrow\pm\infty$. Note we predict a slower convergence rate here.

We compare these decay rates with the numerically computed values in Fig. 2 for the functions $(1+x^{2})^{-1}$ and $(1+x^{4})^{-1}$, respectively. Importantly, the MT coefficients of each function decay exponentially fast, reaching round-off error at around $N=80$ and $N=120$, a modest amount of grid points. The rapid decay of the MT coefficients for slowly decaying rational functions in $L^{2}(\mathbb{R})$ suggest they could form the basis of an efficient spectral method. The next section considers spectral differentiation using the MT expansion. This, as with all spectral methods, benefits from exponentially decaying coefficients.

3 Spectral differentiation

Let us now explore differentiation via MT spectral methods. To begin, consider a function $f(x)$ expressed in terms of MT functions as in (8) whose derivative is

$$f^{\prime}(x)=\sum_{n=-\infty}^{\infty}\check{f}_{n}\phi^{\prime}_{n}(x).$$

We wish to express this derivative function in terms of MT functions itself, that is

$$f^{\prime}(x)=\sum_{n=-\infty}^{\infty}\check{f}_{n}\phi^{\prime}_{n}(x)=\sum_{n=-\infty}^{\infty}\check{f}^{\prime}_{n}\phi_{n}(x),$$

where $\check{f}^{\prime}_{n}$ denotes the derivative function coefficients. Recall that the derivative of the MT functions satisfies the recurrence relation (6). As a result, the MT coefficients of the derivative can be expressed solely in terms of MT functions

$$f^{\prime}(x)=\sum_{n=-\infty}^{\infty}\check{f}^{\prime}_{n}\phi_{n}(x)=\sum_{n=-\infty}^{\infty}\check{f}_{n}\left[-n\phi_{n-1}(x)+i(2n+1)\phi_{n}(x)+(n+1)\phi_{n+1}(x)\right].$$

After a shift of indices the MT coefficients of the derivative are expressed in terms of the the coefficients of the function itself

$$f^{\prime}(x)=\sum_{n=-\infty}^{\infty}\check{f}^{\prime}_{n}\phi_{n}(x)=\sum_{n=-\infty}^{\infty}\left(n\check{f}_{n-1}+i(2n+1)\check{f}_{n}-(n+1)\check{f}_{n+1}\right)\phi_{n}(x).$$ (15)

Notice that the MT derivative involves self $\{n\}$ and nearest {$n-1,n+1$} modes; this is a direct consequence of recurrence relation (6). On the other hand, many spectral methods often involve completely dense interactions [40]. So, the MT differentiation matrices are relatively sparse, and as we shall see below, have the ability to converge exponentially fast for appropriate rational functions.

The derivative MT coefficients given in (15) are related to the coefficients of the original function by

$$\check{\textbf{f}}^{\prime}=\left(\dots,\check{f}^{\prime}_{-1},\check{f}^{\prime}_{0},\check{f}^{\prime}_{1},\dots\right)^{T},~{}~{}~{}~{}\check{\textbf{f}}^{\prime}=\mathbb{D}_{1}\check{\textbf{f}},~{}~{}~{}~{}\check{\textbf{f}}=\left(\dots,\check{f}_{-1},\check{f}_{0},\check{f}_{1},\dots\right)^{T},$$ (16)

where $\mathbb{D}_{1}$ is the banded matrix

$$\mathbb{D}_{1}=\left[\begin{array}[]{*{10}c}\ddots&\ddots&&&&&&&&\\ \ddots&\ddots&\ddots&&&&&&&\\ &-3&-5i&2&&&&&&\\ &&-2&-3i&1&&&&&\\ &&&-1&-i&0&&&&\\ &&&&0&i&-1&&&\\ &&&&&1&3i&-2&&\\ &&&&&&2&5i&-3&\\ &&&&&&&\ddots&\ddots&\ddots\\ &&&&&&&&\ddots&\ddots\end{array}\right].$$

Notice that this matrix is tridiagonal and skew-Hermitian, that is, $\mathbb{D}_{1}^{\dagger}=-\mathbb{D}_{1}$, where ^† denotes the conjugate transpose.

To approximate derivatives via MT expansions in practice, we first truncate the series in (15) to a finite range of $N$ modes: $-N/2,-N/2+1,\dots,N/2-1$, where $N$ is a positive even integer. Then one computes the MT coefficients of the function $f(x)$ using the discrete version of (12), given in (32). Next, the truncated differentiation matrix in (16) is applied. The discrete and finite version of the series (15), given in (33), is summed to give the derivative approximation. This process is summarized in A.

Second-order differentiation matrix

The differentiation matrix for approximating the second derivative can be derived in a similar fashion. Consider a function $f$ expressed in basis of MT functions, i.e. (8), whose second derivative is clearly

$$f^{\prime\prime}(x)=\sum_{n=-\infty}^{\infty}\check{f}_{n}\phi^{\prime\prime}_{n}(x).$$

The goal is to express this function itself in terms of MT functions, that is

$$f^{\prime\prime}(x)=\sum_{n=-\infty}^{\infty}\check{f}_{n}\phi^{\prime\prime}_{n}(x)=\sum_{n=-\infty}^{\infty}\check{f}^{\prime\prime}_{n}\phi_{n}(x),$$

where $\check{f}^{\prime\prime}_{n}$ denote the second derivative coefficients, which we intend to find. Differentiating the first derivative relation in (6) yields the five-term second derivative relation

$$\frac{d^{2}}{dx^{2}}\phi_{n}(x) =-n\frac{d}{dx}\phi_{n-1}(x)+i(2n+1)\frac{d}{dx}\phi_{n}(x)+(n+1)\frac{d}{dx}\phi_{n+1}(x)$$ (17)

$$=n(n-1)\phi_{n-2}(x)-4in^{2}\phi_{n-1}(x)-(2+6n(n+1))\phi_{n}(x)+4i(n+1)^{2}\phi_{n+1}(x)+(n+2)(n+1)\phi_{n+2}(x).$$

Note that this involves five nodal values, $\{n-2,n-1,n,n+1,n+2\}$. Similar to first-order derivative, the second derivative coefficients are obtained by exploiting the orthogonality of the functions and shifting the indices of the series and thus the MT coefficients of the second derivative can be represented as

$$\sum_{n=-\infty}^{\infty}\check{f}^{\prime\prime}_{n}\phi_{n}(x)=\sum_{n=-\infty}^{\infty}\left[n(n-1)\check{f}_{n-2}+4in^{2}\check{f}_{n-1}-(2+6n(n+1))\check{f}_{n}-4i(n+1)^{2}\check{f}_{n+1}+(n+2)(n+1)\check{f}_{n+2}\right]\phi_{n}(x).$$ (18)

Using this information, the MT coefficients of the second derivative can be computed via

$$\check{\textbf{f}}^{\prime\prime}=\left(\dots,\check{f}^{\prime\prime}_{-1},\check{f}^{\prime\prime}_{0},\check{f}^{\prime\prime}_{1},\dots\right)^{T},~{}~{}~{}~{}\check{\textbf{f}}^{\prime\prime}=\mathbb{D}_{2}\check{\textbf{f}},~{}~{}~{}~{}\check{\textbf{f}}=\left(\dots,\check{f}_{-1},\check{f}_{0},\check{f}_{1},\dots\right)^{T},$$ (19)

where $\mathbb{D}_{2}$ is the banded matrix

$$\mathbb{D}_{2}=\left[\begin{array}[]{*{16}c}\ddots&\ddots&\ddots&&&&&&&&&&&&&\\ \ddots&\ddots&\ddots&\ddots&&&&&&&&&&&&\\ \ddots&\ddots&\ddots&\ddots&\ddots&&&&&&&&&&&\\ &30&100i&-122&64i&12&&&&&&&&&&\\ &&20&64i&-74&36i&6&&&&&&&&&\\ &&&12&36i&-38&16i&2&&&&&&&&\\ &&&&6&16i&-14&-4i&0&&&&&&&\\ &&&&&2&4i&-2&0&0&&&&&&\\ &&&&&&0&0&-2&-4i&2&&&&&\\ &&&&&&&0&4i&-14&-16i&6&&&&\\ &&&&&&&&2&16i&-38&-36i&12&&&\\ &&&&&&&&&6&36i&-74&-64i&20&&\\ &&&&&&&&&&12&64i&-122&-100i&30&\\ &&&&&&&&&&&\ddots&\ddots&\ddots&\ddots&\ddots\\ &&&&&&&&&&&&\ddots&\ddots&\ddots&\ddots\\ &&&&&&&&&&&&&\ddots&\ddots&\ddots\end{array}\right].$$

Notice that this matrix is pentadiagonal and Hermitian ($\mathbb{D}_{2}^{\dagger}=\mathbb{D}_{2}$). Similar to $\mathbb{D}_{1}$, this differentiation matrix is banded and fairly sparse; this is a general theme. Based on recurrence relation (6), increasing the order of the derivative increases the bandwidth of the differentiation matrix by one super and one sub diagonal. To implement this in practice, the system in (19) is truncated to some sufficiently large set of modes and the truncated version of $\mathbb{D}_{2}$ is applied (see A).

A natural question to ask is whether to approximate $\partial_{x}^{2}$ by $\mathbb{D}_{1}^{2}=\mathbb{D}_{1}\mathbb{D}_{1}$ (two applications of the first derivative matrix (16)), rather than $\mathbb{D}_{2}$. When the system is infinite (no truncation), $\mathbb{D}_{2}=\mathbb{D}_{1}^{2}$. However, upon truncation the matrices differ at the first and last modes due to different boundary conditions. In our experience, approximation by $\mathbb{D}_{2}$ is slightly more accurate (by about a factor of 2), but one can use $\mathbb{D}_{1}^{2}$ and still expect spectral convergence rates.

Figure 3: Infinity norm error of the (a) first and (b) second MT derivative approximations applied to the rational function $f(x)=(1+x^{2})^{-1}$ as a function of the number of MT modes, $N$. The error is observed to converge exponentially fast. All derivatives were computed using

A typical set of convergence results is highlighted in Fig. 3. For the rational function $f(x)=1/(1+x^{2})$, the MT approximation of derivatives is remarkably effective. The error converges exponentially fast to zero until we reach round-off error effects at around $N=64$ modes. We notice with each order of derivative, the error increases by about one order of magnitude. These rapid convergence results are expected based on the exponentially decay of the MT coefficients shown in Fig. 2. Given that we are using fast modified Fourier transforms in (12) and (13), these calculations are performed quickly.

4 Numerical approximation of PDEs

In this section we develop a method for numerically approximating rational solutions of nonlinear PDEs. An explicit split-step time-stepping approach is implemented using the differentiation matrices described in Sec. 3. These techniques will be used to approximate rational solutions of a well-known PDE, the nonlinear Schrödinger (NLS) equation.

To begin, consider the general nonlinear initial boundary value problem

$$u_{t}=\mathcal{L}u+\mathcal{M}[u],~{}~{}~{}~{}~{}u_{0}(x)=u(x,t_{0}),~{}~{}~{}~{}~{}(x,t)\in\mathbb{R}\times[t_{0},T],$$ (20)

where the linear and nonlinear operators are denoted by $\mathcal{L}$ and $\mathcal{M}$, respectively. The linear operator is assumed to be a constant coefficient differential operator with the form

$$\mathcal{L}=c_{0}+c_{1}\frac{\partial}{\partial x}+c_{2}\frac{\partial^{2}}{\partial x^{2}}+\dots~{},$$

for scalar coefficients $c_{0},c_{1},\dots$ The nonlinear operator is assumed to be local and depend on some product of $u$ and it’s derivatives. If both operators in (20) are linear ($\mathcal{M}[u]=Mu$) and constant coefficient, then the solution is

$$u(x,t)=\exp\left[(t-t_{0})(\mathcal{L}+\mathcal{M})\right]u_{0}(x).$$ (21)

We consider solutions that approach a constant value at the boundaries, i.e. $u(x,t)\rightarrow U_{b}$ as $x\rightarrow\pm\infty$. Ultimately we have in mind rational solutions where the degree of the denominator polynomial is greater than or equal to the degree of the numerator polynomial.

4.1 MT-SS4 approximation of the NLS equation

We now introduce a fourth-order Malmquist-Takenaka split-step (MT-SS4) method for integrating nonlinear PDEs. The idea here is to split the linear and nonlinear parts of (20) into two parts that are easier to solve separately. This method has the benefit of being quite accurate and straightforward to implement while also possessing a relatively large stability region. A thorough derivation of the split-step and its properties can be found in [20].

Consider splitting (20) into two equations

$$\frac{\partial{\bf w}}{\partial t}=\mathcal{L}{\bf w},~{}~{}~{}~{}~{}~{}~{}~{}\frac{\partial{\bf v}}{\partial t}=\mathcal{M}[{\bf v}].$$ (22)

The idea behind the method is to solve the equations in (22) separately and then re-combine them in a specific order to approximate the solution. For example, the first-order approximation of (21) for $t_{0}=0$ is $\exp(t(\mathcal{L}+\mathcal{M}))\approx\exp(t\mathcal{L})\cdot\exp(t\mathcal{M}).$ Only when $\mathcal{L}$ and $\mathcal{M}$ commute does this become an equality. That is, solve one equation (22) and then use its solutions as initial condition for the other equation. The nonlinear version operates in a similar manner [20].

The fourth-order split-step method for one time-step $\Delta t$ is given by [20, 42]

$${\bf u}(t+\Delta t)\approx\exp\left[\alpha_{4}\Delta t\mathcal{M}\right]\exp\left[\beta_{3}\Delta t\mathcal{L}\right]\exp\left[\alpha_{3}\Delta t\mathcal{M}\right]\exp\left[\beta_{2}\Delta t\mathcal{L}\right]\exp\left[\alpha_{2}\Delta t\mathcal{M}\right]\exp\left[\beta_{1}\Delta t\mathcal{L}\right]\exp\left[\alpha_{1}\Delta t\mathcal{M}\right]{\bf u}(t),$$ (23)

in terms of the coefficients

$$\alpha_{1}=\frac{1}{2}c,\hskip 3.61371pt\alpha_{2}=\frac{1}{2}(1-c),\hskip 3.61371pt\alpha_{3}=\alpha_{2},\hskip 3.61371pt\alpha_{4}=\alpha_{1},\hskip 3.61371pt\beta_{1}=c,\hskip 3.61371pt\beta_{2}=1-2c,\hskip 3.61371pt\beta_{3}=\beta_{1},$$

where $c=\frac{1}{2-2^{1/3}}$. More specifically, starting from the solution at time $t$, ${\bf u}(t)$, we first integrate the $\mathcal{M}$-operator equation (22) by a step $\alpha_{1}\Delta t$ to get an intermediate solution ${\bf v}_{1}$. Then we integrate the $\mathcal{L}$-operator equation (22) by a step $\beta_{1}\Delta t$, using ${\bf v}_{1}$ as the initial condition to get another intermediate solution ${\bf w}_{1}$. This process repeats itself until the final step in (23).

Note that the overall accuracy of the method does not depend on the choice of which operator is first and second. That is, the fourth-order method in (23) starts and ends with $\mathcal{M}$, but it is equally valid to begin and end with $\mathcal{L}$, and the order of the method will not be affected. Due to its simplicity, for the NLS equation below we choose to solve the nonlinear part of four times ($\alpha_{n}$ coefficients) and the linear part three times ($\beta_{n}$ coefficients). Whenever possible the exponential operators in (23) should be pre-computed.

From a practical point-of-view, this is an effective technique only when the equations in (22) are solved efficiently and accurately. Typically, the linear part of equation in (22) can be rapidly solved using some spectral approximation, e.g. Fourier or MT. The crucial step is the nonlinear one. For an equation like NLS, this is a nearly trivial solve (see next section) [38]. However, for other equations, like those in the Korteweg-de Vries family, solving this is not so straightforward and this method is not as attractive [39].

Figure 4: A top (left) and side (right) view of the modulus using MT-SS4 method for the NLS equation (24) and Peregrine initial condition (26) with $t_{0}=-2$. These results were obtained using $N=256$ MT modes and a time-step of $\Delta t=0.0001$.

Consider the nonlinear Schrödinger equation

$$iu_{t}+u_{xx}+2(|u|^{2}-1)u=0,~{}~{}~{}~{}~{}u_{0}(x)=u(x,t_{0}),~{}~{}~{}~{}(x,t)\in\mathbb{R}\times[t_{0},T],$$ (24)

with boundary conditions $u\to 1$ as $x\to\pm\infty$. Note that through the phase transformation $u(x,t)=\psi(x,t)e^{-2it}$ this equation can be recast in the more standard form

$$i\psi_{t}+\psi_{xx}+2|\psi|^{2}\psi=0.$$

The form in (24) is more convenient for our purposes. Two well-known solutions of the NLS equaiton are the plane wave solution

$$u_{c}(x,t)=1,$$ (25)

and Peregrine soliton solution

$$u_{p}(x,t)=\frac{4x^{2}-16it+16t^{2}-3}{4x^{2}+16t^{2}+1},$$ (26)

both of which are rational functions with unity boundary conditions. Physically, (25) corresponds to a nonlinear periodic wave train. The Peregrine solution in (26) is a model for rogue waves; a depiction of it’s dynamics is shown in Fig. 4. Unlike most solitons that maintain their profile during propagation, the Peregrine soliton is localized in both space and time.

Consider solving the NLS equation (24) solution via split-step method (23) where $\mathcal{L}=i\partial_{x}^{2}$ and $\mathcal{M}[u]=2i\left(|u|^{2}-1\right)$ in (20). Explicitly, the split-step equations in (22) are

$$\frac{\partial w}{\partial t}=i\frac{\partial^{2}w}{\partial x^{2}},~{}~{}~{}~{}~{}~{}~{}~{}\frac{\partial v}{\partial t}=2i\left(|v|^{2}-1\right)v.$$ (27)

The second equation in (27) is clearly solved by

$$v(x,t_{0}+\Delta t)=e^{2i\Delta t\left(|v(x,t_{0})|^{2}-1\right)}v(x,t_{0}).$$ (28)

To numerically solve the first equation in (27) we wish to approximate spatial derivative using an MT spectral approach. Normally the first step is to expand

Figure 5: Top (left) and side (right) views of the solution modulus using MT-SS4 method. All parameters are the same as Fig. 4 except the time-step $\Delta t=0.1$. The large time-step reveals an instability in the approximation for $t\geq 0$.

the solution $w(x,t)$ in basis of MT functions (see (8)). However, the solution $w(x,t)\rightarrow 1$ as $|x|\rightarrow\infty$ and is not in $L^{2}(\mathbb{R})$; an expansion in terms of the MT functions requires functions decay to zero at infinity. As a result, the functions in (25) and (26) cannot directly be expressed in terms of a MT expansion.

To get around this barrier, we first observe the trivial calculus fact

$$\frac{df}{dx}=\frac{d}{dx}(f-1),$$

for some differential function $f(x)$ such that $f(x)\rightarrow 1$ as $|x|\rightarrow\infty$. Crucially, even though both sides of the equation have the same derivatives, the function on the right-hand side has zero boundary conditions and can be directly expanded in a basis of MT functions. That is, we make the simple change of variables $\widetilde{w}(x,t)=w(x,t)-1$ to the linear (first) equation in (27), solve, and then return the unit boundary conditions. In this way we can approximate derivatives of functions $u(x)$ with constant and nonzero boundary conditions, without resorting to a reformulation of the problem.

After shifting by one, we expand $\widetilde{w}(x,t)$ as (8), and apply the differentiation matrix given in (19). The system is truncated to $N$ MT modes and the linear equation is reformulated as

$$\frac{d\check{\bf{w}}}{dt}=i\mathbb{D}_{2}{\check{\bf w}},$$

where $\mathbb{D}_{2}$ is the truncated $N\times N$ differentiation matrix given in (19) and $\check{{\bf w}}$ consists of the MT coefficients associated with the function $\widetilde{w}(x,t)$. This is a finite system of ODEs and can be solved exactly by

$$\check{{\bf w}}(t_{0}+\Delta t)=e^{i\Delta t\mathbb{D}_{2}}\check{{\bf w}}(t_{0}).$$ (29)

Afterward, the solution in physical space is recovered by computing a truncated sum from (8) (see discrete version in (33)).

Figure 6: (a) Convergence of the MT-SS4 approximation for the Peregrine soliton (26) solved on the interval $-2\leq t\leq 2$. The error at time $t=2$ for $N=256$ MT nodes is shown. The numerical approximation apparently converges at a fourth-order rate when time-step $\Delta t$ is decreased, until other errors become significant. (b) Average CPU runtime as a function of $\Delta t$ and fixed $N$. The runtime clearly increases linearly with $\Delta t$.

4.2 A convergence study

Typical results obtained by solving the NLS equation (24) with the MT-SS4 method for the Pergrine soliton solution (26) are discussed next. The evolution of a typical Pergrine soliton obtained using this method is shown in Fig. 4 over the time interval $[-2,2]$. Notice that the soliton peaks at the origin $(x,t)=(0,0)$.

A benefit of the split-step method is the relatively large stability region. That is, the method yields stable results for quite large values of $\Delta t$.

A numerically induced stability is shown in Fig. 5. This is the result of taking an apparently too large time-step. The derivation of a rigorous stability bound for the MT-SS4 method is an open problem.

Figure 7: Final time snapshot of the real (a) and imaginary (b) parts of the MT-SS4 approximation for the Peregrine soliton shown in Fig. 4; computational parameters are the same. (c) Absolute difference between MT-SS4 approximation and exact solution at $t=2$.

Next, the rate of convergence and runtime as a function of $\Delta t$ is discussed. In Fig. 6(a) the convergence rate of the MT-SS4 method applied to the NLS equation (24) for the Peregrine solution (26) is shown. A clear fourth-order convergence rate is observed until other errors become significant. We have observed some apparent aliasing errors for certain $\Delta t$ values. These anomalies and a padding approach to mitigate them is discussed in Sec. 4.3.

The corresponding runtimes are shown in Fig. 6(b). This timer encapsulates only the split-step solver portion of the code. We observe a roughly linear dependence for this range of computational parameters. For all runs considered in Fig. 6, the total CPU time ranged from 0.5-10 seconds. All runs were performed on a HP Laptop personal computer with an 1.60 GHz 4 Core Intel processor with 8.00 GB of Ram and 15.7 GB of memory. The data values are the average of over 20 runs to mitigate various CPU-related fluctuations in the runtime.

A comparison between the exact and numerical values is highlighted in Fig. 7 at the final time $t=2$. For a time-step size of $\Delta t=10^{-4}$, the difference between the exact and numerical approximation is $\mathcal{O}(10^{-9})$ or better over this time interval. The largest source of the error is observed to occur near the soliton peak located near the origin. We point out that under transformation (9), the spatial points cluster near the origin, even though the points in $\theta$ are evenly spaced. In practice, we observe it essential to cluster points near nontrivial portions of the solution, e.g. rogue wave peak, to achieve maximal accuracy.

Lastly, we do not include any results for the plane wave solution (25) as it is solved exactly by the MT-SS4 method. Examining (28) we see that the nonlinear integrator simply returns the initial plane wave solution, i.e. $v(x,t+\Delta t)=v(x,t)$. For the linear integrator (4.1), we note that $u_{c}$ has boundary condition $U_{b}=1$ and so $\widetilde{w}(x,t)=w(x,t)-1=0$. As a result, all MT coefficients are zero, $\check{\bf w}={\bf 0}$, and the (trivial) solution to (4.1) is zero. Upon returning the boundary condition $U_{b}=1$, we recover the exact constant solution, and so on.

4.3 Aliasing errors

In this section we address an aliasing-type error that can occur with this method. Aliasing effects are a well-known consequence in certain spectral method approximations of nonlinear PDEs [18]. In our case, sometimes when implementing the MT-SS4 algorithm for a given time-step $\Delta t$, we observe anomalously large errors relative to similar time-step sizes. This is due to unexpectedly large spectral coefficient values corresponding to high-frequency MT modes. We propose a spectral padding approach to mitigate these effects.

The error convergence rate for the MT-SS4 method applied to the NLS equation (24) is shown in Fig. 8 for the Peregrine soliton. The only difference between this graph and that of Fig. 6(a) is that we include more values of $\Delta t$. Specifically, in Fig. 6(a) we decrease $\Delta t$ by a factor of 10/7 each time, whereas in Fig. 8(A) we multiply the value of $\Delta t$ by $0.99$ each time (equivalently, decrease by a factor of $100/99$). The overall convergence rate appears fourth-order except for a few problematic $\Delta t$ values. The most egregious example occurs at (b): $\Delta t=0.001(0.99)^{60}$ where the error is ${\mathcal{O}}(10^{-6})$, meanwhile a slightly larger value of (c): $\Delta t=0.001(0.99)^{59}$ or slightly smaller value of (a): $\Delta t=0.001(0.99)^{61}$ have errors that are ${\mathcal{O}}(10^{-7})$. This order of magnitude increase of error is concerning and requires further investigation.

To understand the source of this unexpectedly large error, we examine the evolution of the MT coefficients of the function $u(x,t)-1$. In Fig. 9(a) the modulus of the MT-coefficients at the initial condition $u(x,-2)-1$ for the solution given in (26) is shown. In the other panels we display the MT coefficients at $t=2$ (final time) for different choices of $\Delta t$. In Fig. 9(c) the MT coefficients for the problematic case are observed to be unexpectedly large at large modal values, i.e. at large values of $n$ in (1). The MT functions $\phi_{n}(x)$ at large values of $|n|$ correspond to highly oscillatory functions (see Fig. 1). In contrast, the MT coefficient profiles in Figs. 9(b) and 9(d) appear to remain relatively small at large values of $n$.

The discrepancy in the MT coefficients for $\Delta t=0.001(0.99)^{60}$ compared to the others appears to be an aliasing-type error. For certain values of $\Delta t$, high MT modes (large $n$) are unexpectedly large. This is a result of the nonlinearity in the NLS equation. For Fourier spectral methods, there are a myriad of techniques for de-aliasing, see [18]. One of the simplest techniques, which we will employ, is that of padding. Spectral padding involves setting a certain number of high MT modes (large $n$) to zero every time-step. Explicitly, we apply the following filter once at each linear solve in the MT-SS4

$$\tilde{u}_{n}=\begin{cases}\check{u}_{n},&|n|<\frac{N}{2}-P\\ 0,&\text{otherwise}\\ \end{cases},$$ (30)

for $P\geq 0$ whereby all MT modes with corresponding modenumber $|n|$ greater than or equal to $N/2-P$ ($N$ is even) are effectively truncated.

Figure 8: Convergence of the MT-SS4 approximation for the Peregrine soliton (26). Other than the time-step $\Delta t$, all computational parameters are the same as Fig. 6. Overall, a fourth-order convergence rate is observed for both cases. (A) The numerical approximation converges as expected except for a few problematic $\Delta t$ values. Specific points are indicated: (a) $\Delta t=0.001(0.99)^{61}$, (b) $\Delta t=0.001(0.99)^{60}$, and (c) $\Delta t=0.001(0.99)^{59}$. (B) Convergence rate of the padded MT-SS4 method for the Peregrine soliton (26). The spectral padding defined in (30) with $P=1$ has been applied; otherwise, this is the same as (A).

Figure 9: Modulus of MT coefficients without padding for the exact solution given in (26) with $N=256$ MT modes. (a) MT coefficients at the initial condition $u(x,-2)-1$. Numerical and exact MT coefficients at final time $t=2$ for (b) $\Delta t=0.001(0.99)^{59}$, (c) $\Delta t=0.001(0.99)^{60}$, and (d) $\Delta t=0.001(0.99)^{61}$.

The convergence rate for the padded MT-SS4 method is shown in Fig. 8(B) applied to the NLS equation. Here we take a minimal padding ($P=1$) which zeros out the $n=-N/2,-N/2+1,N/2-1$ modes. This approach is found to prevent the growing tails observed in Fig. 9(c). The convergence rate is fourth-order and should be compared with the convergence rate in Fig. 8(A). Unlike the unpadded case, we do not observe any sharp jumps in error for the $\Delta t$ values considered before. In the event it does appear, the size of pad can be extended by increasing $P$ in (30), thereby truncated more high modenumber MT coefficients.

5 Instability and rogue waves in the NLS equation

Having developed the MT-SS4 method in the previous sections, we now look to apply it on a worthy problem. In this section we explore two models for rogue waves: the Peregrine soliton (26) and instability. The Peregrine soliton is found to be unstable. This will limit its ability as a robust rogue wave mechanism. Next, we examine a perturbation of the constant background. A variety of localized perturbations are found to generate unstable rogue wave, Peregine-like modes.

The instability of the Peregrine soliton

To begin, we examine the error introduced by the MT-SS4 method when approximating the Peregrine soliton. The infinity-norm error as a function of time is shown in Fig. 10. Importantly, the error is observed to grow exponentially fast. For linearly stable problems, often one can bound the error as a function of term that grows algebraically in $t$ [43]. This is a serious error growth. Decreasing the time-step size is found to postpone, but not eliminate this error growth. We have included a set of fits in Fig. 10 that illuminate that the error growth rate, before and after the peak, is on the order of $e^{2t}$. Countless simulations reveal that, regardless of the number of MT nodes $N$, time-step $\Delta t$, or initial time $t_{0}$, the growth rate of the error is eventually $\mathcal{O}(e^{2t})$.

So then, what is the nature of the significant error? Here we attribute it not necessarily to a deficiency in our method per se, but rather to the intrinsically unstable nature of the Peregrine soliton and the NLS equation with nonzero boundary conditions. Several recent works [14, 16, 44] have highlighted the unstable nature of the Peregrine solution. Moreover, in [16, 46] the authors precisely predicted a linear instability growth rate of $e^{2t}$. Hence, the MT-SS4 method has revealed this expected instability growth rate.

Figure 10: Error over time interval $[-4,6]$ using the MT-SS4 for the NLS equation (24) and initial condition $u(x,t_{0})$ where $t_{0}=-4$ for the solution given in (26). The numerical method uses $N=512$ MT modes and $\Delta t=0.0001$.

This instability will most likely always be a challenge in simulation of the Peregrine solution. The reason is that numerics always makes some sort of approximation, e.g. truncation, round-off error. Hence, the moment we try to simulate this solution, regardless of method, we introduce some error. Decreasing the time-step or increasing the number of modes can mitigate these errors, but eventually they will become significant.

In a physical context, we conjecture this instability will make it difficult for Peregrine soliton to form naturally. In a controlled environment, like a wave flume, these waves have been realized [45]. However, in the open ocean a wide variety of disturbances such as wind or collision with any object will certainly introduce perturbations that will affect the formation of the genuine solution in (26). A more likely scenario is a disturbance of a periodic wave train, or the constant solution (25). We explore this next and find that it generates Peregrine-like structures.

Rogue waves generated by instability

The constant background of the NLS equation is unstable to perturbations and so next we examine the rogue wave-like states that form from this instability. Recall that the Peregrine soliton (26) itself approaches the constant solution at large times when $t\rightarrow\pm\infty$. In this section we consider several examples of localized perturbations of the constant solution. This would be akin to an ocean wave that is perturbed by a gust of wind or another similar stimuli.

Starting at an initial time of $t_{0}=-2$ we perturb the constant solution, (25), with initial condition

$$u_{j}(x,-2)=1+\varepsilon f_{j}(x),~{}~{}~{}~{}j=1,2,3,4,$$ (31)

where $0<\varepsilon\ll 1$ and

$$f_{1}=\frac{1}{1+x^{2}},~{}~{}~{}~{}~{}~{}f_{2}=\frac{1}{1+x^{4}},~{}~{}~{}~{}~{}~{}f_{3}=e^{-x^{2}},~{}~{}~{}~{}~{}~{}f_{4}=\text{sech}(x).$$

and then evolve the equation using the MT-SS4 method. In each case considered we choose $\varepsilon=0.1.$ For each case considered, we observe single-hump Peregrine-like structures peak at approximately $t=0$, with a peak height of around 2.6-2.8. A typical evolution is shown in Fig. 11. Notice the peak amplitude is a factor of 26-28 larger than the initial perturbation amplitude ($\varepsilon=0.1$). The structure that forms around $t=0$ has an uncanny resemblance to that of the Pergreine soliton (see Fig. 4). Namely, a localized peak of approximate amplitude 2.7 appears and then suddenly disappears. Moreover, near the peak, we observe a set of adjacent troughs. We do not dwell on them here, but near the final time ($t\approx 1$), a structure consisting of two-humps forms. This is reminiscent of modulational instability on periodic finite domains [47].

To further establish the Peregrine (rational) nature of this peak, in Fig. 12 we show some solution profiles at $t=0$ (near the peak time) generated by each initial condition in (31). In each case considered, the overall decay rate of $|u(x,0)-1|$ is found to decay rationally. Notice that we have subtracted off the background to isolate the decay rate. Furthermore, when we apply a variety of rational fits, we consistently find the tails to decay like $\mathcal{O}\left(\frac{1}{x^{2}}\right)$ as they approach the constant background. This matches the rational structure of (26) at $t=0$.

These results suggest a likely mechanism for generating rogues wave is this instability. We have established the first peak is nearly the Peregrine soliton at $t=0$. All four of initial perturbations considered in (31) appear to be generating roughly the same rational structure (see Fig. 11(left)). These results were obtained by applying the MT-SS4 method developed in the first part of this paper.

Figure 11: Typical results obtained by applying the MT-SS4 using the perturbed initial condition (31) with $f_{1}(x)=\frac{1}{1+x^{2}}$. Top (left) and side (right) views of the solution modulus. The computational parameters used are $N=512$,$\Delta t=0.0001$, and with a pad of $P=1$.

Figure 12: Plots of $\ln|u(x,0)-1|$ obtained by solving the NLS via the MT-SS4 method with the perturbed initial conditions (31). The computational parameters are the same as Fig. 11. These results are obtained from initial perturbations (a) $f_{1}(x)=(1+x^{2})^{-1}$, (b) $f_{2}(x)=(1+x^{4})^{-1}$, (c) $f_{3}(x)=e^{-x^{2}}$, (d) $f_{4}(x)=\text{sech}(x)$. The rational function $\ln|(1+x^{2})^{-1}|$ is indicated by curve (e); it indicates that the decay rate of all the perturbed approximations is rational and decays like $\frac{1}{x^{2}}$ as $x\rightarrow\pm\infty$.

6 Conclusion

In this work we have developed a numerical method for approximating rogue wave models of certain partial differential equations. These rogue wave models are rational in form, so we expand them in a basis of Malmquist-Takenaka rational functions. The spectral coefficients can decay exponentially fast and can be numerically computed in terms of (fast) Fourier methods. A fourth-order split-step method for approximation the nonlinear Schrödinger equation was introduced. The accuracy of the method for rational solutions was observed to be spectral in space and fourth-order in time.

This method was applied to approximating rogue wave models. The Peregrine soliton is linearly unstable and the instability manifests itself in the simulations. Any slight perturbation will excite this instability making the genuine Peregrine soliton solution challenging to see in a deep water, open ocean setting. On the other hand, when a constant background (traveling plane wave) is perturbed, a Peregrine-like rogue wave is found to emerge.

Appendix A Malmquist-Takenaka approximation of derivatives

In this appendix we give the general algorithm for approximating spatial derivatives using the MT spectral approximation. The semi-discrete motivations are shown in (12) and (13). The discrete modified Fourier transforms are computed using standard FFT algorithms. We abbreviate the modified FFT by MFFT and the modified inverse Fourier transform by MIFFT. Their precise definitions are given in B.

$$\text{Input:}~{}~{}~{}N\text{(even)},\hskip 10.84006pt\text{function }f(x),\hskip 14.45377pt\text{differentiation matrix }\mathbb{D}_{m}$$

$$x_{j+\frac{N}{2}+1}=\frac{1}{2}\tan\left(\frac{\pi j}{N}\right),~{}~{}~{}~{}j=-\frac{N}{2},-\frac{N}{2}+1,...,\frac{N}{2}-1$$

$$\textbf{x}=(x_{1},x_{2},...,x_{N})^{T}$$

$$\textbf{f}=(f_{1},f_{2}.,...,f_{N})^{T},\hskip 7.22743ptf_{j}\approx f(x_{j})$$

$$\check{\textbf{f}}=\text{MFFT}[\textbf{f}]$$

$$\check{\textbf{f}}^{(m)}=\mathbb{D}_{m}\check{\textbf{f}}$$

$$\textbf{f}^{(m)}=\text{MIFFT}\left[\check{\textbf{f}}^{(m)}\right]$$

$$\text{Output:}~{}~{}~{}\textbf{f}^{(m)}=(f_{1}^{(m)},f_{2}^{(m)},...,f_{N}^{(m)})^{T}$$

This algorithm gives derivative approximations at the equally spaced angular points $\theta_{j}=\frac{2\pi j}{N}$ for $j=-N/2,-N/2+1,\dots,N/2-1$. Differentiation matrices for the first two derivatives are shown in (16) and (19). The output is the approximation of the $m^{\rm th}$-order derivative evaluated at the grid points, i.e. $f_{j}^{(m)}\approx f^{(m)}(x_{j})$. Once the MT coefficients are obtained, one can approximate $f(x)$ and its derivatives at arbitrary values of $x$ (not just the spatial grid points) by summing the truncated MT series $f^{(m)}(x)\approx\overset{N/2-1}{\underset{n=-N/2}{\Sigma}}\check{f}^{(m)}_{n}\phi_{n}(x),~{}m=0,1,2,$ for the MT functions given in (1).

Appendix B Discrete modified Fourier and inverse Fourier transforms

In this appendix we discuss the practical implementation of (12) and (13) to compute the MT coefficients and sum them via the discrete Fourier transform. We refer to these as the modified fast Fourier transform (MFFT) and modified inverse fast Fourier transform (MIFFT), respectively. First we discretize $\theta\in[-\pi,\pi]$ by $\theta_{j}=jh$, for $j=-\frac{N}{2},-\frac{N}{2}+1,\dots,\frac{N}{2}-1$, where $N$ is an even integer with angular spacing $h=\frac{2\pi}{N}$. There is no point at $j=N/2$ due to the periodic boundary conditions. The discrete MT (Fourier) coefficients are given by

$$\check{f}_{n}=(-i)^{n}\sqrt{\frac{\pi}{2}}\frac{1}{N}\sum_{j=-N/2}^{N/2-1}f\left(\frac{1}{2}\tan\left(\frac{\theta_{j}}{2}\right)\right)\left(1-i\tan\left(\frac{\theta_{j}}{2}\right)\right)e^{-in\theta_{j}},~{}~{}~{}~{}n=-\frac{N}{2},-\frac{N}{2}+1,\dots,\frac{N}{2}-1.$$ (32)

This is the discrete approximation of (11). Summing the MT coefficients yields the corresponding function

$$f_{j}=\sqrt{\frac{2}{\pi}}\frac{1}{1-i\tan\left(\frac{\theta_{j}}{2}\right)}\sum_{n=-N/2}^{N/2-1}(i)^{n}\check{f}_{n}e^{in\theta_{j}},~{}~{}~{}~{}j=-\frac{N}{2},\cdots,\frac{N}{2}-1,$$ (33)

where $f_{j}=f\left(\frac{1}{2}\tan\left(\frac{\theta_{j}}{2}\right)\right)$. This is the discrete approximation of the sum in (8). Once the MT coefficients are known, the corresponding derivative coefficients are found using $N\times N$ truncated versions of the differentiation matrices derived in Sec. 3.

Appendix C Fourth-order split-step integration

All PDEs in this work are solved using the explicit Malmquist-Takenaka split-step (MT-SS4) method described in Sec. 4.1. These approximations tends to be spectrally accurate in space and fourth-order accurate in time. The algorithm pseduo-code is given below.

The temporal interval $[t_{0},T]$ is discretized into $S+1$ equally spaced time levels $t_{l}=t_{0}+l\Delta t$ for $l=0,1,\dots,S$ with time-steps of $\Delta t=\frac{T-t_{0}}{S}$. The real line $(-\infty,\infty)$ is discretized by $x_{j+\frac{N}{2}+1}=\frac{1}{2}\tan\left(\frac{\pi j}{N}\right),$ for $j=-\frac{N}{2},-\frac{N}{2}+1,...,\frac{N}{2}-1$ and ${\bf x}=(x_{1},x_{2},...,x_{N})^{T}$ for even $N$. Denote the numerical approximation of the solution by ${\bf u}(t)=(u_{1},u_{2},...,u_{N})^{T}$ where $u_{j}\approx u(x_{j},t)$. The $N\times N$ second differentiation matrix $\mathbb{D}_{2}$ is defined in (19). The initial condition is given by ${\bf u}_{0}=u({\bf x},t_{0})$.

Consider the NLS equation (24) and splitting (27). The solution of the nonlinear equation was given in (28); notice this is solved in the physical domain. To solve the linear equation, first the solution is shifted, ${\bf u}-1$, to ensure zero boundary conditions. Next, the MFFT transform (32) is taken to find the coefficients and the linear system is solved by (29). Afterward, the function is reconstructed by MIFFT transform (33) and the unity boundary conditions are restored. The split-step method given in (23) gives the appropriate combination of steps to achieve fourth-order accuracy.

One MT-SS4 time-step for the NLS equation

$${\bf v}_{1}=\exp\left[2i\alpha_{1}\Delta t\left(|{\bf u}_{j}|^{2}-1\right)\right]{\bf u}_{j}$$

$$\check{{\bf v}}_{1}=\text{MFFT}[{\bf v}_{1}-1]$$

$$\check{{\bf w}}_{1}=\exp\left[i\beta_{1}\Delta t\mathbb{D}_{2}\right]\check{{\bf v}}_{1}$$

$${\bf y}_{1}=\text{MIFFT}\left[\check{{\bf w}}_{1}\right]+1$$

$${\bf v}_{2}=\exp\left[2i\alpha_{2}\Delta t\left(|{\bf y}_{1}|^{2}-1\right)\right]{\bf y}_{1}$$

$$\check{{\bf v}}_{2}=\text{MFFT}[{\bf v}_{2}-1]$$

$$\check{{\bf w}}_{2}=\exp\left[i\beta_{2}\Delta t\mathbb{D}_{2}\right]\check{{\bf v}}_{2}$$

$${\bf y}_{2}=\text{MIFFT}\left[\check{{\bf w}}_{2}\right]+1$$

$${\bf v}_{3}=\exp\left[2i\alpha_{3}\Delta t\left(|{\bf y}_{2}|^{2}-1\right)\right]{\bf y}_{2}$$

$$\check{{\bf v}}_{3}=\text{MFFT}[{\bf v}_{3}-1]$$

$$\check{{\bf w}}_{3}=\exp\left[i\beta_{3}\Delta t\mathbb{D}_{2}\right]\check{{\bf v}}_{3}$$

$${\bf y}_{3}=\text{MIFFT}\left[\check{{\bf w}}_{3}\right]+1$$

$${\bf u}_{j+1}=\exp\left[2i\alpha_{4}\Delta t\left(|{\bf y}_{3}|^{2}-1\right)\right]{\bf y}_{3}$$

$$t_{j+1}=t_{j}+\Delta t$$