Fabien Caubet
email: fabien.caubet@math.univ-toulouse.fr

Marc Dambrine
email: marc.dambrine@univ-pau.fr

AND Helmut Harbrecht
email: helmut.harbrecht@unibas.ch

A New Method for the Data Completion Problem and Application to Obstacle Detection

Keywords: Cauchy problem, Data completion problem, Inverse obstacle problem, Laplace's equation, Kohn-Vogelius functional AMS subject classifications. 35R30, 35R25, 35N25, 49M15

The present article is devoted to the study of two well-known inverse problems, that is the data completion problem and the inverse obstacle problem. The general idea is to reconstruct some boundary conditions and/or to identify an obstacle or void of different conductivity which is contained in a domain, from measurements of voltage and currents on the outer boundary of the domain. We focus here on Laplace's equation.

Firstly, we use a penalized Kohn-Vogelius functional in order to numerically solve the data completion problem, which consists in recovering some boundary conditions from partial Cauchy data. The functional to be minimized is quadratic, hence we compute its minimum by solving the linearized equation. Secondly, we propose to build an iterative method for the inverse obstacle problem based on the combination of the previously mentioned data completion subproblem and the so-called trial method. The underlying boundary value problems are efficiently computed by means of boundary integral equations and several numerical simulations show the applicability and feasibility of our new approach. For the numerical simulations, we focus on star-shaped domains in the two dimensional case.

1. Introduction. We deal in this article with the inverse obstacle problem which is defined as follows. Let Ω be a bounded connected Lipschitz open set of R d , where d = 2 or d = 3, with a boundary ∂Ω. For some nontrivial Cauchy data (f m , g m ) ∈ H 1/2 (∂Ω) × H -1/2 (∂Ω) on the outer boundary, which corresponds to a pair of current/voltage measurements, we consider the following inverse problem: Then the problem under consideration is a special case (assuming a constant conductivity with a void of connected complement) of the general conductivity reconstruction problem also called Electrical Impedance Tomography (EIT). EIT is used in medical imaging to reconstruct the electric conductivity of a part of the body from measurements of currents and voltages at the surface (see, e.g., [START_REF] Lionheart | Electrical Impedance Tomography: Methods, History and Applications[END_REF]). The same technique is also used in geophysical explorations. An important special case consists in reconstructing the shape of an unknown inclusion or void assuming a constant conductivity. Notice that the continuum model is not a realistic model for EIT and more realistic electrode models exist (see [START_REF] Harrach | Monotonicity-based shape reconstruction in electrical impedance tomography[END_REF][START_REF] Harrach | Combining frequency-difference and ultrasound modulated electrical impedance tomography[END_REF][START_REF] Garde | Staboulis Convergence and regularization for monotonicity-based shape reconstruction in electrical impedance tomography[END_REF][START_REF] Lechleiter | The factorization method applied to the complete electrode model of impedance tomography[END_REF]). However we focus in this paper on this continuum model.

It is known that problem (1.1) admits at most one solution, as claimed by the following identifiability result (see, e.g., [START_REF] Bourgeois | A quasi-reversibility approach to solve the inverse obstacle problem[END_REF]Theorem 1.1] or [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]Theorem 5.1]).

Theorem 1.1. The domain ω and the function u which is assumed to be continuous up to the boundaries that satisfy (1.1) are uniquely defined by the nontrivial Cauchy data (f m , g m ).

It is also well-known that problem (1.1) is severely ill-posed: the problem may fail to have a solution and, even when a solution exists, the problem is highly unstable (see, e.g., [START_REF] Alessandrini | Optimal stability for inverse elliptic boundary value problems with unknown boundaries[END_REF][START_REF] Cristo | Examples of exponential instability for inverse inclusion and scattering problems[END_REF]).

This inverse problem has been intensively investigated by numerous methods. We can cite for example sampling methods [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF], methods based on conformal mappings [START_REF] Afraites | Conformal mapping and inverse conductivity problem with one measurement[END_REF][START_REF] Haddar | Conformal mappings and inverse boundary value problems[END_REF], on integral equations [START_REF] Kress | Nonlinear integral equations and the iterative solution for an inverse boundary value problem[END_REF][START_REF] Rundell | Recovering an obstacle using integral equations[END_REF], methods using the full Dirichlet-to-Neumann map at the outer boundary [START_REF] Brühl | Explicit characterization of inclusions in electrical impedance tomography[END_REF][START_REF] Brühl | Numerical implementation of two noniterative methods for locating inclusions by impedance tomography[END_REF], or level-set methods coupled with quasi-reversibility in the exterior approach [START_REF] Bourgeois | A quasi-reversibility approach to solve the inverse obstacle problem[END_REF]. We also refer for example to [START_REF] Chapko | A hybrid method for inverse boundary value problems in potential theory[END_REF][START_REF] Hettlich | The determination of a discontinuity in a conductivity from a single boundary measurement[END_REF] for numerical algorithms and to [START_REF] Alessandrini | Local uniqueness in the inverse problem with one measurement[END_REF][START_REF] Friedman | On the uniqueness in the inverse conductivity problem with one measurement[END_REF] for particular results concerning uniqueness.

In [START_REF] Roche | Numerical methods for shape identification problems[END_REF], the problem under consideration has been reformulated as a shape optimization problem for the Kohn-Vogelius functional. Then, seeking the unknown inclusion is equivalent to seeking the minimizer of an energy functional. Much attention has been spent on the analysis of this approach (see, e.g., [START_REF] Afraites | On second order shape optimization methods for electrical impedance tomography[END_REF][START_REF] Eppler | A regularized Newton method in electrical impedance tomography using shape Hessian information[END_REF]) and its comparison with a least-squares tracking type functionals. These kind of methods have some advantages such as being adaptable for several partial differential equations, such as the Stokes system (see, e.g., [START_REF] Caubet | On the detection of several obstacles in 2d Stokes flow: topological sensitivity and combination with shape derivatives[END_REF][START_REF] Caubet | A Kohn-Vogelius formulation to detect an obstacle immersed in a fluid[END_REF][START_REF] Guillaume | Topological sensitivity and shape optimization for the Stokes equations[END_REF]), and for obstacles characterized by different boundary conditions, such as Neumann or generalized boundary conditions (see, e.g., [START_REF] Badra | Detecting an obstacle immersed in a fluid by shape optimization methods[END_REF][START_REF] Caubet | Shape optimization methods for the inverse obstacle problem with generalized impedance boundary conditions[END_REF]).

In this article, we propose a new point of view to numerically solve this model inverse obstacle problem. We focus on the two-dimensional case and assume that the inclusion is star-shaped and perfectly conducting. Then, we build an iterative sequence of domains using the combination of a data completion subproblem and the so-called trial method which is used to control the evolution of the inclusion boundary. This means that, given an inclusion ω, we compute an harmonic function u ∈ H 1 (Ω\ω) which admits the Cauchy data (f m , g m ) at the outer boundary ∂Ω. With the help of the Cauchy data at the interior boundary ∂ω, we aim at updating the interior boundary such that the desired Dirichlet condition u = 0 holds at the new interior boundary. The data completion problem is solved by minimizing a Kohn-Vogelius functional thanks to the resolution of a regularized linearized equation, which is efficient since the functional is quadratic.

Organization of the article. The rest of the article is organized as follows.

In Section 2, we introduce the data completion problem which we intend to solve by means of the minimization of a Kohn-Vogelius functional, also introduced in this section. We additionally provide all the general notations and assumptions used throughout the article. In Section 3, we provide some properties concerning the considered boundary value problems and compute the gradient and the Hessian of the Kohn-Vogelius functional. Then, in Section 4, we introduce the Tikhonov regularization term needed to numerically minimize the functional. The minimization of this regularized functional is considered by means of a Newton scheme. Especially, numerical results are presented, which demonstrate that our approach yields an efficient solution of the data completion problem. In Section 5, we focus on the numerical resolution of the inverse obstacle problem. To this end, we present the so-called trial method, used to update an approximated inclusion's boundary from the information obtained from the previous data completion step. Numerical results validate that the proposed approach is feasible to reconstruct inclusions in the inverse obstacle problem at least on a model case a single star-shapes inclusion. Finally, in Section 6, we state concluding remarks. We denote by L 2 (Ω), L 2 (∂Ω), H 1 (Ω), H 1 (∂Ω), the usual Lebesgue and Sobolev spaces of scalar functions in Ω or on ∂Ω. The classical norm and semi-norm on H 1 (Ω) are respectively denoted by

• H 1 and |•| H 1 . Moreover •, • denotes the following product, for all u, v ∈ H 1 (Ω), u, v := Ω ∇u • ∇v so that |u| 2 H 1 (Ω) = u, u .
(f m , g m ) ∈ H 1/2 (∂Ω) × H -1/2 (∂Ω).
The data completion problem consists in recovering data on ∂ω, from the overdetermined data (f m , g m ) on ∂Ω, that is:

(2.1)

Find u ∈ H 1 (Ω\ω) such that    ∆u = 0 in Ω\ω, u = f m on ∂Ω, ∂ n u = g m on ∂Ω.
Of course, the previous problem is not standard. In general, it admits only a local solution (that is defined in the neighborhood of ∂Ω under regularity assumptions) by the Cauchy-Kowalevskaya theorem. Since we are interested in global solutions, this result is not satisfactory and we have to define the notion of compatible data. [START_REF] Fursikov | Optimal control of distributed systems[END_REF], where the following result is given. Lemma 2.2. We have the two following density results.

Definition 2.1. A pair (f m , g m ) ∈ H 1/2 (∂Ω) × H -1/2 (
1. For a fixed f m ∈ H 1/2 (∂Ω), the set of compatible data g is dense in H -1/2 (∂Ω). 2. For a fixed g m ∈ H -1/2 (∂Ω), the set of compatible f is dense in H 1/2 (∂Ω).

Therefore, any numerical scheme should incorporate a regularization step. Several approaches have been considered to solve the data completion problem. Among others, we mention the works of Kozlov et al. [START_REF] Kozlov | An iterative method for solving the Cauchy problem for elliptic equations[END_REF], Cimetière et al. [START_REF] Cimetiere | Solution of the Cauchy problem using iterated Tikhonov regularization[END_REF], Ben Belgacem et al. [START_REF] Belgacem | On Cauchy's problem: I. A variational Steklov-Poincaré theory[END_REF][START_REF] Ben Belgacem | The Lavrentiev regularization of the data completion problem[END_REF].

We shall follow the energy based strategy introduced by Andrieux et al. in [START_REF] Andrieux | Solving Cauchy problems by minimizing an energylike functional[END_REF] by the minimization of the Kohn-Vogelius like functional which admits the solution of problem (2.1) as minimizer, if such a solution exists. Notice that such a functional turns to be quadratic and convex. Surprisingly, to the best of our knowledge, only gradient based numerical schemes have been studied for the resolution of this problem. We will describe in this article the Newton method that is completely suited for quadratic objectives.

In order to deal with the ill-posedness previously mentioned, we consider a Tikhonov regularization of the functional which ensures the existence of a minimizer even for not compatible data thanks to the gained of coerciveness and, in case of compatible data, the convergence towards the exact solution (see, e.g., [16, Proposition 2.5 and Theorem 2.6]).

The Kohn-Vogelius functional.

As previously mentioned, the data completion problem (2.1) can be studied through the minimization of a Kohn-Vogelius cost functional (see [START_REF] Kohn | Determining conductivity by boundary measurements[END_REF]). To this end, we introduce the two maps u DN and u ND defined as follows (2.2)

u DN : H 1/2 (∂Ω) × H -1/2 (∂ω) -→ H 1 (Ω\ω), u ND : H -1/2 (∂Ω) × H 1/2 (∂ω) -→ H 1 (Ω\ω),
where u DN (f, ψ) ∈ H 1 (Ω\ω) solves the boundary values problem

(2.3)    ∆u = 0 in Ω\ω, u = f on ∂Ω, ∂ n u = ψ on ∂ω,
and where u ND (g, φ) ∈ H 1 (Ω\ω) solves the boundary values problem (2.4)

   ∆u = 0 in Ω\ω, ∂ n u = g on ∂Ω, u = φ on ∂ω.
Notice that the indices mean the type of boundary condition, where the first one indicates the outer boundary and the second one the inner boundary.

To tackle the inverse problem, one tries to solve the problem:

Find (φ, ψ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω) such that u DN (f m , ψ) = u ND (g m , φ).
Thanks to the linearity of the maps u ND and u DN , this is equivalent to:

Find (φ, ψ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω) such that u DN (0, ψ) -u ND (0, φ) = -u DN (f m , 0) -u ND (g m , 0) .
Notice that this is a standard linear inverse problem associated to the linear operator

A : H 1/2 (∂ω) × H -1/2 (∂ω) -→ H 1 (Ω\ω) (φ, ψ) -→ u DN (0, ψ) -u ND (0, φ).
Of course, the right hand side b := -u DN (f m , 0) -u ND (g m , 0) may belong or may not belong to the range R(A) of A, depending on the fact that data are compatible or not. We try to solve the previous linear equation in a least squares meaning. We focus on the following optimization problem

(φ * , ψ * ) ∈ argmin (φ,ψ)∈H 1/2 (∂ω)×H -1/2 (∂ω) K(φ, ψ),
where

K : H 1/2 (∂ω) × H -1/2 (∂ω) → R is the non-negative Kohn-Vogelius cost func- tional defined by K(φ, ψ) := 1 2 |A(φ, ψ) -b| 2 H 1 (Ω\ω) .
Notice that this criterion is quadratic as the square of a affine map and also writes (2.5)

K(φ, ψ) = 1 2 |u DN (f m , ψ) -u ND (g m , φ)| 2 H 1 (Ω\ω) = 1 2 Ω\ω |∇u DN (f m , ψ) -∇u ND (g m , φ)| 2 .
The connection between the original inverse problem and the minimization is given by next lemma.

Lemma 2.3. The inverse problem (2.1) has a solution if and only if the functional K has a minimizer (φ, ψ) and K(φ, ψ) = 0. Notice that the minimizer is unique up to an additive constant: if

K(φ 1 , ψ 1 ) = 0 = K(φ 2 , ψ 2 ) then there is a constant c such that φ 1 = φ 2 + c and ψ 1 = ψ 2 .
Proof. For a given and known inclusion ω, if the inverse problem (2.1) has a solution,

then u DN (f m , ψ) = u ND (g m , φ) and K(φ, ψ) = 0. Conversely, if we assume that K(φ, ψ) = 0 then ∇u DN (f m , ψ) = ∇u ND (g m , φ) and there is a constant c such that u DN (f m , ψ) = u ND (g m , φ) + c = u ND (g m , φ + c).
The usual theory of linear inverse problems should be applied. That is, we solve the normal equation A A(φ, ψ) = A b or its regularized version (A A + εB B)(φ, ψ) = A b, where B is a regularization operator and where ε > 0. However it involves the adjoint A of A and hence to manipulate the scalar product in the space H 1/2 (∂ω) × H -1/2 (∂ω).

The main difficulty of this approach is that the scalar product in the spaces H 1/2 (∂ω) and H -1/2 (∂ω) is not tractable with from a practical point of view when one uses the optimize then discretize approach. On the converse, if one first discretizes the equations, and then computes the discrete adjoint, the usual theory can be applied. In the following, we want to remain at the continuous level and reduce the minimization of the Kohn-Vogelius objective to such a linear inverse problem in a simple way thanks to the Newton point of view. Since the objective K is quadratic, its hessian is constant and its minimizer can be computed by a single step in Newton method.

Remark 2.4. Note that the two problems (2.3) and (2.4) are well-posed for any given boundary conditions (φ, ψ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω), without additional compatibility conditions between f m and ψ for the first problem and between g m and φ for the second. This is of particular interest for numerical implementations, as the considered setting allows to consider the classical Sobolev spaces and, therefore, the implementations can be done with classical finite element method softwares without any additional adjustments.

3. Some properties concerning the data completion problem.

3.1. Properties of the maps u DN ant u ND . Let us first emphasize some properties of the maps defined in (2.2) we should use. Proposition 3.1 (Properties of the boundary values problems). We have the following statements.

1. The maps u DN and u ND are linear: for all

(f 1 , ψ 1 ), (f 2 , ψ 2 ) ∈ H 1/2 (∂Ω) × H -1/2 (∂ω), there holds u DN (f 1 + f 2 , ψ 1 + ψ 2 ) = u DN (f 1 , ψ 1 ) + u DN (f 2 , ψ 2 ) = u DN (f 2 , ψ 1 ) + u DN (f 1 , ψ 2 )
and

u ND (g 1 + g 2 , φ 1 + φ 2 ) = u ND (g 1 , φ 1 ) + u ND (g 2 , φ 2 ) = u ND (g 2 , φ 1 ) + u ND (g 1 , φ 2 ).
2. There are positive constants C 1 and C 2 such that

(3.1) C 1 φ H 1/2 (∂ω) u ND (0, φ) H 1 (Ω\ω) C 2 φ H 1/2 (∂ω) , C 1 ψ H -1/2 (∂ω) u DN (0, ψ) H 1 (Ω\ω) C 2 ψ H -1/2 (∂ω) , for all (φ, ψ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω). 3. There holds u DN (f 1 , ψ 1 ), u DN (f 2 , ψ 2 ) = ∂Ω f 1 ∂ n u DN (f 2 , ψ 2 ) + ∂ω u DN (f 1 , ψ 1 )ψ 2 = ∂Ω f 2 ∂ n u DN (f 1 , ψ 1 ) + ∂ω u DN (f 2 , ψ 2 )ψ 1 , u DN (f 1 , ψ 1 ), u ND (g 1 , φ 1 ) = ∂Ω f 1 g 1 + ∂ω u DN (f 1 , ψ 1 )∂ n u ND (g 1 , φ 1 ) = ∂Ω u ND (g 1 , φ 1 )∂ n u DN (f 1 , ψ 1 ) + ∂ω φ 1 ψ 1 , u ND (g 1 , φ 1 ), u ND (g 2 , φ 2 ) = ∂Ω u ND (g 1 , φ 1 )g 2 + ∂ω φ 1 ∂ n u ND (g 2 , φ 2 ) = ∂Ω u ND (g 2 , φ 2 )g 1 + ∂ω φ 2 ∂ n u ND (g 1 , φ 1 ), for all (f 1 , ψ 1 ), (f 2 , ψ 2 ) ∈ H 1/2 (∂Ω) × H -1/2 (∂ω) and all (g 1 , φ 1 ), (g 2 , φ 2 ) ∈ H -1/2 (∂Ω) × H 1/2 (∂ω).
Proof. The first assertion is a trivial consequence of the superposition principle. Let us prove the point

C 1 ψ H -1/2 (∂ω) u DN (0, ψ) H 1 (Ω\ω)
. By definition of the dual norm,

ψ H -1/2 (∂ω) = sup φ H 1/2 (∂ω) =1 φ, ψ H 1/2 (∂ω)×H -1/2 (∂ω) = ∂ω ∂ n u DN (0, ψ)Φ = Ω\ω ∇u DN (0, ψ) • ∇Φ,
where Φ ∈ H 1 (Ω\ω) is the harmonic extension of φ so that it holds

Φ H 1 (Ω\ω) = φ H 1/2 (∂ω) = 1.
We conclude that

ψ H -1/2 (∂ω) u DN (0, ψ) H 1 (Ω\ω) Φ H 1 (Ω\ω) = u DN (0, ψ) H 1 (Ω\ω) .
The others inequalities of the second second assertion comes from usual a priori elliptic estimate. The third assertion is a consequence of the Green formulae.

Let us remark that the third assertion of Proposition 3.1 implies in particular the following formulae which is useful to transform an integral over the outer boundary into an integral over the inner boundary:

u DN (0, ψ), u DN (f, 0) = 0 = ∂Ω f ∂ n u DN (0, ψ) + ∂ω u DN (f, 0)ψ, ∂Ω u ND (0, φ)∂ n u DN (f, 0) = ∂ω u DN (f, 0)∂ n u ND (0, φ), u ND (g, 0), u ND (0, φ) = 0 = ∂Ω u ND (0, φ)g + ∂ω φ∂ n u ND (g, 0).
Another consequence of these integration by parts formulae is the following alternative expression of the Kohn-Vogelius objective:

(3.2) K(φ, ψ) = ∂Ω f m -u ND (g m , φ) ∂ n u DN (f m , ψ) -g m + ∂ω u DN (f m , ψ) -φ ψ -∂ n u ND (g m , φ) . Remark 3.2.
Let us emphasize that we use the usual abusive notation since the boundary integrals are to be understood as duality product between a Dirichlet trace in H 1/2 and a Neumann trace in H -1/2 .

3.2.

Properties of the Kohn-Vogelius objective K. The main properties of the functional defined by (2.5) can be summarized as follow.

Proposition 3.3 (Properties of the Kohn-Vogelius objective). We have the following statements.

1. The functional K is convex, positive and

inf K(φ, ψ) ; (φ, ψ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω) = 0. 2. The functional K is quadratic with gradient DK(φ, ψ) • [δφ, δψ] = ∂ω [∂ n u ND (∂ n u DN (f m , ψ), φ) -ψ] δφ + u DN u ND (g m , φ), ψ -φ δψ,
and constant Hessian

D 2 K(φ, ψ) • ([δφ 1 , δψ 1 ], [δφ 2 , δψ 2 ]) = u DN (0, δψ 1 ) -u ND (0, δφ 1 ), u DN (0, δψ 2 ) -u ND (0, δφ 2 ) .
The first point of this statement is proved in [

(f n m , ψ * n ) = ∇u ND (g m , φ * n )
. Then, we have

K(φ * n , ψ * n ) = 1 2 |u DN (f m , ψ * n ) -u DN (g m , φ * n )| 2 H 1 (Ω\ω) = 1 2 |u DN (f m , ψ * n ) -u DN (f n m , ψ * n )| 2 H 1 (Ω\ω) = 1 2 |u DN (f m -f n m , 0)| 2 H 1 (Ω\ω) C f m -f n m 2 H 1/2 (∂Ω) -→ n→∞ 0,
which concludes the proof. 2. An elementary computation shows that

K(φ + δφ, ψ + δψ) = 1 2 |u DN (f m , ψ) -u ND (g m , φ) + u DN (0, δψ) -u ND (0, δφ)| 2 H 1 (Ω\ω) = K(φ, ψ) + u DN (f m , ψ) -u ND (g m , φ), u DN (0, δψ) -u ND (0, δφ) + 1 2 |u DN (0, δψ) -u ND (0, δφ)| 2 H 1 (Ω\ω)
We identify the linear and quadratic terms in this expression (that are continuous thanks to (3.1)) and obtain a first expression of its derivatives

DK(φ, ψ) • [δφ, δψ] = u DN (f m , ψ) -u ND (g m , φ), u DN (0, δψ) -u ND (0, δφ)
and the quadratic form

D 2 K(φ, ψ) • ([δφ, δψ], [δφ, δψ]) = |u DN (0, δψ) -u ND (0, δφ)| 2 H 1 (Ω\ω) ,
and finally by polarity

D 2 K(φ, ψ) • ([δφ 1 , δψ 1 ], [δφ 2 , δψ 2 ]) = u DN (0, δψ 1 ) -u ND (0, δφ 1 ), u DN (0, δψ 2 ) -u ND (0, δφ 2 ) .
However, this writing is not directly useful since the dependency in (δφ, δψ) is not completely explicit. We therefore expand the expression of the gradient

DK(φ, ψ) • [δφ, δψ] = u DN (f m , ψ), u DN (0, δψ) -u DN (f m , ψ), u ND (0, δφ) -u ND (g m , φ), u DN (0, δψ) + u ND (g m , φ), u ND (0, δφ) ,
and use the properties of the maps u DN and u DN given in Proposition 3.1 to express each of the previous expressions as integrals over the inner boundary.

We first consider products of the same map. One the one hand, we have

u DN (f m , ψ), u DN (0, δψ) = ∂Ω f m ∂ n u DN (0, δψ) + ∂ω u DN (f m , ψ)δψ = - ∂ω u DN (f m , 0)δψ + ∂ω u DN (f m , ψ)δψ = ∂ω [u DN (f m , ψ) -u DN (f m , 0)] δψ = ∂ω u DN (0, ψ)δψ.
On the other hand, we get in a similar manner

u ND (g m , φ), u ND (0, δφ) = ∂Ω u ND (0, δφ)g m + ∂ω δφ ∂ n u ND (g m , φ) = - ∂Ω δφ ∂ n u ND (g m , 0) + ∂ω δφ ∂ n u ND (g m , φ) = ∂Ω δφ ∂ n [u ND (g m , φ) -u ND (g m , 0)] = ∂Ω δφ ∂ n u ND (0, φ).
Then, we consider the mixed products

u DN (f m , ψ), u ND (0, δφ) = ∂Ω u ND (0, δφ)∂ n u DN (f m , ψ) + ∂ω ψδφ = - ∂ω ∂ n u ND (∂ n u DN (f m , ψ), 0)δφ + ∂ω ψδφ = ∂ω [ψ -∂ n u ND (∂ n u DN (f m , ψ), 0)] δφ,
and in the very same manner

u ND (g m , φ), u DN (0, δψ) = ∂Ω u ND (g m , φ)∂ n u DN (0, δψ) ∂ω φδψ = ∂ω [φ -u DN (u ND (g m , φ), 0)] δψ.
Gathering the terms, we obtain the announced expression of the gradient.

Remark 3.4. As seen from the previous proof, we have

D 2 K(φ, ψ) • ([δφ 1 , δψ 1 ], [δφ 2 , δψ 2 ]) = DK(δφ 1 , δψ 1 ).([δφ 2 , δψ 2 ]) = DK(δφ 2 , δψ 2 ).([δφ 1 , δψ 1 ]).
Thus, the Hessian does not depend on (φ, ψ). Especially, it holds

D 2 K(φ, ψ) • ([δφ 1 , δψ 1 ], [δφ 2 , δψ 2 ]) = ∂ω [∂ n u ND (∂ n u DN (0, δψ 1 ), δφ 1 ) -δψ 1 ] δφ 1 + u DN u ND (0, δφ 2 ), δψ 2 -δφ 2 δψ 2 .
4. Numerical resolution of the data completion problem.

4.1. Newton scheme. As mentioned previously, a regularization process is needed for such a severely ill-posed inverse problem. We use here the standard Tikhonov regularization as studied, e.g., in [START_REF] Caubet | On the data completion problem and the inverse obstacle problem with partial Cauchy data for Laplace equation[END_REF]. Therefore, we consider a nonnegative real number ε and introduce the regularized Kohn-Vogelius cost functional

K ε : H 1/2 (∂ω) × H -1/2 (∂ω) → R defined by (4.1) K ε (φ, ψ) = K(φ, ψ) + εT (φ, ψ),
where the regularizing term T (φ, ψ) is defined by

T (φ, ψ) := 1 2 |u ND (0, φ)| 2 H 1 (Ω\ω) + |u DN (0, ψ)| 2 H 1 (Ω\ω) + ∂ω φ 2 = 1 2 ∂ω φ + ∂ n u ND (0, φ) φ + u DN (0, ψ)ψ.
Lemma 4.1. The regularizing functional T is non-negative with the unique minimizer (0, 0).

Proof. Assume that (φ, ψ) is such that |u ND (0, φ)| 2 H 1 (Ω\ω) + |u DN (0, ψ)| 2 H 1 (Ω\ω) = 0.
Then, u ND (0, φ) and u DN (0, ψ) are constant functions. By exploiting the boundary conditions on ∂ω, we check that φ is constant and ψ = 0 on ∂ω. Hence φ = ψ = 0 is a minimizer and the proof is complete.

Remark 4.2. This regularization permits to obtain the existence of a unique minimizer, for all ε > 0, even for not compatible data thanks to the gained of coerciveness and, in case of compatible data, the convergence of these minimizers towards the exact solution when ε goes to 0 (see, e.g., [16, Proposition 2.5 and Theorem 2.6]).

Remark 4.3. Notice that the usual regularizing term would have been the term (φ, ψ) 2 H 1/2 (∂ω)×H -1/2 (∂ω) , but the two quantities are equivalent as claimed in the following lemma. Lemma 4.4. There exist two non-negative constants C 1 and C 2 such that

C 1 ψ 2 H -1/2 (∂ω) |u DN (0, ψ)| 2 H 1 (Ω\ω) C 2 ψ 2 H -1/2 (∂ω)
and

C 1 φ 2 H 1/2 (∂ω) |u ND (0, φ)| 2 H 1 (Ω\ω) + ∂ω φ 2 C 2 φ 2 H 1/2 (∂ω) ,
for all (φ, ψ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω).

Proof. Firstly, thanks to the continuity of the trace operator and thanks to Proposition 3.1 (see equation (3.1)), there exists two non-negative constants C 1 and C 2 such that

C 1 ψ 2 H -1/2 (∂ω) u DN (0, ψ) 2 C 2 ψ 2 H -1/2 (∂ω) ,
for all ψ ∈ H -1/2 (∂ω). Since u DN (0, ψ) = 0 on ∂Ω, Poincaré's inequality ensures that the semi-norm is equivalent to the H 1 -norm. Hence, we obtain the following inequalities (for others constants still denoted by C 1 and C 2 )

C 1 ψ 2 H -1/2 (∂ω) |u DN (0, ψ)| 2 H 1 (Ω\ω) C 2 ψ 2 H -1/2 (∂ω) .
Secondly, using again (3.1) and the continuous embedding H 1/2 (∂ω) → L 2 (∂ω), there exists a non-negative constant C 3 such that

|u ND (0, φ)| 2 H 1 (Ω\ω) + ∂ω φ 2 C 3 φ 2 H 1/2 (∂ω) ,
for all φ ∈ H 1/2 (∂ω).

Let us now show that there exists a non-negative constant C 4 such that

C 4 u ND (0, φ) 2 H 1 (Ω\ω) |u ND (0, φ)| 2 H 1 (Ω\ω) + ∂ω φ 2 ,
for all φ ∈ H 1/2 (∂ω), which will complete the proof by using the continuity of the trace operator. To this end, let us proceed by contradiction assuming that, for all n ∈ N, there exists φ n ∈ H 1/2 (∂ω) such that

v n H 1 (Ω\ω) = 1 and 1 n > |v n | 2 H 1 (Ω\ω) + ∂ω v 2 n ,
where we have set

v n := u ND (0, φ) u ND (0, φ) H 1 (Ω\ω)
. Hence, ∇v n H 1 (Ω\ω) -→ n→+∞ 0 and

∂ω v 2 n -→ n→+∞ 0. Furthermore, (v n ) n is bounded in H 1 (Ω\ω) which is compactly
embedded in L 2 (Ω\ω) and then, up to a subsequence, there exists v ∈ H 1 (Ω\ω) such that

v n v in H 1 (Ω\ω) and v n → v in L 2 (Ω\ω).
Since ∇v n → 0 in L 2 (Ω\ω), we conclude that v n → v in H 1 (Ω\ω) and ∇v = 0. Therefore, v is constant and

∂ω v 2 n -→ 0 = ∂ω v 2 = v 2 |∂ω| . Thus, v = 0 which contradicts the fact that v n = 1 -→ v = 1.
We now compute the gradient of the objectives with respect to (φ, ψ).

Proposition 4.5. The derivatives of T with respect to (φ, ψ) are

DT (φ, ψ) • [δφ, δψ] = ∂ω δφ φ + ∂ n u ND (0, φ) + u DN (0, ψ)δψ and D 2 T (φ, ψ) • [(δφ 1 , δψ 1 ), (δφ 2 , δψ 2 )] = ∂ω δφ 1 δφ 2 + ∂ n u ND (0, δφ 2 ) + u ND (0, δψ 2 )δψ 1 .
Proof. We expand the quadratic quantity T and follow the proof of Proposition 3.3 to get the announced result.

Remark 4.6. The first idea to build a numerical scheme for minimizing the objective is to use a descent method. This leads to a sequence (φ n , ψ n ) by the update rule

φ n+1 ψ n+1 = φ n ψ n + τ n+1 d n ,
where τ n+1 is a descent step and where the descent direction d n is naturally chosen as the formal anti-gradient:

d n = - ∂ n u ND ∂ n u DN (f m , ψ n ), φ n -ψ n u DN u ND (g m , ψ n ), ψ n -φ n + ε φ n + ∂ n u ND (0, φ n ) u DN (0, ψ n ) .
In fact, d n is not the gradient of K. Indeed, the true gradient should be computed with respect to the scalar product on H 1/2 (∂ω) × H -1/2 (∂ω) and not with respect to L 2 (∂ω) × L 2 (∂ω). A nice way to see that is that this updates is not in the right spaces: while (φ n , ψ n ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω), the update d n lies in H -1/2 (∂ω) × H 1/2 (∂ω).

The true gradient is much more complex to compute at least at the continuous level. Therefore, we will not consider the gradient method here and directly solve the linearized equation. This is the idea of Newton method.

Proposition 4.7. The Newton update (δφ, δψ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω) for the regularized Kohn-Vogelius objective K ε (φ n , ψ n ) is given by the linear system

       ∂ n u ND ∂ n u DN (0, δψ), δφ -δψ + ε δφ + ∂ n u ND (0, δφ) = ψ n -∂ n u ND ∂ n u DN (f m , ψ n ), φ n -ε φ n + ∂ n u ND (0, φ n ) , u DN u ND (0, δφ), δψ -δφ + εu DN (0, δψ) = φ n -u DN u ND (g m , ψ n ), ψ n -εu DN (0, ψ n ).
Proof. The Newton scheme leads to the sequence (φ n , ψ n ) such that the increment

(δφ, δψ) = (φ n+1 , ψ n+1 ) -(φ n , ψ n )
satisfies the weak problem

(D 2 K + εD 2 T )(φ n , ψ n ) • ([δφ, δψ], [h, ]) = -(DK + εDT )(φ n , ψ n ) • [h, ] for all [h, ] ∈ H 1/2 (∂ω) × H -1/2 (∂ω).
To solve this equation in the unknowns (δφ, δψ), we need to rewrite the left hand side in order to make explicit the dependency on [h, ]. We proceed as for the gradient: we first expand u DN (0, δψ) -u ND (0, δφ), u DN (0, ) -u ND (0, h) = u DN (0, δψ), u DN (0, ) + u ND (0, δφ), u ND (0, h) -u ND (0, δφ), u DN (0, ) -u DN (0, δψ), u ND (0, h) and then compute each term:

u DN (0, δψ), u DN (0, ) = ∂ω u DN (0, δψ) , u ND (0, δφ), u ND (0, h) = ∂ω ∂ n u ND (0, δφ)h, u ND (0, δφ), u DN (0, ) ' = ∂ω δφ -u DN u ND (0, δφ), 0 , u DN (0, δψ), u ND (0, h) = ∂ω δψ -∂ n u ND ∂ n u DN (0, δψ), 0 h. The weak problem reads: find (δφ, δψ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω) such that ∂ω u DN u ND (0, δφ), δψ -δφ + εu DN (0, δψ) + ∂ n u ND ∂ n u DN (0, δψ), δφ -δψ + ε δφ + ∂ n u ND (0, δφ) h = - ∂ω ∂ n u ND ∂ n u DN (f m , ψ n ), φ n -ψ n + ε φ n + ∂ n u ND (0, φ n ) h - ∂ω u DN u ND (g m , ψ n ), ψ n -φ n + εu DN (0, ψ n )
for all (h, ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω). This implies the assertion. 4.2. Implementation. Our approach to determine the Cauchy data of the functions u DN (f m , ψ) and u ND (g m , φ) relies on a system of boundary integral equations arising from the direct formulation based on Green's fundamental solution.

In this section, in order to simplify the formulae, we use the following notation: Σ := ∂Ω and Γ := ∂ω.

Assuming that u ∈ H 1 (Ω\ω) solves Laplace's equation, then Green's representation formula implies the relation

u(x) = Σ∪Γ G(x, y) ∂u ∂n (y) - ∂G(x, y) ∂n y u(y) dσ y , x ∈ Ω\ω,
where the Green function G is given by

G(x, y) =        - 1 2π log |x -y|, if d = 2, 1 4π|x -y| , if d = 3.
Using the jump properties of the layer potentials, we obtain the direct boundary integral formulation of the problem

(4.2) u(x) = Γ∪Σ G(x, y) ∂u ∂n (y)dσ y + 1 2 u(x)- Γ∪Σ ∂G(x, y) ∂n y u(y)dσ y , x ∈ Γ∪Σ.
We restrict ourselves to the two-dimensional situation d = 2.

Let A, B ∈ {Γ, Σ} which represent the boundaries. Then the expression (4.2) includes the single layer operator defined by

(4.3) S : C(A) → C(B), S AB (ρ)(x) = - 1 2π A log |x -y| ρ(y)dσ y ,
and the double layer operator defined by

(4.4) D : C(A) → C(B), D AB (ρ)(x) = 1 2π A x -y, n y |x -y| ρ(y)dσ y ,
with the densities ρ being the Cauchy data of u at the boundary A. The equation (4.2) in combination with (4.3) and (4.4) amounts to the Dirichlet-to-Neumann map, which is given by the following system of integral equations

S ΓΓ S ΣΓ S ΓΣ S ΣΣ ∂ n u| Γ ∂ n u| Σ = 1 2 I + D ΓΓ D ΣΓ D ΓΣ 1 2 I + D ΣΣ u| Γ u| Σ .
Reordering this system of boundary integral equations yields the missing Cauchy data of the solution u DN (f m , ψ) ∈ H 1 (Ω\ω) of problem (2.3) by (4.5)

1 2 I + D ΓΓ -S ΣΓ -D ΓΣ S ΣΣ u| Γ ∂ n u| Σ = S ΓΓ -D ΣΓ -S ΓΣ 1 2 I + D ΣΣ ψ f m and of the solution u ND (g m , φ) ∈ H 1 (Ω\ω) of Problem (2.4) by (4.6) S ΓΓ -D ΣΓ -S ΓΣ 1 2 I + D ΣΣ ∂ n u| Γ u| Σ = 1 2 I + D ΓΓ -S ΣΓ -D ΓΣ S ΣΣ φ g m .
The boundary integral operators on the left hand side of the systems (4.5) and (4.6) of boundary integral equations are continuous and satisfy a Gårding inequality with respect to L 2 (Γ) × H -1/2 (Σ) and H -1/2 (Γ) × L 2 (Σ), respectively, provided that diam(Ω) < 1. Since injectivity follows from potential theory, these systems of integral equations are uniquely solvable according to the Riesz-Schauder theory. With the help of the solutions given by (4.5) and (4.6), we can compute the Kohn-Vogelius functional (3.2) and its gradient and Hessian. The same holds true for the regularized functional K ε . The next step towards the solution of the boundary value problem is the numerical approximation of the integral operators included in (4.5) and (4.6), respectively, which first requires the parametrization of the integral equations. To that end, we insert parameterizations γ : [0, 2π] → Γ and σ : [0, 2π] → Σ of the boundaries. For the approximation of the unknown Cauchy data, we use the collocation method based on trigonometric polynomials. Applying the trapezoidal rule for the numerical quadrature and the regularization technique along the lines of [START_REF] Kress | Linear integral equations[END_REF] to deal with the singular integrals, we arrive at an exponentially convergent boundary element method provided that the data and the boundaries and thus the solution are arbitrarily smooth.

Since the Hessian of the functional K ε (φ, ψ) does not depend on the argument (φ, ψ), the Newton scheme converges in one iteration. Especially, the Hessian is positive definite since the functional under consideration is convex. We can thus use the conjugate gradient method [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] to compute the Newton update. Therefore, we can efficiently determine the functions (φ, ψ) ∈ H 1/2 (Γ) × H -1/2 (Γ) such that it holds u ND (g m , φ) = u DN (f m , ψ) for the current boundary Γ = ∂ω. In our first example, we consider the ellipse with semi-axis h x = 0.25 and h y = 0.15 as inclusion ω, compare the left plot of Figure 4.1. In the second example, we have a potato shaped inclusion ω, compare the right plot of Figure 4.1. We moreover prescribe the Cauchy data (f m , g m ) of the harmonic function u(x, y) = x 2 -y 2 . The approximation errors for the different geometries and regularization parameters are found in Table 4.1, where the respective Sobolev norms are computed by means of appropriately weighted Fourier series. In particular, we observe that the best choice for the regularization parameter is ε = 0.01 since therefore the approximation errors are the smallest ones. We like to stress that an increase of the number N of boundary elements does not increase accuracy, which reflects the severe ill-posedness of the problem under consideration. We like to mention that, adding noise to the Neumann data g m as described in the next section, does not significantly change the results found in Table 4.1. The approximation errors for the data completion approach in dependence of the regularization parameter.

regularization ε = 0.1 ε = 0.05 ε = 0.01 ε = 0.005 ε = 0.001 parameter ellipse φ-φ N H 1/2 (Γ) φ H 1/2 (Γ) 1.7256 • 10 -1 1.5610 • 10 -1 1.4621 • 10 -1 2.3244 • 10 -1 2.3302 • 10 -1 ψ-ψ N H -1/2 (Γ) ψ H -1/2 (Γ)

Resolution of the inverse obstacle problem with the trial method.

We now aim to numerically solve the inverse obstacle problem (1.1). The general idea is to use the previous data completion step in order to reconstruct u| ∂ω and ∂ n u| ∂ω on an approximation ∂ω of the real inclusion and then to use the so-called trial method in order to update the shape of the inclusion.

We first present the trial method and then illustrate the efficiency of our method with some numerical simulations.

5.1. Background and motivation. The trial method is a fixed-point type iterative method, which is well-known from the solution of free boundary problems (see, e.g., [START_REF] Beurling | On free-boundary problems for the Laplace equation[END_REF][START_REF] Flucher | Bernoulli's free-boundary problem, qualitative theory and numerical approximation[END_REF][START_REF] Tiihonen | On fixed point (trial) methods for free boundary problems[END_REF] and the references therein). In the context of inverse problems, it has been used for example in [START_REF] Serranho | A hybrid method for inverse obstacle scattering problems[END_REF].

We shall assume in the following that the domain ω is starlike. Hence, we can represent the inclusion's boundary ∂ω by a parametrization γ : [0, 2π] → R 2 in polar coordinates, that is ∂ω = γ(s) = r(s)e r (s) ; s ∈ [0, 2π] ,

where e r (s) = cos(s), sin(s) denotes the unit vector in the radial direction. The radial function r(s) is supposed to be a positive function in C per ([0, 2π]), where

C per ([0, 2π]) = r ∈ C([0, 2π]) ; r(0) = r(2π) , such that dist(∂Ω, ∂ω) > 0.
The trial method to solve the conductivity problem problem (1.1) requires an update rule. Suppose that the actual void's boundary is ∂ω k . Then the data completion problem yields a state u k which satisfies

   ∆u k = 0 in Ω\ω k , u k = f m on ∂Ω, ∂ n u k = g m on ∂Ω.
The new boundary ∂ω k+1 is now determined by moving the old boundary into the radial direction, which is expressed by the update rule (5.1)

γ k+1 = γ k + δr k e r .
The computation of the update function δr k is the topic of the next section.

Update rule.

The update function δr k ∈ C per ([0, 2π]) should be constructed in such a way that the desired homogeneous Dirichlet boundary condition will be (approximately) satisfied at the new boundary ∂ω k+1 , i.e.,

u k • γ k+1 ! = 0 on [0, 2π],
where u k is assumed to be smoothly extended into the exterior of Ω\ω k if required.

The traditional update rule is obtained by linearizing u k • (γ k + δr k e r ) with respect to the update function δr. This yields the equation

u k • γ k+1 ≈ u k • γ k + ∂u k ∂e r • γ k δr k .
We decompose the derivative of u k in the direction e r into its normal and tangential components:

∂u k ∂e r = ∂u k ∂n e r , n + ∂u k ∂t e r , t on ∂ω k . Hence, defining F (δr k ) = u k • γ k + ∂u k ∂e r • γ k δr k , we arrive at the update equation (5.2) F (δr k ) = u k • γ k + ∂u k ∂n • γ k e r , n + ∂u k ∂t • γ k e r , t δr k ! = 0.
Remark 5.1. We mention that the solution of the data completion problem according to Section 4 immediately yields the quantities

u k | ∂ω k and ∂ n u k | ∂ω k . Since u k | ∂ω k is expressed in terms of trigonometric polynomials, it is also straightforward to compute ∂u k ∂t • γ k (s) = 1 |γ k (s)| ∂u k ∂s (s).
Consequently, all terms required in (5.2) are available.

Since the trial method is fixed-point type iterative method, its convergence is only ensured if it is a contraction. We refer the reader to [START_REF] Tiihonen | On fixed point (trial) methods for free boundary problems[END_REF] for sufficient conditions on its convergence.

Discretization of the sought boundary.

For the numerical computations, we discretize the radial function r n k associated with the boundary ∂ω k by a finite Fourier series according to

(5.3) r n k (s) = a 0 + n-1 =1 a cos( s) + b sin( s) + a n cos(ns).
This obviously ensures that r n k is always an element of C per ([0, 2π]). To determine the update function δr n k , represented likewise by a finite Fourier series, we insert the N 2n equidistantly distributed points s = 2π /N into the update equations (5.2):

F (δr n k ) ! = 0 in all the points s 1 , . . . , s N .
This is a discrete least-squares problem which can simply be solved by the normal equations.

Numerical results.

We shall illustrate our approach by some numerical test examples. We choose Ω as the ball of radius 0.4, which is centered in 0. Knowing the inclusion ω, we can compute the current measurements g m from a given voltage distribution f m . For the specific voltage distribution f m (x, y) = y| ∂Ω , we try to reconstruct the inclusion ω in the following by using the synthetic data (f m , g m ).

The settings for the data completion algorithm are as in Subsection 4.3, i.e., the regularization parameter ε for the data completion functional (4.1) is set to 0.01, which turned out to be the best choice. We use 200 boundary elements per boundary for the discretization of the boundary integral equations in (4.5) and (4.6), respectively.

The reconstruction is performed by the trial method detailed previously. Precisely, the sought inclusion is represented by a Fourier series with 10 terms, i.e., we have n = 5 in (5.3). The initial guess is always a ball of radius 0.3 (always indicated via the red boundary in the subsequent figures) and we stop the trial method after 20 iterations, where the damping factor h = 0.2 is used. This means that the update from (5.2) is multiplied with h before computing the new iterate in accordance with (5.1).

In our first example, we consider again the ellipse with semi-axis h x = 0.25 and h y = 0.15 as inclusion ω, compare the left plot of Figure 4.1. In the second example, we have again the potato shaped inclusion ω, compare the right plot of Figure 4.1. We add 1% and 5% noise, respectively, to the respective synthetic current measurement g m , where noise g noise is generated as a random vector with on [-1, 1]/ √ n uniformly distributed random variables, such that g δ m := g m + δ g m L 2 (Ω) g noise . In order to better understand the reconstruction algorithm, we draw 25 realizations with noise level 5 % for both obstacles. They can be found in Figure 5.3 for the ellipse inclusion on the left-hand side and for the potato inclusion on the right-hand side. Since we are dealing with starshaped domains, we can compute the mean of the realizations by taking the mean of the respective Fourier coefficients. The mean shape is plotted in green. As one readily infers, the reconstructions vary around the mean shapes which basically coincides with the sought inclusions. We shall finally explore the hard case of a non-convex, small obstacle. The results are presented in Figure 5.4. It underlines the fact that the reconstruction is less performant: the right concavity is not perfectly recovered. Moreover, the reconstruction becomes worse when the void becomes smaller.

Conclusion.

In the present article, we solved the data completion problem for Laplace's equation by the minimization of the Kohn-Vogelius functional, which has been regularized with respect to the associated energy norm. The minimization of the regularized Kohn-Vogelius functional has been performed by a Newton scheme, which results in a direct solver. By employing a collocation method based on trigonometric polynomials, we arrive at a very efficient numerical method for the data completion problem. We then combined this data completion algorithm with the trial method in order to solve an inverse obstacle problem. It updates a given inclusion such that the desired inclusion's homogeneous Dirichlet boundary condition is approximately satisfied. This yields an iterative method for the detection of inclusions in electric impedance tomography.

(1. 1 )

 1 Find a set ω * ∈ D and a solution u ∈ H 1 Ω\ω * of the following overdetermined boundary value problem: Ω\ω * , u = f m on ∂Ω, ∂ n u = g m on ∂Ω, u = 0 on ∂ω * , where D = {ω Ω, ∂ω is Lipschitz and Ω\ω is connected} and n is the outer unit normal to Ω\ω. In other words, we want to detect an inclusion ω * characterized by a homogeneous Dirichlet boundary condition from the knowledge of the Dirichlet and Neumann boundary conditions on the exterior boundary ∂Ω. Notice that a only a single pair of Cauchy data is used.

2 .

 2 Introduction of general notations and of the considered problems. 2.1. Introduction of the general notations. For a bounded open set Ω of R d (d ∈ N \ {0}) with a (piecewise) Lipschitz boundary ∂Ω, we precise that the notation Ω u means Ω u(x)dx which is the classical Lebesgue integral. Moreover we use the notation ∂Ω u to denote the boundary integral ∂Ω u(x)ds(x), where ds represents the surface Lebesgue measure on the boundary. We also introduce the exterior unit normal n of the domain Ω and ∂ n u will denote the normal derivative of u.

2. 2 .

 2 Introduction of the data completion problem. We recall that we consider in this article a bounded connected Lipschitz open set Ω of R d (d = 2 or d = 3) and an inclusion ω ∈ D, with D = {ω Ω, ∂ω is Lipschitz and Ω\ω is connected}. Let us consider some nontrivial Cauchy data

4. 3 .

 3 Numerical validation. We shall validate the data completion approach by some numerical test examples. We choose Ω as the ball of radius 0.4, which is centered in 0. If we take the Cauchy data of a known harmonic function u at the outer boundary ∂Ω as data (f m , g m ), the desired Cauchy data (φ, ψ) ∈ H 1/2 (∂ω) × H -1/2 (∂ω) at the inclusion's boundary ∂ω are just (u| ∂ω , ∂ n u| ∂ω ). The settings for the data completion algorithm are as follows. The regularization parameter ε for the data completion functional (4.1) is set to ε = 0.1, ε = 0.01, or ε = 0.001. We use N = 200 boundary elements per boundary for the discretization of the boundary integral equations in (4.5) and (4.6), respectively. Hence, the computed Cauchy data (φ N , ψ N ) at the inclusion's boundary ∂ω are represented by N = 200 boundary elements each.

Fig. 4 . 1 .

 41 Fig. 4.1. The ball with ellipse inclusion (left-hand side) and with potato inclusion (right-hand side).

Fig. 5 . 1 .

 51 Fig. 5.1. Iterates (in blue) and the final reconstructions (in green) for the ellipse inclusion with 1% noise (left-hand side) and with 5% noise (right-hand side).The final reconstructions, indicated via the green boundary, can be found in Figure5.1 for the first example and in Figure5.2 for the second example; on the left-hand side the reconstructions for the noise level 1% are shown and on the righthand side the reconstructions for the noise level level 5% are shown. The intermediate iterates are indicated by the blue boundaries.

Fig. 5 . 2 .

 52 Fig. 5.2. Iterates (in blue) and the final reconstructions(in green) for the potato inclusion with 1% noise (left-hand side) and with 5% noise (right-hand side).

Fig. 5 . 3 .

 53 Fig. 5.3. 25 reconstructions and mean shape (in green) in case of the ellipse inclusion (left-hand side) and in case of the potato inclusion (right-hand side).

Fig. 5 . 4 .

 54 Fig. 5.4. Iterates (in blue) and the final reconstructions (in gree) for different scalings of a kite shaped void with 1% noise.

  When data is compatible, the unknown u is uniquely determined by either its Dirichlet trace or its Neumann trace on the inner boundary ∂ω. The property of being compatible is not open: indeed the set of compatible couples is dense inH 1/2 (∂Ω) × H -1/2 (∂Ω). Hence, if a given pair (f m , g m ) is not compatible, we may approximate it by a sequence of compatible data. Typical density results have been stated by Andrieux et al. in [1, Lemma 2.1] or Fursikov in

	∂Ω) is said compatible if
	the Cauchy problem (2.1) has a (necessarily unique) solution.

  16, Proposition 2.2] but we recall the proof for reader's convenience. ) = 0. Let us now focus on the non-compatible case. Thanks to the density lemma 2.2, we can approximate f m by a sequence (f n m ) n in a way that the pairs (f n m , g m ) n are compatibles for all n ∈ N. For each n, consider (φ * n , ψ * n ) to be the minimizer of the Kohn-Vogelius function for the data (f n m , g m ) which implies that ∇u DN

	Proof. We prove each statement.

1. Convexity and positiveness are obvious. To prove that inf (φ,ψ) K(φ, ψ) = 0, we have to consider two cases. If the pair (f m , g m ) is compatible, we consider the solution u ex of the Cauchy problem (2.1) and we define φ * := u ex | ∂ω and ψ * := ∂ n u ex | ∂ω and then obtain K(φ * , ψ *

  1.3689 • 10 -2 1.3586 • 10 -2 7.7903 • 10 -3 1.2139 • 10 -1 1.2085 • 10 -1 10 -2 2.6365 • 10 -2 3.5529 • 10 -2 3.6279 • 10 -2 3.7131 • 10 -2

	potato	φ-φ N H 1/2 (Γ) φ H 1/2 (Γ) ψ-ψ N H -1/2 (Γ) ψ H -1/2 (Γ)	8.0964 • 10 -2 5.0993 • 10 -2 8.1978 • 10 -2 8.4939 • 10 -2 8.7713 • 10 -2 4.2829 •

Table 4 . 1

 41


Marc Dambrine

Helmut Harbrecht

Bénédicte Puig

Computing quantities of

Keywords: 60G35, 65N75, 65N99 random domain, second order shape sensitivity, low-rank approximation

come

Introduction

Many practical problems in engineering lead to boundary value problems for an unknown function which needs to be computed to obtain a real quantity of interest. In structural mechanics, the equations of elasticity are usually considered and solved to compute the leading mode of a structure or its compliance. Usually, the input parameters of the model like the geometry or the physical coefficients (typically the value of the Young modulus or Poisson ratio) are not perfectly known. It is therefore important to take these uncertainties into account.

In this work, we consider uncertainties in the geometric definition of the domain. Such uncertainties are motivated by tolerances in the fabrication process or by a damaged boundary during the life of a mechanical device. Manufactured devices are close to a nominal geometry but differ of course from its mathematical definition. Since we are motivated by tolerances, we can make the crucial assumption of the smallness of the random perturbations. This means that we assume that the non-dimensional ratio ε between the characteristic size of the perturbation and the characteristic size of the domain is small. By identifying domains with their boundary, domains close to the nominal domain D can be seen as a normal perturbation of the nominal boundary ∂D. In that case, the random domain D(ω) can be defined thanks to a real valued random field X over ∂D according to ∂D(ω) = {x + X(x, ω)n(x); x ∈ ∂D}, where n(x) denotes the outer unit normal field at a point x in ∂D. In order to take into account the uncertain geometrical definition in the numerical simulation, we have to incorporate the randomness of the computational domain into the underlying model equations. As a consequence, the quantities of engineering interest are also random. We address thus the following problem: given a complete probabilistic description of the random perturbation of the nominal boundary, compute as much information as possible to the distribution of the quantity of interest.

The most common approach to study boundary value problems with stochastic inputs is the Monte-Carlo method. In many situations, this approach is easy to implement since it requires only a large number of samples. However, for boundary values problems on random domains, each sample means a new domain and thus a new mesh, the building of new mass and stiffness matrices, etc. All these steps are mandatory to compute the quantity of interest. Therefore, the Monte-Carlo method is extremely costly and not so easy to implement in our context. This work is a contribution to the development of cheap and deterministic numerical methods to recover statistical information on the distribution of the output quantity of interest. The main advantage of the approach we propose here is to use a single mesh and to build once for all the matrices.

The smallness assumption formalized in the sequel allows us to use a sensitivity analysis with respect to geometrical perturbations of the boundary. We apply shape calculus to perform a second order sensitivity analysis and compute asymptotic expansions of the moments of the output quantity's distribution with respect to the smallness parameter ε of the random perturbation. In particular, we obtain in this work a third order expansion for the expectation and a fourth order expansion for the variance. It turns out that the leading coefficients of these asymptotic expansions do only depend on the autocovariance of the random field which defines the random perturbation. We therefore use the associated integral operator to compute a low-rank approximation of the random field of the type

X(x, ω) = N i=1 α i (ω)f i (x).
For such a random field, we derive an analytic expression of the coefficients of the previous asymptotic expansions that are very cheap to compute: only N + 1 boundary values problems are to be solved on the nominal geometry D in the examples we present in this work. Note that 2N + 2 problems are needed when the derivatives are evaluated thanks to an adjoint state.

We mention that random shape functionals have also been considered in [START_REF] Harbrecht | On output functionals of boundary value problems on stochastic domains[END_REF] by means of a first order perturbation analysis. A rather general framework of the first order perturbation analysis for functionals with random input parameters, particularly random domains, has been presented [START_REF] Chernov | First order k-th moment finite element analysis of nonlinear operator equations with stochastic data[END_REF]. To obtain more precise approximations, the present paper is based on a second order perturbation analysis of the random shape functional under consideration. Thus, the shape Hessian will enter the asymptotic expansions which makes the computations much harder. However, the shape Hessian is meanwhile well understood and has been considered for example in [START_REF] Dambrine | On variations of the shape Hessian and sufficient conditions for the stability of critical shapes[END_REF][START_REF] Dambrine | Stability in shape optimization with second variation[END_REF][START_REF] Eppler | Boundary integral representations of second derivatives in shape optimization[END_REF][START_REF] Eppler | Preconditioning the pressure tracking in fluid dynamics by shape Hessian information[END_REF][START_REF] Henrot | Variation et optimisation de formes. Une analyse géométrique[END_REF][START_REF] Novruzi | Structure of shape derivatives[END_REF][START_REF] Schulz | A Riemannian view on shape optimization[END_REF]. In comparison to [START_REF] Harbrecht | On output functionals of boundary value problems on stochastic domains[END_REF], we derive here a more precise expansion of the random shape functional. This fact is also verified by our numerical experiments.

This work is organized as follows. In Section 2, we introduce our main theoretical tool: the shape calculus of order two. We define the shape derivative of a general shape function and present the structure theorem. In order to illustrate the method, we choose two examples which we will focus on in this work: the first eigenvalue of the Dirichlet-Laplace operator and the Dirichlet energy. We provide the expressions of the shape derivatives of these two quantities. Then, in Section 3, we detail our random model and explain how one can in general obtain asymptotic expansions of the moments of the quantity of interest's distribution. In Section 4, we emphasize the role of the autocovariance of the field X to compute the previously derived Taylor expansions and explain how to obtain a low-rank approximation of the field X itself. In Section 5, we make explicit the expansions of Section 3 and obtain expressions which are directly workable for computations. Finally, we validate our theoretical findings in Section 6 by comparing our results with Monte-Carlo simulations.

2 Second order shape calculus

Definitions and structure theorem

Shape calculus was founded by Hadamard a century ago but was really developed from the seventies on with the works of Murat and Simon as well as Sokolowski and Zolesio. It's objective is to provide a differential calculus for functions which depend on the geometry of a domain D. It is achieved through the action of a family of diffeomorphisms acting on a model domain. We refer readers interested in the precise definitions and useful properties to modern books on the subject such as [START_REF] Delfour | Shapes and geometries. Metrics, analysis, differential calculus, and optimization[END_REF][START_REF] Henrot | Variation et optimisation de formes. Une analyse géométrique[END_REF].

It is well-known since Hadamard's work that, in smooth situations, the shape gradient is a distribution supported on the moving boundary and acting on the normal component of the deformation field that is a vector field defined in R d with values in R d . In this work, we consider a smooth vector field in C k,1 (R d , R d ), the space of C k -functions with Lipschitz continuous k-th derivatives. The second order shape derivative has also a specific structure as stated by Pierre and Novruzi in [START_REF] Novruzi | Structure of shape derivatives[END_REF]. We quote their result from the book [16, Theorem 5-9-2]. (i) If D ∈ O k+1 and J is differentiable at 0, then there exists a continuous linear form L on C k (∂D) such that:

DJ (0)ξ = L(ξ • n) for all ξ ∈ C k,1 (R d , R d ).
(ii) If D ∈ O k+2 and J is twice differentiable at 0, then there exists a continuous symmetric bilinear form

B on C k (∂D) × C k (∂D) such that for all (ξ, ζ) ∈ C k,1 (R d , R d ) 2 D 2 J (0)(ξ, ζ) = B(ξ • n, ζ • n) + L((D τ nζ τ ) • ξ τ -∇ τ (ζ • n) • ξ τ -∇ τ (ξ • n) • ζ τ ),
where ∇ τ is the tangential gradient and ξ τ and ζ τ stands for the tangential components of ξ and ζ. Finally, D τ n is the differential of the Gauss map known as the shape operator or Weingarten map.

The so-called shape derivatives are then the shape gradient, usually noted DJ(D) := DJ (0), and the shape Hessian, usually noted D 2 J(D) := D 2 J (0). With respect to this work, it is important to notice that the shape Hessian is reduced to B for normal deformations fields.

Following the structure theorem, let us consider a C 5 -domain D 0 and consider a C 3 -neighbourhood O of D and a twice (locally) differentiable shape function J. If the size of the neighbourhood is small enough, the local inversion theorem shows that the boundary ∂D of any domain D in O can be represented as a graph over ∂D 0 of the form: there is a real valued function ϕ defined on ∂D 0 such that

∂D = {x + ϕ(x)n(x); x ∈ ∂D 0 }. (1)
In particular, one can restrict oneself to normal perturbations of amplitude ϕ. Moreover, for a given function ϕ defined on ∂D 0 , the interior of the set {x + ϕ(x)n(x); x ∈ ∂D 0 } determines a domain D ϕ . Thus, the Taylor formula for the mapping J defined in Theorem 2.1 for θ = ϕn, where we still have denoted smooth extensions to R d by ϕ and n, reads as

J(D ϕ ) = J(D 0 ) + L[J](D 0 )ϕ + 1 2 B[J](D 0 )(ϕ, ϕ) + R 2 (ϕ). (2)
Here, the reminder R 2 is uniformly in ϕ negligible with respect to ϕ 2 . Note that, in many applications, the shape functions are regular and in particular C 3 (see the discussion in [START_REF] Dambrine | On the shape sensitivity of the first Dirichlet eigenvalue for two-phase problems[END_REF]). In that case, the remainder is of order 3 in the amplitude ε of the considered random perturbations.

Examples of shape derivatives

Shape gradients have been computed for many objectives covering many boundary value problems motivated by applications in civil engineering, inverse problems, aeronautics, etc. A lot of examples can be found in the book [START_REF] Sokolowski | Introduction to Shape Optimization[END_REF]. The shape gradient is widely used in level-set methods for shape optimization (see [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF] among may others). Due to their much higher complexity, shape Hessians are not so often applied and the literature about computations with second order shape derivatives is restricted (see, e.g., [START_REF] Eppler | A regularized Newton method in electrical impedance tomography using shape Hessian information[END_REF][START_REF] Eppler | Second-order shape optimization using wavelet BEM[END_REF][START_REF] Harbrecht | A Newton method for Bernoulli's free boundary problem in three dimensions[END_REF][START_REF] Novruzi | Newton's method in shape optimisation: a three-dimensional case[END_REF]). Nevertheless, the shape Hessian underlies the study of stability issues in shape optimization (see [START_REF] Afraites | On second order shape optimization methods for electrical impedance tomography[END_REF][START_REF] Caubet | Instability of an inverse problem for the stationary Navier-Stokes equations[END_REF][START_REF] Caubet | Stability of critical shapes for the drag minimization problem in Stokes flow[END_REF][START_REF] Dambrine | On variations of the shape Hessian and sufficient conditions for the stability of critical shapes[END_REF][START_REF] Dambrine | On the shape sensitivity of the first Dirichlet eigenvalue for two-phase problems[END_REF][START_REF] Dambrine | Stability in shape optimization with second variation[END_REF][START_REF] Eppler | Preconditioning the pressure tracking in fluid dynamics by shape Hessian information[END_REF][START_REF] Hintermüller | An inexact Newton-CG-type active contour approach for the minimization of the Mumford-Shah functional[END_REF][START_REF] Schmidt | Impulse response approximations of discrete shape Hessians with application in CFD[END_REF][START_REF] Simon | Second variations for domain optimization problems[END_REF][START_REF] Yang | A shape Hessian-based boundary roughness analysis of Navier-Stokes flow[END_REF] for examples in imaging, tomography, fluid mechanics, aircraft construction, etc.). For the sake of readability, we present two academic but representative examples for which the expressions of the shape derivatives remains simple.

We need to precise some geometrical definitions. The mean curvature (understood as the sum of the principal curvatures of ∂D) is denoted by H. For a domain D ⊂ R d , we consider its Dirichlet energy E(D) defined as

E(D) = - 1 2 D |∇u D | 2 ,
where u D is the solution of -∆u = 1 in H 1 0 (D) and λ 1 is the first eigenvalue of the Dirichlet-Laplace operator. The shape derivatives of these functionals are well-known (see [16, Sections 5.9.3 and 5.9.6]).

Lemma 2.2 (Expressions of shape derivatives)

If D is C 2 , one has L[E](D).ϕ = - 1 2 ∂D (∂ n u) 2 ϕ; (3a) B[E](D).(ϕ, ϕ) = -∂ n u ϕ, Λ(-∂ n u ϕ) H 1/2 ×H -1/2 + ∂D ∂ n u + 1 2 H(∂ n u) 2 ϕ 2 ; (3b) L[λ 1 ](D).ϕ = - ∂D (∂ n v) 2 ϕ; (3c) B[λ 1 ](D).(ϕ, ϕ) = ∂D 2w(ϕ) ∂ n w(ϕ) + H(∂ n v) 2 ϕ 2 ; ( 3d
)
where Λ : H 1/2 (∂D) → H -1/2 (∂D) is the Dirichlet-to-Neumann map for the domain D defined as

Λ(ϕ) = -∂ n V (ϕ) with V (ϕ) being the solution of -∆V (ϕ) = 0 in D, V (ϕ) = -ϕ on ∂D, ( 4
)
and v is the associated normalized eigenfunction solution in

H 1 0 (D) of -∆v = λ 1 v with v > 0 in D and w(ϕ) is the solution of                  -∆w(ϕ) = λ 1 w(ϕ) -v ∂D (∂ n v) 2 ϕ in D, w(ϕ) = -ϕ∂ n v on ∂D, D v w(ϕ) = 0. (5)
As consequence of these examples and being the general case for smooth domains to the best of our knowledge, there is a function ℓ, defined on ∂D, so that the shape gradient can be written as

L(D).ϕ = ∂D ℓϕ. (6)
This property to be an integral operator with a nice kernel is not true in general for the second order derivative.

3 Asymptotic expansions of the moments

The stochastic model: random graphs over ∂D

For modeling the stochastic perturbations D(.) of D 0 , we introduce a probability space (Ω, F, P) and consider stochastic functions X : ∂D × Ω → R that define D(ω) as the interior of the graph x + X(x, ω)n ∂D 0 (x). In order to keep a pertinent geometrical description of sets, this of course requires that the domains D(ω) remain close to D 0 at least in the L ∞ sense. In fact, to use the expansion (2), we need closeness in the sense of C 2 . Introducing a small parameter ε > 0 and plugging perturbations of the form x + εX(x, ω)n ∂D 0 (x) into (2), we obtain

J(D(ω)) = J(D 0 ) + εL(X(ω)) + ε 2 2 B(X(ω), X(ω)) + R 2 (εX(ω)). (7)
We make the assumption that X is uniformly (in ω) bounded in the norm in the spatial variable where shape differentiability holds at the second order (typically here the C 2 -norm), i.e., X is a member of the Bochner space L ∞ (Ω, C 2 (∂D)). Then, there exists

C > 0 such that |R 2 (εX(ω))| = O(ε 2 ).
In most applications, the quantity J(D(ω)) is defined through the solution u D(ω) of a boundary value problem set in D(ω). Let B be a fixed large ball such that D(ω) ⊂ B for all ω ∈ Ω. To fix the idea, consider the Dirichlet problem for a second order elliptic partial differential operator L: u D is the solution of Lu = f in H 1 0 (D). We still denote u D for its extension by 0 to H 1 0 (B). Let us remark that the mapping ω → u D(ω) is the composition of the mappings ω → X(ω) ∈ C 2 (∂D) and X → u D(X) ∈ H 1 0 (B) defined through [START_REF] Afraites | On second order shape optimization methods for electrical impedance tomography[END_REF]. Since the second mapping X → u D(X) ∈ H 1 0 (B) is continuous, it is measurable. Moreover, it is strongly measurable by the Orlicz-Pettis theorem since both C 2 (∂D) and H 1 0 (B) are separable (we used the fact that D is bounded so that ∂D is compact).

General preliminary results

In view of (7), we derive the asymptotic expansions of the expectation and variance of J(D(ω)) as our first result. From now, we consider a twice shape differentiable function J.

Proposition 3.1 The expectation and variance of J(D(ω)) admit the asymptotic expansions

E(J(D)) = J(D 0 ) + εE(L(X)) + ε 2 2 E(B(X, X)) + O(ε 2 ) (8)
and

var(J(D)) = ε 2 E [L(X) -E(L(X))] 2 + ε 3 E ([L(X) -E(L(X))] [B(X, X) -E(B(X, X))]) + O(ε 3 ). ( 9
)
Proof of Proposition 3.1.

It suffices to integrate the pointwise Taylor expansion [START_REF] Dambrine | On the shape sensitivity of the first Dirichlet eigenvalue for two-phase problems[END_REF] over the space of probability thanks to the uniform estimate of the remainder. One thus immediately obtains for the expectation

E(J(D)) = J(D 0 ) + εE(L(X)) + ε 2 2 E(B(X, X)) + E(R 2 (εX)) = J(D 0 ) + εE(L(X)) + ε 2 2 E(B(X, X)) + O(ε 2 ).
Likewise, for the variance, one has

var(J(D)) = E [J(D) -E(J(D))] 2 = E ε(L(X) -E(L(X))) + ε 2 2 (B(X, X) -E(B(X, X))) + O(ε 2 ) 2 = ε 2 E [L(X) -E(L(X))] 2 + ε 3 E ([L(X) -E(L(X))] [B(X, X) -E(B(X, X))]) + O(ε 3 ). Remark 3.2 Let us notice that it follows E(L(X)) = L(E(X)) = 0 if X is centered, i.e., if E(X) = 0.
In the same spirit as above, we can also compute asymptotic expansions of higher moments. For example, the centered, normalized moments (as skewness and kurtosis) could be obtained in the following way.

Proposition 3.3 For all k ≥ 2, the centered normalized moment of J(D(ω)) admits the asymptotic expansion

M k (J(D)) := E [J(D) -E(J(D))] k var(J(D)) k = a k + b k ε + O(ε), (10)
where the deterministic coefficients a k and b k are

a k = γ k E [L(X) -E(L(X))] k E [L(X) -E(L(X))] 2 , (11)
b k = γ k k 2 E [L(X) -E(L(X))] k-1 [B(X, X) -E(B(X, X))] E [L(X) -E(L(X))] 2 -E ([L(X) -E(L(X))][B(X, X) -E(B(X, X))]) E [L(X) -E(L(X))] k , (12)
and the normalization constant γ k is

γ k = E [L(X) -E(L(X))] 2 -1-k/2 . ( 13
)
Proof of Proposition 3.3.

The second order Taylor expansion of the real function

x → x -k/2 around a real number x is 1 √ x + h k = 1 √ x k - hk 2 √ x k+2 + O(h 2 ).
Inserting the expansion of the variance (9), we find

1 var(J(D)) k = 1 ε k E [L(X) -E(L(X))] 2 -k/2 - εk 2 E [L(X) -E(L(X)] 2 -1-k/2 E ([L(X) -E(L(X))][B(X, X) -E(B(X, X))]) + O(ε 2 ) .
On the other hand, using the expansion of the expectation (8), we have

E [J(D) -E(J(D))] k = ε k E [L(X) -E(L(X))] k ) + εk 2 E [L(X) -E(L(X))] k-1 [B(X, X) -E(B(X, X))] + O(ε 2 ) .
Then, a tedious but elementary calculus (the product of two asymptotic expansions) leads to the desired result.

4 Efficient computation of the Taylor coefficients -the general case

Direct computation

For the sake of simplicity, we shall assume that the boundary perturbation field X is centered which induces E(L(X)) = 0. We then can compute as in [START_REF] Harbrecht | On output functionals of boundary value problems on stochastic domains[END_REF]]

E [L(X) -E(L(X))] 2 = E L(X) 2 = E ∂D ℓ(x)X(x, ω)dσ(x) 2 = E ∂D ∂D
X(x, ω)X(y, ω)ℓ(x)ℓ(y)dσ(x)dσ(y)

=

∂D ∂D E[X(x, .)X(y, .)]ℓ(x)ℓ(y)dσ(x)dσ(y).

The crucial operator to be studied here is the (two-point) autocovariance function of X defined as Cov X (x, y) = E(X(x, .)X(y, .)).

In fact, the knowledge of the law of X is not needed for the computation we have in mind, we only need the autocovariance function. Notice that various laws for X can provide the same function Cov X . We can proceed in two different ways to compute the expectation of the shape Hessian. The first approach has already been pointed out in e.g. [START_REF] Harbrecht | Sparse second moment analysis for elliptic problems in stochastic domains[END_REF]. We shall explain the proceeding in case of the Dirichlet energy. In that case, according to (3b), we have

E(B(X, X)) = E ∂D ∂ n u + 1 2 H(∂ n u) 2 X 2 dσ + E -∂ n u X, Λ(-∂ n u X) H 1/2 ×H -1/2 .
For the first term, one has easily by Fubini's theorem that

E ∂D ∂ n u + 1 2 H(∂ n u) 2 X 2 dσ = ∂D E(X 2 ) ∂ n u + 1 2 H(∂ n u) 2 dσ = ∂D V(X) ∂ n u + 1 2 H(∂ n u) 2 dσ
since X is centered. To compute the second term, we shall define the following tensor product type boundary value problem:

(id ⊗ (-∆))V = 0 in ∂D × D, V = (id ⊗ (-∂ n u)) Cov X on ∂D × ∂D. (14)
Here, -∂ n u has to be understood in the sense of a multiplication operator H 1/2 (∂D) → H 1/2 (∂D). Note that the differential operator which underlies this boundary value problem is id ⊗ Λ. Due to its linearity and Fubini's theorem, we get thus for the second term

E -∂ n u X, Λ(-∂ n u X) H 1/2 ×H -1/2 = ∂D E (-∂ n u(x) X(x, •)(Λ(-∂ n u X))(y, •)) x=y dσ(x) = ∂D -∂ n u(x) V (x, y) x=y dσ(x).
Boundary value problems like that in ( 14) can be solved in essentially linear complexity (with respect to the number of unknowns to discretize the state equation on the domain D) if a sparse tensor product discretization is employed as proposed in e.g. [START_REF] Harbrecht | Second moment analysis for Robin boundary value problems on random domains[END_REF][START_REF] Harbrecht | Combination technique based k-th moment analysis of elliptic problems with random diffusion[END_REF][START_REF] Harbrecht | Sparse second moment analysis for elliptic problems in stochastic domains[END_REF]. However, the implementation of this approach would be highly intrusive.

Toward a low-rank approximation

The second way, which is much simpler to implement, consists in computing an expansion of the autocovariance function of X of the form

Cov X = k κ k ⊗ κ k . (15)
For example, such a representation can be achieved by a spectral decomposition as an application of Mercer's theorem which ensures the representation of Cov X in the form

Cov X (x, y) = k λ k e k (x)e k (y). (16)
Another way to obtain such a decomposition is an (possibly infinite) Cholesky decomposition of autocovariance function. With the expansion (15) at hand, we will obtain in Section 5 the following expansion of the expectation of the shape Hessian:

E B[J](D).(X, X) = k B[J](D).(κ k , κ k ).
We therefore need only to be able to evaluate the shape Hessian in certain directions {κ i }.

Numerical realization

In general, the expansion (15) will be infinite and has to be appropriately truncated for numerical computations. Consider a suitable ansatz space V n = span{ϕ i : i = 1, 2, . . . , n} ⊂ C 2,1 (∂D). Then, we aim at a low-rank approximation

Cov X (x, y) ≈ m k=1 n i=1 ℓ k,i ϕ i (x) n j=1 ℓ k,j ϕ j (y) ∈ V n ⊗ V n (17)
to the two-point autocovariance Cov X ∈ L 2 (∂D × ∂D) in the tensor product space V n ⊗ V n . The unknown coefficient vectors in ( 17) can be computed as follows. Define the discrete autocovariance matrix C = ∂D ∂D Cov X (x, y)ϕ i (x)ϕ j (y)dσ(x)dσ(y) i,j ∈ R n×n and the mass matrix

G = ∂D ϕ i (x)ϕ j (x)dσ(x) i,j ∈ R n×n
with respect to the ansatz space V n . Then, [START_REF] Harbrecht | Second moment analysis for Robin boundary value problems on random domains[END_REF] amounts to an approximation

C ≈ C m = m k=1 ℓ k ℓ T k with ℓ k = [ℓ k,i ] i = G -1 ℓ k such that the truncation error C -C m is rigorously controllable.
The best low-rank approximation in L 2 (∂D × ∂D) is given by truncating the spectral decomposition [START_REF] Henrot | Variation et optimisation de formes. Une analyse géométrique[END_REF]. The computation requires the knowledge of the eigenpairs (ϕ i , λ i ) of the integral operator (Cu)(x) := ∂D Cov X (x, y)u(y)dσ(y), x ∈ ∂D which is a very demanding task. In particular, the decay of the eigenvalues {λ k } and thus the rank m depend heavily on the smoothness of the autocovariance function Cov X . Related decay rates have been proven in [START_REF] Schwab | Karhunen-Loéve approximation of random fields by generalized fast multipole methods[END_REF].

We suggest using the pivoted Cholesky decomposition to compute a low-rank approximation of Cov X as proposed in [START_REF] Harbrecht | On the low-rank approximation by the pivoted Cholesky decomposition[END_REF]. It is a purely algebraic approach which is quite simple to implement, see e.g. [START_REF] Harbrecht | On the low-rank approximation by the pivoted Cholesky decomposition[END_REF]. It produces a low-rank approximation to C for any given precision ε > 0 where the approximation error is rigorously controlled in the trace norm. A rank-m approximation is computed in O(m 2 n) operations. Exponential convergence rates in m are proven under the assumption that the eigenvalues of C exhibit a sufficiently fast exponential decay, see [START_REF] Harbrecht | On the low-rank approximation by the pivoted Cholesky decomposition[END_REF]. Nevertheless, numerical experiments show that, in general, the pivoted Cholesky decomposition converges optimally in the sense that the rank m is uniformly bounded with respect to the truncation error ε by the number of terms required for the spectral decomposition of Cov X to get the same error ε.

Computing the coefficients for a low-rank approximation of X

Let us apply Proposition 3.1 which provides the Taylor expansion of shape functionals. To that end, let (f i ) i∈N be a sequence of smooth (at least C 2,1 ) real valued functions which are defined on ∂D, possibly with compact support on ∂D if needed, and such that the series i∈N f i C 2,1 is convergent. We consider the domains defined by their boundary

∂D ε (ω) = {x + εX(x, ω)n(x); X ∈ ∂D}, (18)
with

X(x, ω) = ∞ i=1 α i (ω)f i (x) (19)
where the α i are either uncorrelated or independent and identically distributed following the same centered distribution L of a random variable α with finite second order moment. Note that such a model can also be provided by the Karhunen-Loève expansion, even though the f i , given by Mercer's theorem, are only in L 2 (∂D) and the additionally required regularity has to be checked, see e.g. [START_REF] Schwab | Karhunen-Loéve approximation of random fields by generalized fast multipole methods[END_REF].

Proposition 5.1 When the random graph has the form (18) with uncorrelated and centered random coefficients (α i ), one has

E(J(D)) = J(D 0 ) + ε 2 2 E(α 2 ) ∞ i=1 B(f i , f i ) + O(ε 2 ). ( 20
)
If the random coefficients are independent and if the distribution α has finite moments up to the order 3, then

var(J(D)) = ε 2 E(α 2 ) ∞ i=1 (L(f i )) 2 + ε 3 E(α 3 ) ∞ i=1 L(f i )B(f i , f i ) + O(ε 3 ). ( 21
)
Moreover, if the random variable α is symmetric, it even holds

var(J(D)) = ε 2 E(α 2 ) ∞ i=1 (L(f i )) 2 + O(ε 3 ). ( 22
)
Proof of Proposition 5.1.

We first notice that, for any integer k ≥ 1, the series

i∈N f i k C 2,1 is convergent.
1. Expansion of the expectation. By linearity, it comes directly

E(L(X)) = L (E(X)) = 0.
For any integer n, we introduce the truncated series X n defined as

X n (x, ω) = n i=1 α i (ω)f i (x),
so that

E(B(X n , X n )) = E B n i=1 α i f i , n i=1 α i f i = n i,j=1 E(α i α j )B(f i , f j ) = n i=1 E(α 2 )B(f i , f i )
since B is bilinear and the (α i ) are uncorrelated and centered. Since B is continuous, we get

|B(f i , f i )| ≤ B f i 2 C 2
and the series ∞ i=1 B(f i , f i ) converges, so that by continuity and Lebesgue's dominated convergence theorem

E(B(X, X)) = ∞ i=1 E(α 2 )B(f i , f i ).
2. Expansion of the variance. We proceed as before for the second order term:

E([L(X n ) -E(L(X n ))] 2 ) = E L n i=1 α i f i 2 = n i,j=1 E(α i α j )L(f i )L(f j ) = n i=1 E(α 2 )L(f i ) 2 .
We conclude by using the continuity of the map L that

E([L(X) -E(L(X))] 2 ) = ∞ i=1 E(α 2 )L(f i ) 2 .
For the third order term, one thus has

E ([L(X n ) -E(L(X n ))] [B(X n , X n ) -E(B(X n , X n )))]) = E L n i=1 α i f i B n i=1 α i f i , n i=1 α i f i = n i,j,k=1 E(α i α j α k ) L(f i ) B(f j , f k ) = n i=1 E(α 3 ) L(f i ) B(f i , f i )
where we used that the distribution of the random variable α is centered.

Remark 5.2 For all k ≥ 2, if the distribution α has finite moments up to the order k + 1, then the normalized, centered moment of J(D(ω)) admits the asymptotic expansion

M k (J(D)) = a k + b k ε + O(ε),
where the deterministic coefficients a k and b k could be obtained by a similar tedious but easy calculus. The crucial point is the computation of E(Π k j=1 α i j ) for i j ∈ {1, • • • , N } with the help of the multinomial formula. Thanks to the independence of the α i and the fact that they follow a centered distribution, the previous expression cancels as soon as there is an integer i ∈ {1, • • • , N } that is reached only once, i.e., if the equation i = i j for 1 ≤ j ≤ k has only a single solution.

For instance, for k = 3, it holds

M 3 (J(D)) := a 3 + b 3 ε + O(ε)
with

a 3 = E(α 3 ) γ 3/2 ∞ i=1 L(f i ) 3 , b 3 = 3 2γ 3/2   E(α 2 ) 2   i =j L(f i )L(f j )B(f i , f j ) - ∞ i=1 L(f i ) 2 B(f i , f i )   + E(α 4 ) ∞ i=1 L(f i ) 2 B(f i , f i )   - 3 2γ 5/2 E(α 3 ) 2 ∞ i=1 L(f i )B(f i , f i ) ∞ i=1 L(f i ) 3 , where γ = E(α 2 ) ∞ i=1 L(f i ) 2 .
With the help of the pivoted Cholesky decomposition, we can compute a low-rank approximation of the random field X which corresponds to a finite sum of the type

X(x, ω) = N i=1 α i (ω)f i (x)
where the α i are either uncorrelated or independent and identically distributed following the same centered distribution L of a random variable α. In that case, the Taylor coefficients can very easily be computed in accordance with Proposition 5.1: the computational cost reduces to the solution of N + 1 boundary values problems defined on the same domain so that a single mesh can be used. This amount is doubled to 2N + 2 resolutions in case of general shape functionals where also an adjoint state has to be computed to efficiently evaluate the shape derivatives (see [START_REF] Henrot | Variation et optimisation de formes. Une analyse géométrique[END_REF] for an example of the adjoint state in shape optimization).

Finally, using the Bienaymé-Chebychev inequality, we obtain intervals in which the shape functional takes its values with a fixed probability. For example, fixing a desired probability p, we get the bounds

P |J(D) -E(J(D))| ≤ var(J(D)) 1 -p ≥ p.
Thanks to the asymptotic expansions ( 20)-( 21), we obtain then

P J - ε,p (D) ≤ J(D) ≤ J + ε,p (D) ≥ p (23)
where we have set

J ± ε,p (D) = J(D 0 ) ± ε √ 1 -p E(α 2 ) N i=1 (L(f i )) 2 + ε 2 2         E(α 2 ) N i=1 B(f i , f i ) ± E(α 3 ) √ 1 -p N i=1 L(f i )B(f i , f i ) E(α 2 ) N i=1 (L(f i )) 2         + O(ε 2 ).

Numerical illustrations

We now present numerical validations of the proposed asymptotic expansions. We will consider the two-dimensional case for the Dirichlet energy and the first eigenvalue of the Dirichlet-Laplace operator. Numerical resolution of the boundary values problems is made either with the finite element method or with the boundary element method. We will consider both, the uniform and the Beta distribution.

During the Monte-Carlo simulations, the random number generator has been reinitialized after at most 10 000 samples.

The Dirichlet energy around disk in dimension two

Random perturbations defined by a given spectral representation

We consider the unit disk in the plane and random perturbations of the type

D ε (ω) = (r, θ) : 0 ≤ r < 1 + εf (θ, ω) and f (θ) = N i=1 α i (ω)f i (θ) (24)
where the α i are i.i.d. In the following computations, we have taken N = 11 and the coefficient functions f i are the first normalized functions in the Fourier basis, that is Let us present some realizations of such domains with various values of ε and various distributions. In Figure 1, we present some realizations of such random domains for α following the uniform distribution on

f 1 = 1, f 2 (θ) = cos(θ), f 3 (θ) = sin(θ), f 4 (θ) = cos(2θ)/4, f 5 (θ) =
[-1/2, 1/2]. Note that the range of α is [-1/2, 1/2].
In Figure 2, α follows the centered and normalized Beta distribution of parameter [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Chernov | First order k-th moment finite element analysis of nonlinear operator equations with stochastic data[END_REF] in order to use a non-symmetric distribution. Note that the range of this distribution is larger (it can take values greater than 4). Hence, for the same value of the parameter ε, the perturbations can be wider for the Beta distribution than for the uniform distribution. Therefore, we will consider a smaller range of values for the parameter ε. [START_REF] Pironneau | Optimal Shape Design for Elliptic Systems[END_REF] with the Beta distribution.

Monte-Carlo simulations

We shall proceed with the Monte-Carlo method. For each value of ε in the abscissa of the graphs, we perform 4 000 simulations using P 1 finite elements over around 2000 triangles, each of which is made with the FREEFEM++ code. The convergence is illustrated in Figure 3. It seems that a visually correct asymptotic regime is reached for some thousands of simulations.

Comparison with the asymptotic expansions

Finally, we compare with the asymptotic expansions obtained in Proposition 5.1. Let us emphasize that the coefficients appearing in Proposition 5.1 are computed numerically. The first step is to compute the state equation in the reference domain, then the shape gradient is obtained via the computation of an integral over the boundary. Finally, the computation of each term B(f i , f i ) requires the solution of a boundary value problem of the same type than the one solved for the state equation. Therefore, the computation of these coefficients is cheap since it requires only N + 1 resolutions of a boundary value problem. The results for both, the uniform distribution and the Beta distribution, are presented in Figure 4.

We can also use [START_REF] Novruzi | Newton's method in shape optimisation: a three-dimensional case[END_REF] to obtain bounds of intervals of values taken by the Dirichlet energy. In Figure 5, we have plot both, the graphs of the functions J ± ε,p and the empirical quantiles. Notice that the bounds obtained by the asymptotic expansions of the expectation and the variance are better than expected. Such a nice behavior may depend on the specific choice of the Fourier basis, since the shape gradient and the shape Hessian evaluated at the f i decrease fast as shown by the following computed values:

L(f 1 ) = -0.786069 B(f 1 , f 1 ) = -2.35621 L(f 2 ) = -1.74702.10 -6 B(f 2 , f 2 ) = -0.392049 L(f 4 ) = 4.71636.10 -6 B(f 4 , f 4 ) = 0.0246288 L(f 6 ) = 3.00416.10 -7 B(f 6 , f 6 ) = 0.0145671 L(f 8 ) = 4.24149.10 -8 B(f 8 , f 8 ) = 0.00767894 L(f 10 ) = -1.48682.10 -7
B(f 10 , f 10 ) = 0.00440313

Random perturbations known by their mean and covariance

In this second test, we change the type of random perturbations: instead of being decomposed with respect to the Fourier basis, X(ω) is now described by its zero mean and its covariance function. We shall consider models with standard covariance kernels, namely Mátern kernels of smoothness class ν. We still compare a Monte-Carlo simulation with the asymptotic expansions derived in this work for different values ε ≤ 0.2. For the Monte-Carlo simulation, the truncated Karhunen-Lóeve expansion is exploited:

X(ω, θ) = m i=0 α i (ω) λ i ϕ i (θ).
Herein, the same accuracy is used as for the pivoted Cholesky decomposition which yields nearly the same ranks m. Moreover, the random variables are i.i.d. and modeled as uniformly distributed in [-2/ √ 3, 2/ √ 3]. The numerical computation of the shape functional, shape gradient and shape Hessian is now performed by a fast boundary element method as outlined in [START_REF] Eppler | Second-order shape optimization using wavelet BEM[END_REF]. The Monte-Carlo method uses 10 000 samples. As seen in Figures 6,7, and 8, the asymptotic expansions fit well the behaviour of the expectation and variance of the stochastic shape functional under consideration for all covariance kernels under consideration.

The first eigenvalue of the Dirichlet-Laplace operator in a perforated ellipse

To consider another output function, we present simulations for the first eigenvalue of the Dirichlet-Laplace operator. Again, the related shape gradient and shape Hessian are obtained thanks to Lemma 2.2. Here, the reference domain D is the ellipse of semi-axes 4 and 3 where the disk of radius 1 which is centered in the center of the ellipse has been removed. Random perturbations are applied to the disk's boundary. As in the first example, we take

f (θ) = N i=1 α i (ω)f i (θ)
where the α i are i.i.d. In the following computations, we have taken N = 21 and the f i are the first normalized functions in the Fourier basis. We emphasize the local character of the approximation which is clearly seen in Figure 9.

Theorem 2 . 1 (

 21 Structure theorem of first and second order shape derivatives) Let k ≥ 1 be an integer and J a real valued shape function which is defined on O k , the open bounded domains of R d with a C k -boundary. Let us define the function J on C k,1 (R d , R d ) by J (θ) = J((I + θ)(D)).

  sin(2θ)/4, . . . , f 11 (θ) = sin(5θ)/25. the random domain -epsilon=1/4 (b) ε = 1/4.

Figure 1 :

 1 Figure1: Some realizations of domains according to[START_REF] Pironneau | Optimal Shape Design for Elliptic Systems[END_REF] with the uniform distribution.

8 (

 8 the random domain -Beta law -epsilon=1/the random domain -Beta law -epsilon=1/4 (b) ε = 1/4.

Figure 2 :

 2 Figure 2: Some realizations of domains according to (24) with the Beta distribution.

Figure 3 :

 3 Figure 3: Convergence of Monte-Carlo simulations for perturbed disks defined in (24).

Figure 4 :

 4 Figure 4: Comparison of the Monte-Carlo simulation with the asymptotic expansions for perturbed disks defined in (24).

  Specifically, we choose ν = ∞ which yields the Gaussian kernel k ∞ (x, y) = exp(-|x-y| 2 /2), we choose ν = 1/2 which yields the exponential kernel k 1/2 (x, y) = exp(-|x -y|), and we choose ν = 3

Figure 5 :

 5 Figure 5: Comparison between asymptotic expansions' bounds and empirical quantiles for perturbed disks defined in (24).

Figure 6 :

 6 Figure 6: Expectation and variance of the Dirichlet energy -Gaussian covariance kernel k ∞ .

Figure 7 :

 7 Figure 7: Expectation and variance of the Dirichlet energy -Mátern covariance kernel k 3/2 .

Figure 8 :

 8 Figure 8: Expectation and variance of the Dirichlet energy -Exponential covariance kernel k 1/2 .

Figure 9 :

 9 Figure 9: Comparison between asymptotic expansion and Monte-Carlo simulation for the eigenvalue problem -Uniform distribution.

Acknowledgements:

The first author thanks support from the ANR through the grants OPTIFORM and ARAMIS.