On high-order conforming finite element methods for ill-posed Helmholtz problems

Mihai Nechita
(preprint), Finite element (FEM), html, paper, PDEs (numerical)

Abstract


We consider the numerical approximation of the linear ill-posed problem of unique continuation for the Helmholtz equation. We first review the conditional stability of this problem and then discuss high-order conforming finite element methods, with regularization added on the discrete level using gradient jump penalty and Galerkin least squares. The method is shown to converge in terms of the stability of the problem, quantified by the Hölder exponent, and the polynomial degree of approximation. The analysis also takes into account possibly noisy data. Numerical examples illustrating the theory are presented for a Helmholtz version of the classical Hadamard example for two geometric configurations with different stability properties for the continuation problem, including data perturbations of different amplitudes.

Authors

Mihai Nechita
“Tiberiu Popoviciu” Institute of Numerical Analysis, Romanian Academy

Keywords

PDF

Preprint pdf

About this paper

Cite this paper as:
Journal
Publisher Name
DOI

Not available yet.

Print ISSN

Not available yet.

Online ISSN

Not available yet.

Google Scholar Profile

[1] E. Burman, G. Delay, and A. Ern., A hybridized high-order method for unique continuation subject to the helmholtz equation, SIAM J. Numer. Anal., 59(5):2368–2392,2021.[BHL18] E. Burman, P. Hansbo, and M. G. Larson,  Solving ill-posed control problems by stabilized finite element methods: an alternative to Tikhonov regularization, Inverse Problems, 34:035004, 2018.
[2] S. M. Berge and E. Malinnikova, On the three ball theorem for solutions of the Helmholtz equation, Complex Anal. Synerg., 7(2):1–10, 2021.
[3] E. Burman, M. Nechita, and L. Oksanen, Unique continuation for the Helmholtz equation using stabilized finite element methods, J. Math. Pures Appl., 129:1–22, 2019.
[4] I. M. Babuska and S. A. Sauter,  Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers? SIAM J. Numer. Anal., 34(6):2392–2423, 1997.
[5] E. Burman, Stabilized finite element methods for nonsymmetric, noncoercive, and ill-posed problems. Part I: Elliptic equations, SIAM J, Sci. Comput., 35(6):A2752–A2780, 2013.
[6] E. Burman,  Stabilised finite element methods for ill-posed problems with conditional stability. In Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, volume 114 of Lect. Notes Comput. Sci. Eng., pages 93–127. Springer, 2016.
[7] A. Ern and J.-L. Guermond,  Finite elements I: Approximation and interpolation, volume 72 of Texts in Applied Mathematics, Springer Nature, 2021.
[8] T. Hrycak and V. Isakov, Increased stability in the continuation of solutions to the Helmholtz equation, Inverse Problems, 20(3):697–712, 2004.
[9] V. Isakov,  Inverse problems for partial differential equations, volume 127 of Applied Mathematical Sciences, Springer, 3rd edition, 2017.
[10] J. M. Melenk and S. A. Sauter, Wavenumber explicit convergence analysis for Galerkin discretizations of the Helmholtz equation,  SIAM J. Numer. Anal., 49(3):1210–1243, 2011.
[11] M. Nechita, Unique continuation problems and stabilised finite element methods, PhD thesis, University College London, 2020.

Paper (preprint) in HTML form

helmholtz_ho

On high-order conforming finite element methods for ill-posed Helmholtz problems

Mihai Nechita 1 , 2 1 , 2 ^(1,2){ }^{1,2}1,2 1 1 ^(1){ }^{1}1 Commedia, Inria Paris and Laboratoire Jacques-Louis Lions, Sorbonne Université, France 2 2 ^(2){ }^{2}2 Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Cluj-Napoca, Romania

email: mihai.nechita@inria.fr, mihai.nechita@ictp.acad.ro
January 7, 2022

Abstract

We consider the numerical approximation of the linear ill-posed problem of unique continuation for the Helmholtz equation. We first review the conditional stability of this problem and then discuss high-order conforming finite element methods, with regularization added on the discrete level using gradient jump penalty and Galerkin least squares. The method is shown to converge in terms of the stability of the problem, quantified by the Hölder exponent, and the polynomial degree of approximation. The analysis also takes into account possibly noisy data. Numerical examples illustrating the theory are presented for a Helmholtz version of the classical Hadamard example for two geometric configurations with different stability properties for the continuation problem, including data perturbations of different amplitudes.

1 Introduction

Acoustic waves are typically modeled by the Helmholtz equation when their propagation is considered in the frequency domain. For such partial differential equations (PDEs), we are interested in the numerical approximation of linear ill-posed problems appearing in the context of data assimilation, control theory and inverse problems. More precisely, we consider the unique continuation problem, in which no boundary conditions are known, but measurements of the solution are available in a subset of the domain. For an elliptic operator such data can be extended uniquely due to the unique continuation principle: our goal is to numerically approximate the solution in a target domain.
This ill-posed problem has been recently considered in [BNO19], with theoretical results addressing quantitative unique continuation (in norms suitable for numerical analysis) with explicit dependence on the frequency. The stability of the problem is conditional-one must a priori assume a certain bound on the solution-and has Hölder continuity, quantifying the degree of ill-posedness. The sensitivity of the stability estimates with respect to the frequency depends on the geometric configuration and is crucially determined by a convexity condition on the geometry of the target domain relative to the data domain. Roughly speaking, when the target domain is inside the convex hull of the data domain, the stability constant is robust in the frequency, or grows at most linearly, depending on the norm used ( L 2 L 2 L^(2)L^{2}L2 or H 1 H 1 H^(1)H^{1}H1 ). On the other hand, when the solution is estimated outside of the convex hull of the data domain, the stability is very sensitive to the frequency as the constants can grow superpolynomially.
From the computational point of view, [BNO19] introduced and analyzed a low-order finite element method (FEM). In brief, the Helmholtz unique continuation problem-formulated as a PDE-constrained optimization problem-is first discretized using standard piecewise affine elements and then regularization is added using established techniques from the literature on stabilized FEM. Convergence is proven using the conditional stability estimates from the continuous level. This discretize-then-regularize framework represents an alternative to the more standard Tikhonov regularization (see [BHL18] for a comparison) and, in particular, it provides a way of unifying the analysis of the discretization error and the data perturbation error. In this paper we adopt the same framework but consider high-order conforming approximations.
In general, (well-posed) Helmholtz problems are challenging to solve numerically: first, due to the oscillatory nature of the solutions, and secondly, due to the so-called pollution effect [BS97]: the constant bounding the error by the best-approximation error increases with the frequency, preventing the error estimate from being quasi-optimal. High-order methods are typically preferred for such problems since, besides the usual improvement in accuracy, the pollution effect can be reduced or, under some conditions, even eliminated when increasing the polynomial degree of the approximation [MS11]. When it comes to ill-posed problems for elliptic equations, high-order FEM were introduced in [Bur13] and subsequently analyzed in [Bur16] for conditionally stable problems (with a focus on the elliptic Cauchy problem), where it was illustrated numerically that, when dealing with such problems, faster convergence comes with a price: a higher sensitivity to possible perturbations in data. For the Helmholtz equation, a high-order hybridized discontinuous Galerkin method for the unique continuation problem was recently proposed in [BDE21]. Error estimates were proven in terms of the degree of illposedness, the polynomial order and perturbations in data (which get amplified as the order increases). However, a study of standard high-order conforming methods for such ill-posed Helmholtz problems is currently missing from the literature and the aim of this paper is to fill this gap.

2 Ill-posed problem

Figure 1: Domains for the unique continuation problem (1).
Let Ω ⊂ R d , d ∈ { 2 , 3 } Ω ⊂ R d , d ∈ { 2 , 3 } Omega subR^(d),d in{2,3}\Omega \subset \mathbb{R}^{d}, d \in\{2,3\}Ω⊂Rd,d∈{2,3}, be an open, bounded and connected set, and let ω ⊂ Ω ω ⊂ Ω omega sub Omega\omega \subset \Omegaω⊂Ω be a subset (also open, bounded and connected), as in Figure 1. For the Helmholtz equation with wave number k > 0 k > 0 k > 0k>0k>0 (from now on called frequency, for convenience), we consider the unique continuation problem of finding u ∈ H 1 ( Ω ) u ∈ H 1 ( Ω ) u inH^(1)(Omega)u \in H^{1}(\Omega)u∈H1(Ω) such that
{ (1) − Δ u − k 2 u = f in Ω , u = g in ω , (1) − Δ u − k 2 u = f  in  Ω , u = g  in  ω , {[(1)-Delta u-k^(2)u=f" in "Omega","],[u=g" in "omega","]:}\left\{\begin{align*} -\Delta u-k^{2} u=f & \text { in } \Omega, \tag{1}\\ u=g & \text { in } \omega, \end{align*}\right.{(1)−Δu−k2u=f in Ω,u=g in ω,
with source term f ∈ L 2 ( Ω ) f ∈ L 2 ( Ω ) f inL^(2)(Omega)f \in L^{2}(\Omega)f∈L2(Ω) and interior datum g ∈ H 1 ( ω ) g ∈ H 1 ( ω ) g inH^(1)(omega)g \in H^{1}(\omega)g∈H1(ω). Note that no boundary conditions are prescribed, and partial measurements are given in a subset of the domain. If problem (1) has a solution, then its uniqueness is guaranteed by the unique continuation principle for elliptic operators, see for example [Isa17]. We will assume that the data f f fff and g g ggg are compatible, such that the function g g ggg is the restriction on ω ω omega\omegaω of a solution to the Helmholtz equation with source term f f fff.
It is well-known, see for example [Isa17], that problem (1) is ill-posed in the sense of Hadamard: there is no uniform stability with respect to the data f f fff and g g ggg. The reason for this is the fact that the operator mapping a function u u uuu to the corresponding pair ( f , g f , g f,gf, gf,g ) does not have a closed range, and its generalized inverse is, thus, discontinuous. This means that one can find functions in the kernel of the Helmholtz operator that are small in the data set ω ω omega\omegaω, but grow exponentially away from it, as in the example given in Section 4 concerning the Hadamard-type solutions (20). Nonetheless, assuming an additional a priori bound, the solution can be bounded by the data-this is known as conditional stability and in what follows we will present a local version of it.

2.1 Conditional stability and frequency dependence

We denote by B ⊂ Ω B ⊂ Ω B sub OmegaB \subset \OmegaB⊂Ω the target set containing ω ω omega\omegaω such that B ∖ ω ¯ ⊂ Ω B ∖ ω ¯ ⊂ Ω B\\ bar(omega)sub OmegaB \backslash \bar{\omega} \subset \OmegaB∖ω¯⊂Ω, that is B ∖ ω B ∖ ω B\\omegaB \backslash \omegaB∖ω does not touch the boundary of Ω Ω Omega\OmegaΩ. It is well-known, see for example [Isa17], that there exist constants C s t ( k ) > 0 C s t ( k ) > 0 C_(st)(k) > 0C_{s t}(k)>0Cst(k)>0 and α ∈ ( 0 , 1 ] α ∈ ( 0 , 1 ] alpha in(0,1]\alpha \in(0,1]α∈(0,1] such that the following Hölder stability estimate holds
(2) ‖ u ‖ H 1 ( B ) ≤ C s t ( k ) ( ‖ f ‖ L 2 ( Ω ) + ‖ g ‖ H 1 ( ω ) ) α ‖ u ‖ H 1 ( Ω ) 1 − α , (2) ‖ u ‖ H 1 ( B ) ≤ C s t ( k ) ‖ f ‖ L 2 ( Ω ) + ‖ g ‖ H 1 ( ω ) α ‖ u ‖ H 1 ( Ω ) 1 − α , {:(2)||u||_(H^(1)(B)) <= C_(st)(k)(||f||_(L^(2)(Omega))+||g||_(H^(1)(omega)))^(alpha)||u||_(H^(1)(Omega))^(1-alpha)",":}\begin{equation*} \|u\|_{H^{1}(B)} \leq C_{s t}(k)\left(\|f\|_{L^{2}(\Omega)}+\|g\|_{H^{1}(\omega)}\right)^{\alpha}\|u\|_{H^{1}(\Omega)}^{1-\alpha}, \tag{2} \end{equation*}(2)‖u‖H1(B)≤Cst(k)(‖f‖L2(Ω)+‖g‖H1(ω))α‖u‖H1(Ω)1−α,
for any u ∈ H 1 ( Ω ) u ∈ H 1 ( Ω ) u inH^(1)(Omega)u \in H^{1}(\Omega)u∈H1(Ω) satisfying (1), where the stability constant C s t ( k ) C s t ( k ) C_(st)(k)C_{s t}(k)Cst(k) depends on the frequency k k kkk, typically in an implicit way. The exponent α ∈ ( 0 , 1 ] α ∈ ( 0 , 1 ] alpha in(0,1]\alpha \in(0,1]α∈(0,1] encodes the degree of ill-posedness for the continuation problem: note that α = 1 α = 1 alpha=1\alpha=1α=1 would correspond to a well-posed problem with Lipschitz stability, and as α < 1 α < 1 alpha < 1\alpha<1α<1 gets smaller the Hölder stability deteriorates. Both C s t C s t C_(st)C_{s t}Cst and α α alpha\alphaα depend on the geometric configuration in a nontrivial way.
Figure 2: Different frequency dependence for the stability constant. Data set ω ω omega\omegaω (dark grey) and target set B B BBB (light grey). Ω Ω Omega\OmegaΩ is the large box.
How the stability constant C s t ( k ) C s t ( k ) C_(st)(k)C_{s t}(k)Cst(k) depends on the frequency k k kkk is an important aspect for the stability of this ill-posed Helmholtz problem. For example, when there is a straight line that intersects B B BBB but not ω ¯ ω ¯ bar(omega)\bar{\omega}ω¯, as in Figure 2a, it was proven in [BNO19, Example 4] that for any N ∈ N , C s t ( k ) ≤ k N N ∈ N , C s t ( k ) ≤ k N N inN,C_(st)(k) <= k^(N)N \in \mathbb{N}, C_{s t}(k) \leq k^{N}N∈N,Cst(k)≤kN cannot hold uniformly in k k kkk-in other words, the stability constant in (2)
grows superpolynomially in the frequency. Estimates with such a fast growth might be relevant in practice only for low frequencies, and it is thus of interest to consider obtaining estimates with a more favorable frequency dependence. Bounds that have a good behavior with respect to the frequency can be obtained under a convexity condition of the target domain B B BBB relative to the data set ω ω omega\omegaω-essentially that B B BBB is included in the convex hull of ω ω omega\omegaω, as for example in Figure 2b. Such a condition was first considered in [HI04], where is was shown that the stability of the solution in the L 2 L 2 L^(2)L^{2}L2-norm can actually improve in a certain sense as the frequency increases. In this vein, similar results that make use of a convexity condition were proven in [BNO19] for a particular geometric setting prototypical for continuation inside the convex hull of ω ω omega\omegaω. It was shown in [BNO19, Corollary 2] that the bound (2) is robust with respect to the frequency, i.e., it holds with the stability constant C s t C s t C_(st)C_{s t}Cst independent of k k kkk. This robust bound is also valid for other geometric configurations, as for example the one in Figure 2b.
In order to use such conditional stability estimates for the numerical analysis of a finite element method, the norms for measuring the data must be weaken: from L 2 ( Ω ) L 2 ( Ω ) L^(2)(Omega)L^{2}(\Omega)L2(Ω) and H 1 ( ω ) H 1 ( ω ) H^(1)(omega)H^{1}(\omega)H1(ω) to H − 1 ( Ω ) H − 1 ( Ω ) H^(-1)(Omega)H^{-1}(\Omega)H−1(Ω) and L 2 ( ω ) L 2 ( ω ) L^(2)(omega)L^{2}(\omega)L2(ω), respectively. A result of this kind was proven in [BNO19, Corollary 3 and Lemma 2] under a geometric convexity condition. More precisely, there exist constants C > 0 C > 0 C > 0C>0C>0 and α ∈ ( 0 , 1 ] α ∈ ( 0 , 1 ] alpha in(0,1]\alpha \in(0,1]α∈(0,1] such that
(3) ‖ ∇ u ‖ L 2 ( B ) + k ‖ u ‖ L 2 ( B ) ≤ C k ( ‖ f ‖ H − 1 ( Ω ) + ‖ g ‖ L 2 ( ω ) ) α ( ‖ f ‖ H − 1 ( Ω ) + ‖ u ‖ L 2 ( Ω ) ) 1 − α , (3) ‖ ∇ u ‖ L 2 ( B ) + k ‖ u ‖ L 2 ( B ) ≤ C k ‖ f ‖ H − 1 ( Ω ) + ‖ g ‖ L 2 ( ω ) α ‖ f ‖ H − 1 ( Ω ) + ‖ u ‖ L 2 ( Ω ) 1 − α , {:(3)||grad u||_(L^(2)(B))+k||u||_(L^(2)(B)) <= Ck(||f||_(H^(-1)(Omega))+||g||_(L^(2)(omega)))^(alpha)(||f||_(H^(-1)(Omega))+||u||_(L^(2)(Omega)))^(1-alpha)",":}\begin{equation*} \|\nabla u\|_{L^{2}(B)}+k\|u\|_{L^{2}(B)} \leq C k\left(\|f\|_{H^{-1}(\Omega)}+\|g\|_{L^{2}(\omega)}\right)^{\alpha}\left(\|f\|_{H^{-1}(\Omega)}+\|u\|_{L^{2}(\Omega)}\right)^{1-\alpha}, \tag{3} \end{equation*}(3)‖∇u‖L2(B)+k‖u‖L2(B)≤Ck(‖f‖H−1(Ω)+‖g‖L2(ω))α(‖f‖H−1(Ω)+‖u‖L2(Ω))1−α,
for any u ∈ H 1 ( Ω ) u ∈ H 1 ( Ω ) u inH^(1)(Omega)u \in H^{1}(\Omega)u∈H1(Ω) satisfying (1). Note that the bound is robust in the L 2 L 2 L^(2)L^{2}L2-norm, while the frequency dependence is linear for the H 1 H 1 H^(1)H^{1}H1-seminorm. Apart from the good dependence on the frequency, this kind of estimate is particularly suitable for numerical analysis since it can be directly applied to the finite element error equation.
As discussed above, when the domains do not satisfy the convexity condition, the stability constant can be sensitive to the frequency. To mark the difference, in general we write the stability estimate as follows: there exist constants C s t ( k ) > 0 C s t ( k ) > 0 C_(st)(k) > 0C_{s t}(k)>0Cst(k)>0 and α ∈ ( 0 , 1 ] α ∈ ( 0 , 1 ] alpha in(0,1]\alpha \in(0,1]α∈(0,1] such that
(4) ‖ u ‖ H 1 ( B ) ≤ C s t ( k ) ( ‖ f ‖ H − 1 ( Ω ) + ‖ g ‖ L 2 ( ω ) ) α ( ‖ f ‖ H − 1 ( Ω ) + ‖ u ‖ L 2 ( Ω ) ) 1 − α . (4) ‖ u ‖ H 1 ( B ) ≤ C s t ( k ) ‖ f ‖ H − 1 ( Ω ) + ‖ g ‖ L 2 ( ω ) α ‖ f ‖ H − 1 ( Ω ) + ‖ u ‖ L 2 ( Ω ) 1 − α . {:(4)||u||_(H^(1)(B)) <= C_(st)(k)(||f||_(H^(-1)(Omega))+||g||_(L^(2)(omega)))^(alpha)(||f||_(H^(-1)(Omega))+||u||_(L^(2)(Omega)))^(1-alpha).:}\begin{equation*} \|u\|_{H^{1}(B)} \leq C_{s t}(k)\left(\|f\|_{H^{-1}(\Omega)}+\|g\|_{L^{2}(\omega)}\right)^{\alpha}\left(\|f\|_{H^{-1}(\Omega)}+\|u\|_{L^{2}(\Omega)}\right)^{1-\alpha} . \tag{4} \end{equation*}(4)‖u‖H1(B)≤Cst(k)(‖f‖H−1(Ω)+‖g‖L2(ω))α(‖f‖H−1(Ω)+‖u‖L2(Ω))1−α.
Let us also remark that for three-ball inequalities (where ω , B , Ω ω , B , Ω omega,B,Omega\omega, B, \Omegaω,B,Ω are concentric balls) it was recently shown in [BM21] that, when using the maximum norm, C s t ( k ) C s t ( k ) C_(st)(k)C_{s t}(k)Cst(k) grows exponentially in k k kkk and this dependence is optimal.

3 Discretization, regularization and error analysis

3.1 Finite element preliminaries

From now on we will consider for simplicity that Ω ⊂ R d , d ∈ { 2 , 3 } Ω ⊂ R d , d ∈ { 2 , 3 } Omega subR^(d),d in{2,3}\Omega \subset \mathbb{R}^{d}, d \in\{2,3\}Ω⊂Rd,d∈{2,3}, is a convex polygonal/polyhedral domain. Let { T h } h > 0 T h h > 0 {T_(h)}_(h > 0)\left\{\mathcal{T}_{h}\right\}_{h>0}{Th}h>0 be a quasi-uniform family of shape-regular meshes covering Ω Ω Omega\OmegaΩ, each triangulation T h T h T_(h)\mathcal{T}_{h}Th containing elements K K KKK with maximal diameter h h hhh. Let p p ppp denote the polynomial degree of approximation and let P p ( K ) P p ( K ) P_(p)(K)\mathbb{P}_{p}(K)Pp(K) be the set of polynomials of degree at most p p ppp defined on an element K K KKK. We consider the conforming finite element space
V h p := { v h ∈ H 1 ( Ω ) : v h | K ∈ P p ( K ) , ∀ K ∈ T h } , V h p := v h ∈ H 1 ( Ω ) : v h K ∈ P p ( K ) , ∀ K ∈ T h , V_(h)^(p):={v_(h)inH^(1)(Omega):v_(h)|_(K)inP_(p)(K),AA K inT_(h)},V_{h}^{p}:=\left\{v_{h} \in H^{1}(\Omega):\left.v_{h}\right|_{K} \in \mathbb{P}_{p}(K), \forall K \in \mathcal{T}_{h}\right\},Vhp:={vh∈H1(Ω):vh|K∈Pp(K),∀K∈Th},
with continuous, piecewise polynomial functions defined on the mesh T h T h T_(h)\mathcal{T}_{h}Th, and its subspace with homogeneous boundary conditions
W h p := V h p ∩ H 0 1 ( Ω ) . W h p := V h p ∩ H 0 1 ( Ω ) . W_(h)^(p):=V_(h)^(p)nnH_(0)^(1)(Omega).W_{h}^{p}:=V_{h}^{p} \cap H_{0}^{1}(\Omega) .Whp:=Vhp∩H01(Ω).
For a given mesh, we collect all the interior edges/faces of the elements in the set F i F i F_(i)\mathcal{F}_{i}Fi and denote the jump of a quantity across such an edge/face F F FFF by [ [ ⋅ ] ] F [ [ ⋅ ] ] F [[*]]_(F)\llbracket \cdot \rrbracket_{F}[[⋅]]F, omitting the subscript whenever there is no confusion. For a discrete function v h ∈ V h p v h ∈ V h p v_(h)inV_(h)^(p)v_{h} \in V_{h}^{p}vh∈Vhp, we denote the jump of the normal gradient across F ∈ F i F ∈ F i F inF_(i)F \in \mathcal{F}_{i}F∈Fi by
[ [ ∇ v h ⋅ n ] ] F := ∇ v h ⋅ n 1 | K 1 + ∇ v h ⋅ n 2 | K 2 , [ [ ∇ v h ⋅ n ] ] F := ∇ v h ⋅ n 1 K 1 + ∇ v h ⋅ n 2 K 2 , [[gradv_(h)*n]]_(F):= gradv_(h)*n_(1)|_(K_(1))+ gradv_(h)*n_(2)|_(K_(2)),\llbracket \nabla v_{h} \cdot n \rrbracket_{F}:=\left.\nabla v_{h} \cdot n_{1}\right|_{K_{1}}+\left.\nabla v_{h} \cdot n_{2}\right|_{K_{2}},[[∇vh⋅n]]F:=∇vh⋅n1|K1+∇vh⋅n2|K2,
with K 1 , K 2 ∈ T h K 1 , K 2 ∈ T h K_(1),K_(2)inT_(h)K_{1}, K_{2} \in \mathcal{T}_{h}K1,K2∈Th being two elements such that K 1 ∩ K 2 = F K 1 ∩ K 2 = F K_(1)nnK_(2)=FK_{1} \cap K_{2}=FK1∩K2=F, and n 1 , n 2 n 1 , n 2 n_(1),n_(2)n_{1}, n_{2}n1,n2 the outward pointing normals of K 1 , K 2 K 1 , K 2 K_(1),K_(2)K_{1}, K_{2}K1,K2, respectively.
We now recall some classical finite element results, see for example [EG21]. We will use the inverse inequality
(5) ‖ ∇ v h ‖ L 2 ( K ) ≤ C h − 1 ‖ v h ‖ L 2 ( K ) , ∀ v h ∈ P p ( K ) (5) ∇ v h L 2 ( K ) ≤ C h − 1 v h L 2 ( K ) , ∀ v h ∈ P p ( K ) {:(5)||gradv_(h)||_(L^(2)(K)) <= Ch^(-1)||v_(h)||_(L^(2)(K))","quad AAv_(h)inP_(p)(K):}\begin{equation*} \left\|\nabla v_{h}\right\|_{L^{2}(K)} \leq C h^{-1}\left\|v_{h}\right\|_{L^{2}(K)}, \quad \forall v_{h} \in \mathbb{P}_{p}(K) \tag{5} \end{equation*}(5)‖∇vh‖L2(K)≤Ch−1‖vh‖L2(K),∀vh∈Pp(K)
and the trace inequality
(6) ‖ v ‖ L 2 ( ∂ K ) ≤ C ( h − 1 2 ‖ v ‖ L 2 ( K ) + h 1 2 ‖ ∇ v ‖ L 2 ( K ) ) , ∀ v ∈ H 1 ( K ) . (6) ‖ v ‖ L 2 ( ∂ K ) ≤ C h − 1 2 ‖ v ‖ L 2 ( K ) + h 1 2 ‖ ∇ v ‖ L 2 ( K ) , ∀ v ∈ H 1 ( K ) . {:(6)||v||_(L^(2)(del K)) <= C(h^(-(1)/(2))||v||_(L^(2)(K))+h^((1)/(2))||grad v||_(L^(2)(K)))","quad AA v inH^(1)(K).:}\begin{equation*} \|v\|_{L^{2}(\partial K)} \leq C\left(h^{-\frac{1}{2}}\|v\|_{L^{2}(K)}+h^{\frac{1}{2}}\|\nabla v\|_{L^{2}(K)}\right), \quad \forall v \in H^{1}(K) . \tag{6} \end{equation*}(6)‖v‖L2(∂K)≤C(h−12‖v‖L2(K)+h12‖∇v‖L2(K)),∀v∈H1(K).
We denote by π h : L 2 ( Ω ) → V h p π h : L 2 ( Ω ) → V h p pi_(h):L^(2)(Omega)rarrV_(h)^(p)\pi_{h}: L^{2}(\Omega) \rightarrow V_{h}^{p}πh:L2(Ω)→Vhp the L 2 L 2 L^(2)L^{2}L2-projection, that satisfies the orthogonality
( u − π h u , v h ) L 2 ( Ω ) = 0 , ∀ u ∈ L 2 ( Ω ) , ∀ v h ∈ V h p . u − π h u , v h L 2 ( Ω ) = 0 , ∀ u ∈ L 2 ( Ω ) , ∀ v h ∈ V h p . (u-pi_(h)u,v_(h))_(L^(2)(Omega))=0,quad AA u inL^(2)(Omega),AAv_(h)inV_(h)^(p).\left(u-\pi_{h} u, v_{h}\right)_{L^{2}(\Omega)}=0, \quad \forall u \in L^{2}(\Omega), \forall v_{h} \in V_{h}^{p} .(u−πhu,vh)L2(Ω)=0,∀u∈L2(Ω),∀vh∈Vhp.
The L 2 L 2 L^(2)L^{2}L2-projection is stable in the L 2 L 2 L^(2)L^{2}L2-norm
‖ π h u ‖ L 2 ( Ω ) ≤ ‖ u ‖ L 2 ( Ω ) , ∀ u ∈ L 2 ( Ω ) π h u L 2 ( Ω ) ≤ ‖ u ‖ L 2 ( Ω ) , ∀ u ∈ L 2 ( Ω ) ||pi_(h)u||_(L^(2)(Omega)) <= ||u||_(L^(2)(Omega)),quad AA u inL^(2)(Omega)\left\|\pi_{h} u\right\|_{L^{2}(\Omega)} \leq\|u\|_{L^{2}(\Omega)}, \quad \forall u \in L^{2}(\Omega)‖πhu‖L2(Ω)≤‖u‖L2(Ω),∀u∈L2(Ω)
and, since the family { T h } T h {T_(h)}\left\{\mathcal{T}_{h}\right\}{Th} is quasi-uniform, it is also stable in the H 1 H 1 H^(1)H^{1}H1-norm
‖ π h u ‖ H 1 ( Ω ) ≤ C ‖ u ‖ H 1 ( Ω ) , ∀ u ∈ H 1 ( Ω ) π h u H 1 ( Ω ) ≤ C ‖ u ‖ H 1 ( Ω ) , ∀ u ∈ H 1 ( Ω ) ||pi_(h)u||_(H^(1)(Omega)) <= C||u||_(H^(1)(Omega)),quad AA u inH^(1)(Omega)\left\|\pi_{h} u\right\|_{H^{1}(\Omega)} \leq C\|u\|_{H^{1}(\Omega)}, \quad \forall u \in H^{1}(\Omega)‖πhu‖H1(Ω)≤C‖u‖H1(Ω),∀u∈H1(Ω)
We also use the Scott-Zhang interpolant π s z : H 1 ( Ω ) → V h p π s z : H 1 ( Ω ) → V h p pi_(sz):H^(1)(Omega)rarrV_(h)^(p)\pi_{s z}: H^{1}(\Omega) \rightarrow V_{h}^{p}πsz:H1(Ω)→Vhp which preserves vanishing Dirichlet boundary conditions. This interpolant is stable in both the L 2 L 2 L^(2)L^{2}L2 - and the H 1 H 1 H^(1)H^{1}H1-norm and enjoys the same approximation error estimate as the L 2 L 2 L^(2)L^{2}L2-projection
(7) ‖ u − i h u ‖ L 2 ( Ω ) + h ‖ ∇ ( u − i h u ) ‖ L 2 ( Ω ) ≤ C h m | u | H m ( Ω ) , ∀ u ∈ H m ( Ω ) , i h = π h , π s z . (7) u − i h u L 2 ( Ω ) + h ∇ u − i h u L 2 ( Ω ) ≤ C h m | u | H m ( Ω ) , ∀ u ∈ H m ( Ω ) , i h = π h , π s z . {:(7)||u-i_(h)u||_(L^(2)(Omega))+h||grad(u-i_(h)u)||_(L^(2)(Omega)) <= Ch^(m)|u|_(H^(m)(Omega))","quad AA u inH^(m)(Omega)","i_(h)=pi_(h)","pi_(sz).:}\begin{equation*} \left\|u-i_{h} u\right\|_{L^{2}(\Omega)}+h\left\|\nabla\left(u-i_{h} u\right)\right\|_{L^{2}(\Omega)} \leq C h^{m}|u|_{H^{m}(\Omega)}, \quad \forall u \in H^{m}(\Omega), i_{h}=\pi_{h}, \pi_{s z} . \tag{7} \end{equation*}(7)‖u−ihu‖L2(Ω)+h‖∇(u−ihu)‖L2(Ω)≤Chm|u|Hm(Ω),∀u∈Hm(Ω),ih=πh,πsz.

3.2 Discretization

We denote the standard inner product (with its induced norm) on a set S S SSS by
( v h , w h ) S := ∫ S v h w h d x v h , w h S := ∫ S   v h w h d x (v_(h),w_(h))_(S):=int_(S)v_(h)w_(h)dx\left(v_{h}, w_{h}\right)_{S}:=\int_{S} v_{h} w_{h} \mathrm{~d} x(vh,wh)S:=∫Svhwh dx
and define on V h p × W h p V h p × W h p V_(h)^(p)xxW_(h)^(p)V_{h}^{p} \times W_{h}^{p}Vhp×Whp the standard bilinear form corresponding to the Helmholtz equation
a ( u h , w h ) := ( ∇ u h , ∇ w h ) Ω − k 2 ( u h , w h ) Ω . a u h , w h := ∇ u h , ∇ w h Ω − k 2 u h , w h Ω . a(u_(h),w_(h)):=(gradu_(h),gradw_(h))_(Omega)-k^(2)(u_(h),w_(h))_(Omega).a\left(u_{h}, w_{h}\right):=\left(\nabla u_{h}, \nabla w_{h}\right)_{\Omega}-k^{2}\left(u_{h}, w_{h}\right)_{\Omega} .a(uh,wh):=(∇uh,∇wh)Ω−k2(uh,wh)Ω.
The discrete weak form of the continuation problem (1) reads as follows: find u h ∈ V h p u h ∈ V h p u_(h)inV_(h)^(p)u_{h} \in V_{h}^{p}uh∈Vhp such that
(8) a ( u h , w h ) = ( f , w h ) Ω , ∀ w h ∈ W h p and u h = g in ω . (8) a u h , w h = f , w h Ω , ∀ w h ∈ W h p  and  u h = g  in  ω . {:(8)a(u_(h),w_(h))=(f,w_(h))_(Omega)","quad AAw_(h)inW_(h)^(p)quad" and "quadu_(h)=g" in "omega.:}\begin{equation*} a\left(u_{h}, w_{h}\right)=\left(f, w_{h}\right)_{\Omega}, \quad \forall w_{h} \in W_{h}^{p} \quad \text { and } \quad u_{h}=g \text { in } \omega . \tag{8} \end{equation*}(8)a(uh,wh)=(f,wh)Ω,∀wh∈Whp and uh=g in ω.
To take into account the presence of noise in the measurements, we will consider the perturbed datum
g ~ = g + δ g , δ g ∈ L 2 ( ω ) g ~ = g + δ g , δ g ∈ L 2 ( ω ) tilde(g)=g+delta g,quad delta g inL^(2)(omega)\tilde{g}=g+\delta g, \quad \delta g \in L^{2}(\omega)g~=g+δg,δg∈L2(ω)
As a starting point, the continuation problem can be formulated on the discrete level as a constrained minimization:
min u h ∈ V h p 1 2 ‖ u h − g ~ ‖ ω 2 subject to a ( u h , w h ) = ( f , w h ) Ω , ∀ w h ∈ W h p min u h ∈ V h p   1 2 u h − g ~ ω 2  subject to  a u h , w h = f , w h Ω , ∀ w h ∈ W h p min_(u_(h)inV_(h)^(p))(1)/(2)||u_(h)-( tilde(g))||_(omega)^(2)quad" subject to "quad a(u_(h),w_(h))=(f,w_(h))_(Omega),quad AAw_(h)inW_(h)^(p)\min _{u_{h} \in V_{h}^{p}} \frac{1}{2}\left\|u_{h}-\tilde{g}\right\|_{\omega}^{2} \quad \text { subject to } \quad a\left(u_{h}, w_{h}\right)=\left(f, w_{h}\right)_{\Omega}, \quad \forall w_{h} \in W_{h}^{p}minuh∈Vhp12‖uh−g~‖ω2 subject to a(uh,wh)=(f,wh)Ω,∀wh∈Whp
where the L 2 L 2 L^(2)L^{2}L2-distance to the interior datum is minimized under the constraint of the weak form of the differential equation. Introducing a test function as a Lagrange multiplier, one can consider the discrete Lagrangian functional
L h 0 ( u h , z h ) := 1 2 ‖ u h − g ~ ‖ ω 2 + a ( u h , z h ) − ( f , z h ) Ω , ( u h , z h ) ∈ V h p × W h p . L h 0 u h , z h := 1 2 u h − g ~ ω 2 + a u h , z h − f , z h Ω , u h , z h ∈ V h p × W h p . L_(h)^(0)(u_(h),z_(h)):=(1)/(2)||u_(h)-( tilde(g))||_(omega)^(2)+a(u_(h),z_(h))-(f,z_(h))_(Omega),quad(u_(h),z_(h))inV_(h)^(p)xxW_(h)^(p).L_{h}^{0}\left(u_{h}, z_{h}\right):=\frac{1}{2}\left\|u_{h}-\tilde{g}\right\|_{\omega}^{2}+a\left(u_{h}, z_{h}\right)-\left(f, z_{h}\right)_{\Omega}, \quad\left(u_{h}, z_{h}\right) \in V_{h}^{p} \times W_{h}^{p} .Lh0(uh,zh):=12‖uh−g~‖ω2+a(uh,zh)−(f,zh)Ω,(uh,zh)∈Vhp×Whp.
The solution to the above optimization problem is given by the saddle-points of L h 0 L h 0 L_(h)^(0)L_{h}^{0}Lh0, satisfying the Euler-Lagrange equations. However, the ill-posedness of problem (8) has not been alleviated in any way by this reformulation, and the optimality conditions lead to a linear system that might be singular.

3.3 Regularization

To mitigate the ill-posedness, one must regularize the problem by adding some penalty terms that provide well-posedness. The approach that we take herein has been introduced in [Bur13] and is a discretize-then-regularize one. We add the regularization terms on the discrete level, using techniques from stabilized FEM for numerically unstable well-posed problems: gradient jump penalty and Galerkin least squares. Apart from providing the needed stability, the regularizing terms should ideally perturb the problem in a weakly consistent way by converging to zero at optimal rates for smooth enough solutions.
For an approximation u h ∈ V h p u h ∈ V h p u_(h)inV_(h)^(p)u_{h} \in V_{h}^{p}uh∈Vhp and a Lagrange multiplier z h ∈ W h p z h ∈ W h p z_(h)inW_(h)^(p)z_{h} \in W_{h}^{p}zh∈Whp, we consider the regularized Lagrangian
(9) L h ( u h , z h ) := 1 2 ‖ u h − g ~ ‖ L 2 ( ω ) 2 + 1 2 s ( u h − u , u h − u ) − 1 2 s ∗ ( z h , z h ) + a h ( u h , z h ) − ( f , z h ) Ω , (9) L h u h , z h := 1 2 u h − g ~ L 2 ( ω ) 2 + 1 2 s u h − u , u h − u − 1 2 s ∗ z h , z h + a h u h , z h − f , z h Ω , {:(9)L_(h)(u_(h),z_(h)):=(1)/(2)||u_(h)-( tilde(g))||_(L^(2)(omega))^(2)+(1)/(2)s(u_(h)-u,u_(h)-u)-(1)/(2)s^(**)(z_(h),z_(h))+a_(h)(u_(h),z_(h))-(f,z_(h))_(Omega)",":}\begin{equation*} L_{h}\left(u_{h}, z_{h}\right):=\frac{1}{2}\left\|u_{h}-\tilde{g}\right\|_{L^{2}(\omega)}^{2}+\frac{1}{2} s\left(u_{h}-u, u_{h}-u\right)-\frac{1}{2} s^{*}\left(z_{h}, z_{h}\right)+a_{h}\left(u_{h}, z_{h}\right)-\left(f, z_{h}\right)_{\Omega}, \tag{9} \end{equation*}(9)Lh(uh,zh):=12‖uh−g~‖L2(ω)2+12s(uh−u,uh−u)−12s∗(zh,zh)+ah(uh,zh)−(f,zh)Ω,
with s s sss and s ∗ s ∗ s^(**)s^{*}s∗ being the stabilization/regularization terms, to be specified later. The saddlepoints ( u h , z h ) ∈ V h p × W h p u h , z h ∈ V h p × W h p (u_(h),z_(h))inV_(h)^(p)xxW_(h)^(p)\left(u_{h}, z_{h}\right) \in V_{h}^{p} \times W_{h}^{p}(uh,zh)∈Vhp×Whp of L h L h L_(h)L_{h}Lh satisfy the optimality conditions: for any ( v h , w h ) ∈ V h p × W h p v h , w h ∈ V h p × W h p (v_(h),w_(h))inV_(h)^(p)xxW_(h)^(p)\left(v_{h}, w_{h}\right) \in V_{h}^{p} \times W_{h}^{p}(vh,wh)∈Vhp×Whp,
{ 0 = ∂ u h L h v h = ( u h − g ~ , v h ) ω + s ( u h − u , v h ) + a ( v h , z h ) 0 = ∂ z h L h w h = − s ∗ ( z h , w h ) + a ( u h , w h ) − ( f , w h ) Ω 0 = ∂ u h L h v h = u h − g ~ , v h ω + s u h − u , v h + a v h , z h 0 = ∂ z h L h w h = − s ∗ z h , w h + a u h , w h − f , w h Ω {[0=del_(u_(h))L_(h)v_(h)=(u_(h)-( tilde(g)),v_(h))_(omega)+s(u_(h)-u,v_(h))+a(v_(h),z_(h))],[0=del_(z_(h))L_(h)w_(h)=-s^(**)(z_(h),w_(h))+a(u_(h),w_(h))-(f,w_(h))_(Omega)]:}\left\{\begin{array}{l} 0=\partial_{u_{h}} L_{h} v_{h}=\left(u_{h}-\tilde{g}, v_{h}\right)_{\omega}+s\left(u_{h}-u, v_{h}\right)+a\left(v_{h}, z_{h}\right) \\ 0=\partial_{z_{h}} L_{h} w_{h}=-s^{*}\left(z_{h}, w_{h}\right)+a\left(u_{h}, w_{h}\right)-\left(f, w_{h}\right)_{\Omega} \end{array}\right.{0=∂uhLhvh=(uh−g~,vh)ω+s(uh−u,vh)+a(vh,zh)0=∂zhLhwh=−s∗(zh,wh)+a(uh,wh)−(f,wh)Ω
The stabilized finite element method is thus looking for ( u h , z h ) ∈ V h p × W h p u h , z h ∈ V h p × W h p (u_(h),z_(h))inV_(h)^(p)xxW_(h)^(p)\left(u_{h}, z_{h}\right) \in V_{h}^{p} \times W_{h}^{p}(uh,zh)∈Vhp×Whp such that
(10) { ( u h , v h ) ω + s ( u h , v h ) + a ( v h , z h ) = ( g ~ , v h ) ω + s ( u , v h ) a ( u h , w h ) − s ∗ ( z h , w h ) = ( f , w h ) Ω , ∀ ( v h , w h ) ∈ V h p × W h p . (10) u h , v h ω + s u h , v h + a v h , z h = g ~ , v h ω + s u , v h a u h , w h − s ∗ z h , w h = f , w h Ω , ∀ v h , w h ∈ V h p × W h p . {:(10){[(u_(h),v_(h))_(omega)+s(u_(h),v_(h))+a(v_(h),z_(h)),=(( tilde(g)),v_(h))_(omega)+s(u,v_(h))],[a(u_(h),w_(h))-s^(**)(z_(h),w_(h)),=(f,w_(h))_(Omega)]quad,quad AA(v_(h),w_(h))inV_(h)^(p)xxW_(h)^(p).:}:}\left\{\begin{array}{rl} \left(u_{h}, v_{h}\right)_{\omega}+s\left(u_{h}, v_{h}\right)+a\left(v_{h}, z_{h}\right) & =\left(\tilde{g}, v_{h}\right)_{\omega}+s\left(u, v_{h}\right) \tag{10}\\ a\left(u_{h}, w_{h}\right)-s^{*}\left(z_{h}, w_{h}\right) & =\left(f, w_{h}\right)_{\Omega} \end{array} \quad, \quad \forall\left(v_{h}, w_{h}\right) \in V_{h}^{p} \times W_{h}^{p} .\right.(10){(uh,vh)ω+s(uh,vh)+a(vh,zh)=(g~,vh)ω+s(u,vh)a(uh,wh)−s∗(zh,wh)=(f,wh)Ω,∀(vh,wh)∈Vhp×Whp.
The linear system associated to this problem can be represented in block matrix form as
[ M ω + S A T A − S ∗ ] [ u h z h ] = [ g ~ + s u f ] . M ω + S A T A − S ∗ u h z h = g ~ + s u f . [[M_(omega)+S,A^(T)],[A,-S^(**)]][(u_(h))/(z_(h))]=[[ tilde(g)+s_(u)],[f]].\left[\begin{array}{c|c} M_{\omega}+S & A^{T} \\ \hline A & -S^{*} \end{array}\right]\left[\frac{u_{h}}{z_{h}}\right]=\left[\begin{array}{c} \tilde{g}+s_{u} \\ \hline f \end{array}\right] .[Mω+SATA−S∗][uhzh]=[g~+suf].
To address the well-posedness of the discrete problem we must prove that it satisfies an inf-sup condition. For this, let us introduce the bilinear form
A h [ ( u h , z h ) , ( v h , w h ) ] := ( u h , v h ) ω + s ( u h , v h ) + a h ( v h , z h ) − s ∗ ( z h , w h ) + a h ( u h , w h ) , A h u h , z h , v h , w h := u h , v h ω + s u h , v h + a h v h , z h − s ∗ z h , w h + a h u h , w h , A_(h)[(u_(h),z_(h)),(v_(h),w_(h))]:=(u_(h),v_(h))_(omega)+s(u_(h),v_(h))+a_(h)(v_(h),z_(h))-s^(**)(z_(h),w_(h))+a_(h)(u_(h),w_(h)),A_{h}\left[\left(u_{h}, z_{h}\right),\left(v_{h}, w_{h}\right)\right]:=\left(u_{h}, v_{h}\right)_{\omega}+s\left(u_{h}, v_{h}\right)+a_{h}\left(v_{h}, z_{h}\right)-s^{*}\left(z_{h}, w_{h}\right)+a_{h}\left(u_{h}, w_{h}\right),Ah[(uh,zh),(vh,wh)]:=(uh,vh)ω+s(uh,vh)+ah(vh,zh)−s∗(zh,wh)+ah(uh,wh),
for which the system (10) reads as follows: find ( u h , z h ) ∈ V h p × W h p u h , z h ∈ V h p × W h p (u_(h),z_(h))inV_(h)^(p)xxW_(h)^(p)\left(u_{h}, z_{h}\right) \in V_{h}^{p} \times W_{h}^{p}(uh,zh)∈Vhp×Whp such that
A h [ ( u h , z h ) , ( v h , w h ) ] = ( g ~ , v h ) ω + s ( u , v h ) + ( f , w h ) Ω , ∀ ( v h , w h ) ∈ V h p × W h p A h u h , z h , v h , w h = g ~ , v h ω + s u , v h + f , w h Ω , ∀ v h , w h ∈ V h p × W h p A_(h)[(u_(h),z_(h)),(v_(h),w_(h))]=(( tilde(g)),v_(h))_(omega)+s(u,v_(h))+(f,w_(h))_(Omega),quad AA(v_(h),w_(h))inV_(h)^(p)xxW_(h)^(p)A_{h}\left[\left(u_{h}, z_{h}\right),\left(v_{h}, w_{h}\right)\right]=\left(\tilde{g}, v_{h}\right)_{\omega}+s\left(u, v_{h}\right)+\left(f, w_{h}\right)_{\Omega}, \quad \forall\left(v_{h}, w_{h}\right) \in V_{h}^{p} \times W_{h}^{p}Ah[(uh,zh),(vh,wh)]=(g~,vh)ω+s(u,vh)+(f,wh)Ω,∀(vh,wh)∈Vhp×Whp
Notice that a simple choice of the test functions gives that
A h [ ( u h , z h ) , ( u h , − z h ) ] = ‖ u h ‖ ω 2 + s ( u h , u h ) + s ∗ ( z h , z h ) . A h u h , z h , u h , − z h = u h ω 2 + s u h , u h + s ∗ z h , z h . A_(h)[(u_(h),z_(h)),(u_(h),-z_(h))]=||u_(h)||_(omega)^(2)+s(u_(h),u_(h))+s^(**)(z_(h),z_(h)).A_{h}\left[\left(u_{h}, z_{h}\right),\left(u_{h},-z_{h}\right)\right]=\left\|u_{h}\right\|_{\omega}^{2}+s\left(u_{h}, u_{h}\right)+s^{*}\left(z_{h}, z_{h}\right) .Ah[(uh,zh),(uh,−zh)]=‖uh‖ω2+s(uh,uh)+s∗(zh,zh).
To prove an inf-sup condition, it is enough to define the stabilizing bilinear forms s s sss and s ∗ s ∗ s^(**)s^{*}s∗ such that
(11) ‖ ( v h , w h ) ‖ s := ( ‖ v h ‖ ω 2 + s ( v h , v h ) + s ∗ ( w h , w h ) ) 1 2 (11) v h , w h s := v h ω 2 + s v h , v h + s ∗ w h , w h 1 2 {:(11)||(v_(h),w_(h))||_(s):=(||v_(h)||_(omega)^(2)+s(v_(h),v_(h))+s^(**)(w_(h),w_(h)))^((1)/(2)):}\begin{equation*} \left\|\left(v_{h}, w_{h}\right)\right\|_{s}:=\left(\left\|v_{h}\right\|_{\omega}^{2}+s\left(v_{h}, v_{h}\right)+s^{*}\left(w_{h}, w_{h}\right)\right)^{\frac{1}{2}} \tag{11} \end{equation*}(11)‖(vh,wh)‖s:=(‖vh‖ω2+s(vh,vh)+s∗(wh,wh))12
represents a norm on V h p × W h p V h p × W h p V_(h)^(p)xxW_(h)^(p)V_{h}^{p} \times W_{h}^{p}Vhp×Whp. Indeed, we have that
A h [ ( u h , z h ) , ( u h , − z h ) ] = ‖ ( u h , z h ) ‖ s 2 A h u h , z h , u h , − z h = u h , z h s 2 A_(h)[(u_(h),z_(h)),(u_(h),-z_(h))]=||(u_(h),z_(h))||_(s)^(2)A_{h}\left[\left(u_{h}, z_{h}\right),\left(u_{h},-z_{h}\right)\right]=\left\|\left(u_{h}, z_{h}\right)\right\|_{s}^{2}Ah[(uh,zh),(uh,−zh)]=‖(uh,zh)‖s2
and hence we obtain the following inf-sup condition
(12) sup ( v h , w h ) ∈ V h p × W h p A h [ ( u h , z h ) , ( v h , w h ) ] ‖ ( v h , w h ) ‖ s ≥ ‖ ( u h , z h ) ‖ s , (12) sup v h , w h ∈ V h p × W h p   A h u h , z h , v h , w h v h , w h s ≥ u h , z h s , {:(12)s u p_((v_(h),w_(h))inV_(h)^(p)xxW_(h)^(p))(A_(h)[(u_(h),z_(h)),(v_(h),w_(h))])/(||(v_(h),w_(h))||_(s)) >= ||(u_(h),z_(h))||_(s)",":}\begin{equation*} \sup _{\left(v_{h}, w_{h}\right) \in V_{h}^{p} \times W_{h}^{p}} \frac{A_{h}\left[\left(u_{h}, z_{h}\right),\left(v_{h}, w_{h}\right)\right]}{\left\|\left(v_{h}, w_{h}\right)\right\|_{s}} \geq\left\|\left(u_{h}, z_{h}\right)\right\|_{s}, \tag{12} \end{equation*}(12)sup(vh,wh)∈Vhp×WhpAh[(uh,zh),(vh,wh)]‖(vh,wh)‖s≥‖(uh,zh)‖s,
which is equivalent to the well-posedness of the discrete problem (10).
We now present some stabilization operators that satisfy the criteria discussed above. First, we define the gradient jump penalty, also known as continuous interior penalty (CIP),
J h ( u h , v h ) := ∑ F ∈ F i ∫ F h [ [ ∇ u h ⋅ n ] ] F [ [ ∇ v h ⋅ n ] ] F d s , J h u h , v h := ∑ F ∈ F i   ∫ F   h [ [ ∇ u h ⋅ n ] ] F [ [ ∇ v h ⋅ n ] ] F d s , J_(h)(u_(h),v_(h)):=sum_(F inF_(i))int_(F)h[[gradu_(h)*n]]_(F)[[gradv_(h)*n]]_(F)ds,\mathcal{J}_{h}\left(u_{h}, v_{h}\right):=\sum_{F \in \mathcal{F}_{i}} \int_{F} h \llbracket \nabla u_{h} \cdot n \rrbracket_{F} \llbracket \nabla v_{h} \cdot n \rrbracket_{F} \mathrm{~d} s,Jh(uh,vh):=∑F∈Fi∫Fh[[∇uh⋅n]]F[[∇vh⋅n]]F ds,
to which we add the Galerkin least squares and define the bilinear form s : V h p × V h p s : V h p × V h p s:V_(h)^(p)xxV_(h)^(p)s: V_{h}^{p} \times V_{h}^{p}s:Vhp×Vhp by
(13) s ( u h , v h ) := γ J h ( u h , v h ) + γ h 2 ( L u h , L v h ) Ω + h 2 p ( ∇ u h , ∇ v h ) Ω . (13) s u h , v h := γ J h u h , v h + γ h 2 L u h , L v h Ω + h 2 p ∇ u h , ∇ v h Ω . {:(13)s(u_(h),v_(h)):=gammaJ_(h)(u_(h),v_(h))+gammah^(2)(Lu_(h),Lv_(h))_(Omega)+h^(2p)(gradu_(h),gradv_(h))_(Omega).:}\begin{equation*} s\left(u_{h}, v_{h}\right):=\gamma \mathcal{J}_{h}\left(u_{h}, v_{h}\right)+\gamma h^{2}\left(\mathcal{L} u_{h}, \mathcal{L} v_{h}\right)_{\Omega}+h^{2 p}\left(\nabla u_{h}, \nabla v_{h}\right)_{\Omega} . \tag{13} \end{equation*}(13)s(uh,vh):=γJh(uh,vh)+γh2(Luh,Lvh)Ω+h2p(∇uh,∇vh)Ω.
This represents the primal stabilization that will act on the approximation u h ∈ V h p u h ∈ V h p u_(h)inV_(h)^(p)u_{h} \in V_{h}^{p}uh∈Vhp. The constant parameter γ > 0 γ > 0 gamma > 0\gamma>0γ>0 is meant to tune the stabilization empirically and will be omitted from the analysis. The presence of the gradient term in s s sss is motivated by the a priori H 1 H 1 H^(1)H^{1}H1 bound that the discrete solution must satisfy (uniformly in h h hhh ) in order to be able to use the conditional stability from the continuous level, see also [Bur16]. This gradient term is not needed for piecewise affine approximations. The right-hand side term in (10) for the true solution u u uuu is given by
s ( u , v h ) = h 2 ( f , L v h ) . s u , v h = h 2 f , L v h . s(u,v_(h))=h^(2)(f,Lv_(h)).s\left(u, v_{h}\right)=h^{2}\left(f, \mathcal{L} v_{h}\right) .s(u,vh)=h2(f,Lvh).
For the dual stabilizer s ∗ : W h p × W h p s ∗ : W h p × W h p s^(**):W_(h)^(p)xxW_(h)^(p)s^{*}: W_{h}^{p} \times W_{h}^{p}s∗:Whp×Whp acting on z h ∈ W h p z h ∈ W h p z_(h)inW_(h)^(p)z_{h} \in W_{h}^{p}zh∈Whp, we take
(14) s ∗ ( z h , w h ) = ( ∇ z h , ∇ w h ) Ω . (14) s ∗ z h , w h = ∇ z h , ∇ w h Ω . {:(14)s^(**)(z_(h),w_(h))=(gradz_(h),gradw_(h))_(Omega).:}\begin{equation*} s^{*}\left(z_{h}, w_{h}\right)=\left(\nabla z_{h}, \nabla w_{h}\right)_{\Omega} . \tag{14} \end{equation*}(14)s∗(zh,wh)=(∇zh,∇wh)Ω.
With this choice, since the finite element space contains piecewise polynomials we have that s s sss induces a norm on V h p V h p V_(h)^(p)V_{h}^{p}Vhp, and also that s ∗ s ∗ s^(**)s^{*}s∗ induces a norm on W h p W h p W_(h)^(p)W_{h}^{p}Whp. We will denote these norms by
(15) ‖ v h ‖ V h p := s ( v h , v h ) 1 2 , ‖ w h ‖ W h p := s ∗ ( w h , w h ) 1 2 . (15) v h V h p := s v h , v h 1 2 , w h W h p := s ∗ w h , w h 1 2 . {:(15)||v_(h)||_(V_(h)^(p)):=s(v_(h),v_(h))^((1)/(2))","quad||w_(h)||_(W_(h)^(p)):=s^(**)(w_(h),w_(h))^((1)/(2)).:}\begin{equation*} \left\|v_{h}\right\|_{V_{h}^{p}}:=s\left(v_{h}, v_{h}\right)^{\frac{1}{2}}, \quad\left\|w_{h}\right\|_{W_{h}^{p}}:=s^{*}\left(w_{h}, w_{h}\right)^{\frac{1}{2}} . \tag{15} \end{equation*}(15)‖vh‖Vhp:=s(vh,vh)12,‖wh‖Whp:=s∗(wh,wh)12.
Note that (11) is indeed a norm on V h p × W h p V h p × W h p V_(h)^(p)xxW_(h)^(p)V_{h}^{p} \times W_{h}^{p}Vhp×Whp, thus guaranteeing the inf-sup condition (12) is satisfied. This means that the discrete problem is well-posed.
From the standard approximation (7) and the trace inequality (6), we have the approximation bound in the stabilization norm
(16) ‖ u − π h u ‖ V h p ≤ C h p ‖ u ‖ H p + 1 ( Ω ) (16) u − π h u V h p ≤ C h p ‖ u ‖ H p + 1 ( Ω ) {:(16)||u-pi_(h)u||_(V_(h)^(p)) <= Ch^(p)||u||_(H^(p+1)(Omega)):}\begin{equation*} \left\|u-\pi_{h} u\right\|_{V_{h}^{p}} \leq C h^{p}\|u\|_{H^{p+1}(\Omega)} \tag{16} \end{equation*}(16)‖u−πhu‖Vhp≤Chp‖u‖Hp+1(Ω)
Let us now motivate the choice for the gradient jump penalty J h J h J_(h)\mathcal{J}_{h}Jh and the Galerkin least squares stabilization in the primal stabilizer (13) by the following continuity result.
Lemma 1. Let u ∈ H 1 ( Ω ) , u h ∈ V h p , w ∈ H 0 1 ( Ω ) u ∈ H 1 ( Ω ) , u h ∈ V h p , w ∈ H 0 1 ( Ω ) u inH^(1)(Omega),u_(h)inV_(h)^(p),w inH_(0)^(1)(Omega)u \in H^{1}(\Omega), u_{h} \in V_{h}^{p}, w \in H_{0}^{1}(\Omega)u∈H1(Ω),uh∈Vhp,w∈H01(Ω) and the norm ‖ w ‖ # := h − 1 ‖ w ‖ Ω + ‖ ∇ w ‖ Ω ‖ w ‖ # := h − 1 ‖ w ‖ Ω + ‖ ∇ w ‖ Ω ||w||_(#):=h^(-1)||w||_(Omega)+||grad w||_(Omega)\|w\|_{\#}:=h^{-1}\|w\|_{\Omega}+\|\nabla w\|_{\Omega}‖w‖#:=h−1‖w‖Ω+‖∇w‖Ω. Then
a ( u − u h , w ) ≤ C ‖ ( u − u h , 0 ) ‖ s ‖ w ‖ # ⋅ a u − u h , w ≤ C u − u h , 0 s ‖ w ‖ # ⋅ a(u-u_(h),w) <= C||(u-u_(h),0)||_(s)||w||_(#)*a\left(u-u_{h}, w\right) \leq C\left\|\left(u-u_{h}, 0\right)\right\|_{s}\|w\|_{\#} \cdota(u−uh,w)≤C‖(u−uh,0)‖s‖w‖#⋅
Proof. Integrating by parts and using the continuous trace inequality (6)
a ( u − u h , w ) = ( ∇ ( u − u h ) , ∇ w ) Ω − k 2 ( u − u h , w ) Ω = ( L ( u − u h ) , w ) Ω + ∑ F ∈ F i ∫ F [ [ ∇ ( u − u h ) ⋅ n ] ] F w d s = ( h L ( u − u h ) , h − 1 w ) Ω + J h ( u − u h , u − u h ) 1 2 ( h − 1 ‖ w ‖ Ω + ‖ ∇ w ‖ Ω ) ≤ C ‖ u − u h ‖ V h p ( h − 1 ‖ w ‖ Ω + ‖ ∇ w ‖ Ω ) a u − u h , w = ∇ u − u h , ∇ w Ω − k 2 u − u h , w Ω = L u − u h , w Ω + ∑ F ∈ F i   ∫ F   [ [ ∇ u − u h ⋅ n ] ] F w d s = h L u − u h , h − 1 w Ω + J h u − u h , u − u h 1 2 h − 1 ‖ w ‖ Ω + ‖ ∇ w ‖ Ω ≤ C u − u h V h p h − 1 ‖ w ‖ Ω + ‖ ∇ w ‖ Ω {:[a(u-u_(h),w)=(grad(u-u_(h)),grad w)_(Omega)-k^(2)(u-u_(h),w)_(Omega)],[=(L(u-u_(h)),w)_(Omega)+sum_(F inF_(i))int_(F)[[grad(u-u_(h))*n]]_(F)wds],[=(hL(u-u_(h)),h^(-1)w)_(Omega)+J_(h)(u-u_(h),u-u_(h))^((1)/(2))(h^(-1)||w||_(Omega)+||grad w||_(Omega))],[ <= C||u-u_(h)||_(V_(h)^(p))(h^(-1)||w||_(Omega)+||grad w||_(Omega))]:}\begin{aligned} a\left(u-u_{h}, w\right) & =\left(\nabla\left(u-u_{h}\right), \nabla w\right)_{\Omega}-k^{2}\left(u-u_{h}, w\right)_{\Omega} \\ & =\left(\mathcal{L}\left(u-u_{h}\right), w\right)_{\Omega}+\sum_{F \in \mathcal{F}_{i}} \int_{F} \llbracket \nabla\left(u-u_{h}\right) \cdot n \rrbracket_{F} w \mathrm{~d} s \\ & =\left(h \mathcal{L}\left(u-u_{h}\right), h^{-1} w\right)_{\Omega}+\mathcal{J}_{h}\left(u-u_{h}, u-u_{h}\right)^{\frac{1}{2}}\left(h^{-1}\|w\|_{\Omega}+\|\nabla w\|_{\Omega}\right) \\ & \leq C\left\|u-u_{h}\right\|_{V_{h}^{p}}\left(h^{-1}\|w\|_{\Omega}+\|\nabla w\|_{\Omega}\right) \end{aligned}a(u−uh,w)=(∇(u−uh),∇w)Ω−k2(u−uh,w)Ω=(L(u−uh),w)Ω+∑F∈Fi∫F[[∇(u−uh)⋅n]]Fw ds=(hL(u−uh),h−1w)Ω+Jh(u−uh,u−uh)12(h−1‖w‖Ω+‖∇w‖Ω)≤C‖u−uh‖Vhp(h−1‖w‖Ω+‖∇w‖Ω)

3.4 Error analysis

We start by proving that the stabilization/regularization converges optimally for smooth solutions in the norm ‖ ( u h − π h u , z h ) ‖ s u h − π h u , z h s ||(u_(h)-pi_(h)u,z_(h))||_(s)\left\|\left(u_{h}-\pi_{h} u, z_{h}\right)\right\|_{s}‖(uh−πhu,zh)‖s. Note that this is independent of the stability of the continuous problem.
Lemma 2. Let u ∈ H p + 1 ( Ω ) u ∈ H p + 1 ( Ω ) u inH^(p+1)(Omega)u \in H^{p+1}(\Omega)u∈Hp+1(Ω) be the solution to (1) and let ( u h , z h ) ∈ V h p × W h p u h , z h ∈ V h p × W h p (u_(h),z_(h))inV_(h)^(p)xxW_(h)^(p)\left(u_{h}, z_{h}\right) \in V_{h}^{p} \times W_{h}^{p}(uh,zh)∈Vhp×Whp be the finite element solution to (10). Then
‖ ( u − u h , z h ) ‖ s ≤ C ( h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω ) u − u h , z h s ≤ C h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω ||(u-u_(h),z_(h))||_(s) <= C(h^(p)||u||_(H^(p+1)(Omega))+||delta g||_(omega))\left\|\left(u-u_{h}, z_{h}\right)\right\|_{s} \leq C\left(h^{p}\|u\|_{H^{p+1}(\Omega)}+\|\delta g\|_{\omega}\right)‖(u−uh,zh)‖s≤C(hp‖u‖Hp+1(Ω)+‖δg‖ω)
Proof. From the triangle inequality we have that
‖ ( u − u h , z h ) ‖ s ≤ ‖ ( u − π h u , 0 ) ‖ s + ‖ ( u h − π h u , z h ) ‖ s u − u h , z h s ≤ u − π h u , 0 s + u h − π h u , z h s ||(u-u_(h),z_(h))||_(s) <= ||(u-pi_(h)u,0)||_(s)+||(u_(h)-pi_(h)u,z_(h))||_(s)\left\|\left(u-u_{h}, z_{h}\right)\right\|_{s} \leq\left\|\left(u-\pi_{h} u, 0\right)\right\|_{s}+\left\|\left(u_{h}-\pi_{h} u, z_{h}\right)\right\|_{s}‖(u−uh,zh)‖s≤‖(u−πhu,0)‖s+‖(uh−πhu,zh)‖s
The first term is bounded from the approximation estimate (16) and we only need to focus on the second term. Let us start by expanding
(17) ‖ ( u h − π h u , z h ) ‖ s 2 = ( u h , u h − π h u ) ω + s ( u h , u h − π h u ) + s ∗ ( z h , z h ) − ( π h u , u h − π h u ) ω − s ( π h u , u h − π h u ) (17) u h − π h u , z h s 2 = u h , u h − π h u ω + s u h , u h − π h u + s ∗ z h , z h − π h u , u h − π h u ω − s π h u , u h − π h u {:[(17)||(u_(h)-pi_(h)u,z_(h))||_(s)^(2)=(u_(h),u_(h)-pi_(h)u)_(omega)+s(u_(h),u_(h)-pi_(h)u)+s^(**)(z_(h),z_(h))],[-(pi_(h)u,u_(h)-pi_(h)u)_(omega)-s(pi_(h)u,u_(h)-pi_(h)u)]:}\begin{align*} \left\|\left(u_{h}-\pi_{h} u, z_{h}\right)\right\|_{s}^{2}= & \left(u_{h}, u_{h}-\pi_{h} u\right)_{\omega}+s\left(u_{h}, u_{h}-\pi_{h} u\right)+s^{*}\left(z_{h}, z_{h}\right) \tag{17}\\ & -\left(\pi_{h} u, u_{h}-\pi_{h} u\right)_{\omega}-s\left(\pi_{h} u, u_{h}-\pi_{h} u\right) \end{align*}(17)‖(uh−πhu,zh)‖s2=(uh,uh−πhu)ω+s(uh,uh−πhu)+s∗(zh,zh)−(πhu,uh−πhu)ω−s(πhu,uh−πhu)
For the first row of this expansion we use ( 10 ) 1 ( 10 ) 1 (10)_(1)(10)_{1}(10)1 and we have that
( u h , u h − π h u ) ω + s ( u h , u h − π h u ) + a ( u h − π h u , z h ) = ( g ~ , u h − π h u ) ω + s ( u , u h − π h u ) . u h , u h − π h u ω + s u h , u h − π h u + a u h − π h u , z h = g ~ , u h − π h u ω + s u , u h − π h u . (u_(h),u_(h)-pi_(h)u)_(omega)+s(u_(h),u_(h)-pi_(h)u)+a(u_(h)-pi_(h)u,z_(h))=(( tilde(g)),u_(h)-pi_(h)u)_(omega)+s(u,u_(h)-pi_(h)u).\left(u_{h}, u_{h}-\pi_{h} u\right)_{\omega}+s\left(u_{h}, u_{h}-\pi_{h} u\right)+a\left(u_{h}-\pi_{h} u, z_{h}\right)=\left(\tilde{g}, u_{h}-\pi_{h} u\right)_{\omega}+s\left(u, u_{h}-\pi_{h} u\right) .(uh,uh−πhu)ω+s(uh,uh−πhu)+a(uh−πhu,zh)=(g~,uh−πhu)ω+s(u,uh−πhu).
We can use ( 10 ) 2 ( 10 ) 2 (10)_(2)(10)_{2}(10)2 to write
a ( u h − π h u , z h ) = s ∗ ( z h , z h ) + ( f , z h ) Ω − a ( π h u , z h ) = s ∗ ( z h , z h ) + a ( u − π h u , z h ) a u h − π h u , z h = s ∗ z h , z h + f , z h Ω − a π h u , z h = s ∗ z h , z h + a u − π h u , z h {:[a(u_(h)-pi_(h)u,z_(h))=s^(**)(z_(h),z_(h))+(f,z_(h))_(Omega)-a(pi_(h)u,z_(h))],[=s^(**)(z_(h),z_(h))+a(u-pi_(h)u,z_(h))]:}\begin{aligned} a\left(u_{h}-\pi_{h} u, z_{h}\right) & =s^{*}\left(z_{h}, z_{h}\right)+\left(f, z_{h}\right)_{\Omega}-a\left(\pi_{h} u, z_{h}\right) \\ & =s^{*}\left(z_{h}, z_{h}\right)+a\left(u-\pi_{h} u, z_{h}\right) \end{aligned}a(uh−πhu,zh)=s∗(zh,zh)+(f,zh)Ω−a(πhu,zh)=s∗(zh,zh)+a(u−πhu,zh)
thus giving that
( u h , u h − π h u ) ω + s ( u h , u h − π h u ) + s ∗ ( z h , z h ) = ( g ~ , u h − π h u ) ω + s ( u , u h − π h u ) − a ( u − π h u , z h ) u h , u h − π h u ω + s u h , u h − π h u + s ∗ z h , z h = g ~ , u h − π h u ω + s u , u h − π h u − a u − π h u , z h {:[(u_(h),u_(h)-pi_(h)u)_(omega)+s(u_(h),u_(h)-pi_(h)u)+s^(**)(z_(h),z_(h))=(( tilde(g)),u_(h)-pi_(h)u)_(omega)+s(u,u_(h)-pi_(h)u)],[-a(u-pi_(h)u,z_(h))]:}\begin{aligned} \left(u_{h}, u_{h}-\pi_{h} u\right)_{\omega}+s\left(u_{h}, u_{h}-\pi_{h} u\right)+s^{*}\left(z_{h}, z_{h}\right)= & \left(\tilde{g}, u_{h}-\pi_{h} u\right)_{\omega}+s\left(u, u_{h}-\pi_{h} u\right) \\ & -a\left(u-\pi_{h} u, z_{h}\right) \end{aligned}(uh,uh−πhu)ω+s(uh,uh−πhu)+s∗(zh,zh)=(g~,uh−πhu)ω+s(u,uh−πhu)−a(u−πhu,zh)
Returning to the stabilization norm (17), we obtain that
(18) ‖ ( u h − π h u , z h ) ‖ s 2 = ( g ~ − π h u , u h − π h u ) ω + a ( π h u − u , z h ) + s ( u − π h u , u h − π h u ) (18) u h − π h u , z h s 2 = g ~ − π h u , u h − π h u ω + a π h u − u , z h + s u − π h u , u h − π h u {:(18)||(u_(h)-pi_(h)u,z_(h))||_(s)^(2)=(( tilde(g))-pi_(h)u,u_(h)-pi_(h)u)_(omega)+a(pi_(h)u-u,z_(h))+s(u-pi_(h)u,u_(h)-pi_(h)u):}\begin{equation*} \left\|\left(u_{h}-\pi_{h} u, z_{h}\right)\right\|_{s}^{2}=\left(\tilde{g}-\pi_{h} u, u_{h}-\pi_{h} u\right)_{\omega}+a\left(\pi_{h} u-u, z_{h}\right)+s\left(u-\pi_{h} u, u_{h}-\pi_{h} u\right) \tag{18} \end{equation*}(18)‖(uh−πhu,zh)‖s2=(g~−πhu,uh−πhu)ω+a(πhu−u,zh)+s(u−πhu,uh−πhu)
We now discuss the convergence of these three terms. First, by Cauchy-Schwarz and standard approximation (7) we bound
( g ~ − π h u , u h − π h u ) ω ≤ ‖ g ~ − π h u ‖ ω ‖ u h − π h u ‖ ω ≤ C ( h p + 1 ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω ) ‖ u h − π h u ‖ ω . g ~ − π h u , u h − π h u ω ≤ g ~ − π h u ω u h − π h u ω ≤ C h p + 1 ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω u h − π h u ω . {:[(( tilde(g))-pi_(h)u,u_(h)-pi_(h)u)_(omega) <= ||( tilde(g))-pi_(h)u||_(omega)||u_(h)-pi_(h)u||_(omega)],[ <= C(h^(p+1)||u||_(H^(p+1)(Omega))+||delta g||_(omega))||u_(h)-pi_(h)u||_(omega).]:}\begin{aligned} \left(\tilde{g}-\pi_{h} u, u_{h}-\pi_{h} u\right)_{\omega} & \leq\left\|\tilde{g}-\pi_{h} u\right\|_{\omega}\left\|u_{h}-\pi_{h} u\right\|_{\omega} \\ & \leq C\left(h^{p+1}\|u\|_{H^{p+1}(\Omega)}+\|\delta g\|_{\omega}\right)\left\|u_{h}-\pi_{h} u\right\|_{\omega} . \end{aligned}(g~−πhu,uh−πhu)ω≤‖g~−πhu‖ω‖uh−πhu‖ω≤C(hp+1‖u‖Hp+1(Ω)+‖δg‖ω)‖uh−πhu‖ω.
Second, since π h π h pi_(h)\pi_{h}πh is the L 2 L 2 L^(2)L^{2}L2-projection we use orthogonality and again the approximation (7) to write
a ( π h u − u , z h ) = ( ∇ ( π h u − u ) , ∇ z h ) Ω ≤ C h p ‖ u ‖ H p + 1 ( Ω ) ‖ z h ‖ W h p a π h u − u , z h = ∇ π h u − u , ∇ z h Ω ≤ C h p ‖ u ‖ H p + 1 ( Ω ) z h W h p a(pi_(h)u-u,z_(h))=(grad(pi_(h)u-u),gradz_(h))_(Omega) <= Ch^(p)||u||_(H^(p+1)(Omega))||z_(h)||_(W_(h)^(p))a\left(\pi_{h} u-u, z_{h}\right)=\left(\nabla\left(\pi_{h} u-u\right), \nabla z_{h}\right)_{\Omega} \leq C h^{p}\|u\|_{H^{p+1}(\Omega)}\left\|z_{h}\right\|_{W_{h}^{p}}a(πhu−u,zh)=(∇(πhu−u),∇zh)Ω≤Chp‖u‖Hp+1(Ω)‖zh‖Whp
Third,
s ( u − π h u , u h − π h u ) ≤ ‖ u − π h u ‖ V h p ‖ u h − π h u ‖ V h p ≤ C h p ‖ u ‖ H p + 1 ( Ω ) ‖ u h − π h u ‖ V h p , s u − π h u , u h − π h u ≤ u − π h u V h p u h − π h u V h p ≤ C h p ‖ u ‖ H p + 1 ( Ω ) u h − π h u V h p , s(u-pi_(h)u,u_(h)-pi_(h)u) <= ||u-pi_(h)u||_(V_(h)^(p))||u_(h)-pi_(h)u||_(V_(h)^(p)) <= Ch^(p)||u||_(H^(p+1)(Omega))||u_(h)-pi_(h)u||_(V_(h)^(p)),s\left(u-\pi_{h} u, u_{h}-\pi_{h} u\right) \leq\left\|u-\pi_{h} u\right\|_{V_{h}^{p}}\left\|u_{h}-\pi_{h} u\right\|_{V_{h}^{p}} \leq C h^{p}\|u\|_{H^{p+1}(\Omega)}\left\|u_{h}-\pi_{h} u\right\|_{V_{h}^{p}},s(u−πhu,uh−πhu)≤‖u−πhu‖Vhp‖uh−πhu‖Vhp≤Chp‖u‖Hp+1(Ω)‖uh−πhu‖Vhp,
where we used the approximation bound (16). Collecting the bounds for the terms in (18), we conclude that
‖ ( u h − π h u , z h ) ‖ s 2 ≤ C ( h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω ) ‖ ( u h − π h u , z h ) ‖ s u h − π h u , z h s 2 ≤ C h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω u h − π h u , z h s ||(u_(h)-pi_(h)u,z_(h))||_(s)^(2) <= C(h^(p)||u||_(H^(p+1)(Omega))+||delta g||_(omega))||(u_(h)-pi_(h)u,z_(h))||_(s)\left\|\left(u_{h}-\pi_{h} u, z_{h}\right)\right\|_{s}^{2} \leq C\left(h^{p}\|u\|_{H^{p+1}(\Omega)}+\|\delta g\|_{\omega}\right)\left\|\left(u_{h}-\pi_{h} u, z_{h}\right)\right\|_{s}‖(uh−πhu,zh)‖s2≤C(hp‖u‖Hp+1(Ω)+‖δg‖ω)‖(uh−πhu,zh)‖s
and hence
‖ ( u h − π h u , z h ) ‖ s ≤ C ( h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω ) . u h − π h u , z h s ≤ C h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω . ||(u_(h)-pi_(h)u,z_(h))||_(s) <= C(h^(p)||u||_(H^(p+1)(Omega))+||delta g||_(omega)).\left\|\left(u_{h}-\pi_{h} u, z_{h}\right)\right\|_{s} \leq C\left(h^{p}\|u\|_{H^{p+1}(\Omega)}+\|\delta g\|_{\omega}\right) .‖(uh−πhu,zh)‖s≤C(hp‖u‖Hp+1(Ω)+‖δg‖ω). ◻ ◻ ◻\square◻
We continue by proving that the dual norm of the residual can be controlled by the stabilization.
Lemma 3. Let u ∈ H p + 1 ( Ω ) u ∈ H p + 1 ( Ω ) u inH^(p+1)(Omega)u \in H^{p+1}(\Omega)u∈Hp+1(Ω) be the solution to (1) and let ( u h , z h ) ∈ V h p × W h p u h , z h ∈ V h p × W h p (u_(h),z_(h))inV_(h)^(p)xxW_(h)^(p)\left(u_{h}, z_{h}\right) \in V_{h}^{p} \times W_{h}^{p}(uh,zh)∈Vhp×Whp be the finite element solution to (10). Let r ∈ H − 1 ( Ω ) r ∈ H − 1 ( Ω ) r inH^(-1)(Omega)r \in H^{-1}(\Omega)r∈H−1(Ω) be the residual ⟨ r , w ⟩ = a ( u − u h , w ) , w ∈ H 0 1 ( Ω ) ⟨ r , w ⟩ = a u − u h , w , w ∈ H 0 1 ( Ω ) (:r,w:)=a(u-u_(h),w),w inH_(0)^(1)(Omega)\langle r, w\rangle=a\left(u-u_{h}, w\right), w \in H_{0}^{1}(\Omega)⟨r,w⟩=a(u−uh,w),w∈H01(Ω). Then
‖ r ‖ H − 1 ( Ω ) ≤ ‖ ( u − u h , z h ) ‖ s ≤ C ( h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω ) . ‖ r ‖ H − 1 ( Ω ) ≤ u − u h , z h s ≤ C h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω . ||r||_(H^(-1)(Omega)) <= ||(u-u_(h),z_(h))||_(s) <= C(h^(p)||u||_(H^(p+1)(Omega))+||delta g||_(omega)).\|r\|_{H^{-1}(\Omega)} \leq\left\|\left(u-u_{h}, z_{h}\right)\right\|_{s} \leq C\left(h^{p}\|u\|_{H^{p+1}(\Omega)}+\|\delta g\|_{\omega}\right) .‖r‖H−1(Ω)≤‖(u−uh,zh)‖s≤C(hp‖u‖Hp+1(Ω)+‖δg‖ω).
Proof. Let us first recall from (10) that
a ( u h , w h ) − s ∗ ( z h , w h ) = ( f , w h ) Ω = a ( u , w h ) , ∀ w h ∈ W h p , a u h , w h − s ∗ z h , w h = f , w h Ω = a u , w h , ∀ w h ∈ W h p , a(u_(h),w_(h))-s^(**)(z_(h),w_(h))=(f,w_(h))_(Omega)=a(u,w_(h)),quad AAw_(h)inW_(h)^(p),a\left(u_{h}, w_{h}\right)-s^{*}\left(z_{h}, w_{h}\right)=\left(f, w_{h}\right)_{\Omega}=a\left(u, w_{h}\right), \quad \forall w_{h} \in W_{h}^{p},a(uh,wh)−s∗(zh,wh)=(f,wh)Ω=a(u,wh),∀wh∈Whp,
which gives the Galerkin orthogonality
a ( u − u h , w h ) = − s ∗ ( z h , w h ) , ∀ w h ∈ W h p . a u − u h , w h = − s ∗ z h , w h , ∀ w h ∈ W h p . a(u-u_(h),w_(h))=-s^(**)(z_(h),w_(h)),quad AAw_(h)inW_(h)^(p).a\left(u-u_{h}, w_{h}\right)=-s^{*}\left(z_{h}, w_{h}\right), \quad \forall w_{h} \in W_{h}^{p} .a(u−uh,wh)=−s∗(zh,wh),∀wh∈Whp.
Using an interpolant π s z π s z pi_(sz)\pi_{s z}πsz of Scott-Zhang type, we write the residual
⟨ r , w ⟩ = a ( u − u h , w ) = a ( u − u h , w − π s z w ) + a ( u − u h , π s z w ) = a ( u − u h , w − π s z w ) − s ∗ ( z h , π s z w ) ⟨ r , w ⟩ = a u − u h , w = a u − u h , w − π s z w + a u − u h , π s z w = a u − u h , w − π s z w − s ∗ z h , π s z w {:[(:r","w:)=a(u-u_(h),w)=a(u-u_(h),w-pi_(sz)w)+a(u-u_(h),pi_(sz)w)],[=a(u-u_(h),w-pi_(sz)w)-s^(**)(z_(h),pi_(sz)w)]:}\begin{aligned} \langle r, w\rangle=a\left(u-u_{h}, w\right) & =a\left(u-u_{h}, w-\pi_{s z} w\right)+a\left(u-u_{h}, \pi_{s z} w\right) \\ & =a\left(u-u_{h}, w-\pi_{s z} w\right)-s^{*}\left(z_{h}, \pi_{s z} w\right) \end{aligned}⟨r,w⟩=a(u−uh,w)=a(u−uh,w−πszw)+a(u−uh,πszw)=a(u−uh,w−πszw)−s∗(zh,πszw)
From the continuity Lemma 1 and interpolation inequality (7) we bound the first term
a ( u − u h , w − π s z ) ≤ ‖ ( u − u h , 0 ) ‖ s ‖ w − π s z w ‖ # ≤ ‖ ( u − u h , 0 ) ‖ s ‖ w ‖ H 1 ( Ω ) a u − u h , w − π s z ≤ u − u h , 0 s w − π s z w # ≤ u − u h , 0 s ‖ w ‖ H 1 ( Ω ) a(u-u_(h),w-pi_(sz)) <= ||(u-u_(h),0)||_(s)||w-pi_(sz)w||_(#) <= ||(u-u_(h),0)||_(s)||w||_(H^(1)(Omega))a\left(u-u_{h}, w-\pi_{s z}\right) \leq\left\|\left(u-u_{h}, 0\right)\right\|_{s}\left\|w-\pi_{s z} w\right\|_{\#} \leq\left\|\left(u-u_{h}, 0\right)\right\|_{s}\|w\|_{H^{1}(\Omega)}a(u−uh,w−πsz)≤‖(u−uh,0)‖s‖w−πszw‖#≤‖(u−uh,0)‖s‖w‖H1(Ω)
The second term is bounded by Cauchy-Schwarz and the stability of the interpolant
− s ∗ ( z h , π s z ) ≤ ‖ ∇ z h ‖ Ω ‖ π s z w ‖ Ω ≤ ‖ ( 0 , z h ) ‖ s ‖ w ‖ H 1 ( Ω ) , − s ∗ z h , π s z ≤ ∇ z h Ω π s z w Ω ≤ 0 , z h s ‖ w ‖ H 1 ( Ω ) , -s^(**)(z_(h),pi_(sz)) <= ||gradz_(h)||_(Omega)||pi_(sz)w||_(Omega) <= ||(0,z_(h))||_(s)||w||_(H^(1)(Omega)),-s^{*}\left(z_{h}, \pi_{s z}\right) \leq\left\|\nabla z_{h}\right\|_{\Omega}\left\|\pi_{s z} w\right\|_{\Omega} \leq\left\|\left(0, z_{h}\right)\right\|_{s}\|w\|_{H^{1}(\Omega)},−s∗(zh,πsz)≤‖∇zh‖Ω‖πszw‖Ω≤‖(0,zh)‖s‖w‖H1(Ω),
and we obtain that
⟨ r , w ⟩ ≤ ‖ ( u − u h , z h ) ‖ s ‖ w ‖ H 1 ( Ω ) , ⟨ r , w ⟩ ≤ u − u h , z h s ‖ w ‖ H 1 ( Ω ) , (:r,w:) <= ||(u-u_(h),z_(h))||_(s)||w||_(H^(1)(Omega)),\langle r, w\rangle \leq\left\|\left(u-u_{h}, z_{h}\right)\right\|_{s}\|w\|_{H^{1}(\Omega)},⟨r,w⟩≤‖(u−uh,zh)‖s‖w‖H1(Ω),
with the conclusion following from
‖ r ‖ H − 1 ( Ω ) ≤ ‖ ( u − u h , z h ) ‖ s ‖ r ‖ H − 1 ( Ω ) ≤ u − u h , z h s ||r||_(H^(-1)(Omega)) <= ||(u-u_(h),z_(h))||_(s)\|r\|_{H^{-1}(\Omega)} \leq\left\|\left(u-u_{h}, z_{h}\right)\right\|_{s}‖r‖H−1(Ω)≤‖(u−uh,zh)‖s
and Lemma 2. ◻ ◻ ◻\square◻
We will now prove an error estimate on the solution in a target domain B B BBB, using the conditional stability from Section 2. We will write a general version, with implicit dependence on the frequency-this can be made explicit for particular geometric settings satisfying a convexity condition.
Theorem 1. Let u ∈ H p + 1 ( Ω ) u ∈ H p + 1 ( Ω ) u inH^(p+1)(Omega)u \in H^{p+1}(\Omega)u∈Hp+1(Ω) be the solution to (1) and let ( u h , z h ) ∈ V h p × W h p u h , z h ∈ V h p × W h p (u_(h),z_(h))inV_(h)^(p)xxW_(h)^(p)\left(u_{h}, z_{h}\right) \in V_{h}^{p} \times W_{h}^{p}(uh,zh)∈Vhp×Whp be the finite element solution to (10). Then
‖ u − u h ‖ H 1 ( B ) ≤ C s t ( k ) h α p ( ‖ u ‖ H p + 1 ( Ω ) + h − p ‖ δ g ‖ ω ) . u − u h H 1 ( B ) ≤ C s t ( k ) h α p ‖ u ‖ H p + 1 ( Ω ) + h − p ‖ δ g ‖ ω . ||u-u_(h)||_(H^(1)(B)) <= C_(st)(k)h^(alpha p)(||u||_(H^(p+1)(Omega))+h^(-p)||delta g||_(omega)).\left\|u-u_{h}\right\|_{H^{1}(B)} \leq C_{s t}(k) h^{\alpha p}\left(\|u\|_{H^{p+1}(\Omega)}+h^{-p}\|\delta g\|_{\omega}\right) .‖u−uh‖H1(B)≤Cst(k)hαp(‖u‖Hp+1(Ω)+h−p‖δg‖ω).
Proof. We consider the error u − u h u − u h u-u_(h)u-u_{h}u−uh together with the residual ⟨ r , w ⟩ = a ( u − u h , w ) ⟨ r , w ⟩ = a u − u h , w (:r,w:)=a(u-u_(h),w)\langle r, w\rangle=a\left(u-u_{h}, w\right)⟨r,w⟩=a(u−uh,w), for w ∈ H 0 1 ( Ω ) w ∈ H 0 1 ( Ω ) w inH_(0)^(1)(Omega)w \in H_{0}^{1}(\Omega)w∈H01(Ω), and apply the conditional stability estimate (4)
‖ u − u h ‖ H 1 ( B ) ≤ C s t ( k ) ( ‖ u − u h ‖ L 2 ( ω ) + ‖ r ‖ H − 1 ( Ω ) ) α ( ‖ u − u h ‖ L 2 ( Ω ) + ‖ r ‖ H − 1 ( Ω ) ) 1 − α . u − u h H 1 ( B ) ≤ C s t ( k ) u − u h L 2 ( ω ) + ‖ r ‖ H − 1 ( Ω ) α u − u h L 2 ( Ω ) + ‖ r ‖ H − 1 ( Ω ) 1 − α . ||u-u_(h)||_(H^(1)(B)) <= C_(st)(k)(||u-u_(h)||_(L^(2)(omega))+||r||_(H^(-1)(Omega)))^(alpha)(||u-u_(h)||_(L^(2)(Omega))+||r||_(H^(-1)(Omega)))^(1-alpha).\left\|u-u_{h}\right\|_{H^{1}(B)} \leq C_{s t}(k)\left(\left\|u-u_{h}\right\|_{L^{2}(\omega)}+\|r\|_{H^{-1}(\Omega)}\right)^{\alpha}\left(\left\|u-u_{h}\right\|_{L^{2}(\Omega)}+\|r\|_{H^{-1}(\Omega)}\right)^{1-\alpha} .‖u−uh‖H1(B)≤Cst(k)(‖u−uh‖L2(ω)+‖r‖H−1(Ω))α(‖u−uh‖L2(Ω)+‖r‖H−1(Ω))1−α.
Note that ‖ u − u h ‖ L 2 ( ω ) u − u h L 2 ( ω ) ||u-u_(h)||_(L^(2)(omega))\left\|u-u_{h}\right\|_{L^{2}(\omega)}‖u−uh‖L2(ω) is bounded from Lemma 2 and the residual ‖ r ‖ H − 1 ( Ω ) ‖ r ‖ H − 1 ( Ω ) ||r||_(H^(-1)(Omega))\|r\|_{H^{-1}(\Omega)}‖r‖H−1(Ω) is bounded from Lemma 3. Bounding
‖ u − u h ‖ L 2 ( Ω ) ≤ C h − p ‖ ( u − u h , 0 ) ‖ s ≤ C ( ‖ u ‖ H p + 1 ( Ω ) + h − p ‖ δ g ‖ ω ) , u − u h L 2 ( Ω ) ≤ C h − p u − u h , 0 s ≤ C ‖ u ‖ H p + 1 ( Ω ) + h − p ‖ δ g ‖ ω , ||u-u_(h)||_(L^(2)(Omega)) <= Ch^(-p)||(u-u_(h),0)||_(s) <= C(||u||_(H^(p+1)(Omega))+h^(-p)||delta g||_(omega)),\left\|u-u_{h}\right\|_{L^{2}(\Omega)} \leq C h^{-p}\left\|\left(u-u_{h}, 0\right)\right\|_{s} \leq C\left(\|u\|_{H^{p+1}(\Omega)}+h^{-p}\|\delta g\|_{\omega}\right),‖u−uh‖L2(Ω)≤Ch−p‖(u−uh,0)‖s≤C(‖u‖Hp+1(Ω)+h−p‖δg‖ω),
we obtain that
‖ u − u h ‖ H 1 ( B ) ≤ C s t ( k ) h ( α − 1 ) p ( h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω ) ≤ C s t ( k ) h α p ( ‖ u ‖ H p + 1 ( Ω ) + h − p ‖ δ g ‖ ω ) u − u h H 1 ( B ) ≤ C s t ( k ) h ( α − 1 ) p h p ‖ u ‖ H p + 1 ( Ω ) + ‖ δ g ‖ ω ≤ C s t ( k ) h α p ‖ u ‖ H p + 1 ( Ω ) + h − p ‖ δ g ‖ ω {:[||u-u_(h)||_(H^(1)(B)) <= C_(st)(k)h^((alpha-1)p)(h^(p)||u||_(H^(p+1)(Omega))+||delta g||_(omega))],[ <= C_(st)(k)h^(alpha p)(||u||_(H^(p+1)(Omega))+h^(-p)||delta g||_(omega))]:}\begin{aligned} \left\|u-u_{h}\right\|_{H^{1}(B)} & \leq C_{s t}(k) h^{(\alpha-1) p}\left(h^{p}\|u\|_{H^{p+1}(\Omega)}+\|\delta g\|_{\omega}\right) \\ & \leq C_{s t}(k) h^{\alpha p}\left(\|u\|_{H^{p+1}(\Omega)}+h^{-p}\|\delta g\|_{\omega}\right) \end{aligned}‖u−uh‖H1(B)≤Cst(k)h(α−1)p(hp‖u‖Hp+1(Ω)+‖δg‖ω)≤Cst(k)hαp(‖u‖Hp+1(Ω)+h−p‖δg‖ω)
Note that in Theorem 1 the convergence rate increases linearly with the polynomial degree p p ppp of the approximation and in the case of well-posed problems (with α = 1 α = 1 alpha=1\alpha=1α=1 ) we recover the optimal rate p p ppp. In the presence of noise, the error bound can degenerate and perturbations in data get amplified as p p ppp increases, meaning that in practical settings the advantages of high-order methods might be reduced. Also, one should consider a certain threshold after which refining the mesh becomes detrimental.

4 Numerical experiments

We present numerical experiments for the Helmholtz unique continuation problem (1) solved with the high-order conforming method (10). The test case we will focus on has been considered in [Nec20] and [BDE21], and represents the Helmholtz version of the classical Hadamard example for ill-posed elliptic equations. Let n ∈ N n ∈ N n inNn \in \mathbb{N}n∈N and consider the Cauchy problem
{ (19) Δ u + k 2 u = 0 in Ω := ( 0 , π ) × ( 0 , 1 ) , u ( x , 0 ) = 0 for x ∈ [ 0 , π ] , u y ( x , 0 ) = sin ⁡ ( n x ) for x ∈ [ 0 , π ] , (19) Δ u + k 2 u = 0  in  Ω := ( 0 , π ) × ( 0 , 1 ) , u ( x , 0 ) = 0  for  x ∈ [ 0 , π ] , u y ( x , 0 ) = sin ⁡ ( n x )  for  x ∈ [ 0 , π ] , {[(19)Delta u+k^(2)u=0" in "Omega:=(0","pi)xx(0","1)","],[u(x","0)=0" for "x in[0","pi]","],[u_(y)(x","0)=sin(nx)" for "x in[0","pi]","]:}\left\{\begin{align*} \Delta u+k^{2} u & =0 & & \text { in } \Omega:=(0, \pi) \times(0,1), \tag{19}\\ u(x, 0) & =0 & & \text { for } x \in[0, \pi], \\ u_{y}(x, 0) & =\sin (n x) & & \text { for } x \in[0, \pi], \end{align*}\right.{(19)Δu+k2u=0 in Ω:=(0,π)×(0,1),u(x,0)=0 for x∈[0,π],uy(x,0)=sin⁡(nx) for x∈[0,π],
whose solution for n > k n > k n > kn>kn>k is given by
(20) u ( x , y ) = 1 n 2 − k 2 sin ⁡ ( n x ) sinh ⁡ ( n 2 − k 2 y ) (20) u ( x , y ) = 1 n 2 − k 2 sin ⁡ ( n x ) sinh ⁡ n 2 − k 2 y {:(20)u(x","y)=(1)/(sqrt(n^(2)-k^(2)))sin(nx)sinh(sqrt(n^(2)-k^(2))y):}\begin{equation*} u(x, y)=\frac{1}{\sqrt{n^{2}-k^{2}}} \sin (n x) \sinh \left(\sqrt{n^{2}-k^{2}} y\right) \tag{20} \end{equation*}(20)u(x,y)=1n2−k2sin⁡(nx)sinh⁡(n2−k2y)
for n = k n = k n=kn=kn=k by u ( x , y ) = sin ⁡ ( k x ) y u ( x , y ) = sin ⁡ ( k x ) y u(x,y)=sin(kx)yu(x, y)=\sin (k x) yu(x,y)=sin⁡(kx)y, and for n < k n < k n < kn<kn<k by
u ( x , y ) = 1 k 2 − n 2 sin ⁡ ( n x ) sin ⁡ ( k 2 − n 2 y ) u ( x , y ) = 1 k 2 − n 2 sin ⁡ ( n x ) sin ⁡ k 2 − n 2 y u(x,y)=(1)/(sqrt(k^(2)-n^(2)))sin(nx)sin(sqrt(k^(2)-n^(2))y)u(x, y)=\frac{1}{\sqrt{k^{2}-n^{2}}} \sin (n x) \sin \left(\sqrt{k^{2}-n^{2}} y\right)u(x,y)=1k2−n2sin⁡(nx)sin⁡(k2−n2y)
For such Hadamard-type solutions, we will consider the interior datum g = u | ω g = u ω g=u|_(omega)g=\left.u\right|_{\omega}g=u|ω and study the convergence in the target set B B BBB for two geometric configurations of ω ω omega\omegaω and B B BBB, as discussed in Section 2, namely
(21) ω = Ω ∖ [ π 4 , 3 π 4 ] × [ 0 , 0.25 ] , B = Ω ∖ [ π 4 , 3 π 4 ] × [ 0 , 0.95 ] , (21) ω = Ω ∖ π 4 , 3 π 4 × [ 0 , 0.25 ] , B = Ω ∖ π 4 , 3 π 4 × [ 0 , 0.95 ] , {:(21)omega=Omega\\[(pi)/(4),(3pi)/(4)]xx[0","0.25]","quad B=Omega\\[(pi)/(4),(3pi)/(4)]xx[0","0.95]",":}\begin{equation*} \omega=\Omega \backslash\left[\frac{\pi}{4}, \frac{3 \pi}{4}\right] \times[0,0.25], \quad B=\Omega \backslash\left[\frac{\pi}{4}, \frac{3 \pi}{4}\right] \times[0,0.95], \tag{21} \end{equation*}(21)ω=Ω∖[π4,3π4]×[0,0.25],B=Ω∖[π4,3π4]×[0,0.95],
sketched in Figure 2b, and
(22) ω = ( π 4 , 3 π 4 ) × ( 0 , 0.5 ) , B = ( π 8 , 7 π 8 ) × ( 0 , 0.95 ) , (22) ω = π 4 , 3 π 4 × ( 0 , 0.5 ) , B = π 8 , 7 π 8 × ( 0 , 0.95 ) , {:(22)omega=((pi)/(4),(3pi)/(4))xx(0","0.5)","quad B=((pi)/(8),(7pi)/(8))xx(0","0.95)",":}\begin{equation*} \omega=\left(\frac{\pi}{4}, \frac{3 \pi}{4}\right) \times(0,0.5), \quad B=\left(\frac{\pi}{8}, \frac{7 \pi}{8}\right) \times(0,0.95), \tag{22} \end{equation*}(22)ω=(π4,3π4)×(0,0.5),B=(π8,7π8)×(0,0.95),
sketched in Figure 2a. To assess the effect of increasing the frequency, we will take n = 5 , k = 1 n = 5 , k = 1 n=5,k=1n=5, k=1n=5,k=1 and n = 11 , k = 10 n = 11 , k = 10 n=11,k=10n=11, k=10n=11,k=10, having similar values for k 2 − n 2 k 2 − n 2 sqrt(k^(2)-n^(2))\sqrt{k^{2}-n^{2}}k2−n2 in (20). Noisy data is considered by polluting the nodal values of the given restriction g g ggg with uniformly distributed values in [ − h , h ] [ − h , h ] [-h,h][-h, h][−h,h] or [ − h 2 , h 2 − h 2 , h 2 -h^(2),h^(2)-h^{2}, h^{2}−h2,h2 ], for a noise amplitude O ( h ) O ( h ) O(h)\mathcal{O}(h)O(h) or O ( h 2 ) O h 2 O(h^(2))\mathcal{O}\left(h^{2}\right)O(h2), respectively. The results shown below are obtained with the stabilization parameter γ = 10 − 3 γ = 10 − 3 gamma=10^(-3)\gamma=10^{-3}γ=10−3 in the primal stabilizer (13).
For the domains in (21), data is given around three sides of the domain Ω Ω Omega\OmegaΩ and the stability of the continuation problem is expected to be close to Lipschitz. Indeed, for approximations of order p ∈ { 1 , 2 , 3 } p ∈ { 1 , 2 , 3 } p in{1,2,3}p \in\{1,2,3\}p∈{1,2,3} we observe in Figure 3 clear convergence with rate p p ppp, indicating that the exponent α ≈ 1 α ≈ 1 alpha~~1\alpha \approx 1α≈1 in the conditional stability estimate (3). The convergence rates are in agreement with the theoretical error estimates. Also, in this case the target set B B BBB is inside the convex hull of the data set ω ω omega\omegaω, and we notice that increasing the frequency has only a small effect on the size of the errors, showing the robustness of the stability constants.
Figure 3: Convergence rates for domains in (21), good stability.
Figure 4: Noisy data, convergence rates for domains in (21), n = 5 , k = 1 n = 5 , k = 1 n=5,k=1n=5, k=1n=5,k=1.
In Figure 4 we consider noisy data for the domains in (21) with stability exponent α ≈ 1 α ≈ 1 alpha~~1\alpha \approx 1α≈1. We observe that for a noise amplitude of O ( h 2 ) O h 2 O(h^(2))\mathcal{O}\left(h^{2}\right)O(h2) the perturbations have no impact for p ∈ { 1 , 2 } p ∈ { 1 , 2 } p in{1,2}p \in\{1,2\}p∈{1,2}, while for p = 3 p = 3 p=3p=3p=3 one order of convergence is lost, as predicted by the theory. Increasing the data perturbations to O ( h ) O ( h ) O(h)\mathcal{O}(h)O(h) makes no difference for p = 1 p = 1 p=1p=1p=1, while the convergence decreases by a factor of h h hhh for p = 2 p = 2 p=2p=2p=2, and by a factor of h 2 h 2 h^(2)h^{2}h2 for p = 3 p = 3 p=3p=3p=3. This confirms the theoretical results and shows that high-order methods are more sensitive to the size of the data perturbations.
Figure 5: Convergence rates for domains in (22), poor stability.
Figure 5 shows the numerical results for the domains (22), for which the stability of the continuation problem-quantified by the exponent α ≪ 1 − α ≪ 1 − alpha≪1-\alpha \ll 1-α≪1− is expected to be poor. For k = 1 k = 1 k=1k=1k=1, the convergence rate for p = 1 p = 1 p=1p=1p=1 is close to 0.25 and increasing the frequency to k = 10 k = 10 k=10k=10k=10 has a strong effect on the numerics: the convergence rate for p = 1 p = 1 p=1p=1p=1 is now close to 0.1 and the size of the errors increase. This is due to the fact that continuation is done outside of the convex hull of the data region, and the stability constants increase very fast with the frequency. Nevertheless, in both plots we observe that the rates for quadratic and cubic approximations increase linearly with the polynomial order p p ppp, as expected from Theorem 1.
Figure 6: Noisy data, convergence rates for domains in (22), convexity condition not satisfied.
Perturbations in data for the domains (22) are considered in Figure 6. We see that for a noise amplitude of O ( h 2 ) O h 2 O(h^(2))\mathcal{O}\left(h^{2}\right)O(h2) the perturbations have no impact, while for an amplitude of O ( h ) O ( h ) O(h)\mathcal{O}(h)O(h) the method does not seem to converge for p = 3 p = 3 p=3p=3p=3. The noise effect for p = 2 p = 2 p=2p=2p=2 does not seem apparent for such coarse meshes, probably due to the large size of the stability constants. When using high-order methods for problems with poor stability and large perturbations in data, it is clear that one should consider stopping the mesh refinement before the accuracy of the approximation starts to deteriorate.

References

[BDE21] E. Burman, G. Delay, and A. Ern. A hybridized high-order method for unique continuation subject to the helmholtz equation. SIAM J. Numer. Anal., 59(5):2368-2392, 2021.
[BHL18] E. Burman, P. Hansbo, and M. G. Larson. Solving ill-posed control problems by stabilized finite element methods: an alternative to Tikhonov regularization. Inverse Problems, 34:035004, 2018.
[BM21] S. M. Berge and E. Malinnikova. On the three ball theorem for solutions of the Helmholtz equation. Complex Anal. Synerg., 7(2):1-10, 2021.
[BNO19] E. Burman, M. Nechita, and L. Oksanen. Unique continuation for the Helmholtz equation using stabilized finite element methods. J. Math. Pures Appl., 129:1-22, 2019.
[BS97] I. M. Babuška and S. A. Sauter. Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers? SIAM J. Numer. Anal., 34(6):2392-2423, 1997.
[Bur13] E. Burman. Stabilized finite element methods for nonsymmetric, noncoercive, and ill-posed problems. Part I: Elliptic equations. SIAM J. Sci. Comput., 35(6):A2752A2780, 2013.
[Bur16] E. Burman. Stabilised finite element methods for ill-posed problems with conditional stability. In Building Bridges: Connections and Challenges in Modern Approaches
to Numerical Partial Differential Equations, volume 114 of Lect. Notes Comput. Sci. Eng., pages 93-127. Springer, 2016.
[EG21] A. Ern and J.-L. Guermond. Finite elements I: Approximation and interpolation, volume 72 of Texts in Applied Mathematics. Springer Nature, 2021.
[HI04] T. Hrycak and V. Isakov. Increased stability in the continuation of solutions to the Helmholtz equation. Inverse Problems, 20(3):697-712, 2004.
[Isa17] V. Isakov. Inverse problems for partial differential equations, volume 127 of Applied Mathematical Sciences. Springer, 3rd edition, 2017.
[MS11] J. M. Melenk and S. A. Sauter. Wavenumber explicit convergence analysis for Galerkin discretizations of the Helmholtz equation. SIAM J. Numer. Anal., 49(3):1210-1243, 2011.
[Nec20] M. Nechita. Unique continuation problems and stabilised finite element methods. PhD thesis, University College London, 2020.
2022

Related Posts