André É Kabakouala
email: andre.kabakouala@lmpt.univ-tours.fr

Luc Molinet
email: luc.molinet@lmpt.univ-tours.fr

ON THE STABILITY OF THE SOLITARY WAVES TO THE (GENERALIZED) KAWAHARA EQUATION

In this paper we investigate the orbital stability of solitary waves to the (generalized) Kawahara equation (gKW) which is a fifth order dispersive equation. For some values of the power of the nonlinearity, we prove the orbital stability in the energy space H 2 (R) of two branches of even solitary waves of gKW by combining the well-known spectral method introduced by Benjamin [3] with continuity arguments. We construct the first family of even solitons by applying the implicit function theorem in the neighborhood of the explicit solitons of gKW found by Dey et al. [8]. The second family consists of even travelling waves with low speeds. They are solutions of a constraint minimization problem on the line and rescaling of perturbations of the soliton of gKdV with speed 1.

Introduction

We consider the generalized Kawahara equation (gKW):

(1.1)

∂ t u + u p ∂ x u + ∂ 3 x u -µ∂ 5
x u = 0, (t, x) ∈ R * + × R, with initial data u(0) = u 0 ∈ H 2 (R), where p ∈ N * denotes the power of nonlinearity, and µ > 0 the parameter which control the fifth-ordre dispersion term. In the case p = 1, this equation has been derived by Kawahara [11] as a model for water waves with weak amplitude in the long-wave regime approximation for moderate values of surface tension and a Weber number close to 1/3. For such Weber numbers the usual description of long water waves via the Korteweg-de Vries (KdV) equation fails since the cubic term in the linear dispersion relation vanishes and fifth order dispersion becomes relevant at leading order. Note that positive values of the parameter µ in (1.1) correspond to Weber numbers larger than 1/3.

The Cauchy problem associated to (1.1) is locally well-posed in H 2 (R) (see for instance [START_REF] Abdelouhab | Nonlocal models for nonlinear, dispersive waves[END_REF]). The H 2 -solutions of (1.1) satisfy the following two conservation laws in time:

(1.2) E p,µ (u(t)) = R µ 2 (∂ 2 x u) 2 (t) + 1 2 (∂ x u) 2 (t) - 1 (p + 1)(p + 2) u p+2 (t) = E µ (u 0 ) (energy)
and

(1.3) V (u(t)) = 1 2 R u 2 (t) = V (u 0 ) (mass).
These conserved quantities enable to extend the solutions for all positive times so that (1.1) is actually globally well-posed in H 2 (R). We are interested in solitary waves of gKW, i.e., the solutions to equation (1.1) of the form u(t, x) = ϕ c,p,µ (x -ct), traveling to the right with the speed c > 0. Substituting u by ϕ c,p,µ in (1.1), integrating on R with the assumption ∂ k x ϕ c,p,µ (±∞) = 0 for k = 0, . . . 4, we obtain (1.4) µ∂ 4

x ϕ c,p,µ (x) -∂ 2 x ϕ c,p,µ (x) + cϕ c,p,µ (x) = 1 p + 1 ϕ p+1 c,p,µ (x), ∀x ∈ R.

In [START_REF] Dey | Stationary solitons of the fifth order KdV-type. Equations and their stabilization[END_REF], Dey et al. compute explicit solutions to (1.4) that write :

(1.5) ϕ c,p,µ (x) = (p + 1)(p + 4)(3p + 4)c 8(p + 2) 2 µ .

and study their orbital stability. More generally, the existence of solutions to (1.1) with µ > 0 and has been proven in [START_REF] Levandosky | A stability analysis of fifth-order water wave models[END_REF] by solving an associated constrained minimization problem. The orbital stability of the set of ground states to this minimization problem is also studied. Note that (1.6) is a fourth order differential equation and, up to our knowledge, no uniqueness result (up to symmetries) is known even for ground state solutions.

On the other hand, for µ = 0, (1.7) becomes the well-known generalized Korteweg-de Vries equation (gKdV):

(1.7)

∂ t u + u p ∂ x u + ∂ 3 x u = 0, (t, x) ∈ R * + × R and thus, for µ = 0 in (1.4), we recover the equation of solitons of gKdV:

(1.8)

-∂ 2 x ϕ c,p,0 (x) + cϕ c,p,0 (x) = 1 p + 1 ϕ p+1 c,p,0 (x), ∀x ∈ R.

Recall that, the solitons ϕ c,p,0 of gKdV are unique, up to translations, and to the transformation: ϕ c,p,0 → -ϕ c,p,0 if p is even. Moreover, they are explicitly defined by:

(1.9) ϕ c,p,0 (x) = (p + 1)(p + 2)c 2

1/p sech 2/p p √ c 2
x , ∀x ∈ R.

In this paper, we start by constructing two branches of solutions to equation (1.4) with µ = 1 and establishing some uniqueness results on these solutions. Firstly, we fix p ∈ N * , µ = 1 and c p = 2 2 (p+2) 2 (p 2 +4p+8) 2 . By applying the Implicit Function Theorem in the neighborhood of the explicit solution ϕ cp,p , we construct a continuous in H 4 (R) branch {ϕ c,p , c ∈]c p -δ p , c 0 + δ p [}, with 0 < δ p ≪ 1, of even solutions to (1.4). For each c ∈]c 0 -δ p , c 0 + δ p [, ϕ c,p is the unique even H 4 (R) solution of (1.4) in some H 4 -neighborhood of ϕ cp,p . Secondly, for all p ∈ N * , c > 0 and µ > 0, following [START_REF] Levandosky | A stability analysis of fifth-order water wave models[END_REF], we minimize in the even functions of H 2 (R) the functional:

(1.10)

I c,µ (ψ) = R µ 2 (∂ 2 x ψ) 2 + 1 2 (∂ x ψ) 2 + c 2 ψ 2 ,
under the constraint:

(1.11) K p (ψ) = 1 (p + 1)(p + 2) R ψ p+2 = K p (ϕ c,p,0 ).
For µ > 0 small enough we prove the uniqueness of the associated even ground states by using ideas of [START_REF] Kenig | Local well-posedness and blow-up in the energy space for a class of L 2 critical dispersion generalized Benjamin-Ono equations[END_REF] . In this way, we construct a family {ϕ c,p,µ , 0 < µ ≪ 1} of even solutions to (1.4) such that

(1.12) lim µ→0 + ϕ c,p,µ -ϕ c,p,0 H 1 (R) = 0, where ϕ c,p,0 is defined by (1.8)-(1.9). Noticing that u is a solution to (1.4) with c = 1, p = p 0 ≥ 1 and µ > 0 if and only if v = µ 1/p u( √ µ•) is a solution to (1.4) with c = µ, p = p 0 and µ = 1, we obtain a H 1 -continuous branch {ϕ c,p,1 , c ∈]0, δ ′ p [} with 0 < δ ′ p ≪ 1
, of even solutions to (1.4) (with µ = 1) traveling with low speeds. For each c ∈]0, δ ′ p [, ϕ c,p is the unique even solution of the constraint minimizing problem (1.15).

The main result of this paper is the orbital stability of the solitary waves that form these two branches. Note that this improves earlier results (see [START_REF] Levandosky | A stability analysis of fifth-order water wave models[END_REF]) where the stability of the set of ground states is proven. Before stating our main result, let us recall the definition of orbital stability and define what we will call ground state solutions to (1.4). Definition 1.1 (Orbital Stability). Let φ ∈ H 2 (R) be a solution of (1.4). We say that φ is orbitally stable in H 2 (R), if for all ε > 0, there exists δ ε > 0, such that for all initial data u 0 ∈ H 2 (R), satisfying

(1.13) u 0 -φ H 2 (R) ≤ δ ε , the solution u ∈ C(R + , H 2 (R)) of gKW emanating form u 0 satisfies (1.14) sup t∈R+ inf z∈R u(t, • + z) -φ H 2 (R) ≤ ε.
In this paper, for s ≥ 0, we set

H s e (R) = {u ∈ H s (R) : u(-•) = u(•)} with L 2 e (R) = H 0 e (R)
. This space, endowed with the metric of H s (R) is an Hilbert space. Definition 1.2 (Even ground state solution). We say that a solution to (1.4) is an even ground state solution to (1.4), if it is also a solution to the constraint minimization problem

(1.15) S β c,p,µ = inf I c,µ (ψ) : ψ ∈ H 2 e (R), K p (ψ) = β . for some β > 0.
Our main results can be summarized as follows :

Theorem 1.1. Let us fix µ = 1 and set c p = 2 2 (p+2) 2
(p 2 +4p+8) 2 for any p ≥ 1. (i) For p ∈ {1, 2, 3, 4}, there exist δ p > 0 such that for any c ∈]c p -δ p , c p + δ p [, (1.4) has a unique even solution in the ball of H 4 (R) centered in ϕ cp,p,1 (explicitly defined in (1.5)-(1.6)) with radius δ p .

These solutions form a curve of class C 1 in H 4 (R), passing by ϕ cp,p,1 , of even solitary waves to (1.1) which are all stable. (ii) For p ∈ {1, 2, 3}, there exist δ ′ p > 0 such that, for any c ∈]0, δ ′ p [, there exists a unique even ground state ϕ c,p,1 to (1.4). These solutions form a H 1 -continuous curve of even solitary waves of (1.1) which are all stable.

Remark 1.1. If in Definition 1.2, we replace the requirement β > 0 by β = 0 then the uniqueness result for the even ground states to (1.13) with µ = 1 and c > 0 small enough still holds but up to the symmetry v → -v in the case p even.

The stability result will follow from a continuity argument together with a positivity property of the quadratic form associated with the second Fréchet derivative of the action functional at respectively ϕ 1,p,µp and ϕ 1,p,0 where µ p = 2 2 (p+2) 2 (p 2 +4p+8) 2 . Note that to prove this positivity property for ϕ 1,p,µp we check numerically a sign condition on some L 2 -scalar product linked to ϕ 1,p,µp .

We will apply the following classical proposition for equations of form the ∂ t V ′ (u) = ∂ x E ′ (u) where E and V are conservation laws, V being quadratic. Proposition 1.1 (see for instance de Bouard [START_REF] De Bouard | Equations dispersives non linéaires[END_REF]). Let φ c,p,µ ∈ H 2 (R) be a solution of (1.4) traveling with the speed c > 0, and L c,p,µ be the linearized operator associated to the second derivative of the action functional E p,µ + cV at φ c,p , defined by

L c,p,µ v = µ∂ 4 x v -∂ 2 x v + cv -φ p c,p,µ v, for all v ∈ H 4 (R). Assume that the exists δ > 0 such that (1.16) L c,p,µ v, v L 2 ≥ δ v 2 H 2 (R) , for all v ∈ H 2 (R) satisfying the orthogonalities (1.17) v, φ c,p,µ L 2 = v, φ ′ c,p,µ L 2 = 0. Then, φ c,p,µ is stable in H 2 (R).
The proof of Proposition 1.1 follows from the theory developed by Benjamin [START_REF] Benjamin | The stability of solitary waves[END_REF], Bona [START_REF] Bona | On the stability theory of solitary waves[END_REF] and Weinstein [19] which relies on the spectral properties of L c,p,µ (see Section 2 for details).

For other stability results on solitons the similar models, see for instance: Karpman [9]-[10], Dey et al. [START_REF] Dey | Stationary solitons of the fifth order KdV-type. Equations and their stabilization[END_REF], Levandosky [START_REF] Levandosky | A stability analysis of fifth-order water wave models[END_REF]- [START_REF] Levandosky | Stability of solitary waves of a fifth-order water wave model[END_REF], Bridges et al. [START_REF] Bridges | Stability and instability of solitary waves of the fifth-order KdV equation: a numerical framework[END_REF], Pava [START_REF] Angulo | On the instability of solitary-wave solutions for fifth-order water wave models[END_REF] .

This paper is organized as follows. In Section 2, we prove the assertion (i) of Theorem 1.1 and, in Section 3, we prove the assertion (ii) of Theorem 1.1.

Existence and Stability for the branch crossing the explicit solitary waves of gKW

In this section, we prove the point (i) of Theorem 1.1. First, we observe that u is a solution to (1.1) with µ > 0 if and only if u µ (t, x) = µ 1/p u(µ 3/2 t √ µx) is a solution to (1.1) with µ = 1. In particular, ϕ is a solution to (1.4) with c = 1 and µ > 0 if and only if

ϕ µ = µ 1/p ϕ( √ µ •) is a solution to (1.4)
with c = µ and µ = 1. This ensures that it is equivalent to prove the assertion (i) of Theorem 1.1 with µ = µ p = 2 2 (p+2) 2 (p 2 +4p+8) 2 and c close to 1. For the numerical checking of the sign of some L 2 -scalar product (see Subsection 2.4) we prefer to work with this last normalization in this section.

So, let p ∈ N * be fixed and set µ = µ p = 2 2 (p+2) 2 (p 2 +4p+8) 2 , then the explicit solitary wave given by (1.5) travels with the speed c = 1. The strategy is as follows. By applying the Implicit Function Theorem, we will establish for all real c > 0 sufficiently close to 1, the existence of even solitary waves ϕ c,p,µp of gKW traveling with speed c (see Lemma 2.2). Next, following Albert [START_REF] Albert | Positivity properties and stability of solitary-wave solutions of model equations for long waves[END_REF], we use that the hypotheses of Proposition 1.1 are fulfilled for the explicit solitary wave ϕ 1,p,µp whenever a sign condition on some L 2scalar product involving ϕ 1,p,µp is satisfied. This sign condition is checked numerically in Subsection 2.4. Finally, arguing by continuity, we prove that the hypotheses of Proposition 1.1 are still fulfilled for c close enough to 1 which leads to the orbital stability of all the the family ϕ c,p,µp , c ∈]1 -

δ ′ p , 1 + δ ′ p [ ⊂ H 4 (R), with 0 < δ ′ p ≪ 1.
We will divide the proof of Theorem 1.1 part (i) into several lemmas. We first study the spectral properties of the linearized operator L 1,p,µp associated with the second derivative of the action functional at the explicit soliton ϕ 1,p,µp (defined in (1.5)-(1.6)).

2.1. Spectral properties of L 1,p,µp . Lemma 2.1 (Spectral Properties of L 1,p,µp ). Let p ∈ N * . We consider the unbounded operator L 1,p,µp :

L 2 (R) → L 2 (R), defined by: u → µ p ∂ 4 x u -∂ 2 x u + u -ϕ p 1,
(R), since u H 4 (R) u L 2 (R) + L 1,p,µp u L 2 . Moreover, it is a compact perturbation of µ p ∂ 4 x -∂ 2 x + 1, since ϕ p
1,p,µp is smooth and decays exponentially to 0, and thus, its essential spectrum is given by σ

ess (L 1,p,µp ) = [1, +∞[.
According to Albert [START_REF] Albert | Positivity properties and stability of solitary-wave solutions of model equations for long waves[END_REF] (Theorem 3.2 and Lemma 10), to get the properties (P 2)-(P 3), it suffices to prove that:

(2.1) F (ϕ 1,p,µp )(ω) > 0, ∀ω ∈ R, and d 2 dw 2 logF (ϕ p 1,p,µp )(ω) < 0, ∀ω ∈ R * ,
where F denotes the spatial Fourier transform. Thanks to the computations done by Magnus and Oberhettinger [START_REF] Magnus | Formulas and theorems for the special functions of mathematical physics[END_REF] (p. 34), we first recall that

(2.2) F (sech ν (•)) (ω) = 2 ν-1 (Γ(ν)) -1 Γ ν 2 + i ω 2 2 ≥ 0, ∀ω ∈ R, ν ∈ R + ,
where Γ denotes the gamma function. On the other hand, using the usual Fourier transforms table, we establish that

(2.3) F sech 4 (•/2) (ω) = 2 4 π 3! ω(ω 2 + 1)cosech(πω), ∀ω ∈ R and thus (2.4) d 2 dω 2 log F sech 4 (•/2) (ω) = 2 4 π 3! - 1 ω 2 + 2(1 -ω 2 ) (1 + ω 2 ) 2 + π 2 sinh 2 (πω) < 0, ∀ω ∈ R * .
Then, combining (2.2)-(2.4), we get the sufficient condition (2.1). Thus L 1,p,µp has the spectral properties (P 1)-(P 3). Let λ 1,p,µp be the unique negative (simple) eigenvalue of L 1,p,µp , and let χ 1,p,µp be the, normalized in H 2 (R), eigenfunction associated with λ 1,p,µp . One can notice that ϕ 1,p,µp , χ 1,p,µp ⊂ (Ker L 1,p,µp ) ⊥ for the usual inner product on L 2 (R), since χ 1,p,µp , ϕ ′ 1,p,µp L 2 = λ -1 1,p,µp χ 1,p,µp , L 1,p,µp ϕ ′ 1,p,µp L 2 = 0. Moreover, since ϕ 1,p,µp is an even function, then L 1,p,µp χ 1,p,µp (-•) = λ 1,p,µp χ 1,p,µp (-•), and by uniqueness χ 1,p,µp (-•) = χ 1,p,µp (•). This implies that χ 1,p,µp is an even function.

2.2. Construction of the C 1 branch of solitary waves. We construct now the new solitons in a neighborhood of the explicit solitary wave ϕ 1,p,µp .

Lemma 2.2 (Existence of Solitons ϕ 1,p,µp for c close to 1). There exist δ p > 0 and δp > 0 such that for any c > 0 with |c -1| < δ p , there exists a unique H 4 (R) even solution ϕ c,p,µp of (1.4) in the ball of H 4 (R) centered at ϕ 1,p,µp with radius δp > 0. Moreover, the function c → ϕ c,p,µp is of class

C 1 from ]1 -δ p , 1 + δ p [ into H 4 (R).
Proof. We apply of the Implicit Function Theorem, (see for instance [START_REF] Mariş | Global branches of travelling-waves to a Gross-Pitaevskii-Schrödinger system in one dimension[END_REF] for a similar application). Recall that we set

H s e (R) = {u ∈ H s (R) : u(-•) = u(•)}. For γ > 0, we define the map T :]1 -γ, 1 + γ[×H 4 e (R) → L 2 e (R), by: (c, ψ) → µ p ∂ 4 x ψ -∂ 2 x ψ + cψ -1 p+1 ψ p+1 . T is obviously of class C 1 and, since ϕ 1,p,µp satisfies (1.4)-(1.6), it holds (2.5) T (1, ϕ 1,p,µp ) = 0 and ∂ ψ T (1, ϕ 1,p,µp ) = L 1,p,µp H 4 e (R)
,

where L 1,p,µp

H 4 e (R)
denotes the restriction of L 1,p,µp to H 4 e (R).

It is easy to check that the spectrum of L 1,p,µp

H 4 e (R)
:

H 4 e (R) → L 2 e (R)
H 4 e (R)
is an isomorphism.

Therefore, according to the Implicit Function Theorem, there exist δ p > 0 and a

C 1 map R :]1 - δ p , 1 + δ p [→ H 4 e (R) which is uniquely determined such that: T (c, R(c)) = 0 for all (c, R(c)) ∈]1 -δ p , 1 + δ p [×B H 4 e (ϕ 1,p,µp , δ p )
, where B H 4 e (ϕ 1,p,µp , δ p ) denotes the unit ball of H 4 e (R) centered at ϕ 1,p,µp with radius δ p . This proves the lemma by setting ϕ c,p,µp = R(c) for c ∈]1 -δ p , 1 + δ p [. Remark 2.1 (Some Important Properties of ϕ c,p for c ∼ 1). Since ϕ c,p,µp ∈ H 4 (R), using the equation (1.4) of ϕ c,p,µp , and a classical bootstrap argument, we get that ϕ c,p,µp ∈ H k (R) for all k ∈ N. The well-known Sobolev embedding of

H k+1 (R) into C k (R) leads to (2.6) ϕ c,p,µp C k ([n,n+1]) ≤ C S ϕ c,p,µp H k+1 ([n,n+1]) → 0 as n → +∞,
and thus ∂ k x ϕ c,p,µp (±∞) = 0 for all k ∈ N. Now, for c close to 1, in view of the solutions to the linear asymptotic equation: µ p ∂ 4

x ψ -∂ 2 x ψ + cψ = 0, we infer that ϕ c,p,µp and its derivatives satisfy

(2.7) |∂ k x ϕ 1,p,µp (x)| e - √ c|x| for k ∈ {0, 1, 2, 3, 4}.
2.3. Orbital stability result assuming L -1 1,p,µp ϕ 1,p,µp , ϕ 1,p,µp L 2 < 0. Recall that, according to Lemma 2.1, the operator L 1,p,µp possesses the spectral properties (P 1)-(P 3). Then it is well known that if w ∈ H 2 (R) satisfies L 1,p,µp w, w L 2 < 0 then (1.16) holds for any v ∈ H 2 (R) such that v, ϕ ′ c,p L 2 = v, L 1,p,µp w L 2 = 0. Now, since ϕ 1,p,µp is even, it follows from the last subsection that ϕ 1,p,µp ∈ Im L 1,p,µp and applying the above criterium with w ∈ L -1 1,p,µp ϕ 1,p,µp , we obtain that the hypotheses of Proposition 1.1 are satisfied as soon as (2.8)

J p = L -1 1,p,µp ϕ 1,p,µp , ϕ 1,p,µp L 2 < 0 . Note that the above quantity does not depend on the choice of the element of the preimage of ϕ 1,p,µp since ϕ 1,p,µp ∈ (ker L 1,p,µp ) ⊥ . In the next subsection we check numerically that this sign condition is fulfilled for p ∈ {1, 2, 3, 4}. In this subsection, we prove the orbital stability result assuming (2.8).

According to Proposition 1.1 , it suffices to check that for any c close enough to 1 and

any v ∈ H 2 (R) such that v, ϕ c,p L 2 = v, ϕ ′ c,p L 2 = 0, there exists δ p > 0 such that (2.9) L c,p v, v L 2 ≥ δ p v 2 H 2 . First, we notice that setting γ c,p = ϕ c,p,µp -ϕ 1,p,µp H 1 for |c -1| < δ p , it holds (2.10) v, ϕ 1,p,µp L 2 ≤ γ c,p v L 2 (R) and v, ϕ ′ 1,p,µp L 2 ≤ γ c,p v L 2 (R) .
Next, for all v ∈ H 2 (R) satisfying the almost orthogonality conditions (2.10) with γ c,p small enough, let us prove that L 1,p,µp is coercive in H 2 (R). We argue as, for instance, Côte et al ([6], Lemma 2.6).

(2.10) clearly leads to

(2.11) v, ϕ 1,p,µp ϕ 1,p,µp L 2 (R) L 2 + v, ϕ ′ 1,p,µp ϕ ′ 1,p,µp L 2 (R) L 2 γ c,p v H 2 (R) .
Now, we decompose v as follows

(2.12) v = v 1 + a 1 ϕ 1,p,µp ϕ 1,p,µp L 2 + a 2 ϕ ′ 1,p,µp ϕ ′ 1,p,µp L 2 = v 1 + v 2 , with v 1 , ϕ 1,p,µp L 2 = v 1 , ϕ ′ 1,p,µp L 2
= 0. Then, combining (2.11) and (2.12), we infer that (2.13)

|a 1 | + |a 2 | γ c,p v H 2 (R) .
Moreover, for k = 0, 1, 2, using that

∂ k x ϕ 1,p,µp , ∂ k x ϕ ′ 1,p,µp L 2 = 0, it clearly holds (2.14) ∂ k x v 2 2 L 2 (R) = a 2 1 ∂ k x ϕ 1,p,µp 2 L 2 (R) ϕ 1,p,µp 2 L 2 (R) + a 2 2 ∂ k x ϕ ′ 1,p,µp 2 L 2 (R) ϕ ′ 1,p,µp 2 L 2 (R)
.

Then, combining (2.12)-(2.14), we deduce that there exists 0 < ε p < 1 such that

(2.15) v 1 H 2 (R) ∼ v H 2 (R) as soon as γ c,p < ε p .
Next, we compute

(2.16) L 1,p,µp v, v L 2 = L 1,p,µp v 1 , v 1 L 2 + L 1,p,µp v 2 , v 2 L 2 + 2 L 1,p,µp v 1 , v 2 L 2 .
Note that, by construction, v 1 is orthogonal in L 2 (R) with ϕ 1,p,µp and ϕ ′ 1,p,µp . Therefore, according to the discussion in the beginning of this Subsection, there exists δ 0 > 0 such that (2.17)

L 1,p,µp v 1 , v 1 L 2 ≥ δ 0 v 1 2 H 2 (R) δ 0 v 2 H 2 (R)
. On the other hand, using (2.13)-(2.15), it holds

(2.18) L 1,p,µp v 2 , v 2 L 2 |a 1 | 2 + |a 2 | 2 (γ c,p ) 2 v 2 H 2 (R) and (2.19) L 1,p,µp v 1 , v 2 L 2 v 1 H 2 (R) v 2 H 2 (R) v H 2 (R) (|a 1 | + |a 2 |) γ c,p v 2 H 2 (R) .
Therefore, combining (2.16)-(2.19), we deduce that there exists δ ′ p > 0 and

ε p > 0 such that if ϕ c,p - ϕ cp,p L 2 < ε p then (2.20) L 1,p,µp v, v L 2 ≥ δ ′ p v 2 H 2 (R) .
(2.9) follows immediately by noticing that

L c,p,µp v, v L 2 = L 1,p,µp v, v L 2 + (c -1) v 2 L 2 (R) -v 2 , ϕ p c,p,µp -ϕ p 1,p,µp L 2 ≥ δ ′ p v 2 H 2 (R) -(c -1 + Kε p ) v 2 L 2 (R) ≥ δ ′ p 2 v 2 H 2 (R) , (2.21) as soon as |c -1 + Kε p | < δ ′ p /2.
But the continuity of the branch of the ϕ c,p,µp in H 4 (R) ensures that this is true as soon as c is close enough to 1 and we are done.

Numerical checking of the sign condition on

L -1 1,p,µp ϕ 1,p,µp , ϕ 1,p,µp L 2 . 2.4.1. Numerical computing of J 1 = L -1 1,1,µ1 ϕ 1,1,µ1 , ϕ 1,1,µ1 L 2 . Let ρ 1 ∈ L -1 1,1,µ1 ϕ 1,1
Note that, since L 1,1,µ1 is a self-adjoint operator on L 2 (R), the value of J 1 does not depend on the choice of the element ρ 1 . Moreover, using that

L 1,1,µ1 ϕ 1,1,µ1 = -1 2 ϕ 2 1,1,µ1 , we observe that (2.22) ρ 1 , ϕ 2 1,1,µ1 L 2 = -2 ϕ 1,1,µ1 2
L 2 (R) < 0.
From (2.22) and the exponential decay properties of ϕ 1,1,µ1 and ρ 1 (see (2.7)), we infer that ρ 1 takes necessarily some negative values near the origin. Then, we can expect that J 1 is negative.

The goal is to solve numerically the following equation:

(2.23)

∂ 4 x ρ 1 = 13 2 6 2 ∂ 2 x ρ 1 -ρ 1 + ϕ 1,1,µ1 ρ 1 + ϕ 1,1,µ1 .
First, we remark that ρ 1 (•) and ρ 1 (-•) are both solutions of (2.23), since ϕ 1,1,µ1 is an even function.

Then, from (2.22) we deduce that ρ 1 is also an even function. Moreover, we get that ∂ 2k x ρ 1 is even and ∂ 2k+1

x ρ 1 is odd for all k ∈ N. We can restrict ourself to study equation (2.23) on R + . We fix the domain to [0, r max ], with 0 < r max < +∞, and we rewrite equation (2.23) as a system of four first-order ODE. We set:

ρ 1,1 = y 1 , y ′ 1 = y 2 , y ′ 2 = y 3 , y ′ 3 = y 4 and y ′ 4 = 13 2 6 2 (y 3 -y 1 + ϕ 1,1,µ1 y 1 + ϕ 1,1,µ1
). Next, we choose the Robin type boundary conditions at point r max : y 1 (r max )+y 2 (r max ) = 0 and y 3 (r max )+y 4 (r max ) = 0, and at point 0 we use the symmetry of ρ 1,1 : y 2 (0) = 0 and y 4 (0) = 0. Finally, in our numerical scheme, we take into account the exponential decay properties:

∂ k x ϕ 1,1,µ1 (x) ∼ ∂ k x ρ 1 (x) ∼ e -x
, for x ∼ r max , and for all k ∈ N. Therefore, computing with the MATLAB solver, the result is: ρ 1,1 , ϕ 1,1 L 2 ([0,rmax]) ≈ -10.0787 < 0 (see Fig. 1a-1b).

Numerical values of J

p = L -1 1,p,µp ϕ 1,p,µp , ϕ 1,p,µp L 2
for p ∈ {2, 3, 4, 5}. Proceeding as in Subsection 2.3, Step 1, for p ∈ {2, 3, 4, 5}, we compute respectively that: J 2 ≈ -1.9325 < 0, J 3 ≈ -0.5649 < 0, J 4 ≈ -0.1443 < 0 (see Fig. 2a-2b) and J 5 ≈ 0.0252 > 0 (see Fig. 2c-2d). More precisely, we lose the sufficient condition of stability exactly for p crit = 4.84. We observe that the same critical value appears in the stability analysis done by Bridges et al. [START_REF] Bridges | Stability and instability of solitary waves of the fifth-order KdV equation: a numerical framework[END_REF] (see Section 8, p. 209-210) on the similar model, using the Evans functions approach. Finally, we point out that thanks to the term ∂ 5

x u, the family {ϕ c,p,1 , c ∈]c p -δ p , c p + δ p [} ⊂ H 4 (R), with 0 < δ p ≪ 1, of even solitons of gKW remains stable in the critical gKdV case p = 4.

Existence and stability of slow solitons of gKW

In this section, we prove the point (ii) of Theorem 1.1. The idea is as follows. Around the functional equation E ′ µ (ϕ 1,p,µ ) + V ′ (ϕ c,p,µ ) = 0 (equivalent to equation (1.4)), we will construct a minimization problem on the line. By applying the Concentration-Compactness Principe [START_REF] Lions | The concentration-compactness principle in the calculus of variations. The locally compact case[END_REF], we will prove the existence of even solitary waves ϕ 1,p,µ of gKW, traveling with the speed 1. Next, we will establish the strong convergence in H 1 (R) of the family {ϕ 1,p,µ , 0 < µ ≪ 1} to the solitons ϕ 1,p = ϕ 1,p,0 of gKdV as µ tends to 0 + (see Subsection 3.1). Once again, we recall that the solitons ϕ c,p,0 of gKdV satisfy equation (1.4) with µ = 0, i.e.:

(3.1) -∂ 2 x ϕ c,p,0 (x) + cϕ c,p,0 (x) = 1 p + 1 ϕ p+1 c,p,0 (x), ∀x ∈ R.
They are unique, up to the translations, and the transformations: ϕ c,p → -ϕ c,p if p is even. Moreover, they are explicitly given by:

(3.2) ϕ c,p (x) = (p + 1)(p + 2)c 2 1/p sech 2/p p √ c 2
x , ∀x ∈ R.

Finally, using the well-known spectral properties of the linearized operateur L c,p,0 = -∂ 2 x (•) + c(•)ϕ p c,p,0 of gKdV around the solitons ϕ c,p,0 , and arguing as in Subsection 2.3, we will prove the stability in H 2 (R) of the family {ϕ 1,p,µ , 0 < µ ≪ 1} for the sub-critical gKdV nonlinearity p ∈ {1, 2, 3}. Note that we also obtain a uniqueness result on the solution to our minimization problem. After an appropriate rescaling, this enables us to construct a continuous branch of even solitary waves of gKW traveling with low speeds that are all orbitally stable .

3.1. Existence and limit of even ground states of gKW as µ tends to 0 + . In this subsection, we first prove the existence of even solitary waves ϕ 1,p,µ of gKW for µ > 0, by solving a minimization problem on the line. In this step we argue as in Levandosky [START_REF] Levandosky | A stability analysis of fifth-order water wave models[END_REF] (Theorem 2.3). Second, we study the strong convergence in H 1 (R) of the family {ϕ 1,p,µ , 0 < µ ≪ 1} to the explicit solitons ϕ 1,p of gKdV (given in (3.2)) as µ goes to 0 + (see Lemma 3.1). Moreover, we establish the uniqueness of ϕ 1,p,µ for 0 < µ ≪ 1 and p ∈ {1, 2, 3}(see Lemma 3.2).

Lemma 3.1 (Existence and limit of ϕ 1,p,µ ). Let p ∈ N * be fixed. There exists a family {ϕ 1,p,µ , 0 < µ ≪ 1} ⊂ H 4 (R) of even solutions to (1.4) with c = 1, such that

(3.3) lim µ→0 + ϕ 1,p,µ -ϕ 1,p,0 H 1 (R) = 0.
Proof. We split the proof in five steps.

Step 1. Minimization problem on R.

We argue similarly as Levandosky [START_REF] Levandosky | A stability analysis of fifth-order water wave models[END_REF] (Theorem 2.3). First, we set

H 2 e (R) = ψ ∈ H 2 (R) : ψ(-•) = ψ(•) . Clearly, H 2 e (R)
, endowed with the scalar product of H 2 (R), is a Hilbert space since we have a continuous embedding of H 2 (R) into C 1 (R). For all ψ ∈ H 2 e (R), we define the functional associated with the linear part of (1.4) by:

(3.4) I µ (ψ) = I 1,µ (ψ) = R µ 2 (∂ 2 x ψ) 2 + 1 2 (∂ x ψ) 2 + 1 2 ψ 2 ,
and the functional associated with the nonlinear part of (1.4) by:

(3.5) K p (ψ) = 1 (p + 1)(p + 2) R ψ p+2 .
One can see that I µ (•) is coercive in H 2 (R), since for all ψ ∈ H 2 (R), we have

(3.6) I µ (ψ) ≥ min{µ, 1} 2
ψ 2 H 2 (R) .
Also, one can check that K p (•) is locally Lipschitz in L p+2 (R). Indeed, for all (ψ, φ) ∈ L p+2 (R)×L p+2 (R), applying the Hölder inequality, we have

|K p (ψ) -K p (φ)| = 1 (p + 1)(p + 2) R (ψ p+2 -φ p+2 ) ≤ 1 (p + 1)(p + 2) ψ -φ L p+2 (R) p+1 k=0 |ψ| p-k+1 |φ| k L p+2 p+1 (R) , (3.7)
and applying the Young inequality, we have

(3.8) |ψ| p-k+1 |φ| k p+2 p+1 L p+2 p+1 (R) ≤ p -k + 1 p + 1 ψ p+2 L p+2 (R) + k p + 1 φ p+2 L p+2 (R) .
Then combining (3.7) and (3.8), we obtain

(3.9) |K p (ψ) -K p (φ)| ≤ C( ψ L p+2 (R) , φ L p+2 (R) ) ψ -φ L p+2 (R) .
Now, we will study the minimization of I µ (•) in H 2 e (R) subject to the constraint K p (•) = K p (ϕ 1,p,0 ) = β p where ϕ 1,p,0 is defined in (1.9). We chose this constraint to be sure that the family of minimizer (depending on the parameter µ) will converge to the solitons ϕ 1,p,0 of gKdV (in a sense that we specify in Step 5). For β > 0, we define: We claim that there exist a sub-sequence (ψ j ) j≥1 and a function ψ µ ∈ H 2 e (R), such that ψ j → ψ µ strongly in H 2 (R) as j → +∞. Moreover, ψ µ is a minimizer of the problem (3.10).

(3.10) S β p,µ = S β 1,p,µ = inf I µ (ψ) : ψ ∈ H 2 e (R), K p (ψ) = β . Since ϕ 1,p ∈ H 2 e (R
From the convergence of the energy (3.13), we deduce that there exists ε 0 > 0, such that

(3.14) 0 < I c,µ (ψ k ) < S c,p,µ + ε, with 0 < ε < ε 0 ≪ 1, and this implies that the sequence (ψ k ) k≥1 is bounded in H 2 (R) (since I c,µ (•) is equivalent to H 2 (R))
. Since H is reflexive, there exist a sub-sequence (ψ j ) j≥1 ⊂ H and a function ψ µ ∈ H such that the following hold:

(3.15)

ψ j ⇀ ψ µ weakly in H 2 (R), (3.16) ψ j → ψ µ strongly in H 1 loc (R),
(3.17) ψ j → ψ µ a.e. on R, and

(3.18) ψ µ 2 H 2 (R) ≤ lim inf j→+∞ ψ j 2 H 2 (R) .
Now, we define a sequence of positive and even functions:

(3.19) Φ k = ∂ 2 x ψ k 2 + |ψ k | 2 , ∀k ∈ N * .
From (3.14), one can see that (Φ k ) k≥1 is bounded in L 1 (R). After extracting a sub-sequence, we may assume that lim k→+∞ R Φ k = L < +∞. By normalizing, we may assume further that R Φ k = L for all k ∈ N * . Then, by the Concentration-Compactness Lemma [START_REF] Levandosky | A stability analysis of fifth-order water wave models[END_REF], there are three possibilities: (a) Compactness: there exists (y k ) k≥1 ⊂ R, such that for all ε > 0, there exists R ε > 0, such that for all k ∈ N * , (3.20) (c) Dichotomy: there exists l ∈ (0, L), such that for all ε > 0, there exist R > 0, R k → +∞, (y k ) k≥1 and k 0 , such that for all k ≥ k 0 ,

|x-y k |≤Rε Φ k ≥ L -ε ⇔ |x-y k |≥Rε Φ k ≤ ε.
(3.22) |x-y k |≤R Φ k -l ≤ ε and R<|x-y k |<R k Φ k ≤ ε.
First, let us assume that (a) holds, and we will prove that ψ µ is a minimizer of the problem (3.10). One can remark that |y k | ≤ R L/2 for all k ∈ N * , otherwise, there exists k 0 such that {|x -

y k0 | ≤ R L/2 } ⊂ R + .
Then, using that Φ k0 is even, we get a contradiction:

(3.23) R Φ k0 > 2 |x-y k 0 |<R L/2 Φ k0 > L.
Thus, we can assume that y k = 0 for all k ∈ N * in hypothesis (3.20) by taking as radius R L/2 + R ε . Now, by the compact embedding of H 1 loc (R) into L p+2 loc (R), and using (3.16), we obtain (3.24) ψ j → ψ µ strongly in L p+2 loc (R). On the other hand, using that ψ j L p+2 (R) ≤ 1 2 p/2(p+2) ψ j H 1 (R) , and the compactness of (Φ j ) j≥1 , one can easily check that the sequence (|ψ j | p+2 ) j≥1 is also compact, i.e. This clearly leads to

(3.25) ψ j → ψ µ strongly in L p+2 (R).
Finally, by lower semi-continuity (3.18), we get I µ (ψ µ ) ≤ S βp p,µ , and from the fact that K p (•) is locally Lipschitz in L p+2 (R) (see (3.9)), and using (3.25), we get K p (ψ µ ) = β p . Therefore, ψ µ is a minimizer of the problem (3.10).

Next, we suppose that (b) holds and prove that this leads to a contradiction. One can easy estimate that for all y ∈ R,

(3.26) |x-y|≤1 ψ p+2 k |x-y|≤1 Φ k p+2 2
.

Then, using (b), for any fixed ε > 0, there exists

k ε ≥ 0 such that for k ≥ k ε , (3.27) |x-y|≤1 ψ p+2 k ε p+1 2 .
Now, multiplying (3.27) by 1 (p+1)(p+2) and summing over intervals centered at even integers, we obtain

(3.28) K p (ψ k ) ε p-1 2 → 0 as ε → 0 .
This contradicts the fact that K p (ψ k ) = β p for all k ∈ N * .

At last, we assume that (c) holds, and we will rule out this possibility. We define two smooth non negative test-functions (C ∞ (R)) as follows:

(3.29) ξ 1 (x) = 1, |x| ≤ 1, 0, |x| > 2, and ξ 2 (x) =    1, |x| ≥ 1, 0, |x| ≤ 1 2 .
Next, we construct two sequences

ψ k,1 (x) = ξ 1 (|x -y k | /R)ψ k (x) and ψ k,2 (x) = ξ 2 (|x -y k | /R k )ψ k (x),
for all x ∈ R and k ∈ N * . Using the assumption (c) and the continuous embedding from

H 1 (R) into L p+2 (R), one can check that for k ≥ k 0 , (3.30) I µ (ψ k ) = I µ (ψ k,1 ) + I µ (ψ k,2 ) + O(ε) and (3.31) K p (ψ k ) = K p (ψ k,1 ) + K p (ψ k,2 ) + O(ε). Since (ψ k ) k≥1 is uniformly bounded in H 2 (R) and ψ k,i H 2 (R) ≤ ψ k H 2 (R) , for i = 1, 2, it follows that (ψ k,i ) k≥1 , for i = 1, 2
, are also bounded uniformly in H 2 (R). Hence K p (ψ k,i ), for i = 1, 2, are bounded and we may pass to a sub-sequence to define β i (ε) = lim k→+∞ K p (ψ k,i ), for i = 1, 2. Since β i (ε), for i = 1, 2, are bounded uniformly, we can choose a sequence ε j → 0 such that β p,i = lim j→+∞ β i (ε j ) < ∞, for i = 1, 2. We clearly get β p,1 + β p,2 = β p , and there are three cases to consider. Case 1: if β p,1 ∈ (0, β p ). Using (3.11) and (3.30), we compute that

I µ (ψ k ) ≥ S Kp(ψ k, 1 ) p,µ + S Kp(ψ k,2 ) p,µ
+ O(ε j ) = K p (ψ k,1 ) β p 2 p+2 + K p (ψ k,2 ) β p 2 p+2 S βp p,µ + O(ε j ). (3.32)
Letting k → +∞, and using that (ψ k ) k≥1 is a minimizing sequence, we obtain

(3.33) S βp p,µ ≥ β p,1 β p 2 p+2 + β p,2 β p 2 p+2
S βp p,µ + O(ε j ).

Letting j → +∞, we get a contradiction:

(3.34) S βp p,µ ≥ β p,1 β p 2 p+2 + β p,2 β p 2 p+2
S βp p,µ > S βp p,µ .

Case 2: if β 1 = 0 (or equivalently β 1 = β c,p ). Using the coercivity of I µ (•) in H 2 (R), and the assumption (c), we have

I µ (ψ k,1 ) = R ∂ 2 x ψ k,1 2 + |ψ k,1 | 2 = |x-y k |≤2R ∂ 2 x ψ k 2 + |ψ k | 2 + O(ε j ) = l + O(ε j ). ( 3
I µ (ψ k ) ≥ l + K p (ψ k,2 ) β p 2 p+2 S βp p,µ + O(ε j ).
Letting k, j → +∞, we get a contradiction:

(3.37) S βp p,µ ≥ l + S βp p,µ > S βp p,µ . Case 3: if β 1 > β c,p (or equivalently β 1 < 0). Using the non negativity of I µ (•), we have

(3.38) I µ (ψ k ) ≥ I µ (ψ k,1 ) + O(ε j ) ≥ K p (ψ k,1 ) β p 2 p+2 S βp p,µ + O(ε j ).
As previously, letting k, j → +∞, once again, we get a contradiction:

(3.39) S βp p,µ ≥ β 1 β p 2 p+2
S βp p,µ > S βp p,µ .

Therefore, we deduce that the sequence (ψ k ) k≥1 is compact.

Step 2. Euler-Lagrange equation.

At this stage, we can write the Euler-Lagrange equation related to the minimization problem (3.10). Since ψ → I p (ψ) and ψ → K p (ψ) are obviously of class C 1 in H 2 (R), the minimizer ψ p,µ satisfies the following functional equation: For 0 < µ < 1 4 , let us set:

(3.40) µ ∂ 2 x ψ p,µ , ∂ 2 x φ L 2 + ∂ x ψ p,µ , ∂ x φ L 2 + ψ p,µ , φ L 2 = α p,µ p + 1 ψ p+1 p,µ , φ L 2 , ∀φ ∈ H 2 e (R
(3.42)

s 1 = ± 1 - √ 1 -4µ √ 2µ and s 2 = ± 1 + √ 1 -4µ √ 2µ .
One can check that e s1x and e s2x are solutions of the linear asymptotic equation: µ∂ 4 x ψ -∂ 2 x ψ + ψ = 0. Then, we infer the following exponential decay properties: Step 3. Refined estimate on Lagrange multipliers. We will prove that the family {α p,µ , 0 < µ < 1} of Lagrange multipliers is bounded. Multiplying (3.41) by ψ p,µ , and integrating on space, we obtain

(3.44) α p,µ = 2I µ (ψ p,µ ) (p + 2)β p = 2S βp p,µ (p + 2)β p ≤ 2I µ (ϕ 1,p,0 ) (p + 2)β p
where ϕ 1,p,0 is the soliton of the gKdV equation given by (1.9). From the uniqueness of ϕ 1,p,0 (up to the symmetries of the equation) and the definition of β p , we know that ϕ 1,p is the unique solution to the constraint minimization problem S This last estimate together with (3.44) clearly ensure that the family {ψ p,µ , 0 < µ ≪ 1} is bounded in H 1 (R).

Step 4. H 2 -boundedness of the family {ψ p,µ , 0 < µ < 1} . We will prove that the family {ψ p,µ , 0 < µ < 1} is bounded in H 2 (R). Applying the Fourier transformation on equation (3.41), we obtain for all ξ ∈ R,

(3.47) 1 + |ξ| 2 + µ |ξ| 4 ψ p,µ (ξ) = α p,µ p + 1 ψ p+1 µ (ξ).
Since {α p,µ , 0 < µ < 1} is bounded (from Step 3), it suffices to prove that the family { ψ p+1 p,µ , 0 < µ < 1} is bounded in L 2 (R), and by the Plancherel-Parseval identity, this is equivalent to prove that the family {ψ p,µ , 0 < µ < 1} is bounded in L 2p+2 (R). Thus, we compute

ψ p,µ 2p+2 L 2p+2 (R) ≤ ψ p,µ 2p L ∞ (R) ψ p,µ 2 L 2 (R) ≤ 2 -p ψ p,µ 2p+2 H 1 (R) , (3.48)
where we use the Sobolev embedding of H 1 (R) into C 0 (R) (with the best constant). Then, using that {ψ p,µ , 0 < µ < 1} is bounded in H 1 (R) (from Step 1), and (3.48), we obtain the desired result. Note that this boundedness, together with the fact that ϕ 1,p,0 is a solution to the constraint minimization problem S Step 5. Strong convergence in H 1 (R) as µ → 0 and construction of the family of solitary waves {ϕ 1,p,µ , 0 < µ < 1}.

Let (µ n ) n≥1 ⊂ R * + be a sequence which decreases to 0 + , and for simplicity, let us rename ψ p,µn = ψ n and α p,µn = α p,n . Then, there exist a sub-sequence, still denoted by (ψ n ) n≥1 , that converges weakly in H 2 e (R) towards a function ψ ∞ ∈ H 2 e (R). Passing to the limit in (3.41) by making use of (3.50), we get that ψ ∞ satisfies (1.13) with c = 1. Note that, to pass to the limit on the nonlinear part of equation (3.41) we deal as follows: using the pointwise convergence (3.17 Then, applying Lebesgue's Dominated Convergence Theorem, it holds for all φ ∈ H 2 (R),

(3.52) lim n→+∞ ψ p+1 n , φ L 2 = ψ p+1 ∞ , φ L 2 .
Now, the uniqueness result (up to symmetries) of the solutions to (1.13) ensures that ψ ∞ = ϕ 1,p,0 whenever p is even and ψ ∞ = ∓ϕ 1,p,0 whenever p is odd. Since, according to (3.50),

lim n→+∞ I 0 (ψ n ) = I 0 (∓ϕ 1,p,0 ) = I 0 (ψ ∞ ),
it follows that (ψ n ) n≥0 converges actually strongly to ψ ∞ in H 1 (R) and thus K p (ψ ∞ ) = β p > 0. This ensures that in all cases, ψ ∞ = ϕ 1,p,0 and thus

(3.53) ψ µ → ϕ 1,n,0 in H 1 (R) as µ → 0 .
Finally, setting ϕ 1,p,µ = α 1/p p,µ ψ p,µ , ∀n ∈ N * . It is easy to check that for any p ∈ N * and any µ > 0, ϕ 1,p,µ is a solution of (1.4) with c = 1. Moreover, according to Step 4., (3.50) and (3.53), the family {ϕ 1,p,µ , 0 < µ < 1} is bounded in H 2 (R) and

(3.54) ϕ 1,p,µ → ϕ 1,p,0 in H 1 (R) as µ → 0 .
This completes the proof of the lemma.

3.2. Uniqueness of even ground states. Note that the construction of the ϕ 1,p,µ ensures that they are even ground state solutions (in the sense of Definition 1.2) to (1.13) with c = 1. The following lemma ensures that for µ small enough, the even ground state solution to (1.13), with c = 1, is unique.

Lemma 3.2 (Uniqueness of ϕ 1,p,µ for 0 < µ ≪ 1). Let p ∈ {1, 2, 3} be fixed and set c = 1. There exists δ ′ p > 0 such that for all µ ∈]0, δ ′ p [, there exists a unique even ground state solution ϕ 1,p,µ to equation (1.4). Moreover, the map µ → ϕ 1,p,µ is-continuous from ]0, δ ′ p [ into H 1 (R). Proof. We fix p ∈ {1, 2, 3}. We first claim that for any ε > 0 there exists µ ε > 0 such that for any 0 < µ < µ ε , any even ground state solution ϕ to (1.4) with c = 1 satisfies

(3.55) ϕ -ϕ 1,p,0 H 1 < ε .
To prove this claim we proceed by contradiction, assuming that there exist ε 0 > 0, a sequence (µ n ) n≥0 of positive real numbers that converges to 0 and a sequence (ϕ n ) n≥0 of ground state solutions to (1.4) with c = 1 and µ = µ n such that

(3.56) ϕ n -ϕ 1,p,0 H 1 ≥ ε 0 , ∀n ∈ N .
From (1.4) we infer that K p (ϕ n ) = 2 p+2 I µn (ϕ n ) > 0. We set

ψ n = K p (ϕ 1,p,0 ) K p (ϕ n ) 1 p+2 ϕ n so that K p (ψ n ) = K p (ϕ 1,p,0 ) = β p . (ψ n ) n≥0
is thus a sequence of solutions to S βp p,µn and as in the proof of Lemma 3.1 (Step 5.), it follows that (ψ) n≥0 is bounded in H 2 (R) and ψ n → ϕ 1,p,0 in H 1 (R). Therefore,

ϕ n = 2I µ (ψ n ) (p + 2)β p 1/p ψ n → ϕ 1,p,0 in H 1 (R)
which contradicts (3.56). Let ϕ 1,p,µ and φ1,p,µ be two ground ground states of (1.4) with c = 1 and µ > 0. Following the idea of Kenig et al. [START_REF] Kenig | Local well-posedness and blow-up in the energy space for a class of L 2 critical dispersion generalized Benjamin-Ono equations[END_REF] (Proposition 3), we will prove, arguing by contradiction, that w = ϕ 1,p,µ -φ1,p,µ = 0 as soon as µ is small enough. Note that, on account of (3.55), it holds w H 1 ≤ ε(µ) with ε(y) → 0 as y → 0. We set F (x) = 1 p+1 x p+1 so that w satisfies (3.57) µ∂ 4

x w -∂ 2 x w + w = F (ϕ 1,p,µ ) -F ( φ1,p,µ ) Denoting by L 1,p,µ the operator defined in Proposition 1.1 with φ c,p,µ = ϕ 1,p,µ , it holds L 1,p,µ = µ∂ 4

x -∂ 2 x + 1 -F ′ (ϕ 1,p,µ ) and thus L 1,p,µ w = F (ϕ 1,p,µ ) -F ( φ1,p,µ ) + F ′ (ϕ 1,p,µ )(ϕ for all φ ∈ H 1 (R) that satisfies the orthogonality conditions : φ, ϕ ′ 1,p,0 = φ, ϕ 1,p,0 = 0. Now, since w is even and ϕ ′ 1,p,µ is odd, we have the first orthogonality condition: (3.60)

w n , ϕ ′ c,p L 2 = 0. Moreover, using that L 1,p,0 ∂ c|c=1 ϕ c,p,0 = ϕ 1,p,0 we get w, ϕ 1,p,0 L 2 = w, L 1,p,0 ∂ c|c=1 ϕ c,p,0 L 2 = L 1,p,0 w, ∂ c|c=1 ϕ c,p,0 L 2 ≤ L 1,p,µ w, ∂ c|c=1 ϕ c,p,0 L 2 + w∂ c|c=1 ϕ c,p,0 , ϕ p c,p,µ -ϕ p 1,p,0 L 2 + µ ∂ 2 x w, ∂ 2

x ∂ c|c=1 ϕ c,p,0 L 2 (ε(µ) + µ) w H 1 (3.61) so that for µ small enough, w is almost orthogonal to ϕ 1,p,0 . Therefore, proceeding as in Subsection 2.3, we infer that for µ > 0 small enough,

L 1,p,0 w, w ≥ α p 2 w 2 H 1 = α p 2 ,
which contradicts (3.59) and proves the uniqueness result for 0 < µ < δ ′ p . Finally, once the uniqueness is proved, the convergence of ϕ 1,p,µ towards ϕ 1,p,µ0 as µ → µ 0 can be proved exactly as the convergence of ϕ 1,p,µ towards ϕ 1,p,0 as µ → 0. Now, we notice that v satisfies (1.4) with p ∈ N * , c = 1 and µ > 0 if and only v µ = µ 1/p v( √ µ •) satisfies

(1.4) with p ∈ N * , c = µ and µ = 1. The following uniqueness result on the ground states to (1.4) with µ = 1 is thus a direct consequence of Lemma 3.2.

3.3.

End of the proof of Theorem 1.1. Let p ∈ {1, 2, 3}. It remains to prove the orbital stability result. Since u is a solution to (1.1) if and only if u µ (t, x) = µ 1/p u(µ 3/2 t, √ µx) is a solution to (1.1) with µ = 1, the orbital stability of ϕ µ,p,1 is equivalent to the orbital stability of ϕ c,p,µ . To prove the orbital stability of this last ground state for µ > 0 small enough, we rely on a continuity argument as in the preceding section. First, by the H 1 -convergence result (3.3), it is easy to check that for µ > 0 small enough, any v ∈ H 2 (R) satisfying v, ϕ 1,p,µ = v, ϕ ′ 1,p,µ = 0 is almost orthogonal in L 2 (R) to ϕ 1,p,0 and ϕ ′ 1,p,0 . Second, recalling that, for p ∈ {1, 2, 3}, d dc |c=1 ϕ c,p,0 2 L 2 < 0 and arguing exactly as in Subsection 2.3, we obtain that L 1,p,0 is coercive in H 1 (R) under these almost orthogonality conditions. It follows that

L 1,p,µ v, v L 2 = µ ∂ 2 x v 2 L 2 (R) + L 1,p v, v L 2 -v 2 , ϕ p 1,p,µ -ϕ p 1,p,0 L 2 ≥ µ ∂ 2 x v 2 L 2 (R) + α p 2 v 2 H 1 (R) -ε(µ) v 2 L 2 (R) ≥ min{µ, α p 2 -ε(µ)} v 2 H 2 (R) , (3.62)
with α p > 0 and ε(y) → 0 as y → 0. By applying Proposition 1.1, we thus get that there exists δ p > 0 such that for all µ ∈]0, δ p [, ϕ 1,p,µ is orbitally stable. This completes the proof of the theorem.

  is contained in the spectrum of L 1,p,µp and Lemma 2.1 (properties (P1) and (P3)) ensures that σ ess (L 1,p,µp ) = [1, +∞[ and Ker L 1,p,µp = span {ϕ ′ 1,p,µp }. Since ϕ ′ 1,p,µp is an odd function, we thus infer that 0 / ∈ σ(L 1,p,µp | H 4 e (R) ) and thus L 1,p,µp

  ϕ 1,1,µ 1 (x) and ρ 1 (x) profiles.

∂ k x ρ 1

 1 ∂ k x ϕ 1,1 (x), k = 0, . . . , 4, profiles. (x), k = 0, . . . , 4, profiles.

Figure 1 .

 1 Figure 1. Numerical computation of J 1

  ϕ 1,4,µ 4 (x) and ρ 4 ∈ L -1 1,4,µ 4 ϕ 1,4 (x) profiles. µ 4 ϕ 1,4,µ 4 (y))ϕ 1,4,µ 4 (y) profile. ϕ 1,5,µ 5 (x) and ρ 5 = L -1 1,5,µ 5 ϕ 1,5,µ 5 (x) profiles.

  µ 5 ϕ 1,5,µ 5 (y))ϕ 1,5,µ 5 (y) profile.

Figure 2 .

 2 Figure 2. Numerical computation of J p for p = 4, 5.

  ), and I µ (•) is coercive in H 2 e (R), one can see that S βp p,µ > 0. On the other hand, using the definition of I µ (•) and K p (•), one can check that (, ∀β ∈ R * + , and then we infer the crucial sub-additivity property: (3.12) S β p,µ + S βp-β p,µ > S βp p,µ , ∀β ∈ (0, β p ). Let us solve the constraint minimization problem (3.10) with β = β p . We take (ψ k ) k≥1 ⊂ H a minimizing sequence of the problem (3.10), i.e., for all k ∈ N * , we have (3.13) K p (ψ k ) = β p and lim k→+∞ |I µ (ψ k ) -S p,µ | = 0.

  |ψ j | p+2 = 0 .

  .35) Thus, combining (3.30) and (3.35), we get (3.36)

( 3 .

 3 43) |∂ sx ψ p,µ (x)| e -|x| , for |x| ≫ 1, s ∈ N, and 0 < µ ≪ 1.

  α p,µ ≤ 1 + O(µ) .

I

  µ (ψ p,µ ) = lim inf µ→0 I 0 (ψ p,µ ) ≥ I 0 (ϕ 1,p,0 ) . Gathering this last estimate with (3.44)-(3.46), we get (3.50) lim µ→0 S βp p,µ = lim µ→0 I µ (ψ p,µ ) = lim µ→0 I 0 (ψ p,µ ) = I 0 (ϕ 1,p,0 ) and lim µ→0 α p,µ = 1 .

  ), the H 2 -boundedness and the continuity of the function: y → y p+1 , one can check that (3.51) ψ p+1 n → ψ p+1 ∞ a.e. on R.

  ), where α p,µ ∈ R * is a Lagrange multiplier. It is worth noticing that, equation (3.40) also holds with any odd test-function φ ∈ H 2 (R), since then both members do cancel. Therefore ψ p,µ satisfies equation (3.40) for all φ ∈ H 2 (R). Now, by a standard bootstrap argument, we get that ψ p,µ ∈ H s (R) and ∂ sx ψ p,µ (±∞) = 0 for all s ∈ N. Therefore ψ p,µ is a strong solution of (3.40), i.e.:

	(3.41)	µ∂ 4 x ψ p,µ (x) -∂ 2 x ψ p,µ (x) + ψ p,µ (x) =	α p,µ p + 1	ψ p+1 µ (x), ∀x ∈ R.

  1,p,µ -φ1,p,µ ) = F ′′ (W p,µ ) 2 w2 , with, for all x ∈ R, W p,µ (x) ∈ [ϕ 1,p,µ (x), φ1,p,µ (x)]. Therefore, assuming that w = 0 and setting w = w w H 1 , we get(3.58) L 1,p,µ w L 2 (R) 1 w H 1 w L 2 w H 1 w L 2 w H 1 so that L 1,p,0 w, w L 2 = L 1,p,µ w, w L 2 + w 2 , ϕ p 1,p,µ -ϕ p 1,p,0 L 2 -µ ∂ 2 x wOn the other hand, it is well-known that for p ∈ {1, 2, 3}, d dc |c=1 ϕ c,p,0 2 L 2 < 0 which ensures that there exists α p > 0 such that L 1,p,0 φ, φ ≥ α p φ 2

		2
		L 2 (R)
	(3.59)	ε(µ) + µ .
		H 1


Bernard P Serpette

David Janin

Causalité dans les calculs d'événements

Si l'on considère un événement comme une valeur quelconque associée à une date de réception, un calcul d'événements est une fonction produisant un ensemble d'événements à partir d'un ensemble d'événements reçus. Un sous-ensemble particulièrement intéressant de ces fonctions correspond aux fonctions causales dont les événements produits avant une date donnée ne dépendent que des événements reçus avant cette même date. Dans cet article, nous proposons une définition très simple permettant de caractériser ces fonctions causales.

1. un sol pour un do, un la pour un ré, un si pour un mi, etc. . .

Introduction

Nous étudions ici la sémantique des programmes lisant, transformant et produisant des structures temporelles. L'élément essentiel d'une structure temporelle est l'événement : une valeur quelconque, spécifique au domaine que l'on manipule, à laquelle on attribue une date : quelque chose se passe à cette date. Une structure temporelle est alors simplement définie comme un ensemble d'événements.

La musique offre de nombreux exemples de structures temporelles. Les événements sont, par exemple, des débuts ou des fins de note. L'informatique musicale est donc largement utilisatrice de programmes manipulant des structures temporelles. On distingue cependant les transformations « hors temps », applicables, par exemple, lors de la composition d'une oeuvre, des transformations « dans le temps », appliquées, à la volée, pendant l'exécution d'une oeuvre. A l'évidence, les secondes forment une sous-classe des premières.

Prenons, par exemple, le programme réduisant par deux, dans une pièce musicale, la durée de chaque note de musique, sans pour autant déplacer le début de ces notes. Les blanches deviennent des noires, les noires des croches, etc. . .. Ce programme est applicable à une partition. La durée de chaque note étant connue, on peut la diviser par deux. Par contre, cette transformation n'est pas applicable à la volée ou autrement dit « dans le temps ». La durée d'une note n'est connue qu'à la fin de son écoute. Lorsque sa fin arrive, l'événement lié à la fin de la note correspondante aurait déjà dû être émis, cet événement fait déjà partie du passé. On dira que ce programme n'est pas causal.

La propriété de causalité (ou de non causalité) n'est cependant pas toujours si évidente à détecter. Prenons, par exemple, le programme qui ajoute la quinte 1 de toute note jouée, deux unités de temps après le début de cette note et en la jouant pendant la même durée, sous la condition qu'aucune autre note n'a été jouée pendant ces deux unités de temps. Bien que cette transformation soit spécifiée à l'aide de dépendances temporelles plutôt complexes entre les événements d'entrée et les événements de sortie, elle n'en est pas moins causale.

Autrement dit, un programme peut réaliser des calculs complexes sur les intervalles de temps écoulés entre deux événements, comme illustré par l'exemple des quintes retardées, tout en restant causal, c'est à dire exécutable à la volée, dans le temps. Pour un langage de programmation dédié à la conception de systèmes musicaux interactifs, tel que le T-calcul [START_REF] Archipoff | Structured reactive programming with polymorphic temporal tiles[END_REF], l'étude et l'analyse de cette notion de fonctions causales sont importantes. Elles permettraient, par exemple, via une analyse statique appropriée, de garantir que les programmes réalisés appartiennent (à une sous-classe décidable) des programmes effectivement exécutables dans le temps. Pour une étude plus approfondie on se réfèrera à [START_REF] Matsikoudis | The fixed-point theory of strictly causal functions[END_REF].

Dans cet article nous donnerons une définition très simple de la causalité et nous prouverons que cette définition caractérise exactement les programmes pouvant émettre leurs sorties soit en réaction aux événements d'entrées, soit en réaction au temps qui passe.

Notations

Dans toute la suite, nil représente la liste vide. Nous utiliserons la notation x • l pour ajouter un élément x dans une liste l, et l 1 • l 2 pour effectuer la concaténation de deux listes.

Temps, événements et chronologies

Le temps est central dans l'étude que nous allons aborder. La propriété la plus importante n'est pas que le temps soit mesurable, via une valeur numérique, mais qu'il soit comparable ; que l'on puisse parler d'avant et d'après une date donnée. Nous nous plaçons donc dans un domaine de temps T muni d'une relation d'ordre totale.

Un événement est une valeur datée, c'est à dire une valeur quelconque estampillée par une date. Si l'on suppose un domaine de valeur V , un événement est un élément de T × V . Nous noterons donc (t, v) un événement construit sur une date t ∈ T et une valeur v ∈ V .

Si l'on considère plusieurs événements, le domaine le plus naturel à considérer est celui des ensembles d'événements, c'est à dire le domaine P(T × V ). Une implémentation des multi-ensembles, quand les éléments possèdent un ordre total, peut être faite avec des listes triées. Cette implémentation est particulièrement agréable car l'identité ensembliste correspond à l'identité structurelle sur les listes. La différence ensemble versus multi-ensemble n'est pas ici primordiale, la vision multi-ensemble est plus générale en considérant la possibilité de la présence simultanée d'un même événement.

Voir un ensemble d'événements comme une liste triée selon la date des événements est assez naturel. Nous appellerons chronologie une telle structure. Une chronologie est donc une liste d'événements ordonnée par ordre décroissant de temps. Implicitement, la variable c est une chronologie, élément de l'ensemble C suivant :

Définition 1 (Chronologie) C = {l : List(T × V ) | sorted(l)}
Le prédicat sorted caractérise les listes triées par ordre décroissant de temps. La figure 1 L'intérêt que le temps soit comparable est, pour un certain temps t, de pouvoir couper une chronologie en deux, en séparant l'après et l'avant t. Dans une chronologie c, les événements antérieurs à un certain temps t, seront notés c ↓ t, et sont définis par :

Définition 2 [Antériorité dans une chronologie] c ↓ t = match c with | nil ⇒ nil | (t , v) • c ⇒ if t ≥ t then c ↓ t else c
Dans l'intention, et en supposant que les chronologies soient représentées par des ensembles, c ↓ t représente l'ensemble {(t , x) ∈ c | t < t}. On remarquera que la définition 2 prend en compte le fait qu'une chronologie est une liste triée. De plus, on observe qu'un événement au temps t n'est pas considéré antérieur à t.

De manière similaire à l'antériorité, les événements postérieurs à un certain temps t, dans la chronologie c, seront notés c ↑ t, et sont définis par :

Définition 3 [postériorité dans une chronologie] c ↑ t = match c with | nil ⇒ nil | (t , v) • c ⇒ if t ≥ t then (t , v) • c ↑ t else nil
La figure 2 . Une fonction sur les chronologies est donc un programme créant un morceau de musique à partir d'un autre morceau de musique. Ceci sans considération des meilleures manières de produire des données temporelles, ces fonctions se divisent en deux catégories, celles qui sont causales et celles qui ne le sont pas. Dans notre cadre, l'intuition de la causalité est : "l'avenir ne peut pas influencer le passé". Autrement dit une fonction est causale, si, pour tout temps t, elle peut produire les événements antérieurs à t sans considérer les événements postérieurs à t. On voit que, dans ces deux définitions informelles, la causalité est définie en mentionnant, à la fois, le passé et le futur. L'originalité de la définition que nous proposons réside dans le fait qu'elle ne fait mention que du passé. Une fonction est causale si, pour deux chronologies partageant les mêmes événements avant un certain temps t, alors les images de ces deux chronologies par cette fonction doivent aussi avoir les mêmes événements avec ce temps t.

Définition 4 (fonction causale) ∀f, causal(f ) ⇔ ∀c 1 c 2 t, c 1 ↓ t = c 2 ↓ t ⇒ f (c 1 ) ↓ t = f (c 2 ) ↓ t
L'intuition dérivant de cette définition est : "ce qui est fait est fait". Et derrière ces mots, si on considère que le passé est révolu pour une fonction causale, alors il y a sûrement une chance de calculer cette fonction au présent, ou, pour le moins, de façon incrémentale. Pour reprendre l'exemple de la musique et plus particulièrement de l'acoustique : la réverbération, l'écho et les filtres, entre autres, sont des fonctions qui transforment un signal sonore, et donc une chronologie si le signal est échantillonné. Ces transformations d'un signal sont causales car elles peuvent être directement appliquées sur le flux d'entrée. Par contre, la transformation qui consiste à inverser, dans le temps, un signal d'entrée, n'est pas causale. Cette transformation doit attendre la fin du flux d'entrée pour commencer.

Certaines fonctions ne sont pas immédiates à cerner concernant leur causalité. Prenons, par exemple, la fonction qui duplique les événements de son entrée d'un certain quantum de temps, disons ∆, si et seulement si aucun autre événement n'est intervenu dans l'intervalle. Pour la chronologie [9 6 5 4 3 0], en omettant de mentionner les valeurs associées aux événements, et en prenant ∆ = 2, les événements 9, 6 et 0 seront dupliqués et le résultat de la transformation est [11 9 8 6 5 4 3 2 0]. Bien que regardant dans le futur pour prendre une décision, cette fonction est causale.

Etant donné f une fonction sur les chronologies, pour exécuter f de façon incrémentale sur une chronologie c, il faut pouvoir prendre les événements de c par ordre croissant, donc dans l'ordre où ils arrivent dans le temps, les transmettre à f et considérer, événement après événement, que les résultats sont révolus. Si on note [[f ]], l'interprétation incrémentale de f , cette interprétation peut être définie de façon récursive (non terminale) de la manière suivante :

Définition 5 (Calcul incrémental) [[f ]] = λc. match c with | nil| ⇒ f (nil) | (t, v) • c ⇒ (f (c) ↑ t) • ([[f ]](c ) ↓ t)
Autrement dit, pour une chronologie dont le dernier élément dans le temps (tête de la chronologie) se situe au moment t, le comportement de [[f ]] est le même que f après ce temps. Le comportement de [[f ]] avant ce temps pour cette chronologie est le comportement de [[f ]] en considérant que ce dernier événement n'a pas lieu, et oubliant tous les événements produits qui sont postérieurs à t. On remarquera que la concaténation des deux listes, correspondant à l'union ensembliste, est bien fondée sachant que celle de gauche est triée et supérieure à t et que la seconde est triée et inférieure à t.

Calculons symboliquement [[f ]](c) pour la chronologie c = [D C B A].s Si, pour tout événement X, t X est le temps associé à X, on a [[f ]](c) = (f (c) ↑ t D ) • ([[f ]]([C B A]) ↓ t D ). Soit, en dépliant encore une fois, [[f ]](c) = (f (c) ↑ t D ) • ((f ([C B A]) ↑ t C ) • ([[f ]]([B A]) ↓ t C )) ↓ t D . Si l'on continue à déplier et en simplifiant, on arrive à [[f ]](c) = (f (c) ↑ t D ) • (f ([C B A]) ↑ t C ↓ t D ) • (f ([B A]) ↑ t B ↓ t C ) • (f ([A]) ↑ t A ↓ t B )•f (nil) ↓ t A . On voit que [[f ]](c) se calcule avec des tronçons de la forme f (X) ↑ t 1 X ↓ t 2 X où t 1 X
est le plus grand temps de la chronologie X et t 2 X le temps du prochain événement qui adviendra après X. La figure 4 montre ce calcul par tronçons. Les lignes successives montrent les diverses applications de f (X) pour X variant de nil à c. Sur chaque ligne, nous avons partiellement rempli les événements inclus dans le tronçon f (X) ↑ t 1 X ↓ t 2 X . La dernière ligne, correspondant à [[f ]](c), est l'union des tronçons préalablement calculés.

c=[D C B A] f(c) c [[f]](c) A B C D f(nil) f([C B A]) f([B A]) f([B A])
f

([A])

f ((t, v) • c) = f ((t, v) • c) ↑ t • f ((t, v) • c) ↓ t sachant ∀c t, c = c ↑ t • c ↓ t = f ((t, v) • c) ↑ t • f (c) ↓ t par hyp. de causalité et sachant ∀t v c, f ((t, v) • c) ↓ t = c ↓ t = f ((t, v) • c) ↑ t • [[f ]](c) ↓ t par hyp. de rec. = [[f ]](((t, v) • c) par définition. Pour le sens ⇐, sachant c 1 ↓ t = c 2 ↓ t, on a : [[f ]](c 1 ↓ t) ↓ t = [[f ]](c 2 ↓ t) ↓ t ⇔ [[f ]](c 1 ) ↓ t = [[f ]](c 2 ) ↓ t en utilisant le lemme qui suivra. ⇔ f (c 1 ) ↓ t = f (c 2 ) ↓ t par hyp. Lemme 1 ∀f c t, [[f ]](c ↓ t) ↓ t = [[f ]](c) ↓ t

Version opérationelle

Maintenant qu'est établi le bien-fondé de [[f ]] pour les fonctions causales, nous pouvons en donner une version équivalente, totalement réactive :

play(f, c) = let output = f (c) ↑ t 1 c in let event = emitU ntilEvent(output) in play(f, event • c) A chaque

Conclusion

Dans le cadre des fonctions manipulant des données temporelles, nous avons caractérisé celles qui sont causales. Cela nous a amené à prouver que ces fonctions pouvaient se simuler avec une autre fonction qui est réactive. Nous tenons à souligner ici que le fait qu'on peut simuler une fonction causale par une fonction réactive, n'implique pas qu'il faille le faire de cette manière. En effet, la simulation que nous avons donnée impose de conserver l'ensemble des événements déjà rencontrés (effet cumulatif dans la fonction play) 2 .

Le but de notre projet, au-delà de la caractérisation des fonctions causales, est de trouver la meilleure manière de décrire les fonctions manipulant des données temporelles [START_REF] Hudak | Tiled polymorphic temporal media[END_REF] tout en garantissant une certaine efficacité quant à leurs exécutions. Une approche plus théorique du travail proposé dans cet article peut se trouver dans [START_REF] Janin | Timed Denotational Semantics for Causal Functions over Timed Streams[END_REF].

Les preuves ont été effectuées en utilisant le système Coq [START_REF]The coq proof assistant[END_REF]. Ces preuves sont sans difficultés majeures. Au contraire, il se trouve que l'ensemble des lemmes et le théorème final se prouvent sans l'hypothèse que les chronologies soient triées par ordre décroissant de temps. Il y a peut-être matière à réflexion quant à ce résultat. Les sources sont disponibles sur ftp://ftp-sop.inria.fr/indes/rp/jfla2017.v.

Les auteurs remercient Christine Huet pour sa relecture attentive ainsi que les rapporteurs de l'article pour leurs remarques constructives.

Figure 1 -

 1 Figure 1 -Exemple d'une chronologie

  marque avec des triangles pleins l'antériorité d'une chronologie c = [D C B A] pour le temps t, et avec des triangles vides sa postériorité. On a c ↑ t = [D C] et c ↓ t = [B A].

Figure 2 -

 2 Figure 2 -Exemple d'antériorité et de posteriorité

Figure 3 -

 3 Figure 3 -Exemple d'une fonction sur des chronologies

Figure 4 -

 4 Figure 4 -Exemple de la fonction [[]]

Preuve.

  Par induction sur c. Le cas de base est trivial et correspond à f (nil) ↓ t = f (nil) ↓ t. Pour l'induction sachant que l'on ajoute l'événement (t 0 , v) à c, il faut considérer les deux cas t 0 ≥ t et t 0 < t. Le premier cas utilise l'hypothèse de récurrence et le second cas n'utilise que des lemmes simples sur ↑ et ↓. Nous invitons à regarder les sources Coq pour plus de détails.

  appel récursif, la fonction play émet les événements du tronçon f (c) ↑ t 1 c ↓ t 2 c . La variable output contient la partie droite du tronçon, et la fonction EmitU ntilEvent va émettre, dans le temps, les événements de ce tronçon, tant qu'un nouvel événement n'est pas apparu dans le flux d'entrée, présupposé au temps t 2 c . Cette fonction EmitU ntilEvent retourne ce nouvel événement. Pour prendre en compte la dissymétrie du premier tronçon, le premier appel à la fonction play doit se faire avec la chronologie emitU ntilEvent(f (nil)) • nil.

C'est le principal reproche qui a été exprimé envers Fran[START_REF] Elliott | Functional reactive animation[END_REF]


Major in the East India Company's Army. His mother, Ann Sarah Pritchard, died in 1845, and Edward Mockler sen. married Julia Ferryman in 1851. One of their children was Augustus Ferryman Mockler-Ferryman (1856-1930), who became an officer in the British army in West Africa and published extensively on the region. Edward Mockler himself followed his father into the Bombay Army (successor to the East India Company's Army) in 1859. From 1866 onwards, he was used in more of a political capacity, starting with a post in Aden, an important port on the Red Sea (now Yemen) as assistant to the Political resident Colonel William Merewether (1825-1880), while the latter took part in the advance party for the punitive expedition to Abyssinia (1868). During the expedition, Mockler was in charge of relations with the Chiefs of the local tribes at Zoulla (Zula, a port on the coast in Eritrea) in order to keep communication and supply lines open (Times, 26 Oct. 1867 p. 6;18 Dec. 1867 p. 4; 1 July 1868 p. 5).

Since 1873, Mockler was employed on the Makran Coast of Balochistan, which at the time was a division of the Kalat State, itself under the jurisdiction of the Baluchistan Agency (Imp.Gaz. XVII: 44f.), and in 1879 he became Assistant to the Governor General's Agent for Baluchistan.

Throughout the 1880s, Mockler was again engaged in diplomatic missions in the Gulf region, on the one hand in Basra, which became the capital of a vilayet of the Ottoman Empire in 1884 (Enc.Brit. III: 193), but also Political Agent in Muscat, the capital of Oman, which at the time was informally part of British India (ONLEY 2005: 30f.).

Returning from a leave in England, Mockler became Political Resident / Consul General in Turkish Arabia (i.e. the Ottoman provinces of Baghdad and Basra, ONLEY 2005: 37) in 1892, a post he held until his retirement (all dates from India List 1905: 568).

In Baghdad's small community of some 50 Europeans (SCHWANITZ 2009: 2), Mockler must have been a prominent member. The German diplomat and orientalist Max VON OPPENHEIM (1860OPPENHEIM ( -1946) ) provides the following description (1900: 259): "During my time there, the British Consul General undoubtedly played the most important role among the representatives of the European powers; the holder of the post at that time was Colonel Edward Mockler." The post depended on the British embassy in Constantinople and on the government of British India at the same time, so Mockler "receives a substantial income from both sides. He lives in one of the stateliest houses near to the river and commands a guard of 15 Indian sepoys, and a vessel that belongs to the Indian navy and has a crew of some 60 men." The prestige awarded to this post, comments VON OPPENHEIM, helped British political and trading interests in the region in preeminent fashion.

In 1884 Edward married Sarah, daughter of Colonel Edward Ross (Consul General in Persia, a post responsible for British interests in the whole Gulf region, ONLEY 2005: 30). They had six children, three of whom were born at Bushire (Bushehr) in Persia between 1885 and 1888. Edward retired in 1898 with the rank of Major General, and died 4 March 1927 in Guernsey in the Channel Islands.

Academic interests

During his time in posting in Balochistan, Edward Mockler took an interest in the archaeology and history of the region and published two articles in JRAS 1876, (discussing potentially old remains of settlements in the region, and noting the find of an Arsacid coin) and 1879 (suggesting relations between places mentioned in Greek sources and contemporary settlements).

At the same time, he also researched and wrote a grammar of Balochi (A Grammar of the Baloochee Language, dedicated to Col. Merewether, who in the meantime had become Commissioner in Sindh). This grammar, treating the Southern Balochi (Coastal) dialect of Makran, is in all aspects a pioneering work. Mockler mentions two manuals of Eastern Balochi being accessible to him, which he used for dialectal comparisons (cf. p. 1), but these cannot have provided a model for his grammar (GLADSTONE 1874 presents a rather rudimentary sketch of grammar on a mere five pages).

Obviously building on his own investigations, then, Mockler presents a concise grammar, which nonetheless discusses a noteworthy breadth of patterns and situations, including some that are not easily found in later and more extensive works. These include e.g. the combination -excluded in some other dialects of Balochi -of the "indefinite article" -ē with case endings, the former noteworthily preceding the latter (p. 10f.); imperatives with bil (ki) ... lit. 'let X do...' (p.53f.); competing passive constructions (p. 88f.); subjects of verbs such as "laugh", "cry" being in the oblique in the past domain (p. 20f.), and the imperfective prefix a-(which some later works mistook as merely being an epenthetic vowel).

His terminology is of course not the one employed nowadays, so some previous knowledge of Balochi helps to make full use of the multitude of interesting and linguistically very relevant details. For instance, the pronominal clitics (called "pronominal postpositions") in agent function in an ergative construction are treated as the normal case, i.e. even when an overt agent is present, thus assuming functions of verbal agreement; they are affixed either to the patient, the non-verbal element of a complex predicate, or else to the verb itself (p. 35f.).

Mockler's approach is more comprehensive than that of others in his also noting patterns that are not possible (the crucial difference between not mentioning something and saying that it does not exist, a type of information often absent), such as the reflexive pronoun being inflectable in the singular only (p. 40), as well as his lack of data on certain points, e.g. the pronominal clitics of the 1 st and 2 nd person not being frequent enough in his data to make sure they are employed in the same way as those of the 3 rd person (p.37f.).

Combining knowledge of the classical languages with a considerable amount of intuition and surely much hard work, Mockler's grammar matches VON OPPENHEIM's description (1900: 259) of the British Consuls General in Baghdad, always being officers and magistrates who over the years became well acquainted with the region and its inhabitants. "Nearly all of them have rendered outstanding services to scholarship by their own works and the promotion of the research of their compatriots".


A Survey of Parallel Solvers for the Finite Element Method in Computational Electromagnetics

Keywords: Finite element methods, Parallel algorithms, linear systems, electromagnetic analysis

New trends in parallel methods used to solve Finite Element matrix systems are presented: standard iterative and direct solving methods first, and then domain decomposition methods. For example, the current status and properties of two prevailing programming environments (PVM and MPI) are finally given and compared when implemented together with a Finite Element Time Domain formulation.

Introduction

Numerical computation is more and more used in engineering sciences to develop new devices or to optimize existing one. Only parallel computers provide the increase in computing performances necessary to solve today's problems. Two reasons may be highlighted: large memory is required because of a large amount of data, or speed is required to obtain the solution [START_REF] Agarwal | Parallel Computers and Large Problems in Industry[END_REF].

We have previously presented a survey of the parallelization of numerical techniques used in computational electromagnetics [START_REF] Jacques | Le calcul de champs électromagnétiques sur architectures MIMD[END_REF]. In this paper, it is proposed to focus now our attention on numerical solving methods for the parallelized Finite Element Method (FEM) in electromagnetics [START_REF] Butrylo | New trends in parallel electromagnetic fields computation[END_REF]. The objective is actually to describe lately appeared methods which can be useful to solve efficiently linear systems issued from FE discretization. The main characteristic of such systems is that they are generally sparse.

The FEM is attractive for modeling 2D and 3D devices because of its adaptability for designing complex systems with highly heterogeneous materials. Modeling realistic devices requires the solution of a large problem which can be afforded only by parallel computation [START_REF] Vollaire | Parallel Computing for the Finite Element Method[END_REF].

The first section deals with standard numerical solving methods, such as preconditioned conjugate gradient. In the second section, domain decomposition methods (DDM) are described. Since these methods are less known in computational electromagnetics, their basic principle is presented. Special attention is also paid to multilevel methods, which can be seen as a DDM method from a strict mathematical point of view. As example, the last section shows the implementation of a Finite Element Time Domain (FETD) algorithm used for a wave propagation problem. A comparison between the two message passing libraries PVM and MPI is also performed.

Standard Parallel Solving Methods for the FEM

The Conjugate Gradient Method

The Conjugate Gradient (CG) method is one of the most popular linear solvers for Finite Element (FE) analysis, since it is well adapted to solve sparse matrix systems. In order to improve the solving, preconditioning techniques have to be used together with the solving method. The basic idea consists in

replacing the initial system A x = b by a new system b M b ~ and A M A with b x A ~ b Ax 1 1 - - = = = ⇔ = ( 1
)
where M is the preconditioning matrix. Two conditions are required for efficiency:

-the conditioning number of the new matrix à has to be lower than the initial conditioning number, -and the matrix system M has to be to not expensive to build or to solve.

Several preconditioning techniques exist: Incomplete Cholesky (ICP), SSOR, Diagonal (DP) preconditioning. The following example shows the efficiency of the preconditioning method. It is used for the modeling of radiofrequency hyperthermia [START_REF] Siauve | 3D Modeling of the SAR Distribution for RF Hyperthermia with External Waveguide Applicator[END_REF]. The problem is computed with more than 130 000 degrees of freedom. From Table I, by comparing both DP and ICP, it is shown that better results are obtained with the ICP in terms of number of iterations. On the other hand, since the cost of one iteration is greater than when using the DP, the total CPU time of computation is larger. In this case, the best results are obtained with a SSOR preconditioning: both number of iterations and CPU time are divided by two. Several papers published recently propose new parallel-processing techniques of the preconditioned CG solver.

In [START_REF] Iwashita | Parallel Processing of 3-D Eddy Current Analysis with Moving Conductor Using Parallelized nICCG Solver with Renumbering Process[END_REF], a parallelized transient eddy current FE analysis is presented. Due to the presence of moving conductors, the right-hand side of the system of equations is changed at each time step. First-order bricktype edge elements are used for the FE discretization. Parallel solving is performed using the ICCG, leading to difficulties due to the forward-backward substitutions. Two parallel implementations are compared. First, on N p processors, the block ICCG (BICCG) method divides the global matrix into N p submatrices and performs the IC factorization for the local submatrix in each processor, with the entries between different processors being ignored in the factorization. Second, a renumbering process allows to reorder the matrix into a form similar to the dissection ordering case (PICCG-RP) (fig. 1). Since no entries are neglected during the preconditioning step, a better convergence rate is observed with the PICCG-RP method, leading to a better parallel efficiency (fig. 2).

P1

Interior This reordering technique is extended in [START_REF] Iwashita | Algebraic Multicolor Ordering for Paralleized ICCG Solver in Finite-Element Analyses[END_REF] for the same 3D eddy-current analyses. The key of the method is coloring of the unknowns. Several colors are assigned to the unknowns, such that the unknowns belonging to the same color must have no data dependency on each other. Due to the presence of edge elements, the assignment is algebraic: it is carried out by using only information obtained from a coefficient matrix. This method is actually a complete black-box solver, and its implementation appears to be easy. Parallel performances obtained on a distributed memory Fujitsu VPP-800 are also better than when using a classical BICCG (fig. 3). It is also shown that the use of many colors improves the convergence rate, but it causes an increase in communication costs and a decline in parallelism. For the example of transient eddy currents analysis, made of more than 10 6 unknowns, best results are obtained with 60 colors. The Conjugate Gradient Squared (CGS) is also used for the modeling of 2D open-region frequencydomain scattering problems by the FEM in [START_REF] Maheswaran | MCGS: a Modified Conjugate Gradient Squared Algorithm for Nonsymmetric Linear Systems[END_REF]. By modifying the computation sequence in the algorithm, the synchronization overhead is reduced by a factor of two. The basic idea of this modification is to merge the inner products present in the CGS algorithm so that they can be evaluated using a single global summing operation per iteration of the "while" loop. From both theoretical and experimental analyses it is found that this Modified CGS (MCGS) algorithm becomes slightly faster as the number of processors increases. A set of algorithms is then proposed, where either CGS or MCGS is selected depending on the number of unknowns and number of processors.

Other parallel standard iterative methods

Other parallel iterative solvers may be used for the FE method and can be more efficient than the preconditioned CG. As example, the Quasi-Minimal Residual (QMR) method and the BI-Conjugate Gradient STAbilized (BICGSTAB) method are compared in [START_REF] Iwashita | Parallel Solvers for Electromagnetic Field Computations on Printed Circuit Boards[END_REF] for the computation of 2D axisymmetric magnetic and electric fields on printed circuit board (PCB) containing magnetic films. Both solvers and preconditioners are parallelized by matrix and vector partitioning, the mesh partition being obtained by using METIS [START_REF] Karypis | Analysis of Multilevel Graph Partitioning[END_REF]. The results obtained on different parallel machines show the efficiency of this approach (Table II).

Solver Speedup

ILUT perconditioned BICGSTAB 14.4

Diagonaly scaled BICGSTAB 22.9

Diagonaly scaled QMR 31.9

Direct methods

Iterative solvers are not suitable for any type of problem. For example, the matrix equations obtained from the FEM with perfect matched layers (PML) can be rather ill-conditioned and can lead to very slow or non-convergent results. Direct solvers have then to be employed, since they are much more robust, reliable and efficient. As example, a sparse LU decomposition solver is presented in [START_REF] Sheng | Parallel Electromagnetic Modeling of 3D Microstrip Discontinuities using FEMand PML[END_REF]. It is used to compute 3D microstrip discontinuities. Several implementation techniques are used to speedup the computation: dynamic partial pivoting, reordering scheme, supernodal approach. The parallel version of the package is implemented on a PC-based parallel platform, leading to decent speedups (fig. 4).

Domain Decomposition Methods

A natural way to create solvers for partial differential equations on parallel computers is the decomposition technique based on a decomposition of the spatial domain of the problem into subdomains. Most domain decomposition methods can be classified as either overlapping or nonoverlapping subdomain approach. For storage and time efficiency, iterative algorithms are preferred to direct solvers. To make the methods numerically scalable, a global coupling between subdomains must be introduced, using a coarse global mesh for example.

Overlapping domains methods

Overlapping domains methods are easy to implement as additive Schwarz iteration. In this method, the computation is first performed on each subdomain with known boundary condition and with unknown boundary condition on the interface with other subdomains. From this computation, new boundary conditions are obtained on the interfaces and transferred to the other subdomains. This process is made iteratively until convergence (fig. 6). On each subdomain, an exact solver (LU) or an iterative solver can be used. Of course, this method is easy to parallelize, since each subdomain can be computed in parallel. So it is well adapted to coarse grain parallelism.

Ω 2 Ω 1 Ω 2 Ω 1 Ω 2 Ω 1

Compute with exact and unknown BC 2. Get BC from previous results

Compute with approximated BC

unknown BC iteration Fig. 6. Basic idea of the additive Schwarz iteration.

On the other hand, the performance of multiplicative Schwarz iteration is often better, since the unknown boundary conditions are faster approximated (fig. 7). However this method seems to be more difficult to parallelize, since it appears to be sequential by nature. The technique of multicoloring must then be used to improve its parallel efficiency .

Ω

Non overlapping domain methods

Algorithms with non overlapping domains Schur complement systems. If the system matrix is made of four blocks, the Schur complement matrix S is defined by:

B A B' - C S 1 - = ⇒ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ C B' B A (2)
It is then shown that solving the global system leads first to solve the Schur system:

( ) ⎩ ⎨ ⎧ - = - = ⇒ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - - - 1 1 2 1 - 1 1 1 2 1 b A B' b y B A B' - C By A b A x b b y x C B' B A (3)
The unknowns located on the interface separating two neighbouring subdomains give the Schur complement system. Since the subdomains Ω i become then disconnected, the corresponding block matrices A i become also disconnected, allowing an easy parallelization (fig. 8). The Schur matrix system is solved by the CG method for symmetric definite matrix or more generally by Krylov subspace method such as GMRES. Parallel interface preconditioners are often introduced to improve the convergence.

Many ways for preconditioning exist in the literature for 2D problems. Preconditioners which not destroy the original scalability and which are not expensive are complicated to obtain for complex geometries and irregular discretizations. A smart manner to obtain non overlapping domain decomposition algorithm showing a good parallel efficiency is to enforce the continuity of boundary conditions by introducing Lagrange multipliers [START_REF] Zyserman | Parallel finite element algorithm with domain decomposition for three-dimensional magnetotelluric modelling[END_REF].

The FE Tearing and Interconnecting method is an important class of such hybrid methods [START_REF] Farhat | A Method of Finite Element Tearing and Interconnecting and its Parallel Solution Algorithm[END_REF]. It can be viewed as a dual Schur Complement method . The Lagrange multipliers are solution of a small-size linear system compared to the global problem size. When solved with an iterative method (CG or BiCG method), the corresponding matrix does not need to be explicitly assembled. Each local field problem can be treated independently on each processor. In [START_REF] Navsariwala | An Efficient Implementation of the Finite-Element Time Domain algorithm on parallel Computers using a Finite-Element Tearing and Interconnecting Algorithm[END_REF] the resulting algorithm for the solution of the vector wave equation discretized with tetrahedral elements is presented . It is used to modelize 3D PEC wedge in a rectangular cavity. The continuity of tangential electric fields on the interface between subdomains is used to build the Schur system. This system is calculated using parallel conjugate gradient, while the local matrices A i factorized only once. From Fig. 9, this algorithm is shown to be highly scalable.

An other interesting advantage of the FETI method follows from [START_REF] Wolfe | A parallel Finite-Element Tearing and Interconnecting Algorithm for Solution of the Vector Wave Equation with PML Absorbing Medium[END_REF]. For the solution of Maxwell's equation in the frequency domain, its performance is less affected by the presence of perfectly matched layers solution than a direct iterative method (a divide-and-conquer-type preconditioned BiCGM). On the other hand, the convergence of the method is sensitive to irregular subdomains. The question of constructing efficient parallel preconditionners is also a remaining problem for such a method. The domain decomposition approach offers many opportunities for the coupling of boundary element method (BEM) and FEM. An unified variational formulation has been proposed in [START_REF] Costabel | Symmetric methods for the coupling of finite elements and boundary elements[END_REF]. In [START_REF] Heise | Nonlinear simulation of electromagnetic fields with domain decomposition methods on MIMD parallel computers[END_REF][START_REF] Kuhn | The application of coupled BE/FE formulations in technical magnetic field computations[END_REF], the z-component of the vector potential is determined from the solution of linear, symmetric but indefinite system of coupled equations which are reformulated as symmetric and positive definite systems. A different way to exploit the BEM-FEM coupling is presented in [START_REF] Rischmüller | 3D Transient Analysis of Electromagnetic Devices Using Parallel BEM Coupled to FEM[END_REF] for 3D transient problem. In this case, the coupling with the BEM is made necessary due to the presence of moving part into the electromechanical device. The developed method is based on a Dirichlet-Neumann algorithm: sequential solutions of the FEM part are linked within a preconditioned GMRES method together with parallel solution of the BEM part (fig. 10). On a NEC-SX4 with 16 processors, the obtained speedup is about 6.

Multilevel methods

Multilevel methods can be included in the class of decomposition method, since these algorithms have similarities with Schwarz type domain decomposition method with inexact solvers. Different grid levels or subspaces play the role of subregions. Schematically, each iteration may be seen as performed into 3 steps (fig. 11). A few iterations are first performed on a fine grid, leading to an approximative solution u h . This leads to a residual vector r h which is transferred to the coarse grid and then exactly computed (second step). The solution is finally re-transferred onto the fine grid, allowing to modify the previous values of u h . It is followed by a post-smoothing step.

Multilevel methods can be generalized to several grid sizes. Two types of multigrid methods exist [START_REF] Ruge | Algebraic multigrid[END_REF]:

-the Geometric Multigrid Method (GMM), based on the effective use of mesh grids, -the Algebraic Multigrid Method (AMG), where there is no mesh, and the building of the operators is based on the values of the coefficients of the matrix. Parallel implementation of the AMG is presented in [START_REF] Mifune | A Fast Solver for FEM Analyses Using the Parallelized Algebraic Multigrid Method[END_REF] for eddy-current analyses and magnetostatics analyses. The parallel performance of the algorithm is improved in the setup phase by the use of longrange interpolation instead of the conventional direct interpolation. The numerical results obtained with a rectangular domain show that the AMG is a fast solver, since the CPU time is divided by 10 compared to a classical ICCG (Table III). As shown in Fig. 12, this intrisic property is not destroyed by the parallelization of the method, since both speedups remain comparable. It would be of great interest to test the AMG method with more complicated domains, since the GMG requires a hierarchy of coarse grids which is actually difficult to manage in 3D. An experience of implementation of the GMG has been reported in [START_REF] Haase | Parallel multigrid 3D Maxwell solvers[END_REF] for 3D magnetostatics. Note that all these methods (domain decomposition or multigrid methods) can be used as preconditioners for the Conjugate Gradient algorithm [START_REF] Xu | Iterative methods by space decomposition and subspace correction[END_REF].

An example: the Finite Element Time Domain method for wave propagation

Flexibility of the finite element technique and right physical sense of the edge elements make this formulation useful in modeling of electromagnetic phenomena. High performance simulation of the time domain problem enables to reduce either memory cost or time of computation. As example of implementation, the properties of the vector finite element time domain algorithm are evaluated in the known distributed multi-computer environments MPI and PVM [START_REF] Buyya | High Performance Cluster Computing[END_REF][START_REF] Geist | PVM and MPI: a Comparison of Features[END_REF][START_REF] Pacheco | Parallel Programming with MPI[END_REF][START_REF] Sunderam | Heterogeneous Parallel and Distributed Computing[END_REF].

Finite Element formulation

The distribution of electric field is given by time dependent vector wave equation , 0

= ∂ ∂ + ∂ ∂ + ∂ ∂ + × ∇ × ∇ t t t imp J E E E ε σ μ 1 2 2
where μ, σ, ε are the parameters of the media and J imp is the source current. The domain of analysis is truncated with first order Engquist-Majda absorbing boundary condition (ABC) [START_REF] Sunderam | Heterogeneous Parallel and Distributed Computing[END_REF]. In the case of the excitation by a plane wave, according to general Galerkin scheme, the weak form of eq. ( 4) is given by expression ( )( )

, 1 1 1 2 2 ∫ ∫ ∫ ∫ ∫ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ × ∇ × - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ × ∂ ∂ = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ × ∂ ∂ + ∂ ∂ + ∂ ∂ + × ∇ × ∇ ABC ABC S i i S V V V t c dS t c t dV t dV E n n E W n E W E W E W W E r r r μ μ ε σ μ (5)
where W is the test vector function, S ABC is the external surface of the model and E i is the incident field. This equation is discretized in time domain and in space domain to yield a system of linear equations which must be solved. Space discretization is achieved using incomplete first order H(curl, Ω) tetrahedral edge elements. Considering the central Euler difference approximations of the first and the second order derivatives, the final form of the equation is ( )

, 2 2 2 2 1 2 1 n n n n t t t t f E T R E S T E R T Δ - ⋅ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - Δ + ⋅ Δ - = ⋅ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Δ + - + (6)
where R, T, S are mass, damping and stiffness matrices. The f n vector represents the dynamic load in the analyzed model.

Distributed implementation of the formulation

The Single Program Multiple Data (SPMD) paradigm is used in the elaborated parallel code. Because the unknowns are connected with the edges, the set of edges is decomposed in the distributed algorithm.

Since the matrices are assembled with respect to the degrees of freedom, there are no geometrical restrictions of the data decomposition.

Each computing unit assembles the matrices only for its local subset of degrees of freedom. During the assembling stage, the processing units have to exchange information only about the spatial distribution of boundary conditions (fig. 13). The numerical integration of ( 5) requires the solution of the linear system A x n = b n (fig. 14). Conjugate gradient method with diagonal preconditioning is used. The basic matrix and vector operations in the solver procedure are parallelized. The interdependence of the processes is different in the assembling stage and in the time integration loop. The assembling thread in a processing unit is loosely connected with the others threads, while the subtasks in the distributed implementation of the CG algorithm are tightly coupled.

Processing unit P i+1

Create and init the data structures: ♦ unknowns: 3

Comparison of PVM and MPI performances

Both implementations of the FEM algorithm with PVM and MPI have been tested on an heterogeneous cluster of 4 PC. The modelized problem is the propagation of a plane wave in free space. These results indicate better numerical performance of PVM version of the presented FETD algorithm. However, the configuration of the hardware seems to have significant influence on the performance of the parallel algorithm.

Conclusion

It appears from this review that the parallelization of standard solving methods, such as Conjugate Gradient method preconditioned with classical technique (SSOR, Incomplete Cholesly, …), has reached a sufficient maturity level. On the other hand, progress have to be made regarding the use of new methods, such as domain decomposition techniques or multigrid solvers. These methods may be used themselves as preconditioning techniques for Krylov subspace methods. Domain decomposition methods may be classified into overlapping and non-overlapping methods. Regarding the overlapping methods, additive Schwarz method is easy to parallelize, and it is well adapted to coarse grain parallelism. On the other hand, multiplicative Schwarz iteration is often better, but it appears more difficult to parallelize. Non overlapping domains involve Schur complement systems, and a good parallel efficiency may be obtained by introducing Lagrange multipliers to enforce the continuity of boundary conditions. Multilevel methods can be included in the class of decomposition method. These methods are generally efficient when implemented on sequential computer, and only some parallel implementations have been reported yet. However more complex geometries and domains have to solved in order to obtain more significant results.

Like in many other disciplines, computational electromagnetics requires a computing power considerably higher than that offered by today's conventional supercomputers, and parallel computing is destined to have a large influence in the next scientific developments.

Fig. 1 .

 1 Fig.1. Basic principle of the renumbering process[START_REF] Iwashita | Parallel Processing of 3-D Eddy Current Analysis with Moving Conductor Using Parallelized nICCG Solver with Renumbering Process[END_REF].

Fig. 2 .

 2 Fig.2. Comparison of several preconditioning techniques for the CG solver[START_REF] Siauve | 3D Modeling of the SAR Distribution for RF Hyperthermia with External Waveguide Applicator[END_REF]. Computation is performed on a Hitachi SR2201.

Fig. 3 .

 3 Fig. 3. Comparison of multi-color method for the CG solver and BICCG [6] for a transient eddy currents problem (1011920 unknowns).

Fig. 4 .

 4 Fig.[START_REF] Vollaire | Parallel Computing for the Finite Element Method[END_REF]. LU decomposition for the computation of 3D microstrip discontinuities (21408 dof)[START_REF] Sheng | Parallel Electromagnetic Modeling of 3D Microstrip Discontinuities using FEMand PML[END_REF].

1 Fig. 5 .

 15 Fig.5. Overlapping (left) and nonoverlapping (right) domains.

Fig. 8 .

 8 Fig. 8. Basic idea of the Schur complement system method.

Fig. 9 .

 9 Fig.9. Speedup of the FETI algorithm for the computation of a PEC wedge into a cavity[START_REF] Farhat | A Method of Finite Element Tearing and Interconnecting and its Parallel Solution Algorithm[END_REF]. Computation is performed on a 32-node Intel iSPC/860 hypercube.

Fig. 10 .

 10 Fig. 10. BEM-FEM coupling (left) to modelize electromechanical device (right) with moving part [19].

Fig. 11 .

 11 Fig. 11. Basic idea of the multilevel method (example with 2 grids).

Fig. 14 .

 14 Fig. 14. Relation between two processes in the time integration loop.

  Several mesh sizes have used for the comparison. The total time of execution of the PVM version is about 20% longer than the MPI implementation of the algorithm. However the relative parameters of PVM implementation are better (fig.[START_REF] Farhat | A Method of Finite Element Tearing and Interconnecting and its Parallel Solution Algorithm[END_REF][START_REF] Navsariwala | An Efficient Implementation of the Finite-Element Time Domain algorithm on parallel Computers using a Finite-Element Tearing and Interconnecting Algorithm[END_REF]. The speedup of MPI and PVM implementations depends on the numbers of unknowns in the analysed model. If the relation between time of computation and time of data transfer is small (i.e. for relatively small model) the PVM and MPI implementations do not cause the speedup of the calculation process. The performance of the MPI version is shown to be a non-monotone function. The model of medium size gives the best speedup of the computation (about 2.7 for 4 processing units). The speedups of PVM and MPI are approximately equal only for this type of the FEM model. The efficiency and speedup of MPI version decrease when the limit of RAM memory in the cluster is achieved. The efficiency of MPI version converges to 0.5 for the largest FEM model. In this case the speedup and efficiency of PVM implementation achieve the ideal values.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Speedup of the FETD algorithm-A: 6930 DOF, B: 77832 DOF, C: 187488 DOF.

TABLE I

 I

			25761 s
	SSOR	3481	12225 s
	: COMPARISON OF THE EFFICIENCY OF SEVERAL PRECONDITIONING TECHNIQUES FOR THE CG
		METHOD
	Preconditioning method	Number of iterations	CPU time
	Diagonal	13561	22356 s

TABLE II :

 II SPEEDUP ON 48 PROCESSORS FOR A 16593 DOF MATRIX[START_REF] Iwashita | Parallel Solvers for Electromagnetic Field Computations on Printed Circuit Boards[END_REF]. COMPUTATION ON A NEC CENJU-4.

TABLE III :

 III COMPARISON OF CPU TIMES (IN S) OBTAINED WITH THE AMG AND ICCG SOLVERS ON A 80601DEGREES OF FREEDOM PROBLEM.

	Number of processors		AMG	ICCG
		1		513	5219
		2		273	3163
		4		153	1748
		8		97.2	998
		16		67.5	630
		16,0
			AMG
		12,0	ICCG
	Speedup	8,0	Linear
		4,0
		0,0
		0	4	8	12	16
			Number of processors
		Fig. 12. Speedup of AMG and ICCG methods [19].