\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \usepackage[hang]{subfigure} \usepackage{booktabs,amssymb} \usepackage{pythonhighlight} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2017 (2017), No. 62, pp. 1--23.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2017/62\hfil Wave models with SpecTraVVave] {A numerical study of nonlinear dispersive wave models with SpecTraVVave} \author[H. Kalisch, D. Moldabayev, O. Verdier \hfil EJDE-2017/62\hfilneg] {Henrik Kalisch, Daulet Moldabayev, Olivier Verdier} \address{Henrik Kalisch \newline Department of Mathematics, University of Bergen, P.O. Box 7800, 5020 Bergen, Norway} \email{henrik.kalisch@math.uib.no} \address{Daulet Moldabayev \newline Department of Mathematics, University of Bergen, P.O. Box 7800, 5020 Bergen, Norway} \email{daulet.moldabayev@math.uib.no} \address{Olivier Verdier \newline Department of Mathematics and Statistics, University of Umea, Sweden.\newline Department of Computing, Mathematics and Physics, Western Norway University of Applied Sciences, Norway} \email{olivier.verdier@hvl.no} \thanks{Submitted November 7, 2016. Published March 2, 2017.} \subjclass[2010]{35C07, 35Q53, 45J05, 65M70} \keywords{Traveling Waves; nonlinear dispersive equations; bifurcation; \hfill\break\indent solitary waves} \dedicatory{Communicated by Jerry Bona} \begin{abstract} In nonlinear dispersive evolution equations, the competing effects of nonlinearity and dispersion make a number of interesting phenomena possible. In the current work, the focus is on the numerical approximation of traveling-wave solutions of such equations. We describe our efforts to write a dedicated \textsf{Python} code which is able to compute traveling-wave solutions of nonlinear dispersive equations in a very general form. The \textsf{SpecTraVVave} code uses a continuation method coupled with a spectral projection to compute approximations of steady symmetric solutions of this equation. The code is used in a number of situations to gain an understanding of traveling-wave solutions. The first case is the Whitham equation, where numerical evidence points to the conclusion that the main bifurcation branch features three distinct points of interest, namely a turning point, a point of stability inversion, and a terminal point which corresponds to a cusped wave. The second case is the so-called modified Benjamin-Ono equation where the interaction of two solitary waves is investigated. It is found that two solitary waves may interact in such a way that the smaller wave is annihilated. The third case concerns the Benjamin equation which features two competing dispersive operators. In this case, it is found that bifurcation curves of periodic traveling-wave solutions may cross and connect high up on the branch in the nonlinear regime. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{remark}[theorem]{Remark} \newtheorem{definition}[theorem]{Definition} \allowdisplaybreaks \newcommand{\setc}[2]{\lbrace #1:#2\rbrace} \section{Introduction} \label{intro} This article concerns traveling wave solutions for a class of nonlinear dispersive equations of the form \begin{equation}\label{nonlocal} u_t + [f(u)]_x + \mathcal{L} u_x = 0, \end{equation} where $\mathcal{L}$ is a self-adjoint operator, and $f$ is a real-valued function with $f(0)=0$ and $f'(0)=0$, and which satisfies certain growth conditions. Equations of this form arise routinely in the study of wave problems in fluid mechanics and many other contexts. A prototype of such an equation is the KdV equation that appears if $\mathcal{L}= I + \frac{1}{6} \partial_x^2$ and $f(u) = \frac{3}{4} u^2$. In the current work, the operator $\mathcal{L}$ is considered to be given as a Fourier multiplier operator, such as for instance in the Benjamin--Ono equation, which arises in the study of interfacial waves. In this case, the Fourier multiplier operator is given by $\mathcal{L} = I - \mathcal{H} \partial_x$, where the Hilbert transform $\mathcal{H}$ is defined as \begin{equation} \mathcal{H} u(x) = \frac{1}{\pi} \operatorname{p.v.} \int_{-\infty}^{\infty} \frac{u(x-y)}{y} \, \mathrm{d} y, \quad \widehat{\mathcal{H} u}(k) = -\mathrm{i} \operatorname{sgn}(k)\widehat{u}(k). \end{equation} We also study in detail traveling wave solutions of the Whitham equation, which appears when $\mathcal{L}$ is given by convolution with the integral kernel $K_{h_0}$ in the form \begin{equation}\label{symbol} \mathcal{L} u(x) = \int_{-\infty}^{\infty} K_{h_0} (y) u(x-y) \, \mathrm{d} y, \quad \widehat{K_{h_0}} (k) = \sqrt{\frac{g\tanh (h_0 k)}{k}}, \end{equation} and $f$ is the same function as in the KdV equation. The particular form of equation~\eqref{nonlocal} exhibits the competing effects of dispersion and nonlinearity, which gives rise to a host of interesting phenomena. The most well known special phenomenon is the existence of solitary waves and of periodic traveling waves containing higher Fourier modes. Indeed, note that in the purely dispersive model $u_t + \mathcal{L} u_x = 0$, the only possible permanent progressive waves are simple sinusoidal waves, while in the nonlinear model~\eqref{nonlocal} higher Fourier modes must be considered to obtain solutions. The order of the operator $\mathcal{L}$ appearing in~\eqref{nonlocal} has a major effect on the types of solutions that may be found. A higher-order operator, such as in the Korteweg--de Vries equation, acts as a smoothing operator because of its effect of spreading different frequency components out due to a strongly varying phase speed \cite{Kato}. Lower-order operators such as the operator $K_{h_0}$ in \eqref{symbol} appearing in the Whitham equation may allow solutions to develop singularities, such as derivative blow-up (see \cite{Hur2015,HurTao}) and formation of cusps (see \cite{EW2016}). On the other hand, highly nonlinear functions $f(u)$ may lead to $L^{\infty}$-blow-up. For instance, the generalized KdV equation which is written in normalized form as \begin{equation} \label{kdv} u_t + u^{p}u_x + u_x + u_{xxx} = 0, \end{equation} features global existence of solutions for $p=1,2,3$, but the solutions blow-up in the critical case $p=4$ (the case $p \ge 5$ is open). In the case of the generalized Benjamin--Ono equation \begin{equation*} u_t + u^{p}u_x + u_x - \mathcal{H} u_{xx} = 0, \end{equation*} where $\mathcal{H}$ is the Hilbert transform, numerical evidence points to singularity formation for $p > 2$ \cite{BoKa}, but no proofs are available at this time. To study different phenomena related to equations of the form~\eqref{nonlocal} and their traveling wave solutions, a \textsf{Python}-based solver package \textsf{SpecTraVVave} was developed by the authors \cite{github}. The general idea behind the solver is to use a numerical continuation method \cite{Keller} implemented with a pseudo-spectral algorithm. Similar previous projects include \textsf{AUTO} \cite{Doed97} and \textsf{Wavetrain} \cite{Sherratt}. \textsf{AUTO} is written in \emph{C}, whereas \textsf{Wavetrain} is written in \emph{Fortran}. Both programs are efficient and very general, as they are able to cover a wide range of problems involving bifurcation analyses. However, from a user's perspective, such a generality coupled with low level programming languages may lead to some difficulties in utilizing these programs efficiently. \textsf{SpecTraVVave} is designed to provide researchers with a simple yet effective tool for investigating problems on traveling waves. The package is flexible, and its functionality can be easily expanded. The availability of the \textsf{IPython} notebook \cite{IPython} makes the solver very interactive, so that it should be easier for new users to get started. To maximize ease of use, \textsf{SpecTraVVave} was designed to find even solutions of~\eqref{nonlocal}. Symmetry of steady solutions can be proved for some of the models in the form~\eqref{nonlocal}, but not for all \cite{CB}. Some of these equations also admit non-smooth solutions, for instance as termination points of a bifurcation branch. This happens for example for the Whitham equation, which features bifurcation curves which terminate in a solution with a cusp \cite{EW2016}. One of the goals of the present paper is to investigate the precise nature of the termination of the bifurcation curve. The content of this article is structured as follows. A mathematical description of the numerical method of \textsf{SpecTraVVave} is given in Section \ref{sec:spectral}. Section \ref{sec:experiments} presents results of different experiments carried out with the package. Concluding remarks are given in Section \ref{sec:conclusion}. A method for finding initial guesses for the solver is described in Section \ref{sec:stokes}. Section \ref{sec:implementation} contains a schematic of program and a description of its workflow. \section{Spectral scheme and construction of nonlinear system.} \label{sec:spectral} \subsection{Cosine collocation method.} To compute traveling wave solutions to \eqref{nonlocal} the following ansatz is employed: \begin{equation*}%\label{ansatz} u(x,t) = \phi(x-ct). \end{equation*} Thus, the equation takes the form \begin{equation*}%\label{phi-nonlocal} \phi' + \left[f(\phi) \right]' + \mathcal{L} \phi' = 0, \end{equation*} which can be integrated to give \begin{equation}\label{int_gen1} -c \phi + f(\phi) + \mathcal{L} \phi = B. \end{equation} The constant $B$ is a priori undetermined. One may set the $B$ equal to zero as a way of normalizing the solutions. Another option is to impose an additional condition, for example that the integral of $\phi$ over one wavelength be zero. In this case, $B$ will be found along with the solution $\phi$. We consider $\mathcal{L}$ as a Fourier multiplier operator with symbol $\alpha(k)$. We also assume that $f$ is at least twice differentiable, and we have $f(0)=0$ and $f'(0)=0$. When computing traveling-wave solutions we focus on even periodic solutions. While it can be proved in some cases that solutions of~\eqref{int_gen1} must be even, this is not known for a general operator $\mathcal{L}$. Nevertheless, we make this assumption here in order to make the numerical procedure as uniform as possible. For even periodic solutions, one may use a cosine collocation instead of a Fourier method. In particular, using the cosine functions as basis elements automatically removes the inherent symmetries due to reflective and translational symmetry. Moreover, the number of unknowns is reduced by a factor of $2$, and the problem of the asymmetric arrangement of nodes in the FFT is circumvented. Of course, all these problems could also be dealt with a collocation method based on the Fourier basis, but the cosine basis does all of the above automatically. In addition, the \textsf{Python} cosine transform is based on an integrated algorithm, which relies on an optimized version of the discrete cosine transform (DCT). The following description of computation scheme was presented in detail in \cite{EK2}, but we will briefly repeat it here for consistency of the manuscript. For the purpose of clarity, we will refer to full wavelength $L$ of a solution as the (full) wavelength, and half of fundamental wavelength will be called half-wavelength. Such a definition is required because the method computes a half of a solution profile, the other half is automatically constructed due to symmetry. Traveling wave solutions to \eqref{int_gen1} are to be computed in the form of a linear combination of cosine functions of different wave-numbers, i.e., in the space \begin{equation} \mathcal{S}_N = \operatorname{span}{}_{\mathbb{R}}\setc{\cos(lx)}{0 \leq l \leq N-1 }. \end{equation} This is a subspace of $L^2(0,2\pi)$, and the collocation points $x_n = \pi \frac{2n -1}{2N}$ for $n = 1,\dots, N$ are used to discretize the domain. If the required full wavelength of solutions is to be $L \neq 2\pi$, one can use a scaling on the $x$-variable. Defining the new variable \begin{equation*} x' = \frac{L}{2\pi}x, \end{equation*} yields collocation points $x_n'$ and wavenumbers $\kappa_l$ defined by \begin{equation*} x_n' = \frac{L}{2} \frac{2n -1}{2N}, \quad\quad \kappa_l = \frac{2\pi}{L}l. \end{equation*} We are seeking a function $\phi_N \in \mathcal{S}_N$ that satisfies the equations \begin{equation} -c\phi_N(x_n') + f(\phi_N)(x_n') + \mathcal{L}^N\phi_N(x_n') = 0, \label{int_gen-disc} \end{equation} at the collocation points $x_n'$. The operator $\mathcal{L}^N$ is the discrete form of the operator $\mathcal{L}$, and $\phi_N$ is the discrete cosine representation of $\phi$ which is given by \begin{gather*} \phi_N(x') = \sum_{l = 0}^{N-1} \omega(\kappa_l)\Phi_N(\kappa_l)\cos(\kappa_{l}x'),\\ \omega(\kappa_l) = \begin{cases} \sqrt{1/N}, \quad \kappa_l = 0, \\ \sqrt{2/N}, \quad \kappa_l > 0, \end{cases} \\ \Phi_N(\kappa_l) = \omega(\kappa_l) \sum_{n = 1}^{N} \phi_N(x'_n) \cos(\kappa_l x'_n), \end{gather*} where $\kappa_l = 0,\frac{2\pi}{L},\ldots, \frac{2\pi}{L}(N-1)$ are the scaled wavenumbers, and $\Phi_N(\cdot)$ are the discrete cosine coefficients. As equation \eqref{int_gen-disc} is enforced at the collocation points $x_n'$, one may evaluate the term $\mathcal{L}^N\phi_N$ using the matrix $\mathcal{L}^N(i,j)$ defined by \begin{gather*} \mathcal{L}^N\phi_N(x'_i) = \sum_{j=1}^{N}\mathcal{L}^N(i,j) \phi_N(x'_j),\\ \mathcal{L}^N(i,j) = \sum_{l=0}^{N-1} \omega^2(\kappa_l) \alpha(\kappa_l) \cos(\kappa_l x'_i) \cos(\kappa_l x'_j), \end{gather*} where $\alpha(\cdot)$ is the Fourier multiplier function of the operator $\mathcal{L}$. \subsection{Construction of nonlinear system.} Equation \eqref{int_gen-disc} enforced at $N$ collocation points yields a nonlinear system of $N$ equations in $N$ unknowns, which can be written in shorthand as \begin{equation*} F(\phi_N) = 0. \end{equation*} This system can be solved by a standard iterative method, such as Newton's method. In this system, the value of phase speed $c$ has to be fixed for computing one particular solution. Such an approach becomes impractical when a turning point on the bifurcation curve appears. In \textsf{SpecTraVVave} a different approach is employed: both the amplitude $a$ and the phase speed $c$ of a solution are treated as functions of a parameter $\theta$: $ a = a(\theta)$, $c = c(\theta)$. The parameter $\theta$ is to be computed from the system~\eqref{system-N3}. This construction makes it possible to follow turning points on the bifurcation branch with relative ease. Having computed two solutions, i.e., two points on the bifurcation curve $P_1 = (c_1, a_1)$ and $P_2 = (c_2, a_2)$, one may find a direction vector $\mathbf{d}=(d^c, d^a)$ of the line that contains these points: \[ \mathbf{d}: \quad d^c = c_2 - c_1, \quad d^a = a_2-a_1. \] Then the point $P_3 = (c_3,a_3)$ is fixed at some (small) distance $s$ from the point $P_2$ in the direction $\mathbf{d}$. \[ P_3: \quad c_3 = c_2 + s \cdot d^c, \quad a_3 = a_2 + s \cdot d^a. \] The point $P_3$ plays the role of the initial guess for velocity and amplitude when computing the next solution $P_\ast = (c_\ast, a_\ast)$. The solution point $P_\ast$ is required to lay on the line with direction vector $\mathbf{d_{\bot}}=(d^c_{\bot}, d^a_{\bot})$, which is orthogonal to the vector $\mathbf{d}$. \begin{gather*} \mathbf{d_{\bot}}: \quad d^c_{\bot} = - d^a, \quad d^a_{\bot} = d^c, \\ P_\ast: \quad c_\ast = c_3 + \theta d^c_{\bot}, \quad a_\ast = a_3 + \theta d^a_{\bot}. \end{gather*} A schematic sketch of finding a new solution $P_\ast$ is given in Figure \ref{navig}. The variable $\theta$ is computed by Newton's method from the extended system \begin{equation} F\begin{pmatrix} \phi_N(x_1)\\ \vdots \\ \phi_N(x_N) \\ B \\ \theta \end{pmatrix} = \begin{pmatrix} (-c+\mathcal{L}_N)\phi_N(x_1)+ f(\phi_N(x_1)) - B \\ \vdots \\ (-c+\mathcal{L}_N)\phi_N(x_N)+ f(\phi_N(x_N)) - B \\ \Omega(\phi_N, c, a, B) \\ \phi_N(x_1)-\phi_N(x_N)-a \end{pmatrix} =\begin{pmatrix} 0 \\ \vdots \\ 0 \\ 0 \\ 0 \end{pmatrix}. \label{system-N3} \end{equation} Here a nonhomogeneous problem ($B \neq 0$) is considered. The equation \begin{equation*} \phi_N(x_1)-\phi_N(x_N)-a = 0, \end{equation*} makes the waveheight of the computed solution to be that of $a$. The equation \begin{equation*} \Omega(\phi_N, c, a, B)=0, \end{equation*} is called the \emph{boundary condition}. It allows to enforce different specifications on the computed traveling wave solution. For example, if one sets \begin{equation*} \Omega(\phi_N, c, a, B)= \phi_N(x_1)+\dots+\phi_N(x_N), \end{equation*} then \emph{the mean of the computed wave over a period will have to be equal to zero}. One may also experiment with \[ \Omega(\phi_N, c, a, B) = B, \] to consider the homogeneous problem $(B=0)$. It can be also interesting to set \begin{equation} \Omega(\phi_N, c, a, B) = \phi_N(x_N). \label{solitary} \end{equation} This enables us to compute traveling wave solutions that mimic solitary wave solutions. \begin{figure}[ht] \begin{center} \includegraphics[width=0.35\textwidth]{fig1} % navig \end{center} \caption{Navigation on the bifurcation curve.} \label{navig} \end{figure} \subsection{Convergence} To test the numerical implementation of the discretization, the method is applied to a case where the solution is known. One such case is the KdV equation \[ u_t + u_x + \frac{3}{2} u u_x + \frac{1}{6}u_{xxx} = 0, \] which has a known solution, given in the form \[ u_{\mathrm{exact}}(x,t) = a \operatorname{sech}^2\Big( \sqrt{\frac{3a}{4}} (x-c t)\Big), \] with $c=1+a/2$. Using the boundary equation~\eqref{solitary}, \textsf{SpecTraVVave} is capable of computing approximations to solitary wave solutions of nonlinear wave equations. Solitary wave solutions are treated as traveling waves with sufficiently long wavelength that have the wave trough at zero. In case of the KdV equation solitary wave solutions have exponential decay, and therefore, considering the symmetry of solitary solutions, the half-wavelength of $30$ is considered for the comparison. Approximation errors are summarized in Table \ref{table1}. \begin{table}[ht] \begin{center} \begin{tabular}{cccc} \toprule grid points & $\log_{10}(\|u_{\mathrm{exact}} - u\|_{L^{\infty}})$ & $\log_{10}(\|u_{\mathrm{exact}} - u\|_{L^2})$ & Ratio of $L^2$-errors\\ \midrule 32 & $-1.389$ & $-2.092$ & \\ 64 & $-3.705$ & $-4.549$ & $286.8$\\ 128 & $-8.809$ & $-9.508$ & $90935.0$\\ 256 & $-15.353$ & $-16.144$ & $4329670.9$\\ 512 & $-15.353$ & $-16.087$ & $0.9$\\ \bottomrule \end{tabular} \end{center} \caption{Estimates of error between the exact and computed solitary wave solutions for the KdV equation. Half-wavelength $30$, waveheight $a=1.2651$} \label{table1} \end{table} \begin{figure}[ht] \begin{center} \includegraphics[width=0.47\textwidth]{fig2} % convergence \end{center} \caption{Graph of error estimates given in Table \ref{table1}.} \label{L2convergence} \end{figure} \section{Experiments with SpecTraVVave.} \label{sec:experiments} \subsection{Termination of waveheight-velocity bifurcation curve of the Whitham equation} \label{sec-W-term} The waveheight-velocity bifurcation curve of the Whitham equation \begin{align} u_t + \frac{3}{2} u u_x + Ku_x = 0, \quad \widehat{Ku}(k) = \sqrt{\frac{\tanh(k)}{k}}, \label{W} \end{align} was studied numerically in \cite{EK2}. An attempt was made to identify the termination point of the Whitham bifurcation curve. The investigation was limited by computational tools and complete results were not obtained. In particular, the authors could not confirm that traveling wave solutions do not exist past the point where the authors, based on pioneering work of Whitham \cite{Wh1} suspected a cusped solution. In this section a number of numerical results on nature of the bifurcation curve for the Whitham equation are presented. Solutions to \eqref{W} are computed in the form of traveling waves $u(x,t) = \phi(x - c t)$ and the homogeneous $(B = 0)$ integrated version the equation is considered: \begin{equation} -c \phi + \frac{3}{4} \phi^2 + K \phi = 0. \label{W-int} \end{equation} Special attention is given to relation between stability of solutions and their waveheight and velocity parameters, i.e., their position on the bifurcation curve. The following questions are under study: \begin{enumerate} % chktex 9 chktex 10 \item[(a)] Where does the bifurcation curve terminate? \item[(b)] Where on the bifurcation curve do solutions change their stability? \item[(c)] Is there any role that the turning point on the bifurcation curve plays? \end{enumerate} The results presented here focus on $2\pi$-periodic solutions to \eqref{W-int}, i.e., solutions of system~\eqref{system-N3}. Figure \ref{bif-curve-all} presents Whitham bifurcation curves with numbers of grid points $N = 512$, $N=1024$ and $N=2048$. The current implementation of the \textsf{SpecTraVVave} package allows fixing the number of grid points $N$ and a so-called doubling parameter $\mathcal{D}$, i.e., the number by which $N$ is doubled as computations are made. This allows us to get sets of solutions with $N, 2N, \ldots, 2^{\mathcal{D}} N$ grid points. If $\mathcal{D} = 1$ then only two sets of solutions are computed and they are regarded as lower grid (lower resolution) and higher grid (higher resolution) solutions. While system \eqref{system-N3} is processed by Newton solver, lower grid solutions are taken as initial guesses for higher grid solutions. All curves shown in this manuscript have been produced after tests with a number of resolutions were run, and the curves shown did not change significantly under further refinement. \begin{figure}[ht] \begin{center} \includegraphics[width=0.48\textwidth]{fig3} % 512_1024_2048 \end{center} \caption{Whitham bifurcation curve in different grid resolutions.} \label{bif-curve-all} \end{figure} Figure \ref{bif-L2-1024}\subref{sfig:bif-ch} presents the Whitham bifurcation curve computed by \textsf{SpecTraVVave} with $N=1024$ and $\mathcal{D}=1$. There are three solutions which deserve to be singled out: \begin{enumerate} \item Traveling wave solution with minimum velocity (rhombus); \item Traveling wave solution with maximum $L^2$-norm (circle); \item Cusped traveling wave solution (square). \end{enumerate} Profiles of the above listed solutions are given in Figure \ref{3waves}\subref{sfig:prof-sing}. The solution with minimum velocity corresponds to the turning point of the bifurcation curve. The solution with maximum $L^2$-norm is very close to the latter one, although it has a higher waveheight and a different velocity. The solution marked by a square is called here the \emph{terminal solution}. As already mentioned, previous studies, such as \cite{EK1,EK2} did not provide any conclusive analysis on the part of the bifurcation curve past the turning point. In particular, it was not clear whether solutions ceased to exist at or after the turning point, or whether solutions were stable or unstable after the turning point. Let us first focus on the stability of solutions. Note that \textsf{SpecTraVVave} has an evolution integrator routine, which enables one to check the stability of computed solutions. The current version of the package uses the fourth-order method developed in \cite{dFSS}. In addition one may use a more refined analysis, resting on the evaluation of invariant functionals. This analysis is based on the observation that the traveling waves can be thought of as solutions of a constrained minimization problem. This analysis is based on ideas developed by Boussinesq, first exploited in \cite{B2}, and later used in \cite{BSS,TJB,NgK}, and many other works. \begin{figure}[ht] \begin{center} \subfigure[][Bifurcation curve in the wavespeed-wave\-height parameter space]{\includegraphics[width=0.48\textwidth]{fig4a} % 1024_2048 {}\label{sfig:bif-ch}} \subfigure[][Bifurcation curve in a wavespeed-$L^2$-norm diagram] {\includegraphics[width=0.48\textwidth]{fig4b} % L2_1024_2048 {}} \end{center} \caption{Whitham bifurcation curves for $2\pi$-periodic solutions.} \label{bif-L2-1024} \end{figure} Let us define two functionals $V$ and $E$: \[ V( \phi) = \frac{1}{2} \int_{-\infty}^{+\infty} \phi^2\,\mathrm{d}\zeta, \quad E( \phi) = \int_{-\infty}^{+\infty} \left\{ \frac{1}{2} \phi^3 - \phi \, K\! \phi \right\}\, \mathrm{d}\zeta. \] Equation \eqref{W-int} can be then written in terms of variational derivatives of $E$ and $V$ as \begin{equation} E' ( \phi ) - c V ' ( \phi ) = 0.\label{zero} \end{equation} It is known from \cite{BSS} that under certain conditions, the stability of solitary wave solutions depends on convexity of the function $d(c) = E(\phi) - c V(\phi)$. Solutions with values of $c$ for which $d''(c) > 0$ are stable solutions, and solutions with wave speeds for which $d''(c) <0$ are unstable solutions. The current numerical investigation can therefore be interpreted as an indication of the stability properties of traveling wave solutions of the Whitham equation. Note that differentiation of $d(c)$ yields \[ d '(c) = \underbrace{E'(\phi) - c V'(\phi)}_{=0} - V(\phi). \] Using~\eqref{zero} as indicated yields \[ d '(c) = - V(\phi) = - \frac{1}{2} \int_{-\infty}^{+\infty} \phi ^2\,\mathrm{d}\zeta = - \frac{1}{2}\|\phi\|_{L^2}^2 . %\label{d'c} \] Therefore, to understand the convexity of $d(c)$, it is sufficient to find points of maximum $L^2$-norm on the curve in the right panel of Figure \ref{bif-L2-1024}. It is straightforward to see that $d''(c)$ changes sign in the neighbourhood of the maximum point of this curve, i.e., around the solution with maximum $L^2$-norm. In particular, $d''(c) > 0$, i.e., solutions are stable to the left of the maximum point, and $d''(c) < 0$, i.e., solutions are unstable to the right of the maximum point. In addition, the solutions were tested with the evolution integrator to confirm their stability/instability in time. The solution with maximum $L^2$-norm and those on the left-hand side were always found to be stable in the time-dependent computations. Solutions on the right-hand side do not preserve their shape and thus are unstable. Examples are given in Figure \ref{time-evo}. This analysis confirms that the point corresponding to the minumum wave speed (the turning point), and the point of stability inversion are two distinct points on the bifurcation curve. Moreover, the point of stability inversion is a little further up the branch from the turning point. \begin{figure}[ht] \begin{center} \subfigure[][Profiles of the solutions singled out in Figure \ref{bif-L2-1024} \subref{sfig:bif-ch}] {\includegraphics[width=0.48\textwidth]{fig5a} % 1024_waves_1 {}\label{sfig:prof-sing}} \subfigure[][Profile after terminal solution] {\includegraphics[width=0.48\textwidth]{fig5b} % 1024_waves_2 {}\label{sfig:prof-term}} \end{center} \caption{Profiles of specific traveling wave solutions.} \label{3waves} \end{figure} Next, we turn our attention to the analysis of the terminal point. There are two main questions. Does the branch terminate, and if so, does the terminal point on the branch correspond to a cusped traveling wave. First of all, note that the solution, which is computed by \textsf{SpecTraVVave}, past the terminal solution has two crests, no matter how small the stepping on the bifurcation branch is taken. (see Figure \ref{3waves}\subref{sfig:prof-term}). Secondly, as will be explained presently, the relation \[ \frac{c}{\sup_{x} \phi(x)} = \frac{3}{2} \] holds for the terminal solution with a good degree of approximation. For the most accurate runs, we obtain $c / \sup_{x \in \mathbb{R}} \phi(x)$ $\approx$ $1.51$. To explain how this relation comes about note that the steady integrated form of the Whitham equation can be written as \begin{equation}\label{Whitham_Alt} \Big( \frac{c}{\sqrt{3}} - \frac{\sqrt{3}}{2}\phi \Big)^{2} = \frac{1}{3} c^2 - K \phi. \end{equation} It is clear that for any $\phi < 3c/2$, the relation~\eqref{Whitham_Alt} can be used in a bootstrap argument to show that any continuous solution must be in fact smooth. However for the case $\phi = 3c/2$ this bootstrap argument fails since the left-hand side vanishes. It can be concluded that a solutions containing a cusp will have a maximum value of $3c/2$. As an additional check, the discrete cosine coefficients of the solutions were examined, and fitted to the following models: \[ \mathcal{E}(k) = \nu_1 e^{-\nu_2 k^n}, \quad \mathcal{P}(k) = \frac{\mu_1}{\mu_2 + \mu_3 k^m}, \] where $\nu_1$, $\nu_2$, $\mu_1$, $\mu_2$, $\mu_3$, $n$ and $m$ are fitting parameters. A smooth function is known to have discrete cosine coefficients with exponential decay in $k$. On the other hand, if a function is not smooth, the discrete cosine coefficients feature polynomial decay. To identify the best fit, two parameters were used: $L^2$-norm of the residual and the Akaike information criterion (AIC) measure. From the data given in the Table \ref{table2}, one can deduce that for solutions with minimum speed and maximum $L^2$-norm exponential fit is better than polynomial. That is not the case for the terminal solution. Thus, the first two solutions are smooth and the terminal solution is nonsmooth. In fact, the polynomial fit is better than exponential for solutions that are between the maximum $L^2$-norm solution and the terminal solution. \begin{figure}[ht] \begin{center} \subfigure[][Minimum speed solution]{ \includegraphics[width=0.30\textwidth]{fig6a} % 1024_full_3p_1 } \subfigure[][Maximum $L^2$-norm solution]{ \includegraphics[width=0.30\textwidth]{fig6b} % 1024_full_3p_2 } \subfigure[][Terminal solution]{ \includegraphics[width=0.30\textwidth]{fig6c} % 1024_full_3p_3 } \end{center} \caption{Evolution of specific solutions in time. Time range for each solution is three periods.} \label{time-evo} \end{figure} The numerical evidence brought forward supports the conclusion that the Whitham bifurcation branch terminates at the terminal point indicated in Figure \ref{bif-curve-all}. Of course, as already mentioned, this conclusion has now also been reached using tools of mathematical analysis \cite{EW2016}. \begin{table}[ht] \begin{center} \begin{tabular}{l|rr|rr|rr} \cline {1-7} & \multicolumn{2}{ c| }{Min speed solution} & \multicolumn{2}{ c |} {Max $L^2$-norm solution}& \multicolumn{2}{ c }{Terminal solution} \\ \multicolumn{1}{l|} {Model} &$\mathcal{E}(k)$&$\mathcal{P}(k)$ &$\mathcal{E}(k)$&$\mathcal{P}(k)$&$\mathcal{E}(k)$&$\mathcal{P}(k)$ \\ \midrule \multicolumn{1}{l|} {\small Res.~$L^2$-norm} & $5\times 10^{-5}$ &$4\times 10^{-3}$ & $7 \times 10^{-5}$ & $3\times 10^{-3}$ & $6\times 10^{-3}$& $6\times 10^{-4}$ \\ \multicolumn{1}{l|} {AIC} & -543 & -321 & -529 & -333 & -298 & -416 \\ \bottomrule \end{tabular} \end{center} \caption{Results for measures of fit.} \label{table2} \end{table} \subsection{Interaction of solitary wave solutions of modified Benjamin--Ono equation} \label{mB-O-interact} In this section, we utilize the \textsf{SpecTraVVave} package to obtain high-precision approximations to solitary-wave solutions of the modified Benjamin--Ono \[ u_t + u^2u_x + u_x - \mathcal{H} u_{xx} = 0, \] which is a special case of the generalized Benjamin--Ono equation, with $p=2$. This case corresponds to the critical scaling, i.e., invariance of the energy norm under the natural invariant scaling. The Benjamin--Ono equation was found by Benjamin \cite{B1} as a model for long small-amplitude interfacial waves in deep water. The validity of approximating the more physically correct configuration of a continuous density distribution by the two-layer approximation has recently been justified mathematically \cite{ChenWalsh2016}. Solitary-wave solutions of the modified Benjamin--Ono equation with $p=3$, $p=4$ and $p=5$ were approximated in \cite{BoKa} with a standard Newton scheme. The solutions in \cite{BoKa} were not very accurate, but since singularity formation of the evolution equations were under study, the accuracy of the solitary-wave approximation was not an important issue. The problem with the method of \cite{BoKa} and some other works was that the fft used there was not purged of possible symmetries (translational and reflective). In the current code, since a cosine formulation is chosen, these symmetries are automatically eliminated, and the resulting computations are able to to render more accurate approximations. Solitary-wave solutions of these equation could be computed with higher accuracy using a type of Petviashvili method in \cite{Pelinovsky}, but here we employ the \textsf{SpecTraVVave} package using the boundary equation~\eqref{solitary}, and treating solitary waves as traveling waves with sufficiently long wavelength that have wave trough at zero. Once these high-accuracy solutions are found, they are aligned in an evolution code using a high-order time integrator, and the interaction of two waves is studied. Two questions are investigated. First, the interaction is investigated for evidence of integrability. Second, we are looking for possible annihilation of one of the waves, such as may happen in some other evolution equations \cite{Tjon}. A possible approach to studying the question of complete integrability is analyzing the interaction of two solitary wave solutions of the equation, such as carried out in \cite{HK,HKJB} for other nonlocal equations. In Figure \ref{solitons-full}, snapshots of interaction of two solitary waves at different times are shown. The time difference between two consecutive snapshots is constant. As it may be observed, during the process of interaction, the two initial solitary waves combine into a single wave, and an additional oscillation is produced. This leads us to the conclusion that the interaction of solitary waves is not elastic and the modified Benjamin--Ono equation may not be integrable. In addition, it appears that the smaller wave disappears as most of its mass is acquired by the larger wave. Thus one may argue that the small wave is annihilated by the larger wave. It can also be observed that the larger wave starts growing, and it is likely that this growth will lead to finite-time blow-up. This question was not investigated further since blow-up phenomena have already been studied closely in \cite{HKJB}. Indeed, very recently, the finite-time blow-up has been proven mathematically in \cite{MP}. \begin{figure}[ht] \begin{center} \includegraphics[width=0.5\textwidth]{fig7} % soliton_inter \end{center} \caption{Interaction of two solitary waves of the modified Benjamin-Ono equation.} \label{solitons-full} \end{figure} \subsection{Effect of competing dispersion in the Benjamin equation} \label{Benjamin_eq} The Benjamin equation was found by Benjamin \cite{B3} as a model for two-layer flow in the case when the interface is subject to surface tension. The approximation may not be a good model for a stratified situation, but more applicable to the case where two fluids are separated by a sharp interface. The equation is \begin{equation}\label{benjamin} u_t + u_x + uu_x - \mathcal{H} u_{xx} - \tau u_{xxx} = 0, \end{equation} where $\tau$ is a parameter similar to the inverse of the Bond number in free surface flow \cite{B3,KalischBenjamin,Walsh2014a}. In this section, a study relating to the effects of competing dispersion operators on the shape of periodic traveling waves in the Benjamin equation is presented. An in-depth study of solitary waves was carried out in \cite{DDM}. As will come to light, the periodic case features some new phenomena, such as secondary bifurcations, connecting and crossing branches. For the purpose of this study, we fix the parameter $\tau=0.1$, so that the dispersion relation for the linearized equation is \begin{equation} c(k) = 1 - |k| + 0.1 k^2. \label{Ben_disp_rel} \end{equation} Traveling wave solutions with full wavelengths $L_1 = \pi/5$, $L_2 = 4\pi /19$ and $L_3 = 4\pi$ are computed for \eqref{benjamin}. The corresponding wavenumbers are $k_1 = 10$, $k_2 = 19/2$, and $k_3 = 0.5$, respectively. A plot of the dispersion relation~\eqref{Ben_disp_rel} is given in Figure \ref{disp-rel}. Bifurcation branches of traveling wave solutions with the selected wavelengths are given in Figure \ref{fig:ben-bifur}. \begin{figure}[ht] \begin{center} \includegraphics[width=0.47\textwidth]{fig8} % disp_rel_plot \end{center} \caption{Dispersion relation \eqref{Ben_disp_rel}} \label{disp-rel} \end{figure} \begin{figure} \begin{center} \includegraphics[width=0.47\textwidth]{fig9} % Benj_bifur_plot \end{center} \caption{Bifurcation curves of \eqref{benjamin} with different wavelengths, $\blacktriangle$ -- points of bifurcation, $\bullet$ -- selected solutions (see Figure \ref{long-solutions}).} \label{fig:ben-bifur} \end{figure} The branch denoted by $L_1$ originates at the bifurcation point located at $c = 1$ and zero waveheight. The branches denoted by $L_2$ and $L_3$ originate from the same bifurcation point, located at $c = 0.525$ and zero waveheight. These two branches continue in different directions, due to differences in wavelength. In particular, the $L_3$ branch contains waves with shorter wavelengths, and falls into the capillary regime. On the other hand, the $L_2$ branch falls in the gravity regime. As the waveheight grows, solutions on the $L_3$ branch first cross the $L_1$ branch without connecting. Additional oscillations develop in the solutions until a new fundamental wavelength $\pi/5$ is reached, and the branch terminates as it connects to the $L_1$ branch. The situation is depicted in Figure \ref{long-solutions}. The point where the $L_1$ and $L_3$ branches meet is approximately $(c^*,a^*) = (0.945, 5.938)$. The corresponding profiles essentially overlap, as shown on Figure \ref{last-solutions}. This point can also be interpreted as a secondary bifurcation point of the $L_1$ branch, where solutions with wavelengths that are multiples of $\pi/5$ develop. We should note that similar phenomena concerning crossing and connecting branches were previously observed in \cite{RK} for the Whitham equation with surface tension which was introduced in \cite{HurJohnson}. \begin{figure}[ht] \begin{center} \begin{subfigure}[] {\includegraphics[width=0.47\textwidth]{fig10a} % left_long_solutions }\end{subfigure} \begin{subfigure}[] {\includegraphics[width=0.47\textwidth]{fig10b} % right_long_solutions }\end{subfigure} \end{center} \caption{Selected solutions of \eqref{benjamin}. Labels preserved as shown in Figure \ref{fig:ben-bifur}.} \label{long-solutions} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.47\textwidth]{fig11} % last2 \end{center} \caption{Solution profiles at the point $(c^*, a^*)$ where the $L_1$ and $L_3$ branches meet (see Figure \ref{fig:ben-bifur}).} \label{last-solutions} \end{figure} \section{Conclusions and future work} \label{sec:conclusion} The numerical algorithm of \textsf{SpecTraVVave} features ample flexibility for researching different aspects of nonlocal dispersive wave equations and their traveling wave solutions. The solver package is simpler in use when compared with programs such as \textsf{AUTO} and \textsf{Wavetrain}, however it does not have the same level of generality. Moreover, \textsf{AUTO} and \textsf{Wavetrain} are programmed in low-level programming languages and will therefore run more efficiently. \textsf{SpecTraVVave} is implemented in an object-oriented fashion \cite{PyBook}, which makes the program easily expandable. \textsf{IPython} provides means for interactive work with the package, and enables users to create convenient notebook-programs. A parametric approach in defining amplitude and phase speed makes it possible to follow turning points on bifurcation curves. Specification of different boundary conditions allows computing solutions with certain features, such as traveling waves with zero mean, or approximations to solitary waves. In this work, the \textsf{SpecTraVVave} package has been put to use for the study of a number on nonlinear evolution equations: the Whitham equation, the modified Benjamin--Ono equation and the Benjamin equation. For the chosen set of parameters, experiments on the Whitham equation resulted in numerical confirmation of the conjecture on cusped solutions. It was also possible to identify the point of stability inversion of traveling wave solutions of the equation and the termination point of its bifurcation curve. In case of the modified Benjamin--Ono equation, the study on solitary wave solutions lead us to conclude that interaction process ended with annihilation of one of the two waves. The experiment on the Benjamin equation showed one more example of the effect of competing dispersion. As the amplitude increased, traveling wave solutions of wavelength $4\pi$ developed additional oscillations, and later connected up with a branch of solutions with wavelength $\pi/5$. Future work on the \textsf{SpecTraVVave} package will be focused on development of its functionality and broadening the range of problems that can be studied. Possible extensions may include implementation of algorithms based on the Petviashvili method \cite{Duran2014CNSNS,Duran2014JCAM} and generalization to systems of equations. \section{Appendix: Computing initial guesses from Stokes expansion.} \label{sec:stokes} \begin{equation} u_t+\left[f(u)\right]_{x} + \mathcal{L}u_x=0, \label{gen} \end{equation} The goal of this section is to explain how the idea of Stokes's approximation works in providing the initial data (guess) on wave and phase velocity for solving \eqref{gen} numerically. We will consider $\mathcal{L}$ being linear and self-adjoint Fourier multiplier operator, and a function $f$ that has degree of zeros $p \geq 2$: \begin{gather*} \widehat{\mathcal{L} u}(k) = \alpha(k) \cdot \widehat{u}(k),\\ \langle\mathcal{L} u, v \rangle_{L^2(0,L)} = \langle u, \mathcal{L} v\rangle_{L^2(0,L)}, \\ p = \min_{k \in \mathbb{N}}{f^{(k-1)}(0) =0 \text{ and } f^{(k)}(0)\neq 0} \end{gather*} Consider equation \eqref{gen} and its solution in the form $u(x,t) = \phi(x-ct)$, which is a traveling wave solution. Inserting $\phi(x-ct)$ into~\eqref{gen} leads to the equation \[ -c\phi' + f'(\phi)\phi'+\mathcal{L}\phi'= 0, \] which can be integrated to give \begin{equation} -c\phi + f(\phi)+\mathcal{L}\phi= B, \quad B = \mathrm{const}. \label{int_gen-B} \end{equation} Consider $B = 0$ in equation~\eqref{int_gen-B}, and expansions of $\phi$ and $c$: \begin{gather} \phi = \xi = \varepsilon\xi_1 + \varepsilon^2\xi_2+ \dots, \label{fi_exp}\\ c = c_0 + \varepsilon c_1 + \varepsilon^2 c_2+\dots. \label{c_exp} \end{gather} The next step is to insert~\eqref{fi_exp} and~\eqref{c_exp} to \eqref{int_gen-B} and write out the terms at powers of $\varepsilon$. The function $f(\phi)$ is expanded around zero and, therefore, will appear only in $\varepsilon^p$ terms. Thus, the term at the first power of $\varepsilon$ reads \begin{equation} \varepsilon: \quad -c_0\xi_1 + \mathcal{L}\xi_1 = 0, \label{epsilon1} \end{equation} Hence, $c_0$ is an eigenvalue of the operator $\mathcal{L}$, regarded as defined on $L$-periodic functions. Taking the Fourier transform of \eqref{epsilon1} gives: \begin{equation} -c_0\widehat{\xi_1}(k) + \alpha(k)\widehat{\xi_1}(k) = 0 . \label{eps1-transformed} \end{equation} Equation~\eqref{eps1-transformed} has two trivial solutions: either $\xi_1(k) \equiv 0$ or $\alpha(k) \equiv c_0$. If we assume non-trivial $\xi_1$ and $\alpha(k) \neq \text{const}$, the following solves the problem \begin{equation} \widehat{\xi_1}(k) = 2\pi\delta(k - k_0), \quad \text{and} \quad c_0 = \alpha(k_0), \label{solution-xi} \end{equation} for some $k_0 \in \mathbb{R}$. Since $\xi_1$ is the first-order approximation to $\phi$, the corresponding wave number should be equal to 1. The $L$-periodicity condition entails that $k_0 = 2\pi/L\cdot 1$. The spatial variable $x$ has to be scaled to $x'= L/2\pi x$, accordingly. From the solutions in~\eqref{solution-xi} we have \begin{equation} \xi_1(x') = e^{ik_0x'} = \cos(k_0x') + \sin(k_0x'). \end{equation} Considering the projection onto the space $\mathcal{S}_N$, we are led to choose $\xi_1(x') = \cos(k_0x')$. For further analysis, let us define an operator $\mathcal{A}$, \[ \mathcal{A} := - c_0 \mathcal{E} + \mathcal{L}, \] where $\mathcal{E}$ is the identity operator. The operator $\mathcal{A}$ inherits the property of being self-adjoint from $\mathcal{L}$. Moreover, it follows from~\eqref{epsilon1} that $\mathcal{A} \xi_1 = 0$ and $\xi_1 \in \ker(\mathcal{A})$. If $p>2$ then $f''(0)=0$ and the terms at $\epsilon^2$ are: \begin{equation} \mathcal{A} \xi_2 - c_1\xi_1 = 0. \end{equation} Taking scalar multiplication of the latter with $\xi_1$, one obtains \begin{gather*} \langle\xi_1, \mathcal{A}\xi_2 \rangle_{L^2(0,L)} = c_1 \| \xi_1\|^2_{L^2(0,L)}, \quad \langle\xi_1, \mathcal{A}\xi_2\rangle_{L^2(0,L)} = \langle\mathcal{A} \xi_1, \\ \xi_2\rangle_{L^2(0,L)} = \langle0, \xi_2\rangle_{L^2(0,L)} = 0. \end{gather*} As a result, one has $c_1 \| \xi_1\|^2_{L^2(0,L)} = 0$ and, hence, $c_1 = 0$. Repeating the same argument, it becomes clear that $c_k = 0$ for any $k \leq p-1$. Besides that, $\xi_2$ is in the kernel of $\mathcal{A}$, so it may be assumed to be proportional to $\xi_1$. The terms at order $\varepsilon^p$ are: \begin{align} \mathcal{A} \xi_p - c_{p-1}\xi_1 + \frac{f^{(p)}(0)}{p!} \xi_1^p = 0. \label{epsp} \end{align} Let us denote for brevity \[ f_p := \frac{f^{(p)}(0)}{p!}. \] Pairing~\eqref{epsp} with $\xi_1$ (and assuming $\| \xi_1 \|_{L^2(0,L)} = 1$) gives \begin{equation} c_{p-1} = f_p \cdot \langle\xi_1^p, \xi_1\rangle_{L^2(0,L)}, \label{cp-1} \end{equation} which gives us the value of $c_{p-1}$. It only remains to solve the following problem numerically in order to obtain $\xi_p$: \begin{equation} \mathcal{A} \xi_p = c_{p-1}\xi_1 - f_p\xi_1^p \label{epsilonp} \end{equation} For the last equation to be solved, the operator $\mathcal{A}$ has to be invertible. It is also required that $\langle\xi_1, \xi_p\rangle_{L^2(0,L)}=0$. Therefore the solution is sought in the space orthogonal to $\ker(\mathcal{A})$. Since $\mathcal{A}$ is still a Fourier multiplier operator, one can take the Fourier transform of \eqref{epsilonp} to find \begin{gather*} \widehat{\mathcal{A}}(k) \widehat\xi_p(k) = c_{p-1}\widehat\xi_1(k) - f_p \widehat{\xi_1^p}(k), \\ \widehat\xi_p(k) = \widehat{\mathcal{A}}(k)^{-1} \left(c_{p-1}\widehat\xi_1(k) - f_p \widehat{\xi_1^p}(k)\right). \end{gather*} Taking the inverse Fourier transform of $\widehat\xi_p(k)$ gives $\xi_p$. Since only even solutions of the problem are considered the cosine part of the Fourier transforms will be required. It is sufficient to use $\xi_1$ and $c_0$ as the initial guesses for the Newton method. However, it should be noted that for different values of $p$ the pair of parameters $\xi_p$ and $c_{p-1}$ are computed in different ways. \begin{enumerate} \item[(a)] If $p = 2$, then $\xi_2$ is computed from~\eqref{epsilonp}, but $c_{p-1}$ here becomes zero. Therefore one has to consider the next level of the expansion $\varepsilon^p$. \item[(b)] For odd values of $p$ the parameter $c_{p-1}$ can be computed from~\eqref{cp-1} and $\xi_p$ from~\eqref{epsilonp}. \item[(c)] For even values of $p \geq 4$ the parameter $\xi_p$ can be computed, but $c_{p-1}$ may not be non-zero in general. In such cases a different strategy of fixing the initial guess should be used. \end{enumerate} \section{Appendix: Presentation of \textsf{SpecTraVVave} and its workflow} \label{sec:implementation} \subsection{Overview} There are several classes in the \textsf{SpecTraVVave} package. An over\-view of the program is shown in Figure \ref{schematic}. The workflow begins with defining a flux function $f$ and the Fourier multiplier function $\alpha$ to set up an equation. The traveling wave solution is characterized by the wavelength $L$ and a boundary condition $\Omega(c, a, \phi_N, B)$. These parameters are fixed for a given problem. The defined equation is then discretized. The \pyth{Discretization} object contains all required elements such as grid points, wave-numbers and the discrete linear operator. The initial guess and the equation's residual are passed from the \pyth{Discretization} to the \pyth{Solver} object. The \pyth{Navigation} object is responsible for finding good initial guesses for $c$ and $a$ that are passed to the \pyth{Solver} object. The \pyth{Solver} object applies Newton's method to find a solution to system of equations~\eqref{system-N3}. The new solution is sent back to the \pyth{Discretization} and \pyth{Navigation} objects, where variables get updated. All computed solutions are stored for further analysis. This finishes one iteration. For the next iteration the updated variables are used and a new solution is found. The process may be continued as long as the Jacobian of the problem is non-singular. \begin{figure}[ht] \begin{center} \includegraphics[width=0.8\textwidth]{fig12} % block-scheme \end{center} \caption{Overview of the \textsf{SpecTraVVave} package.} \label{schematic} \end{figure} \subsection{Class Description} We present an overview of the classes used in \textsf{SpecTraVVave} package. Note that, since the package is under continuous modification and development, we describe here only the basic classes and functions the package. We refer to the package repository \cite{github} for up-do-date tutorials and installation instructions. The \pyth{Equation} class is the general class for all model equations. Its only role is to store a parameter $L$, the wavelength: \begin{python} class Equation(object): def __init__(self, L): self.length = L \end{python} A subclass of the \pyth{Equation} class has to implement two functions, \pyth{compute_kernel} and \pyth{flux}. The KdV model equation \[ u_t + \frac{3}{2} u u_x + u_x + \frac{1}{6} u_{xxx} = 0 \] with $f(u) = \frac{3}{4} u^2$ and $\widehat{\mathcal{L} u}(k) = (1 - \frac{1}{6} k^2)\widehat{u}(k)$ is presented in the program as a subclass of the \pyth{Equation} class: \begin{python} class KDV(Equation): def compute_kernel(self, k): return 1-1/6*k*k def flux(self, u): return 3/4*u*u \end{python} One can then create an object of the class \pyth{KDV} with the command: \begin{python} kdv = KDV(L=np.pi) \end{python} % The solver will compute only a half of a solutions profile. The full wavelength of the solutions of the defined equation will be equal to $2\pi$. In order to find solutions with specific features, boundary conditions are introduced as separate classes. For instance, the boundary condition specifying a constant of integration is implemented as follows: \begin{python} class Const(object): """ The boundary condition under which the constant of integration (B) is always set to zero. """ def __init__(self, level=0): self.level = level def enforce(self, wave, variables, parameters): """ Enforces the Const boundary condition. """ return np.hstack([variables[0] - self.level]) def variables_num(self): """ The number of additional variables that are required to construct the Const boundary condition. """ return 1 \end{python} A \pyth{Const} boundary condition object is created as follows: \begin{python} boundary_condition = Const() \end{python} The next step is to create an object of \pyth{Discretization} class, which is initialized with a model equation such as \pyth{kdv_model} and the number of grid points. The main parts of the class are the following: \begin{python} class Discretization(object): def __init__(self, model_equation, grid_size): self.equation = model_equation self.size = grid_size def operator(self, u): u_ = scipy.fftpack.dct(u, norm='ortho') Lv = self.fourier_multiplier()*u_ result = scipy.fftpack.idct(Lv, norm='ortho') return result def residual(self, u, wavespeed, const_B): residual = - wavespeed*u + self.equation.flux(u) + self.operator(u) - const_B return residual \end{python} The call \pyth{Discretization.operator(u)} computes $\mathcal{L} u $ as the inverse transform of a transformed convolution \[ \mathcal{L} u = \mathcal{F}^{-1}[ \widehat{\mathcal{L} u}(k)] = \mathcal{F}^{-1}\left[ \alpha(k)\cdot\widehat{u}(k)\right]. \] The result of the call \pyth{Discretization.residual(u, wavespeed, const_B)} is then used in the \pyth{Solver} class. An object of the \pyth{Solver} class is initialized with an object of the \pyth{Discretization} class, and a boundary condition object. \begin{python} class Solver(object): def __init__(self, discrete_problem, boundary_condition): self.discretization = discrete_problem self.boundary = boundary_condition def solve(self, guess_wave, parameter_anchor, direction): """ Runs a Newton solver on a system of nonlinear equations once. Takes the residual(vector) to solve. parameter_anchor is the initial guess for (c,a) values and it is taken from the Navigation class. """ size = len(guess_wave) self.discretization.size = size def residual(vector): """ Contructs a system of nonlinear equations. First part, main_residual, is from given wave equation; second part, boundary_residual, comes from the chosen boundary conditions. """ . . . return np.hstack([main_residual, boundary_residual, amplitude_residual]) . . . return new_wave, new_boundary_variables, new_parameter \end{python} % Some omitted parts in the above script are substituted by \textsf{'. . .'} sign. Each iteration on a \pyth{Solver} object is run from a \pyth{Navigation} object, which takes the \pyth{Solver} object for initialization. \begin{python} class Navigation(object): """ Runs the Solver and stores the results. """ def __init__(self, solver_object, size=32): self.solve = solver_object.solve # function to run Newton method self.size = size # size for navigation . . . def run_solver(self, current_wave, pstar, direction): new_wave, variables, p3 = self.solve(current_wave, pstar, direction) return new_wave, variables, p3 \end{python} All the above classes can be modified and developed further, new classes may be defined as well. \subsection{Detailed Workflow} The workflow with the package consists of three basic steps: \begin{enumerate} \item Once the necessary classes have been imported in the current namespace, generate all necessary objects: \begin{python} equation = KDV(L=np.pi) boundary_condition = Const() discretization = Discretization(equation, grid_size=64) solver = Solver(discretization, boundary_condition) navigator = Navigation(solver) \end{python} \item Choose a number of iterations, i.e., the number of solutions to compute, and run the solver: \begin{python} n_iter = 50 navigator.run(n_iter) \end{python} \item All computed solutions are stored in \pyth{navigation_object} \begin{python} last_computed = -1 wave_profile = navigator[last_computed]['solution'] wave_speed = navigator[last_computed]['current'][0] wave_amplitude = navigator[last_computed]['current'][1] \end{python} \end{enumerate} For up-to-date instructions on how to run the code, we refer the reader to the code repository \url{https://github.com/olivierverdier/SpecTraVVave}. \subsection*{Acknowledgments} The authors would like to thank Mats Ehrnstr\"om and Erik Wahl\'en for fruitful discussions on the subject of the current paper. This research was supported by the Research Council of Norway on grant no. 213474/F20 and by the J C Kempe Memorial Fund on grant no. SMK-1238. \begin{thebibliography}{00} \bibitem{ABFS} Abdelouhab, L.; Bona, J. L.; Felland, M. Saut, J.-C.; {\em Nonlocal models for nonlinear dispersive waves}. Physica D \textbf{40},360-392 (1989). \bibitem{APR} Albert, J. P.; Bona, J. L.; Restrepo, J.M.; {\em Solitary-wave solutions of the Benjamin equation}. SIAM J. Appl. Math. \textbf{59},2139--2161 (1999). \bibitem{Duran2011JCAM} \'{A}lvarez, J.; Dur\'{a}n, A.; {\em A numerical scheme for periodic travelling-wave simulations in some nonlinear dispersive wave models}. J. Comput. Appl. Math. \textbf{235}, 1790--1797 (2011). \bibitem{Duran2014CNSNS} \'{A}lvarez, J.; Dur\'{a}n, A.; {\em An extended Petviashvili method for the numerical generation of traveling and localized waves}. Commun. Nonlinear Sci. Numer. Simul. \textbf{19}, 2272--2283 (2014). \bibitem{Duran2014JCAM} \'{A}lvarez, J.; Dur\'{a}n, A.; {\em, Petviashvili type methods for traveling wave computations: I. Analysis of convergence}. J. Comput. Appl. Math. \textbf{266}, 39--51 (2014). \bibitem{Duran2014JCAMCor} \'{A}lvarez, J.; Dur\'{a}n, A.; {\em Corrigendum to "Petviashvili type methods for traveling wave computations: I. Analysis of convergence}, [J. Comput. Appl. Math. 266 (2014) 39--51]. J. Comput. Appl. Math. \textbf{277}, 215--216 (2015). \bibitem{B1} Benjamin, T.B .; {\em Internal Waves of permanent form in fluids of great depth}. J. Fluid Mech. \textbf{29}, 559--592 (1967). \bibitem{B2} Benjamin, T. B.; {\em The stability of solitary waves}. Proc. Roy. Soc. London A 328, 153-183 (1972). \bibitem{B3} Benjamin, T. B.; {\emph{A new kind of solitary wave}}. J. Fluid Mech. \textbf{245}, 401--411 (1992). \bibitem{BDKM} Bona, J. L.; Dougalis, V.A.; Karakashian, O. A.; McKinney, W.R.; {\em Conservative, high-order numerical schemes for the generalized Korteweg--de Vries equation}, Philos. Trans. Royal Soc. London Ser. A, \textbf{351}, 107--164 (1995). \bibitem{BoKa} Bona, J. L.; Kalisch, H.; {\em Singularity formation in the generalized Benjamin--Ono equation}. Discrete Contin. Dyn. Syst. \textbf{11}, 779--785 (2004). \bibitem{BSS} Bona, J. L.; Souganidis, P. E.; Strauss, W. A.; {\em Stability and instability of solitary waves of Korteweg--de Vries type}. Proc. R. Soc. Lond. A \textbf{411}, 395-412 (1987). \bibitem{HKN} Borluk, H.; Kalisch, H.; Nicholls, D. P.; {\em A numerical study of the Whitham equation as a model for steady surface water waves}. J. Comput. Appl. Math. \textbf{296}, 293--302 (2016). \bibitem{TJB} Bridges, T. J.; {\em Superharmonic instability, homoclinic torus bifurcation and water-wave breaking}. J. Fluid Mech. \textbf{505}, 153-162 (2004). \bibitem{Chen} Chen, H.; {\em Existence of periodic travelling-wave solutions of nonlinear, dispersive wave equations}. Nonlinearity \textbf{17}, 2041--2056 (2004). \bibitem{CB} Chen, H.; Bona, J. L.; {\em Existence and asymptotic properties of solitary-wave solutions of Benjamin-type equations}. Adv. Diff. Eq. \textbf{3}, 51-84 (1998). \bibitem{ChenWalsh2016} Chen, R. M.; Walsh, S.; {\em Continuous dependence on the density for stratified steady water waves}. Arch. Ration. Mech. Anal. \textbf{219}, 741--792 (2016). \bibitem{CC2} Choi, W.; Camassa, R.; {\em Fully nonlinear internal waves in a two-fluid system}. J. Fluid Mech. \textbf{396}, 1-36 (1999). \bibitem{Tjon} Courtenay Lewis, J.; Tjon, J. A.; {\em Resonant production of solitons in the RLW equation}. Phys. Lett. A \textbf{73}, 275--279 (1979). \bibitem{dFSS} de Frutos, J.; Sanz-Serna, J. M.; {\em An easily implementable fourth-order method for the time integration of wave problems}. J. Comp. Phys. \textbf{103}, 160--168 (1992). \bibitem{Doed97} Doedel, E. J.; Champneys, A. R.; Fairgrieve, T. F.; Kuznetsov, Y.A.; Sandstede, B. and Wang, X.-J.; {\em AUTO97 : Continuation and bifurcation software for ordinary differential equations}, available at \url{http://indy.cs.concordia.ca/auto/}, Department of Computer Science and Software Engineering, Concordia University, Montreal, Canada, 1997. \bibitem{DDM} Dougalis, V. A.; Duran, A.; Mitsotakis, D.; {\em Numerical solution of the Benjamin equation}. Wave Motion \textbf{52}, 194--215 (2015). \bibitem{EK1} Ehrnstr{\"o}m, M.; Kalisch, H.; {\em Traveling waves for the Whitham equation}. Differential Integral Equations \textbf{22}, 1193--1210 (2009). \bibitem{EK2} Ehrnstr{\"o}m, M., Kalisch, H.; {\em Global bifurcation for the Whitham equation}. Math. Mod. Nat. Phenomena \textbf{8}, 13--30 (2013). \bibitem{EW2016} Ehrnstr{\"o}m, M.; Wahl\'en, E.; {\em On Whitham's conjecture of a highest cusped wave for a nonlocal dispersive equation}. arXiv:1602.05384 (2016). \bibitem{FW} Fornberg, B.; Whitham, G. B.; {\em A numerical and theoretical study of certain nonlinear wave phenomena}. Phil. Trans. Roy. Soc. A \textbf{289}, 373--404 (1978). \bibitem{PyBook} F\"uhrer, C.; Solem, J. E.; Verdier, O.; Computing with Python: An introduction to Python for science and engineering, Pearson (2014). \bibitem{GK} Gruji\'{c}, Z.; Kalisch, H.; {\em Gevrey regularity for a class of water-wave models}. Nonlinear Analysis \textbf{71}, 1160--1170 (2009) \bibitem{Hur2015} Hur, V. M.; {\em Breaking in the Whitham equation for shallow water waves}. arXiv:1506.04075 (2015). \bibitem{HurJohnson} Hur, V. M.; Johnson, M.; {\em Modulational instability in the Whitham equation of water waves}. Studies in Applied Mathematics \textbf{134}, 120--143 (2015). \bibitem{HurTao} Hur, V. M.; Tao, L.; {\em Wave breaking for the Whitham equation with fractional dispersion}. arXiv:1410.1570, 2014. \bibitem{HK} Kalisch, H.; {\em Error analysis of a spectral projection of the regularized Benjamin-Ono equation}. BIT Numerical Mathematics \textbf{45}, 69--89 (2005). \bibitem{KalischBenjamin} Kalisch, H.; {\em Derivation and comparison of model equations for interfacial capillary-gravity waves in deep water}. Math. Comput. Simulation \textbf{74}, 168--178 (2007). \bibitem{HKJB} Kalisch, H.; Bona, J. L.; {\em Models for internal waves in deep water}. Discrete and Continuous Dynamical Systems \textbf{6}, 1--20 (2000). \bibitem{Kato} Kato, T.; {\em On the Cauchy problem for the (generalized) Korteweg--de Vries equation}. Studies in applied mathematics, 93--128, Adv. Math. Suppl. Stud., \textbf{8} (Academic Press, New York, 1983). \bibitem{Keller} Keller H. B.; {\em Numerical solution of bifurcation and nonlinear eigenvalue problems}. Applications of Bifurcation Theory, Rabinowitz, P. H., Academic Press, (1977), pp. 359--384. \bibitem{LS} Lannes, D.; Saut, J.-C.; {\em Remarks on the full dispersion Kadomtsev--Petviashvli equation}. Kinet. Relat. Models \textbf{6}, 989--1009 (2013). \bibitem{MM} Martel, Y.; Merle, F.; {\em Blow up in finite time and dynamics of blow up solutions for the L2-critical generalized KdV equation}. J. Amer. Math. Soc. \textbf{15}, 617--664 (2002). \bibitem{MP} Martel, Y.; Pilod, D.; {\em Construction of a minimal mass blow up solution of the modified Benjamin-Ono equation}. arXiv:1605.01837 (2016). \bibitem{MKD2015} Moldabayev, D.; Kalisch, H.; Dutykh, D.; {\em The Whitham Equation as a model for surface water waves}. Phys. D \textbf{309}, 99--107 (2015) \bibitem{github} Moldabayev, D.; Verdier, O.; Kalisch, H.; {\em SpecTraVVave}, available at \newline \url{https://github.com/olivierverdier/SpecTraVVave}. \bibitem{NgK} Nguyen, N. T.; Kalisch, H.; {\em Orbital stability of negative solitary waves}. Mathematics and Computers in Simulation \textbf{80}, 139-150 (2009). \bibitem{O} Ono, H.; {\emph{Algebraic Solitary Waves in Stratified Fluids}}. J. Phys. Soc. of Japan \textbf{39}, 1082-1091 (1975). \bibitem{Pelinovsky} Pelinovsky, D. E.; Stepanyants, Y. A.; {\em Convergence of Petviashvili's iteration method for numerical approximation of stationary solutions of nonlinear wave equations}. SIAM J. Numer. Anal. \textbf{42}, 1110--1127 (2004). \bibitem{IPython} P\'erez, F.; Granger, B. E.; \emph{IPython: A System for Interactive Scientific Computing}, Computing in Science and Engineering, vol. 9, no. 3, pp. 21-29, May/June 2007. URL: \url{http://ipython.org} \bibitem{Petviashvili} Petviashvili, V.I.; {\em Equation of an extraordinary soliton}, Sov. J. Plasma Phys. \textbf{2}, 257--258 (1976). \bibitem{RK} Remonato, F.; Kalisch, H.; {\em Numerical bifurcation for the capillary Whitham equation}, Physica D: Nonlinear Phenomena, 343, 51--62 (2017). \bibitem{Sanford} Sanford, N.; Kodama, K.; Carter, J. D.; Kalisch, H.; {\em Stability of traveling wave solutions to the Whitham equation}. Phys. Lett. A \textbf{378}, 2100--2107 (2014). \bibitem{Sherratt} Sherratt, J. A.; {\em Numerical continuation methods for studying periodic travelling wave (wavetrain) solutions of partial differential equations}. Appl. Math. Comp. \textbf{218}, 4684--4694 (2012). \bibitem{Walsh2014a} Walsh, S.; {\em Steady stratified periodic gravity waves with surface tension I: Local bifurcation}. Discrete Contin. Dyn. Syst. \textbf{34}, 3241--3285 (2014). \bibitem{Walsh2014b} Walsh, S.; {\em Steady stratified periodic gravity waves with surface tension II: global bifurcation}. Discrete Contin. Dyn. Syst. \textbf{34}, 3287--3315 (2014). \bibitem{Wh1} Whitham, G. B.; {\em Variational methods and applications to water waves}. Proc. Roy. Soc. London A, \textbf{299}, 6--25 (1967). \bibitem{Wh2} Whitham, G. B.; Linear and Nonlinear Waves. Wiley, New York, (1974). \end{thebibliography} \end{document}