FACTORIZATION METHOD FOR THE INVERSE STOKES PROBLEM

. We propose an imaging technique for the detection of porous inclusions in a stationary ﬂow based on the Factorization method. The stationary ﬂow is described by the Stokes-Brinkmann equations, a standard model for a ﬂow through a (partially) porous medium, involving the deformation tensor of the ﬂow and the permeability tensor of the porous inclusion. On the boundary of the domain we prescribe Robin boundary conditions that provide the free- dom to model, e.g., in- or outlets for the ﬂow. The direct Stokes-Brinkmann problem to ﬁnd a velocity ﬁeld and a pressure for given boundary data is a mixed variational problem lacking coercivity due to the indeﬁnite pressure part. It is well-known that indeﬁnite problems are diﬃcult to tackle theoreti- cally using Factorization methods. Interestingly, the Factorization method can nevertheless be applied to this non-coercive problem, as long as one uses data consisting only of velocity measurements. We provide numerical experiments showing the feasibility of the proposed technique.


1.
Introduction. The inverse problem that we consider in this paper is to image penetrable inclusions modeled by the Stokes-Brinkmann equations in a bounded Lipschitz domain Ω ⊂ R d , d = 2 or 3, with Robin boundary conditions. These equations describe a fluid in steady-state in terms of the velocity field u = (u 1 , . . . , u d ) : Ω → R d and the pressure p : Ω → R. Both quantities are related through the deformation tensor D(u) of the fluid,  (2) div u = 0 in Ω.
We should mention already at this point that our notation does not respect some of the usual notational conventions in the physics literature for viscous flows. Most differences are intended to reduce the number of required symbols and allow for a compact notation. The Stokes-Brinkmann equations are particularly important, at various scales, for modeling flow through porous and partially porous media. Possible applications include for instance the modeling of flow of liquids or gas through the ground, see, e.g., [18,23] that investigate flow through carbonate karst reservoirs in detail. The problem of detecting inclusions in flows could moreover be applied to nondestructive testing problems of isolation materials by streaming air into a block of the material and measuring the resulting outflow.
On the boundary Γ := ∂Ω we prescribe Robin boundary conditions for the normal and tangential components of the velocity field. These conditions go back to Navier, see [22]; they will be introduced in (8) below. These conditions offer considerably more flexibility to model boundary phenomena than pure Dirichlet (or no-slip) or Neumann boundary conditions. The Robin conditions for instance allow to crudely model outlets on the boundary of the domain where liquid or gas streams in or out. However, since these boundary conditions necessarily link the pressure, the deformation tensor, and the velocity field, it is impossible to eliminate the pressure in (1) by, e.g., seeking merely for a velocity field in spaces of divergence-free functions (the latter technique would be applicable, e.g., when considering no-slip boundary conditions).
The explicit appearance of the pressure in the boundary conditions complicates the analysis of the direct and the inverse problem. First, the variational problem for the direct Stokes-Brinkmann problem necessarily yields a mixed, non-coercive variational formulation that is coercive in the velocity variable but indefinite in the pressure variable. Of course, this mixed formulation can be tackled using standard inf-sup conditions. Second, it is well-known that the analysis of the Factorization method requires, roughly speaking, coercivity of the direct problem. In our framework, the indefiniteness of the direct problem does not perturb the analysis of the Factorization method for the inverse problem, as long as one considers (difference) measurements that are merely based on the velocity field. The analysis that we present would indeed not work for any (reasonable) measurements involving pressure.
The number of publications on inverse problems for Stokes(-Brinkmann) equations and related partial differential equations is relatively small when compared to, e.g., the Laplace, the Helmholtz, or the Maxwell's equations, and most papers appeared quite recently. Among the work that is most related to ours is a paper on a point source method for the Oseen equation (a linearization of the Navier-Stokes equations), see [26], and a Factorization method for the Stokes equations, see [25], as well as non-linear integral equations for reconstructing obstacles in a Stokes flow, see [2]. In contrast to our paper, these works consider flows in an unbounded domain for impenetrable obstacles with no-slip (Dirichlet) boundary conditions for the velocity field. Usually, the viscosity µ outside the obstacle is constant, and the term div(µD(u)) becomes a multiple of the (vectorial) Laplacian. For such problems one can avoid to solve an indefinite problem involving the velocity field and the pressure, for example by solving an integral equation of the second kind for a single unknown density.
2. Stokes flow and variational formulation. The Stokes equations describe incompressible viscous fluids at low Reynolds number in the steady state. The governing equations for the velocity u : Ω → R d and the pressure p : with suitable boundary conditions of Robin type that we introduce below. The coefficient µ : Ω → R denotes the fluid's viscosity and f : Ω → R d is a source term. The second equation in (3) is the so-called incompressibility condition for the fluid. In contrast to many papers and textbooks on the Stokes equations, we do not assume here that the viscosity µ is constant. Note that a direct computation shows that for constant viscosity, the term div(µD(u)) reduces (in the distributional sense) to µ∆u /2, where the Laplacian is applied component-wise.
All differential equations in this text will be formulated variationally. To this end, the next lemma provides an essential integration-by-parts formula.
The right-most term has to be interpreted as a duality product of v ∈ H 1/2 (Γ) d with the co-normal derivative µD(u)n ∈ H −1/2 (Γ) d , see [21]. Seeking for a solution u ∈ H 1 (Ω) d to (3) we multiply the first equation with a test function v ∈ H 1 (Ω) d , integrate over Ω, and apply Lemma 2.2 as well as the divergence theorem, Let us now consider the boundary terms in (5). Every vector field u on the boundary can be decomposed into its tangential component v and its normal component v ⊥ , The same decomposition can be done for the vector field D(u)n, Remark 1. Strictly speaking, we would have to restrict a vector field u ∈ H 1 (Ω) d to the boundary Γ before taking its normal or tangential component. To simplify notation, we will nevertheless write u ⊥ and u instead of (u| Γ ) ⊥ and (u| Γ ) for u ∈ H 1 (Ω) d , respectively.
We want the tangential and normal components of u and of the deformation tensor D(u) to satisfy Robin-type boundary conditions, (8) µD (u) = g − α u on Γ and pn − µD ⊥ (u) = −g ⊥ + α ⊥ u ⊥ on Γ for given boundary data g ∈ L 2 (Γ) d . Plugging these relations into (7) yields that For given f ∈ L 2 (Ω) d and g ∈ L 2 (Γ) d , the variational formulation to find u ∈ H 1 (Ω) d and p ∈ L 2 (Ω) will be presented below, under the assumption that µ ∈ L ∞ (Ω) and α ⊥, ∈ L ∞ (Γ) are positive. Before stating this mixed formulation, let us introduce the Stokes-Brinkmann equations. A penetrable porous inclusion inside the viscous liquid can be modeled by the Stokes-Brinkmann equations: Using the above notation, these equations describe the steady-state flow via an additional absorption term M : Ω → R d×d that is typically supported inside the inclusions, −div(µD(u)) + M u + ∇p = f in Ω, (10) in Ω, The material parameters µ and M will typically take different values in the inside and the outside of a porous inclusion. In particular, outside the inclusion the flow is free, that is, M vanishes. To state existence and uniqueness results for the Stokes(-Brinkmann) equations, we suppose that µ ∈ L ∞ (Ω) is strictly positive, µ ≥ c > 0, that α ⊥, ∈ L ∞ (Γ) with α ⊥, ≥ c > 0, and that M ∈ L ∞ (Ω) d×d takes values in the symmetric and positive semidefinite matrices, Define the bounded bilinear forms for u, v ∈ H 1 (Ω) d , p ∈ L 2 (Ω), f ∈ L 2 (Ω) d and g ∈ L 2 (Γ) d . Repeating exactly the same integrations by parts that lead us to (9) we obtain the following mixed variational formulation of the Stokes-Brinkmann equations (10), which obviously include the Stokes equations by setting M = 0: Seek u ∈ H 1 (Ω) d and p ∈ L 2 (Ω) such that Theorem 2.3. For any bounded linear form F : Proof. This is an application of a standard result on mixed variational formulations, see [6]. For the problem under investigation, one uses Korn's inequality (see Lemma A.2(b)) to show coercivity of the bilinear form a(·, ·), and an inf-sup condition for the form b, see, e.g., [20,6].
Remark 2. Different physically senseful choices for the boundary conditions on Γ are possible, e.g., u = 0 and (µD ⊥ (u) − pn) + α ⊥ u ⊥ = g ⊥ on Γ. The ansatz space for the solution u is then H 1 ⊥ (Ω) d = {u ∈ H 1 (Ω) d : u = 0 on Γ}. The formulation (7) has to be adapted by replacing H 1 (Ω) d by H 1 ⊥ (Ω) d and skipping the term involving α . It is also possible to replace all Robin boundary conditions by Neumann conditions. Then, however, the solution space has to be adapted, typically yielding a quotient space where rigid motions are factored out.
3. Inverse Stokes problem and factorization of the data operator. In this section we consider an inverse Stokes problem which is to reconstruct an inclusion D ⊂ Ω where both viscosity and permeability differ from the background parameters from boundary measurements taken on Γ = ∂Ω. To this end, we denote by µ 0 the known viscosity of the background medium, while µ 1 denotes the unknown viscosity which differs from µ 0 inside an unknown inclusion D ⊂ Ω. In the entire paper we assume that Ω\D is connected. We also assume that the viscosity inside the inclusion is higher than in the background medium, and that the permeability inside D is positive. As above, we describe the permeability effect by a symmetric matrix M : Ω → R d×d such that M (x) = 0 a.e. in Ω\D and y M (x)y ≥ c|y| 2 > 0 for a.e. x ∈ D and all 0 = y ∈ R d . In the following, we abbreviate this property as Later on, it will be convenient to work with the viscosity contrast κ, The reconstruction of the shape of D is based on the difference of two Robin-to-Dirichlet operators. Roughly speaking, the first of these operators, Λ 0 : L 2 (Γ) d → L 2 (Γ) d , maps Robin boundary values to the trace of the solution u 0 to the Stokes problem with (background) viscosity µ 0 (and without the Brinkmann term M ). The second of these operators, Λ 1 : L 2 (Γ) d → L 2 (Γ) d , does the same for viscosity µ 1 and permeability M . More precisely, Remark 3. (a) In (15) we added the incompressibility condition from the second line of (13) to the first line of (13). Of course, choosing v = 0 in (15) shows that this condition is still satisfied. (b) The operator Λ 0 can in principle be computed (numerically) without knowing D, M , or µ 1 ; it is hence not necessarily required to measure Λ 0 experimentally. If, however, such measurements are available, then the modeling error of the difference Λ 0 − Λ 1 can be reduced significantly.
Analogously, the Robin-to-Dirichlet operator Λ 1 maps g to u 1 | Γ , where the pair Theorems 2.3 and A.1 show that Λ 0,1 are well-defined and compact operators on A crucial tool to arrive at a Factorization method for the inverse Stokes-Brinkmann problem is a factorization of the relative Robin-to-Dirichlet operator Λ 0 − Λ 1 . To construct this factorization we set Λ 0 g = u 0 | Γ and Λ 1 g = u 1 | Γ and write w = u 0 − u 1 and r = p 0 − p 1 . We take the difference of (15) and (16), use the fact that µ 0 = µ 1 − κ, and find that To derive a first factorization of Λ 0 − Λ 1 we define the (evaluation) operator H : (15). Additionally, we define the solution operator Comparing the last mixed variational formulation with (17) directly shows the factorization Λ 0 − Λ 1 = GH. In the following we will need to work with the adjoint H * of H. We compute this adjoint with respect to the (special) inner product Inverse Problems and Imaging Volume 7, No. 4 (2013), 1271-1293 Lemma 3.1. If κ ≥ c > 0 and if y M (x)y ≥ c|y| 2 > 0 a.e. in D and for all 0 = y ∈ R d , then the inner product from (19) is equivalent to the usual inner product on H 1 (D) d .
Proof. The proof follows from the coercivity of the bilinear form in (19) due to Lemma A.2(a).
By definition of the inner product (19) and by definition of H, (20).
Proof. This follows from the trace theorem, see Theorem A.1.
Proof. We exploit that µ 1 = µ 0 + κ to rewrite the variational formulation (18) that defines G(h) = w| Γ by taking all integrals containing κ or M on the right-hand side, Since we already know that Λ 0 − Λ 1 = GH, it remains to use the definition of H * in (20) to show that H * T = G to conclude that (21) holds true. Hence, we need to Obviously, (24) is exactly the same variational problem for the unknowns ( w, s) as (23) for (w, s). Since this problem is uniquely solvable (see Theorem 3.5. The operator T : (22), is self-adjoint on (H 1 (D), ·, · * ); see (19) for a definition of the inner product ·, · * .
Proof. For h 1,2 ∈ H 1 (Ω) d let (w 1,2 , r 1,2 ) ∈ H 1 (Ω) d × L 2 (Ω) be the solution to (22). We need to show that The variational formulation (22) for (w 2 , r 2 ) with (w 1 , r 1 ) as test functions shows that the left-hand side of (25) equals The variational formulation (22) for (w 1 , r 1 ) with (w 2 , r 2 ) as test functions shows that the right-hand side in (25) equals the same expression, up to exchanging w 1 and w 2 . Hence, (25) holds true and T is self-adjoint. Proof. The definition of ·, · * in (19) and the definition of T in Theorem 3.4 imply that Obviously, Finally choosing the test function q to be r, (26) becomes The task is now to estimate the terms of the right of (27). First, Second, the coercivity of the bilinear form a(·, ·) on H 1 (Ω) d implies Consequently, Proof. The factorization Λ 0 − Λ 1 = H * T H and the compactness of H * imply that the difference Λ 0 − Λ 1 is compact; Theorem 3.5 and the factorization yield the selfadjointness of the relative Robin-to-Dirichlet operator. The coercivity of T implies that (Λ 0 − Λ 1 )g, g L 2 (Γ) d ≥ 0 for all g ∈ L 2 (Γ) d . Hence, Λ 0 − Λ 1 is non-negative.
(This follows from the same computation as in the last part). Obviously, Hg = 0, since u 0 vanishes on Ω. Since H and Λ 0 − Λ 1 are related by the factorization of Theorem 3.4, we deduce that (Λ 0 − Λ 1 )g = 0.
Since the relative Robin-to-Dirichlet operator is a compact and self-adjoint operator, there exists a complete orthonormal eigensystem (λ l , π l ) l∈N of Λ 0 − Λ 1 such that λ l → 0 as l → ∞. Since Λ 0 − Λ 1 is non-negative, it is obvious that λ l ≥ 0 and the characterization of the kernel of Λ 0 − Λ 1 shows that precisely one eigenvalue vanishes. For simplicity, we define λ 1 = 0 to be the zero eigenvalue, and order the remaining eigenvalues λ 2 ≥ λ 3 ≥ · · · > 0 in decreasing order.

Fundamental solution and Green's function of the Stokes equations.
The second important ingredient of the Factorization method is a characterization of the inclusion D via the Green's function of the Stokes equations. In this section, we construct this Green's function via a correction of the boundary values of the fundamental solution of the Stokes equations. To this end, we will assume that the viscosity µ 0 of the background medium is constant. Without loss of generality, we choose the simplest case, µ 0 = 2 in Ω, since then div(µ 0 D(u)) = ∆u for any divergence-free function due to Lemma 2.1. It is not difficult to see that any other (positive) constant viscosity can be treated analogously by a simple scaling argument. The treatment of varying background viscosities requires more involved adaptions of the fundamental solution to get the Green's function for varying viscosity, see, e.g., [13] for the case of the Laplace equation.
The fundamental solution of the Stokes equations can be found in many references on boundary integral equations for the Stokes equations, see, e.g., [10,19,15]. We will concentrate in this section on the case d = 2 and announce the results for d = 3 in the end of the section without proof. For x = y ∈ R 2 , we set The Green's function is now defined to correct the inhomogeneous boundary values of the fundamental solution to homogeneous Robin boundary conditions. For y ∈ Ω, define correction terms where R ij ∈ H 1 (Ω), R i = (R 1i , R 2i ) , and r i ∈ L 2 (Ω). If, for i = 1, 2 and y ∈ Ω, the pair (R i , r i ) ∈ H 1 (Ω) 2 × L 2 (Ω) is a variational solution to the Stokes problem (30) − ∆R i (·, y) + ∇r i (·, y) = 0 and divR i (·, y) = 0 in Ω, subject to the boundary conditions (31) 2D (R i (·, y)) + α R i (·, y) = −[2D (S i (·, y)) + α S i (·, y) ] on Γ and (r i (·, y)n−2D ⊥ (R i (·, y))n) − α ⊥ R i (·, y) ⊥ (32) = −[(s i (·, y)n − 2D ⊥ (S i (·, y))n) − α ⊥ S i (·, y) ⊥ ] on Γ, then (E, e) = (S, s) + (R, r) is called the Green's function for the Stokes equations in Ω. Note that all right-hand sides in the above boundary conditions for (R i , r i ) belong to L 2 (Γ) 2 , since the singularity y ∈ Ω of the fundamental solution is away from the boundary Γ. The explicit variational formulation defining (R i , r i ) can be found from (13), choosing µ ≡ 2, Again, we emphasize that the 2 × 2 matrix-valued function E consists of the two columns E 1 and E 2 . Obviously, the Green's function (E, e) solves, for i = 1, 2, The following sampling function in two dimensions plays a main role for the construction of the Factorization method in the next section. For every z ∈ Ω and some t ∈ [0, 2π) we define The parameter t could in principle be chosen in dependence on the sampling point z.
However, in all our numerical examples, we will always work with a fixed t, for simplicity. The sampling function φ z characterizes the inclusion D via the range R(H * ) of the adjoint H * : L 2 (D) 2×2 → L 2 (Γ) 2 , see Theorem 3.2.
Theorem 4.1. For z ∈ Ω it holds that z ∈ D if and only if φ z ∈ R(H * ).
Due to the divergence constraint for the velocity field, the proof of Theorem 4.1 is substantially more involved than the corresponding result for, e.g., the Laplace equation. The easy direction of the proof is based on the following auxiliary result that we will not prove in detail since the arguments are rather standard, see, e.g., [16]. In a nutshell, the proof of Lemma 4.2 follows from the definition of E i = S i + R i , from the fact that R i ∈ H 1 (Ω) 2 by definition, and from the singularity of S i at x = y, which prevents ∇S i from being square integrable in Ω. Admitting this result, we next prove Theorem 4.1.
The two-dimensional vector potential used in the proof of Theorem 4.1 has to be replaced by a corresponding three-dimensional construction that can be found in, e.g., [24]. The rest of the proof essentially remains the same, up to adapting the space dimension.

5.
The factorization method. Theorem 4.1 explicitly characterizes the inclusion D ⊂ Ω using the range of the operator H * . This operator is, however, unknown because it cannot be computed without knowing the inclusion D. Fortunately, the factorization Λ 0 − Λ 1 = H * T H from Theorem 3.4 relates H * with the known data Λ 0 − Λ 1 . A range identity will then provide to a characterization of the inclusion D in terms of Λ 0 − Λ 1 in the main result of this paper, see Theorem 5.1. We exploit the eigendecomposition λ l g, π l L 2 (Γ) d π l , g ∈ L 2 (Γ) d , of Λ 0 − Λ 1 from Theorem 3.7 to define the (positive semi-definite) square root This square root is again a bounded self-adjoint operator on L 2 (Γ) d and Now we have two factorizations of the relative Robin-to-Dirichlet operator on L 2 (Γ) d . Since T is self-adjoint and coercive, we can apply the range characterization from [17,Lemma 5.15] (see Lemma A.5 in the appendix) to the first factorization. The result is that for any 0 = φ ∈ L 2 (Γ) d there holds Applying the same result to the second factorization Together, (41) and (42) imply the range identity Since (λ l , π l ) l∈N is an eigensystem of Λ 0 − Λ 1 and since λ l ≥ 0, it is obvious that ( √ λ l , π l ) l∈N is an eigensystem of (Λ 0 −Λ 1 ) 1/2 . Picard's range criterion [17, Theorem A.54] then implies the following characterization of D via the sampling functions φ z defined in (34) and in (40) for the case of dimension two and three, respectively.
Proof. Theorem 4.1 and the range identity (43) imply that z ∈ D if and only if φ z ∈ R((Λ 0 −Λ 1 ) 1/2 ). Recall that in last paragraph of Section 3 we have set λ 1 = 0, such that the first eigenfunction π 1 spans the kernel of (Λ 0 −Λ 1 ) 1/2 . Moreover, since φ z is the trace of a linear combination of the divergence-free columns E 1,...,d (·, z) of the Green's function, the divergence theorem states that Γ n · φ z dS = 0. In consequence, the characterization of the kernel of Λ 0 − Λ 1 implies that φ z , π 1 L 2 (Γ) d = 0. Thus, Hence, Picard's lemma shows that the series on the right of (44) is finite if and only if φ z ∈ R((Λ 0 − Λ 1 ) 1/2 ), that is, if and only if z ∈ D.

Numerical examples.
In this section we present numerical experiments illustrating the theoretical results of the last sections and the feasibility of the method. We present images that are computed from relatively few normal measurements that are additionally perturbed by additive noise. For simplicity, all examples are two-dimensional. Further, all experiments are computed using synthetic data computed by a mixed finite element method based on Taylor-Hood (P2/P1) elements, see, e.g., [8,6]. All finite element simulations are done using the C++-based open source finite element code FreeFem++ (see www.freefem.org/ff++ ). For our numerical examples we decompose the boundary of the domain Ω into N connected pieces Γ j , j = 1, . . . , N of equal length and approximate functions on the boundary by piecewise constant (vector-valued) functions. We will denote this approximation space by X (N ) ; later on, we will always choose N = 16. This rather low-dimensional choice is motivated by the difficulty to experimentally measure many degrees of freedom of a boundary flow in practice. Valves to measure or control the in-and outflow require sufficiently large holes in the boundary of the domain and cannot easily be shifted. Note that the mesh of the domain that we use to compute the synthetic data has many degrees of freedom on each boundary piece such that mean-value projections have to be used to project the boundary values of general finite-element functions into X (N ) . It is by the way not necessary that the boundary pieces cover the entire boundary, even if in all our examples this will be the case.
To approximate the difference Λ 0 − Λ 1 , we firstly discretize (15) using Taylor-Hood elements to evaluate Hg = u 0 | D on D, for all g ∈ X (N ) . Secondly, we exploit that (Λ 0 −Λ 1 )g = GHg, inserting (Hg)| D in the right-hand side of the discretization of (18) The background viscosity equals µ 0 ≡ 2 and µ = 100 and 2 inside and outside the inclusion D, respectively. The (scalar) material parameter M equals 50 inside the inclusion and vanishes outside. At first glance, these values seem rather big. However, Figure 1(a) shows that the resulting inclusion is well penetrated by a flow that is generated by an inflow on the right and an outflow on the left. The boundary parameters for this example, as well as for all the following ones are α = 5 and α ⊥ = 10. Moreover, the plot of the first component of this flow in Figure 1(b) shows that the velocity field is influenced by the inclusion. Of course, this influence decreases if one decreases the contrast; when setting, e.g., µ = 10 and M = 5 inside D, the flow is no longer visibly influenced.
For the configuration explained in the last paragraph, the eigenvalues of the numerical approximation to (Λ 0 − Λ 1 ) (16) that we use for our inversion experiments are shown (in decreasing order) in Figure 1(c). The fact that one eigenvalue of Λ 0 − Λ 1 vanishes is reflected by the large gap in magnitude between the smallest eigenvalue (roughly equal to 5e − 18) and all other eigenvalues. Of course, the precision of the computed eigenvalues depends on the spatial discretization of Ω. We test this precision by plotting the largest eigenvalue for different meshes in Figure 1(d). The parameter on the horizontal axis, say k, describes the maximal mesh width h k = π/(8k). We observe convergence of the largest eigenvalue as k increases: 2, 3, and 4 correct digits are attained for k = 6, 11, and 22, respectively. We use k = 12 for generating the synthetic data later on.
Since the Green's function of the Stokes equations in Ω has in general no explicit representation we first compute the correction term (R, r) by discretizing (30-32). Second, we compute the Green's function by adding the (analytic evaluation) of the fundamental solution function and the numerical approximation to the correction term to construct the sampling function φ z from (34). Computing the Green's function requires a finer mesh than computing the synthetic data, at least when the sampling point is close to the boundary. We found that k = 39 gives sufficient accuracy for reasonable images. Despite we could in principle compute testfunctions φ z for several dipole directions (parametrized by t in (34)) we restrict ourselves to fix t = 0. Improving the reconstruction quality by, e.g., averaging or optimizing over dipole directions might be possible; our aim here is merely to give a proof of concept for the method.  (16) (no artificial noise was added to the data). Red diamonds: 32 singular values of a random additive perturbation of (Λ 0 − Λ 1 ) (16) , the relative noise level equals 1%. (d) Blue dots: the largest eigenvalue of (Λ 0 − Λ 1 ) (16) computed on different meshes. The eigenvalues are plotted against an integer k determining the mesh width h k = π/(8k).
The sampling points z determining the sampling function φ z are chosen on a uniform grid G of Ω. We compute the sampling function φ z from (34) for every sampling point (without projecting the boundary data of the fundamental solution into X (N ) ). Note that this requires to solve the auxiliary problem (30-32) once for every sampling point. However, the sampling function does not depend on the inclusion and can be precomputed and re-used for several reconstructions in the same setting. After assembling the sampling function φ z we project the result into X (N ) using a mean-value projection and define Tangential and normal components φ are then defined componentwise as in (6).
Finally, to identify an inclusion D in a domain Ω we approximate the reciprocal of the function from (44) by where 1 ≤ N − < N + ≤ dN − 1 are truncation indices. Roughly speaking, if the approximation quality of our discrete data is good enough, then Theorem 5.1 motivates that this procedure yields an image of the inclusion, since the function in (45) should be small outside D while taking large values inside D. Of course, the quality of the discrete data determines the resolution and the contrast of the resulting image. In our numerical experiments, we noted that the tangential component of the measurements Λ 0 −Λ 1 is fairly small, usually about two orders of magnitude smaller than the normal component. This observation holds both for tangential and normal boundary excitations. To this end, in our subsequent examples we merely take the "normal" part of Λ 0 − Λ 1 as data, that is, we only use normal excitations and only measure the normal component of the trace of the resulting flow. Denote by X (N ) ⊥ = {g ∈ X (N ) , g · n = 0 on Γ} the space of tangential vector fields in X (N ) . Replacing the space X (N ) in the above construction of (Λ 0 − Λ 1 ) (N ) we obtain a corresponding matrix representation (Λ 0 − Λ 1 ) ⊥,j , the corresponding eigenvalues and vectors, the indicator function that we plot is hence In our experiments, the reconstructions obtained merely from normal measurements gave better results than those obtained from all measurements. This is presumably due to the large difference in magnitude of the normal and tangential data, which causes the tangential measurements to be inexact even when no artificial noise has been added. It made relatively little difference whether the excitations were chosen normally or tangentially, as long as we took normal measurements of the velocity field. Such normal measurements are, by the way, much simpler to measure reliably in a practical experiment. In our implementation, we actually replace the factor 1/λ (N ) ⊥,l in (46) by 1/(λ (N ) ⊥,l + 2 ) where > 0 corresponds to the relative noise level, to avoid division by zero.
(There are definitely more sophisticated numerical techniques for regularizing the Picard series, see, e.g., [11,12].) In our numerical experiments we noted that the image quality can be significantly improved by taking the freedom to vary both indices N ± . For the images shown in Figure 2 we always chose N − = 10 and N + = 15. This choice produced the best overall results. Figures 2(b,c,d) show the images of the inclusion D using W ⊥ (z) from synthetic data (Λ 0 − Λ 1 ) (16) ⊥⊥ generated in the setting explained above. In Figures 2(c,d) we added a random matrix (generated using uniformly distributed random variables) to the discretization of Λ 0 − Λ 1 . The noise was scaled such that the relative noise level equals 1% in (c) and 5% in (d). The experiments show that the method is able to separate the two inclusions up to a noise level of about 1%. However, even without artificial noise, the shape of the two inclusions is not reconstructed correctly, in contrast to position and size. This might be surprising when compared to reconstructions in impedance tomography, see, e.g., [11], but is probably due to the coarse approximation space X (16) that we employ for approximating functions on the boundary. The separation of the two inclusions is lost when the noise level increases to 5%, but geometric information on the inclusions can still be deduced from the corresponding image.
Appendix A. Auxiliary results. As in the main part of this text, Ω is a Lipschitz domain with boundary Γ. The following trace theorem is a well-known analytic tool, see, e.g., [21,9].
Note that we do not denote the trace operator explicitly in this text, but rather write u| Γ . Often, we even omit to denote the trace operation completely if it follows from the context (as, e.g., inside boundary integrals). For the following Korn's inequalities, we refer to, e.g., [20,21,6].  for all v ∈ H 1 (Ω) d and q ∈ L 2 (Ω) with compact support in Ω. This variational formulation can of course also be formulated using ∇u instead of the deformation tensor D(u) due to Lemma 2.1.
The next lemma extends the analytic continuation result from [17, Lemma 6.13] from the Helmholtz equation to the Stokes equations. The main tool to generalize the latter result are Green's first and second identity for the Stokes equations, see, e.g., [15,Chapter 2.3].
Then u can be extended to an analytic function in Ω and (u, p) satisfies the Stokes equations in Ω.
The next lemma is a well-known characterization of certain operator ranges, see, e.g., [17,Lemma 5.15].
Lemma A.5. (Inf-criterion) Let X and Y be Hilbert spaces, B : X → X, A : X → Y , and T : Y → Y linear and bounded such that B = A * T A. Furthermore, let T be self-adjoint and coercive, i.e., T y, y Y ≥ c y 2 Y for all y ∈ Y . Then it holds for any 0 = ϕ ∈ X that ϕ ∈ R(A * ) ⇔ inf{ Bx, x X : x ∈ X, x, ϕ X = 1} > 0.