Computer Algebra Applied to a Solitary Waves Study

Didier Clamond, Denys Dutykh, André Galligo

► To cite this version:

Didier Clamond, Denys Dutykh, André Galligo. Computer Algebra Applied to a Solitary Waves Study. ISSAC'2015, Jul 2015, Bath, United Kingdom. 10.1145/2755996.2756659 . hal-01255412

HAL Id: hal-01255412

https://inria.hal.science/hal-01255412

Submitted on 13 Jan 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Computer Algebra Applied to a Solitary Waves Study

Didier CLAMOND*; Denys DUTYKH†; André GALLIGO‡

## Abstract

We apply computer algebra techniques, such as algebraic computations of resultants and discriminants, certified drawing (with a guaranteed topology) of plane curves, to a problem in fluid dynamics: We investigate “capillary-gravity” solitary waves in shallow water, relying on the framework of the Serre-Green-Naghdi equations. So, we deal with two-dimensional surface waves, propagating in a shallow water of constant depth. By a differential elimination process, the study reduces to describing the solutions of an ordinary non linear first order differential equation, depending on two parameters. The paper is illustrated with examples and pictures computed with the computer algebra system Maple.

## 1 Introduction

In this paper we apply computer algebra techniques such as algebraic computations of resultants and discriminants, certified drawing (with a guaranteed topology) of plane curves, to a problem in fluid dynamics. More precisely, we investigate “capillary-gravity” solitary waves in shallow water, relying on the framework of the Serre-Green-Naghdi equations. So, we deal with two-dimensional surface waves, propagating in a shallow water of constant depth (without assuming small wave amplitudes), see [6, 7, 11, 23]. Figure 1 represents two successive positions of a water wave, the bottom level and the still water level are emphasized (in our study the bottom will be horizontal), the velocities are indicated by small arrows. Beside the kinematic energy, the governing equations take into account an effect of capillarity pressure quantified by a Bond number $B$ (depending on a tension surface coefficient) and an effect of gravity force, jointly quantified by a (squared) Froude number $F$, which also indicates the wave phase velocity. To ease the presentation and for physical reasons, we can concentrate on the cases where $F \in [0, 2]$ and $B \in [-1, 2]$. The variable $x$ indicates the abscissa in the channel, the principal unknowns are the

<footnote>*Mathématiques, University of Nice-Sophia Antipolis</footnote>

<footnote>†LAMA, UMR 5127 CNRS</footnote>

<footnote>‡corresponding author, INRIA and University of Nice-Sophia Antipolis</footnote>
<figure><img src="image_1.png" /><figcaption>Figure 1: Water wave</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2: x-velocity</figcaption></figure>

elevation $h(x)$ of the free surface above $x$, and $\bar{u}(x)$ the depth-averaged horizontal velocity. The height of the surface at rest is denoted by $d$; $h = d$ is therefore a trivial solution. To illustrate some typical features of a water wave we have shown, in Figure 2, a color representation of the (normalized) variations of the horizontal component of the velocity $u(x, y)$, the coordinates $x$, and $y$ are also normalized. The bottom surface is assumed horizontal, and we only consider a (vertical) section of the wave. One can also see the graph of $y = h(x)$ which forms the elevation curve. It is called free curve (or more generally free surface) because unlike the bottom it is not prescribed.

The primitive equations of the general motion were established already in the 18th century by the great mathematicians Euler, Bernoulli and d'Alembert relying on physical conservative laws: conservations of mass and energy. However these non linear partial differential equations are difficult to solve and need to be contextualized (e.g. in the shallow water regime). Usually, the physicists consider that some quantity is small and replace the equations by the first terms of a Taylor approximation, see e.g. [20] for a mathematical general presentation. Here we follow the path opened by Lord Kelvin and his followers, and choose the equations of gravity waves first derived by Serre [23] in 1953, and later rediscovered by Green and Naghdi [11], which uses the depth-averaged horizontal velocity. As in [8], we generalize these equation by adding a term which takes into account the tension forces induced by capillarity. So the PDEs depend on the gravity $g$ and a surface tension coefficients $\tau$. In the second section, we present these non linear differential equations with boundary conditions and we explain that by considering a mobile frame, we are reduced to investigate their steady waves. Then, by a first integration and a simple substitution eliminating the variable $\bar{u}(x)$, the study reduces to describing the solutions of a system of ordinary differential equations depending on two parameters, $F$ and $B$. By a differential elimination process, inspired by physics considerations, we were able to reduce it further to a non linear first order differential equation. In this paper, we concentrate on the description of solitary waves, the steady waves which decay to a constant value at infinity. Solitary waves play a central role in nonlinear sciences, see [21]. For some special PDEs, such as KdV, NLS and
classical Serre-Green-Naghdi equations, the solitary waves can be found ana-
lytically; this is not the case for the equations, with capillarity effects, that we
consider. But we will get certified qualitative results. In the third section, we
review the different kinds of solutions of differential equations and emphasize the
interest of singular solutions, first investigated algebraically by Hamburger in
1893 [12] and studied more recently in computer algebra by Hubert [13, 14, 15] in
a more general framework. The angular solitons were first considered by Stokes
in 1880 [24] and received a renewed interest with the many studies on Camassa–
Holmes “peakons” [3], see also [19]. Then, we explain our phase plane analysis:
to the first order non linear ODE, depending on two parameters, we associate
a family of plane curves; and we translate the solution of the fluid mechanics
problem into a geometrical problem. In the fourth section, we make judicious
use of resultants and discriminants computations to get a cell decomposition
of the parameter space to classify the different behaviors of the (generalized)
solutions.

## 2 Formulation of the problem

The flow is assumed to be two-dimensional and irrotational, we neglect viscosity and compressibility of the fluid. The bottom is assumed to be horizontal. We will take a coordinate system (x, y) moving with the wave at the same speed. The x-coordinate is taken horizontally to the right and the y-coordinate is taken vertically upward, the bottom equation is y = 0. In this moving frame, the wave profile is stationary and there is an underlying flow traveling in the opposite direction. Let us recall that following Newton’s law of motion, the description of the flow is governed by the Euler equations in the plane which express mass and energy conservations. They are constraints on the velocity (u, v) and the pressure p. The pressure can be replaced by p = p₀ + τK, where p₀ is the pressure of the air above the free boundary, τ is the surface tension coefficient and K is the curvature of the free boundary.

2.1 Approximation in shallow water regime

In shallow water regime, following Serre, Green and Naghdi, we assume that
the horizontal velocity $u(x,y)$ can be approximated by its averaged value $\bar{u}(x)$
over the water column, so in Euler equations $u(x,y)$ is replaced by $\bar{u}(x)$; $\bar{u}_x$
denotes the derivative with respect to $x$ (more generally a subscript denotes
a derivative). Applying the incompressibility hypothesis, we can also replace
the vertical velocity $v(x,y)$ by $-y \bar{u}_x$. Finally, the new variables characterizing
the flow are the elevation of the free surface $h(x)$ and the averaged horizon-
tal velocity $\bar{u}(x)$. Time dependent Euler equations are also approximated and
become the following two PDEs, called (gravity-capillarity) Serre-Green-Naghi
equations:

$$
h_t + [h \bar{u}]_x = 0,
$$

$$
[h \bar{u}]_t + \left[ h \bar{u}^2 + \frac{1}{2} g h^2 + \frac{1}{3} h^2 \gamma - \tau R \right]_x = 0, \quad (1)
$$
$g$ being the acceleration due to gravity and

$$R = h h_{xx} (1 + h_x^2)^{-3/2} + (1 + h_x^2)^{-1/2},$$ $$\gamma = h (\bar{u}_x^2 - \bar{u}_{xt} - \bar{u} \bar{u}_{xx}),$$

$\gamma$ being the vertical acceleration of the fluid at the free surface.

In the second equation, which expresses energy conservation, we note that $\gamma$ contains a derivative in $t$. To concentrate all such derivatives in the first bracket, we consider two conservation laws of other momentums. Formally this amounts to add and subtract in equation (1) a third of $(h^3 \bar{u}_x)_{xt}$ to $h\bar{u}$ in order to get the following equation (2).

$$\left[ h\bar{u} - \frac{(h^3 \bar{u}_x)_x}{3} \right]_t + \left[ h\bar{u}^2 + \frac{gh^2}{2} - \frac{2h^3 \bar{u}_x^2}{3} - \frac{h^3 \bar{u} \bar{u}_{xx}}{3} - h^2 h_x \bar{u} \bar{u}_x - \tau R \right]_x = 0,$$

$$\left[ \bar{u} - \frac{(h^3 \bar{u}_x)_x}{3h} \right]_t + \qquad (3)$$

$$\left[ \frac{\bar{u}^2}{2} + gh - \frac{h^2 \bar{u}_x^2}{2} - \frac{\bar{u} (h^3 \bar{u}_x)_x}{3h} - \frac{\tau h_{xx}}{(1 + h_x^2)^{3/2}} \right]_x = 0.$$

Similarly considering the derivative of the similar quantity divided out by $h$, we get equation (3). Equations (2) and (3) also imply energy conservation, so we can replace the second equation of (1) by equations (2) and (3).

Since our aim is to study steady waves, for now on, we set to zero all derivatives with respect to $t$. Integrating the first equation of (1), we obtain a constant that we set to $-cd$. For a physical interpretation, we consider a $2L$-periodic solution, we have:

$$d = \frac{1}{2L} \int_{-L}^{L} h \, dx, \quad -cd = \frac{1}{2L} \int_{-L}^{L} h \bar{u} \, dx, \qquad (4)$$

thus $c$ is the wave phase velocity observed in the frame of reference without mean flow. This extends to the case of solitary waves, by letting $L$ tend to infinity. The values $c, g$ and $\tau$ serve to define three constants used, in fluid mechanics, to classify the waves: the (squared) Froude number $F = c^2/gd$, the Bond number $B = \tau/gd^2$ and the Weber number $W = \frac{B}{F} = \tau/c^2d$.

The mass conservation of (1) yields

$$\bar{u} = -cd/h, \qquad (5)$$

hence

$$\bar{u}_x = \frac{cd}{h^2}, \quad \bar{u}_{xx} = -2 \frac{cd}{h^3}. \qquad (6)$$
Then, replacing $R$ by its value, substitutions of $\bar{u}$ and its derivatives into (2) and (3), followed by an integration, give

$$ \frac{\mathcal{F}d}{h} + \frac{h^2}{2d^2} + \frac{\gamma h^2}{3gd^2} - \frac{Bh h_{xx}}{(1+h_x^2)^{\frac{3}{2}}} - \frac{B}{(1+h_x^2)^{\frac{1}{2}}} = \mathcal{F} + \frac{1}{2} - \mathcal{B} + K_1. \quad (7) $$

$$ \frac{\mathcal{F}d^2}{2h^2} + \frac{h}{d} + \frac{\mathcal{F}d^2h_{xx}}{3h} - \frac{\mathcal{F}d^2h_x^2}{6h^2} - \frac{Bd h_{xx}}{(1+h_x^2)^{\frac{3}{2}}} = \frac{\mathcal{F}}{2} + 1 + \frac{\mathcal{F}K_2}{2}, \quad (8) $$

with

$$ \gamma/g = \frac{\mathcal{F}d^3 h_{xx}}{h^2} - \frac{\mathcal{F}d^3 h_x^2}{h^3}, \quad (9) $$

the $K_1$ and $K_2$ are integration constants.

In the case of solitary waves, which interest us here, all these derivations remain correct with $K_1 = K_2 = 0$.

## 2.2 Solitary waves

Letting $K_1 = K_2 = 0$, in equations (7) and (8), we obtain two nonlinear differential equations governing the solitary waves solutions.

A key remark is that once we replace $\gamma$ by its value, the terms in $h_{xx}$ of these two equations are similar.

Computing $(7)-(h/d) \times (8)$ in order to eliminate $h_{xx}$ between (7) and (8), one obtains

$$ \frac{\mathcal{F}d}{2h} - \frac{h^2}{2d^2} - \frac{\mathcal{F}dh_x^2}{6h} - \frac{B}{(1+h_x^2)^{\frac{1}{2}}} = \mathcal{F} + \frac{1}{2} - \mathcal{B} - \frac{(\mathcal{F}+2)h}{2d}, \quad (10) $$

that is a first-order differential equation for $h$.

Normalizing $d$ to $d=1$, writing $h' = \frac{dh(x)}{dx}$ and multiplying by $h$, equation (10) becomes

$$ F(h', h) := \frac{\mathcal{F}h'^2}{3} + \frac{2Bh}{(1+h'^2)^{\frac{1}{2}}} - \mathcal{F} + (2\mathcal{F}+1-2\mathcal{B})h - (\mathcal{F}+2)h^2 + h^3 = 0. \quad (11) $$

Notice that, by taking the derivative of this equation with respect to $x$, assuming that $h$ is twice derivable, we obtain a system equivalent to the system formed by equations (7) and (8). In other words, starting from the first equations (1) in their steady solitary wave form, our manipulations amount to perform a differential elimination but directed by a physical intuition.

Now, we are reduced to a mathematical problem: study the solutions of the nonlinear first order differential equation $F(h', h) = 0$.

In a strict physical interpretation only the twice derivable functions $h(x)$ (such that the curvature of the wave, thence the capillarity force are well defined) should be considered. However in the next section, we will see that it is worthwhile to also consider less regular solutions of this equation; they are sometimes called singular solutions.
# 3 Generalized solutions

## 3.1 Corner waves

The theory of corner waves goes back to G.G. Stokes [24] who, in 1880, conjectured that finite amplitude crest gravity waves are limited by a wave with a corner at its crest; he proved that this corner has a 120 degrees angle. It also can be proved that the two branches at the crest of this highest corner wave have vanishing curvature. In other words, if $y = h(x)$ is the equation of the free surface and the crest occurs at $x = 0$ then $h'(0_{-}) = -h'(0_{+}) = \frac{\sqrt{3}}{3}$ and $h''(0_{-}) = h''(0_{+}) = 0$ (see Figure 3). For a presentation of Stokes waves see [10] and Wikipedia or another encyclopedia for a brief introduction.

<figure><img src="image_3.png" /><figcaption>Figure 3: Stokes' limit wave</figcaption></figure>

More recently, a famous equation was introduced by Camassa and Holm [3] to provide another simple model for waves in shallow water,

$$u_t + 2\kappa u_x - u_{xxt} + 3uu_x = 2u_x u_{xx} + uu_{xxx}.$$

It depends on a parameter $\kappa$ which, in this context, should be positive so that the solitary wave solutions are smooth solitons.

When $\kappa$ is zero, the Camassa-Holm this equation has so-called “peakon” solitary wave solutions. The wave slope at the peak does not vanish but jumps from one value to the opposite value. In other words, they are shaped like the graph of the function $h(x) = \exp(-|x|)$. Moreover, explicit formulas for the peakon interactions are known, [22]. There is an extensive literature on Camassa-Holm equation.

In 2013, a unified wave model for progressive gravity waves in finite water depth was proposed by Liao [19]. Based on the symmetry and the exact wave equations, it admits not only all traditional smooth periodic solitary waves but also the peaked solitary waves including the famous peaked solitary waves of Camassa-Holm equation mentioned above.

In all these cases, the corner wave can be viewed as a limit of a sequence of regular functions $h_{\kappa}(x)$ solutions of a differential equation depending on $\kappa$.
## 3.2 Singular solutions

Singular solutions of an ordinary differential equation of the form $G(x, y, dy/dx) = 0$ are not always clearly defined. As written in [17], in the classical treaties (see [16]) the answer is informal: a “general solution” is defined to be a one parameter family of solutions and a singular solution is a solution which is not contained in that family, see also [5]. Clearly if $G$ depends on a parameter $\lambda$ the envelope or the limit of the graphs of a family of solutions will be the graphs of general solutions. Note that these functions need not have singular points. These are the generalized solutions usually admitted under the vocable of singular solutions. See [17] for a formalization of these concepts and also the book of V. Arnol'd [1].

In computer algebra, the works of E. Hubert [13, 14, 15] present a computational view on general as well as singular solutions, that have enveloping properties of polynomial ordinary differential equations. We can also cite the work initiated by M. Hamburger already in 1893 [12].

Here, we consider another kind of generalized solutions which could also be called singular.

## 3.3 Our generalized solutions

The differential equations governing the steady waves of the (gravity-capillarity) Serre-Green-Naghi equations, (7) and (8) presented in Section 2, involve the derivative $h_x$ only by its square. Since these equations express the conservative laws, the non regular solutions $h(x)$ such that the derivative jumps at some points, but such that $h_x(x)^2$ and $h_{xx}(x)$ are continuous have a physical interest.

More precisely, in the sequel of this paper, we will investigate both the regular solutions and the “peakon” or angular solutions $h(x)$ such that $h_x$ is differentiable except for a finite number of values $x_i$, $i = 1..N$, at which only $h_x(x)^2$ is continuous and differentiable and $h_{xx}(x_i) = 0$.

Unlike the singular solutions of the previous subsection, these peakon solutions are not isolated and we do not know yet if they are limit solutions of a system of differential equation generalizing the Serre-Green-Naghi equations and involving more parameters.

## 3.4 Phase plane analysis

In order to investigate the number and behavior of the solutions and peakon solutions (defined in the previous subsection) of the differential equation $F_{\mathcal{F},\mathcal{B}}(h, h') = 0$, with respect to the pair of parameters $(\mathcal{F}, \mathcal{B}) \in [-1, 2]^2$, our approach is to describe graphically the variations of $h$ and $h'$, through the corresponding family of real algebraic curves $C_{\mathcal{F},\mathcal{B}} \in \mathbb{R}^2$. The algebraic implicit equation of $C_{\mathcal{F},\mathcal{B}}$ is $F_{\mathcal{F},\mathcal{B}}(h, h') = 0$. We will decompose the parameters space $(\mathcal{F}, \mathcal{B})$, with our restriction to the square $[0, 2]^2$, into sub-domains where the plane curves $C_{\mathcal{F},\mathcal{B}}$ have the same “shape” (in particular the same topology). Some sub-domains can be very small, to overcome this computational difficulty, the study of the
curves delimiting the sub-domains will be performed with a certified topology
(see, e.g. [9] or [4] and the references therein).

Let us illustrate our approach with an example: the phase plane description
for some values of the parameters that we chose randomly, $\mathcal{F} = 0.4$ and $\mathcal{B} = 0.9$.
The corresponding curve $C_{0.4,0.9}$ is shown on the left of Figure 4. It can be
plotted by simply using a discretization, e.g. using Matlab or drawn with a
guaranteed topology using specialized programs such as the a algebraic curves
package of Maple or Axel [2]. In particular, this means that ovals of whatever
size are preliminary detected relying on exact (given certain precision) numerical
and algebraic computations.

After computing the points of the curve with a horizontal tangent and the
singular points, we can decompose the curve into a finite number *N* of portions
of graphs $\Gamma_j, j = 1..N$ of type $h = \phi_j(h')$, for some differentiable functions
$\phi_j, j = 1..N$. Moreover, since when $h'$ is positive (resp. negative), $h$ should
increase (resp. decrease), so each $\Gamma_j$ can be oriented. See Figure 4 where there
are *N* = 6 oriented portions of graphs. The points with a horizontal tangent or
a singular point satisfy the condition $\partial F / \partial h' = 0$.

Then, on each branch, the corresponding differential equation $dx = \frac{dh}{\phi(h)}$ can
be solved numerically. It provides, on the corresponding interval, an approxi-
mate solution $h(x)$ which is qualitatively correct.

In this example, there are four points such that $\partial F / \partial h' = 0$: $A_1 = (0, 1)$, $A_2 = (0, 0.4)$, $A_3 \approx (2.1, 1.955)$ and $A_4 \approx (-2.1, 1.955)$. Call $B_1$ and $B_2$ the intersections of the curve with the $h'$-axis. Then, we get the following going up branches: $\Gamma_1$ from $B_1$ to $A_3$, $\Gamma_2$ from $A_2$ to $A_1$ (on the right side) and $\Gamma_3$ from $A_1$ to $A_3$. Symmetrically, the down-going branches $\Gamma_j, j = 4, 5, 6$ are in the left half-plane $h' < 0$.

Constrained by the asymptotic conditions, we only consider loops starting
from $A_1$ and arriving at $A_1$. In this example, the only continuous one is going
down on $\Gamma_5$ from $A_1$ to $A_2$ (on the left side), then going up on $\Gamma_2$ from $A_2$ to
$A_1$ (on the right side). It will correspond to the unique differentiable solution
$h(x)$ of the differential equation.

As we said in the previous subsection, we will also consider other solutions
that we call **angular solutions**.

Note that there are 6 points of the curve $C_{0.4,0.9}$ with a vertical tangent,
they are symmetric with respect to the $h$-axis. We call $V_1, V_2, V_3$ three of these 6
points, their coordinates are approximately $V_1 = (0.3, 0.5)$, $V_2 = (3.5, 1.5)$, $V_1 =
(-0.3, 0.5)$. The two semi second derivatives at these points are zero. The points
$V_2$ belongs to $\Gamma_1$ which is not connected to $A_1$, so we do not consider it; but
$V_1$ and $V_3$ are connected to $A_1$ by respectively $\Gamma_2$ and $\Gamma_5$. So, we can form
a discontinuous symmetrical loop going down from $A_1$ to $V_3$ (on $\Gamma_5$), jumping
from $V_3$ to $V_1$, then going up from $V_1$ to $A_1$ (on $\Gamma_2$). It corresponds to an
interesting new solitary wave with an angle, a “peakon”, $h(x)$ which satisfies
all the conservative laws. We can as well consider a loop going down from $A_1$
to $A_2$ (on the left side), then going up from $A_2$ to $V_1$, then jumping from $V_1$
to $V_3$, then cycling through $V_3, A_2, V_1, V_3$, then going up from $A_2$ to $A_1$. The
<figure><img src="image_4.png" /><figcaption>Figure 4: Curve $C_{0.4,0.9}$ of Example 1</figcaption></figure>

corresponding wave has several “peaks”, see Figure 5. The angular points are emphasized by green solid disks.

## 4 A family of plane curves

By the discussion in the previous section, we are led to study and classify the family of curves $C_{\mathcal{F},\mathcal{B}}$, in the real $(k, h)$ plane, defined by the algebraic (but non polynomial) equation $F(k, h) = 0$ with exact coefficients

$$ F(k, h) = \frac{\mathcal{F}k^2}{3} + \frac{2\mathcal{B}h}{(1+k^2)^{\frac{1}{2}}} - \mathcal{F} + (2\mathcal{F} + 1 - 2\mathcal{B})h - (\mathcal{F} + 2)h^2 + h^3 = 0, \quad (12) $$

To compute the points with horizontal or vertical tangent and the singular points, we will use resultants after changes of coordinate.

We have the following partial derivatives of $F$ with respect to $h$ and $k$

$$ F_h = \frac{2\mathcal{B}}{(1+k^2)^{\frac{1}{2}}} + (2\mathcal{F} + 1 - 2\mathcal{B}) - 2(\mathcal{F} + 2)h + 3h^2, $$

$$ F_k = \frac{2\mathcal{F}k}{3} + \left( \frac{2\mathcal{B}kh}{(1+k^2)^{\frac{3}{2}}} \right). $$
<figure><img src="image_5.png" /><figcaption>Figure 5: Multiangular wave</figcaption></figure>

## 4.1 Points on the h-axis

The points on the h-axis are important; to determine them it suffices to substitute $k = 0$ and solve the univariate equation $F_{\mathcal{F},\mathcal{B}}(0, h) = 0$, which gives either $h = 1$ with multiplicity 2, or $h = \mathcal{F}$. Call $A_1$ and $A_2$ the points of coordinates $(0, 1)$ and $(0, \mathcal{F})$.

In order to investigate the local behavior at $A_1$, we compute the second-order Taylor expansion of $F$ at that point

$$ (\mathcal{F} - 3\mathcal{B}) k^2 - 3(\mathcal{F} - 1)(h - 1)^2 = 0. $$

Therefore, if $(\mathcal{F} - 1)(\mathcal{F} - 3\mathcal{B}) < 0$, $A_1$ is an isolated point of the curve $C_{\mathcal{F},\mathcal{B}}$; this implies that the only solution is the trivial one $h = 1$. If $(\mathcal{F} - 1)(\mathcal{F} - 3\mathcal{B}) > 0$, $A_1$ is a double point corresponding to the crossing of two branches of $C_{\mathcal{F},\mathcal{B}}$; hence there is no obstruction (at the level of this preliminary local analysis) for the existence of a solitary wave. If $\mathcal{F} = 1$, then $A_1 = A_2$ and we must consider the third-order Taylor expansion

$$ (1 - 3\mathcal{B}) k^2 + 3(h - 1)^3 - (h - 1) k^2 = 0. $$

Therefore if $3\mathcal{B} > 1$, (resp. $3\mathcal{B} < 1$), $C_{\mathcal{F},\mathcal{B}}$ admits at $A_1$ a cusp above (resp. below) $A_1$; while if $3\mathcal{B} = 1$, the Taylor expansion corresponds to three lines including $h = 1$; i.e., no obstruction (at this local analysis) for the existence of a (non constant) solitary wave.

Now, a Taylor expansion of $F$ at the point $A_2$, i.e. $(k = 0, h = \mathcal{F})$, yields

$$ 3(\mathcal{F} - 1)^2 (h - \mathcal{F}) = \mathcal{F}(3\mathcal{B} - 1) k^2. $$
So, when {$\mathcal{F} \neq 1, 3\mathcal{B} \neq 1$}, we get a point with an horizontal tangent; it is of convex type if $3\mathcal{B} > 1$ and concave if $3\mathcal{B} < 1$. Therefore, if $\mathcal{F} < 1$ and $3\mathcal{B} > 1$, or if $\mathcal{F} > 1$ and $3\mathcal{B} < 1$, there is no obstruction (at this local analysis level) for the existence of a solitary wave. While if $\mathcal{F} < 1$ and $\mathcal{F} < 3\mathcal{B} < 1$, or if $\mathcal{F} > 1$ and $\mathcal{F} > 3\mathcal{B} > 1$, this local analysis only permits an angular solitary wave.

When $\mathcal{F} \neq 1$ and $3\mathcal{B} = 1$, the Taylor expansion of $F$ at the point $A_2$ yields

$$ \mathcal{F}(h')^4 \sim 4(\mathcal{F}-1)^2(h-d). $$

In other words, the curve has a flat point (the curvature vanishes).

Notice that at the end of this first local analysis, our discussion partitioned the parameter plane $(\mathcal{F}, \mathcal{B})$ into domains delimited by the three lines $\mathcal{F} = 1$, $\mathcal{B} = \frac{1}{3}$ and $\mathcal{F}-3\mathcal{B} = 0$, and discarded the domain defined by $(\mathcal{F}-1)(\mathcal{F}-3\mathcal{B}) < 0$.

### 4.1.1 Points with a horizontal tangent

These points satisfy the two equations $F = 0$ and $\partial F / \partial k = 0$. For $k \neq 0$, the second equation vanishes when $\mathcal{F}^2(1 + k^2)^3 = (3\mathcal{B}h)^2$. From which we express $k^2$ as a function in $h$ and replace this expression into the equation $F = 0$. We get an equation depending only on $h$, but with a cubic root operator. To get rid of it, and deal with polynomials, we introduce a new variable $Y$ such that $h = (\mathcal{F}/3\mathcal{B})Y^3$, thence $k^2 = Y^2 - 1$ with $Y \ge 1$. We are led to the following polynomial equation in $Y$

$$ \begin{aligned} f(Y) = \mathcal{F}^2 Y^9 & - (3\mathcal{F}-2)\mathcal{F}\mathcal{B}Y^6 \\ & + 9\mathcal{B}^2(1+2\mathcal{F}-2\mathcal{B})Y^3 + 27\mathcal{B}^3Y^2 - 36\mathcal{B}^3 = 0. \end{aligned} \quad (13) $$

The adapted tool to discuss the number of real roots of $f$ is its discriminant $D_1(\mathcal{F}, \mathcal{B})$.

The real trace of the discriminant locus $D_1$ is an algebraic curve in the $(\mathcal{F}, \mathcal{B})$-plane, it characterizes the values of $(\mathcal{F}, \mathcal{B})$ where $f(Y)$ has (real or complex) multiple roots. A key property, for our analysis, is that this curve $D_1$ divides the parametric plane into domains where $f$ has a constant number of real roots. However, since we do not have a necessary and sufficient condition, two adjacent domains delimited by $D_1$ may correspond to the same number of real roots, they should be merged to form only one domain. The sub-domains of these domains delimited by the implicit curve defined by $f_{\mathcal{F},\mathcal{B}}(1) = 0$, in the parametric plane, have a constant number of real roots greater than 1. Also, two adjacent sub-domains may have the same number of real roots greater than 1 and should be merged. Eventually, we explicitly construct a partition into connex domains. In each such domain the number of local extrema of $h$ on $C_{\mathcal{F},\mathcal{B}}$ is fixed, and we are able to classify the possible shapes of $C_{\mathcal{F},\mathcal{B}}$.

### 4.1.2 Points with a vertical tangent

We follow roughly the same method to study the points of $C_{(\mathcal{F},\mathcal{B})}$ with a vertical tangent. They satisfy the two equations $F = 0$ and $\partial F / \partial h = 0$. Combining
these equations, and introducing a new variable $Z$ such that $k^2 = Z^2 - 1$ thence $\mathcal{F}(Z^2 - 1) = 3(h-1)(2h^2 - h - 1)$. Eliminating $h$ by the computation of a resultant, we get a polynomial $g(Z, \mathcal{F}, \mathcal{B})$, of degree 6 in $Z$, which plays the same role than the polynomial $f(Y)$ in the previous paragraph.

Then we compute the discriminant polynomial of $g$ with respect to $Z$, we find a polynomial $D_2(\mathcal{F}, \mathcal{B})$. As explained above, it allows to decompose the parameters space into domains where the number of points with a vertical tangent is constant. We also require that these points satisfy $h \ge 0$.

Then we take the intersection of the two families of sub-domains defined by $D_1(\mathcal{F}, \mathcal{B})$ and by $D_2(\mathcal{F}, \mathcal{B})$.

Eventually, we explicitly construct a “cell decomposition”. In each such cell the number of local extrema of $h$ and the number of local extrema of $k$ on $C_{\mathcal{F},\mathcal{B}}$ is fixed, and we are able to classify the possible shapes of $C_{\mathcal{F},\mathcal{B}}$, with respect to admissible paths connecting $A_1$ either to a point with $k=0$ or to a point with a vertical tangent.

### 4.1.3 Pictures

Using the commands `discrim()` and `factor()` of the computer algebra system Maple, we decompose the discriminant polynomial $D_1$ into the product of a square $(\mathcal{F} - 3\mathcal{B})^2$ by another polynomial $D$ of degree 10 in $(\mathcal{F}, \mathcal{B})$. The zero locus of $D$ (in red) together with the lines $\mathcal{B} = \frac{1}{3}, 3\mathcal{B} = \mathcal{F}$ and $\mathcal{F} = 1$ are shown in Figure 6.

We also decompose the other discriminant polynomial $D_2$, with respect to $Z$. We find the product of the same polynomial $D$ by powers of $\mathcal{F}$, powers of $\mathcal{B}$, and by the cube of another polynomial $D_3(\mathcal{F}, \mathcal{B})$. The zero locus of $D_3$ is shown in green in Figure 6.

The line defined by the equation $\mathcal{F} - 3\mathcal{B} = 0$ is of special interest since it corresponds to a double root at a fixed value $Y = 1$, which corresponds to the double point at $h = 1$ and $k = 0$, and another simple solution $Y(\mathcal{B})$.

### 4.1.4 Classification and special behaviors

Considering each of the delimited domains, allow to classify all the possible curves $C_{\mathcal{F},\mathcal{B}}$, by the relative location of their points with horizontal tangent, vertical tangent, singularity. For each such type of curve, we analyzed all the possible admissible paths (with respect to the orientation) departing from $A_1$ and arriving to a point with a vertical tangent, or to a point of the line $k = 0$. To each of them is associated a regular or an angular solitary wave.

We found only two types of regular solitary waves: a crest when $\mathcal{F} > 1$ and a through when $\mathcal{F} < 1$, they appear for some sub-domain of the square. For $(\mathcal{F}, \mathcal{B})$ in another sub-domain, we may have a multi angular solitary wave, they correspond to curves similar to the one in Figure 4.

Then in another parameter domain, we find curves similar to $C_{-0.2,0.32}$, see Figure 8 which allow only one angular solitary wave as in Figure 15. Here is a list of additional extremal behaviors.
1. The limiting case where $(\mathcal{F}, \mathcal{B})$ belongs to a branch of the discriminant locus e.g. when $(\mathcal{F} = 0.8, \mathcal{B} = 0.3538557)$ approximately, we obtain the curve shown in Figure 12. There is a a connex path but with an angle going from $A_1$ to a lower point on $h' = 0$. It gives rise to a weak regular wave shown in Figure 16, which has two points where the regularity is only $C^1$ (we represented them by changing the colors of the branches).

2. A tiny branch of the discriminant with $\mathcal{F} > \infty$ gives rise to another interesting type of “semi extremal curve” e.g. $C_{1.01,0.33341215}$, shown in Figure 13, it corresponds to (multi) angular weak crests, shown in Figure 17. The angles are indicated with a red bullet and the weak defect of regularity is indicated with a green bullet.

3. The most extremal curve is $C_{1,\frac{1}{3}}$ (see Figure 11) since it has a branch with an horizontal tangent passing through the point $A_1$ which means that the decreases at infinity is algebraic and not exponential. It admits a symmetric single angular solitary wave.

# 5 Conclusion

In this paper, we have presented for the computer algebra community, an actual problem coming from the fluid mechanics community. Indeed, the discovery of giant ocean waves, sometimes called freak or rogue waves, renewed the interest of that community in classical water wave problems and their “mathematical” aspects. We explained how the deep understanding of the solitary capillary-gravity waves in shallow water, including the investigation of singular or “angular” waves, could benefit from the mathematical description of the shapes of a family of plane real algebraic curves depending on two real parameters.

Then, relying on computer algebra techniques (adapted coordinates changes, geometric interpretation and computations of resultants and discriminants, certified graph drawing) we were able to list all the possible shapes of the phase diagram curves and all possible regular or “angular” solitary waves.

Of course, much work remains to be done. In particular, it would be worthwhile to investigate all the periodic steady waves. It would also be very interesting to derive from the primitive equation of motion, a set of equations depending on more (physical) parameters such that the angular solutions described here appear as limit of regular solutions of this new set of equations. Then, our angular solutions will become singular solutions in the more usual sense.

# References

[1] ARNOLD, V. I., Geometrical Methods In The Theory Of Ordinary Differential Equations, Springer-Verlag (1988).

[2] See HTTP://AXEL.INRIA.FR.
[3] CAMASSA, R.; HOLM, D.D., An integrable shallow water equation with peaked solitons, *Phys. Rev. Lett.* **71** (11): 1661–1664, (1993),

[4] CHENG, J., ET AL., On the topology of real algebraic plane curves. Mathematics in Computer Science, Springer, 4 (1), pp.113–137. (2010), https://hal.inria.fr/inria-0051717.

[5] COHN, R. The general solution of a first order differential polynomial. Proceedings of the American Mathematical Society 55,1, 14–16. (1976).

[6] DUTYKH, D., CLAMOND, D., MILEWSKY, P. & MITSOTAKIS, M., An implicit-explicit finite volume scheme for fully nonlinear Serre equations. Eur. J. App. Math. **24**, 761–787, (2013).

[7] DUTYKH, D. AND CLAMOND, D., Efficient computation of steady solitary gravity waves. *Wave Motion* **51**, 86–99, (2014).

[8] CLAMOND, D., DUTYKH, D. AND GALLIGO, A, Extreme solitary capillarity-gravity waves in shallow water. *In preparation*, (2015).

[9] GONZALEZ-VEGA, L. & NECULA, I., Efficient topology determination of implicitly defined algebraic plane curves. *Computer Aided Geometric Design* **19**, 9, 719–743, (2002).

[10] GRANT, M.A., Standing Stokes waves of maximum height, *Journal of Fluid Mechanics* 60 (3): 593604, (1973).

[11] GREEN, A.E AND NAGHDI, P.M, A derivation of equations for wave propagation of water of variable depth. *J. of Fluid Mechanics*, 78, 237–246, (1976).

[12] HAMBURGER, M. Ueber die singulären lösungen der algebraischen differentialgleichungen erster ordnung. *J. Reine Ang. Math.*, 112, 205–246, (1893).

[13] HUBERT, E., The general solution of an ordinary differential equation. *In ISSAC'96*. ACM Press. (1996).

[14] HUBERT, E., Detecting Degenerate Behaviors in Non-Linear Differential Equations of First order Algebraic Differential Equations. *Theoretical Computer Science*, vol 187 (1-2), pages 7–25, (1997).

[15] HUBERT, E. (1998), Essential components of an algebraic differential equation. *Journal of Symbolic Computations*, volume 28(4-5), 657–680, (1999).

[16] INCE, E., Ordinary Differential Equations. Dover (1956).

[17] IZUMIYA, S. AND YU, J. How to define singular solutions. Kodai Math. J.16, 227–234, (1993).

[18] LANG, S, Algebra. Chapter 5. em Addison-Wesley, (1984).
[19] LIAO, S.J., Do peaked solitary water waves indeed exist?, *Communications in Nonlinear Science and Numerical Simulation*, 19:1792–1821, (2014).

[20] OKAMOTO, I. AND SHOJI, M., The mathematical theory of permanent progressive water-waves. *Adv. Ser. Nonlin. Dyn.*, 20. Worlds Scientific. (2001).

[21] OSBORN, A., Nonlinear ocean waves and the inverse scattering transform. *Academic Press.*, 97, Elsevier (2010).

[22] PARKER, A., On the Camassa–Holm equation and a direct method of solution. III. N-soliton solutions, *Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.* 461: 3893–3911, (2005).

[23] SERRE, F., Contribution à l'étude des écoulements permanents et variables dans les canaux. *Houille Blanche*, 8, 374–388, (1953).

[24] STOKES, G.G., Supplement to a paper on the theory of oscillatory waves, *Mathematical and Physical Papers, Volume I, Cambridge University Press*, pp. 314–326, (1880).
<figure><img src="image_6.png" /><figcaption>Figure 6: Discriminant loci</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7: Zoom on a detail</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 8: $C_{-02,032}$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9: $C_{01,032}$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 10: $C_{096,032}$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 11: Extremal curve</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12: Curve $C_{08,03538}$ from the discriminant locus</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 13: Semi extremal curve</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 14: Regular</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 15: Angular</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 16: Weak regular wave</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 17: Weak multi angular crest wave</figcaption></figure>