Unique continuation for the Helmholtz equation using stabilized finite element methods

In this work we consider the computational approximation of a unique continuation problem for the Helmholtz equation using a stabilized finite element method. First conditional stability estimates are derived for which, under a convexity assumption on the geometry, the constants grow at most linearly in the wave number. Then these estimates are used to obtain error bounds for the finite element method that are explicit with respect to the wave number. Some numerical illustrations are given.


Introduction
We consider a unique continuation (or data assimilation) problem for the Helmholtz equation and introduce a stabilized finite element method (FEM) to solve the problem computationally. Such methods have been previously studied for the Laplace equation in [5], [6] and [8], and for the heat equation in [9]. The main novelty of the present paper is that our method is robust with respect to the wave number k, and we prove convergence estimates with explicit dependence on k, see Theorems 1 and 2 below.
An abstract form of a unique continuation problem is as follows. Let ω ⊂ B ⊂ Ω be open, connected and non-empty sets in R n and suppose that u ∈ H 2 (Ω) satisfies (1) in Ω. Given u in ω and f in Ω, find u in B.
In general, the constants C and α depend on k, as can be seen in Example 3 given in Appendix A. However, under suitable convexity assumptions on the geometry and direction of continuation it is possible to prove that in (2) both the constants C and α are independent of k, see the uniform estimate in Corollary 2 below, which is closely related to the so-called increased stability for unique continuation [16]. Obtaining optimal error bounds in the finite element approximation crucially depends on deriving estimates similar to (2), with weaker norms in the right-hand side, as in Corollary 3 below, or in both sides, by shifting the Sobolev indices one degree down, as in Lemma 2 below.
In addition to robustness with respect to k, an advantage of using stabilized FEM for this unique continuation problem is that, when designed carefully, its implementation does not require information on the constants C and α in (2), or any other quantity from the continuous stability theory, such as a specific choice of a Carleman weight function. Moreover, unlike other techniques such as Tikhonov regularization or quasi-reversibility, no auxiliary regularization parameters need to be introduced. The only asymptotic parameter in our method is the size of the finite element mesh, and in particular, we do not need to saturate the finite element method with respect to an auxiliary parameter as, for example, in the estimate (34) in [4].
Throughout the paper, C will denote a positive constant independent of the wave number k and the mesh size h, and which depends only on the geometry of the problem. By A B we denote the inequality A ≤ CB, where C is as above.
For the well-posed problem of the Helmholtz equation with the Robin boundary condition (3) ∆u + k 2 u = −f in Ω and ∂ n u + iku = 0 on ∂Ω, the following sharp bounds (4) ∇u L 2 (Ω) + k u L 2 (Ω) ≤ C f L 2 (Ω) and (5) hold for a star-shaped Lipschitz domain Ω and any wave number k bounded away from zero [3]. The error estimates that we derive in Section 3, e.g. u − u h H 1 (B) ≤ C(hk) α u * in Theorem 2, contain the term which corresponds to the well-posed case term k f L 2 (Ω) . It is well known from the seminal works [2,17,18] that the finite element approximation of the Helmholtz problem is challenging also in the well-posed case due to the so-called pollution error. Indeed, to observe optimal convergence orders of H 1 and L 2 -errors the mesh size h must satisfy a smallness condition related to the wave number k, typically for piecewise affine elements, the condition k 2 h 1. This is due to the dispersion error that is most important for low order approximation spaces. The situation improves if higher order polynomial approximation is used. Recently, the precise conditions for optimal convergence when using hp-refinement (p denotes the polynomial order of the approximation space) were shown in [22]. Under the assumption that the solution operator for Helmholtz problems is polynomially bounded in k, it is shown that quasi-optimality is obtained under the conditions that kh/p is sufficiently small and the polynomial degree p is at least O(log k).
Another way to obtain absolute stability (i.e. stability without, or under mild, conditions on the mesh size) of the approximate scheme is to use stabilization. The continuous interior penalty stabilization (CIP) was introduced for the Helmholtz problem in [24], where stability was shown in the kh 1 regime, and was subsequently used to obtain error bounds for standard piecewise affine elements when k 3 h 2 1. It was then shown in [10] that, in the one dimensional case, the CIP stabilization can also be used to eliminate the pollution error, provided the penalty parameter is appropriately chosen. When deriving error estimates for the stabilized FEM that we herein introduce, we shall make use of the mild condition kh 1. To keep down the technical detail we restrict the analysis to the case of piecewise affine finite element spaces, but the extension of the proposed method to the high order case follows using the stabilization operators suggested in [5] (see also [7] for a discussion of the analysis in the ill-posed case).
From the point of view of applications, unique continuation problems often arise in control theory and inverse scattering problems. For instance, the above problem could arise when the acoustic wave field u is measured on ω and there are unknown scatterers present outside Ω.

Continuum stability estimates
Our stabilized FEM will build on certain variations of the basic estimate (2), with the constants independent of the wave number, and we derive these estimates in the present section. The proofs are based on a Carleman estimate that is a variation of [16, Lemma 2.2] but we give a self-contained proof for the convenience of the reader. In [16] the Carleman estimate was used to derive the so-called increased stability estimate where F = u H 1 (ω) + ∆u + k 2 u L 2 (Ω) and the constants C and α are independent of k.
Here F can be interpreted as the size of the data in the unique continuation problem and the H 1 -norm of u as an a priori bound. As k grows, the first term on the right-hand side of (7) dominates the second one, and the stability is increasing in this sense.
As our focus is on designing a finite element method, we prefer to measure the size of the data in the weaker norm Taking u to be a plane wave solution to (1) suggests that , could be the right analogue of (7) when both the data and the a priori bound are in weaker norms. We show below, see Lemma 2, a stronger estimate with only the second term on the right-hand side.
Lemma 1 below captures the main step of the proof of our Carleman estimate. This is an elementary, but somewhat tedious, computation that establishes an identity similar to that in [21] where the constant in a Carleman estimate for the wave equation was studied. For an overview of Carleman estimates see [20,23], the classical references are [14,Chapter 17] for second order elliptic equations, and [15, Chapter 28] for hyperbolic and more general equations. In the proofs, the idea is to use an exponential weight function e ℓ(x) and study the expression ∆(e ℓ w) = e ℓ ∆w + lower order terms, or the conjugated operator e −ℓ ∆e ℓ . A typical approach is to study commutator estimates for the real and imaginary part of the principal symbol of the conjugated operator, see e.g. [20]. This can be seen as an alternative way to estimate the cross terms appearing in the proof of Lemma 1. Sometimes semiclassical analysis is used to derive the estimates, see e.g. [20]. This is very convenient when the estimates are shifted in the Sobolev scale, and we will use these techniques in Section 2.2 below.
Proof. We employ the equality in Lemma 1 with ℓ = τ φ and σ = ∆ℓ. With this choice of σ, it holds that a = 0. As the two first terms on the right-hand side of the equality are positive, it is enough to consider The claim follows by combining this with and τ |∇(∆φ)||∇v||v| ≤ C(|∇v| 2 + τ 2 |v| 2 ).
The above Carleman estimate implies an inequality that is similar to the three solid balls inequality, see e.g. [1]. The main difference is that here the foliation along spheres is followed in the opposite direction, i.e. the convex direction.
For simplicity we consider here a specific geometric setting, illustrated in Figure 1. Figure 1. Geometric setting. We use the following notation for a half space Corollary 2. Let r > 0, β > 0, R > r and r 2 + β 2 < ρ < R 2 + β 2 . Define y = (β, 0, . . . , 0) and Then there are C > 0 and α ∈ (0, 1) such that for all u ∈ C 2 (Ω) and k ≥ 0 Then φ is smooth and strictly convex in Ω, and it does not have critical points there.

2.2.
Shifted three balls inequality. In this section we prove an estimate as in Corollary 2, but with the Sobolev indices shifted down one degree, and our starting point is again the Carleman estimate in Corollary 1. When shifting Carleman estimates, as we want to keep track of the large parameter τ , it is convenient to use the semiclassical version of pseudodifferential calculus. We write > 0 for the semiclassical parameter that satisfies = 1/τ . The semiclassical (pseudo)differential operators are (pseudo)differential operators where, roughly speaking, each derivative is multiplied by , for the precise definition see Section 4.1 of [25]. The scale of semiclassical Bessel potentials is defined by and the semiclassical Sobolev spaces by . Then a semiclassical differential operator of order m is continuous from H m+s scl (R n ) to H s scl (R n ), see e.g. Section 8.3 of [25].
We will give a shifting argument that is similar to that in Section 4 of [11]. To this end, we need the following pseudolocal and commutator estimates for semiclassical pseudodifferential operators, see e.g. (4.8) and (4.9) of [11]. Suppose that ψ, χ ∈ C ∞ 0 (R n ) and that χ = 1 near supp(ψ), and let A, B be two semiclassical pseudodifferential operators of orders s, m, respectively. Then for all p, q, N ∈ R, Both these estimates follow from the composition calculus, see e.g. [25,Theorem 4.12].
Let φ be as in Corollary 1 and set ℓ = φ/ and σ = ∆ℓ in Lemma 1. Then where Ω ′ ⊂ R n is open and bounded, and Ω ⊂ Ω ′ . Then, rescaling by 4 , and for small enough > 0 we obtain . Now the conjugated operator P is a semiclassical differential operator, where we used the pseudolocality (11) to absorb the second term on the right-hand side by the left-hand side. We have and using the commutator estimate (12), we have . This can be absorbed by the left-hand side of (13). Thus √ . Take now s = −1, and let φ and χ be as in the proof of Corollary 2. In fact, we need to be slightly more careful and choose χ so that there is ψ ∈ C ∞ 0 (B(y, s) ∪ ω) satisfying ψ = 1 in supp([P, χ]).
Let u ∈ C ∞ (R n ) and set v = e φ/ χu. Write also w = e φ/ u.
Similarly with the proof of Corollary 2, , and, fixing small > 0, we arrive to . Note that for fixed , the norms of H −1 scl (Ω) and H −1 (Ω) are equivalent. Hence we have proved the following result.

Stabilized finite element method
We aim to solve the unique continuation problem for the Helmholtz equation (14) ∆u where ω ⊂ Ω ⊂ R 1+n are open, f ∈ H −1 (Ω) and q ∈ L 2 (ω) are given. Following the optimization based approach in [5,6,8,9] we will make use of the continuum stability estimates in Section 2 when deriving error estimates for the finite element approximation.
3.1. Discretization. Consider a family T = {T h } h>0 of triangulations of Ω consisting of simplices such that the intersection of any two distinct ones is either a common vertex, a common edge or a common face. Also assume that the family T is quasi-uniform. Let be the H 1 -conformal approximation space based on the P 1 finite element and let and the Scott-Zhang interpolator π h : H 1 (Ω) → V h , that preserves vanishing Dirichlet boundary conditions. Both operators have the following stability and approximation properties, see e.g. [12, Chapter 1] where i = π, Π, k = 1, 2 and m = 0, k − 1.
The regularization on the discrete level will be based on the jump stabilizer where F h is the set of all internal faces, and the jump over F ∈ F h is given by with K 1 , K 2 ∈ T h being two simplices such that K 1 ∩ K 2 = F , and n j the outward unit normal of K j , j = 1, 2.
Proof. See [9, Lemma 2] when the interpolator is π h . Since this proof uses just the approximation properties of π h , it holds verbatim for Π h .
Adopting the notation a(u, z) = (∇u, ∇z) L 2 (Ω) , G f (u, z) = a(u, z) − k 2 (u, z) L 2 (Ω) − f, z , G = G 0 we have that for u ∈ H 1 (Ω) G f (u, z) = 0, z ∈ H 1 0 (Ω), is the weak formulation of ∆u + k 2 u = −f . Our approach is to minimize the Lagrangian functional where s and s * are symmetric bilinear forms-the primal and dual stabilizers-and · ω denotes the norm · L 2 (ω) . We make the choice s(u, u) = J(u, u) + hk 2 u 2 L 2 (Ω) , s * = a, and define on V h and W h , respectively, the norms where A is the symmetric bilinear form Since A[(u, z), (u, −z)] = u 2 ω + u 2 V + z 2 W we have the following inf-sup condition (20) sup that guarantees a unique solution in V h × W h for (19).

3.2.
Error estimates. We start by deriving some lower and upper bounds for the norm · V . For u h ∈ V h , z ∈ H 1 0 (Ω), we use (17) and hence For u ∈ H 2 (Ω), from (18) and the stability of the L 2 -projection where u * is defined as in (6).

Lemma 4.
Let u ∈ H 2 (Ω) be the solution to (14) and (u h , z h ) ∈ V h × W h be the solution to (19). Then there exists C > 0 such that for all h ∈ (0, 1) Proof. Due to the inf-sup condition (20) it is enough to prove that The weak form of (14) implies that Using (16) we bound the first term to get For the second term we use the L 2 -orthogonality property of Π h , and (16) to obtain while for the last term we employ (22) to estimate Theorem 1. Let ω ⊂ B ⊂ Ω be defined as in Corollary 2. Let u ∈ H 2 (Ω) be the solution to (14) and (u h , z h ) ∈ V h × W h be the solution to (19). Then there are C > 0 and α ∈ (0, 1) such that for all k, h > 0 with kh 1 Proof. Consider the residual r, w = G(u h − u, w) = G(u h , w) − f, w , w ∈ H 1 0 (Ω). Taking v = 0 in (19) we get G(u h , w) = f, w + s * (z h , w), w ∈ W h which implies that . Using (21) and (16) we estimate the first term

since, due to Lemma 4 and (22)
The second term is bounded by using (16) and the last term by using Lemma 4 and the H 1 -stability (15) Hence the following residual norm estimate holds Using the continuum estimate in Lemma 2 for u − u h we obtain the following error estimate By (16) and Lemma 4 we have the bounds ≤ Ck −2 u * thus leading to the conclusion.

Theorem 2.
Let ω ⊂ B ⊂ Ω be defined as in Corollary 2. Let u ∈ H 2 (Ω) be the solution to (14) and (u h , z h ) ∈ V h × W h be the solution to (19). Then there are C > 0 and α ∈ (0, 1) such that for all k, h > 0 with kh 1 Proof. We employ a similar argument as in the proof of Theorem 1 with the same estimates for the residual norm and the L 2 -errors in ω and Ω, only now using the continuum estimate in Corollary 3 to obtain which ends the proof.
Let us remark that if we make the assumption k 2 h 1 then the estimate in Theorem 2 becomes and combining Theorems 1 and 2 we obtain the following result.

Corollary 4.
Let ω ⊂ B ⊂ Ω be defined as in Corollary 2. Let u ∈ H 2 (Ω) be the solution to (14) and (u h , z h ) ∈ V h × W h be the solution to (19). Then there are C > 0 and α ∈ (0, 1) such that for all k, h > 0 with k 2 h 1 Comparing with the well-posed boundary value problem (3) and the sharp bounds (4) and (5), we note that the k −1 u * term in the above estimate is analogous to the well-posed case term f L 2 (Ω) .
The critical points (u, z) ∈ V h × W h of the perturbed Lagrangian Lq ,f satisfy Lemma 5. Let u ∈ H 2 (Ω) be the solution to the unperturbed problem (14) and (u h , z h ) ∈ V h × W h be the solution to the perturbed problem (23). Then there exists C > 0 such that for all h ∈ (0, 1) Proof. Proceeding as in the proof of Lemma 4, the weak form gives We bound the perturbation terms by and we conclude by using the previously derived bounds for the other terms.
Theorem 3. Let ω ⊂ B ⊂ Ω be defined as in Corollary 2. Let u ∈ H 2 (Ω) be the solution to the unperturbed problem (14) and (u h , z h ) ∈ V h × W h be the solution to the perturbed problem (23). Then there are C > 0 and α ∈ (0, 1) such that for all k, h > 0 with kh 1 Proof. Following the proof of Theorem 1, the residual satisfies Bounding the first term in the right-hand side by Lemma 5 and (22) and the third one by Lemma 5 again, we obtain The continuum estimate in Lemma 2 applied to u − u h gives where u − u h L 2 (ω) was bounded by using Lemma 5 and (16). Then the bound concludes the proof.

Theorem 4.
Let ω ⊂ B ⊂ Ω be defined as in Corollary 2. Let u ∈ H 2 (Ω) be the solution to the unperturbed problem (14) and (u h , z h ) ∈ V h × W h be the solution to the perturbed problem (23). Then there are C > 0 and α ∈ (0, 1) such that for all k, h > 0 with kh 1 Proof. Following the proof of Theorem 3, we now use Corollary 3 to derive which ends the proof.
Analogous to the unpolluted case, if k 2 h 1 the above result becomes and combining Theorems 3 and 4 gives the following.

Corollary 5.
Let ω ⊂ B ⊂ Ω be defined as in Corollary 2. Let u ∈ H 2 (Ω) be the solution to the unperturbed problem (14) and (u h , z h ) ∈ V h × W h be the solution to the perturbed problem (23). Then there are C > 0 and α ∈ (0, 1) such that for all k, h > 0 with k 2 h 1

Numerical examples
We illustrate the above theoretical results for the unique continuation problem (14) with some numerical examples. Drawing on previous results in [5], we adjust the primal stabilizer in (19) with a fixed stabilization parameter γ > 0 such that s(u, v) = γJ(u, v) + γh 2 k 4 (u, v) L 2 (Ω) . The error analysis stays unchanged under this rescaling. Various numerical experiments indicate that γ = 10 −5 is a near-optimal value for different kinds of geometries and solutions. The implementation of our method and all the computations have been carried out in FreeFem++ [13]. Ω is the unit square, and the triangulation is uniform with alternating left and right diagonals, i.e. of Union Jack type, as shown in Figure 2. The mesh size is taken as the inverse square root of the number of nodes.
In the light of the convexity assumptions in Section 2, we shall consider two different geometric settings: one in which the data is continued in the convex direction, inside the convex hull of ω, and one in which the solution is continued in the non-convex direction, outside the convex hull of ω. We will give results for two kinds of solutions.
In the convex setting, given in Figure 3a, we take and Figure 4a shows that, for k = 10, the numerical results strongly agree with the convergence rates expected from Theorems 1 and 2, and Lemma 4, i.e. sub-linear convergence for the relative error in the L 2 and H 1 -norms, and quadratic convergence for J(u h , u h ). Although in Figure 4b we do obtain smaller errors and better than expected convergence rates when k = 50, various numerical experiments indicate that this example's behaviour when increasing the wave number k is rather a particular one. For oscillatory solutions, such as those in Example 2, with fixed n, or the homogeneous u = sin(kx/ √ 2) cos(ky/ √ 2), we have noticed that the stability deteriorates when increasing k.
In the non-convex setting we let           Let us recall that stability estimates for the unique continuation problem are closely related to those for the notoriously ill-posed Cauchy problem, see e.g. [1]. It is of interest to consider the following variation of a well-known example due to Hadamard, since this example can be used to show that conditional Hölder stability is optimal for the unique continuation problem.
It can be seen in Figure 7a that the convergence rates agree with the ones predicted for the convex setting i.e. sub-linear convergence for the relative error in the L 2 and H 1 -norms, and quadratic convergence for J(u h , u h ), although one can notice that the values of the jump stabilizer J(u h , u h ) visibly increase compared to Example 1.
Choosing large n we see that C depends on k, and for any N ∈ N, C ≤ k N cannot hold.