ON THE EXISTENCE AND COMPUTATION OF PERIODIC TRAVELLING WAVES FOR A 2D WATER WAVE MODEL

. In this work we establish the existence of at least one weak periodic solution in the spatial directions of a nonlinear system of two coupled diﬀerential equations associated with a 2D Boussinesq model which describes the evolution of long water waves with small amplitude under the eﬀect of sur- face tension. For wave speed 0 < | c | < 1, the problem is reduced to ﬁnding a minimum for the corresponding action integral over a closed convex subset of the space H 1 k ( R ) ( k -periodic functions f ∈ L 2 k ( R ) such that f (cid:48) ∈ L 2 k ( R )). For wave speed | c | > 1, the result is a direct consequence of the Lyapunov Center Theorem since the nonlinear system can be rewritten as a 4 × 4 system with a special Hamiltonian structure. In the case | c | > 1, we also compute numerical approximations of these travelling waves by using a Fourier spectral discretization of the corresponding 1D travelling wave equations and a Newton-type iteration.

1. Introduction. In this paper, we study the existence of non-trivial periodic travelling waves in both spatial variables for the Boussinesq type system derived in [11] to describe the evolution of three dimensional long water waves with small amplitude: where is the amplitude parameter (nonlinearity coefficient), µ = (h 0 /L) 2 (with h 0 and L being the characteristic depth and characteristic length, respectively) is the long-wave parameter (dispersion coefficient), σ −1 is the Bond number (associated with the surface tension), and p is a rational number. The variable Φ is the rescaled nondimensional velocity potential on the bottom z = 0, and the variable η is the rescaled free surface elevation.
One of the main features of this model analyzed by J. Quintero in [11] is that the generalized Benney-Luke equation and the generalized Kandomtsev-Petviashvili equation can be obtained formally from this Boussinesq-type model (1). J. Quintero also proved in [11] that in an appropriate scaling there are sequences of solitons for this Boussinesq system that after a translation converge as , µ → 0 to a soliton of the generalized Kandomtsev-Petviashivili equation, showing that (1) has physically meaningful solitons. Further in [12], the existence and uniqueness of solutions of the Cauchy problem associated to system (1) is considered using Kato's approach [5,4,6] by exploiting that this system has a 2D KdV structure. Stability of solitary waves of system (1) is also analyzed in [12]. Other advantage of this bi-directional system is that equations in system (1) are not reduced to a unique scalar equation for the potential velocity Φ by eliminating the wave elevation η. This prevents us to neglect additional nonlinear and dispersive terms as in the scalar water wave models derived in [9,15] and [10], for instance. Furthermore, J. Quintero also considered in [13] the existence of x-periodic travelling wave solutions of system (1) by using a variational approach for wave speed |c| < 1.
J. Quintero in [14] in the regime |c| > 1 and σ > 1 2 (big enough) described all small waves η(x, y, t) = u(x − ct, y), Φ(x, y, t) = v(x − ct, y) (that translate steadily with supercritical speed c and that are periodic in the direction of translation (or orthogonal to it). In this regime the travelling-wave equations (after rescaling , µ) for (1) satisfied by u, v have mixed type, when looking for x-periodic travelling waves as solutions for the initial-value problem for (3) considered as an evolution equation in y. The associated Cauchy problem turns out to be linearly ill-posed. More exactly, at the linear level there are finite many central modes (pure imaginary eigenvalues) and infinitely many hyperbolic modes. So, the result involves using an invariant center manifold of finite dimension and infinite codimension. It was also shown that this center manifold contains all globally defined small-amplitude solutions of the travelling wave equation for the Boussinesq system that are periodic in the direction of propagation. In this paper, in the case of wave speed |c| < 1 and σ > 1 2 , we use the Arzelà-Ascoli Theorem and the characterization of bounded functional on convex sets to establish the existence of periodic solutions of the form for the travelling wave system (3), where ψ and ϕ are periodic functions of the variable ξ which must satisfy the system The existence of periodic solutions of system (5)-(6) is reduced to finding a minimum for the corresponding action integral over a closed convex subset of the space H 1 k (R) (k-periodic functions f ∈ L 2 k (R) such that f ∈ L 2 k (R)). On the other hand, in the regime in |c| > 1 and σ > 1 2 , we use the Lyapunov Center Theorem [8] (chapter 9) to establish the existence of periodic solutions for system (5)- (6). Furthermore, an explicit expression for the period of such solutions is obtained.
It is important to point out that we already obtain periodic solutions u and v for the travelling wave system (3) in both spatial variables (not considered in previous works on system (1)). In fact, let T be the period of the function ψ. Then, if we choose T 1 and T 2 such that T = T 1 + βT 2 , then u (x, y) = ψ (ξ) and v (x, y) = ϕ (ξ) with ξ = x + βy are periodic functions of period (T 1 , T 2 ). In fact, u(x+T 1 , y+T 2 ) = ψ(x+T 1 +β(y+T 2 )) = ψ(x+βy+T 1 +βT 2 ) = ψ(x+βy) = u(x, y), and analogously for the component v.
The second aim of the present paper is to compute some numerical approximations of solutions of system (5)- (6) in the case that |c| > 1, p = 1 by using a Newton-type iteration. We explain in detail how to compute an appropriate initial step from an exact periodic solution of a KdV equation formally equivalent (up to some order in , µ) to system (1). The accuracy of the travelling wave solutions computed is checked by means of an implicit-explicit (IMEX) second-order strategy for time stepping of the corresponding initial value problem. This type of temporal discretization is described for instance in [1] and it has been used in conjunction with spectral methods [7,2], for the time integration of spatially discretized PDEs of diffusion-convection type. IMEX schemes have also been successfully applied to the incompressible Navier-Stokes equations [3] and in environmental modelling studies [16].
The plan of the rest of this paper is now outlined. Section 2 is devoted to discuss existence of periodic travelling waves in the form (2), (4) of system (1) for wave speed 0 < |c| < 1 via the Arzelà-Ascoli Theorem together with a variational strategy applied to the corresponding action integral. In Section 3 for the complementary regime |c| > 1, we use Lyapunov Center Theorem and the Hamiltonian structure of the travelling wave system to establish existence of such type of solutions of system (1). Finally in Section 4, we present the Galerkin-Fourier spectral numerical schemes employed for approximating a family of solutions of systems (5)- (6) and (1). Some numerical simulations are conducted for particular values of the modelling parameters.
2. Existence of periodic travelling waves via the Arzelà-Ascoli Theorem for 0 < |c| < 1. Hereafter we set L k = [− k 2 , k 2 ] and define the space V k as the closure of the C ∞ k (R) (periodic C ∞ functions of period k) with respect to the norm given by We define the functional space as . Note that V k and X k are Hilbert spaces with inner products given respectively by In inequalities below, C denotes a generic constant whose value may change from instance to instance. Hereafter p will denote a rational number with odd numerator and denominator. More concretely, p = p1 p2 with p i being odd positive integer and gcd(p 1 , p 2 ) = 1.
As we mention in the introduction, we are looking for periodic solutions of system (3) of the form with functions ψ and ϕ being periodic in the variable ξ. A direct computation shows that weak k-periodic solution of system (5)-(6) can be characterized as critical points for the functional where functional I k , G 1,k and G 2,k are defined as It is easy to see that the functionals I k , G 1,k and G 2,k are smooth maps from X k to R. For instance, if f ∈ H 1 k has mean zero in [0, k], then for q ≥ 1 we have that On the other hand, from Cauchy's and Young's inequalities, we get that since ϕ ∈ H 1 k has trivially mean zero on [0, k]. Now, if we set Σ k = I k + G 1,k , we have from (8) that Thus, for σ > 1/2, β ∈ R and 0 < |c| < 1 there is some positive constant C 1 (σ, c, β) > 1 such that We note that critical points of the functional J c,k satisfy the travelling wave system (5)- (6). Hereafter, we will say that weak solutions for (5)-(6) are critical points of the functional J c,k . A direct computation shows that Thus on any critical point w, we have that Lemma 2.1. Assume that the sequence (ψ n , ϕ n ) n ⊂ X k converges weakly in X k to (ψ 0 , ϕ 0 ) ∈ X k . Then there is a subsequence (ψ nj , ϕ nj ) j of (ψ n , ϕ n ) n such that we have that lim inf Proof. Recall that J c,k = Σ k + G 2,k and that Σ k is like a norm in X k , so is convex. More exactly, for λ ∈ (0, 1) and any u, v ∈ X k , we have that On the other hand, we have that From previous remark, we have that Note using that the sequences (ψ n ) n and (ϕ n ) n converge weakly in H 1 k to ψ 0 and In other words, we have that Now, we need to observe that On the other hand, since H 1 k ⊂ L q k is compact for q ≥ 2, then there is a subsequence (ψ nj , ϕ nj ) j of (ψ n , ϕ n ) n such that ψ nj → ψ 0 and ϕ nj → ϕ 0 in L q k for q ≥ 2. Moreover, we also have that which means, after recalling that the sequence (ψ nj , ϕ nj ) j is bounded, that As a consequence of previous remarks, we conclude that We must recall that H 1 k is the space of absolutely continuous functions u which are k-periodic and such that u ∈ L 2 loc (R). For a > 0, we consider the weakly closed subset of X k Lemma 2.2. 1.-There are positive constants C 1 and C 2 such that for any Φ ∈ X k , we have that 2.-There exists a 0 > 0 such that for 0 < a < a 0 the functional J c,k is coercive on X k,a . More exactly, there is C 3 > 0 such that for Φ ∈ X k,a , Proof. 1.-From inequalities (9) and (10), there are positive constants C 1 and C 2 such that So, using inequality (10) and previous one, we have that Our goal now is to show the existence of a non trivial critical point for J c,k . The result will be a direct consequence of the coerciveness of J c,k and that J c,k is (sequentially) weakly lower semi-continuous on X k,a for 0 < a < a 0 .
Proof. It is straightforward to check that X k,a is weakly closed subset of X k . In fact, let (ψ n , ϕ n ) n ⊂ X k,a be a sequence that converges weakly to (ψ 0 , ϕ 0 ). Then we have that the sequence (ψ n , ϕ n ) n is bounded in X k . Now, since ϕ n ∈ H 1 k has mean zero on [0, k], we know that From Arzelá-Ascoli Theorem we have for some subsequence (denoted equal) that (ϕ n ) n converges uniformly to ϕ 0 on [0, k], since we have that |ϕ n (ξ)| ≤ a for a.e. ξ ∈ R and for all n ∈ N. From this fact and the uniform convergence of (ϕ n ) n , we conclude that |ϕ 0 (ξ)| ≤ a for a. e. ξ ∈ R. Moreover, as done in the proof of Lemma (2.1), we also have that G 2,k (ψ 0 , ϕ 0 ) = 1. In other words, (ψ 0 , ϕ 0 ) ∈ X k,a , meaning that X k,a is weakly closed subset of X k . On the other hand, from previous Lemma we have that J c,k is coercive and from inequality (21) that J c,k is bounded below on X k,a . So, let (ψ n , ϕ n ) n be a minimizing sequence for J c,k on X k,a such that We easily see from inequality (21) that (ψ n , ϕ n ) n is bounded on X k , and so there (ψ 0 , ϕ 0 ) ∈ X k,a , such that (ψ n , ϕ n ) n (ψ 0 , ϕ 0 ) (weakly) in X k,a (passing to a subsequence, if necessary). Moreover, by Lemma (2.1) we conclude passing to a subsequence, if necessary, that meaning that (ψ 0 , ϕ 0 ) minimizes J c,k over X k,a , as desired.
3. Periodic travelling waves via the Lyapunov Center Theorem for |c| > 1. In this section, we will see that periodic solutions for |c| > 1 can be obtained by using the Lyapunov Center Theorem when we look for solutions of the form where ψ and ϕ are periodic functions of the variable ξ. The result will be a consequence of the Lyapunov Center Theorem which can be obtained as an application of the Hopf Bifurcation Theorem.
where f is a smooth function which vanishes along with its first partial derivatives at x = 0. Assume that the system admits a first integral of the form where S is a n × n real symmetric matrix with det S = 0. Let the matrix A have eigenvalues . . , n, then system (22) has a one parameter family of periodic solutions emanating from the origin with period 2π λ1 .
If we set z = ϕ , we see after integrating that (z, ψ) must satisfy the system So, in the variable V = (z, ψ, ρ, w), where z = ρ and ψ = w, we obtain an appropriate model to apply the Lyapunov Center Theorem where A and F are defined by We see directly that this system has a first integral H (Hamiltonian) given by where S is a 4 × 4 real symmetric matrix with det S = 0 defined as Moreover, system (23) has the Hamiltonian structure where J is the skew symmetric matrix given by, From this fact, we see that if V is a solutions for system (23), then H (V ) is conserved in time. In fact, since J is skew symmetric we have that Now, a direct computation shows that the equilibrium points for system (23) satisfy that We find that z satisfies the equation, Then we have that z = 0 or Clearly, in the regime c 2 > 1 + β 2 , we have that z ± (p) and c have the same sign. In fact, c 2 p 2 + 4 (p + 1) (1 + β 2 ) ≤ |c| (p + 2) , which implies that −|c| (p + 2) ≤ c 2 p 2 + 4 (p + 1) (1 + β 2 ) ≤ |c| (p + 2) .
On the other hand, a direct computation shows that (26), we conclude in the regime c 2 > 1 + β 2 that z + (p) and ψ + (p) have opposed signs and that z − (p) and ψ − (p) have equal signs meaning that where P e is an equilibrium point for system (23), Then we have that Moreover, for P e = 0, we have that 3.1. Eigenvalue analysis for. p ≥ 1. It is clear that in this case, there are three equilibrium points given by where z ± (p) is given by (28).
Spectral Analysis for P 0 (p). In order to compute the eigenvalues at the point P 0 , we must find the eigenvalues of the linearization of system (23) around the point P 0 . In other words, we need to determine the eigenvalues of the matrix A given by . It is straightforward to see that the eigenvalues λ of A must satisfy the equation λ 4 − d + 1 + β 2 a λ 2 + 1 + β 2 − c 2 ad = 0, whose roots in terms of λ 2 are given by

PERIODIC TRAVELLING WAVES FOR A 2D WATER WAVE MODEL 567
Now we have that In this range for the wave speed c, we have that the eigenvalues of A are In terms of σ, β and c, the eigenvalues of A are given by Clearly, λ 1 , λ 2 ∈ iR and λ 3 , λ 4 ∈ R. Moreover, λ3 λ1 / ∈ Z. In other words, we have verified the hypotheses of the Lyapunov Center Theorem for system (23). More concretely, Theorem 3.2. Let σ > 1 2 , c = 0 and β be such that c 2 > 1 + β 2 . Then for p ≥ 1, system (23) has a one parameter family of periodic solutions (z, ϕ, ρ, w) emanating from the origin with period T 0 > 0 defined as Spectral Analysis for P + (1). Now, we turn out the attention in the analysis for the equilibrium point P + (1) in the regime c 2 > 1 + β 2 . We recall that if and only if Then we conclude that where numbers a, b and d are defined as .
4. Numerical results. We wish to analyze the periodic travelling waves for the supercritical case |c| > 1, σ > 1/2, p = 1 studied in the previous sections for the original system from a numerical point of view. In fist place, we will compute some travelling wave solutions not covered by the theoretical results in section 3 since their wave speeds do not satisfy that c 2 > 1 + β 2 , and their corresponding fundamental periods do not coincide with the ones predicted by Theorems 3.2 and 3.3. Assume we look for a solution (η,Φ) of the formη(x, y, t) = η(x + βy, t) and Φ(x, y, t) = Φ(x + βy, t). Then we see directly that (η, ϕ) with Φ = ϕ satisfies the system We point out that the derivatives in the equations above are with respect to the variable x+βy. Since it is difficult to come out with explicit travelling wave solutions for the previous system, in this section we will describe a numerical strategy to approximate these solutions at least in the case that p = 1. First note that we have formally Now, let ρ > 0 be such that ρ 2 (1 + β 2 ) = 1. Thus multiplying equation (38) by ρ 2 , we get So, in order to these equations be consistent to leading order truncation, we may assume that This assumption corresponds to searching waves moving to the right, as mentioned for instance in [17]. Using this, we see from equations (38)-(39) that η satisfies the following equations, Multiplying the first equation by ρ and using that ρ 2 (1 + β 2 ) = 1, we conclude that and therefore we get that Using this system, it follows that Therefore a direct computation shows that Then we obtain that and also that η satisfies, up to order O( 2 , µ 2 ), the normalized GKdV equation which can be transformed in the generalized KdV equation by defining η(x, t) = γw νx − νρ(1 + β 2 )t, νt , where ν and γ are defined as It is straightforward to see that w(x, t) = −v(−x, t) satisfies the generalized KdV We note that for p = 1, the KdV equation has explicit periodic travelling wave solutions v(x, t) = χ(x − c 0 t) of mean zero on [0, T 0 ] of the cnoidal type where k is known as the modulus being defined by and the constants β 1 , β 2 , β 3 and the period T 0 are given by where K is the complete elliptic integral of the first kind given by We see directly that equation (48) Observe that this periodic wave has period and propagates with speed equal to c = ρ(1 + β 2 ) − c 0 = 1 + β 2 − c 0 . We emphasize that the wave speed of the periodic waves computed does not satisfy the requirement c 2 > 1 + β 2 and the fundamental period T does not coincide in general with the ones given in equations (31)-(34).

Approximation of periodic travelling waves.
To improve the accuracy of the approximation given in (51)-(52) of a travelling wave with period T = 2l of system (36)-(37), we correct it by using a Newton's iteration. We recall that a travelling wave solution η(x, t) = ζ(x − ct), ϕ(x, t) = u(x − ct) for system (36)-(37) must satisy the following equations: Let us introduce truncated cosine expansions for ζ and u: where and analogous expressions for u. We point out that this strategy can also be used for approximating solutions decaying to zero at infinity, provided that the period 2l is taken large enough. By substituting expressions (55) into the travelling wave equations (53)-(54) and evaluating them at the N/2 + 1 collocation points we obtain a system of N + 2 nonlinear equations in the form where the N +2 coefficients ζ n , u n are the unknowns. The nonlinear system (57) can be solved by Newton's iteration. Computation of the cosine series in (55) and the integrals in (56) is performed using the FFT (Fast Fourier Transform) algorithm. The Jacobian of the vector field F : R N +2 → R N +2 is approximated by the secondorder accurate formula where e j = (0, ..., 1, ..., 0) and h = 0.01. We stop the iteration procedure when the relative error between two successive approximations and the value of the vector field F are smaller than 10 −12 .

Computation of solutions of the initial value problem.
To check the accuracy of the travelling wave solutions computed above, we describe the numerical scheme for computing spatially periodic solutions with period T of the initial value problem associated to system (36)-(37). In the numerical solver to be introduced, the spatial computational domain [0, T ] is discretized by N equidistant points, with spacing ∆x = T /N , and the unknowns η and ϕ are expanded as truncated Fourier series in space with time-dependent coefficients: with w j = 2πj T , j = −N/2 + 1, ...0, ...N/2.
The time-dependent coefficientsφ j (t), j = −N/2 + 1, ..., 0, ..., N/2 are calculated by means of the equationφ and analogously forη j (t). Substituting these expressions into equations (36)-(37) and projecting the resulting equations with respect to the L 2 -orthonormal basis φ j = T −1/2 e iwj x and the inner product it follows thatη where and P j [.] denotes the operator 4.3. Temporal discretization. Note that equations (60)-(61) can be seen as a system of ordinary differential equations where the unknowns are the time-dependent Fourier coefficients of the solutions. To solve it, we use an implicit-explicit scheme (IMEX) in the form Here η n , ϕ n denote the approximations of the unknowns η(x, t), ϕ(x, t), respectively, at time t = n∆t, where ∆t is the time step of the method and N n 1 , N n 2 denote the approximations of the functions N 1 , N 2 evaluated at time n∆t. Similarly,η n j ,φ n j denote the approximations to the Fourier transforms of the functions η and ϕ, respectively, with respect to the variable x, evaluated at time n∆t. Observe that the linear dispersive terms are approximated by using an implicit strategy, in contrast to the nonlinear terms which are treated in explicit form. Implicit-explicit schemes (IMEX) were already applied for instance in [1] to scalar dispersive evolution equations and to the incompressible Navier-Stokes equations [3]. The main advantage of the numerical scheme described here is that at each time step we can solve explicitly the approximationsη n+1 j ,φ n+1 j , j = −N/2 + 1, ..., N/2 of the Fourier coefficients of the unknowns η(x, t) and ϕ(x, t) from equations (63)-(64), without using implicit Newton-type iterations. Thus, the scheme results to be cheap and its computer implementation is easier.
In the scheme, the spatial derivatives ϕ x and η x are computed by exact differentiation of the truncated Fourier series. For instance, The numerical calculations presented in this paper were carried out in double precision by using Matlab R2012b on a Mac platform. The Fourier-type integral appearing in the operator P j [.] (see equation (62)) is approximated through the well-known Fast Fourier Transform (FFT) routine.

4.4.
Description of the numerical experiments. In this section we set p = 1 in equations (36)-(37). In first place, we compute some travelling wave solutions of this system by correcting the profiles defined by equations (51)-(52) through the Newton's iteration described above using N = 2 7 FFT collocation points. The travelling waves obtained are displayed in Figures 1 and 2 for two different sets of parameters in the initial step (51)-(52). To verify that the profiles obtained correspond really to travelling waves of system (36)-(37), we employ the numerical scheme given in equations (63)-(64) for solving the corresponding initial value problem. The results are also displayed in Figures 1, 2. We point out that the difference in the norm of maximum between the profiles shown in each figure is of order 1e − 4. Thus we observe that the approximate travelling wave solutions propagate without significant change of shape and with the expected wave speed c. We point out that the numerical solver given in (63)-(64) can also been employed to explore the orbital stability of periodic travelling wave solutions of system (36)-(37) in a future research.
Finally in Figures 3 and 4, we further display the surface plots of the wave elevationη(x, y, t) = η(x + βy, t) of the original system (35) evaluated at t = 0, for the parameters used in Figures 1 and 2, respectively.
As mentioned above, the approximations of travelling wave solutions of system (36)-(37) obtained in the previous computer simulations do not satisfy the condition on the wave speed c 2 > 1+β 2 required in Theorems 3.2, 3.3. Therefore our numerical simulations give evidence on the fact that the existence of travelling wave solutions in the supercritical case |c| > 1 is even possible when this technical requirement is not satisfied.
On the other hand, in Figures 5 and 6, we display periodic solutions of system (53)-(54) which are included in the family of solutions guaranteed by the theoretical results presented in section 3. We recall that η(x, t) = ζ(x − ct), ϕ(x, t) = u(x − ct) corresponds to a periodic travelling wave solution of the original system (36)-(37) with speed c. Observe that the periods T 0 and T + (1) of these solutions are in accordance with the predicted values given by equations (31), (34) in Theorems 3.2 and 3.3, respectively. To obtain these travelling wave solutions, we use 2 8 FFT collocation points and the initial data for Newton's procedure are ζ 0 (x) = 1 + cos 2πx 2T 0 , u 0 (x) = 1 + cos 2πx T 0 , and ζ 0 (x) = 1 + cos (2)(3)πx T + (1) , u 0 (x) = 1 + cos (2)(2)πx T + (1) , for the experiments in Figures 5 and 6, respectively.    Figure 4. Surface plot of the wave elevationη(x, y, t) = η(x + βy, t) in the original system (35) at t = 0, with the parameters used in Figure 2. Observe that this solution satisfies the condition on the wave speed c 2 > 1 + β 2 as required in Theorem 3.3.