Numerical bifurcation and stability for the capillary-gravity Whitham equation
aa r X i v : . [ phy s i c s . f l u - dyn ] F e b NUMERICAL BIFURCATION AND STABILITY FORTHE CAPILLARY-GRAVITY WHITHAM EQUATION
EFSTATHIOS G. CHARALAMPIDIS AND VERA MIKYOUNG HUR
Abstract.
We adopt a robust numerical continuation scheme to examine the global bifurcation ofperiodic traveling waves of the capillary-gravity Whitham equation, which combines the dispersionin the linear theory of capillary-gravity waves and a shallow water nonlinearity. We employ a highlyaccurate numerical method for space discretization and time stepping, to address orbital stabilityand instability for a rich variety of the solutions. Our findings can help classify capillary-gravitywaves and understand their long-term dynamics. Introduction
Whitham in his 1967 paper [30] (see also [31]) put forward(1) u t + c ww ( | ∂ x | ) u x + 32 r gh uu x = 0to 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 ww ( | ∂ x | ) is a Fourier multiplier operator, defined as(2) V c ww ( | ∂ x | ) f ( k ) = r g tanh( kh ) k b f ( k ) , where g is the gravitational constant and h the undisturbed fluid depth. Notice that c ww ( k ) is thephase speed of 2 π/k periodic waves in the linear theory of water waves [31].For relatively shallow water or, equivalently, relatively long waves for which kh ≪
1, one canexpand the right hand side of (2) to obtain c ww ( k ) = p gh (cid:16) − k h (cid:17) + O ( k h ) , whereby arriving at the celebrated Korteweg–de Vries (KdV) equation:(3) u t + p gh (cid:16) u x + 16 h u xxx (cid:17) + 32 r gh uu x = 0 . Therefore (1) and (2) can be regarded as augmenting (3) to include the full dispersion in thelinear 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 parwith or better than the KdV equation and other shallow water models in a wide range of amplitudeand wavelength parameters.
Mathematics Department, California Polytechnic State University San Luis Obispo, CA 93407-0403,USADepartment of Mathematics, University of Illinois at Urbana-Champaign Urbana, IL 61801, USA
E-mail addresses : [email protected], [email protected] . Date : February 4, 2021.VMH is supported by the NSF through award DMS-2009981. he KdV equation admits solitary and cnoidal waves but no traveling waves can ‘peak’. Bycontrast, the so-called extreme Stokes wave possesses a 120 ◦ peaking at the crest. Also no solutionsof (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 ∗ ofthe 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] provedbreaking, and Ehrnstr¨om and Wahl´en [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 > . . . . , comparing with the well-known Benjamin–Feir instability of a Stokeswave. 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] proposedto replace (2) by(4) V c ww ( | ∂ x | ; T ) f ( k ) = r ( g + T k ) tanh( kh ) k b f ( k ) , where T is the surface tension coefficient. Notice that c ww ( k ; T ) is the phase speed in the lineartheory of capillary-gravity waves [31]. When T = 0, (4) becomes (2). Since c ww ( k ; T ) = p gh (cid:16) − (cid:16) − Tgh (cid:17) k h (cid:17) + O ( k h ) as kh → , one arrives at the KdV equation for capillary-gravity waves:(5) u t + p gh (cid:16) u x + 12 (cid:16) − Tgh (cid:17) h u xxx (cid:17) + 32 r gh uu x = 0in the long wavelength limit unless T /gh = 1 /
3. Notice that (3) and (5) behave alike, qualitatively,possibly after a sign change. By contrast, whenever
T > c ww ( k ; T ) → p T | k | as kh → ∞ , sothat (1) and (4) become(6) u t + √ T | ∂ x | / u x + 32 r gh uu x = 0to leading order whereas when T = 0, u t + √ g | ∂ x | − / u x + q gh uu x = 0.In recent years, the capillary-gravity Whitham equation has been a subject of active research[3, 5, 8, 16, 27] (see also [14, 18, 19, 21, 26]). Particularly, Remonato and Kalisch [27] combined aspectral collocation method and a numerical continuation approach to discover a rich variety oflocal bifurcation of periodic traveling waves of (1) and (4), notably, crossing and connecting solutionbranches. Here we adopt a robust numerical continuation scheme to corroborate the results andproduce 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 withfourth-order time differencing, to experiment with (nonlinear) orbital stability and instability for aplethora of periodic traveling waves of (1) and (4). (Spectral) modulational stability and instabilitywere investigated numerically [3] and also analytically for small amplitude [16]. To the best of theauthors’ knowledge, however, nonlinear stability and instability have not been addressed. Our novelfindings include, among many others, orbital stability for the k = 1 branch whenever T > k > , ∈ N branches for great wave height when 0 < T /gh / ∗ The phase speed is √ gh (1 − k h ), poorly approximating c ww ( k ) when kh ≫ . Preliminaries
We rewrite (1) and (4) succinctly as(7) u t + c ww ( | ∂ x | ; T ) u x + ( u ) x = 0 , where V c ww ( | ∂ x | ; T ) f ( k ) = r (1 + T k ) tanh( k ) k b f ( k ) . Seeking a periodic traveling wave of (7), let u ( x, t ) = φ ( z ) , z = k ( x − ct ) , where c ∈ R ( = 0) is the wave speed, k > φ satisfies, by quadrature andGalilean invariance,(8) − cφ + c ww ( k | ∂ z | ; T ) φ + φ = 0 . We assume that φ is 2 π periodic and even in the z variable, so that 2 π/k periodic in the x variable.Notice(9) c ww ( k | ∂ z | ; T ) (cid:26) cossin (cid:27) ( nz ) = c ww ( kn ; T ) (cid:26) cossin (cid:27) ( nz ) , n ∈ Z . For any T > k > c ∈ R , clearly, φ ( z ) = 0 solves (8) and (9). Suppose that T and k arefixed. A necessary condition for nontrivial solutions to bifurcate from such trivial solution at some c is that c ww ( k | ∂ z | ; T ) φ − cφ = 0 admits a nontrivial solution , if and only if c = c ww ( kn ; T ) for some n ∈ N , by symmetry, whencecos( nz ) ∈ ker( c ww ( k | ∂ z | ; T ) − c ww ( kn ; T )) in some space of even functions . When T = 0, c ww ( · ; 0) monotonically decreases to zero over (0 , ∞ ), whence for any k > c ww ( kn ; T ) > c ww ( kn ; T ) whenever n , n ∈ N and n > n . Therefore, for any k > n ∈ N , ker( c ww ( k | ∂ z | ; 0) − c ww ( kn ; 0)) = span { cos( nz ) } in spaces of even functions.When T > / c ww ( · ; T ) monotonically increases over (0 , ∞ ) and unbounded from above,whereby for any k > n ∈ N , ker( c ww ( k | ∂ z | ; T ) − c ww ( kn ; T )) = span { cos( nz ) } .When T < /
3, to the contrary, c ww ( kn ; T ) = c ww ( kn ; T ) for some k > n , n ∈ N and n = n , provided that T = T ( kn , kn ), where(10) T ( kn , kn ) = 1 kn kn n tanh( kn ) − n tanh( kn ) n tanh( kn ) − n tanh( kn ) , and ker( c ww ( k | ∂ z | ; T ) − c ) = span { cos( n z ) , cos( n z ) } , where c = c ww ( kn ; T ) = c ww ( kn ; T );otherwise, ker( c ww ( k | ∂ z | ; T ) − c ww ( kn ; T )) = span { cos( nz ) } .Suppose that T > k >
0, and T = T ( kn, kn ′ ) for any n, n ′ ∈ N and n = n ′ , particularly, either T = 0 or T > /
3. We assume without loss of generality n = 1. There exists a one-parameter curveof nontrivial, 2 π periodic and even solutions of (8) and (9), denoted by(11) ( φ ( z ; s ) , c ( s )) , | s | ≪ , in some function space (see [8, 9], for instance, for details), and(12) φ ( z ; s ) = s cos( z ) + s (cid:18) c ww ( k ; T ) − c ww (0; T ) − c ww ( k ; T ) − c ww (2 k ; T ) cos(2 z ) (cid:19) + s c ww ( k ; T ) − c ww (2 k ; T ))( c ww ( k ; T ) − c ww (3 k ; T )) cos(3 z ) + O ( s ) ,c ( s ) = c ww ( k ; T ) + s (cid:18) c ww ( k ; T ) − c ww (0; T ) −
12 1 c ww ( k ; T ) − c ww (2 k ; T ) (cid:19) + O ( s ) s | s | →
0. Moreover, subject to a ‘bifurcation condition’ (see [8, 9], for instance, for details), (11)extends to all s ∈ [0 , ∞ ). See [8, 9] and references therein for a rigorous proof.When T = 0 and, without loss of generality, k = 1, ( φ ( z ; s j ) , c ( s j )) → ( φ ( z ) , c ) as j → ∞ forsome s j → ∞ as j → ∞ for some φ ∈ C / ([ − π, π ]) T C ∞ ([ − π, S (0 , π ]) and 0 < c < ∞ such that(13) φ (0) = c . See [9] and references therein for a rigorous proof. Indeed, such a limiting solution enjoys φ ( z ) ∼ c − p π | z | / as | z | → T > /π , so that the Fourier transform of c ww ( k | ∂ z | ; T ) − is ‘completely monotone’ [8],and k = 1, on the other hand, k ( φ ( · ; s ) , c ( s )) k → ∞ as s → ∞ in some † function space; or( φ ( z ; s ) , c ( s )) is periodic in the s variable.See [8], for instance, for a rigorous proof. Our numerical findings suggest min z ∈ [ − π,π ] φ ( z ; s ) → −∞ as s → ∞ . But such a limiting scenario would be physically unrealistic, for the capillary-gravityWhitham equation models water waves in the finite depth. We say that a 2 π periodic solution of(8) and (9) is admissible if φ ( z ) − π ˆ π − π φ ( z ) dz > − z ∈ [ − π, π ] , so that the fluid surface would not intersect the impermeable bed, and we stop numerical continu-ation once we reach a limiting admissible solution , for which(14) min z ∈ [ − π,π ] φ ( z ) − π ˆ π − π φ ( z ) dz = − . Section 4 provides examples.Suppose, to the contrary, that T = T ( kn , kn ) (see (10)) for some k > n , n ∈ N and n < n , so that c ww ( kn ; T ) = c ww ( kn ; T ) =: c . Suppose that n does not divide n . Thereexists a two-parameter sheet of nontrivial, periodic and even solutions of (8) and (9), and(15) φ ( z ; s , s ) = s cos( n z ) + s cos( n z ) + O (( | s | + | s | ) ) ,c ( s , s ) = c + O (( | s | + | s | ) ) ,k ( s , s ) = k + O (( | s | + | s | ) )for | s | , | s | ≪
1. We emphasize that k is a bifurcation parameter. Otherwise, n divides n , and(15) holds for | s | , | s | ≪ s > s > s . In other words, there cannot exist 2 π/n periodic and ‘unimodal’ waves, whose profile monotonically decreases from its single crest to thetrough 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 addressnumerically.Turning the attention to the (nonlinear) orbital stability and instability of a periodic travelingwave of (7), notice that (7) possesses three conservation laws:(16) E ( u ) = ˆ (cid:18) uc ww ( | ∂ x | ; T ) u + 13 u (cid:19) dx (energy) ,P ( u ) = ˆ u dx (momentum) ,M ( u ) = ˆ u dx (mass) , † a H¨older-Zygmund space, for instance [8] nd so does the KdV equation with fractional dispersion:(17) u t + | ∂ x | α u x + ( u ) x = 0 , < α , where \ | ∂ x | α f ( k ) = | k | α b f ( k ) , for which E ( u ) = ´ ( u | ∂ x | α u + u ) dx rather than the first equation of (16). We pause to remarkthat (7) becomes (17), α = 1 /
2, for high frequency (see (6)), and (3) and (5) compare with α = 2. Asolitary wave of (17) is a constrained energy minimizer and orbitally stable, provided that α > / P c >
0. See [24] and references therein for a rigorous proof. See also [17] for ananalogous result for periodic traveling waves. It seems not unreasonable to expect that the orbitalstability 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 ofconstrained energy minimization. Indeed, numerical evidence [20] supports the conjecture when T = 0 and k = 1. Here we take matters further to T > k > , ∈ N . Also we numericallyelucidate the instability scenario when T = 0 and k = 1.3. Methodology
We begin by numerically approximating 2 π periodic and even solutions of (8) and (9) by meansof a spectral collocation method [2, 10]. See [20], among others, for nonlinear dispersive equationsof the form (7) and (17). We define the collocation projection as a discrete cosine transform as(18) φ N ( z ) = N − X n =0 ω ( n ) b φ ( n ) cos( nz ) , where the discrete cosine coefficients are(19) b φ ( n ) = ω ( n ) N X m =1 φ N ( z m ) cos( nz m )and ω ( n ) = (p /N for n = 0 , p /N otherwise;the collocation points are z m = π m − N , m = 1 , , . . . , N. Therefore φ ( z ) ≈ φ N ( z ) = N X m =1 N − X n =0 ω ( n ) cos( nz m ) cos( nz ) φ ( z m ) . We can compute (18) and (19) efficiently using a fast Fourier transform (FFT) [10]. For T > k >
0, likewise, c ww ( k | ∂ z | ; T ) φ ( z ) ≈ ( c ww ( k | ∂ z | ; T ) φ ) N ( z ) := N X ℓ =1 N − X n =0 ω ( n ) c ww ( kn ; T ) cos( nz ) cos( nz ℓ ) φ N ( z ℓ ) , and we can compute ( c ww ( k | ∂ z | ; T ) φ ) N ( z m ), m = 1 , , . . . , N , via an FFT.Suppose that T > k > N = 1024. We numerically solve(20) − cφ N ( z m ) + ( c ww ( k | ∂ z | ; T ) φ ) N ( z m ) + φ N ( z m ) = 0 , m = 1 , , . . . , N, by means of Newton’s method. We parametrically continue the numerical solution over c by meansof a pseudo-arclength continuation method [7] (see [20], among others, for nonlinear dispersiveequations of the form (7) and (17)), for which the (pseudo-)arclength of a solution branch is thecontinuation parameter, whereas a parameter continuation method would use c as the continuationparameter and vary it sequentially. The pseudo-arclength continuation method can successfully ypass a turning point of c , at which a parameter continuation method fails because the Jacobianof (20) becomes singular, whence Newton’s method diverges. See [27] for more discussion.We say that Newton’s method converges if the residual vuut N X m =1 (cid:0) − cφ N ( z m ) + ( c ww ( k | ∂ z | ; T ) φ ) N ( z m ) + φ N ( z m ) (cid:1) − (see (20)), and this is achieved provided that an initial guess is sufficiently close to a true solutionof (8) and (9). To this end, we take a small amplitude cosine function and c ≈ c ww ( k ; T ) as aninitial guess for the local bifurcation at φ = 0 and c = c 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 andrelative tolerance of 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 applicabilityof the code used herein) with a fixed arclength stepsize, to numerically locate and trace solutionbranches. We have additionally run AUTO [6], a robust numerical continuation and bifurcationsoftware, to corroborate the results. All the results in Section 4 have been obtained by employingAUTO with strict tolerances ( EPSL = 1e-12 and
EPSU=1e-12 in the AUTO’s constants file). AUTOhas the advantage, among many others, of making variable arclength stepsize adaptation (using theoption
IADS=1 ) and of detecting branch points (using the option
ISP=2 for all special points). Forthe former, provided with minimum and maximum allowed absolute values of the pseudo-arclengthstepsize, AUTO adjusts the arclength stepsize to be used during the continuation process.Let φ and c numerically approximate a 2 π/k periodic traveling wave of (7), and we turn tonumerically experimenting with its (nonlinear) orbital stability and instability. After a 2048-point(full) Fourier spectral discretization in x ∈ [ − π, π ], (7) leads to(21) u t = L u + N ( u ) , where L = − c ww ( | ∂ x | ; T ) ∂ x and N ( u ) = − ( u ) x , and we numerically solve (21), subject to(22) u ( x,
0) = φ ( kx ) + 10 − max x ∈ [ − π,π ] | φ ( kx ) | U ( − , , by means of an integrating factor (IF) method (see [22], for instance, and references therein), where U ( − ,
1) is a uniform random distribution. In other words, at the initial time, we perturb φ byuniformly distributed random noise of small amplitude, depending on the amplitude ‡ of φ . We pauseto remark that the well-posedness for the Cauchy problem of (7) can be rigorously established atleast for short time in the usual manner via a compactness argument.Following the IF method, let v = e − L t u , so that (21) becomes(23) v t = e − L t N ( e L t v ) . 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-stageRunge–Kutta (RK4) method [11] with a fixed time stepsize ∆ t : v n +1 = v n + a + 2 b + 2 c + d , where v n = v ( t n ) and t n = n ∆ t , n = 0 , , , . . . , a = f ( v n , t n )∆ t, b = f ( v n + a/ , t n + ∆ t/ t,c = f ( v n + b/ , t n + ∆ t/ t, d = f ( v n + c, t n + ∆ t )∆ t ‡ We treat k φ k L ∞ as the amplitude of φ although φ is not of mean zero, because the difference is insignificant. nd f ( v, t ) = e − L t N ( e L t v ) (see (23)). We take ∆ t = 10 − × π/c , where c is the wave speed of φ . While such value of ∆ t ensures numerical stability, some computations, depending on c , requiresmaller § values, for instance, ∆ t = 10 − × π/c or 10 − × π/c . At each time step, we removealiasing errors by applying the so-called 3 / 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 − and M ( t ) to machine precision. Wehave also corroborated our numerical results using two other methods: a higher-order Runge–Kuttamethod, such as the Runge–Kutta–Fehlberg (RKF) method [11] with time stepsize adaptationand a two-stage (and, thus, fourth-order) Gauss–Legendre implicit ¶ Runge–Kutta (IRK4) method[12]. For stable solutions, the results obtained from the RK4 method correctly match those fromthe RKF and IRK4 methods. However, and for unstable solutions, the instability manifests atslightly different times for the same initial conditions. This is somewhat expected because thelocal truncation error (LTE) of each method and the convergence error of the IRK4 method act asperturbations. Recall that the RK4 and IRK4 methods have an LTE of the order of (∆ t ) whereasthe RKF method has (∆ t ) . 4. Results
We begin by taking T = 0 and k = 1. Figure 1 shows the wave height(24) H = max z ∈ [ − π,π ] φ ( z ) − min z ∈ [ − π,π ] φ ( z )and the momentum(25) P = 12 ˆ π − π φ ( z ) dz from our numerical continuation of 2 π periodic and even solutions of (8) and (9). The resultagrees, qualitatively, with [20, Figure 4] and others. We find that H monotonically increases to ≈ . P increases and thendecreases, making one turning point. Also we find that the crest becomes sharper and the troughflatter 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 travelingwaves 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 turningpoint of P , remains unchanged for 10 time periods, after translation of the x axis, implyingorbital stability, whereas the inset reveals that wave (b), past the turning point, suffers from crestinstability. Indeed, our numerical experiments point to transition from stability to instability atthe turning point of P . But there is numerical evidence [28] that waves (a) and (b) are (spectrally)modulationally unstable.When T = 4 /π and k = 1, Figure 3 shows H and P versus c . Our numerical findings suggestthat H, P → ∞ monotonically. Also min z ∈ [ − π,π ] φ ( z ) → −∞ . The solution branch discontinuesonce (14) holds, highlighted with the red square, though, because the solution would be physically § In most cases, 10 − × π/c = O (10 − ) but in Figure 7(b), for instance, c = 0 .
36, whence 10 − × π/c ≈ . − × π/c = O (10 − ). ¶ An implicit method, by construction, enjoys wider regions of absolute stability. This allows us to assess all ourresults, while avoiding numerical instability that might be observed in explicit methods such as RK4 and RKFmethods. However, the IRK4 method requires a fixed point iteration (in the form of Newton’s method) at each timestep. nrealistic, for the capillary-gravity Whitham equation models water waves in the finite depth. Ournumerical continuation works well past the limiting admissible solution, nevertheless. We find thatthe crest becomes wider and flatter, the trough narrower and more peaked, as H increases. But allsolutions must be smooth [8].The left panel of Figure 4 shows an example wave. The right panel shows the profiles perturbedby small random noise at t = 0 and of the solution of (7) after 10 time periods, after translationof the x axis. Our numerical experiments suggest orbital stability for all wave height. But there isnumerical evidence of modulational instability when H is large [3].Numerical evidence [3, 27] points to qualitatively the same results whenever T > / T = T (1 , ≈ . φ = 0 and c = c ww (1 , T ) = c ww (2 , T ) ≈ . H versus c for the k = 1 and k = 2branches, all the way up to the limiting admissible solutions. There are no k 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 π periodic andunimodal waves, whose profile monotonically decreases from a single crest to the trough over theperiod [8]. See also Section 2. For small wave height, the fundamental mode seems dominant, sothat there is one crest over the period = 2 π , 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 separatingthe troughs becomes higher. See the left of Figure 6(b). Observe, to the contrary, π periodic andunimodal waves in the k = 2 branch. See the left of Figure 6(c,d). We find that the crests becomewider and flatter, the troughs narrower and more peaked, as H increases in the k = 1 and k = 2branches. See the left of Figure 6(b,d).Our numerical experiments suggest orbital stability for the k = 1 branch (see the right panels ofFigure 6(a,b)) versus instability for k = 2 (see the right of Figure 6(c,d)).We take matters further to T = T (1 , k = 1 and k = 3 branches are similar to those when T = T (1 ,
2) and k = 1 ,
2. We pause to remark that inthe k = 1 branch, for small H , a smaller wave breaks up the trough into two, and a much smallerwave breaks up the crest into two. See the left panel of Figure 8(a). As H increases, the waveseparating the crests becomes lower, transforming into one wide and flat crest, whereas the waveseparating the troughs become higher (see the left of Figure 8(b)), whereby resembling those when T = T (1 ,
2) and k = 1. Observe π periodic and unimodal waves in the k = 2 branch, orbitallystable for small H versus unstable for large H . See Figure 9(e,f).When T = T (2 , π and 2 π/ k = 2and 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 andmore peaked, as H increases. See the left column of Figure 11. Our numerical experiments (seethe right of Figure 11) suggest orbital instability for the k = 2 branch for large H and for k = 3for all H .When T = T (2 , k = 2 branch. The local bifurcation theorem [8] dictates π periodic and unimodal waves, but wenumerically find that they are not in the k = 2 branch. For small H , for wave (a), for instance, asmaller wave breaks up the trough into two over the half period. As H increases, for wave (b), forinstance, the troughs become narrower and more peaked. Our numerical experiments (see the right k We observe a turning point of P for greater wave height, but the solution is inadmissible. f Figures 13 and 14) suggest orbital stability for the k = 2 branch for small H versus instabilityfor k = 2 for large H and for k = 5 for all H .To proceed, when T = T (1 , k = 1 branch lies above and to the right of the k = 4 branch atleast for small H (not shown), but as T increases, c ww (4; T ) increases more rapidly than c ww (1; T ),so that when T = T (1 , . k = 1 branch crosses the k = 4 branch. Figure 15shows H and P versus c in the k = 1 and k = 4 branches for all admissible solutions. The smallheight result agrees with [27, Figure 10(a)]. We find that in the k = 1 branch, H and P turn toconnect the k = 4 branch, whereas in the k = 4 branch, H and P monotonically increase to thelimiting 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 π 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 = 1and k = 4 branches almost connect, wave (e) in the k = 1 branch and wave (f) for k = 4 are almostthe same. In the k = 4 branch, wave (g), for instance, is π/ 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 . 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 thatthe k = 1 branch crosses and connects the k = 5 branch, like when T = T (1 ,
4) + 0 . k = 1 ,
4, but it continues after connecting all the wave up to the limiting admissible solution. Seethe insets. The left panels of Figures 20, 21 and 22 show several profiles along the k = 1 and k = 5branches. The results for the k = 1 branch up to connecting and for k = 5 are similar to thosewhen T = T (1 ,
4) + 0 . k = 1 ,
4. In the k = 1 branch, after it connects the k = 4 branchat 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 orbitalstability for the k = 1 branch for all H and for k = 5 for small H , versus instability for k = 5 forlarge H .We emphasize orbital stability for the k = 1 branch for all T >
Discussion
Here we employ efficient and highly accurate numerical methods for computing periodic travelingwaves of (7) and experimenting with their (nonlinear) orbital stability and instability. Our findingssuggest, among many others, stability whenever
T > k = 1 branch versus instability when0 < T < / k > , ∈ N branches at least for large wave height. Currently under investigationis to take matters further to all T > k , to classify the orbital stability and instability ofall periodic traveling waves. It will be interesting to devise other numerical continuation methods,for instance, deflated continuation techniques [1, 4], for detecting disconnected solution branchesand others. Also interesting will be to numerically investigate (spectral) modulational stability andinstability of orbitally stable waves. Our methodology will be useful for exploring the nonlineardynamics of modulationally unstable waves. Also it can help tackle the capillary-gravity waveproblem and other nonlinear dispersive equations. Of course, it is of great importance to rigorouslyprove the numerical results. cknowledgement. The authors are grateful to Henrik Kalisch for helpful discussions. EGCacknowledges the hospitality of the Department of Mathematics at the University of Illinois atUrbana-Champaign where early stages of this work took place.
References
1. N. Boull´e, E. G. Charalampidis, P. E. Farrell, and P. G. Kevrekidis,
Deflation-based identification of nonlinearexcitations of the three-dimensional Gross-Pitaevskii equation , Phys. Rev. A (2020), no. 5, 053307, 8 pp.–53314. MR 41900032. John P. Boyd,
Chebyshev and Fourier spectral methods , second ed., Dover Publications, Inc., Mineola, NY, 2001.MR 18740713. John D Carter and Morgan Rozman,
Stability of periodic, traveling-wave solutions to the capillary Whithamequation , Fluids (2019), no. 1, 58.4. E. G. Charalampidis, N. Boull´e, P. E. Farrell, and P. G. Kevrekidis, Bifurcation analysis of stationary solutions oftwo-dimensional coupled Gross-Pitaevskii equations using deflated continuation , Commun. Nonlinear Sci. Numer.Simul. (2020), 105255, 24. MR 41019685. Evgueni Dinvay, Daulet Moldabayev, Denys Dutykh, and Henrik Kalisch, The Whitham equation with surfacetension , Nonlinear Dynam. (2017), no. 2, 1125–1138. MR 36283766. Eusebius Doedel, AUTO , http://indy.cs.concordia.ca/auto/ .7. Eusebius Doedel and Laurette S. Tuckerman (eds.), Numerical methods for bifurcation problems and large-scaledynamical systems , The IMA Volumes in Mathematics and its Applications, vol. 119, Springer-Verlag, New York,2000. MR 17683548. Mats Ehrnstr¨om, Mathew A. Johnson, Ola I. H. Maehlen, and Filippo Remonato,
On the Bifurcation Diagramof the Capillary–Gravity Whitham Equation , Water Waves (2019), no. 2, 275–313. MR 41768709. Mats Ehrnstr¨om and Erik Wahl´en, On Whitham’s conjecture of a highest cusped wave for a nonlocal dispersiveequation , Ann. Inst. H. Poincar´e Anal. Non Lin´eaire (2019), no. 6, 1603–1637. MR 400216810. David Gottlieb and Steven A. Orszag, Numerical analysis of spectral methods: theory and applications , Society forIndustrial and Applied Mathematics, Philadelphia, Pa., 1977, CBMS-NSF Regional Conference Series in AppliedMathematics, No. 26. MR 052015211. E. Hairer, S. P. Nø rsett, and G. Wanner,
Solving ordinary differential equations. I , second ed., Springer Seriesin Computational Mathematics, vol. 8, Springer-Verlag, Berlin, 1993, Nonstiff problems. MR 122798512. E. Hairer and G. Wanner,
Solving ordinary differential equations. II , second ed., Springer Series in ComputationalMathematics, vol. 14, Springer-Verlag, Berlin, 1996, Stiff and differential-algebraic problems. MR 143950613. Vera Mikyoung Hur,
Wave breaking in the Whitham equation , Adv. Math. (2017), 410–437. MR 368267314. ,
Shallow water models with constant vorticity , Eur. J. Mech. B Fluids (2019), 170–179. MR 390748115. Vera Mikyoung Hur and Mathew A. Johnson, Modulational instability in the Whitham equation for water waves ,Stud. Appl. Math. (2015), no. 1, 120–143. MR 329887916. ,
Modulational instability in the Whitham equation with surface tension and vorticity , Nonlinear Anal. (2015), 104–118. MR 341492217. ,
Stability of periodic traveling waves for nonlinear dispersive equations , SIAM J. Math. Anal. (2015),no. 5, 3528–3554. MR 339742918. Vera Mikyoung Hur and Ashish Kumar Pandey, Modulational instability in the full-dispersion Camassa-Holmequation , Proc. A. (2017), no. 2203, 20171053, 18. MR 368546919. ,
Modulational instability in a full-dispersion shallow water model , Stud. Appl. Math. (2019), no. 1,3–47. MR 389726220. Henrik Kalisch, Daulet Moldabayev, and Olivier Verdier,
A numerical study of nonlinear dispersive wave modelswith SpecTraVVave , Electron. J. Differential Equations (2017), Paper No. 62, 23. MR 362594221. Henrik Kalisch and Didier Pilod,
On the local well-posedness for a full-dispersion Boussinesq system with surfacetension , Proc. Amer. Math. Soc. (2019), no. 6, 2545–2559. MR 395143122. Aly-Khan Kassam and Lloyd N. Trefethen,
Fourth-order time-stepping for stiff PDEs , SIAM J. Sci. Comput. (2005), no. 4, 1214–1233. MR 214348223. C. T. Kelley, Solving nonlinear equations with Newton’s method , Fundamentals of Algorithms, vol. 1, Society forIndustrial and Applied Mathematics (SIAM), Philadelphia, PA, 2003. MR 199838324. Felipe Linares, Didier Pilod, and Jean-Claude Saut,
Remarks on the orbital stability of ground state solutions offKdV and related equations , Adv. Differential Equations (2015), no. 9-10, 835–858. MR 336039325. Daulet Moldabayev, Henrik Kalisch, and Denys Dutykh, The Whitham equation as a model for surface waterwaves , Phys. D (2015), 99–107. MR 3390078
6. Ashish Kumar Pandey,
The effects of surface tension on modulational instability in full-dispersion water-wavemodels , Eur. J. Mech. B Fluids (2019), 177–182. MR 395264427. Filippo Remonato and Henrik Kalisch, Numerical bifurcation for the capillary Whitham equation , Phys. D (2017), 51–62. MR 360620028. Nathan Sanford, Keri Kodama, John D. Carter, and Henrik Kalisch,
Stability of traveling wave solutions to theWhitham equation , Phys. Lett. A (2014), no. 30-31, 2100–2107. MR 322608429. J. Sullivan, E. G. Charalampidis, J. Cuevas-Maraver, P. G. Kevrekidis, and N. I. Karachalios,
Kuznetsov–Mabreather-like solutions in the Salerno model , The European Physical Journal Plus (2020), no. 7, 607.30. G. B. Whitham,
Variational methods and applications to water waves , Proceedings of the Royal Society of London.Series A, Mathematical and Physical Sciences (1967), no. 1456, 6–25.31. G. B. Whitham,
Linear and nonlinear waves , Pure and Applied Mathematics (New York), John Wiley & Sons,Inc., New York, 1999, Reprint of the 1974 original, A Wiley-Interscience Publication. MR 1699025 .76 0.78 0.8 0.82 0.84 0.86 0.8800.20.40.6 0.76 0.78 0.8 0.82 0.84 0.86 0.8800.020.040.060.08 Figure 1. T = 0, k = 1. Wave height (left, see (24)) and momentum (right, see(25)) vs. wave speed. The red square corresponds to the limiting solution (see (13)),for which c ≈ . - - /2 0 /2-0.2-0.100.10.20.3 - - /2 0 /2-0.2-0.100.10.20.3 - - /2 0 /2-0.200.20.4 -0.1 0 0.10.20.3 - - /2 0 /2-0.200.20.40.6 Figure 2. T = 0, k = 1. Left column: almost limiting waves at the points labelledwith (a) and (b) in the inset of the right panel of Figure 1, prior to (a) and past (b)the turning point of P , for which c = 0 . − × k φ k L ∞ at t = 0 (dash-dotted),and of the solutions of (7) at later instants of time (solid, see the legends), aftertranslation of the x axis (a). The insets in the bottom panels are close-ups near thecrests. Figure 3. T = 4 /π , k = 1. H (left) and P (right) vs. c . The red square corre-sponds to the limiting admissible solution (see (14)), for which c ≈ . - - /2 0 /2-0.4-0.200.20.40.6 - - /2 0 /2-0.4-0.200.20.40.6 Figure 4. T = 4 /π , k = 1. On the left, the profile at the point labelled with (a)in Figure 3, for which c = 1 .
5, and on the right, the profiles perturbed by uniformlydistributed random noise of amplitude 10 − × k φ k L ∞ at t = 0 (dash-dotted) and ofthe solution of (7) at t = 10 × π/c (solid), after translation of the x axis. Figure 5. T = T (1 .
2) (see (10)). H vs. c for k = 1 (blue) and k = 2 (orange). Thered squares correspond to the limiting admissible solutions for the k = 1 and k = 2branches, for which c ≈ . c ≈ . - /2 0 /2-2024 10 -3 - - /2 0 /2-2024 10 -3 - - /2 0 /2-1.2-0.9-0.6-0.30 - - /2 0 /2-1.2-0.9-0.6-0.30- - /2 0 /2-6-3036 10 -3 - - /2 0 /2-1.2-0.600.61.2 10 -2 - - /2 0 /2-1.2-0.8-0.400.4 - - /2 0 /2-1.6-1.2-0.8-0.400.40.8 Figure 6. T = T (1 , k = 1 ,
2. Left column: profiles at the points labelled with(a),(b) ( k = 1) and (c),(d) ( k = 2) in the insets of Figure 5, for which c = 0 .
973 (a),0 .
97 (c), and 0 .
36 (b,d). Right column: profiles perturbed by small random noiseat t = 0 (dash-dotted) and of the solutions at later instants of time (solid, see thelegends), after translation of the x axis (a,b). Notice different vertical scales. .2 0.4 0.6 0.8 100.20.40.60.81 Figure 7. T = T (1 , H vs. c for k = 1 (blue), k = 3 (yellow), and k = 2(orange). The red squares correspond to the limiting admissible solutions, for which c ≈ . k = 1), c ≈ . k = 3), and c ≈ . k = 2).See Figures 8 and 9 for the profiles at the points labelled with (a)-(f) in the insets. - - /2 0 /2-6-4-202 10 -2 - - /2 0 /2-6-4-202 10 -2 - - /2 0 /2-1.2-0.8-0.400.4 - - /2 0 /2-1.2-0.8-0.400.4 Figure 8.
T=T(1,3), k = 1. Left column: profiles at the points labelled with (a)and (b) in the insets of Figure 7, for which c = 0 .
92 and 0 .
35. Right column: profilesperturbed by small random noise at t = 0 (dash-dotted) and of the solutions after10 periods (solid), after translation of the x axis (a). - /2 0 /2-8-404 10 -2 - - /2 0 /2-8-4048 10 -2 - - /2 0 /2-1.2-0.8-0.400.4 - - /2 0 /2-8-4048 10 -2 - - /2 0 /2-2-1012 10 -2 - - /2 0 /2-2-1012 10 -2 - - /2 0 /2-1.2-0.8-0.400.4 - - /2 0 /2-1.2-0.9-0.6-0.300.30.6 Figure 9. T = T (1 , k = 3 (c,d) and k = 2 (e,f), c = 0 .
92 (c,e) and 0 .
35 (d,f). On the right, solid curves are the solution profiles atlater times (see the legends), after translation of the x axis (e). .2 0.4 0.6 0.8 100.40.81.2 Figure 10. T = T (2 , H vs. c for k = 2 (blue) and k = 3 (orange). The redsquares correspond to the limiting admissible solutions, for which c ≈ . k = 2) and c ≈ . k = 3). See Figure 11 for the profiles at the pointslabelled with (a) and (b) in the inset. - - /2 0 /2-1.2-0.8-0.400.4 - - /2 0 /2-1.2-0.9-0.6-0.300.30.6 - - /2 0 /2-1.2-0.8-0.400.4 - - /2 0 /2-1.2-0.9-0.6-0.300.30.6 Figure 11. T = T (2 , k = 2) and (b) ( k = 3) in Figure 10, for which c = 0 .
3. Right column: profilesperturbed by small random noise at t = 0 (dash-dotted) and of the solutions at latertimes (solid, see the legends). .2 0.4 0.6 0.8 100.20.40.60.811.2 Figure 12. T = T (2 , H vs. c for k = 2 (blue) and k = 5 (inset, orange).The red squares correspond to the limiting admissible solutions, for which c ≈ . k = 2) and c ≈ . k = 5). See Figures 13 and 14 forthe profiles at the points labelled with (a)-(d). - - /2 0 /2-1012 10 -2 - - /2 0 /2-1012 10 -2 - - /2 0 /2-0.8-0.6-0.4-0.200.2 - - /2 0 /2-1.2-0.9-0.6-0.300.30.6 Figure 13. T = T (2 , k = 2. Left column: profiles at the points labelled with(a) and (b) in Figure 12, for which c = 0 .
813 (a) and 0 . t = 0 (dash-dotted) and of the solutions at latertimes (solid), after translation of the x axis (a). - /2 0 /2-8-4048 10 -2 - - /2 0 /2-0.8-0.6-0.4-0.200.2 - - /2 0 /2-1.2-0.9-0.6-0.300.30.6 Figure 14. T = T (2 , k = 5, c = 0 .
805 (c) and 0 . -2 -5 Figure 15. T = T (1 ,
4) + 0 . H (left) and P (right) vs. c for k = 1 (blue) and k = 4 (orange). The red square in the insets corresponds to the limiting admissiblesolution for k = 4, for which c ≈ . - /2 0 /2-10123 10 -3 - - /2 0 /2-10123 10 -3 - - /2 0 /2-8-4048 10 -3 - - /2 0 /2-8-4048 10 -3 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 Figure 16. T = T (1 ,
4) + 0 . k = 1. Left column: profiles at the pointslabelled with (a)-(d) in the right panel of Figure 15, prior to (a) and past (b)crossing the k = 4 branch, and prior to (c) and past (d) the turning point of P , forwhich c = 0 . .
939 (b), and 0 . t = 0 (dash-dotted) and of the solutions after 10 periods(solid), after translation of the x axis. - /2 0 /2-6-3036 10 -3 - - /2 0 /2-6-3036 10 -3 Figure 17. T = T (1 ,
4) + 0 . k = 1. Similar to Figure 16 but almostconnecting the k = 4 branch, for which c ≈ . - - /2 0 /2-6-3036 10 -3 - - /2 0 /2-6-3036 10 -3 - - /2 0 /2-1.2-0.8-0.400.4 - - /2 0 /2-1.2-0.8-0.400.4 Figure 18. T = T (1 ,
4) + 0 . k = 4, almostconnecting the k = 1 branch (f) and almost the end of the numerical continuation(g), for which c ≈ . c = 0 . .9282 0.9284 0.9286 0.9288 0.92900.511.52 10 -2 -5 Figure 19. T = T (1 ,
5) + 0 . H and P vs. c for k = 1 (blue) and k =5 (orange). The red squares in the insets correspond to the limiting admissiblesolutions, for which c ≈ . k = 1) and c ≈ . k = 5). SeeFigures 20, 21, 22 for the profiles at the points labelled with (a) through (i). - - /2 0 /2-4-2024 10 -3 - - /2 0 /2-4-2024 10 -3 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 Figure 20. T = T (1 ,
5) + 0 . k = 1. Left column: profiles at the pointslabelled with (a) and (b) in the right panel of Figure 19, prior to (a) and past (b)crossing the k = 5 branch, for which c = 0 . . t = 0 (dash-dotted) and of the solutionsafter 10 periods (solid), after translation of the x axis (a). - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.8-0.6-0.4-0.200.2 - - /2 0 /2-1-0.8-0.6-0.4-0.200.2 Figure 21. T = T (1 ,
5) + 0 . k = 1. Similar to Figure 20, almost connectingthe k = 5 branch (c), prior to (d) and past (e) crossing the k = 1 branch itself, andalmost the end of the numerical continuation (f), for which c = 0 . . . . - /2 0 /2-6-3036 10 -3 - - /2 0 /2-6-3036 10 -3 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-1-0.500.51 10 -2 - - /2 0 /2-0.9-0.6-0.300.3 - - /2 0 /2-1.2-0.9-0.6-0.300.30.6 Figure 22. T = T (1 ,
5) + 0 . k = 5. Left column: profiles at the pointslabelled with (g)-(i) in the right panel of Figure 20, almost crossing the k = 1 branch(g), almost connecting the k = 1 branch (h), and almost the end of the numericalcontinuation (i), for which c = 0 . . . t = 0 (dash-dotted) and of thesolutions at lager times (solid).= 0 (dash-dotted) and of thesolutions at lager times (solid).