On Whitham and related equations

TL;DR

This paper investigates the dynamics of Whitham and related equations through theoretical and numerical methods, establishing bounds with KdV solutions and exploring behaviors beyond traditional water wave modeling regimes.

Contribution

It provides rigorous bounds between Whitham and KdV solutions and offers new insights into the equation's dynamics outside standard water wave regimes.

Findings

Established bounds between Whitham and KdV solutions

Numerical simulations reveal dynamics outside typical water wave regimes

Insights into the validity range of the Whitham equation

Abstract

The aim of this paper is to study, via theoretical analysis and numerical simulations, the dynamics of Whitham and related equations. In particular we establish rigorous bounds between solutions of the Whitham and KdV equations and provide some insights into the dynamics of the Whitham equation in different regimes, some of them being outside the range of validity of the Whitham equation as a water waves model.

Tables

Table 1. Table 1. Critical times τ c = ϵ ​ t c subscript 𝜏 𝑐 italic-ϵ subscript 𝑡 𝑐 \tau_{c}=\epsilon t_{c} and exponents according to ( 7.4 ) for the solutions to the Whitham equation for initial data u ​ ( x , 0 ) = 10 ​ exp ⁡ ( − x 2 ) 𝑢 𝑥 0 10 superscript 𝑥 2 u(x,0)=10\exp(-x^{2}) for different values of ϵ italic-ϵ \epsilon .

$ϵ$	${\tau }_c$	$\mu$
1	0.1175	0.3635
0.8	0.1178	0.3649
0.6	0.1184	0.3641
0.4	0.1197	0.3619
0.2	0.1254	0.3499
0.1	0.1482	0.3855
0.08	0.1674	0.4133
0.06	0.2163	0.4274

Equations

$$u_t+{\mathcal{L}}_{ϵ}{u}_x+ϵu{u}_x=0,$$

$${\mathcal{L}}_{ϵ}=l(\sqrt{ϵ}D):={\left(\frac{\tanh \sqrt{ϵ}|D|}{\sqrt{ϵ}|D|}\right)}^{1/2}\text{and}D=-i\nabla =-i\frac{\partial }{\partial x}.$$

$$u_t+u_x+ϵu{u}_x+\frac{ϵ}{6}{u}_{xxx}=0.$$

$${\partial }_{t}u+{\tilde{c}}_{WW}(\sqrt{ϵ}|D^{ϵ}|){\left(1+ϵ\frac{D_2^2}{D_1^2}\right)}^{1/2}{u}_x+ϵ\frac{3}{2}u{u}_x=0,$$

$${\tilde{c}}_{WW}(\sqrt{ϵ}k)=(1+\beta ϵk^2{)}^{\frac{1}{2}}{\left(\frac{\tanh \sqrt{ϵ}k}{\sqrt{ϵ}k}\right)}^{1/2},$$

$$|D^{ϵ}|=\sqrt{D_1^2+ϵD_2^2},D_1=\frac{1}{i}{\partial }_x,D_2=\frac{1}{i}{\partial }_y.$$

$$ϵ\sim \frac{a}{h}\sim {\left(\frac{{\lambda }_x}{{\lambda }_y}\right)}^2\sim {\left(\frac{h}{{\lambda }_x}\right)}^2\ll 1.$$

$${\partial }_{t}u+c_{WW}(\sqrt{ϵ}|D^{ϵ}|){\left(1+ϵ\frac{D_2^2}{D_1^2}\right)}^{1/2}{u}_x+ϵ\frac{3}{2}u{u}_x=0,$$

$$c_{WW}(\sqrt{ϵ}k)={\left(\frac{\tanh \sqrt{ϵ}k}{\sqrt{ϵ}k}\right)}^{1/2},$$

$$u_t+{\tilde{\mathcal{L}}}_{ϵ}{u}_x+ϵu{u}_x=0,$$

$${\tilde{\mathcal{L}}}_{ϵ}=(1+\beta ϵ|D{|}^2{)}^{1/2}{\left(\frac{\tanh \sqrt{ϵ}|D|}{\sqrt{ϵ}|D|}\right)}^{1/2},\beta >0.$$

$$u_t+u_x+\frac{ϵ}{2}\left(\frac{1}{3}-\beta \right)u_{xxx}+ϵu{u}_x=0.$$

$$u_t+u_x+u{u}_x-|D{|}^{\alpha }{u}_x=0.$$

$$\left\{\begin{array}[]{l}\eta_{t}+\mathcal{L}_{\epsilon}^{2}u_{x}+\epsilon(\eta u)_{x}=0\\ u_{t}+\eta_{x}+\epsilon uu_{x}=0,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\eta_{t}+\mathcal{L}_{\epsilon}^{2}\nabla\cdot{\bf u}+\epsilon\nabla\cdot(\eta{\bf u})=0\\ {\bf u}_{t}+\nabla\eta+\frac{\epsilon}{2}\nabla|{\bf u}|^{2}=0,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\eta_{t}+u_{x}+\frac{\epsilon}{3}u_{xxx}+\epsilon(\eta u)_{x}=0\\ u_{t}+\eta_{x}+\epsilon uu_{x}=0,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\eta_{t}+\nabla\cdot{\bf u}+\frac{\epsilon}{3}\Delta\nabla\cdot{\bf u}+\epsilon\nabla\cdot(\eta{\bf u})=0\\ {\bf u}_{t}+\nabla\eta+\frac{\epsilon}{2}\nabla|{\bf u}|^{2}=0,\end{array}\right.$$

$${\partial }_t\left(\begin{array}{c}\eta \\ \mathbf{u}\end{array}\right)+J\text{grad}{H}_{ϵ}(\eta ,\mathbf{u})=0$$

$$J=\left(\begin{array}{ccc}0& {\partial }_x& {\partial }_y\\ {\partial }_x& 0& 0\\ {\partial }_y& 0& 0\end{array}\right),$$

$$H_{\epsilon}(U)=\frac{1}{2}\int_{\mathbb{R}^{2}}\big{(}|\mathcal{L}_{\epsilon}{{\bf u}}|^{2}+\eta^{2}+\epsilon\eta|{{\bf u}}|^{2}\big{)}dxdy,$$

$$U=\left(\begin{array}{c}\eta \\ \mathbf{u}\end{array}\right),$$

$${\partial }_t\left(\begin{array}{c}\eta \\ u\end{array}\right)+J\text{grad}{H}_{ϵ}(\eta ,u)=0$$

$$J=\left(\begin{array}{cc}0& {\partial }_x\\ {\partial }_x& 0\end{array}\right)$$

$$H_{ϵ}(\eta ,u)=\frac{1}{2}{\int }_{\mathbb{R}}(|{\mathcal{L}}_{ϵ}u{|}^2+{\eta }^2+ϵu^2\eta )dx,$$

$$\left\{\begin{array}[]{l}\eta_{t}+u_{x}+\epsilon(\eta u)_{x}=0\\ u_{t}+\mathcal{P}_{\epsilon}\eta_{x}+\epsilon uu_{x}=0,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\eta_{t}+\nabla\cdot{\bf u}+\epsilon\nabla\cdot(\eta{\bf u})=0\\ {\bf u}_{t}+\mathcal{P}_{\epsilon}\nabla\eta+\frac{\epsilon}{2}\nabla|{\bf u}|^{2}=0,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\eta_{t}+u_{x}+\epsilon(\eta u)_{x}=0\\ u_{t}+\eta_{x}+\epsilon(\frac{1}{3}-\beta)\eta_{xxx}+\epsilon uu_{x}=0,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\eta_{t}+\nabla\cdot{\bf u}+\epsilon\nabla\cdot(\eta{\bf u})=0\\ {\bf u}_{t}+\nabla\eta+\epsilon(\frac{1}{3}-\beta)\nabla\Delta\eta+\frac{\epsilon}{2}\nabla|{\bf u}|^{2}=0,\end{array}\right.$$

$$u_t+u^4{u}_x+u_{xxx}=0,$$

$$u_t+u^2{u}_x-\mathcal{H}{u}_{xx}=0,$$

Full text

1. introduction

The Whitham equation was introduced formally in [60] as an alternative to the KdV equation, by keeping the exact dispersion of the water waves system. No rigorous derivation from the water waves system is known, and there is no consistent scaling allowing to derive it.

Actually the idea in [60] was to propose a model describing also the occurrence of waves of greatest height with the Stokes 120 degrees angle at the crest which is clearly impossible with the KdV equation.

As was conjectured in [60] and later proved in [24], this is actually the case (with a different angle though) but this “wave breaking”  phenomena seems to be irrelevant in the (KdV) regime where the Whitham equation can be really viewed as a consistent water wave model, as will be discussed in the present paper.

In short, and as was noticed by Whitham, “the desired qualitative effect is in”  but this wave breaking effect cannot be quantitatively linked with the water waves system and thus to “real”  waves. From the modeling point of view, the advantage of the Whitham equation with respect to KdV seems to be the enlargement of the range of validity of this asymptotic model in terms of frequencies though this is not so easy to quantify.

On the other hand, due to the very different behavior of the dispersion at low and large frequencies, the Whitham equation is a fascinating mathema-tical object since, as will be discussed later, it has three interesting different asymptotic regimes.

One should thus distinguish between the usefulness of the Whitham equation as a relevant water waves model which seems to be poor (see [15, 14] for some preliminary comparisons with experiments) and its mathematical properties that are rich, in particular as a useful and relevant toy model to provide some insights into the effect of adding a “weak”  dispersion term to a nonlinear hyperbolic equation.

We will restrict in this paper to the Whitham equation posed on the real line. Many interesting aspects of the periodic problem, in particular concerning periodic solitary waves, can be found in [22, 24, 28, 53, 55]. However our results concerning the Cauchy problem and its link with the KdV one are straightforwardly valid in the periodic case.

In order to make the link with the KdV equation (and thus with the full water wave system) we will write the Whitham equation as in [40] :

[TABLE]

where the non local operator $\mathcal{L}_{\epsilon}$ is related to the dispersion relation of the (linearized) water waves system and is defined by

[TABLE]

The (small) parameter $\epsilon$ measures the comparable effects of nonlinearity and dispersion (see below). The Whitham equation is supposed to be an approximation of the full water waves system on time scales of order $\frac{1}{\epsilon},$ see [40] and Section 2 below.

Taking the formal limit $\sqrt{\epsilon}|\xi|\to 0$ in $\mathcal{L}_{\epsilon},$ (1.1) reduces to the KdV equation

[TABLE]

We do not know of a complete rigorous justification of the Whitham equation from the water waves system. More precisely no correct scaling seems to exist allowing to connect directly and rigorously the Whitham equation and the water wave system. See however [40] pages 213-214 and Section 2 below where a comparison with the KdV equation is displayed, justifying thus the Whitham equation via the KdV approximation of weakly nonlinear long surface water waves.

On the other hand the Whitham equation can be viewed as the one-dimensional restriction of the full dispersion KP equation introduced in [40] to overcome the “bad”  behavior of the dispersion relation of the usual KP equations at low frequencies in $x$ (see also the analysis in [41]). We refer to [44] for a further study of the Cauchy problem and to [20] for the existence of localized solitary waves, “close”  to the usual KP I ones in the case of strong surface tension):

[TABLE]

with

[TABLE]

where $\beta\geq 0$ is a dimensionless coefficient measuring the surface tension effects and

[TABLE]

The link with the full water wave system is via the choice of parameters (see [40]). Denoting by $h$ a typical depth of the fluid layer, $a$ a typical amplitude of the wave, $\lambda_{x}$ and $\lambda_{y}$ typical wave lengths in $x$ and $y$ respectively, the relevant regime here is when

[TABLE]

For purely gravity waves, $\beta=0,$ (1.3) becomes

[TABLE]

with

[TABLE]

which reduces to the Whitham equation (1.1) when $u$ does not depend on $y.$

In presence of surface tension ($\beta>0$) (1.3) reduces when $u$ does not depend on $y$ to the capillary Whitham equation (see [53, 16])

[TABLE]

where

[TABLE]

Taking the formal limit $\sqrt{\epsilon}|\xi|\to 0$ in $\widetilde{\mathcal{L}}_{\epsilon},$ (1.5) reduces to the KdV equation

[TABLE]

Note also that the Whitham equation with surface tension looks, for high frequencies, like the following fractional KdV (fKdV) equation with $\alpha=1/2$

[TABLE]

We refer to [37, 42, 43] for some properties of the fKdV equations viewed as toy models to study the influence of a “weak”  dispersive perturbation on the dynamics of the Burgers equation.

There is a “Boussinesq like”, system version of (1.1) for waves propagating in both directions. As for the Whitham equation it cannot been derived directly from the water wave system by a consistent asymptotic analysis but by a rather heuristic and formal argument.

Actually, one obtains a full dispersion system when in the Boussinesq regime one keeps (formally) the original dispersion of the water waves system (see [40], [18], and [1] where interesting numerical simulations of the propagation of solitary waves are performed).111As noticed in [1], the use of nonlocal models for shallow water waves is also suggested in [61].

Setting again $\mathcal{L}_{\epsilon}=\left(\frac{\tanh\sqrt{\epsilon}|D|}{\sqrt{\epsilon}|D|}\right)^{1/2},$ we get with $D=-i\nabla$ or $-i\partial_{x}:$

[TABLE]

when $d=1$ and

[TABLE]

when $d=2.$

Taking the limit $\sqrt{\epsilon}|\xi|\to 0$ in $\mathcal{L}_{\epsilon},$ (1.9) reduces formally to

[TABLE]

while in the two-dimensional case, (1.10) reduces in the same limit to

[TABLE]

that is to the (linearly ill-posed) system one gets by expanding to first order the Dirichlet to Neumann operator with respect to $\epsilon$ in the full water wave system (see [40]). On the other hand, (1.9) and (1.10) are linearly well-posed. As will be seen below the Cauchy problem for the nonlinear system is locally well-posed under a non physical condition on the initial data though.

System (1.11) is also known in the Inverse Scattering community as the Kaup-Kupperschmidt system (see [32, 39]). It is completely integrable though linearly ill-posed since the eigenvalues of the dispersion matrix are $\pm i\xi(1-\frac{\epsilon}{3}\xi^{2})^{1/2}$ (we refer to [4] for an analysis of the ill-posedness of the nonlinear Kaup and related systems). It has explicit solitary waves (see [32]).

The Boussinesq system (1.9) can therefore be seen as a (well-posed) regularization of the Kaup-Kupperschmidt system. It is not known to be completely integrable (this is unlikely, see in particular the simulations in Section 7).

The full dispersion Boussinesq system has the following Hamiltonian structure

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

when $d=2$ and

[TABLE]

where

[TABLE]

and

[TABLE]

when $d=1.$

We will see in the next section that the full dispersion Boussinesq system has mathematical properties that make it doubtful as a relevant water wave model.

When surface tension is taken into account, one should replace the operator $\mathcal{L}_{\epsilon}^{2}$ by $\mathcal{P}_{\epsilon}=(I+\beta\epsilon|D|^{2})\left(\frac{\tanh(\sqrt{\epsilon}|D|)}{\sqrt{\epsilon}|D|}\right)$ where again the parameter $\beta>0$ measures surface tension (see [40]), yielding a more dispersive full dispersion Boussinesq system. When $\beta>\frac{1}{3},$ this full dispersion Boussinesq system yields in the limit $\sqrt{\epsilon}|\xi|\to 0$ in $\mathcal{P}_{\epsilon},$ Boussinesq systems of the class $a<0,b=c=d=0$ (see [7]) for which long time (that is on time scales of order $1/\epsilon$) well-posedness is established in [58], Theorem 4.5.

When $\beta<\frac{1}{3},$ the full dispersion Boussinesq system reduces in the formal limit $\sqrt{\epsilon}|\xi|\to 0$ in $\mathcal{P}_{\epsilon},$ to an ill-posed system, analogous to the Kaup system in dimension $1.$

Remark 1.1*.*

Another Whitham-Boussinesq system is introduced in [29].222This system seems a bit artificial contrary to (1.9), (1.10) which are the exact counter parts of the Whitham equation with respect to the original Boussinesq system It writes

[TABLE]

when $d=1$ and

[TABLE]

when $d=2.$

In the limit $\sqrt{\epsilon}|\xi|\to 0$ in $\mathcal{L}_{\epsilon},$ those systems reduce formally to

[TABLE]

in dimension one and to

[TABLE]

in the two-dimensional case. Both systems are ill-posed when $0\leq\beta<\frac{1}{3}$ while when $\beta>\frac{1}{3}$ they belong to the class of $a=b=d=0,c<0$ of classical (a,b,c,d) Boussinesq systems for which existence on time scales of order $1/\epsilon$ is established in [58], Theorems 4.6 and 4.7. The local well-posedness for (1.13) is established in [29].

According to previous theoretical results and numerical simulations, one expects at least three different regimes for the Whitham equation (1.1) (without surface tension):

1. Scattering for “small”  initial data (see the simulations in [37] and in the last Section 7 of the present paper).

2. Finite time blow-up (cusplike), see [37] for various simulations displaying the structure of a shock like blow-up (blow-up of the gradient with bounded sup-norm of the solution). The occurrence of such phenomena is rigorously proven in [30] for a class of fractional KdV equations and in [26] for the Whitham equation itself (see also the numerical simulations in [37]) but these phenomena have probably nothing to do with the breaking of real water waves. In fact, when one keeps the small parameter $\epsilon$ in the equation (that is not done in [26, 30]), the blow-up should occur on time scales much larger than $1/\epsilon,$ the time scale on which the Whitham equation is supposed to approximate the full water wave system via the KdV equation. We refer again to the simulations in Section 7.

3. A KdV, long wave regime. In fact it is shown in [21] that (1.1) possesses specific solitary waves, close to those of KdV and formally stable. In this regime one can expect a “KdV like”  behavior, namely the soliton resolution, at least on sufficiently long time scales. Those solitary waves and their perturbations are investigated numerically in Section 7 of the present paper.

One aim of the present paper is to give further evidence of the relevance of those conjectures.

The dynamics of the Whitham equation with surface tension (1.5) should be different because of the different behavior of the dispersion at high frequencies. This makes the equation more dispersive and the expected dynamics is that of $L^{2}$ critical KdV type equations. In particular the (expected) finite time blow-up should be similar to that of the $L^{2}$ critical generalized KdV equation

[TABLE]

or of the modified Benjamin-Ono equation (also $L^{2}$ critical)

[TABLE]

where $\mathcal{H}$ denotes the Hilbert transform.

The rigorous analysis of blow-up for those equations can be found respectively in [46] and [48].

As aforementioned one aim of the present paper is to give some evidence for the above conjectures via mathematical analysis and mainly by careful numerical simulations.

We will also give some hints on the qualitative behavior of the full-dispersion Boussinesq systems (1.9), (1.10).

The paper is organized as follows. In a first section we give the expected error estimates on the correct time scales, between the solutions of the Whitham and KdV equation for smooth initial data. Together with the classical results on the KdV approximation of surface water waves (see [40]) this implies a rigorous justification of the Whitham equation in the Boussinesq-KdV regime.

The next section concerns the Cauchy problem for the Whitham equation with surface tension. Contrary to (1.1) dispersive estimates can be used here to enlarge the space of resolution to the Cauchy problem.

We then comment on the local well-posedness of the Cauchy problem for the Boussinesq full dispersion systems (1.9) and (1.10) and of their capillary waves versions. We will see that the well-posedness of the system for pure gravity waves is obtained under a very restrictive, non physical condition (positivity of the wave elevation). When this condition is not satisfied the system is (Hadamard) ill posed and those facts invalidate it as a relevant model for water waves.

The presence of surface tension on the other hand prevents the appearance of Hadamard unstable modes when the initial elevation is not positive, the possible unstable modes being bounded then.

The two next sections review (and comment on) known results concerning finite time blow-up and solitary wave solutions.

Finally the last two sections display many accurate numerical simulations aiming to illustrate and to detail various properties of the Whitham equations and systems, allowing to propose convincing conjectures on their dynamics.

Section 7 is devoted to the Whitham equation. We first construct numerically the solitary waves to the Whitham equation without surface tension and simulate their perturbations for various values of the small parameter $\epsilon.$

Then we solve numerically the Cauchy problem for Gaussian initial data $\lambda\exp(-x^{2}).$ Depending on the size of $\epsilon$ and $\lambda$ a finite time blow-up may occur, at a time outside the physically relevant time scales $O(1/\epsilon)$ though.

The situation is quite different when surface tension is included. Actually one shows a finite time blow-up very similar to the one of the $L^{2}$ critical KdV equation

[TABLE]

In Section 8 we consider the one-dimensional Whitham-Boussinesq systems, with and without surface tension. We construct numerically solitary waves and study their stability. Then we simulate solutions of the Cauchy problem with initial data satisfying or not the well-posedness condition.

2. Comparison between the Whitham and KdV equations

We will compare here the solutions $v$ and $u$ of respectively the Whitham equation (1.1) and the KdV equation (1.2) with the same initial data $v_{0}=u_{0}=\phi\in H^{\infty}(\mathbb{R})$.

Here is the main result of this section.

Theorem 2.1**.**

Let $\phi\in H^{\infty}(\mathbb{R})$. Then, for all $j\in\mathbb{N}$, $j\geq 0$, there exists $M_{j}=M_{j}(\|\phi\|_{H^{j+8}})>0$ such that

[TABLE]

for all $0\leq t\lesssim\epsilon^{-1}$.

Remark 2.1*.*

The implicit constant in the notation $t\lesssim\epsilon^{-1}$ depends on $\|\phi\|_{H^{2}}$ for $j=0$ and $1$ and on $\|\phi\|_{H^{j+1}}^{-1}$ for $j\geq 2.$

Remark 2.2*.*

A similar resul holds mutatis mutandi for the periodic problem since we use energy type methods.

Recall that such a theorem was proved by Bona, Pritchard and Scott [8] for the comparison between BBM and KdV. We also refer to Albert and Bona [2] for other comparison results in the long wave regime.

2.1. The linear case

In this section, we compare the linear versions of (1.2) and (1.1), i.e. the Airy equation

[TABLE]

and the linear Whitham equation

[TABLE]

where $l(\sqrt{\epsilon}D)$ denotes the Fourier multiplier of symbol $l(\sqrt{\epsilon}\xi)$ defined by

[TABLE]

Then, we have the following result

Theorem 2.2**.**

Let $\phi\in H^{\infty}(\mathbb{R})$. Then, for all $j\in\mathbb{N}$, $j\geq 0$, there exists $N_{j}=N_{j}(\|\phi\|_{H^{j+7}})>0$ such that the solutions $u$ of (2.2) and $v$ of (2.3) associated to the same initial datum $\phi$ satisfy

[TABLE]

for all $t\in\mathbb{R}_{+}$.

Remark 2.3*.*

Note that we could also obtain bounds for the differences of $u$ and $v$ in $L^{2}$ arguing as in the proof of Theorem 2.1.

The following technical result which compares the symbol of $l(\sqrt{\epsilon}D)\partial_{x}$ and $\partial_{x}+\frac{\epsilon}{6}\partial_{x}^{3}$ will be needed below.

Lemma 2.3**.**

Assume that $\sqrt{\epsilon}|\xi|\leq 1$. Then,

[TABLE]

Proof.

The proof follows directly from the expansion

[TABLE]

∎

Proof of Theorem 2.2.

The solutions of (2.2) and (2.3) are respectively given by the unitary groups

[TABLE]

and

[TABLE]

Then, it follows from the mean value inequality and Lemma 2.3 that

[TABLE]

Then, we deduce from the Cauchy-Schwarz inequality that

[TABLE]

which concludes the proof of the theorem. ∎

Remark 2.4*.*

Note that fundamental solutions of both the linear Whitham and KdV equations have a quite different behavior.

The KdV one is given by the Airy function

[TABLE]

Set

[TABLE]

333Note that this function is different from the one used in [24] to study the properties of the solitary waves of the Whitham equation.

The fundamental solution of the Whitham equation is thus

[TABLE]

One can easily establish that contrary to the Airy function, K is an unbounded function.

Here we follow the analysis in [9] in a different setting. Actually, we write

[TABLE]

where $r$ is a continuous function exponentially decaying to zero at $\pm\infty.$ Using the elementary identity

[TABLE]

one obtains the decomposition

[TABLE]

where

[TABLE]

Since $f_{t}$ decays exponentially to zero when $|\xi|\to\infty,$ Riemann-Lebesgue lemma implies that for every $t>0,$ $I_{t}^{2}$ is a continuous function of $x$ decaying to zero at infinity.

On the other hand, following the analysis in Section 3 of [9] one can prove that $I_{t}^{1}$ decays algebraically to [math] when $|x|\to\infty$ for t fixed while for instance $I_{1}^{1}(x)\sim|x|^{-3/2}\exp(\frac{i\xi}{x})$ when $x\to 0$ for some non zero $\xi$.

Although we will not use it, we recall for the sake of completeness a dispersive estimate derived in [50], Theorem 2.5, on the free Whitham group (see also [49], Lemma 2.4). Note the difference with the classical $L^{1}-L^{\infty}$ estimates on the Airy group. We denote $S_{\epsilon}(t)=e^{itl(\sqrt{\epsilon}D)\partial_{x}}.$

Theorem 2.4**.**

There exists $C>0$ independent of $\epsilon$ such that for any $\phi\in\mathcal{S}(\mathbb{R}^{2})$ the following estimates hold :

[TABLE]

[TABLE]

[TABLE]

2.2. A priori estimates on $u$ and $v$

It is well-known that the KdV equation is well-posed in $H^{\infty}(\mathbb{R})$. Moreover, by using the complete integrability of KdV and in particular the fact that KdV possesses an infinite number of conserved quantities, one can get global bounds at the $H^{j}$-level for any $j\geq 0$. We refer for example to Saut [57], Bona and Smith [10] and Bona, Pritchard and Scott [8] and the references therein.

Proposition 2.5**.**

Let $\phi\in H^{\infty}(\mathbb{R})$. Then there exists a unique solution $u\in C([0,+\infty):H^{\infty}(\mathbb{R}))$ to (1.2) such that $u(\cdot,0)=\phi$. Moreover, the flow map data-solution $\phi\mapsto u$ is continuous from $H^{\infty}(\mathbb{R})$ into $C([0,+\infty):H^{\infty}(\mathbb{R}))$.

Furthermore, the following bounds hold true. For every $j\geq 0$, there exists $C_{j}=C_{j}(\|\phi\|_{H^{j}})$ (note that the $C_{j}$ can be chosen to be non-increasing functions of their arguments) such that

[TABLE]

Being a skew-adjoint perturbation of the Burgers equation, the Whitham equation is also trivially well-posed on $H^{\infty}(\mathbb{R})$ (and also in $H^{s}(\mathbb{R})$ for $s>\frac{3}{2}$)444Similarly to the Burgers equation, (see eg [42]) the Cauchy problem for the Whitham equation is expected to be ill-posed in $H^{\frac{3}{2}}(\mathbb{R})$ (see [27] for weaker results on related fractional KdV equations). but on a time of interval of length $1/\epsilon$.

Proposition 2.6**.**

Let $\phi\in H^{\infty}(\mathbb{R})$. Then there exist a positive time $T\sim\epsilon^{-1}$, a unique solution $v\in C([0,T]:H^{\infty}(\mathbb{R}))$ to (1.1) such that $v(\cdot,0)=\phi$. Moreover, the flow map data-solution $\phi\mapsto v$ is continuous from $H^{\infty}(\mathbb{R})$ into $C([0,T]:H^{\infty}(\mathbb{R}))$.

Furthermore, the following bounds hold true. For every $j\in\mathbb{N}$, $j\geq 2$

[TABLE]

Note that the implicit constant in (2.9) depends on $\|\phi\|_{H^{j}}$ as in (2.10).

Proof.

We only explain how to prove (2.9). Let $v\in C([0,T]:H^{\infty}(\mathbb{R}))$ be a solution of (1.1) with initial datum $v(\cdot,0)=\phi$.

Let $J^{s}$ denote the Bessel potential of order $-s$, i.e.

[TABLE]

Then, it follows from the Kato-Ponce commutator estimates [33] and integrations by parts that

[TABLE]

for any $s>0$. If $s>\frac{3}{2}$, we deduce from the Sobolev embedding $H^{s-1}(\mathbb{R})\hookrightarrow L^{\infty}(\mathbb{R})$ that

[TABLE]

Hence, we deduce from a classical ODE argument that

[TABLE]

This finishes the proof of Proposition 2.6. ∎

Finally, we also need a bound on $\pi_{\epsilon}(D)u$ where $u$ is the solution of (1.2) and $\pi_{\epsilon}(D)$ is the Fourier multiplier of symbol $\pi_{\epsilon}(\xi)$ defined by

[TABLE]

Proposition 2.7**.**

Let $\phi\in H^{\infty}(\mathbb{R})$ and let $u$ be the solution of (1.2) evolving from $\phi$ obtained in Proposition 2.5. Then, for all $j\in\mathbb{N}$, $j\geq 0$, there exists $A_{j}=A_{j}(\|\phi\|_{H^{j+8}})>0$ such that

[TABLE]

for all $0\leq t\leq\epsilon^{-1}$.

Proof.

We apply the operator $\partial_{x}^{j}\pi_{\epsilon}(D)$ to (1.2), multiply the equation by $\partial_{x}^{j}\pi_{\epsilon}(D)u$ integrate in the space variable over $\mathbb{R}$ and integrate by parts and use the Cauchy-Schwarz inequality to deduce that

[TABLE]

Now, we introduce a cut-off function $\eta\in C_{0}^{\infty}(\mathbb{R})$ such that $\text{supp}\,(\eta)\subset[-1,1]$ and $\eta=1$ over $[-1/2,1/2]$ and we define the Fourier multiplier $P_{<1/\sqrt{\epsilon}}$ of symbol $\eta(\sqrt{\epsilon}\cdot)$ and $P_{\geq 1/\sqrt{\epsilon}}=1-P_{<1/\sqrt{\epsilon}}$. Then, the triangle inequality yields

[TABLE]

and we need to estimate both terms on the right-hand side of (2.14).

To control the first one, we use Plancherel’s identity, Lemma 2.3 and the fact that $H^{s}(\mathbb{R})$ is a Banach algebra for $s>\frac{1}{2}$ to deduce that

[TABLE]

which together with Proposition 2.5 gives that

[TABLE]

To control the second one, we use again Plancherel’s identity to obtain that

[TABLE]

Now, observe that $1\lesssim\epsilon^{4}|\xi|^{8}$ on the support of the integral, so that

[TABLE]

since $H^{8}(\mathbb{R})$ is a Banach algebra and where we used Proposition 2.5 on the last inequality.

Then, we deduce gathering (2.13)–(2.16) that

[TABLE]

Therefore, we deduce integrating between [math] and $t$ that

[TABLE]

Finally, arguing as above, we get from Lemma 2.3 and Plancherel’s identity that

[TABLE]

Then, we deduce combining (2.18) and (2.19) that

[TABLE]

as soon as $0\leq t\leq t/\epsilon$. ∎

2.3. Proof of Theorem 2.1

Let $z=u-v$ and $j\in\mathbb{N}$, $j\geq 0$. It is deduced from equations (1.1) and (1.2) that $z$ solves the initial value problem

[TABLE]

where $\pi_{\epsilon}(D)$ is the Fourier multiplier of symbol $\pi_{\epsilon}(\xi)$ defined in (2.11).

Note that the equation of $z$ is well defined for $0<t\lesssim 1/\epsilon$.

Differentiate $j$ times, multiply the equation (2.21) by $\partial_{x}^{j}z$, integrate in the space variable over $\mathbb{R}$ and integrate by parts to deduce that

[TABLE]

Thus, the Cauchy-Schwarz inequality, the Leibniz rule and integration by parts yield

[TABLE]

in the cases $j=0$ and $j=1$ and

[TABLE]

in the cases where $j\geq 2$.

Now, on the one hand, we get by using the Sobolev embedding and Propositions 2.5 and 2.6 that

[TABLE]

for all $0\leq t\lesssim\epsilon^{-1}$ (where the implicit constant depends on the $\|\cdot\|_{H^{\kappa}}$ norm of $\phi$ as explained in the proof of proposition 2.6). Here, we used the notation $\kappa=2$ for $j=0$ or $1$ and $\kappa=j+1$ for $j\geq 2$. On the other other hand, by applying Proposition 2.7, it follows that

[TABLE]

for all $0\leq t\leq\epsilon^{-1}$.

Therefore, we deduce gathering (2.23), (2.25), (2.26) and using Gronwall’s inequality that

[TABLE]

whenever $0\leq t\lesssim\epsilon^{-1}$, which concludes the proof of Theorem 2.1.

Remark 2.5*.*

A similar comparison result can be established between the Whitham equation with surface tension (1.5) and the KdV equation (1.7) by using the expansion

[TABLE]

3. The Cauchy problem for the capillary Whitham equation

As was already noticed the Cauchy problem for the Whitham equations (1.1) and (1.5) is trivially well-posed in $H^{s}(\mathbb{R}),\;s>\frac{3}{2}$ on time scales of order $\epsilon^{-1}.$

This result can be improved (by enlarging the space of resolution) for the capillary Whitham equation (1.5) by using its dispersive properties. Actually one gets, see [43] where general fractional KdV (fKdV) equations are considered:

Theorem 3.1**.**

Assume that $s>\frac{21}{16}$. Then, for every $u_{0}\in H^{s}(\mathbb{R})$, there exist a positive time $T_{\epsilon}=T_{\epsilon}(\|u_{0}\|_{H^{s}})=0(1/\sqrt{\epsilon})$ (which can be chosen as a non-increasing function of its argument), and a unique solution $u$ to (1.5) satisfying $u(\cdot,0)=u_{0}$ such that

[TABLE]

Moreover, for any $0<T^{\prime}<T_{\epsilon}$, there exists a neighborhood $\mathcal{U}$ of $u_{0}$ in $H^{s}(\mathbb{R})$ such that the flow map data-solution

[TABLE]

is continuous.

Remark 3.1*.*

1. Since the value $\alpha=\frac{1}{2}$ is $L^{2}$ critical for the fKdV equation we conjecture that the Cauchy problem for (1.5) is globally well posed for initial data in the energy space $H^{\frac{1}{4}}(\mathbb{R})$ having a sufficiently small $L^{2}$ norm.

2. This result was recently improved in [51] where the local well-posedness is obtained for $s>\frac{9}{8}.$

3. For “large enough”  initial data, one expects a finite time blow-up silmilar to that of the $L^{2}$-critical generalized KdV equation or to the cubic-Benjamin-Ono equation proven respectively in [46] and [48], as displayed in the simulations of Section 7.

Remark 3.2*.*

Using the same energy methods one can prove a result similar to that of Theorem 2.1 between solutions of (1.5) and (1.7).

4. The Cauchy problem for the Boussinesq-Whitham systems

In both cases (without or with surface tension) the Boussinesq-Whitham systems are linearly well-posed since the linearized systems write, say in dimension one

[TABLE]

where the Fourier transform of the matrix operator A has real eigenvalues

[TABLE]

in absence of surface tension and

[TABLE]

in presence of surface tension.

The theory of the Cauchy problem for the nonlinear system without surface tension is relatively straightforward and we focus on the case of space dimension one. Since the operator $\mathcal{L}_{\epsilon}^{2}\partial_{x}$ has order zero, (1.9) is an order zero perturbation of a first order system, more precisely it writes

[TABLE]

where $U=\begin{pmatrix}\eta\\ u\end{pmatrix}$, $\mathcal{A}_{\epsilon}(U)=\begin{pmatrix}\epsilon u&\epsilon\eta\\ 1&\epsilon u\end{pmatrix},$ $\mathcal{H}_{\epsilon}U=\begin{pmatrix}\mathcal{L}_{\epsilon}^{2}\partial_{x}u\\ 0\end{pmatrix}.$

The matrix $\mathcal{A}_{\epsilon}(U)$ has eigenvalues $\lambda$ satisfying

[TABLE]

so that the system

[TABLE]

is hyperbolic in the regions where

[TABLE]

and Hadamard ill-posed if this condition is violated.

Actually, under the condition

[TABLE]

the system can be symmetrized via the positive definite symmetrizer $S(U)=\begin{pmatrix}1&0\\ 0&\epsilon\eta\end{pmatrix}.$

The same process can be applied to (1.9) yielding a symmetric hyperbolic system perturbed by the order zero operator $\mathcal{H}_{\epsilon}$ namely

[TABLE]

Setting $\zeta=\eta+N_{0}$ where $N_{0}>0$ is a fixed constraint, the standard theory of symmetrizable hyperbolic systems (see eg [45]) imply the local well-posedness of the Cauchy problem for $(\zeta,u)\in H^{s}(\mathbb{R})\times H^{s}(\mathbb{R}),\;s>\frac{3}{2}$ for initial data $\zeta_{0}$ such that $\eta_{0}$ is sufficiently small. Similar arguments in the two-dimensional case yield well-posedness in $H^{s}(\mathbb{R}^{2})\times H^{s}(\mathbb{R}^{2}),\;s>2.$

Note that the condition (4.2) implies that the wave is always of elevation which seems to be the case of solitary wave solutions (see Section 7). On the other hand, condition (4.3) implies that the wave cannot tend to zero at infinity and thus the perturbations of the solitary wave solutions are excluded from the range of well-posed initial data. Those facts are not physically realistic, invaliding the Boussinesq-Whitham system as a relevant water waves model.

Remark 4.1*.*

We do not know of a local well-posedness result for (1.9) or (1.10) under the assumption (4.2).

To illustrate the (Hadamard) instabilities occurring when condition (4.3) is violated, let consider the linearization around $(\eta,u)=(-c,0)$ where $c$ is a positive constant. The linearized system has eigenvalues $\lambda_{\pm}$ where

[TABLE]

All modes are unstable when $\epsilon>1/c$ while when $\epsilon<1/c$ all modes corresponding to $|\xi|\geq\frac{x_{c,\epsilon}}{\sqrt{\epsilon}}$ are unstable where $x_{c,\epsilon}$ is the unique positive solution of $\frac{\tanh x}{x}=c\epsilon.$ Note that $x_{c,\epsilon}\to\infty$ as $\epsilon\to 0.$ Thus the set of stable modes get larger and larger when for a fixed $c,$ the small parameter $\epsilon$ tends to $0.$

When the nonlocal term is removed, all modes are unstable, for any $c>0.$ The effect of dispersion is thus here to create a range of linearly stable modes when $c<1,$ getting larger and larger when $\epsilon\to 0.$

Things are different for Boussinesq-Whitham systems in presence of surface tension since dispersive effects play a significant role here.

In fact, surface tension prevents the appearance of (Hadamard) unstable modes for the linearized system at $(\eta,u)=(-c,0)$.

Actually the linearized system eigenvalues are now

[TABLE]

Whatever the values of $c$ and $\epsilon$ the possible unstable modes are bounded and there are no more Hadamard instabilities.

The effect of surface tension is thus to suppress the Hadamard instabilities.

A convenient way to get an idea of the nature of the full system is to derive an equivalent system by diagonalizing the linear part of the system. More precisely, we define

[TABLE]

the Fourier transform of the dispersion matrix with eigenvalues

[TABLE]

Setting

[TABLE]

and

[TABLE]

the linear part of the system is diagonalized as

[TABLE]

where

[TABLE]

that is two dispersive equations of order $\frac{1}{2}.$

In the $W$ variable the complete system writes now

[TABLE]

where

[TABLE]

While the linear part is of order $\frac{1}{2},$ the nonlinear terms are formally of order $1.$ However the structure of the nonlocal nonlinear terms might prevent the implementation of the methods used in [43] for fractional type KdV equations to obtain the local well-posedness of the Cauchy problem on time scales of order $1/\sqrt{\epsilon}.$

On the other hand one might think of proving the local well-posedness by using an elementary energy method, reminiscent of the one used in [58] for the $(-1,0,0,0)$ Boussinesq system (see [58] Theorem 3.1) or a symmetrization technique as in [58] Theorem 4.5. However those methods do not seem to work since a control of the $H^{s}\times H^{s+1/2},\quad s>\frac{1}{2}$ norm on $(\eta,u)$ is not enough to control the nonlinear terms (in case of the $(-1,0,0,0)$ Boussinesq system one gets a control on the $H^{s}\times H^{s+1}$ norm).

5. Finite time blow-up for the Whitham equation

It has been established in [26] that solutions of the Whitham equation (1.1) may blow-up in finite time, the blow-up being "shock-like", that is blow-up of the gradient with bounded $L^{\infty}$ norm. The numerical simulations in [37] display the cusp nature of the singularity depending on the sign of the initial data. On the other hand those works consider the Whitham equation (1.1) with $\epsilon=1$, thus they do not address the long-wave, KdV limit. In particular the crucial dependence of the blow-up time with respect to $\epsilon$ is not obtained. Note that this blow-up time should be larger that $1/\epsilon,$ the time scale on which the KdV and Whitham equations are close and asymptotic models for the propagation of long, weakly nonlinear surface water waves for which no singularities are expected.

Note also that a quite different finite time blow-up is expected for the Whitham equation with surface tension (1.5) that cannot be anymore viewed as a "weak" dispersive perturbation of the Burgers equation. In fact, as was already noticed, it is reminiscent for large frequencies of the $L^{2}$ critical fractional KdV equation

[TABLE]

for which a finite time blow-up à la Martel-Merle [46, 48] is expected (but not yet proven), see the numerical simulations in [35, 36, 37].

6. Solitary waves

6.1. Whitham without surface tension

We will focus here on non periodic solitary wave solutions of the Whitham equation (1.1), that is solutions of the form $u(x-ct).$

We refer to [12, 22, 23, 28, 55] for interesting (theoretical and numerical) studies on periodic traveling waves. In particular the existence of a global bifurcation branch of $2\pi-$ periodic smooth traveling wave solutions is established in [23].

The existence of solitary wave solutions to the Whitham equation (1.1) decaying to zero to infinity and close to the KdV soliton has been proven in [21] by variational methods. Symmetry and decay properties of Whitham solitary waves are established in [13]. The analysis in [21] is made on (1.1) with $\epsilon=1$ under the scaling

[TABLE]

where $2\alpha-\beta=1$ so that $\frac{1}{2}\int_{\mathbb{R}}u^{2}=\epsilon$.

An interesting issue is that of the transverse stability of the Whitham solitary wave in the framework of the full dispersion KP equation (1.3). Considering the similar problem for the usual KP I/II equations (see [56] and the references therein) one can conjecture that the Whitham solitary wave is transversally stable when the surface tension parameter $\beta$ is less than $\frac{1}{3}$ and unstable otherwise. We intend to go back to this issue in a subsequent paper.

We do not know of any rigorous result on the existence of solitary wave solutions to the full dispersion Boussinesq system.

The following computations give some evidence to the existence of solitary wave solutions "close" to the KdV soliton.

Recall that the system reads

[TABLE]

where $\hat{T}_{\epsilon}=(\tanh(k\sqrt{\epsilon})/(k\sqrt{\epsilon}))^{1/2}$. In the formal limit $\epsilon\to 0$, the system reduces to the wave equation.

We are here interested in localized travelling wave solutions to the system (6.1). Putting $u(x,t)=U(x-ct)$ and $\eta(x,t)=N(x-ct)$, we find after integration

[TABLE]

Eliminating $N$ from the first equation of (6.2) via the second, we find

[TABLE]

Writing $c=1+\alpha\epsilon$ with $\alpha>0$, we get, by performing a formal expansion in $\epsilon$ and neglecting terms of order $\epsilon$ or higher in (6.3)

[TABLE]

which gives after integration

[TABLE]

which has the KdV soliton

[TABLE]

as a solution. Note that for a solution $(N,U)$ with a given velocity $c$ of the above equations, $(-N,U)$ provides a solution for the same system with velocity $-c$.

Similarly, a solitary wave $u(x,t)=U(x-ct)$ of the Whitham equation (1.1) satisfies the equation

[TABLE]

Remark 6.1*.*

Following the method used in [13] one can prove that solutions of (6.7) that tend to [math] as $|x|\to\infty$ decay to [math] exponentially in the sense that for some $\nu>0,$ $e^{\nu|\cdot|}\phi\in L^{1}(\mathbb{R})\cap L^{\infty}(\mathbb{R}).$

Writing $c=1+\delta\epsilon$ where $\delta>0$ is a constant independent of $\epsilon$, we get, neglecting terms of order $\epsilon$ or higher in the formal expansion in $\epsilon$,

[TABLE]

which gives after integration

[TABLE]

Thus we find again the KdV soliton

[TABLE]

Remark 6.2*.*

The above formal considerations suggest that one could obtain the existence of solitary wave solutions to both the Whitham equation and the Whitham system by perturbations arguments starting from the KdV soliton. In particular in the case of the Whitham equation this would give an alternative proof to the one in [21].

6.2. Whitham with surface tension

The situation is quite different for the Whitham equation with surface tension (1.5).

A solitary wave solution $u(x,t)=U(x-ct)$ of (1.5) should satisfy the equation

[TABLE]

where $\tilde{T}_{\epsilon}$ is defined by (1.6).

The existence of localized non trivial solutions of (6.11) results from Theorem 2.1 in [6]. Since (1.5) (as the fKdV equation with $\alpha=\frac{1}{2}$) is $L^{2}$ critical one expects the instability of those solitary waves by blow-up, similarly to the generalized KdV equation with $p=4$, see [46] or to the modified Benjamin-Ono equation ([48]).

On the other hand, a solitary wave solution $(\eta(x-ct),u(x-ct))$ of the full-dispersion Boussinesq system should satisfy

[TABLE]

where $\mathcal{P}_{\epsilon}=(I+\beta\epsilon|D|^{2})\left(\frac{\tanh(\sqrt{\epsilon}|D|)}{\sqrt{\epsilon}|D|}\right).$

Eliminating $\eta$ from the second equation in (6.11) yields the single equation

[TABLE]

One cannot apply directly the result of [6] to prove the existence of non trivial solutions to (6.13) since the nonlinear term is not homogeneous.We plan to come back to this issue in a next paper.

7. Numerical simulations for the Whitham equation

The simulations in this section will illustrate several aspects of the dynamics of the Whitham equation alluded to in the previous sections and will also give evidence for some facts not yet rigorously proven.

7.1. Numerical construction of solitary waves for the Whitham

equation

In this subsection we construct numerically solitary wave solutions to the Whitham equations, i.e., localized solutions of (6.7). To this end we use the same technique as in [37]: since the wanted solutions are expected to be rapidly decreasing, it is convenient to formulate the problem as an essentially periodic one with a sufficiently large period such that the solutions and the numerically interesting derivatives vanish with numerical accuracy at the domain boundaries $\pm L\pi$ (we typically choose $L=5$ in the following). The Fourier transforms are approximated in standard way via discrete Fourier transforms conveniently implemented via the Fast Fourier transform (FFT).

With this approach equation (6.7) is approximated via a finite dimensional (we use $N$ Fourier modes) nonlinear equation system formally written as

[TABLE]

where $\hat{U}$ denotes the discrete Fourier transform of $U$. This nonlinear system is solved with a Newton-Krylov iteration. This means the action of the inverse of the Jacobian of $\mathcal{F}$ on $\mathcal{F}$ in the standard Newton iteration is determined iteratively via GMRES [54]. As the initial iterate we choose the KdV soliton written as $U=3\delta\mbox{sech}^{2}\left(\sqrt{\frac{3\delta}{2}}(x-ct)\right)$. Note that the reality of $U$ has to be enforced during the iteration.

As an example we study the case $\epsilon=0.01$ and vary the constant $\delta$. This is mainly done in order to have for small $\delta$ an initial iterate close to the Whitham soliton. For larger $\delta$, we use as an initial iterate the solution for a slightly smaller value of $\delta$ constructed before. Since we vary $\delta$, the choice of $\epsilon$ is, however, no restriction for the general case since the only important quantity is the velocity $c$. It is to be expected that there will be solitons for $c\sim 1$ because of the analogy to KdV, but that for larger values of $c$ the dispersion of the Whitham equation becomes in contrast to the KdV equation too weak to compensate the nonlinearity. Thus we expect that in contrast to KdV there might be an upper limit of $c$ for solitons to the Whitham equation. In fact we obtain the situation shown in Fig. 1. For $c=1.01,1.02,1.05,1.1,1.2$ we use $N=2^{14}$ Fourier modes, for $c=1.22$ $N=2^{16}$.

It can be seen that maxima of the solution grow as expected, but that the solitons become more peaked and more compressed compared to the corresponding KdV case. In fact we did not succeed to construct solitons for $c$ much greater than $1.22$. The failure of an iteration to converge obviously does not imply that there will be no longer a soliton, it just means that the numerical approach no longer can be used. However the fact that the iteration did not converge even for larger resolutions, for initial iterates being the numerical solution for a slightly smaller $\delta$ and for an iteration with relaxation is a strong indication that there might not be a Whitham solitary waves for much larger velocities.

There are additional differences between Whitham and KdV solitons. The mass (the $L^{2}$ norm) of the KdV soliton (6.10) is given by $M=4\sqrt{6}\delta^{3/2}$, the mass of the Whitham solitons can be seen on the left of Fig. 2. As expected the mass is identical for small $\delta$ to the one of the KdV soliton. But for larger speeds $c$, the mass grows more slowly. The energy of the KdV soliton (6.10) is proportional to $\delta^{5/2}$, the energy of the Whitham soliton can be seen on the right of Fig. 2. The energy curve appears to flatten near the reachable maximal velocities.

7.2. Numerical study of perturbed Whitham solitons

In this subsection, we study the time evolution of perturbations of the Whitham solitons constructed above. To this end we write the Whitham equation for the Fourier transform in a commoving frame, i.e.,

[TABLE]

and approximate the Fourier transform via a discrete Fourier transform as before. For the time integration we use as in [37] an implicit Runge-Kutta method of fourth order with a fixed point iteration, see [37] for details. The accuracy of the solution is controlled via the Fourier coefficients which should decrease exponentially to machine precision (we work here in double precision which allows a maximal accuracy of $10^{-16}$) for smooth functions and via the conservation of the energy

[TABLE]

Due to unavoidable numerical errors, this energy does numerically evolve with time, but thus provides an estimate of the numerical error. As shown in [34], conserved quantities typically overestimate the $L^{\infty}$ error (the maximum of the difference between numerical and exact solution) by 1-2 orders of magnitude.

As a first test of the code, we take the numerically constructed soliton with $\epsilon=0.01$ and $c=1.2$ of the previous subsection as initial data. In the used commoving frame, the solitons should correspond to a stationary solution. We use $N_{t}=10^{4}$ time steps for $t\in[0,10]$ and find that the difference between the evolved solution and the initial data is of the order of $10^{-13}$, the order of the error with which the equation (6.7) had been solved in the previous section. The numerically computed relative energy $1-E(t)/E(0)$ is of the order of $10^{-15.5}$.

We first consider as initial data the soliton with $\epsilon=0.01$ $c=1.2$ plus a Gaussian perturbation which is a perturbation of the order of $1\%$ of the maximum of the soliton. We use $N_{t}=2*10^{4}$ time steps for $t\in[0,20]$. The Fourier coefficients in this case decrease to the order of $10^{-14}$ during the whole computation, and the relative energy is conserved to the order better than $10^{-14}$. In this case there appears a slightly faster soliton travelling to the right in the commoving frame, as can be seen in Fig. 3. But in addition to radiation travelling to the left, a slower soliton travelling to the left in the frame commoving with the unperturbed soliton seems to emanate from the initial data. This is even more visible in the plot on the right of Fig. 3 where a close up of the smaller soliton is shown for $t=20$. The soliton appears to be stable, the perturbation leads to a slightly larger soliton and possibly a smaller soliton plus radiation.

The situation changes if the same soliton as in Fig. 3 is considered, but this time with a perturbation $3\exp(-x^{2})$ corresponding to roughly $4\%$ of the maximum of the soliton. We use $N_{t}=10^{4}$ time steps for $t\in[0,8]$. In this case the solution appears to develop a cusp. We show in Fig. 4 the solution and the modulus of its Fourier coefficients at $t=7.52$. It can be seen that the solution clearly runs out of resolution in Fourier space, already at $t\sim 7$, the Fourier coefficients decrease only to $10^{-2}$.

In Fig. 5 it can be seen that the $L^{\infty}$ norm of the solution grows clearly beyond the initial data, but it is not clear whether it diverges. The same behavior can be seen for the $L^{2}$ norm of $u_{x}$ on the right of the same figure. A fit of the Fourier coefficients as in [59, 37] reveals that a singularity appears to approach the real axis in the complex $x$-plane. The idea of this approach is that an essential singularity of the form $(x-x_{s})^{\mu}$, $\mu\notin\mathbb{Z}$, in the complex plane for a function $u$ leads to an asymptotic behavior of the Fourier transform

[TABLE]

where $\delta={\rm Im}x_{s}$. For $t=7.547$, this singularity seems to hit the real axis ($\delta\to 0$) indicating the formation of a cusp (the fitted coefficient $\mu$ of the essential singularity is positive as in [37]). Thus there seems to be a hyperbolic blow-up in this case. The soliton for $c=1.2$ appears therefore to be unstable against blow-up in the form of a cusp for sufficient size of the perturbation.

Since the solitons of the Whitham equation are rapidly decreasing, it is possible to study soliton interactions. In Fig. 6, we show the solitons with $c=1.05$ and $c=1.1$, the latter being centered at $x=-4$ by simply multiplying its Fourier coefficients by $e^{4ik}$. The sum of the solutions gives two soliton initial data since they both vanish with numerical precision where the other takes values above the numerical error. On the left of Fig. 6 we show the solution to the Whitham equation for $t\in[0,100]$. The modulus of the Fourier coefficients of the solution decreases to $10^{-12}$ during the whole computation. It can be seen that the soliton interaction resembles the KdV two-soliton: the solitons have almost the same shape after the collision, there is just a phase shift (note that the computation is performed in a frame commoving with the smaller soliton at $c=1.05$). The close-up of the solution at the final time on the right of Fig. 6 reveals, however, that this is not an exact two-soliton since there is dispersive radiation propagating to the left. This indicates that the Whitham equation is as expected not integrable.

7.3. Numerical study of Gaussian initial data for different values

of $\epsilon$

In this section we study numerically the time evolution of Gaussian initial data $u(x,0)=10\exp(-x^{2})$ for different values of the parameter $\epsilon$. Since we are interested in studying the solutions on time scales of order $1/\epsilon$, we introduce the rescaled time $\tau=t\epsilon$ in which the Whitham equation (1.1) reads

[TABLE]

where we have used a frame commoving with velocity 1. This equation is solved with the same numerical approach as in the previous subsection.

For $\epsilon=1$ we use $N=2^{14}$ Fourier modes and $N_{t}=10^{4}$ time steps for $\tau\in[0,0.2]$. The solution appears to develop a cusp for $\tau>0.1175$ as can be seen in Fig. 7 on the left. The Fourier coefficients on the right of the same figure indicate a loss of resolution. A fitting of the Fourier coefficients to (7.4) indicates that for $\tau_{c}=0.1175$, the parameter $\delta$ indicates the distance between a singularity in the complex $x$-plane vanishes, and the corresponding parameter $\mu\sim 0.3635$ indicating the formation of a cusp proportional to $(x-x_{s})^{1/3}$.

Note that the critical time $\tau_{c}$ above does not indicate the exact blow-up time, this just indicates that the singularity in the complex plane leading to the asymptotic behavior (7.4) of the Fourier coefficients is too close to the axis to be numerically distinguished from 0. If one considers the same initial data as in Fig. 7 for smaller values of $\epsilon$, one obtains the critical times and exponents shown in the table in Fig. 1. It can be seen that the critical times $\tau_{c}$ are always greater than the $\tau_{c}$ for $\epsilon=1$. In fact the critical times appear to grow with decreasing $\epsilon$. This implies, however, that the solutions do not have a blow-up on time scales of order $1/\epsilon$ (recall that $\tau=t\epsilon$).

For values of $\epsilon$ smaller than $0.04$, there does not appear to be a blow-up at all. In Fig. 8 we show the solution to the same initial data as in Fig. 7 for $\epsilon=0.01$. We use $N=2^{12}$ Fourier modes for $x\in L[-\pi,\pi]$ with $L=10$ and $N_{t}=10^{4}$ time steps for $\tau\in[0,2]$. It can be seen that two stable solitons appear which gives support to the soliton resolution conjecture in the context of the Whitham equation: stable solitons seem to appear in the long time behavior of the solution. In addition there is the usual dispersive radiation also known from KdV solutions. The $L^{\infty}$ norm of the solution on the right of Fig. 8 also seems to correspond to a soliton with speed greater than 1 (we are in a commoving frame with $c=1$). Neither the Fourier coefficients nor the $L^{2}$ norm of $u_{x}$ indicate the formation of a singularity. The latter again indicates the appearance of a soliton.

The solution in Fig.8 for small $\epsilon$ becomes obviously closer to the solution to the KdV equation $u_{\tau}+uu_{x}+u_{xxx}/6$, which can be obtained from the Whitham equation in the formal limit $\epsilon\to 0$ for the same initial data. This KdV solution is shown in Fig.9

7.4. Numerical study of the Whitham equation with surface

tension

In this subsection we will study the Whitham equation with surface tension (1.5) for similar situations as above for the Whitham equation without surface tension: solitons, their stability and general initial data in the Schwartz class. The goal is mainly to highlight differences with respect to the case without surface tension.

Arnesen [6] showed that solitons exist in the case of non-vanishing surface tension for all values of $\epsilon$. For small $\epsilon$, a similar calculation as for the Whitham equation in (6.10) leads to the result that the Whitham soliton should be in some sense close to the KdV soliton given by

[TABLE]

This means that for $\beta<1/3$, the solitons should have positive amplitude, whereas they are depression waves (negative amplitude) for $\beta>1/3$. Numerically we construct the solutions as before. Since the behavior for the high wave numbers of the operator $\hat{\mathcal{L}}_{\epsilon}$ (1.6) is different (linear in $|k|$) than for the case $\beta=0$ (constant in $|k|$), it is useful to divide by $\hat{\mathcal{L}}_{\epsilon}$ in the equation for the solitary wave in Fourier space. This is straight forward because the operator is diagonal in Fourier space. For small $\epsilon$ and $\delta$ the solitons are close to the KdV soliton. All figures are created with $N=2^{12}$ Fourier modes.

We first address the case $\beta=0.1$ and $\delta=2$ which can be seen on the left of Fig. 10. The solitary wave for the Whitham equation with surface tension is slightly more peaked than the KdV soliton in red in this case. If we propagate the exact solitary wave in a commoving frame by the time evolution code (as in [36] it is useful to apply a simplified Newton iteration instead of a fixed point iteration), the difference between numerical solution and initial data for $N_{t}=10^{4}$ time steps for $\tau<1$ is of the order of the numerical error for the solitary wave ($10^{-12}$).

If we perturb this solution slightly by considering the initial data $u(x,0)=0.99U(x,c=1.02)$, i.e., initial data close to the solitary wave with a slightly smaller mass, one can see on the right of Fig. 10 that the soliton is dispersed to infinity. In contrast to the case without surface tension, there are no visible oscillations, the initial pulse just gets broader and decreases in height.

If we consider for the same values of $\epsilon$ and $\beta$ as in Fig. 10 the initial data $u(x,0)=1.01U(x,c=1.02)$, i.e., initial data in the vicinity of the solitary wave with slightly larger mass, one can see in Fig. 11 that the solution has an $L^{\infty}$ blow-up in finite time which is clear from the right figure. In fact the divergence of the $L^{\infty}$ norm after $\tau=5.96$, for which the solution is shown on the left, is so rapid that the fitted distance (7.4) changes from positive to negative values too rapidly to obtain a sensible result for $\mu$ (the Fourier coefficients deteriorate too rapidly). Thus one finds the same behavior for almost solitary initial data as for the $L^{2}$ critical gKdV equation [47]: exact solitary initial data are not affected, initial data with smaller mass will be dispersed to infinity, and initial data with larger mass lead to an $L^{\infty}$ blow-up.

For $\beta=1$ and $\delta=-20$ the soliton has negative amplitude and can be seen on the left of Fig. 12. The solitary wave for the Whitham equation with surface tension is again slightly more peaked than the KdV soliton in red. If we perturb this solution by considering the initial data $u(x,0)=0.99U(x,c=0.8)$, one can see on the right of Fig. 12 that the soliton is once more dispersed to infinity, this time to the right, but again without visible oscillations.

If we consider for the same values of $\epsilon$ and $\beta$ as in Fig. 12 the initial data $u(x,0)=1.01U(x,c=0.8)$, one can see in Fig. 13 that the solution has again an $L^{\infty}$ blow-up in finite time which is clear from the right figure. The divergence of the $L^{\infty}$ norm after $\tau=4.75$, for which the solution is shown on the left, is so rapid that the iteration stops converging. Once more one finds the same behavior for almost solitary initial data as for the $L^{2}$ critical gKdV equation.

To address the case of general initial data in the Schwartz class, we consider the same initial data as in the case $\beta=0$ in Fig. 7. The numerical results indicate that there is a hyperbolic blow-up in this case. In Fig. 14 we show in contrast the corresponding situation for $\beta=1$. We use $N=2^{14}$ Fourier modes on $x\in[-2\pi,2\pi]$ and $N_{t}=50000$ time steps for $\tau\in[0,0.2]$. For small times, the dynamics of Burgers’ equation is again dominant leading to a steepening of the right front towards the formation of a shock. Close to a potential gradient catastrophe, the stronger dispersion (compared to the Whitham equation with $\beta=0$) takes over and appears to generate a dispersive shock as known from the KdV equation. However, as becomes clear from the right figure in Fig. 14, instead of a dispersive shock wave an $L^{\infty}$ blow-up as for generalized KdV equations is observed, see e.g. [46, 36] and references therein.

A fitting of the Fourier coefficients in Fig. 15 according to (7.4) indicates indeed that a singularity in the complex plane approaches the real axis for $\tau\approx 0.1648$. The critical exponent is found to be $\mu=-0.5518$ which confirms an $L^{\infty}$ blow-up. This is in accordance with the $L^{\infty}$ norm of the solution shown in Fig. 15 on the right.

8. Numerical simulations for the Boussinesq system

In this section we numerically construct solitary waves for the Boussinesq system and study their stability as well as solutions to more general initial data.

8.1. Numerical construction of solitary waves for the

Boussinesq system

As in the case of solitary waves for the Whitham equation, we construct the solitary waves by solving equation (6.3) with FFT techniques in Fourier space with a Newton-GMRES iteration. We choose again $\epsilon=0.01$ and use the KdV soliton (6.6) as initial iterate for small $\alpha$. For larger values of $\alpha$ (which implies larger values of the speed $c=1+\alpha\epsilon$), we use the numerical solution to (6.3) for a slightly smaller $\alpha$ as an initial iterate. The iterations are carried out on the interval $x\in[-5\pi,5\pi]$ with $N=2^{14}$ Fourier modes. For $c=1.16$ we use $2^{16}$ Fourier modes, but even increasing $c$ in small increments, we do not observe convergence of the iteration. Again this does not prove the non-existence of solitary waves to the system (6.1) at higher speeds, but is an indication that there might be an upper limit to the speed of the travelling waves for this equation.

In Fig. 16, we show on the left the function $U$ for different values of $c$. It can be seen that with increasing speed, the solitons become again more peaked and localized. Note that all shown solutions are numerically well resolved in the sense that the modulus of the Fourier coefficients decreases to machine precision for large wave numbers. On the right of Fig. 16 we show the corresponding functions $N$ computed for given $U$ via (6.2). Due to the term proportional to $\epsilon$ in (6.2), the functions $N$ are less peaked than the corresponding function $U$ for large $c$.

8.2. Numerical study of perturbed solitary waves

In this section we study numerically solutions to the Cauchy problem for the Boussinesq system (6.1) for several examples. We use the same approach as for the Whitham equation, a Fourier spectral method in $x$ and an implicit Runge-Kutta method of fourth order with a fix point iteration in $t$. Since the system can be seen as a dispersive regularization of the ill-posed Kaup system and since it is ill posed for negative $\eta$, Krasny filtering [38] has to be used in order to stabilize the code as in the case of for instance the focusing nonlinear Schrödinger equation in the semiclassical limit, see e.g. [34, 17]; this means the Fourier coefficients with a modulus smaller than some threshold, typically $10^{-12}$, are put equal to zero. The accuracy of the solution is controlled via the decrease of the modulus of the Fourier coefficients for large wave numbers and the conservation of the numerically computed energy.

If we choose the numerically constructed soliton solution as initial data in a commoving frame, the solution is for $c=1.05$ numerically evolved for $t\in[0,2]$ with an accuracy of $10^{-14}$ ($L^{\infty}$ norm of the difference between numerical solution and initial data). This shows the accuracy both of the time evolution code and the code for the solitons. Note that a rescaling of the time as in (7.5) is not straight forward for the Boussinesq system. Therefore we will always consider the time $t$ in this section.

Using initial data of the form $u(x,0)=U(x)+\exp(-x^{2})$, $\eta(x,0)=N(x)$, one gets the solution shown in Fig. 17. It can be seen that a slightly larger soliton emerges, and that the remaining energy is radiated to the left.

As in the case of the Whitham equation in Fig. 6, the rapid decrease of the solitons to the Boussinesq system allows to study numerically soliton interactions. We consider as initial data the solitons with $c=1.05$ and $c=1.1$, the latter being centered at $x=-3$. The solution to the Boussinesq system for these initial data can be seen in Fig. 18. Visibly both $u$ and $\eta$ show the behavior of the KdV two-soliton. Note that the solution is computed in a frame commoving with $c=1.05$.

However the solution at $t=100$ in Fig. 19 shows upon closer inspection, in particular the close-up on the right, that there is dispersive radiation propagating to the left. Thus despite a similarity to the KdV two-solitons, the Boussinesq systems does not appear to be completely integrable.

8.3. Numerical study of the Boussinesq system for more general

initial data

In Fig. 20 we show the solution to the Boussinesq system (6.1 for the initial $u(x,0)=0$ and $\eta(x,0)=10\exp(-x^{2})$. The solution breaks at $t=0.4115$ since the fit of the Fourier coefficients to (7.4) indicates that a singularity in the complex plane hits the real axis. The fitted coefficient $\mu=0.345$ implies the formation of a cusp. The solution at the critical time is shown in Fig. 20. The behavior is similar to the cusp formation in solutions to the semiclassical NLS equation, see e.g. [17].

For $\epsilon=0.1$ and the same initial data, the solution breaks at $t=2.2262$ with $\mu=.423$, i.e., again a cusp as shown in Fig. 21. But this time the behavior is different from Fig. 20: the solution follows for a certain time the underlying wave equations and two humps are forming from the initial hump. But dispersion is too weak to overcome the nonlinearity and cusps form in finite time. Note that there is no clear scaling for the blow-up time in dependence of $\epsilon$ here since the mechanisms for the blow-up appear to be different.

For even smaller values of $\epsilon$, the solution for the same initial data appears to be global in time. As can be seen in Fig.22 for $\epsilon=0.01$, just two ‘solitons’ appear to emerge from the initial hump. For sufficiently small $\epsilon$, the system appears to be close to the wave equation.

For negative $w$, the system (6.1) is ill-posed. For the numerical studies, this implies that the Krasny filter must be used at a higher level (we choose $10^{-10}$) in order to control numerical instabilities. If we consider the initial data $u(x,0)=0$, $\eta(x,0)=-10\exp(-x^{2})$, i.e., the same data which led to Fig. 20 except for the sign of $w$, we find that the solution has a blow-up for $t=0.17$ as shown in Fig. 23. The fitting of the Fourier coefficients according to (7.4) yields $\mu=0.2029$, i.e., again a cusp.

For the same initial data, but $\epsilon=0.1$, a blow-up very similar to Fig. 23 is observed for $t=0.579$ with $\mu=1.15$. In the case of even smaller $\epsilon=0.01$, the solution appears to be again close to the solution of the wave equation as in Fig. 22.

The effects of the ill-posedness of the system for negative initial data are more visible in the presence of oscillations. In Fig. 24 we show the solution to the system for the initial data $u(x,0)=0$ and $\eta=\sin(10x)e^{-x^{2}}$. The code breaks at $t=0.0898$ since a singularity in the complex plane appears to hit the real axis. The fitting of the Fourier coefficients according to (7.4) is not conclusive ($\mu=0.0619$) whether this is a cusp or a pole.

8.4. Numerical study of the Boussinesq system with surface

tension

If we study the same initial data as in the previous subsection for the Boussinesq system with surface tension, the effects of the stronger dispersion are clearly visible.

The initial data of Fig. 23 lead in the presence of surface tension with $\beta=1$ to the solution in Fig. 25. The solution appears to exist for all times, instead of a shock a dispersive shock wave is observed.

For negative $\eta(x,0)$, the situation in Fig. 23 on the other hand is not really changed. Again a blow-up appears in finite time despite a surface tension of $\beta=1$ as can be seen in Fig. 26. The code breaks at $t=0.2568$ with $\mu=-0.05$ after fitting the Fourier coefficients according to (7.4).

9. Conclusion

We have tried in this paper to answer some key questions on the Whitham equation and systems. As far as modeling of water waves is concerned, the Whitham equation is relevant in the KdV (Boussinesq) regime and then it is probably not better than the KdV equation itself (and actually it is never used in realistic water waves modeling). On the other hand, because of its dispersion relation that behaves drastically differently for large and small frequencies, it has a fascinating variety of dynamical behaviors, most of them displayed in numerical simulations, that deserve further mathematical investigations.

The Boussinesq-Whitham systems, however, appear to have almost no interest for the modeling of water waves, since the local Cauchy problem (in absence of surface tension) can be solved only for initial data on the elevation that satisfies a rather unphysical condition. Its mathematical interest is also limited since its long wave limit is an ill-posed Boussinesq system.

The Cauchy problem for the system with surface tension while having similar shortcomings to model water waves, leads however to an interesting open mathematical question.

Acknowledgements**.**

The Authors were partially supported by the Brazilian-French program in mathematics and the MathAmSud program. J.-C. S. acknowledges support from the project ANR-GEODISP of the Agence Nationale de la Recherche. F.L and D.P. were partially supported by CNPq and FAPERJ/Brazil. J.-C. S. thanks David Lannes for useful discussions related to this work.