Numerical bifurcation and stability for
the capillary-gravity Whitham equation

Whitham in his 1967 paper [30] (see also [31]) put forward

to address, qualitatively, breaking and peaking of waves in shallow water. Here and throughout, $u(x,t)$ is related to the displacement of the fluid surface from the rest state at position $x$ and time $t$, and $c_{\rm{\tiny ww}}(|\partial_{x}|)$ is a Fourier multiplier operator, defined as

where $g$ is the gravitational constant and $h$ the undisturbed fluid depth. Notice that $c_{\text{ww}}(k)$ is the phase speed of $2\pi/k$ periodic waves in the linear theory of water waves [31].

For relatively shallow water or, equivalently, relatively long waves for which $kh\ll 1$, one can expand the right hand side of (2) to obtain

whereby arriving at the celebrated Korteweg–de Vries (KdV) equation:

Therefore (1) and (2) can be regarded as augmenting (3) to include the full dispersion in the linear theory of water waves and, thus, improving over (3) for short and intermediately long waves. Indeed, numerical studies (see [25], for instance) reveal that the Whitham equation performs on par with or better than the KdV equation and other shallow water models in a wide range of amplitude and wavelength parameters.

The KdV equation admits solitary and cnoidal waves but no traveling waves can ‘peak’. By contrast, the so-called extreme Stokes wave possesses a $120^{\circ}$ peaking at the crest. Also no solutions of (3) can ‘break’. That means, the solution remains bounded but its slope becomes unbounded. See [31, Section 13.14] for more discussion. This is perhaps not surprising because the dispersion^∗*∗*$*$The phase speed is $\sqrt{gh}(1-\frac{1}{6}k^{2}h^{2})$, poorly approximating $c_{\text{ww}}(k)$ when $kh\gg 1$. of the KdV equation is inadequate for explaining high frequency phenomena of water waves. Whitham [30, 31] conjectured breaking and peaking for (1) and (2). Recently, one of the authors [13] proved breaking, and Ehrnström and Wahlén [9] proved peaking. Johnson and one of the authors [15] proved that a small amplitude and periodic traveling wave of (1) and (2) is modulationally unstable, provided that $kh>1.145\dots$, comparing with the well-known Benjamin–Feir instability of a Stokes wave. By contrast, all cnoidal waves of (3) are modulationally stable.

When the effects of surface tension are included, Johnson and one of the authors [16] proposed to replace (2) by

where $T$ is the surface tension coefficient. Notice that $c_{\text{ww}}(k;T)$ is the phase speed in the linear theory of capillary-gravity waves [31]. When $T=0$, (4) becomes (2). Since

one arrives at the KdV equation for capillary-gravity waves:

in the long wavelength limit unless $T/gh^{2}=1/3$. Notice that (3) and (5) behave alike, qualitatively, possibly after a sign change. By contrast, whenever $T>0$, $c_{\text{ww}}(k;T)\to\sqrt{T|k|}$ as $kh\to\infty$, so that (1) and (4) become

to leading order whereas when $T=0$, $u_{t}+\sqrt{g}|\partial_{x}|^{-1/2}u_{x}+\frac{3}{2}\sqrt{\frac{g}{h}}uu_{x}=0$.

In recent years, the capillary-gravity Whitham equation has been a subject of active research [16, 27, 5, 3, 8] (see also [19, 18, 14, 26, 21]). Particularly, Remonato and Kalisch [27] combined a spectral collocation method and a numerical continuation approach to discover a rich variety of local bifurcation of periodic traveling waves of (1) and (4), notably, crossing and connecting solution branches. Here we adopt a robust numerical continuation scheme to corroborate the results and produce convincing results for global bifurcation. Our findings support local bifurcation theorems (see [8], for instance) and help classify all periodic traveling waves.

We employ an efficient numerical method for solving stiff nonlinear PDEs implemented with fourth-order time differencing, to experiment with (nonlinear) orbital stability and instability for a plethora of periodic traveling waves of (1) and (4). (Spectral) modulational stability and instability were investigated numerically [3] and also analytically for small amplitude [16]. To the best of the authors’ knowledge, however, nonlinear stability and instability have not been addressed. Our novel findings include, among many others, orbital stability for the $k=1$ branch whenever $T>0$ versus instability for $k\geqslant 2,\in\mathbb{N}$ branches for great wave height when $0<T/gh^{2}\leqslant 1/3$.

The methodology here is potentially useful for tackling the capillary-gravity wave problem and other nonlinear dispersive equations.

We rewrite (1) and (4) succinctly as

Seeking a periodic traveling wave of (7), let

where $c\in\mathbb{R}$ ($\neq 0$) is the wave speed, $k>0$ the wave number, and $\phi$ satisfies, by quadrature and Galilean invariance,

We assume that $\phi$ is $2\pi$ periodic and even in the $z$ variable, so that $2\pi/k$ periodic in the $x$ variable. Notice

For any $T\geqslant 0$, $k>0$ and $c\in\mathbb{R}$, clearly, $\phi(z)=0$ solves (8) and (9). Suppose that $T$ and $k$ are fixed. A necessary condition for nontrivial solutions to bifurcate from such trivial solution at some $c$ is that

if and only if $c=c_{\text{ww}}(kn;T)$ for some $n\in\mathbb{N}$, by symmetry, whence

When $T=0$, $c_{\text{ww}}(\cdot\,;0)$ monotonically decreases to zero over $(0,\infty)$, whence for any $k>0$, $c_{\text{ww}}(kn_{1};T)>c_{\text{ww}}(kn_{2};T)$ whenever $n_{1},n_{2}\in\mathbb{N}$ and $n_{1}>n_{2}$. Therefore, for any $k>0$ for any $n\in\mathbb{N}$, $\ker(c_{\text{ww}}(k|\partial_{z}|;0)-c_{\text{ww}}(kn;0))=\text{span}\{\cos(nz)\}$ in spaces of even functions.

When $T\geqslant 1/3$, $c_{\text{ww}}(\cdot\,;T)$ monotonically increases over $(0,\infty)$ and unbounded from above, whereby for any $k>0$ for any $n\in\mathbb{N}$, $\ker(c_{\text{ww}}(k|\partial_{z}|;T)-c_{\text{ww}}(kn;T))=\text{span}\{\cos(nz)\}$.

When $T<1/3$, to the contrary, $c_{\text{ww}}(kn_{1};T)=c_{\text{ww}}(kn_{2};T)$ for some $k>0$ for some $n_{1},n_{2}\in\mathbb{N}$ and $n_{1}\neq n_{2}$, provided that $T=T(kn_{1},kn_{2})$, where

and $\ker(c_{\text{ww}}(k|\partial_{z}|;T)-c)=\text{span}\{\cos(n_{1}z),\cos(n_{2}z)\}$, where $c=c_{\text{ww}}(kn_{1};T)=c_{\text{ww}}(kn_{2};T)$; otherwise, $\ker(c_{\text{ww}}(k|\partial_{z}|;T)-c_{\text{ww}}(kn;T))=\text{span}\{\cos(nz)\}$.

Suppose that $T\geqslant 0$, $k>0$, and $T\neq T(kn,kn^{\prime})$ for any $n,n^{\prime}\in\mathbb{N}$ and $n\neq n^{\prime}$, particularly, either $T=0$ or $T\geqslant 1/3$. We assume without loss of generality $n=1$. There exists a one-parameter curve of nontrivial, $2\pi$ periodic and even solutions of (8) and (9), denoted by

in some function space (see [9, 8], for instance, for details), and

as $|s|\to 0$. Moreover, subject to a ‘bifurcation condition’ (see [9, 8], for instance, for details), (11) extends to all $s\in[0,\infty)$. See [9, 8] and references therein for a rigorous proof.

When $T=0$ and, without loss of generality, $k=1$, $(\phi(z;s_{j}),c(s_{j}))\to(\phi(z),c)$ as $j\to\infty$ for some $s_{j}\to\infty$ as $j\to\infty$ for some $\phi\in C^{1/2}([-\pi,\pi])\bigcap C^{\infty}([-\pi,0)\bigcup(0,\pi])$ and $0<c<\infty$ such that

See [9] and references therein for a rigorous proof. Indeed, such a limiting solution enjoys $\phi(z)\sim\frac{c}{2}-\sqrt{\frac{\pi}{8}}|z|^{1/2}$ as $|z|\to 0$.

When $T\geqslant 4/\pi^{2}$, so that the Fourier transform of $c_{\text{ww}}(k|\partial_{z}|;T)^{-1}$ is ‘completely monotone’ [8], and $k=1$, on the other hand,

1. $\|(\phi(\cdot;s),c(s))\|\to\infty$ as $s\to\infty$ in some^††††$†$a Hölder-Zygmund space, for instance [8] function space; or

2. $(\phi(z;s),c(s))$ is periodic in the $s$ variable.

See [8], for instance, for a rigorous proof. Our numerical findings suggest $\min_{z\in[-\pi,\pi]}\phi(z;s)\to-\infty$ as $s\to\infty$. But such a limiting scenario would be physically unrealistic, for the capillary-gravity Whitham equation models water waves in the finite depth. We say that a $2\pi$ periodic solution of (8) and (9) is admissible if

so that the fluid surface would not intersect the impermeable bed, and we stop numerical continuation once we reach a limiting admissible solution, for which

Section 4 provides examples.

Suppose, to the contrary, that $T=T(kn_{1},kn_{2})$ (see (10)) for some $k>0$ for some $n_{1},n_{2}\in\mathbb{N}$ and $n_{1}<n_{2}$, so that $c_{\text{ww}}(kn_{1};T)=c_{\text{ww}}(kn_{2};T)=:c$. Suppose that $n_{1}$ does not divide $n_{2}$. There exists a two-parameter sheet of nontrivial, periodic and even solutions of (8) and (9), and

for $|s_{1}|,|s_{2}|\ll 1$. We emphasize that $k$ is a bifurcation parameter. Otherwise, $n_{1}$ divides $n_{2}$, and (15) holds for $|s_{1}|,|s_{2}|\ll 1$ and $s_{2}>s_{0}>0$ for some $s_{0}$. In other words, there cannot exist $2\pi/n_{1}$ periodic and ‘unimodal’ waves, whose profile monotonically decreases from its single crest to the trough over the period. See [8], for instance, for a rigorous proof. The global continuation of (15) and limiting configurations have not been well understood analytically, though, and here we address numerically.

Turning the attention to the (nonlinear) orbital stability and instability of a periodic traveling wave of (7), notice that (7) possesses three conservation laws:

and so does the KdV equation with fractional dispersion:

for which $E(u)=\int(\frac{1}{2}u|\partial_{x}|^{\alpha}u+\frac{1}{3}u^{3})~{}dx$ rather than the first equation of (16). We pause to remark that (7) becomes (17), $\alpha=1/2$, for high frequency (see (6)), and (3) and (5) compare with $\alpha=2$. A solitary wave of (17) is a constrained energy minimizer and orbitally stable, provided that $\alpha>1/2$, if and only if $P_{c}>0$. See [24] and references therein for a rigorous proof. See also [17] for an analogous result for periodic traveling waves. It seems not unreasonable to expect that the orbital stability and instability of a periodic traveling wave of (7) change likewise at a critical point of $P$ as a function of $c$ although, to the best of the authors’ knowledge, there is no rigorous analysis of constrained energy minimization. Indeed, numerical evidence [20] supports the conjecture when $T=0$ and $k=1$. Here we take matters further to $T\geqslant 0$ and $k\geqslant 1,\in\mathbb{N}$. Also we numerically elucidate the instability scenario when $T=0$ and $k=1$.

We begin by numerically approximating $2\pi$ periodic and even solutions of (8) and (9) by means of a spectral collocation method [2, 10]. See [20], among others, for nonlinear dispersive equations of the form (7) and (17). We define the collocation projection as a discrete cosine transform as

where the discrete cosine coefficients are

and

the collocation points are

Therefore

We can compute (18) and (19) efficiently using a fast Fourier transform (FFT) [10]. For $T\geqslant 0$ and $k>0$, likewise,

and we can compute $(c_{\text{ww}}(k|\partial_{z}|;T)\phi)_{N}(z_{m})$, $m=1,2,\dots,N$, via an FFT.

Suppose that $T\geqslant 0$ and $k>0$ are fixed and we take $N=1024$. We numerically solve

by means of Newton’s method. We parametrically continue the numerical solution over $c$ by means of a pseudo-arclength continuation method [7] (see [20], among others, for nonlinear dispersive equations of the form (7) and (17)), for which the (pseudo-)arclength of a solution branch is the continuation parameter, whereas a parameter continuation method would use $c$ as the continuation parameter and vary it sequentially. The pseudo-arclength continuation method can successfully bypass a turning point of $c$, at which a parameter continuation method fails because the Jacobian of (20) becomes singular, whence Newton’s method diverges. See [27] for more discussion.

We say that Newton’s method converges if the residual

(see (20)), and this is achieved provided that an initial guess is sufficiently close to a true solution of (8) and (9). To this end, we take a small amplitude cosine function and $c\approx c_{\text{ww}}(k;T)$ as an initial guess for the local bifurcation at $\phi=0$ and $c=c_{\text{ww}}(k;T)$, or (12) so long as it makes sense. We have also solved (20) using a (Jacobian-free) Newton–Krylov method [23] with absolute and relative tolerance of $10^{-10}$. The results correctly match those obtained from Newton’s method, corroborating our numerical scheme.

We have run a pseudo-arclength continuation code (see [29], for instance, for the applicability of the code used herein) with a fixed arclength stepsize, to numerically locate and trace solution branches. We have additionally run AUTO [6], a robust numerical continuation and bifurcation software, to corroborate the results. All the results in Section 4 have been obtained by employing AUTO with strict tolerances (EPSL = 1e-12 and EPSU=1e-12 in the AUTO’s constants file). AUTO has the advantage, among many others, of making variable arclength stepsize adaptation (using the option IADS=1) and of detecting branch points (using the option ISP=2 for all special points). For the former, provided with minimum and maximum allowed absolute values of the pseudo-arclength stepsize, AUTO adjusts the arclength stepsize to be used during the continuation process.

Let $\phi$ and $c$ numerically approximate a $2\pi/k$ periodic traveling wave of (7), and we turn to numerically experimenting with its (nonlinear) orbital stability and instability. After a $2048$-point (full) Fourier spectral discretization in $x\in[-\pi,\pi]$, (7) leads to

and we numerically solve (21), subject to

by means of an integrating factor (IF) method (see [22], for instance, and references therein), where $U(-1,1)$ is a uniform random distribution. In other words, at the initial time, we perturb $\phi$ by uniformly distributed random noise of small amplitude, depending on the amplitude^‡‡‡‡$‡$We treat $\|\phi\|_{L^{\infty}}$ as the amplitude of $\phi$ although $\phi$ is not of mean zero, because the difference is insignificant. of $\phi$. We pause to remark that the well-posedness for the Cauchy problem of (7) can be rigorously established at least for short time in the usual manner via a compactness argument.

Following the IF method, let $v=e^{-\mathbf{L}t}u$, so that (21) becomes

This ameliorates the stiffness of (21). Notice that (23) is of diagonal form, rather than matrix, in the Fourier space. We advance (23) forward in time by means of a fourth-order four-stage Runge–Kutta (RK4) method [11] with a fixed time stepsize $\Delta t$:

where $v_{n}=v(t_{n})$ and $t_{n}=n\Delta t$, $n=0,1,2,\dots$,

and $f(v,t)=e^{-\mathbf{L}t}\mathbf{N}(e^{\mathbf{L}t}v)$ (see (23)). We take $\Delta t=10^{-4}\times 2\pi/c$, where $c$ is the wave speed of $\phi$. While such value of $\Delta t$ ensures numerical stability, some computations, depending on $c$, require smaller^§§§§$§$In most cases, $10^{-4}\times 2\pi/c=O(10^{-4})$ but in Figure 7(b), for instance, $c=0.36$, whence $10^{-4}\times 2\pi/c\approx 0.00174533$, whereas $10^{-5}\times 2\pi/c=O(10^{-4})$. values, for instance, $\Delta t=10^{-4}\times\pi/c$ or $10^{-5}\times 2\pi/c$. At each time step, we remove aliasing errors by applying the so-called $3/2$-rule so that the Fourier coefficients well decay for high frequencies (see [2], for instance).

We assess the fidelity of our numerical scheme by monitoring (16) for numerical solutions. For unperturbed solutions, $E(t)$, $P(t),M(t)$ have all been conserved to machine precision, whereas for (randomly) perturbed ones, $E(t)$ and $P(t)$ to the order of $10^{-7}$ and $M(t)$ to machine precision. We have also corroborated our numerical results using two other methods: a higher-order Runge–Kutta method, such as the Runge–Kutta–Fehlberg (RKF) method [11] with time stepsize adaptation and a two-stage (and, thus, fourth-order) Gauss–Legendre implicit^¶¶¶¶$¶$An implicit method, by construction, enjoys wider regions of absolute stability. This allows us to assess all our results, while avoiding numerical instability that might be observed in explicit methods such as RK4 and RKF methods. However, the IRK4 method requires a fixed point iteration (in the form of Newton’s method) at each time step. Runge–Kutta (IRK4) method [12]. For stable solutions, the results obtained from the RK4 method correctly match those from the RKF and IRK4 methods. However, and for unstable solutions, the instability manifests at slightly different times for the same initial conditions. This is somewhat expected because the local truncation error (LTE) of each method and the convergence error of the IRK4 method act as perturbations. Recall that the RK4 and IRK4 methods have an LTE of the order of $(\Delta t)^{4}$ whereas the RKF method has $(\Delta t)^{6}$.

We begin by taking $T=0$ and $k=1$. Figure 1 shows the wave height

and the momentum

from our numerical continuation of $2\pi$ periodic and even solutions of (8) and (9). The result agrees, qualitatively, with [20, Figure 4] and others. We find that $H$ monotonically increases to $\approx 0.58156621$, highlighted with the red square, for which (13) holds, whereas $P$ increases and then decreases, making one turning point. Also we find that the crest becomes sharper and the trough flatter as $H$ increases. The limiting solution must possess a cusp at the crest [9].

The left column of Figure 2 shows almost limiting waves. The inset is a close-up near the crest, emphasizing smoothness. The right column shows the profiles of (22), namely periodic traveling waves of (7) perturbed by uniformly distributed random noise of small amplitudes at $t=0$ (dash-dotted), and of the solutions of (7) at later instants of time (solid). Wave (a), prior to the turning point of $P$, remains unchanged for $10^{3}$ time periods, after translation of the $x$ axis, implying orbital stability, whereas the inset reveals that wave (b), past the turning point, suffers from crest instability. Indeed, our numerical experiments point to transition from stability to instability at the turning point of $P$. But there is numerical evidence [28] that waves (a) and (b) are (spectrally) modulationally unstable.

When $T=4/\pi^{2}$ and $k=1$, Figure 3 shows $H$ and $P$ versus $c$. Our numerical findings suggest that $H,P\to\infty$ monotonically. Also $\min_{z\in[-\pi,\pi]}\phi(z)\to-\infty$. The solution branch discontinues once (14) holds, highlighted with the red square, though, because the solution would be physically unrealistic, for the capillary-gravity Whitham equation models water waves in the finite depth. Our numerical continuation works well past the limiting admissible solution, nevertheless. We find that the crest becomes wider and flatter, the trough narrower and more peaked, as $H$ increases. But all solutions must be smooth [8].

The left panel of Figure 4 shows an example wave. The right panel shows the profiles perturbed by small random noise at $t=0$ and of the solution of (7) after $10^{3}$ time periods, after translation of the $x$ axis. Our numerical experiments suggest orbital stability for all wave height. But there is numerical evidence of modulational instability when $H$ is large [3].

Numerical evidence [27, 3] points to qualitatively the same results whenever $T\geqslant 1/3$.

We turn the attention to $T=T(1,2)\approx 0.2396825$ (see (10)). There exists a two-parameter sheet of nontrivial, periodic and even solutions of (8) and (9) in the vicinity of $\phi=0$ and $c=c_{\text{ww}}(1,T)=c_{\text{ww}}(2,T)\approx 0.97166609$ [8]. See also Section 2. Figure 5 shows $H$ versus $c$ for the $k=1$ and $k=2$ branches, all the way up to the limiting admissible solutions. There are no^∥∥∥∥$∥$We observe a turning point of $P$ for greater wave height, but the solution is inadmissible. turning points of $P$. The left column of Figure 6 shows waves in the $k=1$ and $k=2$ branches for small and large $H$. The small height result agrees with [27, Figure 6].

Observe ‘bimodal’ waves in the $k=1$ branch. Indeed, there cannot exist $2\pi$ periodic and unimodal waves, whose profile monotonically decreases from a single crest to the trough over the period [8]. See also Section 2. For small wave height, the fundamental mode seems dominant, so that there is one crest over the period $=2\pi$, but the fundamental and second modes are resonant, whereby a much smaller wave breaks up the trough into two. See the left panel of Figure 6(a). As $H$ increases, the effects of the second mode seem more pronounced, so that the wave separating the troughs becomes higher. See the left of Figure 6(b). Observe, to the contrary, $\pi$ periodic and unimodal waves in the $k=2$ branch. See the left of Figure 6(c,d). We find that the crests become wider and flatter, the troughs narrower and more peaked, as $H$ increases in the $k=1$ and $k=2$ branches. See the left of Figure 6(b,d).

Our numerical experiments suggest orbital stability for the $k=1$ branch (see the right panels of Figure 6(a,b)) versus instability for $k=2$ (see the right of Figure 6(c,d)).

We take matters further to $T=T(1,3)$. See Figures 7, 8 and 9. The results for the $k=1$ and $k=3$ branches are similar to those when $T=T(1,2)$ and $k=1,2$. We pause to remark that in the $k=1$ branch, for small $H$, a smaller wave breaks up the trough into two, and a much smaller wave breaks up the crest into two. See the left panel of Figure 8(a). As $H$ increases, the wave separating the crests becomes lower, transforming into one wide and flat crest, whereas the wave separating the troughs become higher (see the left of Figure 8(b)), whereby resembling those when $T=T(1,2)$ and $k=1$. Observe $\pi$ periodic and unimodal waves in the $k=2$ branch, orbitally stable for small $H$ versus unstable for large $H$. See Figure 9(e,f).

When $T=T(2,3)$, to the contrary, we find $\pi$ and $2\pi/3$ periodic unimodal waves in the $k=2$ and $k=3$ branches, respectively, corroborating a local bifurcation theorem (see [8], for instance). See also Section 2. We find that the crests become wider and flatter, the troughs narrower and more peaked, as $H$ increases. See the left column of Figure 11. Our numerical experiments (see the right of Figure 11) suggest orbital instability for the $k=2$ branch for large $H$ and for $k=3$ for all $H$.

When $T=T(2,5)$, on the other hand, the left column of Figure 13 shows bimodal waves in the $k=2$ branch. The local bifurcation theorem [8] dictates $\pi$ periodic and unimodal waves, but we numerically find that they are not in the $k=2$ branch. For small $H$, for wave (a), for instance, a smaller wave breaks up the trough into two over the half period. As $H$ increases, for wave (b), for instance, the troughs become narrower and more peaked. Our numerical experiments (see the right of Figures 13 and 14) suggest orbital stability for the $k=2$ branch for small $H$ versus instability for $k=2$ for large $H$ and for $k=5$ for all $H$.

To proceed, when $T=T(1,4)$, the $k=1$ branch lies above and to the right of the $k=4$ branch at least for small $H$ (not shown), but as $T$ increases, $c_{\text{ww}}(4;T)$ increases more rapidly than $c_{\text{ww}}(1;T)$, so that when $T=T(1,4)+0.001$, for instance, the $k=1$ branch crosses the $k=4$ branch. Figure 15 shows $H$ and $P$ versus $c$ in the $k=1$ and $k=4$ branches for all admissible solutions. The small height result agrees with [27, Figure 10(a)]. We find that in the $k=1$ branch, $H$ and $P$ turn to connect the $k=4$ branch, whereas in the $k=4$ branch, $H$ and $P$ monotonically increase to the limiting admissible solution, highlighted by the red square in the insets.

The left panels of Figures 16, 17 and 18 show several profiles along the $k=1$ and $k=4$ branches. In the $k=1$ branch, for small $H$, wave (a), for instance, is $2\pi$ periodic and unimodal. After the $k=1$ branch crosses the $k=4$ branch, on the other hand, wave (b), for instance, becomes bimodal, resembling those when $T=T(1,4)$ and $k=1$. Continuing along the $k=1$ branch, for waves (c) and (d), for instance, high frequency ripples of $k=4$ ride a carrier wave of $k=1$. When the $k=1$ and $k=4$ branches almost connect, wave (e) in the $k=1$ branch and wave (f) for $k=4$ are almost the same. In the $k=4$ branch, wave (g), for instance, is $\pi/2$ periodic and unimodal.

Our numerical experiments (see the right panels of Figures 16, 17 and 18) suggest orbital stability for the $k=1$ branch for all $H$ and for $k=4$ for small $H$, versus instability for $k=4$ for large $H$. Particularly, stability and instability do not change at the turning point of $P$.

Last but not least, when $T=T(1,5)+0.0001$, Figure 19 shows $H$ and $P$ versus $c$ for the $k=1$ and $k=5$ branches. The small height result agrees with [27, Figure 10(b)]. We find that the $k=1$ branch crosses and connects the $k=5$ branch, like when $T=T(1,4)+0.0001$ and $k=1,4$, but it continues after connecting all the wave up to the limiting admissible solution. See the insets. The left panels of Figures 20, 21 and 22 show several profiles along the $k=1$ and $k=5$ branches. The results for the $k=1$ branch up to connecting and for $k=5$ are similar to those when $T=T(1,4)+0.0001$ and $k=1,4$. In the $k=1$ branch, after it connects the $k=4$ branch at the point (c), the results are similar to those when $T=T(1,5)$ and $k=1$. See waves (d), (e) and (f). Our numerical experiments (see the right panels of Figures 20, 21 and 22) suggest orbital stability for the $k=1$ branch for all $H$ and for $k=5$ for small $H$, versus instability for $k=5$ for large $H$.

We emphasize orbital stability for the $k=1$ branch for all $T>0$ for all wave height throughout our numerical experiments.

Here we employ efficient and highly accurate numerical methods for computing periodic traveling waves of (7) and experimenting with their (nonlinear) orbital stability and instability. Our findings suggest, among many others, stability whenever $T>0$ for the $k=1$ branch versus instability when $0<T<1/3$ for $k\geqslant 2,\in\mathbb{N}$ branches at least for large wave height. Currently under investigation is to take matters further to all $T\geqslant 0$ to all $k$, to classify the orbital stability and instability of all periodic traveling waves. It will be interesting to devise other numerical continuation methods, for instance, deflated continuation techniques [1, 4], for detecting disconnected solution branches and others. Also interesting will be to numerically investigate (spectral) modulational stability and instability of orbitally stable waves. Our methodology will be useful for exploring the nonlinear dynamics of modulationally unstable waves. Also it can help tackle the capillary-gravity wave problem and other nonlinear dispersive equations. Of course, it is of great importance to rigorously prove the numerical results.