Analysis of a Helmholtz preconditioning problem motivated by uncertainty quantification

This paper analyses the following question: let $\mathbf{A}_j$, $j=1,2,$ be the Galerkin matrices corresponding to finite-element discretisations of the exterior Dirichlet problem for the heterogeneous Helmholtz equations $\nabla\cdot (A_j \nabla u_j) + k^2 n_j u_j= -f$. How small must $\|A_1 -A_2\|_{L^q}$ and $\|{n_1} - {n_2}\|_{L^q}$ be (in terms of $k$-dependence) for GMRES applied to either $(\mathbf{A}_1)^{-1}\mathbf{A}_2$ or $\mathbf{A}_2(\mathbf{A}_1)^{-1}$ to converge in a $k$-independent number of iterations for arbitrarily large $k$? (In other words, for $\mathbf{A}_1$ to be a good left- or right-preconditioner for $\mathbf{A}_2$?). We prove results answering this question, give theoretical evidence for their sharpness, and give numerical experiments supporting the estimates. Our motivation for tackling this question comes from calculating quantities of interest for the Helmholtz equation with random coefficients $A$ and $n$. Such a calculation may require the solution of many deterministic Helmholtz problems, each with different $A$ and $n$, and the answer to the question above dictates to what extent a previously-calculated inverse of one of the Galerkin matrices can be used as a preconditioner for other Galerkin matrices.

as r := x 2 → ∞ (uniformly in x := x/r), where k > 0 is the wavenumber; in this paper, we are interested in the case when the wavenumber is arbitrarily large. In Equation (1.1) the obstacle D − and the coefficients A j and n j satisfy the natural conditions for the scattering problem to make sense; see Assumption 3.1 and Definition 3.2 (to present the main results as close to the beginning of the paper as possible, we postpone to §3 the precise definitions of the Helmholtz problem, the finite-element method, and GMRES).
Let A j , j = 1, 2, be the Galerkin matrices corresponding to the h-version finite-element discretisations (with decreasing mesh size h and any fixed polynomial degree p) of (1.1) truncated to D := D + ∩ D, where D is some open set containing D − (the simplest case is when D = B R := {x : x 2 < R} for some sufficiently large R). We consider the cases where either the radiation condition is realised exactly on ∂ D via the exact Dirichlet-to-Neumann (DtN) map or the DtN map is approximated on ∂ D by an impedance boundary condition (a.k.a. a first-order absorbing boundary condition). See (3.10) for the precise definitions of A j , j = 1, 2, and Definition 3.2 and Remark 3.3 for the truncation procedure.
This paper considers the following question: Q1. How small must A 1 −A 2 L q and n 1 −n 2 L q be (in terms of k-dependence) for the generalised minimum residual method (GMRES) applied to either (A 1 ) −1 A 2 or (A 2 ) −1 A 1 to converge in a k-independent number of iterations for arbitrarily large k?
Before stating results answering Q1 in §2 we describe our motivation for studying Q1.

Motivation from uncertainty quantification of the Helmholtz equation
Consider the Helmholtz equation where A(x; ω) and n(x; ω) are random fields, and ω is an element of an underlying probability space. Suppose Q(u) is a functional of the solution of (1.3) (usually called a quantity of interest).
In the forward problem of uncertainty quantification (UQ), one task is to compute E[Q(u)], and the arguably simplest way to do this uses sampling, i.e. using the approximation Q(u(ω ℓ )), (1.4) where the ω ℓ are elements of the sample space Ω. To calculate the right-hand side of (1.4), one must solve many deterministic Helmholtz problems, corresponding to different samples ω ℓ , i.e. corresponding to different realisations of the coefficients A(·, ω) and n(·, ω). Solving all these deterministic problems is a very computationally-intensive task because linear systems arising from discretisations of the Helmholtz equation are notoriously difficult to solve; this difficulty is due to the following three reasons: (a) The solutions of the homogeneous Helmholtz equation ∆u + k 2 u = 0 oscillate on a scale of 1/k, and so to approximate them accurately with piecewise-polynomial functions one needs the total number of degrees of freedom, N , to be proportional to k d as k increases.
(b) The pollution effect means that, for fixed-order finite-element methods, having N ∼ k d is still not enough to keep the relative finite-element error bounded independently of k as k increases (see the references in §4.3 for how N must depend on k).
(c) When iterative methods such as GMRES are applied to the linear system Au = f , the number of iterations grows with k. Part of the explanation for this is that the standard variational formulation of the Helmholtz equation is not coercive (i.e. it is sign-indefinite) when k is sufficiently large, and this indefiniteness is inherited by the Galerkin linear system. The design of good preconditioners for discretisations of the Helmholtz equation is therefore a very active area of research; see, e.g., the literature reviews in [35, §1.3], [27, §4].
The question therefore arises of how to compute Q(u(ω ℓ )), ℓ = 1, . . . , N , in an efficient way. To simplify the discussion, initially suppose that one calculates the LU factorisation of the Galerkin matrix for one realisation of the coefficients A and n. The answer to Q1 above makes k-explicit the extent to which this factorisation can be used as an effective preconditioner for Galerkin matrices arising from different realisations of A and n.
The benefit of "reusing" an LU factorisation in this way can be seen by recalling that direct solvers involving the LU decomposition of the linear system have a computational cost of the order O N 3/2 in 2-d [15,Section 1] and O N 2 in 3-d [15,Equation 3], with this analysis assuming sufficient grid regularity. However, if one already has computed the LU decomposition, then the cost of applying backsolves using it is much cheaper; O(N log N ) in 2-d [15,Section 1] and O N 4/3 in 3-d [15,Equation 4], hence the LU decomposition of A could be efficiently used as a preconditioner for matrices "near" A.
More generally, for any preconditioner for the Helmholtz equation, one could ask when the preconditioner corresponding to one realisation of (1.3) can be re-used for other realisations. The answer to Q1 does not answer this more-complicated scenario, but forms the first step in analysing it.
In §2.4 we discuss the implications of our answer to Q1 on such a preconditioning strategy. However, we highlight that to perform a complete k-explicit analysis of this type of preconditioning strategy applied to computing E[Q(u)], in addition to the answer to Q1, one also needs the answers to the following two questions.
Q2 How must the number of samples, N , depend on k for the error Q(u(ω ℓ )) (1.5) to be controllably small, independent of k, for arbitrarily large k?
We do not address these two questions in this paper, but remark that (i) the recent paper [28] gives results related to Q2 for a non-standard variational formulation of the heterogeneous Helmholtz equation given in [29] (based on the formulation for the homogeneous Helmholtz equation introduced in [48]), and (ii) the answer to Q3 follows from the specific structure of the randomness assumed in the coefficients A and n.

Novelty of the main results
To our knowledge, this paper is the first time that preconditioning a discretisation of the heterogeneous Helmholtz equation with a discretisation corresponding to a Helmholtz problem with "nearby" coefficients has been analysed. This paper is part of a growing body of work on the analysis and numerical analysis of the heterogeneous Helmholtz equation, with the analysis explicit in both the wavenumber and the coefficients [12,25,31,32,33,43,49,53,62] Furthermore, we believe that the idea of preconditioning with "nearby" coefficients will be relevant to other PDEs with random coefficients (this technique can then then be seen as a generalisation of the mean-based preconditioning discussed in §2.4). We expect that the ideas in our analysis -in particular Lemmas 5.7 and 5.8 which combine bounds on the continuous problem and on the finite-element error to prove bounds on preconditioned matrices -would then be useful in analysing such strategies.
Outline of the paper. §2 states the main results and gives supporting numerical experiments.
§3 gives the precise definitions of the Helmholtz boundary value problem, its finite-element solution, and weighted GMRES. §4 give the precise definitions of the conditions under which the main results hold. §5 and §6 prove the main results.

Statement of the main results
We state in §2.1 the results obtained when the differences A 1 − A 2 and n 1 − n 2 are measured in the L ∞ norm (Theorems 2.2 and 2.3). We discuss in §2.2 the sharpness of these results. We then state in §2.3 generalisations of Theorems 2.2 and 2.3 where the differences are measured in the L q norm for 2 < q ≤ ∞ (Theorems 2.9 and 2.10). We discuss in §2.4 the implications of these results for calculating quantities of interest of the Helmholtz equation with random coefficients (following on from the discussion in §1.2).
We consider both standard GMRES (which uses the Euclidean norm · 2 on vectors) and weighted GMRES, a short description of which is given in §3. In weighted GMRES, we use the weighted vector norms · D k and · D −1 with S I is the standard stiffness matrix for the finite-element discretisation of the Laplacian, and M 1 is the mass matrix -see (3.9) for a precise definition. The key point is that, for a finite-element where v is the vector with ith entry v i ; i.e., · D k is the norm on the finite-element space induced by the weighted H 1 norm · H 1 k (D) in which the PDE analysis of the Helmholtz equation naturally takes place. (Results about convergence of domain-decomposition methods in the norms (2.1) were recently obtained for the Helmholtz equation in [34,35,31] and for the time-harmonic Maxwell equations in [7].) We further define A L ∞ (D;op) := ess sup x∈D A(x) 2 , where here · 2 denotes the spectral/operator norm on matrices, induced by the Euclidean norm · 2 on vectors.
To make the statements of the main results more concise, we make the following definition for weighted/standard GMRES. Definition 2.1 We say that GMRES applied to Cu = f , with zero starting guess and in a norm · on C n , converges in a k-independent number of iterations if, given ε > 0 and k 0 > 0, there exists C 1 (ε, k 0 ) > 0 and C 2 (k 0 ) > 0, both independent of k, h, and p, such that In this paper we only consider GMRES in one of the three norms · , · D k , and · D −1 k (with the last two defined in (2.1)).

Results on Q1 involving
These results are proved under Conditions 4.1 and 4.2. These conditions can be informally stated, respectively, as • the obstacle D − and the coefficients A j and n j are such that u j exists, is unique, and the problem is nontrapping (in the sense described in §4.2), and • the meshsize h and polynomial degree p in the finite-element method are chosen to depend on k to ensure that the finite-element approximation to the solution of the problem with coefficients A j and n j exists, is unique, and has bounded error in the H 1 k -norm (defined in (2.12)) as k → ∞.
Given k 0 > 0, there exist C 1 , C 2 > 0 independent of h and k (but dependent on D − , A 1 , n 1 , p, and k 0 ) such that if for all k ≥ k 0 , then both weighted GMRES working in the vector norm · D k (and the associated inner product) applied to and weighted GMRES working in the vector norm · (D k ) −1 (and the associated inner product) applied to with initial guess the zero vector converge in a k-independent number of iterations (in the sense of Definition 2.1) for all k ≥ k 0 .
Given k 0 > 0, let C 1 and C 2 be as in Theorem 2.2, and let s + and m ± be as in Lemma 5.1 (note that all of these constants are independent of k and h). Then if for all k ≥ k 0 , then standard GMRES (working in the Euclidean norm and inner product) applied to either of the equations (2.4) or (2.5) with initial guess the zero vector converge in a k-independent number of iterations (in the sense of Definition 2.1) for all k ≥ k 0 .
We make the following remarks.
FEM2 given in Condition 4.2, -n 1max , n 1min , and A 1min given in (3.1) and (3.2) (with n replaced by n 1 and A replaced by A 1 ), (b) The constant 1/2 on the right-hand sides of (2.3) and (2.6) can be replaced by any number < 1 and the overall result that GMRES converges in a k-independent number of iterations still holds, although the actual number of GMRES iterations depends on this number (but not on k). Lemma 2.4 (Bound on standard GMRES from condition on weighted GMRES) Given k 0 > 0, assume that (2.3) holds (so that weighted GMRES with initial guess the zero vector converges in a k-independent number of iterations). If h = Ck −1−δ with δ ≥ 0, then there exists C > 0, depending only on s + , m ± , C, δ, and k 0 , such that if standard GMRES is applied to the linear system (2.4) with initial guess the zero vector, and This result is not of the same form as Theorems 2.2 and 2.3, which give sufficient conditions for weighted/standard GMRES to converge in a k-independent number of iterations (in the sense of Definition 2.1). To better understand Lemma 2.4, we observe that the proof of Theorem 2.2 (see §5.3) shows that, given k 0 > 0, if (2.3) holds, and  .3)).
Remark 2.5 (Sufficient conditions for the relative residual to be controllably small) Theorems 2.2 and 2.3 give sufficient conditions for GMRES to converge in a k-independent number of iterations in the sense of Definition 2.1. These conditions are also sufficient for the GMRES relative residual to be controllably small, independent of k, in the sense that, given ε > 0 and k 0 > 0, there exists C 1 (ε, k 0 ) > 0 and C 2 (k 0 ) > 0, both independent of k, h, and p, such that if m ≥ C 1 then r m r 0 ≤ C 2 ε, for all k ≥ k 0 . is not only an answer to Q1 but is also the natural answer to an analogue of Q1 at the level of PDEs, namely Q1 ′ Let u j , j = 1, 2, be the Helmholtz solution with coefficients A j and n j ; how small must A 1 − A 2 L ∞ (D;op) and n 1 − n 2 L ∞ (D;R) be (in terms of k) for the relative error in approximating u 2 by u 1 to be bounded independently of k for arbitrarily-large k?
(The claim that Q1 ′ is an analogous question to Q1 at the level of PDEs is made precise in Remark 6.2.) Furthermore Lemma 2.7 shows that the condition "k n 1 − n 2 L ∞ (D;R) sufficiently small" is the provably-sharp answer to Q1 ′ when A 1 = A 2 = I.
Before stating these PDE results, we recall from (2.2) the definition of the weighted for v ∈ H 1 0,D (D), (2.12) where the space H 1 0,D (D), defined by (3.6), is the natural space containing the solution of the exterior Dirichlet problem. Finally, let (H 1 0,D (D)) * denote the dual space of H 1 0,D (D).
Theorem 2.6 (Answer to Q1 ′ ) Given F ∈ (H 1 0,D (D)) * and coefficients A j , n j , j = 1, 2, let u j , j = 1, 2, be the solutions of the variational problem (3.5) with coefficients A j , n j . Assume that D, A 1 , and n 1 satisfy Condition 4.1 and that D, A 2 , and n 2 are such that u 2 exists.
Then, given k 0 > 0, there exists C 3 > 0, independent of k and F and given explicitly in terms of D − , A 1 , and n 1 in (6.2), such that for all k ≥ k 0 .
In stating the next result (and elsewhere in the paper), for a, b > 0 we write a b when a ≤ Cb for some C > 0, independent of k and h. We also write a ∼ b if a b and b a.
there exists an annulus between D − and ∂ D).
Then there exist particular choices of F, n 1 , and n 2 (with n 1 = n 2 , and both continuous functions on D + ) such that the corresponding solutions u 1 and u 2 of the variational problem (3.5) with A 1 = A 2 = I exist, are unique, and, given k 0 > 0, for all k ≥ k 0 .
Remark 2.8 (Physical interpretation of the condition (2.11)) Recall that the wavelength of the wave u (at least when A = I and n = 1) is 2π/k. Since this is the natural length scale associated with u, one expects from physical considerations that perturbations of magnitude ≤ c/k, for c > 0 sufficiently small, are 'unseen' by the PDE or numerical method. This is indeed what we see in Theorems 2.2, 2.3, and 2.6: perturbations of this form (i.e. A j and n j satisfying (2.11)) give bounded relative difference (for Q1 ′ ) and bounded GMRES iterations for the nearby-preconditioned linear system (for Q1).

Numerical experiments to test the sharpness of Theorems 2.2 and 2.3
We now present experiments testing the sharpness of the condition op) sufficiently small" as a sufficient condition for the kindependent convergence of standard GMRES, we recall from the discussion below Lemma 2.4 that we expect (2.15) to be sufficient. While we do not present any experiments with weighted GMRES in this paper, we highlight that in [34] the number of iterations for standard GMRES and weighted GMRES with the weight D k were found to be almost identical when applied to preconditioned linear systems arising from the Helmholtz FEM.
Description of the set-up. We consider the solution of the variational problem (3.5) with D − = ∅, D + = [0, 1] 2 , and T = ik. These choices correspond to the interior impedance problem on the unit square in 2-d. We set A 1 = I, n 1 = 1, f = 0, and define g I so that the exact solution u 1 is a plane wave incident from the bottom left passing through the homogeneous medium given by coefficients A 1 and n 1 . We solve this variational problem via the finite-element method defined in (3.8); we choose the family of finite-dimensional subspaces (V h,p ) h>0 to be first-order (i.e. p = 1) continuous finite elements on regular grids with h = k −3/2 ; this choice of h means that the finite-element approximation u h is uniformly accurate as k → ∞ (see §4.3 and the references therein). All finite-element calculations were carried out using the Firedrake software library [59,44], which itself uses the PETSc [3,13], Chaco [38], and MUMPS [1,2] software packages. The code used in this paper can be found at https://github.com/orpembery/running-nbpc. The exact versions of the code used in this paper, and the corresponding computational data can be found in the releases [54,55,56].
Description of the four different A 2 and n 2 considered.
1. A 2 = I and 100 realisations of randomly-chosen n 2 that are piecewise constant with respect to a 10 × 10 square grid on D (experiments in Figure 1).
2. n 2 = 1 and 100 realisations of randomly-chosen A 2 that are piecewise constant with respect to a 10 × 10 square grid on D (experiments in Figure 2).

3.
A 2 = I and deterministic n 2 that are piecewise constant in a checkerboard pattern on a 10 × 10 square grid on D (experiments in Figure 3). 4. A 2 = I and deterministic n 2 that are piecewise constant in a checkerboard pattern on a 2 × 2 square grid on D (experiments in Figure 4).
Regarding 1: Each instance of n 2 is piecewise constant with respect to a 10 × 10 square grid on D : Regarding 2: On each element of the grid we set where a, c ∈ Unif(0, α) and α as given in (2.16). For each draw of a, c we choose b ∈ Unif(0, δ), where δ = min{α, (1 + a)(1 + c)}. It follows that each draw of A 2 is positive definite almost surely. Since the elements of A 1 − A 2 = I − A 2 are all bounded above by α, direct computation shows that Regarding 3: As in 1., n 2 is piecewise constant with respect to a 10 × 10 square grid on D, but now n 2 alternates between the values 1 − α and 1 + α (with α as in (2.16)) in a checkerboard pattern. In this case, Regarding 4: This set up is exactly the same as in 3., except that now the square grid on D is 2 × 2. The motivation for this checkerboard is that now there are just under 16 wavelengths in each piece of the board (of width 0.5) when k = 100 (i.e., at the largest k considered), whereas in the 10 × 10 checkerboard there is only just over one wavelength in each piece of the board at k = 100.
In all four cases, to obtain the action of (A 1 ) −1 , we calculate the exact LU decomposition of A 1 . We then solve the linear system (2.4) using standard GMRES (i.e. with residual minimisation in the Euclidean norm) and zero initial guess, and record the number of GMRES iterations taken to achieve a relative residual of 10 −5 . In the first and second experiments, we consider k between 20 and 100, and in the third and fourth experiments we consider k between 20 and 200.
Grey markers indicate GMRES reached 500 iterations without convergence.
Description and discussion of results. In all cases we see the number of GMRES iterations being bounded (over the ranges considered) for β ∈ {0.8, 0.9, 1.0}, and growing for β ∈ {0, . . . , 0.5}, with the growth worse for the two deterministic cases in This behaviour is better than expected from the condition (2.15), which covers β = 1. Because of the results about Q1 ′ at the PDE level in Theorem 2.6 and Lemma 2.7, we expect that the number of iterations will eventually grow for β < 1 for k sufficiently large.

Results involving
The results in Theorems 2.2 and 2.3 involving A 1 − A 2 L ∞ and n 1 − n 2 L ∞ give an answer to Q1 that fits well with both PDE considerations (via Theorem 2.6) and physical intuition (via Remark 2.8), and is consistent with the numerical experiments for moderate k in §2.2.2. However, the following example shows the disadvantage of measuring the difference in coefficients in the L ∞ norm. Let D = {x : x 2 = 2}, A 1 = A 2 = I, and    for some 0 < α < 1/2. Then n 1 − n 2 L ∞ (D;R) = 1/2 for all α, but one may expect that GMRES applied to (A 1 ) −1 A 2 would converge in a k-independent number of iterations for small α. In this subsection, we therefore state analogues of Theorems 2.2 and 2.3, with the differences in n 1 − n 2 and A 1 − A 2 measured in weaker norms than the L ∞ norm.
Statement of the results. To state these results, we need the definition that    Similar to in §2.1, Theorem 2.10 gives results for standard GMRES (i.e., using the vector Euclidean norm · 2 ), whereas Theorem 2.9 gives results in the weighted vector norms · D k and · D −1 k . Theorem 2.9 (Answer to Q1: k-independent weighted GMRES iterations) Assume that • D − , A j , and n j , j = 1, 2, satisfy Condition 4.1, and • D − , A j , n j , j = 1, 2, h and p, satisfy Condition 4.2.
Given k 0 > 0 and q > 2, there exist constants C 1 and C 2 independent of h and k (but dependent on D − , A 1 , n 1 , p, and k 0 ) such that if for all k ≥ k 0 , then both weighted GMRES working in · D k (and the associated inner product) applied to and weighted GMRES working in · (D k ) −1 (and the associated inner product) applied to with initial guess the zero vector converge in a k-independent number of iterations (in the sense of Definition 2.1) for all k ≥ k 0 .
Given k 0 > 0, and q > 2, let C 1 and C 2 be as in Theorem 2.9, and let s + > 0 and m ± > 0 be as in Lemma 5.1 (note that all these constants are independent of k, h, and p). Then if for all k ≥ k 0 , then standard GMRES (working in the Euclidean norm and inner product) applied to either (2.4) or (2.5) with initial guess the zero vector converges in a k-independent number of iterations (in the sense of Definition 2.1) for all k ≥ k 0 .
We make the following remarks.
• The constants C 1 and C 2 are given explicitly in (5.6) and (5.7) in terms of q and -C -C inv,s given in (5.1).
• Theorems 2.2 and 2.3 can be viewed as special cases of Theorems 2.9 and 2.10, respectively, with q formally set to ∞.
• In the conditions (2.20) and (2.21) there is a trade-off between the norm used to measure n 1 − n 2 and the restriction on the magnitude of this norm. Indeed, as q ց 2, n 1 − n 2 is measured in a weaker norm, but the condition on n 1 − n 2 L q in (2.20) and (2.21) becomes more restrictive since h −d/q becomes larger. Conversely, as q ր ∞, n 1 − n 2 is measured in a stronger norm, but the condition on n 1 − n 2 L q in (2.20) and (2.21) becomes less restrictive since h −d/q becomes smaller.
• The condition (2.3) contains no powers of h, but the analogous condition (2.20) does. The reason is that, in the proofs of Theorems 2.2 and 2.3, we use the inequality ( where v h ∈ V h,p , to "pull out" the L ∞ norm of n 1 − n 2 . In Theorems 2.9 and 2.10 we want to "pull out" the L q norm of n 1 − n 2 , so we use the inequality (n 1 − n 2 )v h L 2 (D) ≤ n 1 − n 2 L q (D;R) v h L s (D) where 1/2 = 1/q + 1/s (see (5.23)). We then need to convert the L s norm of v h into an L 2 norm (essentially because the PDE (1.1) is posed in the L 2 -based Sobolev space H 1 0,D (D)), and we do this using the inverse inequality (5.1), introducing powers of h. Similar considerations hold for A 1 − A 2 .

Implications of these results for uncertainty quantification for the Helmholtz equation
The results of Theorems 2.2, 2.3, 2.9, and 2.10 and the numerical results in §2.2.2 show that A 1 − A 2 and n 1 − n 2 must decrease as k increases for A 1 to be a good preconditioner for A 2 . This suggests that the idea of "reusing" preconditioners in calculating Q(u(ω ℓ )), ℓ = 1, . . There is therefore a trade-off between how quickly N increases with k and how quickly the region in coefficient space for which a preconditioner is effective decreases with k.
In the rest of this subsection we give some more details about the experiments in [53, §4.6.3] about "reusing" preconditioners in calculating Q(u(ω ℓ )). Before doing this, however, we highlight that this strategy is related to mean-based preconditioning, where a single preconditioner is calculated corresponding to the mean of the random coefficient. The initial computational work on mean-based preconditioning for the equation −∇·(A∇u) = f was carried out by [30], [52], and [42], with theory following from [58] and [20]. In the Helmholtz context, mean-based preconditioning has been used in [41] and [65] (see the literature review in [53, §4.7] for more discussion).
Preliminary experiments in [53, §4.6.3] show that, at least for moderate k, the idea of "reusing nearby preconditioners" can give a significant speedup in the implementation of Quasi-Monte-Carlo (QMC) methods for Helmholtz problems. One of the experiments in [53, §4.6.3] concerned the UQ of various quantities of interest for the solution of (1.1) on the unit square D = (0, 1) 2 with an impedance boundary condition, and with coefficients taken to be A = I and n(x, y) = 1 + s j=1 y j ψ j (x), with x ∈ D, y j ∈ Unif(−1/2, 1/2) i.i.d. and ψ j (x) = j −2 cos(jπx 1 /4) cos((j + 1)πx 2 /4). For the UQ, a randomly shifted QMC rule taken from [51] was used.
The first part of the experiment uses extrapolation to determine how the number of QMC points N = N (k) need to be chosen to ensure that the error in the expected value of each of the quantities of interest is bounded with respect to k. These experiments suggest that for k ∈ [10,60], N (k) needs to grow somewhere between linearly and quadratically in k to keep the error bounded; thus providing an empirical answer in this case to Q2 in §1.2 (see (1.5)).
Since N (k) instances of the boundary value problem have to be solved for each k, [53, §4.6.2] looks at speeding this computation up by computing a set of 'preconditioning points' on a (coarse) sublattice of the QMC points. By the theory above, for the sublattice to be guaranteed useful, for each point y on it, there should be another coarse lattice point y ′ satisfying approximately n(·, y) − n(·, y ′ ) L ∞ (D) k −1 . In practice, it is hard to compute such a lattice (see Q3 in §1.2), so in [53, §4.6.3] a coarse sublattice of N c = N c (k) points is computed that instead equidistributes the right-hand side of the upper bound with the value of N c again determined by extrapolation from small values of k. Results for this procedure are given in [53, Table 4.5] and these show that the ratio N c /N is never bigger than 0.016 for k ∈ [10,60] in dimension s = 10. Therefore, more than 98% of the required linear systems can be solved by preconditioned GMRES. Moreover, the average number of GMRES iterations needed for the solution of all the required linear systems varies between 6.46 to 7.41 as k increases from 10 to 60, showing the good quality and "k−robustness" of the set of preconditioning points. In the experiments in [53], the restarted version of GMRES (the default in Firedrake) was used with restarts after 30 iterations. This preconditioning strategy for Helmholtz UQ will be investigated in more detail elsewhere.
3 Definitions of the Helmholtz boundary value problem, its finite-element solution, and weighted GMRES Function-space notation: Let SPD be the set of all symmetric-positive-definite matrices in  (ii) n ∈ L ∞ (D; R) is such that supp(1 − n) ⊂⊂ D and there exist 0 < n min ≤ n max < ∞ such that n min ≤ n(x) ≤ n max for almost every x ∈ D,  where ·, · ΓI denotes the duality pairing on Γ I that is linear in the first argument and antilinear in the second argument. When A = A j and n = n j , for either j = 1 or j = 2, we denote the sesquilinear form a(·, ·) by a j (·, ·). A simple approximation to DtN is the operator of multiplication by ik. The boundary condition (3.4) then becomes the impedance boundary condition, which is also known as a first-order absorbing boundary condition (see, e.g, [5], [40, §3.3]). See [24,Equation 5.3] for analysis of the approximation DtN ≈ ik in the case when D is a ball, and [23] for upper and lower bounds on the relative error between the solution of the TEDP and the exterior Dirichlet problem for more general D.

Remark 3.4 (Existence and uniqueness)
With either of the choices T = ik and T = DtN, a(·, ·) satisfies a Gårding inequality (for T = ik this is straightforward; for T = DtN see, e.g., [50, Theorem 2.6.6]). Therefore, once uniqueness is established the Fredholm alternative then gives existence. The unique continuation principle (UCP) holds (and therefore gives uniqueness) when n ∈ L ∞ and (i) d = 2 and A ∈ L ∞ , (ii) d = 3 and A ∈ C 0,1 ; see the references in [32, §1], [33, §2] (where the latter paper applies these results specifically to the Helmholtz equation). The UCP does not hold in general when d = 3 and A is rougher than Lipschitz (see [22]), but uniqueness is proved for a particular class of A ∈ L ∞ in [32, Theorem 2.7]. For such D it is not possible to fit ∂D exactly with simplicial elements, and some analysis of non-conforming error is therefore necessary; since this is standard (see, e.g., [8, Chapter 10]), we ignore this issue here. Definition 3.5 (Finite-element approximation of the solution of the TEDP) With a(·, ·) and F given in (3.7), find u h ∈ V h,p such that (3.8) We say that u h ∈ V h,p is the finite-element approximation of u.
Conditions on k, h, and p under which u h exists and is unique are discussed in §4.3. (Regarding the notation u h : although the finite-element approximation depends on both h and p, since we consider decreasing h with p fixed in this paper, for brevity the notation only reflects the h-dependence.) The matrices associated with our finite-element discretisation are defined as follows. Let {φ i , i = 1, . . . , N } be a nodal basis for V h,p with each φ i real-valued. Let with S I and M 1 the standard stiffness and mass matrices defined by (3.9) above (thus D k is independent of the coefficients A and n). We then define the weighted norm · D k by (2.1).
Weighted GMRES. We now give the set-up for weighted GMRES, introduced in [21] (with GMRES introduced in [60]). Consider the abstract linear system Cx = d in C N , where C ∈ C N ×N is invertible. Let x 0 be the initial guess, and define the initial residual r 0 := Cx 0 − d and the standard Krylov spaces: K m (C, r 0 ) := span C j r 0 : j = 0, . . . , m − 1 .
Analogously to the definition of · D k above, let (·, ·) D denote the inner product on C n induced by some Hermitian positive-definite matrix D, i.e. (v, w) D := (Dv, w) 2 , and let · D be the induced norm. For m ≥ 1, define the mth GMRES iterate x m to be the unique element of K m such that its residual r m := Cx m − d satisfies the minimal residual property: The algorithm stops when r m 2 / r 0 2 ≤ ε for some prescribed tolerance ε. Observe that when D = I this is the standard GMRES algorithm. We also note that in general, weighted GMRES requires the use of weighted Arnoldi process, also introduced in [21], see also the alternative implementations of the weighted Arnoldi process in [37].  bound . Condition 4.2 is a condition on h and p in the finite-dimensional subspaces V h,p and implicitly also a condition on D, A, and n in the TEDP (since D, A, and n must satisfy certain properties for the finite-element approximation to exist). Condition 4.2 involves the constants C FEM1 and C FEM2 ; when D, A j , and n j satisfy this, we write the corresponding constants as C u, exists, is unique, and, given k 0 > 0, u satisfies the bound where C bound is independent of k, but dependent on A, n, D, and k 0 .
To state Condition 4.2 we need the definition of the norm on (H 1 k (D)) * , .

(4.3)
Condition 4.2 (k-independent accuracy of the FE solution for a(·, ·)) Given k 0 > 0, 1. D, A, n, h and p are such that, if f = n j α j φ j for some α j ∈ C and n ∈ L ∞ (D; R) (i.e. f is an arbitrary element of V h,p multiplied by n), then for all k ≥ k 0 , • the solution u h of (3.8) with F (v) given by (4.1) exists and is unique, and • the error bound holds, where C FEM1 is independent of k and h, but dependent on A, n, D, k 0 , and p.
2. D, A, n, h and p are such that, if F (v) = ( A∇ f , ∇v) L 2 (D) , where A ∈ L ∞ (D; SPD) and f := j α j φ j with α j ∈ C (i.e. f is an arbitrary element of V h,p ), then for all k ≥ k 0 , • the solution u h of (3.8) exists and is unique, and • the error bound holds, where C FEM2 is independent of k and h, but dependent on A, n, D, k 0 , and p.

When does Condition 4.1 hold?
The summary is the following: 1. Condition 4.1 holds when D, A, and n are such that the problem is nontrapping, in the sense that all rays starting in D hit Γ I after some uniform time.
2. The definition of nontrapping is subtle in the case when either D − = ∅ or A and n are discontinuous.
3. Expressions for C bound that are explicit in A 1 and n 1 have recently been obtained [25], [32], [49]; essentially C bound is a constant multiple of the length of the longest ray in D.

Regarding 1:
In the case when D − = ∅, and A and n are C ∞ , the rays of the Helmholtz equation ∇ · (A∇u) + k 2 nu = −f are defined as the projections in x of the solutions (x(s), ξ(s)) ∈ R d × R d of the Hamiltonian system where the Hamiltonian H(x, ξ) given by (observe that H is the semiclassical principal symbol of the Helmholtz equation; see, e.g., [32, §7] for an introductory discussion on this).
When T = DtN, the propagation of singularities results of [16,§VI] (see also, e.g., [39,Chapter 24], [68, §12.3]) combined with the parametrix argument of [64] then imply that Condition 4.1 is satisfied. These arguments, however, do not give an explicit expression for C bound in terms of (properties of) A and n.
When T = ik the results [4, Theorems 5.5 and 5.6 and Proposition 5.3] imply that Condition 4.1 is satisfied, but do not give an explicit expression for C bound in terms of (properties of) A and n.
Regarding 2: When A and n are C ∞ , but D − = ∅, defining how the rays interact with the boundary Γ is subtle, and requires the notion of the Melrose-Sjöstrand generalized-bicharacteristic flow [39, §24.3], [46], [47]. However, under a suitable nontrapping hypothesis [47,Definition 7.20], the results of [64] then imply that Condition 4.1 is satisfied when T = DtN. When T = ik, the analogous result is given in [4] (with the specific case of a nontrapping Dirichlet obstacle covered by [4,Equation 5.2]). When A and n are discontinuous, the signs of the jumps dictate whether the problem is nontrapping or trapping; see [57], [10], and [49].
Regarding 3: When T = DtN, D, A, and n are such that the problem is nontrapping, both n 1 and A 1 are globally C 1,1 and C ∞ in a neighbourhood of D − , and no rays are tangent to Γ D to infinite order then [25, Theorem 2 and Equation 6.32] proves that Condition 4.1 holds with where L(A, n, D) is the length of the longest ray in a fixed neighbourhood of D; this result also holds when D − = ∅ and both n 1 and A 1 are globally C 1,1 . (When A 1 = I and n 1 = 1 the rays are straight lines, and so no rays being tangent to Γ D to infinite order means that the boundary is never flat to infinite order.) When (i) T = DtN and D − is star-shaped with respect to a point, or (ii) T = ik and D − and D are both star-shaped with respect to the same point, then expressions for C bound are given in [32] when A and n satisfy certain monotonicity conditions in the radial direction. These monotonicity conditions allow for A and n to be discontinuous (with nontrapping jumps), and the results of [32] thus recover the earlier results of [49] for piecewise constant A and n when T = DtN.

When does Condition 4.2 hold?
Part 1 of Condition 4.2: When D, A, and n satisfy both Condition 4.1 and additional smoothness requirements, Part 1 of Condition 4.2 holds (at least when T = ik) when k(hk) 2p ≤ C 1 for some C 1 depending on C bound .
Indeed, by [53,Theorem 2.39], when D, A, and n satisfy Condition 4.1, T = ik, ∂D ∈ C p,1 , A ∈ C p−1,1 , and n ∈ H max{p−1,⌈d/2⌉+1} there exist C 2 , C 3 > 0, depending on A and n (but not on this property is also ensured by requiring n ∈ C min{p−2,0},1 (D).) Furthermore, the arguments in [43] show that, when D, A, and n satisfy Condition 4.1, T = DtN, p = 1, ∂D ∈ C 0,1 , A ∈ C 0,1 , and n ∈ L ∞ , there exist C 2 , C 3 > 0, depending on A and n (but not on C bound ), such that if (4.6) holds, then Part 1 of Condition 4.1 holds with C FEM1 = C 3 C bound . Indeed, while [43, Theorem 1.5] obtains a bound on the relative error of the finiteelement approximation (rather than the error in terms of the data as in (4.4)), using the bound u H 2 k C bound f L 2 in [43, Equation 2.5] gives (4.4), with this bound on the H 2 norm following from (4.2) by elliptic regularity; see, e.g., [32, Remark 2.14] and the references therein.
Part 2 of Condition 4.2: When D, A, and n satisfy both Condition 4.1 and additional smoothness requirements, Part 2 of Condition 4.2 holds (at least when T = ik) when k(hk) p ≤ C 4 for some C 4 depending on C bound . Indeed, by [12, Theorem 2.15 and Remark 2.2], when D, A, and n satisfy Condition 4.1, T = ik, ∂D ∈ C p,1 , A ∈ C p−1,1 , and n ∈ C p−1,1 then there exists C 4 , C 5 > 0, depending on A, n (but not i.e. the finite-element approximation is quasi-optimal. Lemma 5.10 shows that, if Condition 4.1 holds, then u H 1 k F (H 1 k (D)) * , and then, setting v h = 0 in (4.7), we see that (4.5) holds.

Lemma 5.1 (Norm equivalences of FE functions)
There exist m ± > 0 and s + > 0, independent of h (but dependent on p), such that

2)
and Sketch proof The inequalities (5.2) can be proved by direct computation, and then (5. Remark 5.2 (Relaxing the assumption of quasi-uniformity) We assume that (T h ) h>0 is a quasi-uniform family of meshes so that we can use the inverse inequalities (5.1) and (5.4). Nevertheless, (i) We see below that the proof of Theorem 2.2 does not use either (5.1) or (5.4), and so Theorem 2.2 holds without the quasi-uniformity assumption.
(ii) We expect that analogues of some of the results in Theorems 2.3, 2.9, and 2.10 can be obtained for shape-regular meshes, following the techniques in the proofs [26, Section 3.4 and 4.1.2] (see Remark 5.11 for discussion on how the results of [26] are related to those in the present paper). Indeed, [26] contains bounds on preconditioned mass matrices, analogous to those in Lemma 5.7. From the way Lemma 5.7 is used in the proofs of Theorems 2.3, 2.9, and 2.10, we expect that analogues of these results with A 1 = A 2 hold in the case of shape-regular meshes. However, [26] does not contain any results on preconditioned stiffness matrices analogous to Lemma 5.8. Therefore, from the way Lemma 5.8 is used in the proofs of Theorems 2.3, 2.9, and 2.10, it is unclear at this point whether analogues of these results with A 1 = A 2 can be obtained in the case of shape-regular meshes.

Proofs of Theorems 2.9 and 2.10
In this subsection we prove Theorems 2.9 and 2.10; in the next subsection we indicate how the proofs can be modified slightly to prove Theorems 2.2 and 2.3.
The main ingredient of the proofs of Theorems 2.9 and 2.10 is the following lemma, which we prove in §5.2.2. Let m ± and s + be given as in Lemma 5.1. Given k 0 > 0 and q > 2, let where bound n 1max , (5.6) and C inv,s is defined by (5.1) with 1/s = 1/2 − 1/q. Then, for all k ≥ k 0 ,  where c := 2 √ 2/3 < 1. Therefore, given ε > 0, if m ≥ log(1/ε)/ log(1/c), then r m D ≤ ε r 0 D . By definition of r m and the fact that x 0 = 0, Therefore, Proof of Theorem 2.10 This follows from Lemma 5.3 in an analogous way to the proof of Theorem 2.9, except with D = I and using (5.9) instead of (5.8).
Remark 5.6 (The improvement of the Elman estimate (5.11) in [6]) A stronger result than (5.11) is given for standard (unweighted) GMRES in [6, Theorem 2.1], and then converted to a result about weighted GMRES in [7,Theorem 5.3]. In the bound of [6, Theorem 2.1], the convergence factor sin β is replaced by a function of β strictly less than sin β for β ∈ (0, π/2). Using this stronger result, however, does not improve the k-dependence in the bounds of Theorems 2.9 or 2.10.

Proof of Lemma 5.3
To prove Lemma 5.3 we need the following two lemmas.
Lemma 5.7 (Bounds on (A 1 ) −1 M n and M n (A 1 ) −1 ) Assume that D, A 1 , and n 1 satisfy Condition 4.1 and that D, A 1 , n 1 , h, and p satisfy Part 1 of Condition 4.2. Let m ± be as in Lemma 5.1 and let C 1 and C 2 be given by (5.5). Then, for all n ∈ L ∞ (D; R), q > 2, and k ≥ k 0 , and Lemma 5.8 (Bounds on (A 1 ) −1 S A and S A (A 1 ) −1 ) Assume that D, A 1 , and n 1 satisfy Condition 4.1 and that D, A 1 , n 1 , h, and p satisfy Part 2 of Condition 4.2. Let m ± and s + be as in Lemma 5.1 and let C 1 and C 2 be given by (5.5). Then for all A ∈ L ∞ (D, SPD), q > 2, and k ≥ k 0 , and Proof of Lemma 5.3 using Lemmas 5.7 and 5.8 Using the definition of the matrices A j , S A , and M n in (3.10) and (3.9), we have (5.18) and similarly The bounds in ( These proofs rely on the fact that, from the definitions of · D k and · (D k ) −1 in (2.1), for where C † is the conjugate transpose of C (i.e. the adjoint with respect to (·, ·) 2 ). This result follows from the fact that, for all v j ∈ C N , where w j := (D k ) −1 v j , j = 1, 2,. Therefore, since M n and S A are real, symmetric matrices The matrix (A 1 ) † is the Galerkin matrix corresponding to the adjoint variational problem to (3.5).
The key point is that if Conditions 4.1 and 4.2 are satisfied for the variational problem (3.5), then one can also show that they are satisfied for the adjoint problem. Therefore, the bound in (5.14) on M n (A 1 ) −1 (D k ) −1 follows from the bound on (A 1 ) −1 M n D k , and the bound in (5.16) on (A 1 ) −1 S A D k also holds for ((A 1 ) † ) −1 S A D k . This type of duality argument was first used in [26] and [34].
Proof of the bounds on (A 1 ) −1 M n in Lemma 5.7 We first concentrate on proving (5.14). Given f ∈ C N and n ∈ L ∞ (D; R), we create a variational problem whose Galerkin discretisation leads to the equation A 1 u = M n f . Indeed, let f := j f j φ j ∈ H 1 0,D (D). Define u to be the solution of the variational problem and let u h be the solution of the finite-element approximation of (5.20), i.e., Hölder's inequality implies that, if q, s > 2 such that 1/2 = 1/q + 1/s, then Since f ∈ V h,p , we can use the inverse inequality (5.1) in (5.23), and combining the result with (5.22) we obtain that . The bound on (A 1 ) −1 M n D k in (5.14) then follows from the definition (2.1) of · D k and the definition (5.5) of C 2 .
To prove the bound on (A 1 ) −1 M n 2 in (5.15), we use the consequences of (5.2) Observe that the definitions (4.3) and (2.12) of the norms · (H 1 k (D)) * and · H 1 k (D) , respectively, and the Cauchy-Schwarz inequality imply that . With · L q (D,op) defined by (2.19), Hölder's inequality implies that, with q, s > 2 such that 1/2 = 1/q + 1/s, Combining this with the inverse inequality (5.1) we then have that Let u h be the solution of the finite-element approximation of (5.26), i.e., and let u be the vector of nodal values of u h . The definition of f then implies that (5.28) is equivalent to Similar to the proof of the bounds on (A 1 ) −1 M n in Lemma 5.7, using the triangle inequality, the bound (4.5) from Condition 4.2, the bound (5.25) from Lemma 5.10, the bound (5.27), and the definition (5.6) of C 1 , we find and the bound on (A 1 ) −1 S A D k in (5.16) follows since 1/q = 1/2 − 1/s. To prove the bound on (A 1 ) −1 S A 2 in (5.17), we use the bounds Remark 5.11 (Link to the results of [26]) A result analogous to the Euclidean-norm bounds in Lemma 5.3 was proved in [26,Theorem 1.4] for the case that A 1 = A 2 = I, n 2 = 1, and n 1 = 1 + iε/k 2 , with the 'absorption parameter' or 'shift' ε satisfying 0 < ε k 2 . This is equivalent to approximating the discrete inverse of ∆ + k 2 by the discrete inverse of ∆ + k 2 + iε. The motivation for proving this result was that the so-called 'shifted Laplacian preconditioning' of the Helmholtz equation (introduced in [19]) is based on preconditioning, with these choices of parameters, A 2 with an approximation of A 1 . Similar to the proofs of Theorems 2.9 and 2.10 in §5.2.1, bounds on I − (A 1 ) −1 A 2 2 and I − A 2 (A 1 ) −1 2 then give upper bounds on how large the 'shift' ε can be for GMRES for (A 1 ) −1 A 2 to converge in a k-independent number of iterations in the case when the action of (A 1 ) −1 is computed exactly.
The main differences between [26] and this work are that (i) [26] focuses only on the case when T = ik, (ii) the theory in [26] focuses on the particular case that D − is star-shaped with respect to a ball, finding a k-and ε-explicit expression for C (1) bound in this case using Morawetz identities, (iii) the proof of [26, Theorem 1.4] requires bounds on (A 1 ) −1 M n , M n (A 1 ) −1 (analogous to the bounds in Lemma 5.7), (A 1 ) −1 N, and N(A 1 ) −1 , but not on (A 1 ) −1 S A and S A (A 1 ) −1 , (iv) [26,Theorem 1.4] holds under the assumption that the finite-element approximation is quasi-optimal, whereas here we prove results under Condition 4.2, which holds under less restrictive conditions on h and p than those for quasioptimality (see §4.3), and (v) [26] only proves bounds on I − (A 1 ) −1 A 2 and I − A 2 (A 1 ) −1 in the · 2 norm whereas here we also prove bounds in weighted norms.

Proofs of Theorems 2.2 and 2.3 from Theorems 2.9 and 2.10
Theorems 2.2 and 2.3 follow from repeating the proofs of Theorems 2.9 and 2.10 but now, instead of using Hölder's inequality (5.23), using the simpler inequality The result is that the inverse inequality (5.1) is no longer needed in (5.24) to obtain f L 2 (D) from f L s (D) . Then Lemmas 5.7, 5.8, and 5.3, and Theorems 2.9 and 2.10 all hold without the factors h −d/q on the right-hand sides of the inequalities, and with the constants C j , j = 1, 2, replaced by C j , j = 1, 2. This is equivalent to these results holding with q set formally as ∞, since C inv,s = 1 when s = 2 (which occurs when q = ∞), as noted immediately after ( . Since Condition 4.1 holds, we can then apply the result of Lemma 5.10, i.e. the bound (5.25), to the solution of the variational problem (6.1) to find that Proof of Lemma 2.7 We actually prove the stronger result that, given any function c(k) such that c(k) > 0 for all k ≥ k 0 , there exist f, n 1 , n 2 , with n 1 = n 2 and n 1 − n 2 L ∞ (D;R) ∼ c(k), such that the corresponding solutions u 1 and u 2 of the exterior Dirichlet problem with A 1 = A 2 = I exist, are unique, and satisfy (2.14). The heart of the proof of (2.14) is the equation (∆ + k 2 ) e ikr χ(r) = e ikr ik d − 1 r χ(r) + 2ik dχ dr (r) + ∆χ(r) =: − f (r), (6.3) where χ(r) is chosen to have supp χ ⊂ D. This example proves the sharpness of the nontrapping resolvent estimate (4.2) as k → ∞, since both the L 2 (D) norm of f and the H 1 k (D) norm of e ikr χ(r) are proportional to k, and hence each to other (see, e.g., [11,Lemma 3.10], [63,Lemma 4.12]).
The overall idea of the proof of (2.14) is to set things up so that (u 1 − u 2 )(x) = e ikr χ(r), the rationale being that (6.3) proves the sharpness of (4.2), and (4.2) and its corollary (5.25) (applied to u 1 − u 2 ) are the main ingredients in the proof of Theorem 2.6.
Then, with f defined by (6.3), we formally define u 2 and f by , and f (x) := − ∆ + k 2 n 2 (x) u 2 (x), (6.5) and we define u 1 to be the solution of the TEDP (in the sense of Definition 3.2) with n 1 = 1, g = 0, and f defined as in (6.5). Since χ(r) has compact support in D, we need to tie both the support of χ and how fast χ vanishes near the boundary of its support to the definition of χ for both u 2 and f in (6.5) to be well defined. We ignore this issue for the moment and show that these definitions achieve our goal of having u 1 (x) − u 2 (x) = e ikx1 χ(r). (6.6) First observe that since supp f is a compact subset of D, so is supp u 2 . Therefore u 2 is the solution of the TEDP with g = 0, f defined in (6.5), and coefficient n 2 := 1 + c(k) χ. Combining (6.4) with the definitions of n 2 and u 2 , we see that (∆ + k 2 )(u 1 − u 2 ) = − f . (6.7) Since the solution of the TEDP is unique (see Remark 3.4), (6.3) and (6.7) imply that (6.6) holds. We therefore have that u 1 − u 2 L 2 (D) ∼ 1 and u 1 − u 2 H 1 k (D) ∼ k.
Furthermore, the definitions (6.5) and (6.3) of u 2 and f , respectively, imply that and and, since n 1 − n 2 L ∞ (D;R) = c(k), (2.14) holds. Therefore, to complete the proof, we only need to show that there exists a choice of χ and χ for which u 2 and f defined by (6.5) are in H 1 (D) and L 2 (D) respectively (in fact, we prove that they are in W 1,∞ (D) and L ∞ (D) respectively). Since χ and χ ∈ C ∞ (D), and χ is positive everywhere in the interior of its support, the only issue is what happens at the boundary of supp χ, where u 2 could be singular.
It is now much harder than in (6.4) to set things up so that u 1 (x) − u 2 (x) = e ikr χ(r) (so that one can then use (6.3)).
Remark 6.2 (The relationship between Q1 and Q1 ′ ) In §2.2, we stated that Q1 ′ is an analogue of Q1 on the PDE level; we now make this statement more precise.
Since C 3 defined by (6.2) is independent of F , the relative-error bound (2.13) of Theorem 2.6 can be written as On the other hand, the main ingredient to proving Theorem 2.2 is the bound which implies op) , C 2 n 1 − n 2 L ∞ (D;R) . (6.10) The left-hand side of (6.10) is the discrete analogue of the left-hand side of (6.9). The right-hand sides of these bounds are (up to constants independent of k) identical, hence why we say that the condition (2.11) (i.e. k A 1 − A 2 L ∞ (D;op) and k n 1 − n 2 L ∞ (D;R) both sufficiently small) is the answer to both Q1 and Q1 ′ .