Thierry Daudé
email: thierry.daude@u-cergy.fr

Niky Kamran
email: nkamran@math.mcgill.ca

François Nicoleau
email: francois.nicoleau@math.univ-nantes.fr

Stability in the inverse Steklov problem on warped product Riemannian manifolds

Keywords: Inverse Calderón problem, Steklov spectrum, moment problems, Weyl-Titchmarsh functions, local Borg-Marchenko theorem. 2010 Mathematics Subject Classification. Primaries 81U40, 35P25; Secondary 58J50

In this paper, we study the amount of information contained in the Steklov spectrum of some compact manifolds with connected boundary equipped with a warped product metric. Examples of such manifolds can be thought of as deformed balls in R d . We first prove that the Steklov spectrum determines uniquely the warping function of the metric. We show in fact that the approximate knowledge (in a given precise sense) of the Steklov spectrum is enough to determine uniquely the warping function in a neighbourhood of the boundary. Second, we provide stability estimates of logtype on the warping function from the Steklov spectrum. The key element of these stability results relies on a formula that, roughly speaking, connects the inverse data (the Steklov spectrum) to the Laplace transform of the difference of the two warping factors.

In this paper, we shall consider a class of d-dimensional manifolds

M = (0, 1] × S d-1 , (1.1)
where d ≥ 3 and with a connected boundary consisting in a copy of S d-1 , i.e. ∂M = {1} × S d-1 . We shall assume these manifolds are equipped with Riemannian warped product metrics of the form

g = c 4 (r)[dr 2 + r 2 g S ]. (1.2)
In what follows, g S denotes a fixed smooth Riemannian metric on the sphere S d-1 and c is a positive C m -function, (with m ≥ 2), of the variable r only. Without loss of generality, we may assume c(0) = 1.

We emphasize that under these general assumptions, the metrics (1.2) are not necessarily regular. Actually, one proves that the metrics g are regular if and only if the odd-order derivatives c (2k+1) (0) = 0 and g S = dΩ 2 where dΩ 2 is the round metric on S n-1 , (see [START_REF] Petersen | Riemannian Geometry, Third Edition[END_REF], Section 4.3.4). Nevertheless, in this paper, we consider general metrics of the form (1.2) which could be singular at r = 0, but emphasizing that in the regular case, the proofs of our results are much simpler.

It is convenient to use the change of coordinates x = -log r ∈ [0, +∞[. In these new coordinates, the metric g has the form g = f 4 (x)[dx 2 + g S ],

(

where f (x) = c(e -x )e -x 2 . Using the Taylor expansion of c(r) at r = 0, we see that the conformal factor f (x) satisfies the following asymptotic expansion:

f (x) = e -x 2 + m k=1 c k e -(k+ 1 2 )x + o(e -(m+ 1 2 )x ) , x → +∞ , (1.4)
for some suitable real constants c k .

A typical example of Riemannian manifold belonging to our class is the usual Euclidean metric on the unit ball of R d , obtained by simply taking c(r) = 1, (i.e f (x) = e -x 2 ), and g S = dΩ 2 . Hence, under our hypothesis, the Riemannian manifold (M, g) can be viewed topologically as a unit ball of R n whose metric is a deformation of the usual Euclidean metric both in the radial direction r through the warping function c(r), and also in the transversal directions through the general metric g S .

In this paper, we are interested in studying the amount of information contained in the Steklov spectrum associated to the class of Riemannian manifolds (M, g). Recall [START_REF] Girouard | Spectral geometry of the Steklov spectrum[END_REF] that the Steklov spectrum is defined as the spectrum of the Dirichlet-to-Neumann map Λ g (abbreviated later by DN map) associated to (M, g). More precisely, consider the Dirichlet problem

-g u = 0, on M, u = ψ, on ∂M, (1.5)
where ψ ∈ H 1 2 (∂M ). Using the separation of variables, we shall prove in Section 2 that we can solve uniquely (1.5) even in the case where the metric g is singular. Of course, when the metric is smooth, this result is well-known (see [START_REF] Salo | The Calderón problem on Riemannian manifolds, Inverse problems and applications: inside out[END_REF][START_REF] Taylor | Partial Differential Equations, I. Basic theory[END_REF]) and the unique solution u of (1.5) belongs to the Sobolev space H 1 (M ).

In the latter case, we define the DN map Λ g as the operator Λ g from H 1/2 (∂M ) to H -1/2 (∂M ) as

Λ g ψ = (∂ ν u) |∂M , (1.6)
where u is the unique solution of (1.5) and (∂ ν u) |∂M is the normal derivative of u with respect to the outer unit normal vector ν on ∂M . Here (∂ ν u) |∂M is interpreted in the weak sense as an element of H -1/2 (∂M ) by

Λ g ψ|φ = M
du, dv g dV ol g , for any ψ ∈ H 1/2 (∂M ) and φ ∈ H 1/2 (∂M ) such that u is the unique solution of (1.5) and v is any element of H 1 (M ) such that v |∂M = φ. If ψ is sufficiently smooth, we can check that

Λ g ψ = g(ν, ∇u) |∂M = du(ν) |∂M = ν(u) |∂M ,
so that the expression in local coordinates for the normal derivative is thus given by

∂ ν u = ν i ∂ i u. (1.7)
The DN map is a pseudo-differential operator of order 1 and is self-adjont on L 2 (∂M, dS g ) where dS g denotes the metric induced by g on the boundary ∂M . Therefore, the DN map has a real and discrete spectrum accumulating at infinity. We shall thus denote the Steklov eigenvalues (counted with multiplicity) by

0 = σ 0 < σ 1 ≤ σ 2 ≤ • • • ≤ σ k → ∞. (1.8)
We refer the reader to the nice survey [START_REF] Girouard | Spectral geometry of the Steklov spectrum[END_REF] and references therein. The main results of this paper are the following. First, we obtain the complete asymptotics of the Steklov spectrum σ k as k → ∞ in terms of the Taylor series of the effective potential

q f (x) = (f d-2 ) (x) f d-2 (x) - (d -2) 2 4 , (1.9)
and of the eigenvalues of the positive Laplace-Beltrami operator of (S d-1 , g S ). Note that, in the radial coordinate r, the effective potential is given by:

q f (x) = (d -2) (d -3)r 2 c (r) c(r) 2 + r 2 c (r) c(r) + (d -1)r c (r) c(r) . (1.10)
For later use, we denote the spectrum ofg S (counting multiplicities) by

0 = µ 0 < µ 1 ≤ µ 2 ≤ • • • ≤ µ k → ∞. (1.11)
For each k ≥ 0, the eigenvalue µ k is associated to the normalized eigenfunction

Y k ∈ L 2 (S d-1 , dV ol g S ) such that -g K Y k = µ k Y k .
Let us finally introduce

κ k = µ k + (d -2) 2 4 , ∀k ≥ 0, (1.12)
and recall the Weyl asymptotics ( [START_REF] Safarov | The asymptotic distribution of eigenvalues of partial differential operators[END_REF], Theorem 1.2.1 and Remark 1.2.2)

κ k = c d-1 k 1 d-1 + O(1), k → ∞. (1.13)
Here

c d-1 = 2π (ω d-1 V ol(S d-1 )) -1 d-1 where ω d-1 is the volume of the unit ball in R d-1 , while V ol(S d-1
) is the volume of S d-1 for the metric g S . Our first result is the following: Theorem 1.1. Let (M, g) be a Riemannian manifold given by (1.1)-(1.2) and assume that c ∈ C ∞ ([0, 1]). Then the Steklov spectrum (σ k ) k≥0 satisfies for all N ∈ N,

σ k = (d -2)f (0) f 3 (0) + κ k f 2 (0) + N j=0 β j (0) f 2 (0) κ -j-1 k + O(κ -N -2 k ), k → ∞, (1.14)
where

     β 0 (x) = 1 2 q f (x), β j+1 (x) = 1 2 β j (x) + 1 2 j l=0 β l (x)β j-l (x). (1.15)
We make the following observations:

1. Using to (1.13), we find that the Steklov spectrum satisfies the usual Weyl law (see [START_REF] Girouard | Spectral geometry of the Steklov spectrum[END_REF]):

σ k = c d-1 f 2 (0) k 1 d-1 + O(1) = 2π k ω d-1 V ol(∂M ) 1 d-1 + O(1), k → ∞.
(1.16)

2. Note that the coefficients β j (0), j ≥ 0 only depend on the derivatives q (l) f (0), l = 0, . . . , j up to order j. Hence the values f (0), f (0) and the Taylor series at 0 of the effective potential q f give the leading terms of the asymptotics of the Steklov spectrum in inverse powers of κ k . But without additional knowledge of the asymptotics of the κ k , we are not able to determine the warping function f from this information. 3. In the exceptional case where (S d-1 , g S ) = (S d-1 , dΩ 2 ), the spectrum ofg K is given by {λ

k = k(k + d -2), k ≥ 0} where each eigenvalue λ k has multiplicity 2k+d-2 d-2 k k+d-3 .
Hence, if we order the Steklov spectrum σ k without counting mutiplicities, i.e.

0 = σ 0 < σ 1 < σ 2 < • • • < σ k → ∞, 4 where each σ k has multiplicity 2k+d-2 d-2 k k+d-3 and denoting κ k = λ k + (d-2) 2 4 = k + d-2
2 , we get the more precise asymptotics (without multiplicity):

σ k = (d -2)f (0) f 3 (0) + k + d-2 2 f 2 (0) + N j=0 β j (0) f 2 (0) (k + d -2 2 ) -j-1 + O(k -N -2 ), k → ∞. (1.17)
As a consequence, the asymptotics of the Steklov spectrum allow to recover inductively the values f (0), f (0) and the Taylor series at 0 of the effective potential q f . Using then (1.9), we immediately see that in fact the Steklov spectrum determines the Taylor series of f . Hence we have proved:

Corollary 1.2. If (S d-1 , g S ) = (S d-1 , dΩ 2
) and if the warping function c is analytic on [0, 1], then the Steklov spectrum of (M, g) determines uniquely the warping function c.

This leads to the question: does the Steklov spectrum determine uniquely the warping function c without any assumption of analyticity on c and on the transversal metric g S ? The answer is yes. Precisely, we have: Theorem 1.3. Let (M, g) be a Riemannian manifolds given by (1.1)-(1.2). Then the Steklov spectrum (σ k ) k≥0 determines uniquely the warping function c.

In dimension 2, the problem of recovering the metric from the Steklov spectrum on the unit disk was studied in a paper by Jollivet and Sharafutdinov [START_REF] Jollivet | On an inverse problem for the Steklov spectrum of a Riemannian surface[END_REF]. They show that if the DN maps associated to the metrics g j , j = 1, 2 are intertwined (which is a stronger hypothesis than the Steklov isospectrality), then the metrics g j are equal (up to a diffeomorphism and a conformal factor which is equal to 1 on the boundary of the unit disk).

In fact, we prove an even better uniqueness result than Theorem 1.3: the subsequence (σ k d-1 ) determines uniquely the warping function c. The proof is based on the following local uniqueness property: Theorem 1.4. Let (M, g) and (M, g) be Riemannian manifolds given by (1.1)-(1.2). Then, for a positive constant a > 0, the two following assertions are equivalent:

σ k d-1 -σk d-1 = O(e -2aκ k d-1 ), k → ∞.
(1.18)

c(r) = c(r), ∀r ∈ [e -a , 1]. (1.19)
This theorem asserts that the knowledge of a subsequence of the Steklov spectrum up to an error that decays exponentially as k → ∞ (see (1.18)) determines uniquely the warping function c in a neighbourhood of the boundary r = 1 provided that (S d-1 , g S ) is known. Of course, Theorem 1.3 is an immediate consequence of Theorem 1.4. The question of whether the Steklov spectrum determines uniquely the warping function without the additional knowledge of the transversal Riemannian manifold (S n-1 , g S ) remains open.

Notice that the above result can be viewed as a weak form of a stability result. It stresses the important fact that the asymptotic behaviour of the Steklov spectrum allows one to determine the warping function c in a neighbourhood of the boundary r = 1. In contrast, we shall show that the first Steklov eigenvalues determine in a stable way the warping function c on [0, 1]. Before stating our main result, let us define our set C(A) of admissible warping functions c(r), where A is any positive constant, m ≥ 3 and 2 ≤ p ≤ m-1:

C(A) = {c ∈ C m ([0, 1]) s. t. ||c|| C m ([0,1]) + || 1 c || C m ([0,1]) ≤ A and for k = 1, ..., p -1, c (k) (0) = 0}. (1.20)
Roughly speaking, the conditions on the derivatives of c(r) at r = 0 tell us that the metric g is relatively flat in a neighborhood of the origin r = 0. These assumptions are only introduced to take into account the regularity of the warping function in the stability estimates obtained by integration by parts. We only assume m ≥ 3 for the simplicity of the proofs. In particular, using (1.10) we see that, if c ∈ C(A), the effective potential q f satisfies the uniform estimates:

|q (k) f (x)| ≤ C A e -px , ∀x ≥ 0, ∀k = 0, ..., m -2, (1.21)
where the constant C A depends only on A.

In Section 4, we shall prove two log-type estimates results. The first one concerns the case where the metric g is smooth, i.e we assume that the odd order derivatives of the warping functions c (2k+1) (0) = 0, and we assume that the transversal metric g S = dΩ 2 is the usual Euclidean metric on S n-1 . In this regular case, we have the following estimate: Theorem 1.5. Let (M, g) and (M, g) be smooth Riemannian manifolds given by (1.1)-(1.2) with c, c ∈ C(A) where A > 0 is fixed. Assume that for > 0 small enough, one has:

sup k≥0 |σ k -σk | ≤ . (1.22)
Then, there exists a positive constant C A , depending only on A such that,

||c -c|| L ∞ (0,1) ≤ C A 1 log( 1 ) p-1 (1.23)
This result is relatively close in spirit to the logarithmic stability estimates obtained by Alessandrini [START_REF] Alessandrini | Stable determination of conductivity by boundary measurements[END_REF], and then improved by Novikov [START_REF] Novikov | New global stability estimates for the Gelfand-Calderón inverse problem[END_REF], for the Schrödinger equation on a bounded domain M in R d . Indeed, it is well-known there is a close connection between the DN map which we have defined in the introduction for the metric g = c 4 (r) g e , (g e being the Euclidean metric on R d ), and the DN map associated with the Schrödinger equation (-∆ ge + V )u = 0 where the potential V is defined below. It is induced by the transformation law for the Laplace-Beltrami operator under conformal changes of metric:

-∆ c 4 ge u = c -(d+2) (-∆ ge + V ) c d-2 u , (1.24)
where

V = c -d+2 ∆ ge c d-2 . (1.25)
We recall that if 0 is not an eigenvalue of -∆ ge + V , then for all ψ ∈ H 1/2 (∂M ), there exists a unique u solution of (-∆ ge + V )u = 0, on M, u = ψ, on ∂M.

(1.26)

The DN map for this Schrödinger equation is then given by Λ V (ψ) = ∂ ν u |∂M , where (∂ ν u) |∂M is its normal derivative with respect to the unit outer normal vector ν on ∂M .

If the warping function c(r) satisfies c |∂M = 1 and ∂ ν c |∂M = 0, (or equivalently c(1) = 1 and c (1) = 0 in our case), one can easily show (see for instance [START_REF] Salo | The Calderón problem on Riemannian manifolds, Inverse problems and applications: inside out[END_REF]) that:

Λ ge = Λ V . (1.27)
For the Schrödinger equation, Alessandrini [START_REF] Alessandrini | Stable determination of conductivity by boundary measurements[END_REF] has obtained the following stability estimate : assume that, for j = 1, 2, the potentials V j belong to the Sobolev space W m,1 (M ) with ||V j || W m,1 (M ) ≤ A. Then there exists a constant C A > 0 such that

||V 1 -V 2 || L ∞ (M ) ≤ C A log(3 + |||Λ V1 -Λ V2 ||| -1 -α , (1.28)
with α = 1 -d m and where |||B||| is the norm of the operator B : L ∞ (∂M ) → L ∞ (∂M ). This stability estimate was then improved by Novikov in [START_REF] Novikov | New global stability estimates for the Gelfand-Calderón inverse problem[END_REF], where he proved one can take α = m -d.

In Theorem 1.5, under some additional assumptions on the derivatives of the warping function at r = 0, we obtain a rather similar stability estimate with, roughly speaking, α replaced by p -1. Note that our potential V is only C m-2 and that p ≤ m -1. We emphasize that we do not assume that c(1) = 1 and c (1) = 0, thus the connection between Λ g and Λ V is not so clear and depends implicitly on the unknown values c(1) and c (1) which we are trying to estimate.

In the singular case, we also obtain a log-type estimate result, but which is much less precise. Indeed, in this case, we do not have explicit formula for the angular eigenvalues κ k , so we need to use the Weyl's law to reasonably approximate the κ k . We can prove the following theorem: Theorem 1.6. Let (M, g) and (M, g) be singular Riemannian manifolds given by (1.1)-(1.2) with c, c ∈ C(A) where A > 0 is fixed. Assume that for > 0 small enough, one has:

sup k≥0 |σ k -σk | ≤ . (1.29)
Then, there exists θ ∈ (0, 1) and a positive constant C A , depending only on A such that,

||c -c|| L ∞ (0,1) ≤ C A 1 log( 1 ) (p-1)θ (1.30)
Remark 1.7. In the definition of our admissible warping functions set C(A), we impose that the integer p ≥ 2. In particular, the first derivative at r = 0 of the warping function c(r) must vanish. When the metric g is smooth, this condition is automatically satisfied. When the metric g is singular, we only make this assumption for simplicity of exposition in order to use the same strategy as for the regular case. This technical assumption could be certainly relaxed with a little effort, and one would obtain a logarithmic stabilily estimate close to the one stated in Theorem 1.6.

The rest of our paper is organized as follows. In Section 2, using the separation of variables, we recall the standard construction of the Dirichlet to Neumann map (see [START_REF] Daudé | Non uniqueness results in the anisotropic Calderón problem with Dirichlet and Neumann data measured on disjoint sets[END_REF], [START_REF] Daudé | On the hidden mechanism behind non-uniqueness for the anisotropic Calderón problem with data on disjoint sets[END_REF], [START_REF] Daudé | The anisotropic Calderón problem for singular metric of warped product type : the borderline between uniqueness and invisibility[END_REF] for details), and some basics facts on the Steklov spectrum. In Section 3, we prove that the warping function c(r) is uniquely determined by the knowledge of the Steklov spectrum. We also obtain a local uniqueness result from the approximate knowledge of the Steklov spectrum. In Section 4, we prove our stability results: the proofs are based on the Müntz-Jacskon approximation, starting from the knowledge of the Hausdorff moments.

The Steklov spectrum

In this section, we construct explicitly the Steklov spectrum by using the warped product structure of the manifold (M, g). More precisely, we use the underlying symmetry of the warped product in order to diagonalize the DN map onto the Hilbert basis of harmonics {Y k } k≥0 , i.e. the normalized eigenfunctions, ofg S . On each harmonic, the DN map acts as an operator of multiplication by essentially the Weyl-Titchmarsh function associated to the countable family of Schrödinger operators arising from the separation of variables procedure. The Weyl-Titchmarsh theory will then allow us to prove the asymptotic of theorem 1.1 as well the (local) uniqueness results of Theorems 1.3 and 1.4 in Section 3.

Let us first solve the Dirichlet problem (1.5). In the coordinate system (x, ω), the Laplace equation

-g u = 0 reads [-∂ 2 x -g S + q f (x)]v = - (d -2) 2 4 v, (2.1)
where v = f d-2 u and q f is given by (1.9). Observe that under the hypothesis (1.4), the effective potential q f is a smooth function on [0, +∞) that satisfies the asymptotics

q f (x) = O(e -px ), x → ∞. (2.2)
We now use the warped structure to separate variables. We thus look for solutions of (2.1) of the form

v = ∞ k=0 v k (x)Y k . (2.3)
Hence, for each k ≥ 0, the functions v k satisfy the Schrödinger equation on the half-line

-v k + q f (x)v k = -κ 2 k v k , x ∈ [0, ∞), (2.4)
where κ k is given by (1.12) for all k ≥ 0. Actually, it is very useful to consider complex spectral parameter z ∈ C and we are interested by some special solutions of the Sturm-Liouville equation:

-v + q f (x)v = zv. (2.5)
We denote by {C 0 (x, z), S 0 (x, z)} the fundamental system of solutions of (2.5) with a spectral parameter z ∈ C that satisfy Neumann and Dirichlet conditions at x = 0 respectively, given by C 0 (0, z) = 1, C 0 (0, z) = 0, S 0 (0, z) = 0, S 0 (0, z) = 1.

(2.6)

Note that W (C 0 (x, z), S 0 (x, z)) = 1, (2.7)
where the Wronskian is defined by W (u, v) = uv -u v. Moreover, the functions z → C 0 (x, z), S 0 (x, z) are entire in z. Note also that under the hypothesis (2.2), the Schrödinger operator H = -d 2 dx 2 + q f is in the limit point case at x = ∞. In consequence, for all z ∈ C, there exists a unique (up to constant factor) solution S ∞ (x, z) of (2.5) that is L 2 in a neighbourhood of x = ∞, (see [START_REF] Reed | Methods of modern mathematical physics-Scattering theory[END_REF], Theorem XI.57 where our spectral parameter z = k 2 ). We write this function as

S ∞ (x, z) = A(z) (C 0 (x, z) -M (z)S 0 (x, z)) .
(2.8)

Using (2.7), we thus get the following expressions for the function A(z) and the Weyl-Titchmarsh function M (z)

A(z) = W (S ∞ (x, z), S 0 (x, z)), M (z) = - W (C 0 (x, z), S ∞ (x, z)) W (S 0 (x, z), S ∞ (x, z)) = S ∞ (0, z) S ∞ (0, z) . (2.9)
For later use, we also introduce the characteristic functions

∆(z) = W (S 0 (x, z), S ∞ (x, z)) = -A(z), d(z) = W (C 0 (x, z), S ∞ (x, z)).
(2.10)

We thus have:

M (z) = - d(z) ∆(z) . (2.11)
Notice that the characteristic functions ∆(z) and d(z) are analytic on C, and the WT function M is analytic on C\[β, +∞[ with -β sufficiently large. The zeros (-α 2 j ) j≥0 of the function z → ∆(z) are precisely the Dirichlet eigenvalues of the self-adjoint operator H (and thus are real) whereas the zeros (-γ 2 j ) j≥0 of the function z → d(z) are the Neumann eigenvalues of the self-adjont operator H (and thus are real too). Moreover, we know that σ ess (H) = [0, +∞) and that the essential spectrum contains no embedded eigenvalues ( [START_REF] Reed | Methods of modern mathematical physics-Analysis of operator[END_REF], Thm XIII.56). In consequence of the spectral theorem, the eigenvalues (-α 2 j ) j≥0 and (-γ 2 j ) j≥0 must then satisfy

-α 2 j , -γ 2 j ∈ [min(q f ), 0] , ∀j ≥ 0, (2.12)
with the usual convention that [min(q f ), 0] = ∅ if min(q f ) > 0.

We have a more precise result on the location of the eigenvalues -α 2 j of H thanks to

Lemma 2.1. The discrete Dirichlet spectrum of H = -d 2 dx 2 + q f is finite and contained in (-(d-2) 2 4 , 0). Proof. Introduce the operators L = H + (d-2) 2 4
and K = f -d+2 Lf d-2 . An easy calculation shows that

K = - 1 f d-2 d dx f d-2 d dx .
We remark that K is selfadjoint on L 2 (R + ; f 2d-4 (x)dx) if we impose Dirichlet boundary condition at x = 0. Moreover, the operator K is clearly positive. Hence we have σ pp (K) ⊂ (0, +∞). But λ is an eigenvalue of K if and only if λ

-(d-2) 2 4
is an eigenvalue of H. In particular, we obtain that an eigenvalue

-α 2 j of H always satisfy -α 2 j > -(d-2) 2 4
. Together with (2.12), this proves the result.

Corollary 2.2. 0 does not belong to the Dirichlet spectrum of -∆ g .

Proof. Assume the converse. Then there exists u = 0 such that -∆ g u = 0 on M and u = 0 on ∂M . Using separation of variables, this means that there exists k ≥ 0, v k = 0 and v k (0) = 0 such that Hv k = -κ 2 k v k . In other words, -κ 2 k is an eigenvalue of H with Dirichlet boundary condition. According to Lemma 2.1, we then must have -κ

2 k > -(d-2) 2 4 . But since, µ k ≥ 0, we always have -κ 2 k ≤ -(d-2) 2 4 . Contradiction.
Finally, the DN map can be clearly diagonalized onto the Hilbert basis of harmonics {Y k } k≥0 . If we represent the Dirichlet data as ψ = k≥0 ψ k Y k , then the global DN map has the expression

Λ g ψ = ∞ k=0 (Λ k g ψ k )Y k , (2.13)
where the diagonalized DN map are defined by

Λ k g ψ k = (d -2)f (0) f d+1 (0) v k (0) - v k (0) f d (0) . (2.14) Using that v k (x) = α k C 0 (x, κ k ) + β k S 0 (x, κ k )
, a straightforward calculation shows that the partial DN map Λ k g acts as an operator of multiplication, precisely

Λ k g ψ k = (d -2)f (0) f 3 (0) - M (-κ 2 k ) f 2 (0) ψ k . (2.15)
We see immediately from (2.15) that the partial DN map Λ k g acts essentially by an operator of multiplication by the WT function M (-κ 2 k ) associated to (2.4) up to some boundary values of the warping function f ans its first derivative f . In consequence, we infer that the Steklov spectrum of (M, g), that is the set of eigenvalues of Λ g , is precisely given by

{σ k , k ≥ 0} = (d -2)f (0) f 3 (0) - M (-κ 2 j ) f 2 (0)
, j ≥ 0 .

(2.16)

In fact, we have the following exact identification:

Lemma 2.3. 1. The function y ∈ R → M (y) is strictly increasing on R \ {-α 2 j , j ≥ 0}. 2. Then ∀k ≥ 0, σ k = (d-2)f (0) f 3 (0) - M (-κ 2 k ) f 2 (0) .
Proof. 1. Let y ∈ R\{-α j 2 }. Notice then that C 0 (x, y), S 0 (x, y), S ∞ (x, y) are real and thus ∆(y) and d(y) where ˙denotes the derivative with respect to y. We conclude from this that 1. holds. Observe that between two Dirichlet eigenvalues -α 2 j < -α 2 j+1 , the WT function M (y) is strictly increasing and goes from from -∞ to +∞.

are real too. Let y * ∈ R. Since (S ∞ (x, y)S ∞ (x, y * ) -S ∞ (x, y)S ∞ (x, y * )) = (y -y * )S ∞ (x,
2. From Lemma 2.1, we know that the Dirichlet eigenvalues of H satisfy

-(d-2) 2 4
< -α 2 j ≤ 0 for all j ≥ 0. In other words, we have

∀j ≥ 0, 0 ≤ α j < d -2 2 . Since κ k ≥ d-2
2 for all k ≥ 0 (see (1.12)), we deduce from 1. that the function κ k → M (-κ 2 k ) is strictly decreasing for k ≥ 0. Hence 2. follows from (2.16) and the ordering of the Steklov spectrum (σ k ) k≥0 .

Having obtained the exact expression of the Steklov spectrum in terms of the WT functions (2.9) evaluated at the -κ 2 k , we can prove easily Theorem 1.1 by using the well-known asymptotic of a WT function due to Danielyan and Levitan [START_REF] Danielyan | Asymptotic behavior of the Weyl-titchmarsh m-function[END_REF]. We refer for instance to [START_REF] Simon | A new approach to inverse spectral theory, I. Fundamental formalism[END_REF], Section 4 and the references therein for a proof and more asymptotic results.

Proof of Theorem 1.1. Recall the asymptotic of the WT function from [START_REF] Simon | A new approach to inverse spectral theory, I. Fundamental formalism[END_REF], Thm 4.5. If q f ∈ C N ([0, δ)) with δ > 0, then as κ → ∞, we have

M (-κ 2 ) = -κ - N j=0 β j (0)κ -j-1 + O(κ -N -2 ), (2.17)
where the constants β j (0) can be calculated inductively by (1.15). Since the potential q f is assumed to be smooth, the asymptotic (2.17) together with lemma 2.3 lead to the result.

Uniqueness and local uniqueness

In this section, we prove local uniqueness, that is Theorem 1.4, under the assumption that the transversal Riemanniann manifold (S n-1 , g S ) is known. The general idea of local uniqueness inverse results arises in the different versions of the local Borg-Marchenko Theorem stated first by Simon in [START_REF] Simon | A new approach to inverse spectral theory, I. Fundamental formalism[END_REF] and proved differently or extended to singular settings in [START_REF] Bennewitz | A proof of the local Borg-Marchenko Theorem[END_REF][START_REF] Eckhardt | Uniqueness results for Schrödinger operators on the line with purely discrete spectra[END_REF][START_REF] Gesztesy | A new approach of inverse spectral theory, II. General potentials and the connection to the spectral measure[END_REF][START_REF] Gesztesy | On local Borg-Marchenko uniqueness results[END_REF][START_REF] Kostenko | Weyl-Titchmarsh theory for Schrödinger operators with strongly singular potentials[END_REF]. We shall follow here the initial version due to Simon that makes an intensive use of a representation of the WT function M as the Laplace transform of what Simon called the A-function. Precisely, Simon showed in [START_REF] Simon | A new approach to inverse spectral theory, I. Fundamental formalism[END_REF], Thm 2.1, that there exists a function A on [0, ∞) with A -q f continuous, obeying

|A(α) -q f (α)| ≤ Q(α) 2 e αQ(α) , Q(α) = α 0 |q f (s)|ds, (3.1)
such that, if κ > 1 2 q f L 1 , then M (-κ 2 ) = -κ - ∞ 0 A(α)e -2κα dα. (3.2)
We also have

|A(α, q f ) -A(α, qf )| ≤ q f -qf L 1 [Q(α) + Q(α)]e α[Q(α)+ Q(α)] . (3.3)
Finally, Simon also proved the local uniqueness result Theorem 3.1 ([29], Thm 1.5). The potential q f on [0, a] is a function of A on [0, a]. Explicitly, if q f and qf are two potentials, let A and à be their A-functions. Then

A(α) = Ã(α), ∀α ∈ [0, a] ⇐⇒ q f (x) = qf (x), ∀x ∈ [0, a].
We shall use these results as follows.

Proof of theorem 1.4. Suppose that the assumption (1.18) is satisfied. Using Theorem 1.1 and Lemma 2.3, we see immediately that

f (0) = f (0), f (0) = f (0). (3.4)
Hence the assumption (1.18) can be equivalently read as

M (-κ 2 k d-1 ) -M (-κ 2 k d-1 ) = O(e -2aκ k d-1 ), k → ∞. (3.5)
We use then the representation of the WT functions (3.2), the behavior (2.2) of the potentials q f , q f and the estimate (3.3) to show that (3.5) entails

a 0 A(α) -Ã(α) e -2κ k d-1 α dα = O(e -2aκ k d-1 ), k → ∞. (3.6)
Now, we need the following proposition which is a slight generalization to the case of noninteger (κ k ) k≥0 of Proposition 2.4 in [START_REF] Horvath | Partial identification of the potential from phase shifts[END_REF].

Proposition 3.2. Let f ∈ L 1 (0, a). Assume that a 0 e -κ k d-1 t f (t) dt = O(e -aκ k d-1 ) , k → +∞. (3.7)
Then, f = 0 almost everywhere.

Proof. Thus, Proposition 3.2 implies in particular that

Setting λ k = 1 c d-1 κ k d-
A(α) = Ã(α), ∀α ∈ [0, a],
from which we infer using Theorem 3.1 that

q f (x) = q f (x), ∀x ∈ [0, a].
Recalling the definition (1.9) of q f , we thus see that

(f d-2 ) (x) = ( f d-2 ) (x) f d-2 (x) f d-2 (x), ∀x ∈ [0, a] (3.12)
We finish the proof seeing that (3.12) can be viewed as a linear second-order ODE for f d-2 . Recalling from (3.4) that the Cauchy data (3.4) of f d-2 are equal to those of f d-2 , we conclude that the unique solution of (3.12) on [0, a] is f d-2 = f d-2 . Whence the asserted result.

Conversely, assume that (1.19) holds. In particular, q f (x) = q f (x) for all x ∈ (0, a). So, using Theorem 3.1, we get A(α) = Ã(α) for all α ∈ (0, a). Then, using the same arguments as in the first part of the proof, we obtain M (-κ 2 ) = M (-κ 2 ) + O(e -2aκ ) , κ → +∞.

(3.13)

Then, using Lemma 2.3 and noting that f (0), f (0) are also known, we get immediately (1.18).

4 Proof of Theorems 1.5 and 1.6

A Volterra type integral operator

Let us begin by an elementary lemma:

Lemma 4.1. For any κ large enough, we have:

S ∞ (0, -κ 2 ) S∞ (0, -κ 2 ) M (-κ 2 ) -M (-κ 2 ) = +∞ 0 (q f (x) -q f (x)) S ∞ (x, -κ 2 ) S∞ (x, -κ 2 ) dx.
Proof. To simplify the notation, we set z = -κ 2 and we integrate over the interval (0, ∞) the obvious equality:

S ∞ (x, z) S ∞ (x, z) -S ∞ (x, z) S∞ (x, z) = (q f (x) -q f (x))S ∞ (x, z) S∞ (x, z). (4.1)
For k = 0 or 1, S

∞ (x, z) → 0 as x → +∞, (see for instance [START_REF] Berezin | The Schrödinger equation[END_REF]), so we get immediately

S ∞ (0, z) S∞ (0, z) -S ∞ (0, z) S ∞ (0, z) = +∞ 0 (q f (x) -q f (x))S ∞ (x, z) S∞ (x, z) dx, (4.2)
which implies the Lemma, thanks to (2.9).

In the following lemma, we express the Weyl solution S ∞ (x, z) with the help of the well-known Marchenko's representation, (we refer to [START_REF] Marchenko | Sturm-Liouville operators and their applications[END_REF], Chapter III for details): Lemma 4.2. Assume that c ∈ C(A). Then, there exists a C m-1 function K(x, t) for 0 ≤ x ≤ t < ∞, satisfying the properties:

S ∞ (x, -κ 2 ) = e -κx + +∞ x K(x, t)e -κt dt , κ > 0. (4.3) K(x, x) = 1 2 +∞ x q f (t) dt. (4.4)
Moreover, there exists a constant C A > 0 depending only on A such that,

|∂ k x ∂ l t K(x, t)| ≤ C A e -p 2 (x+t) , ∀k, l ≤ m -1. (4.5)
Proof. The existence of the Marchenko's kernel K(x, t) is proved in ( [START_REF] Marchenko | Sturm-Liouville operators and their applications[END_REF], Lemma 3.1.1), and we have the following estimate:

|K(x, t)| ≤ 1 2 σ( x + t 2 ) e σ1(x)-σ1( x+t 2 ) , (4.6)
where

σ(x) = +∞ x |q f (s)| ds , σ 1 (x) = +∞ x σ (
= l = 0. Now, let us define H(u, v) = K(u -v, u + v) for 0 ≤ v ≤ u. Thanks to ([22], Lemma 3.1.2), H obeys: ∂H ∂u (u, v) = - 1 2 q f (u) - v 0 q f (u -s)H(u, s) ds, (4.8)
∂H ∂v (u, v) = +∞ u q f (s -v)H(s, v) ds. (4.9)
Then, (4.5) follows from a straightforward calculation.

As a by-product, we get:

Corollary 4.3. Let c ∈ C(A) be a warping function for the metric (1.2). Then, there exists a constant

C A such that, |S ∞ (0, -κ 2 ) -1| ≤ C A κ + 1
, for all κ ≥ 0.

Proof. Using (4.5) for k = p = 0, we get:

S ∞ (0, -κ 2 ) -1 = +∞ 0 K(0, t) e -κt dt ≤ C A κ + 1 .
Now, let us introduce a new kernel K 1 (x, t) for 0 ≤ t ≤ x < ∞, by the formula:

K 1 (x, t) = 2K(t, 2x -t) + 2 K(t, 2x -t) + 2 2x-t t K(t, u) K(t, 2x -u) du, (4.10)
where K(x, t) is the Marchenko's kernel associated with the potential q f , (see Lemma 4.2). We have the following estimate which follows immediately from (4.5):

Lemma 4.4. Assume that c ∈ C(A). Then, for all α < p, there exists a constant C A,α > 0 depending only on A and α such that,

|∂ k x ∂ l t K 1 (x, t)| ≤ C A,α e -αx , ∀k, l ≤ m -1. (4.11)
Finally, we consider the corresponding Volterra type integral operator B given by:

Bh(x) = h(x) + x 0 K 1 (x, t)h(t) dt. (4.12)
This operator B is crucial to our analysis because it links the Steklov spectrum to the difference of the potentials q f and q f . More precisely, one has the following result:

Lemma 4.5. For κ sufficiently large, we have:

S ∞ (0, -κ 2 ) S∞ (0, -κ 2 ) M (-κ 2 ) -M (-κ 2 ) = +∞ 0 e -2κx B[q f -q f ](x) dx. (4.13)
Proof. Thanks to ([18], Lemma 2.5), we have:

+∞ 0 (q f (x) -q f (x))S ∞ (x, -κ 2 ) S∞ (x, -κ 2 ) dx = +∞ 0 e -2κx B[q f -q f ](x) dx. (4.14)
Then, using Lemma 4.1 and (4.14), we get (4.13).

This Volterra operator B also possesses good properties on some L 2 -spaces equipped with exponential weights, and which are defined by:

H δ = {q : ||q|| 2 H δ := +∞ 0 |q(x)| 2 e δx dx < ∞}. (4.15)
If c ∈ C(A), then it results from (1.21) that q f ∈ H δ for any δ < 2p. Moreover, there exists a constant C A,δ depending only on A and δ such that

||q f || H δ ≤ C A,δ for all c ∈ C(A). (4.16)
In the following Proposition, which is close to [START_REF] Horvath | Inverse scattering with fixed energy and an inverse eigenvalue problem on the half-line[END_REF], Lemmas 2.5 -2.6, we give a uniform estimate on the norm of B : H δ → H δ , and its inverse, when the warping functions belong to the admissible set C(A). This result will be very useful to estimate the difference of the potentials q f -q f , (for the topology of H δ ).

Proposition 4.6. Let c, c be warping functions belonging to C(A). Then, for any 0 < δ < p, we have:

B : H δ → H δ is an isomorphism,
and there exists a constant C A,δ depending only on A and δ such that

||B|| + ||B -1 || ≤ C A,δ . ( 4
K n+1 (x, t) = x t K 1 (x, s)K n (s, t) ds , t ≤ x, (4.20)
and we have the following estimates which can be easily proved by induction: for δ < 2,

|K n (x, t)| ≤ C A,δ e -δx n (x -t) n-1 (n -1)! , t ≤ x. (4.21)
For n ≥ 2, using the rough bound ||C n || ≤ ||C n || HS , we get:

||C n || ≤ +∞ 0 x 0 C A,δ e -δx 2n (x -t) n-1 (n -1)! 2 e δ(t+x) dt dx ≤ (C A,δ ) n (n -1)! 2 +∞ 0 x 2n-1 e -2(n-1)δx dx ≤ (C A,δ ) n (n -1)! 2 1 2(n -1)δ 2n +∞ 0 y 2n-1 e -y dy ≤ C A,δ 2(n -1)δ 2n Γ(2n) Γ(n) 2 .
So, thanks to Stirling's formula, we get

∞ n=0 ||C n || ≤ 1 + ||C|| + +∞ n=2 C A,δ 2(n -1)δ 2n Γ(2n) Γ(n) 2 ≤ C A,δ . (4.22)
This means that the Neumann series

B -1 = +∞ n=0 (-1) n C n (4.23)
is convergent in the operator norm and

||B -1 || ≤ C A,δ .
Now, let us assume again that the warping functions c, c ∈ C(A) and that for all k ≥ 0,

|σ k -σk | ≤ . (4.24)
By making k → +∞ in (4.24), and thanks to Theorem 1.1, we deduce that f (0) = f (0), and also

(d -2)f (0) f 3 (0) - (d -2) f (0) f 3 (0) ≤ . (4.25)
Thus, recalling that for sufficiently large k, we have

σ k = (d -2)f (0) f 3 (0) - M (-κ 2 k ) f 2 (0) , (4.26)
we obtain:

| M (-κ 2 ) -M (-κ 2 )| ≤ 2f 2 (0) ≤ 2A 2 , (4.27)
since c ∈ C(A). So, plugging (4.27) into (4.13) and using Corollary 4.3, we easily get the following result:

Lemma 4.7. Let c, c ∈ C(A).
Assume that for all k ≥ 0, |σ k -σk | ≤ . Then, there exists a positive constant C A which does not depend on such that

+∞ 0 e -2κ k x B[q f -q f ](x) dx ≤ C A . (4.28)
Note that this kind of estimates fits into the so-called moment theory, (see for instance [START_REF] Ang | Moment theory and some inverse problems in potential theory and heat conduction[END_REF] for a nice exposition and references therein), and this is the object of the next sections.

A Müntz-Jackson's theorem.

Let Λ ∞ = {0 ≤ λ 0 < λ 1 < ... < λ n < ...}, λ n → +∞ be a sequence of positive real numbers. The classical Müntz-Szász's Theorem characterizes the sequences (λ k ) for which all functions in C 0 ([0, 1]), or in L 2 ([0, 1]), can be approximated by "Müntz polynomials" of the form:

P (x) = n k=0 a k x λ k , (4.29)
with real coefficients a k . More precisely, one has:

Theorem 4.8. Let Λ ∞ be a sequence of positive real numbers as above. Then, span

{x λ0 , x λ1 , ...} is dense in L 2 ([0, 1]) if and only if +∞ k=1 1 λ k = ∞. (4.30)
Moreover, if λ 0 = 0, the denseness of the Müntz polynomials in C 0 ([0, 1]) in the sup norm is also characterized by (4.30).

Now, for n ≥ 1, let us consider the finite sequence

Λ := Λ n : 0 ≤ λ 0 < λ 1 < ... < λ n . (4.31)
We define the subspace of the "Müntz polynomials of degree λ n " as:

M(Λ) = {P : P (x) = n k=0 a k x λ k }. (4.32)
The error of approximation from M(Λ) of a function

f in C 0 ([0, 1]) or in L 2 ([0, 1]
) is given by:

E(f, Λ) p := inf P ∈M(Λ) ||f -P || p = ||f -P 0 || p , (4.33)
for some P 0 ∈ M(Λ) depending on whether p = 2 or p = ∞. Clearly, one has

E(f, Λ) 2 ≤ E(f, Λ) ∞ and E(f, Λ) 2 = ||f -π n (f )|| 2 where π n (f )
is the orthogonal projection of f on the subspace M(Λ). An estimation from above of E(f, Λ) p in terms of the smoothness of f is called a Müntz-Jackson's theorem.

To estimate this error of approximation for the uniform norm, let us consider the so-called Blaschke product B(z), z ∈ C,

B(z) := B(z, Λ) = n k=0 z -λ k z + λ k . (4.34)
The index of approximation of Λ in C 0 ([0, 1]) is defined as :

∞ (Λ) = max y≥0 B(1 + iy) 1 + iy . (4.35)
Its relevance is justified by the following result ([21], Theorem 2.6, Chapter 11) when λ 0 = 0: Proposition 4.9. Let Λ : 0 = λ 0 < λ 1 < ... < λ n be a finite sequence. Then, for each

f ∈ C 1 ([0, 1]), E(f, Λ) ∞ ≤ 20 ∞ (Λ) ||f || ∞ . (4.36)
As a by-product, we can easily prove:

Corollary 4.10. Let Λ * : 0 < λ 1 < ... < λ n be a finite sequence. Then, for each f ∈ C 1 ([0, 1]) with f (0) = 0, one has:

E(f, Λ * ) ∞ ≤ 40 ∞ (Λ * ) ||f || ∞ . (4.37)
Proof. Consider the finite sequence Λ : 0 = λ 0 < λ 1 < ... < λ n . Thanks to Proposition 4.9, there exists

P ∈ M(Λ) such that E(f, Λ) ∞ = ||f -P || ∞ ≤ 20 ∞ (Λ) ||f || ∞ . (4.38) Since f (0) = 0, we deduce that |P (0)| ≤ 20 ∞ (Λ) ||f || ∞ . Thus, setting Q = P -P (0) ∈ M(Λ * ), we obtain: ||f -Q|| ∞ ≤ 40 ∞ (Λ) ||f || ∞ , (4.39)
which concludes the proof since ∞ (Λ) = ∞ (Λ * ).

In the same way, we can recursively approximate C r -differentiable functions, (r ≥ 1). To this end, let Λ * : r -1 < λ 1 < ... < λ n be a finite sequence. For 0 ≤ k ≤ r -1, we set:

Λ (k) : λ (k) 1 = λ 1 -k, ... , λ (k) n = λ n -k, (4.40)
and the indices of approximation

(k) ∞ = ∞ (Λ (k)
), k = 0, 2, ..., r -1. We have the following result: Corollary 4.11. For each f ∈ C r ([0, 1]) such that f (k) (0) = 0 for all k = 0, ..., r -1, one has:

E(f, Λ * ) ∞ ≤ 40 r r-1 k=0 (k) ∞ ||f (r) || ∞ . (4.41)
Proof. For r = 1, this estimate is nothing but Corollary 4.10. Now, assume that r = 2 and consider f ∈ C 2 ([0, 1]) with f (0) = f (0) = 0. Using Corollary 4.10 for the function f and the finite sequence Λ (1) , we see there exists 1) ) such that

P (x) = n k=1 a k x λ k -1 ∈ M(Λ (
E(f , Λ (1) ) ∞ = ||f -P || ∞ ≤ 40 (1) ∞ ||f || ∞ . (4.42)
We set

F (x) = f (x) - x 0 P (t) dt = f (x) - n k=1 a k λ k x λ k .
Since F (0) = 0, using again Corollary 4.10 for the finite sequence Λ * , we see there exists Q ∈ M(Λ * ) such that

E(F, Λ * ) ∞ = ||F -Q|| ∞ ≤ 40 ∞ ||F || ∞ , (4.43)
or equivalently

||f -(

x 0 P (t) dt + Q)|| ≤ 40 ∞ ||f -P || ∞ . (4.44)
Observing that

x 0 P (t) dt + Q(x) = n k=1 a k λ k x λ k + Q(x) ∈ M(Λ * )
, and using (4.42), we obtain:

E(f, Λ * ) ∞ ≤ ||f -( x 0 P (t) dt + Q)|| ≤ 40 2 ∞ (1) ∞ ||f || ∞ , (4.45)
which proves Corollary 4.11 in the case r = 2. For r ≥ 3, the proof is identical.

For special finite sequences Λ, the index of approximation ∞ (Λ) can be replaced by a much simpler expression. For instance, we have the following result, ( [START_REF] Lorentz | Makovoz Y Constructive approximation, A series of comprehensive studies in Mathematics[END_REF], Theorem 4.1, Chapter 11): Theorem 4.12. Let Λ : 0 = λ 0 < λ 1 < λ 2 < ... < λ n be a finite sequence. Assume that λ k+1 -λ k ≥ 2 for k ≥ 0. Then,

∞ (Λ) = |B(1, Λ)| = n k=1 λ k -1 λ k + 1 (4.46)
In particular, if λ k = 2k + b for k = 1, ..., n where b > 0, one has As an application of the previous section, we shall give an approximation of the L 2 -norm of a function f from the approximately knowledge of a finite number of these moments:

∞ (Λ) = b + 1 2n + b + 1 . ( 4
m k = 1 0 t λ k f (t) dt , k ∈ X ⊂ N. (4.48) Thanks to Theorem 4.8, if Λ ∞ = {0 ≤ λ 0 < λ 1 < ... < λ n < ...}, λ n → +∞ is a sequence of positive real numbers such that +∞ k=1 1 λ k = ∞, (4.49)
the system {x λ0 , x λ1 , ...} is complete in L 2 ([0, 1]). Thus, the knowledge of the complete sequence (m k ) k≥0 , uniquely determines the function f . But, in practice, one has available only a finite set m 0 , ..., m n of moments, and furthermore these moments are usually corrupted with noise. It is well-known that this problem is severely ill-posed, (see for instance [START_REF] Ang | Moment theory and some inverse problems in potential theory and heat conduction[END_REF] and references therein, for further details).

Let us briefly explain the approach given in [START_REF] Ang | Moment theory and some inverse problems in potential theory and heat conduction[END_REF]. Using the Gram-Schmidt process, we define the polynomials (L m (x)) as L 0 (x) = 1, and for m ≥ 1,

L m (x) = m j=0 C mj x λj , (4.50)
where we have set

C mj = 2λ m + 1 m-1 r=0 (λ j + λ r + 1) m r=0,r =j (λ j -λ r ) . ( 4

.51)

The family (L m (x)) defines an orthonormal Hilbert basis of L 2 ([0, 1]). For instance, if λ k = k, the polynomials (L m (x)) are the Legendre polynomials, and in this case, the latter coefficients C mj are given by:

C 0 mj := √ 2m + 1 (-1) m-j (m + j)! (m -j)! j! 2 (4.52)
Note that, the generalized binomial theorem gives easily the upper bound: We denote π n (f ) the orthogonal projection on on the subspace M(Λ):

|C 0 mj | ≤ √ 2m + 1 3 m+j . ( 4
π n (f ) = n k=0 < f, L k > L k . (4.56)
Thus, we deduce that:

||π n (f )|| 2 2 = n k=0 | < f, k p=0 C kp x λp > | 2 ≤ 2 n k=0 k p=0 |C kp | 2 , (4.57)
thanks to our hypothesis on the moments (4.55). So, we get immediately:

||f || 2 2 = ||π n (f )|| 2 2 + ||f -π n (f )|| 2 2 = ||π n (f )|| 2 2 + E(f, Λ) 2 2 ≤ 2 n k=0 k p=0 |C kp | 2 + E(f, Λ) 2 ∞ ≤ 2 n k=0 k p=0 |C kp | 2 + E(f, Λ * ) 2 ∞ ,
where Λ * : 0 < λ 1 < ... < λ n . In particular, if f ∈ C r ([0, 1] with f (k) (0) = 0 for k = 0, ..., r -1, and if r -1 < λ 1 , we get using Corollary 4.11:

||f || 2 2 ≤ 2 n k=0 k p=0 |C kp | 2 + 40 r r-1 k=0 (k) ∞ ||f (r) || ∞ 2 . (4.58)
At this stage, it is important to make the following remark: on the one hand, the double sum appearing in the (RHS) of (4.58) can be very large with respect to n. Indeed, in the case λ k = k, the coefficients C 0 k0 are equal to √ 2k + 1 thanks to (4.52). On the other hand, if λ k = 2k, Theorem 4.12 suggests that the second term in the (RHS) of (4.58) is equal to O(n -r ), n → +∞. Thus, if we want to control reasonably (with respect to ) the L 2 -norm of the function f , we have to choose a suitable n = n( ) in the equation (4.58). Of course, this choice will depend heavily of the behaviour of the coefficients C kp . This will be done in the next two sections where we separate the simpler case where the metric in regular, and the case where the metric is singular at r = 0.

4.4 The regular case.

Stability estimates for the potentials.

In this case, we recall that necessarily, the metric g S = dΩ 2 and κ

k = k + d-2
2 , (we order the Steklov spectrum without counting multiplicity). So, making the change of variables t = e -x in (4.28), we obtain:

| 1 0 t 2k+d-3 B[q f -q f ](-log t) dt| ≤ C A , ∀k ≥ 0.
(4.59) Thus, introducing δ ∈]0, 1[, we get:

| 1 0 t λ k h(t) dt| ≤ C A , ∀k ≥ 0, (4.60)
where we have set h

(t) = t -δ+1 2 B[q f -q f ](-log t) and λ k = 2k + d -3 + δ+1 2 . Note that ||h|| L 2 (0,1) = ||B[q f -q f ]|| H δ . (4.61)
Using Lemma 4.4, we see that h(t) ∈ C m-2 ((0, 1]) and we have the following result:

Lemma 4.13. The function t → h(t) extends in a C p-1 -differentiable function on [0, 1] with h (k) (0) = 0 for k = 0, ..., p -1. Moreover, there exists a constant C A > 0 such that, for all k ≤ p -1,

||h (k) || ∞ ≤ C A . ( 4

.62)

Proof. We only sketch the proof since the arguments are straightforward. Recalling that B = Id + C, we write :

h(t) = t -δ+1 2 [q f -q f ](-log t) + t -δ+1 2 C[q f -q f ](-log t) := h 1 (t) + h 2 (t).
First, let us estimate the derivatives of h 1 (t). Setting Q(x) := [q f -q f ](x), we deduce immediately from (1.21) that:

|Q (k) (x)| ≤ C A e -px , ∀x ≥ 0, ∀k = 0, ..., m -2. (4.63)
Then, we get easily:

|h (k) 1 (t)| ≤ C A t p-k-δ+1 2 , ∀k = 0, ..., m -2, ∀t ∈]0, 1]. (4.64)
In the same way, using again Lemma 4.4, a tedious calculation shows that for all α < p, there exists a constant C A,α > 0 such that

|h (k) 2 (t)| ≤ C A,α t α-k-δ+1 2 , ∀k = 0, ..., m -2. ( 4
C mj = √ 4m + 2b + 1 m-1 r=0 (2j + 2r + 4b + 1) m r=0,r =j (2j -2r) . ( 4

.66)

Setting M = max(2, 4b + 1), we get easily:

|C mj | ≤ 4m + M + 1 2 M 2 m m-1 r=0 (j + r + 1) m r=0,r =j (j -r) , (4.67)
or equivalently

|C mj | ≤ 4m + M +1 2 2m + 1 M 2 m |C 0 mj |. (4.68)
Thus, there exists a universal constant B > 0 such that

|C mj | ≤ B M 2 m |C 0 mj |. (4.69)
Thanks to (4.69), we can estimate the L 2 -norm of the ortogonal projection π n (h) on the subspace M(Λ) where Λ = Λ n : 0 < λ 0 < λ 1 < ... < λ n , where

λ k = 2k + d -3 + δ+1 2 : ||π n (h)|| 2 2 ≤ 2 n k=0 k p=0 |C kp | 2 ≤ 2 n k=0 k p=0 M 2 k |C 0 kp | 2 ≤ B 2 2 n k=0 M 2 k k p=0 |C 0 kp | 2 .
Using (4.53), we see that

k p=0 |C 0 kp | ≤ √ 2k + 1 k p=0 3 k+p ≤ 3 2 √ 2k + 1 3 2k . (4.70)
We deduce that

n k=0 M 2 k k p=0 |C 0 kp | 2 ≤ 9 4 n k=0 M 2 2k (2k + 1)3 4k ≤ 9 4 ( 2n 1
)

n k=0 9M 2 2k ≤ 9 4 (2n + 1) 9M 2 2n+2 9M 2 2 -1 := g(n) 2 ,
where g : [0, +∞[ is the strictly increasing function defined for t ∈ [0, +∞[ as

g(t) = 3 2 1 9M 2 2 -1 √ 2t + 1 9M 2 t+1 . ( 4
E(h, Λn( ) ) ∞ ≤ 40 p-1 p-2 k=0 (k) ∞ ||h (p-1) || ∞ , (4.75)
where

(k) ∞ = ∞ ( Λ(k) n( ) ). Since λ (k) j+1 -λ (k) j
= 2, we can use Theorem 4.12, and thanks to (4.62), we get easily: As a conclusion, thanks to (4.61), (4.73) and (4.78), we have obtained for small enough:

E(h, Λn( ) ) ∞ ≤ C A 1 n( ) p-1 . ( 4
||B[q f -q f ]|| H δ ≤ C A 1 log( 1 ) p-1 .
(4.79) By Proposition 4.6, B : H δ → H δ is an isomorphism and ||B -1 || ≤ C A , then we get a stability estimate between the two potentials q f and q f for the topology of H δ : First, let us prove a stability estimate for the warping functions f and f in the variable x ∈ (0, +∞). In order to simplify the notation, we set F = f d-2 , F = f d-2 . We we write: We integrate again this equality on (x, ∞) and we get: and the proof of Theorem 1.5 is now complete.

||q f -q f || H δ ≤ C A 1 log( 1 )
( F F -F F ) (y) = F F (q f -q f )(y) ≤ C A e -(d-2)y |(q f -q f )(y)|,
|1 - F F (x)| ≤ C A e -δ

The singular case.

We only sketch the proof since it is very similar to the previous one in the regular case. The real difference lies in the fact that we have no explicit formula for the angular eigenvalues κ k , and thus we have to use the Weyl asymptotics (1.13). First, let us consider the sub-sequence ν k = κ k d-1 . Thanks to Lemma 4.7, we get:

+∞ 0 e -2ν k x B[q f -q f ](x) dx ≤ C A , ∀k ≥ 0. (4.84) Thus, making the change of variables t = e -x , we obtain immediately:

1 0 t 2ν k -1 B[q f -q f ](-log t) dt ≤ C A , ∀k ≥ 0. (4.85)
Now, let us choose α ∈] 1 2 , 3 2 [ and let N be an integer large enough which we shall specify below. We deduce easily from (4.85) that:

1 0 t λ k h(t) dt ≤ C A , ∀k ≥ 0, (4.86)
where we have set h(t) = t -α B[q f -q f ](-log t) and λ k = 2ν N k + α -1. Note that ||h|| L 2 (0,1) = ||B[q f -q f ]|| H δ with δ = (2α -1) ∈ (0, 1[.

First, let us verify that we are in the separate case in the framework of the Müntz-Jackson approximation, i.e λ k+1 -λ k ≥ 2. It follows from the Weyl's law (1.13) that there exists C ≥ 1 such that: Now, following exactly the same approach as in the regular case, and taking θ = 1 B ∈ (0, 1[, we get Theorem 1.6. The details are left to the reader.

4 4 . 1 A 25 1

 4125 Proof of Theorems 1.5 and 1.6 13 Volterra type integral operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.2 A Müntz-Jackson's theorem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.3 A Hausdorff moment problem with non-integral powers. . . . . . . . . . . . . . . . . . . . 20 4.4 The regular case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.4.1 Stability estimates for the potentials. . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.4.2 Stability estimates for the warping functions. . . . . . . . . . . . . . . . . . . . . . 24 4.5 The singular case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The model and statement of the results

0 S 0 S

 00 y)S ∞ (x, y * ), and using that for k = 0, 1, lim x→∞ S (k) ∞ (x, y) = 0 and (2.10), we have:(y -y * ) ∞ ∞ (x, y)S ∞ (x, y * )dx = ∆(y)d(y * ) -∆(y * )d(y).Hence letting y * → y, we get ∞ ∞ (x, y) 2 dx = d(y) ∆(y) -∆(y) ḋ(y) = ∆ 2 (y) Ṁ (y),

1 ,

 1 we deduce from (1.13) there exists C > 0 such that |λ k -k| ≤ C. For a fixed N ∈ N large enough, we set ν k = λ kN N , and we have |ν k -k| ≤ C N < 1 4 . Using (3.7), a straightforward calculation gives: b 0 g(y)e -ν k y dy = O(e -bν k ), k → +∞, (3.8) where we have set b = c d-1 N a, and g(y) = f ( y c d-1 N ). Now, let us define for z ∈ C, F (z) = e bz b 0 g(y)e -zy dy.(3.9)Clearly, F (z) is an entire function which obeys|F (z)| ≤ ||g|| 1 e bRe+(z) ,(3.10) where Re + (z) is the positive part of Re z. Moreover, from (3.8), we get F (ν k ) = O(1). It follows from a theorem of Duffin and Schaeffer ([5], Theorem 10.5.1) that F (x) is bounded for x > 0, or equivalently b 0 g(y)e -xy dy = O(e -bx ) , x → +∞. (3.11) Then, using ([29], Lemma A.2.1), we have g = 0 almost everywhere in (0, b) which concludes the proof of this Proposition.

.47) 4 . 3 A

 43 Hausdorff moment problem with non-integral powers.

1 0f

 1 .53) Now, let us consider the finite sequenceΛ := Λ n : 0 ≤ λ 0 < λ 1 < ... < λ n . (4.54) Assume that the (n + 1) first moments of a function f ∈ L 2 ([0, 1]) are equal to zero up to noise, i.e there exists > 0 such that | (t) t λ k dt| ≤ , ∀k = 0, ..., n. (4.55)

2

 2 Stability estimates for the warping functions.

-δ 2 y e δ 2 y 1 ,

 21 since for instance, f (y) = c(e -y )e -y 2 and c ∈ C(A). Integrating this last inequality on the interval (x, ∞), and setting G = F F -F F , we get:|G(x)| ≤ C A e -(d-2)x +∞ x |(q f -q f )(y)| dy ≤ C A e -(d-2)x +∞ x e |(q f -q f )(y)| dy ≤ C A e -(d-2+ δ 2 )x ||q f -q f || H δ ≤ C A e -(d-2+ δwhere we have used (4.80) and the Cauchy-Schwarz inequality. Writing F F = G F 2 , and since c ∈ C(A), we deduce from (4.81) that F F (y) ≤ C A e (d-2)y e -(d-2+ δ

  deduce easily that, for all x ≥ 0,| f (x) -f (x)| ≤ C A e -,in the variable r = e -x ∈ (0, 1], we get:|c(r) -c(r)| ≤ C A r

  |ν k -c d-1 k| ≤ C , ∀k ≥ 0. (4.87) Now, let us fix an integer N large enough such that B := c d-1 N > 3C. One deduces thatλ k+1 -λ k = 2[ν (k+1)N -ν kN ] ≥ 2[(B(k + 1) -C) -(Bk + C)] ≥ 2B -4C ≥ 2C ≥ 2.Second, let us estimate the index of approximation ∞ (Λ). We recall that have to estimate the Müntz coefficients C mj associated with our sequence λ k as in the previous section in the regular case. We recall thatC mj = 2λ m + 1 m-1 r=0 (λ j + λ r + 1) m r=0,r =j (λ j -λ r ) . (4.90)Clearly, for j > r (for instance), one hasλ j -λ r = 2(ν kN -ν jN ) ≥ 2[B(j -r) -2C]. Thus, choosing b ∈ (0, B -2C],we get immediately λ j -λ r ≥ b(j -r). We deduce that:| m r=0,r =j (λ j -λ r )| ≥ (2b) m m r=0,r =j |j -r| (4.91)In the same way, one has:|λ j + λ r | = |ν kj + ν kr + 2α -1| ≥ |B(j + r) + 2C + 2α -1| ≥ M |j + r + 1|,where M = max{B, 2C + 2α -1}. It follows there exists D > 0 such that |C mj | ≤ D M 2b m |C 0 mj |. (4.92)

  ||Ch|| H δ ≤ C A,δ ||h|| H δ . Then, B : H δ → H δ is bounded and we have ||B|| ≤ C A,δ .Now, let us prove that C : H δ → H δ is a Hilbert-Schmidt operator. To this end, we calculate: It follows that C is a Hilbert-Schmidt operator, and a fortiori C is a compact operator. So, if B is not an isomomorphism, then -1 must be an eigenvalue of C. But this is impossible for a Volterra operator with continuous kernel.

	+∞	+∞
	x |Ch(x)| ≤ C A,δ e -αx ||h|| H δ , 0 K 1 (x, t) h(t) dt For h ∈ H δ , using (4.11), we get immediately HS := 0 0 |K 1 (x, t)1 {t≤x} | 2 e δ(x+t) dx dt = +∞ 0 x 0 |K 1 (x, t)| 2 e δ(x+t) dt dx ≤ C A,δ +∞ 0 x 0 e -2αx e δ(x+t) dt dx ≤ C A,δ +∞ 0 xe -2αx e 2δx dt dx ≤ C A,δ , which clearly implies ||C|| 2 since α ∈]δ, p[.	(4.18) (4.19)

.17) Proof. By convention, in what follows, C A or C A,δ denote constants depending only on A, (or only on A and δ), which can differ from one line to the other. Let α ∈]δ, p[ be fixed. We split the operator B as B = Id + C where C is the Volterra operator given by Ch(x) = Now, let us estimate the norm of the inverse operator B : H δ → H δ . We denote by K n (x, t) the integral kernel of the operator C n , n ≥ 1. Clearly, these kernels satisfy the relation

  It remains to estimate ||h -π n( ) h|| 2 = E(h, Λ n( ) ) 2 . First, let us introduce for small enough,Λn( ) : λ k0 < λ k0+1 < ... < λ n( ) ,(4.74)where λ k0 > p -2. Obviously, E(h, Λ n( ) ) 2 ≤ E(h, Λn( ) ) 2 ≤ E(h, Λn( ) )) ∞ . So using Corollary 4.11 with r = p -1, we obtain:

		.71)
	At this stage, we have obtained
	||Π n h|| 2 2 ≤ B 2 2 g(n) 2 .	(4.72)
	Now, let us choose a suitable integer n to control properly the norm of the projection ||π n h|| 2 2 . We set n( ) := E (g -1 ( 1 √ )). Clearly, since g is an increasing function, one has g(n( )) ≤ 1 √ , and thanks to
	(4.72) we get immediately:
	||π n( ) h|| 2 2 ≤ B 2 .	(4.73)

* Research supported by the French National Research Projects AARG, No. ANR-12-BS01-012-01, and Iproblems, No. ANR-13-JS01-0006 † Research supported by NSERC grant RGPIN 105490-2011 ‡ Research supported by the French GDR Dynqua 1 2 The Steklov spectrum 7


A Déprez
email: deprez@univ-grenoble-alpes.fr

N Cotte

A Walpersdorf

Toward the generation of EPOS-GNSS products: Massive GNSS processing & evaluation of a statistical approach for the velocity field computation

1.Massive data set & processing method

The GNSS solution presented here is the result of a massive GNSS processing based on a double difference method (using GAMIT / GLOBK software, Herring et al., 2010). Our data set is composed of rinex files from more than 1500 permanent stations over 18 years (2000-2017 included). This represents more than 4 millions of rinex files being an amount of data of more than 2000G. Such a processing justifies the use of a high performance computing platform (CIMENT hosted at The massive processing implementation requires innovative strategies in order satisfy the efficiency, robustness & automation criteria for the computation of GNSS velocity field in a "big data" context. MIDAS software appears to be a good, fast and the most appropriate method (Blewitt et al., 2016).

We obtained a "MIDAS velocity" for 1362 stations. In order to validate this "unusual" appoach on our data set, we compared these velocities with the ones obtained using a more "traditional" approach (using PYACS, Nocquet et al. in prep.) I.e. : estimation of a linear trend on position time series after an automatic or semiautomatic cleanage.

Horizontal velocity field computed with MIDAS Vertical velocity field computed with MIDAS

We presented here the statistics associated with the comparison between the MIDAS evaluation and the linear regression of the time series after spending a lot of time in order to carrefully visualize the time series data set and to correct them using human eyes and knowledge.

The standard deviations of the differences were 0.182, 0.240 and 0.509 mm/ yr respectively for the east, north and vertical component.

There were still 74 remaining stations with a difference > 0.5 mm/yr for the horizontal velocity and 21 stations with a difference > 1.5 mm/yr for the vertical velocity. However they are mainly located in active zones and are associated with area for which a constant velocity is meaningless, regardless the chosen method. Some differences in the vertical component were also associated with a high level of seasonal noise du to snow.

-

Differences in terms of horizontal velocities between the statistical approac (MIDAS) and the linear regression (PYACS), with respect to

Eurasia fixe.

Differences in terms of vertical velocities between the statistical approac (MIDAS) and the linear regression (PYACS).

.

  the data download to the position time series generation -Network of the European permanent stations used in our solution 2.Approaches tested to evaluate the velocity field 3.Comparison of the results in terms of velocity

  References: Herring, T. A., King, R. W. & McClusky, S. C. Reference Manual (GPS Analysis at MIT) Gamit-Globk Blog http://www-gpsg.mit.edu/∼simon/ gtgk/GAMIT_Ref.pdf (2010) Blewitt, G., Kreemer, C., Hammond, W. C. & Gazeaux J., MIDAS robust trend estimator for accurate GPS station velocities without step detection, J. Geophys. Res. Solid Earth, 121, 2054-2068 (2016), doi:10.1002/2015JB012552 Nocquet et al., in preparation PYACS, a framework (i.e. a set of tools and a way to think) to analyze and model geodetic data References: RENAG: doi:10.15778/RESIF.RG RGP: http://rgp.ign.fr ORPHEON: http://reseau-orpheon.fr RING (EPOS-glass node): http://ring.gm.ingv.it NOA (EPOS-glass node): http://194.177.194.238


Selima Bennaceur
email: selima.bennaceur@iup.univ-evry.fr.

Naoufel Azouz

Azgal Abichou

AN EFFICIENT MODELING OF FLEXIBLE BLIMPS: EULERIAN APPROACH

Unmanned Aerial Vehicles (U.A.V.) have a need of a greater autonomy in their new missions. Autonomous U.A.V. flight control systems require a precise modeling of the dynamic behavior taking into account the effect of the flexibility and the interaction with the surrounding fluid. In this paper, we present an efficient modeling of the autonomous flexible blimps. These flying objects are assumed to undergo large rigid-body motion and small elastic deformations. The formalism used is based on the Newton-Euler approach. This one is frequently used for flying rigid objects. In this study we develop a method to generalize the existing Newton-Euler "rigid body" formalisms by including the effect of the flexibility without destroying the global methodology. The method is hybrid. It uses the Lagrange equations and the Eulerian variables. The flexibility appears in the global dynamical system by the way of few supplementary degrees of freedom. This method has the advantage of making easier the elaboration of algorithms of control, stabilization or generation of trajectories. The added mass phenomenon is also taken into account in the dynamical system. This phenomenon is important for big and light objects moving in a fluid such as airships. As validation we use the parameters of an AS-200 blimp belonging to the University of Evry.

INTRODUCTION

Autonomous aerial vehicles such as airships and drones have recently gained importance. The airships began these last decades to retrieve a new youth after a half century of hibernation. Capability of airships are expanding rapidly now, and the range of missions they designed to support is growing.

We can give the examples of climate research, surveillance or even infiltration in war scenarios, and more closer to us we always saw them as means of advertising. A common denominator in all these situations is the impossibility or unwillingness to have human presence at the scene. In that case a precise dynamic model should be elaborated to permit an easy control, stabilization or navigation task of this autonomous object. In order to fulfill this requirement, it is necessary to introduce the effect of the structural flexibility in the dynamic model. It is important to note that several kinds of airships, usually called blimps, are mainly constituted of a balloon filled with gas. The only solid parts are the careen and the tail fins (for more details see [START_REF] Khoury | Airship Technology[END_REF]). The integration of the structural flexibility in the dynamic analysis is then useful; however it is now in an embryonic state and is only just emerging. Several researches were done by using the assumption of rigid body behavior for airships [START_REF] Fossen | Guidance and control of ocean vehicles[END_REF], [START_REF] Bestaoui | Dynamic modeling of small autonomous blimps[END_REF]. The flexibility effects are sometimes modeled as a perturbation. However in other flying objects, such as light aircrafts, the introduction of the flexibility in the dynamic model becomes essential [START_REF] Calise | Adaptive compensation for flexible dynamics[END_REF][START_REF] Bianchin | State space reduced order models for static aeroelasticity and flight mechanics of flexible aircrafts[END_REF] . Thus we try to contribute to the study of the deformation of the airships by introducing the effect of flexibility as non controlled supplementary degrees of freedom. The deformation of the blimp is not considered as a perturbation but rather acting on the motion of the airship. The influence of structural flexibility on the dynamics of mechanical systems has become increasingly important in classical robotics [START_REF] Simo | The role of non-linear theories in transient dynamic analysis of flexible aircrafts[END_REF][START_REF] Pascal | Vibrations analysis of flexible multibody systems[END_REF][START_REF] Schiehlenw | Multibody dynamics in computational mechanics and engineering applications[END_REF], and recently in flying robots (i.e. Airships, drones…).

Several approaches, to study the problems of flexible bodies, have been proposed in the literature. These approaches can be classified into two groups. The first one uses the total Lagrangian method [START_REF] Shabana | Dynamics of multibody systems[END_REF]. This consists in defining the motion relatively to a fixed reference frame, but this often leads to complex relations when describing stresses and strains in the flexible body subject to large displacements and small deformations. An Updated Lagrangian Method (U.L.M.) was proposed by Bathe & al. [START_REF] Bathe | Finite Elements for large deformation dynamic analysis[END_REF] and developed for deformable bodies that undergo large translational and rotational displacements. The resolution of the dynamic problem is incremental. The configuration and the motion of the body are identified using a moving reference configuration representing the position of the deformable body in the preceding step. Bennaceur & al[11] propose as reference configuration a rigid body configuration which follows the motion of the body without coinciding with it. This approach is convenient for a flying body with small deformations. The motion is given by coupled sets of rigid and elastic variables. The nonlinear equations are formulated in terms of a set of time invariant matrices expressed in a reference configuration (i.e. mass and stiffness matrices). Time-variant quantities appear in the nonlinear terms that represent the dynamic coupling between the rigid body modes and the elastic deformation. A suitable technique to actualize these terms using matrix partitioning and canonical decomposition is proposed, and the dynamic system is reduced and solved via a modal synthesis.

The Newton-Euler description [START_REF] Boyer | Generalization of Newton-Euler model for flexible manipulators[END_REF], which is an interesting method in regard of the time computation, was extensively studied by another important group of searchers in the case of rigid airships. This choice is mainly motivated by the facility to build control or stabilization algorithms based on this model. However the use of the Newton-Euler approach in the dynamic analysis of flexible complex structures is rare [START_REF] Boyer | Generalization of Newton-Euler model for flexible manipulators[END_REF]. We propose through this paper to extend the classical rigid bodies' model to the deformable bodies, without destroying the general formalism obtained. However this requires a total reformulation of the spatial flexible structure modelization by the use of a hybrid method based on both Lagrangian and Eulerian description. The model obtained should be acceptable for control community and closer to the one used for the rigid bodies.

On the other hand airships are also governed by the aerodynamic forces that have to be modeled. The basis to analyze the motion of a rigid body in a perfect fluid has been established in the 19th century and has been described by Lamb [START_REF] Lamb | On the motion of solids through a liquid[END_REF] . In his work, Lamb considered the case of simple displacement in a big infinite mass of fluid where the movement of this last is entirely due to the motion of a solid, and it is irrotational and acyclic. He proved that the kinetic energy of the fluid can be expressed as a quadratic shape of the six velocities of translation and rotation of the vehicle. The derivations given by Lamb will be used in the description of the airship, in a stationary uniform atmosphere. The terms depending on the acceleration (i.e. the added masses) come from the fact that the fluid considered perfect is accelerated. When an ellipsoid body moves in an incompressible and infinite inviscid fluid in order that the external flow is everywhere irrotational and continuous, the kinetic energy of the fluid produces an effect equivalent to an important increase of the mass and of the moments of inertia of the body [START_REF] Thomasson | Equation of motion of a vehicle in a moving fluid[END_REF], [START_REF] Ortega | On the motion of a rigid body immersed in a bidimensional incompressible perfect fluid[END_REF], [START_REF] Fossen | Guidance and control of ocean vehicles[END_REF]. In this paper, we introduce the added mass phenomenon in our model to have a complete dynamic model of a flying flexible object.

EULERIAN DESCRIPTION OF A RIGID BODY:

Kinematics:

We present an hybrid approach based first on an energetic method and using the Lagrange equations. The analysis of the motion is made with respect to an earth-fixed reference frame 0 0, 0 0 ( , , ) = R O X Y Z . A parallel is then done with the Eulerian approach, and the dynamical system is expressed with the Eulerian variables. The origin of the reference frame can be chosen in an arbitrary point such as, for example, the initial position on the ground of the airship. Its axes are selected as it follows: [ , ] ν ν ν

t t t = such as: 1 2 [ , , ], [ , , ] ν ν u v w p q r = =
are the linear and rotational speeds expressed in m R , ν ɺ are the relative acceleration expressed in m R .

Rotation matrices:

Commonly in aeronautics, a parameterization in yaw, pitch and roll is used to describe the orientation of the airship in the inertial reference frame. The configuration of the blimp is described by means of three rotations defined by three angles of orientation.

Let us note H ψ , H θ , H φ the orthogonal matrices of rotation.

The whole transformation between 0 R and m R is the combination of these preceding matrices and is given by: 1 cos cos sin cos sin cos sin sin sin sin cos cos ( ) cos sin cos cos sin sin sin cos sin cos sin sin sin cos sin cos cos

J ψ θ ψ φ φ ψ θ φ ψ θ ψ φ θ ψ ψ φ θ ψ φ ψ φ φ θ ψ θ θ φ θ φ - + +     = + - +     -   2 η such as: 1 1 1 1 3
( ). ( ) ( ). ( )

2 2 2 2 η η η η T T J J J J I = =
I 3 : the identity matrix.

Transformations of velocities:

Using the rotation matrix 1 ( ) J 2 η , the expression of the linear velocity in 0 R is given by:

1 ( ). J = 1 2 1 η η ν ɺ (1)
On the other hand, the angular velocity of the blimp 2 ν is the combination of the angular velocity around the three axes of yaw, pitch and roll. It can be written related to ɺ 2 η as:

1 0 sin 0 cos sin cos . 0 sin cos cos

p q r θ φ φ φ θ θ φ φ θ ψ   -           = =             -       2 ν ɺ ɺ ɺ (2) or:
( )

2     -       = = -         -       2 η
It is noticed that the Euler angles parameterization have a singularity in

2 k π θ π = +
.This parameterization is acceptable because it is impossible for a blimp to reach this orientation of 90 degrees pitching angle.

The global kinematics equation is then:

1 2 ( ) 0 0 ( ) J J    =       2 1 2 2 η ν η η ν ɺ (4)

Notations:

Subsequently we can note 1 ( ) J 2 η by 1 J and 2 ( ) J 2 η by 2 J .

Dynamics:

In order to be clearer when developing the flexible body formalism, we present here the analysis of a rigid airship by a hybrid methodology based on both Lagrangian and Euler description. The analysis of the flexible blimps is presented as an extension of the one developed for rigid flying objects.

We assumed that the airship is a rigid body with a mass m, and G as centre of gravity.

The position of an arbitrary point P expressed in the total reference frame, of the airship is given by the following vectorial relation: The position of a given point P in the blimp is:

OP = OG + GP

1 J = + 1 0 q η u (5)
Let us note u 0 the local position vector of the point P in R m .

By deriving ( 5) we find the following velocity's expression:

1 1 J J + + 1 0 0 0 q = η u u ɺ ɺ ɺ ɺ (6)
After computation we obtain:

1 1 1 . . ( ) . ( ) = ∧ = - ∧ ɺ 0 2 0 0 2 u ν u u ν J J J 1 1 0 1 0 2 J u J u J - = - = - 2 2 ν η ɺ ɶ ɶ (7)
0 u ɶ is the skew-matrix associated to the vector 0 u . ∧ is the vectorial product of two vectors. Thus the velocity vector could be expressed in a compact form as:

1 3 1 2 -     =       ɺ ɺ ɶ ɺ 1 0 2 η q u η I J J (8)
We will use this last form of the velocity vector thereafter in the development of the kinetic energy.

The following definition of the kinetic energy is used:

1 2 c V E d V ρ = ∫ T q q ɺ ɺ (9)
ρ and V are respectively the mass density and the volume of the body. This leads to the following relation:

1 2 c E M = T η η ɺ ɺ ( 10
)
Where M is recognized as the symmetrical mass matrix and defined as follows:

1 3 1 2 1 1 1 1 0 2 1 2 1 2 ( ) T V I J J M dV J u J J J J J ρ - - - -   -   =   - - -   ∫ 0 0 0 u u u ɶ ɶ ɶ ɶ (11)
In a more symbolic form the mass matrix is:

TT TR RT RR m m M m m   =     With: TT V m I dV ρ = ∫ 1 1 2 ( ) TR RT V m J J dV m ρ - = - = ∫ 0 u ɶ , 1 1 1 2 1 2 ( ) ( ) T RR V m J J J J dV ρ - - = - - ∫ 0 0 u u ɶ ɶ
We note that the term TT m , associated with the translation motion, is invariant. However, the other matrices depend on the generalized co-ordinates and consequently they are implicit functions of time. As long as the centre of gravity of the blimp coincides with the centre of the reference frame we can consider that the two extra-diagonal sub-matrices of m are null. Thus:

3*3 3*3 0 0 TT RR m M m   =     (12)

Lagrange's equations:

Let us use the Lagrange's equations to define the motion of the blimp. This leads to the following relation:

T ( ) ∂ ∂ - = ∂ ∂ ɺ τ η η C C E E d dt (13)
τ represents the external forces and moments.

First of all let us evaluate the first left hand term:

( ) [ ( )] C TT RR RR E d m m m dt ∂ = + ∂ T T T 1 2 2 η η η η ɺɺ ɺɺ ɺ ɺ ɺ (14)
The inertia matrix RR m is expressed in 0 R . If we note by RR I the constant symmetric inertia matrix in m R , we obtain:

1 1 2 2 ( ) T RR RR m J I J - - = (15)
By using the following relations [START_REF] Goldestein | Classical Mechanics[END_REF] :

1 2 0 J - = 2 η ɺ ɺ (16) 1 1 2 2 J J - - = - 2 2 η η ɺ ɺ (17)
and (3) :

1 2 J - = 2 2 ν η ɺ
and while differentiating with respect to time the equation [START_REF] Ortega | On the motion of a rigid body immersed in a bidimensional incompressible perfect fluid[END_REF] , we obtain:

1 2 RR RR m I J - = T T 2 2 η ν ɺ ɺ ɺ (18)
When replacing the last relation in ( 14) it gives:

1 2 ( ) [ ( )] C TT RR RR E d m m I J dt - ∂ = + ∂ T T T 1 2 2 η η ν η ɺ ɺɺ ɺɺ ɺ (19)
The second left term of Lagrange's equation is:

( ) 15), this gives:

3 1 3 2 1 0 ( ) 2 0 T C RR T RR E m I J -   ∂ ∂ =   ∂ ∂   = T 2 2 2 T 2 η η η η -ν ɺ ɺ ɺ (
1 1 1 2 2 2 1 1 1 2 2 2 ( ) 2 ( ) T R R T R R J J I J J J J I - - - - - - = - 2 2 2 η τ ν ɺɺ ɺ (22)
The following relations are proved (see for example [START_REF] Shabana | Dynamics of multibody systems[END_REF]):

1 1 2 2 3 ( ) 4 T J J I - - = (23) 1 1 2 2 ( ) 2 T R R R R J J I I - - = ∧ 2 2 2 ν ν ν ɺ (24)
Thus, relation (22) becomes:

1 2 1 2 4 4 R R R R I J J I - - = - ∧ 2 2 2 2 η τ ν ν ɺɺ (25)
Using relations (3) and ( 16) one can easily find:

1 2 J - = 2 2 ν η ɺ ɺɺ (26)
And we obtain the dynamic relation of rotation expressed with Eulerian variables:

1 2 1 4 - = -∧ ɺ 2 2 2 2 ν τ ν ν RR RR I J I (27)
After recalling the classical differential relation:

1 1 J - = + ∧ 1 1 2 1 η ν ν ν ɺɺ ɺ (28)
The dynamical system of the blimp expressed with the Euler variables becomes:

1 1 1 2 ( ) 0 1 0 ( ) 4 - -   - ∧       =      -∧         ɺ ɺ ɺ 1 2 1 1 2 2 2 2 ν τ ν ν ν ν τ ν ν E TT TT RR RR M J m m I J I (29)
Or in compact form:

E G M Q = + ν τ ɺ ( 30
) E M
is the constant mass matrix in the Eulerian description,

-1 1 -1 2 J 1 J 4     =       1 2 τ τ τ
is the vector of external forces, and

finally TT G RR m ( ) Q (I ) - ∧   =   -∧   2 1 2 2
ν ν ν ν is a vector containing the gyroscopic and Coriolis forces.

EULERIAN DESCRIPTION OF A FLEXIBLE BODY

In the flexible case, the motion analysis will be similar to that used for rigid body, but with the adjunction of supplementary degrees of freedom (d.o.f). Indeed the body's flexibility will be represented by a limited number (nd) of lower frequency deformation modes, resulting from the Rayleigh-Ritz decomposition. The choice of the number of modes to be retained depends on the structure and the solicitations applied (see for example [START_REF] Gibert | Vibrations des structures. Interactions avec les fluides[END_REF], [START_REF] Clough | Dynamics of Structures McGraw[END_REF]).

Dynamics:

Figure 3. Definitions of the position vectors.

The position of a point in the blimp is given by: (see figure 3)

' OP' = OG+GP+ PP (31)
Knowing that P' is a point of the deformed configuration.As in the rigid case, we will start with a Lagrangian description. Consequently all the dynamic parameters are expressed in 0 R . Consequently:

1 1 J J = + = + 1 0 d 1 U q η (u +u ) η U (32)
By deriving the relation (32) with respect to time, and when using similar developments to those presented in (5), we obtain:

1 1 2 1 J J J - = - + 1 2 d q η U η u ɶ ɺ ɺ ɺ ɺ (33)
The

Y

: the associated modal amplitude. This writing can be condensed in the following way: The modes of deformation S being constant compared to time, we can thus deduce that the speed of a given point from the blimp will be:

S = d d u Y (35)
1 1 2 1 J U J J S - = - + 1 2 d q η η Y ɺ ɶ ɺ ɺ ɺ (36)
The preceding form of the speed will be used for the development of the kinetic energy of a deformable body.

Expression of the kinetic energy of an arbitrary point P':

Similarly to the developments presented in § 2, we will have:

1 2 η η = ɺ ɺ T c E M (37)
M is now defined by: ( )

1 1 1 2 1 2 1 1 ρ - -     = - -         ∫ ɶ ɶ V I M J U J I J U J J S dV J S ( 38
)
The mass matrix can be written in a symbolic form as follows:

    =         TT TR TD T TR RR RD T T TD RD DD m m m M m m m m m m with: TT V m IdV ρ = ∫ ; 1 1 2 TR V m J UJ dV ρ - = -∫ ɶ ; T DD V m S SdV ρ = ∫ ; 1 TD V m J SdV ρ = ∫ , 1 1 2 2 ( ) T T RR V m J U UJ dV ρ - - = ∫ ɶ ɶ 1 2 1 ( ) T T RD V m J U J dV ρ - = ∫ ɶ .
It is clear that the two sub-matrices respectively associated with the translation and the elastic co-ordinates are constant.

The other matrices depend on the generalized co-ordinates of the system, and consequently they are implicit functions of time.

By coinciding the centre of gravity of the machine and the centre of m R , and by assuming that the position of this last does not vary significantly in m R in presence of small deformations, we can thus consider that the sub-matrix m TR is null. We will also assume here that the terms of coupling translation-deformation m TD are void. The kinetic energy is then:

1 ( 2 2
)

C TT RR RD DD E m m m m = + + + T T 1 1 2 2 T T 2 d d d η η η η η Y Y Y ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ (39)
By using the Lagrange's equations, we obtain the following relation:

G M K Q = - + η τ η ɺɺ ( 40
)
Where τ is the vector of the external forces and moments, and K is the matrix of stiffness given by: 0 0 0 0 0 0 0 0

DD K K     =      
The relation (40) can be given in a matrix form as follows:

1 2 0 0 0 0 0 0 0 0 0 0 0 0 TT RR RD T RD DD G G DD Gd m m m m m Q Q K Q η η                                       + = +                         1 2 d 1 1 2 2 d η η Y η τ η τ Y 0 ɺɺ ɺɺ ɺɺ (41) Then 0 0 G DD M Q K τ     = - +       d η Y ɺɺ (42)
Let us now analyze the vector of gyroscopic and Coriolis forces Q G .

From the Lagrange's equation, we have:

1 ( ) 2 T G Q M M   ∂ = - +   ∂   T η η η η ɺ ɺ ɺ ɺ (43)
Or in more explicit form:

1 2 1 ( ) 2 1 ( ) 2 1 ( ) 2 T G T G G RR RD Gd T RD DD M Q Q Q m m M Q m m M η η     ∂     ∂               ∂ = = - - +       ∂             ∂   - - +     ∂     T 1 T 2 d 2 T 2 d d η η η η Y η η η η Y η η Y ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ (44)
The development of this relation leads to these following results:

• The first component could be written as:

1 1 ( ) 0 2 T G Q M η   ∂ = =   ∂   T 1 η η η ɺ ɺ (45)
• For the second component

2 G Q η we have: 2 1 ( ) 2 T G RR RD Q m m M η   ∂ = - - +   ∂   T 2 d 2 η Y η η η ɺ ɺ ɺ ɺ ɺ ɺ (46)
The first term is:

.

1 1 2 2 ( ) T RR RR m J I J - - - = - 2 2 η η ɺ ɺ ɺ (47)
With RR I the matrix of rotation inertia expressed in m R . Basing on the fact that this analysis is intended for the airships subjected to small deformations, therefore we will keep the same assumptions used for the rigid case. Consequently the inertial terms expressed in the pointer are invariants.

Notations

We note RR I defined by According to this and to the equations ( 16) (3) and (47) we have:

.

1 1 1 2 2 2 ( ) ( ) T T RR RR J I J J I - - - - = - 2 2 η ν ɺ ɺ (48)
and we can say for the second term that:

1 2 ( ) T RD RD m J I - - = - d d Y Y ɺ ɺ ɺ ɺ (49)
For the third term:

1 1 2 2 1 1 ( ) ( ) 2 2 1 ( ( ) ) 2 RR T RR M m J I J - - ∂ ∂ = ∂ ∂ ∂ = ∂ T T 2 2 2 2 T 2 2 2 η η η η η η η η η ɺ ɺ ɺ ɺ ɺ ɺ (50)
By using the relation ( 16), and after differentiating, we will have:

1 2 1 ( ) ( ) 2 T T RR M J I -   ∂ = -   ∂   T 2 2 η η ν η ɺ ɺ ɺ (51) Thus: 2 1 1 2 2 2( ) ( )
T T G RR RD Q J I J I η - - = - - 2 d ν Y ɺ ɺ ɺ (52)
• For the third component G d Q we have:

0 1 ( ) 2 T G d RD T Q m M = -   ∂ +   ∂   2 T d η η η Y ɺ ɺ ɺ ɺ (53)
According to the preceding notations we can obtain the following relation:

1 2 0 0 T G d RD Q I J - = - = 2 η ɺ ɺ (54)
When collecting all these developments, the equation ( 41) could be expressed in this form:

1 1 2 2 0 0 0 0 0 2( ) ( ) 0 TT RR RD T dd RD DD T T RR RD I m m K m I J I J I - -             = +             -             - -         1 1 2 2 d d 2 d η τ η τ Y Y ν Y ɺɺ ɺɺ ɺɺ ɺ ɺ ɺ (55)
We will now transform this system equation ( 55) by the use of the Eulerian variables. The methodology is similar to that developed in the rigid case, and by using the relations (24), (25). This leads to the following system:

1 1 1 2 ( ) 0 0 1 0 ( ) ( ) 4 0 0 TT TT RR RD RR RD T dd RD DD J I I I I J I I K I I - - τ ν ν ν ν τ ν ν ν Y Y Y ɺ ɺ ɺ ɺɺ (56) Or in compact form E G M Q ν = + τ ɺ (57)
We thus obtain a Newton-Euler system equation similar to that defined in equation ( 30) but extended to the flexible blimps.

AERODYNAMIC CONTRIBUTION:

We show here the influence of the surrounding air on the general behavior of the flexible blimps. Let us consider a simple dynamical model for the action of the air on a body. To present this model, we assume that the flow is quasi-steady, i.e. the distribution of the velocities of particles of the medium coincides with the distribution corresponding to the steady motion of the body. Thus the medium responds only to the current motion of the body and forgets its initial conditions. Therefore, within the framework of this hypothesis, the resultant force and torque acting on the body can be represented in the form of a function of the instantaneous distribution of velocities in this body. Thus we reach at the statement of the problem of the motion of a body in a dragging medium as a problem of classical dynamics.

Flow representation:

To take into account the interaction of the airship with the surrounding fluid medium, a model of the flow is needed.

Here, we rely on the potential flow theory corresponding to the following hypothesis:

• the air can be considered as a perfect fluid with uniform density ρ air , i.e. an incompressible gas with vanishing viscosity, • the flow is irrotational Only the flow outside the airship contributes significantly to the aerodynamics. Denoting by v the velocity field in the fluid domain "air", the incompressibility and irrotational assumptions leads to:

. 0 ; 0

∇ = ∇ ∧ = v v (58)
∇ is the gradient symbol, and the flow field may be described in terms of a potential F such as:

= ∇ Φ v (59)
From the incompressibility constraint, it is easy to show that the potential obeys to the homogeneous Laplace equation: (see [START_REF] Katz | Low-speed Aerodynamics, from Wings Theory to Panel Methods[END_REF] for details on the numerical treatment).

Fluid forces:

For this assumption, the pressure at any point in the fluid domain (including ( )

t C ∂
) is given by Bernoulli theorem:

air air

1 1 P . P . . 2 t 2 ∞ ∞ ∞   ∂Φ +ρ + = + ρ   ∂   vv v v (62)
The subscript ∞ denotes the undisturbed conditions far from the airship. This pressure distribution over the airship surface can be integrated to compute the resulting forces and torques. At the end, and with the linear property of the Laplace equation, the generalized fluid forces vector can be rewritten as:

ad f M B = - - f F ν ν ɺ (63)
Where M ad is the matrix of the added masses (virtual masses), and Bf is a damping due to the flexibility of the hull. Taking into account the previous developments, the dynamic equation (57) becomes:

E G M Q ′ ′ = + ν τ ɺ (64) we note E E ad M M M ′ = + , G G f Q Q B ′ = -ν

SIMULATIONS RESULTS:

The method for computation is based on the following flow chart:

To illustrate this Eulerian approach we study the blimp belonging to the LSC-IBISC and having the following characteristics:

-The envelope: Length: 6.25 m. Diameter: 1.52 m. Volume: 7.48 m3.

-Mass of the airship: 5.8 Kg.

-Payload: 1.58 Kg. The blimp is thrust by two contrarotating propellers. In the first step we apply to the blimp two opposite forces in the propellers to generate a yaw motion (Fig 5,6). The flexible airship should rotate about 90° around the z-axis. A Proportional-Derivative (P.D) controller will impose this task, and we study the behavior of the airship during this maneuver.

The number of deformable modes kept is nd = 2. This number seems to give an acceptable approximation of the flexible behavior. We were guided for that by the modal masses of these two modes which represent roughly 70% of the total mass of the airship. We applied the force on the point P (2, 0, 0) of the blimp.

In this simulation we applied around the yaw angle a torque consigns :

With d K = 800, v K = 50, d ψ : The desired angular position and ψ : The angular position of the airship. One notices in figure 5 and 6 that the amplitude of the deviation of the yaw angle decreases significantly in a few seconds. it stabilizes itself while merging with the instruction. In figure 7 we superimpose total displacement along the Xaxis of the rigid airship and the flexible device. It is noticed that the flexible device continues to oscillate which proves the impact of flexibility on displacement. For the rigid behavior we eliminate the two deformable modes and keep only the rigid motion. The deformations (in figure 8) are about 0.20m what is more or less significant considering one has small deformations. In this simulation we visualized the prow of the airship, and we see well the oscillation of this later with non-negligible amplitude.

One classical motion of an airship is the helicoidally motion. This usually represents the optimal trim trajectories. This represents a complete aerial motion of the airship. The propellers were oriented adequately to assure a combined motion along the moving x-axis and z-axis and the tails were oriented to give a homogeneous yaw motion. We can see (in figure 9) the influence of the deformation on the overall motion of the airship. This proves that in our model, the deformation is not a sample perturbation around the main rigid motion, but it acts on this motion. Experimental data will be available in few months to validate the results of our scheme. Finally we show the influence of the added masses on the total displacement of the airship. It is noticed (in figure 10) that the immersed airship has a delay compared to the body alone and this is because of the addition of the virtual masses representing the influence of the displacement of the mass of air around the airship when this last accelerates.

CONCLUSION

In this paper we presented a description of the flexible airships by an Eulerian method. We introduced the flexibility of the airship as an extension of a rigid body algorithm based on Eulerian description and taking into account the coupling between the rigid body motion and the deformation. A hybrid methodology is used, based on both Lagrangian and Eulerian theory. The final goal of developing this model is to enable an easy implementation of control and stabilization algorithms for the flexible flying objects. The model built presents also an interesting ratio precision /computational time. Simulation results prove that the integration of the flexibility in the dynamic system of the airship is very important and could not be neglected. The added mass phenomenon is also taken into account. However we neglected in this study the added masses issued from the vibration of the airship in the air and their coupling with those issued from the overall rigid body motion. This point is under investigation and will be presented in future works.

  fixed at the centre of gravity of the airship G, its axes are selected as follows: m X : aligned with longitudinal axis, m Y : The transverse axis, m Z : The normal axis directed downwards.

Figure 1 .

 1 Figure 1. Definitions of the reference frames. The airship's motion is described by the following parameters: 0 0 0 [ , , ] T x y z = 1 η : Vector position of the origin expressed in 0 R , [ , , ] T φ θ ψ = 2 η : The Euler angles. The kinetic and dynamic Eulerian parameters are those used usually in aeronautics ν and ν ɺ . Where1 2

Figure 2 .

 2 Figure 2. Definition of the position vectors.

20) 3 0 2 . 6 2 J

 3262 T is the null vector of three dimensions. The Lagrange's equation becomes then: Eulerian description: In this part, we will define the motion of the blimp with the Eulerian variables starting from (21.a), (21.b) and presenting the correlation between the Lagrangian and Eulerian description. When multiplying by 1 -the relation (21.a), and according to the equation (

  vector of an arbitrary point of the flexible airship.

  of inertial coupling translation-deformation.As announced in the beginning, we will suppose that RR I , RD I , TD I are constant matrices with respect to time through the assumption of small deformations. matrices of translation and deformation.

  of the most important characteristic of this representation is that v only depends on the current boundary conditions, and not on the history of the flow: the model is quasi-steady. To solve the potential equation, we use the boundary integral representation of the Laplace equation, together with standard boundary element method. It consists in the determination of a piecewise constant distribution of singularities over ( ) t C ∂

  . The effect of the fluid on the structure is then represented mainly by the adjunction of the added masses matrix Mad to the mass matrix of the structure. For a quasi-ellipsoid airship the extra-diagonal terms of ad M can be neglected (for more details about the constitutive terms of ad M , the reader can see[START_REF] Fossen | Guidance and control of ocean vehicles[END_REF] ,[START_REF] Lamb | On the motion of solids through a liquid[END_REF]).

Figure 4 .

 4 Figure 4. The LSC-IBISC's airship.

Figure 5 .Figure 6 .

 56 Figure 5. Superposition of the deviation of the yaw angle and the instruction.

Figure 7 .

 7 Figure 7. Superposition of the motions of the flexible and rigid airship.

Figure 8 .

 8 Figure 8. The deformation's displacement.

Figure 9 .

 9 Figure 9. Rise of the airship.

Figure 10 .

 10 Figure 10. The impact of the fluid on the Blimp.