SMAI-JCM SMAI Journal of Computational Mathematics A reconstruction method for the inverse gravimetric problem

. We propose a reconstruction method to solve the inverse gravimetric problem with constant mass density. The method is based on the computation of the harmonic moments of the unknown domain. Convergence results are proved and numerical experiments are provided to illustrate the method and show its efficiency.

The gravitational potential generated by a uniform mass distribution in ω is given by the Newtonian potential where G(x) := −1/(2π) ln |x| is the fundamental solution of the Laplacian in R 2 and m is the Lebesgue measure.We can easily verify that the function U ω satisfies the Laplace equation and admits the following behavior at infinity (|ω| stands for the area of ω): 2) The inverse problem we are interested in is often referred to as the inverse gravimetric problem and reads as follows: how to reconstruct ω from the knowledge of ∇U ω on Γ?
Remark 1.1.Two comments are in order.First, it is natural to assume that the available measurement is ∇U ω (and not U ω ) on Γ since this quantity represents the gravitational field, the actual quantity measured in real experiments, typically in geophysics.Second, it is possible to recover U ω on Γ from the knowledge of ∇U ω .Indeed, integrating the identity (1.1) over B and applying Green's formula leads to the expression of the area of ω: Then, the expression of U ω in B + = R 2 \B follows by remarking that U ω = u+|ω|G in B + where the function u is the unique solution of the following exterior Neumann problem (see [29,Theorem 8.18] for the well-posedness):

Background on uniqueness and stability
The inverse problem stated above is known to be severely ill-posed as proved in [20,Section 3.4], since even uniqueness is not guaranteed in general.Nonetheless, uniqueness results can be established under either one of the following assumptions: (S) ω is star-shaped with respect to its center of gravity.
(C) ω is convex in one direction (say x 1 ), i.e the intersection of any straight line parallel to the x 1 -axis with ω is an interval.
It is worth mentioning that these uniqueness results have been improved by Gardiner and Sjödin [9,10] that proved that for solid domains, uniqueness is ensured if we know that one of the two domains is convex.
A stability result can also be found in the literature, at least in the specific class of star-shaped domains having the same center of gravity.This common center of gravity is assumed to be 0 for the rest of this subsection for the sake of simplicity.In polar coordinates, any such a domain ω can be described by means of a function ρ ω : [0, 2π[−→ R + defined by: ρ ω (θ) = max{ρ ⩾ 0 : (ρ cos θ, ρ sin θ) ∈ ω} for all θ ∈ [0, 2π[.
Using this notation we can state the following logarithmic-type stability estimate [20,Theorem 3.6.1]: A reconstruction method for the inverse gravimetric problem Theorem 1.3.Let M, ρ − , ρ + and α be four positive constants and let: Then, there exists κ > 0 such that for all star-shaped domains ω 1 and ω 2 satisfying ρ ω 1 , ρ ω 2 ∈ R we have

Overview of the method
The aim of this paper is to develop an algorithm that generates a sequence of domains (ω N ) N converging (in some sense) to the unknown domain ω, or at least whose gravitational fields fit the targeted one (i.e, is graviequivalent to ω).To achieve this goal, we only make use of the harmonic moments of ω, quantities that can easily be deduced from the measurements.Indeed, for every function v harmonic in B, Green's formula yields: where we recall that ∇U ω and U ω are both supposed to be known on Γ (see Remark 1.1).Thus, we are led to study a general shape-from-moments problem which can be stated as follows: how to reconstruct a domain ω from the knowledge of its harmonic moments?Throughout the paper, any complex number z = x 1 + ix 2 will be identified with the point (x 1 , x 2 ) of the plane.Let us denote by z ℓ (ℓ ⩾ 0) the harmonic monomials in the plane.Given an integer N ⩾ 1, the main idea of this work is to construct a sequence of domains (ω N ) N such that: for all ℓ = 0, . . ., 2N − 1.
To achieve this goal we proceed in two steps.For any positive integer N : Step 1 Construction of a quadrature formula associated with the data, i.e. computation of weights c 1 , . . ., c N ∈ C and nodes z 1 , . . ., z N ∈ C such that: The values of the c j and z j (j = 1, . . ., N ) will be obtained by means of the so-called Prony's method.
Step 2 Construction of a quadrature domain associated with the previous quadrature formula, i.e. determination of a domain ω N (if any) satisfying identity (1.5) but for any harmonic polynomial this time: Such a domain is called a "quadrature domain" and for given c k and z k (k = 1, . . ., N ), its existence is not guaranteed in general.We shall test two approaches to reconstruct this domain: the first one rests on purely complex analysis tools while the second one, called "partial balayage", is formulated as a convex minimization problem.
Notice that our method may be applied to any shape-from-moments problem, though we are able to obtain convergence results only for the gravimetric problem.
Our method is based on an original combination of classical notions.The relation between gravitational potential and quadrature domains is natural (see for example [12]).On the other hand, although Prony's method (and related approaches) is a classical tool in signal processing (particularly in exponential fitting), its use for the reconstruction of quadrature domains is not very common.

Related works
Inverse gravimetric problems may play an important role in future practical applications (primarily in geophysics), more especially since recent advances in gravitational field measurements [39].As far as theoretical questions are concerned, the book by Isakov [20] remains an essential reference.Let us also mention the monograph of Cherednichenko [6].Following Potthast in his survey paper [32], it is relevant to classify the methods for solving geometric inverse problems in two main categories: direct and iterative ones.
Our strategy belongs to the first category, and has similarities with those presented in [11,13].In these articles, the authors address the problem of reconstructing an unknown domain from its moments.Only polygonal shapes are considered in the former while all the moments (and not only the harmonic ones) are required for the reconstruction in the latter.
Moments methods were also used in [33] to recover the Fourier expansion of the boundaries' parametrizations of domains having the property (S).Finally, let us also mention [5], in which a ball fitting approach, evoking our ball representation of quadrature domains, is used.
In the class of iterative methods, a first glance shows the predominance of those based on levelset optimization [21,22,26,27,28,40].Before the spreading of this powerful general tool, Hettlich and Rundell [18] tackled the problem using shape derivatives coupled to a Newton-type method, while Zidarov implemented a least squares approximation, as described in his classic book [41].Finally, Kress and Rundell [25] proposed an original rewriting of this problem using boundary integral equations.

Outline of the paper
The paper is organized as follows: Step 1 and Step 2 are elaborated in Section 2 and 3 respectively.Next, we establish convergence results in Section 4 and finally we present some numerical experiments and give some details about their implementation in Section 5.

Construction of the quadrature formula
Recall that Step 1 of our method consists in seeking, for any positive integer N , weights c 1 , . . ., c N ∈ C and nodes z 1 , . . ., z N ∈ C such that: where the quantities in the right hand side are the harmonic moments of the domain ω defined by: 2) The non-linear system (2.1) turns out to be a so-called Prony's system.Prony's systems have been extensively studied, as they appear in several fields of theoretical and applied mathematics, such as Padé's approximants, signal processing or error correction codes.
A reconstruction method for the inverse gravimetric problem Remark 2.1.Observe that (2.1) can be rewritten as: for all ℓ = 0, . . ., 2N − 1, ( where the measure µ N is equal to the linear combination of Dirac masses N k=1 c k δ z k .Intuitively, identity (2.3) suggests that we are trying to approximate the measure 1 ω with point masses.It is thus natural to expect the weights c k to be positive for all k.However we do not add this physical constraint, for the resulting constrained system may have no solution.
We formulate the problem of solving (2.1) as follows: Let N be a positive integer and τ 0 , . . ., τ 2N −1 ∈ C be given.Determine for some positive integer p, z 1 , . . ., z p ∈ C pairwise distincts and c 1 , . . ., c p ∈ C \ {0} such that (2.4) Let us emphasize that in this statement, the τ ℓ are arbitrary complex numbers and do not need to be harmonic moments as in (2.2).

Existence and uniqueness
When p is allowed to be as large as 2N , Problem (2.4) admits always a solution, as proved for instance in the first theorem in [31, Section 2].However, this does not meet our requirement in (2.1), which is p ⩽ N .Moreover, the nodes z k of this solution lie on the unit circle, which is not satisfactory for our purpose.Indeed, we want to reconstruct the unknown domain ω (compactly contained in the unit disk by assumption), so we expect the nodes to be located in ω (this intuition will become obvious in Section 3).
The next results stated in this subsection might have been already proved, but we have not been able to find them in the literature.The proofs are given in Appendix A.
We begin with introducing the matrix: The statement of the existence result requires introducing the following polynomial: A. Gerber-Roth, A. Munnier, et al.
Theorem 2.4.Problem (2.4) admits a solution for which p = N if and only if the polynomial P N admits N simple roots.In this case, this solution is unique and the nodes z 1 , . . ., z N are the roots of P N .

Numerical solutions
A review of available methods to solve (2.1) can be found in [19].Roughly speaking (see, e.g., [38,Section 4]), one can distinguish direct nonlinear minimization methods (typically nonlinear least squares methods), recurrence-based methods (such as the original Prony's method), and subspace methods (such as Pisarenko's method, MUSIC, ESPRIT and matrix pencil methods).In this work, we apply the matrix pencil method, as described in [11].
Let us introduce another Hankel matrix of size N (obtained by applying a shift to the indices of H (N ) 0 ): First we denote by z 1 , . . ., z N the (complex) eigenvalues that verify the following generalized eigenvalue problem: , admits solutions c 1 , . . ., c N , one can verify that z 1 , . . ., z N and c 1 , . . ., c N solve (2.1) (see [11,Section 3]).Since the inverse gravimetric problem is ill-posed, it is not surprising that it leads to Prony's systems which are generally ill-conditioned as well (see [3]).Actually, the above numerical method involves solving two ill-conditioned problems: a generalized eigenvalue problem with Hankel matrices [35] and a Vandermonde system.To improve the overall condition number, we use in our implementation the regularization techniques detailed in [11,Section 4].It is worth mentioning also that numerical instabilities increase along with N .

An example of Prony's system without solution
In this section, we provide an example of domain whose moments lead to a Prony's system that does not have a solution (and hence for which the first step of our reconstruction algorithm fails).Consider the domain ω whose boundary is described in polar coordinates by ρ(θ) = 1 + 2a cos(θ) with 0 < a < 1/2 (see Figure 2.1).
For every N ⩾ 1, let the Hankel matrix H (N ) 0 and the polynomial P N be defined as in Section 2.1, where the measurements (τ ℓ ) ℓ⩾0 are given by (2.2).
(2) rank H A reconstruction method for the inverse gravimetric problem Proof.For all ℓ ⩾ 0, direct computations lead to: and since: we end up with: Consequently, det H which proves the first point.Concerning the second point, we have: .
Multiplying the k-th row by a −(k−1) and the j-th column by a −(j−1) , we obtain that this matrix has the same rank as , which has the same rank as , as it can be seen by replacing column C j by C j − C j−1 for all j ⩾ 2.
We are now in position to prove: Proposition 2.6.For every N ⩾ 2, Problem (2.4) has no solution for which p ⩽ N .
Proof.The result is proved by induction on the index N .For N = 2, H 0 is invertible, so it follows from Proposition 2.2 that every solution for which p ⩽ 2 would actually be such that p = 2.However, since P 2 has a double root, Theorem 2.4 asserts that Problem (2.4) has no such solution.Assume now that the claim of the Proposition holds true for some N ′ ⩾ 2. We first notice that the coefficient in front of the monomial z Consequently, P N ′ +1 cannot have N ′ + 1 simple roots and it follows that Problem (2.4) for N = N ′ + 1 has no solution for which p = N ′ + 1.If a solution with p < N ′ + 1 existed, it would be also a solution to Problem (2.4) for N = N ′ , yielding a contradiction with the assumption.
Pathological domains, like the example provided in this subsection, seem to be interesting mainly from a theoretical point of view.In practice, we have never encountered this type of situation.However, the cancellation of det H (N ) 0 for some values of N (at least in the floating-point precision sense) is common and arises typically when N is too large (see also Remark 4.8 for a related discussion).

Construction of the quadrature domain
We adress now the second step of our method summarized in Section 1.3.This step consits in constructing a domain associated with a given quadrature formula.More precisely, given N nodes z 1 , . . ., z N ∈ C and N weights c 1 , . . ., c N ∈ C (those computed in the first step of the reconstruction algorithm detailed in Section 2), we want to construct a domain ω N such that for all harmonic functions v on ω N : Obviously, such a domain will satisfy in particular: The domains satisfying quadrature identities in the genus of (3.1) are called quadrature domains.The simplest quadrature domain is provided by the disk D(a, R) of center a ∈ C and radius R > 0 for which, according to the mean value property: for every function v harmonic in D(a, R).Let us now collect some of the main properties of quadrature domains that will be needed in the sequel (for a complete introduction to quadrature domains, we refer the interested reader to [16] and references therein).
A reconstruction method for the inverse gravimetric problem

Definitions and properties
Definition 3.1.A planar domain Ω is called a quadrature domain associated to a complex Radon measure µ when: for all complex analytic integrable function v in Ω.The set of all domains satisfying this property is denoted by Q(µ, AL 1 ).Replacing "complex analytic" by "harmonic" in this definition leads to the definition of harmonic quadrature domains which form the set Q(µ, HL 1 ).Eventually, a domain Ω satisfying: for all subharmonic integrable function v, is called a subharmonic quadrature domain.The set of all the subharmonic quadrature domains is denoted by Q(µ, SL 1 ).
Remark 3.2.The following inclusions hold: (the latter is obvious and for the former, notice that when v is a harmonic function, v and −v are both subharmonic).
Returning to our main concern, that is to construct domains satisfying (3.1), we shall focus in the sequel on the particular case where µ is an atomic measure.More precisely, we assume from now on that: where N is a positive integer, z 1 , . . ., z N are pairwise distinct complex numbers (the nodes), and c 1 , . . ., c N are complex numbers (the weights).
Remark 3.3.Depending on the type of quadrature domain considered, the weights c k must satisfy additional constraints: (1) If Ω ∈ Q(µ N , HL 1 ), then all the weights c k are real.Indeed, taking the imaginary part in (3.3) leads to: for all harmonic function v, and thus Im(c k ) = 0.
As already mentioned, the disks are quadrature domains and it can be shown that they are the only ones for which N = 1 [34,Example 1.1].Rather counter-intuitively, there are many quadrature domains.Actually, every domain whose boundary consists in non-intersecting C ∞ Jordan curves is arbitrarily close to a quadrature domain (see [4,Section 2]).
Given a set of nodes and weights, existence of a subharmonic quadrature domain is asserted by Gustafsson in [12, Theorem 2.4 (vi)], providing that the weights are positive (recall that every subharmonic quadrature domain is also a harmonic quadrature domain by (3.5)).The construction of this domain from µ N is strongly related to partial balayage of measures [14] and free boundary problems [15].Things are more complicated when it comes to the question of uniqueness.Focusing again on subharmonic quadrature domains, uniqueness is proved in the already mentioned paper [12, Theorem 2.2].In view of (3.4), uniqueness means up to a set of null Lebesgue measure and in [12, Theorem 2.2], the author claims that this domain can be chosen open.This uniqueness property is lost when we move to harmonic quadrature domains.Indeed, even in the simplest case of an atomic measure like µ N (with positive weights), uniqueness does not hold in general.An example is provided in the paper [1], to which we refer for general questions related to uniqueness.
We shall now present two numerical methods allowing the construction of a quadrature domain from the knowledge of a given set of nodes and weights: a conformal mapping method and a convex minimization method.It is worth noticing that these methods won't always succeed if we do not assume the c k to be positive.In particular, characterizing nodes and weights for which there exists a quadrature domain is an open problem for general complex weights c k .

Conformal mapping
Let N be a positive integer, and z 1 , . . ., z N be complex numbers pairwise distincts and located in the open unit disk B. Consider in addition complex weights c 1 , . . ., c N .The method used in this subsection to recover a quadrature domain from these numbers is taken from the recent article [1].It allows for complex weights but is able to yield only quadrature domains ω N that are connected and "solid" (i.e. an open set such that ω e N = R 2 \ ω N is connected and ∂ω e N = ∂ω N ).More precisely, let us consider µ N as introduced in (3.6) and assume that at least one solid domain ω N containing the origin belongs to Q(µ N , AL 1 ).Then it can be obtained as the image of the unit disk B by a conformal map φ N taking the form: where x 1 , . . ., x N and λ 1 , . . ., λ N are complex parameters (with |λ k | < 1 for evey k).These parameters are related to the nodes z 1 , . . ., z N and the weights c 1 , . . ., c N through the identities: for all k = 1, . . ., N. ( To solve this system we use the pre-implemented Matlab function fsolve, which is based on an iterative method.Since the nodes and weights are supposed to be the harmonic moments of the domain ω N , it follows that N k=1 c k = |ω N | is real and positive.We initialize the method with the parameters: where r k := |c k |/π.With these values for the parameters, the image of the unit disk by the conformal mapping φ N is a disk of area |ω N | i.e. with the same first harmonic moment as ω N .

Convex minimization
Consider now in addition that the weights c 1 , . . ., c N are real numbers and define the measure µ N as in (3.6).We aim to construct a solid element of Q(µ N , HL 1 ).For this purpose we introduce a balayage A reconstruction method for the inverse gravimetric problem operator that will provide, under some conditions, the characteristic function of the quadrature domain ω N .This balayage operator can be evaluated by solving a convex minimization problem.This will be done using Finite Element Method (FEM).
In the following, we denote by M c the set of real Radon measures with compact support and by L ∞ c the set of compactly supported functions in L ∞ (R N ).As already mentioned, we will identify two domains when they differ only by a set of measure zero.
It will be more convenient to deal with the following regularized version of the measure µ N : for given α 1 , . . ., α N > 0. The next result shows that µ N and µ N yield the same quadrature domains.
Proposition 3.4.We have: Proof.By the mean value property, for every harmonic function v we can write and the result follows.
We define now the balayage operator F : Definition 3.5.For every µ ∈ M c , we set and we define: (1) Some properties of the balayage operator F are collected in the next proposition (a more complete analysis can be found in [12,14]).Proposition 3.6.Let µ ∈ M c , then A. Gerber-Roth, A. Munnier, et al.
(3).From, [14, Example 2.2], we obtain Ω = Ω(µ) up to a null set, and The cornerstone of our method consists in noticing that when F ( µ N ) is a characteristic function of a domain ω N : then we can deduce that: If in addition ω N is solid, it is a solution to our problem.Let us emphasize that it is crucial to work with the regularized measure µ N because results from [12,14] do not apply to the singular measure µ N .
A reconstruction method for the inverse gravimetric problem

Basic elements from potential theory
For reader's convenience, we recall in this subsection some basic properties from potential theory used in the sequel.We shall focus on the single layer potential and we refer to the monograph [29] for a thorough presentation of potential theory and boundary integral operators.
Let D be a bounded Lipschitz domain whose boundary C is a Jordan curve.Let n be the unit normal directed towards the exterior of D. We denote by γ − and γ + the one-sided trace operators respectively from the interior and the exterior of D. We use a similar convention for the normal derivatives ∂ ± n .Definition 4.1.The single layer potential is the weakly singular integral operator defined for any q ∈ L 2 (C ) by: and extended by density as a bounded operator S C : According to [29, p. 261], the single layer potential admits the following asymptotic expansion as |x| goes to infinity: where ⟨ • , • ⟩ C stands for the duality pairing on C ) extending the L 2 scalar product.For every q ∈ H −1/2 (C ), the outer trace γ + S C q and inner trace γ − S C q coincide on C .We denote by S C q this common value and we recall that the operator S C : Regarding the inner and outer normal traces, the single layer potential satisfies the so-called "jump relation", namely: Let us recall the following result involving the operator S C , referred to as Theorem 8.16 in [29].
Let us emphasize that the assumption on the diameter of D is not restrictive (otherwise, it suffices to rescale the problem) and simply ensures that the capacity of C is less than 1 (see, for instance, [37, p. 143] and references therein).We shall also require the following uniqueness result (see [29,Lemma 8.14]):

Main results
The main steps of our algorithm are summarized in Algorithm 4.1.
At the end of this process we then obtain a domain ω N such that We denote by U ω N the Newtonian potential G * 1 ω N and we claim: and next, by continuity of the trace operator, there exists a constant C 1 > 0 such that: Now, since x → G(x, y) is harmonic when y ∈ Γ r , we can write (recall that c k,N ∈ R from Remark 3.3) On the other hand, according to Theorem 4.3, there exists a unique (q ω , c ω ) and ⟨q ω , 1⟩ Γr = |ω|.Comparing (4.1) and (1.2) we deduce that: ) and thanks to [29, Theorem 8.10] (uniqueness result for the exterior Dirichlet problem), we conclude that c ω = 0 and for all |x| > r.
It follows that for every harmonic function v in H 1 (B): Integrating by parts over B r the first term in the integral, and taking into account the jump relation (4.2), we end up with: The very same computations can be carried out replacing ω with ω N and we get as well: where ℓ for every positive integer ℓ and every (x 1 , x 2 ) ∈ R 2 , we obtain, combining (1.4) and (4.3): for all ℓ = 0, . . ., 2N − 1. Returning to (4.4) and invoking the compact embedding of H 3/2 (Γ r ) into H 1/2 (Γ r ), we can assume that, up to a subsequence extraction, the traces of the functions U ω N strongly converge in H 1 2 (Γ r ) as N goes to +∞.Since S Γr has a bounded inverse (as asserted by Theorem 4.2 above), the sequence (q ω N ) N converges in H − 1 2 (Γ r ) to some density q and for every ℓ ⩾ 0 we have: Letting N go to +∞, we get: The subspace of trigonometric polynomials being a dense subset in H 1 2 (Γ r ) (see [24,Theorem 8.2]), the identity (4.7) leads to q = q ω .Notice now that the function G( • , • ) and all its derivatives are bounded on the compact set O × O r , which entails the existence of a constant C 2 > 0 such that: A. Gerber-Roth, A. Munnier, et al.
By continuity of the trace operator, it follows that there exists a constant C 3 > 0 such that: . The conclusion follows by letting N −→ +∞ and noticing, by a classical argument, that the uniqueness of the limit ensures that the convergence holds for the whole sequence (U ω N ) N and not only for a subsequence.As expected, the previous theorem gives a convergence result only for the gravitational fields.However, adding the hypotheses of Theorem 1.3 we can prove the following result regarding the domains.
Corollary 4.6.Under the hypotheses of Theorem 4.4 assume in addition that ω and ω N both satisfy (S) from Section 1.2 and that their polar representations ρ ω and ρ ω N belong to the set R for all N ⩾ 1, where R is defined in Theorem 1.3.Then, (ω N ) N converges to ω in the sense of Proof.Note that, thanks to (4.3), ω and all the ω N have the same center of gravity and then apply Theorems 1.3 and 4.4.
In the special case where the unknown domain is a subharmonic quadrature domain the above result can be strongly improved, as an exact reconstruction of ω can be achieved in a finite number of steps.Theorem 4.7.Assume that ω is a subharmonic quadrature domain associated to a finite sum of point masses.Let (ω N ) N be the sequence of subharmonic quadrature domains constructed with the algorithm 4.1, using partial balayage.Then there exists a positive integer M such that ω M = ω.According to [7,Lemma 2], the value of M can be determined by noticing that: rank where H (n) 0 is defined in (2.5a).Remark 4.8.Theoretically, increasing the number N of nodes in our algorithm should increase the accuracy of the reconstruction.In practice, it also increases the numerical errors (mainly in the resolution of Prony's system).Hence taking N as large as possible does not always result in a successful strategy.Inspired by (4.8), we could take for N the number of "significant" singular values of the matrix H (n) 0 , with n large enough.
A reconstruction method for the inverse gravimetric problem

Numerical experiments
We present in this section some numerical experiments.For every simulation, we choose a domain ω to be reconstructed and compute its gravitational potential U ω and gravitational field ∇U ω to which we apply a white gaussian noise.These data are then taken as inputs in Algorithm 4.1 for different values of N .
The harmonic moments are computed from these boundary data using the Fast Fourier Transform (FFT).In order to reduce the impact of noise on the measurements, we use a large number of quadrature points to evaluate the moments.
From the harmonic moments, we construct the Hankel matrices (2.5) and solve the corresponding Prony's system which provides the (complex) nodes z 1,N , ..., z N,N and (complex) weights c 1,N , ..., c N,N .
Finally we compute a quadrature domain using one of the two methods previously introduced (in Section 3.2 and Section 3.3).The results are presented in the following two subsections.

Conformal mapping
As this method presents convergence issues, we restrict our experiments to the cases where ω satisfies Property (S) and we use noiseless data.The boundary of the target domain is parameterized in polar coordinates by a function ρ : ] − π, π] −→ R + taking the form of a Fourier series with a finite number of modes: We test the conformal mapping reconstruction method with the domains depicted on ) ⩽ 0, we color the corresponding disk in red).On the second row, we plot the image of the unit disk by the computed conformal mapping φ N , as introduced in (3.7), for several values of N .We can notice that the quality of the reconstruction improves with N and that we obtain a fairly accurate approximation of the target shape.However, let us emphasize that the nonlinear solver is quite sensitive to the shape of the target domain and to the initial guess.For example, on Figure 5.2d and Figure 5.3e, the algorithm failed to converge.It is not clear whether this failure is due to the lack of existence of solutions to System (3.8) or to a poor choice of initial guess.Nevertheless, one may notice that a partial balayage applied to the measure obtained by replacing the complex weights by their real parts, would provide a fairly good approximation of the target shape.This can be seen by considering the disks representation (see next subsection).Besides its sensitivity to the initial guess, this algorithm is limited to the cases of solid domains.Its main advantage is that it allows for non-real (and non-positive) weights.However, in our tests, the imaginary parts of the weights were systematically close to zero.For all these reasons, in most of the cases, the partial balayage method should be preferred.

Convex minimization
We focus hereafter on the partial balayage method (described in Section 3.3).Since we need the computed weights c k,N to be real, we replace them by their real part if necessary.In our tests, though, this choice turned out to be of little importance because we observed that the imaginary part was always very small with respect to the real part.We recall that the problem to be solved is summarized in (3.16).We follow the approach described in [23] to solve the variational inequality involved in this convex minimization problem and we implemented this algorithm using FreeFEM++ [17].We use P 2 finite elements in order to compute the laplacian of the minimizer, and hence evaluate an approximation of the operator F .In the tests of this subsection, we use the regularization parameter α = 0.99.
We begin with a domain having the property (S) (see Figure 5.4).On each picture is displayed the targeted domain (in red), the disks (B(z k , r k )) k associated with the solution of Prony's system (the centers are materialized by small markers), and a finite element approximation of F (μ α N ), for different values of N .These tests are performed without noise and allow a fairly accurate reconstruction of the domain.One may notice that some centers lay outside the domain (for N = 15), but they correspond to nodes for which the weight is zero.We wish now to investigate the influence of noise through some examples: simple ones (the previous star-shaped domain in Figure 5.5, a domain convex in one direction in Figure 5.6 and a set of disjoint disks in Figure 5.7) and more sophisticated ones (a rabbit-shaped domain in Figure 5.8 and Figure 5.9, and a plane-shaped one in Figure 5.10 and Figure 5.11).For these two last examples, we represent only the tests with a reasonable ratio between the truncation parameter N and the noise.For all these cases, even in the presence of noise, we observe that the algorithm generates very few positive mass outside the targeted domain.Although the noise has a deteriorating effect on the details, the raw reconstruction of the domain contour seems to be relatively robust.In Figure 5.7, note the good quality of the reconstruction on the last image, which is unusual at this noise level.
One limitation of our method is the condition number of the first step, as invoked in Section 2.3.Indeed, for some values of N the computations may conduct to totally wrong reconstructions.We illustrate this fact in Figure 5.12, for which we deal with the same unknown domain as in Figure 5.6.We remark that in Figure 5.12a, all the weights are almost equal to zero.Of course, searching for an associated quadrature domain do not yield a satisfactory approximation of ω.However, such situations can be avoided by choosing another value for N , as visible in Figure 5.12b.

Conclusion
In this work, we propose an algorithm to solve the planar gravitational inverse problem.The approach is based on an original combination of Prony's method and generation of quadrature domains.Under some technical hypotheses, we are able to prove that the sequence of reconstructed domains converges A. Gerber-Roth, A. Munnier, et al.
These equations can be reformulated as the following linear system: Since H (N ) 0 is invertible, the vector q is uniquely determined and hence this is also true for Q and for its roots (z 1 , . . ., z N ).Uniqueness for the weights follows because they are solution of a linear system with an invertible Vandermonde matrix.
Proof of Corollary 2.3.Assume that there exists a solution with p = N to Problem (2.4).In that case, we deduce from (A.1) that H (defined in (2.5a)) and solve the linear system: .
(2) Find the N roots z 1 , . . ., z N of the polynomial: β j z j .
Proof of Theorem 2.4.Assume that Problem (2.4) admits a solution with p = N and denote by z 1 , . . ., z N the nodes.We introduce β 0 , β 1 , . . ., β N −1 such that: In the other hand, let for all z ∈ C.

A reconstruction method for the inverse gravimetric problem
One can verify that = 0, for all k = 1, . . ., N and thus, det(M (z k )) = 0 = P N (z k ), which proves that z k are the roots of P N (remark that the degree of P N is N ).Conversely suppose that P N admits N simple roots z 1 , . . ., z N .Keeping the same notation, we have: Observe that deg(P N ) = N and since its leading coefficient is (−1) N det H (N ) 0 , we deduce that det H (N ) 0 ̸ = 0.In addition, ker(M (z k )) ̸ = {0} , for all k = 1, . . ., N, and we can verify that every nonzero element of ker(M (z k )) is proportional to (β 0 , . . ., β N −1 , 1) with: . . .
Now, writing = 0, and using the last line of this system, we obtain that z 1 , . . ., z N are the roots of the polynomial Since the coefficients of this polynomial satisfy (A.2), we can apply Algorithm 1.1 to construct a solution for which p = N , whose nodes are exactly z 1 , . . ., z N .Finally, since det(H (N ) 0 ) ̸ = 0, the uniqueness follows from Proposition 2.2.

Remark 4 . 5 .
Condition 3 in Theorem 4.4 holds, for example, when all the c k,N are positive, since N k=1 c k,N = |ω|, for all N ⩾ 1.

Figure 5 . 1 .
We represent on the first row of Figure 5.2 and Figure 5.3, the union of the disks centered at the nodes z k,N and with radii |Re(c k,N )|/π (if Re(c k,N