The Helmholtz equation with uncertainties in the wavenumber

We investigate the Helmholtz equation with suitable boundary conditions and uncertainties in the wavenumber. Thus the wavenumber is modeled as a random variable or a random field. We discretize the Helmholtz equation using finite differences in space, which leads to a linear system of algebraic equations including random variables. A stochastic Galerkin method yields a deterministic linear system of algebraic equations. This linear system is high-dimensional, sparse and complex symmetric but, in general, not hermitian. We therefore solve this system iteratively with GMRES and propose two preconditioners: a complex shifted Laplace preconditioner and a mean value preconditioner. Both preconditioners reduce the number of iteration steps as well as the computation time in our numerical experiments.


Introduction
The Helmholtz equation is a linear partial differential equation (PDE), whose solutions are time-harmonic states of the wave equation, see [14,19].Important applications of this model are given in acoustics and electromagnetics [2].The Helmholtz equation includes a wavenumber, which is either a constant parameter or a space-dependent function.Furthermore, boundary conditions are imposed on the spatial domain.
We consider uncertainties in the wavenumber.Thus the wavenumber is replaced by a random variable or a spatial random field to quantify the uncertainties.The solution of the Helmholtz equation changes into a random field, which can be expanded into the (generalized) polynomial chaos, see [30].We employ the stochastic Galerkin method to compute approximations of the unknown coefficient functions.Stochastic Galerkin methods were used for linear PDEs of different types including random variables, for example, see [12,32] on elliptic type, [13,22] on hyperbolic type, and [21,31] on parabolic type.Wang et al. [29] applied a multi-element stochastic Galerkin method to solve the Helmholtz equation including random variables.We investigate the ordinary stochastic Galerkin method, which is efficient if the wavenumbers are not close to resonance.
The stochastic Galerkin method transforms the random-dependent Helmholtz equation into a deterministic system of linear PDEs.Likewise, the original boundary conditions yield boundary conditions for this system.We examine the system of PDEs in one and two space dimensions.A finite difference method, see [15], produces a high-dimensional linear system of algebraic equations.When considering absorbing boundary conditions, the coefficient matrices are complex-valued and non-hermitian.
We focus on the numerical solution of the linear systems of algebraic equations.The dimension of these linear systems rapidly grows for increasing numbers of random variables.Hence we use iterative methods like GMRES [24] in the numerical solution.The efficiency of an iterative method strongly depends on an appropriate preconditioning of the linear systems.We propose two preconditioners in the general case where the wavenumber can depend on space and on multiple random variables: a complex shifted Laplace preconditioner, see [6,8], and a mean value preconditioner, see [10,29].Statements on the location of spectra and estimates of matrix norms are shown.Furthermore, results of numerical computations are presented for both settings.
The article is organized as follows.The stochastic Helmholtz equation is introduced in Section 2 and discretized in Section 3. We discuss the complex shifted Laplace preconditioner in Section 4 and the mean value preconditioner in Section 5.Sections 6 and 7 contain numerical experiments in one and two spatial dimensions, respectively, which show the effectiveness of the preconditioners.An appendix includes the detailed formulas of the discretizations in space.

Problem Definition
We illustrate the stochastic problem associated to the Helmholtz equation.

Helmholtz equation
The Helmholtz equation is a PDE of the form Often homogeneous Dirichlet boundary conditions, i.e., u = 0 on ∂Q, ( are applied for simplicity.Alternatively, absorbing boundary conditions read as where ∂ n denotes the derivative with respect to the outward normal of Q and i = √ −1 is the imaginary unit.

Stochastic modeling
We consider uncertainties in the wavenumber.A simple model to include a variation of the wavenumber is to replace the constant k by a random variable on a probability space (Ω, A, P ).We write k = k(ξ), where ξ : Ω → R is some random variable with a traditional probability distribution.More generally, the wavenumber can be a spacedependent function on Q including a multidimensional random variable ξ : Ω → Ξ with Ξ ⊆ R s .We assume ξ = (ξ 1 , . . ., ξ s ) with independent random variables ξ for = 1, . . ., s.Now the wavenumber becomes a random field with given functions k : Q → R for = 0, 1, . . ., s, as in [29].A truncation of a Karhunen-Loève expansion, see [11, p. 17], also yields a random input of the form (2.4).Consequently, the solution of the deterministic Helmholtz equation (2.1) changes into a random field u : Q × Ξ → K.We write u(x, ξ) to indicate the dependence of the solution on space as well as the random variables.We assume that each random variable ξ has a probability density function ρ .Since the random variables are independent, the product ρ = ρ 1 • • • ρ s is the joint probability density function.Without loss of generality, let ρ(ξ) > 0 for almost all ξ ∈ Ξ.The expected value of a measurable function f : Ξ → K depending on the random variables is (2.5) In the following, L 2 (Ξ, ρ) denotes the Hilbert space of square-integrable functions.The associated norm is f L 2 (Ξ,ρ) = f, f .Later we will focus on uniformly distributed random variables ξ : Ω → [−1, 1].In this case, the joint probability density function is constant, i.e., Ξ = [−1, 1] s and ρ ≡ 2 −s .

Polynomial chaos expansions
We assume that there is an orthonormal polynomial basis (φ i ) i∈N 0 in L 2 (Ξ, ρ).Thus it holds that φ i , φ j = δ i,j = 1 for i = j 0 for i = j with the inner product (2.5).In the case of uniform probability distributions, the multivariate functions φ i are products of the (univariate) Legendre polynomials.We assume that φ 0 ≡ 1.The number m + 1 of multivariate polynomials in s variables up to a total degree r is see [30, p. 65].This number grows fast for increasing r or s.
with (a priori unknown) coefficient functions The series (2.7) converges in L 2 (Ξ, ρ) pointwise for x ∈ Q.If the wavenumber k is an analytic function of the random variables, then the rate of convergence is exponentially fast for traditional probability distributions.

Discretization of the stochastic Helmholtz equation
We consider the stochastic Helmholtz equation with given source term f : Q → R and random wavenumber k : Q × Ξ → R + , together with either homogeneous Dirichlet boundary conditions or with absorbing boundary conditions All derivatives are taken with respect to x.We discretize this boundary value problem in two steps, with a finite difference method (FDM) in space and the stochastic Galerkin method in the random-dependent part.The steps can be done in any order.We first give an overview of the procedure when beginning with the FDM in Section 3.1.In Section 3.2, we discuss the discretization when beginning with the stochastic Galerkin method.

FDM and stochastic Galerkin method
A spatial discretization of the boundary value problem with a FDM leads to a (stochastic) linear algebraic system with S(ξ) ∈ K n,n for ξ ∈ Ξ (see Section A for details) and constant vector F 0 ∈ R n .In a second step, we consider a PC approximation of U (ξ) of the form and φ i are polynomials as in Section 2.3.The coefficient vectors V i are determined by the orthogonality of the residual to the subspace span{φ 0 , φ 1 , . . ., φ m } with respect to the inner product •, • in (2.5), i.e., by R m (ξ), φ i (ξ) = 0 for i = 0, 1, . . ., m.Here the inner product is taken componentwise.The orthogonality condition is equivalent to due to φ 0 ≡ 1.This leads to a (deterministic) linear algebraic system where the stochastic Galerkin projection A ∈ K (m+1)n,(m+1)n is a block matrix with m + 1 blocks of size n × n, and F i = 0 ∈ R n for i = 1, . . ., m.
Remark 3.1.The Galerkin approximation (3.5) can be interpreted as a spatial discretization of a Galerkin approximation As it turns out, the matrix S(ξ) is a (complex) linear combination of symmetric positive (semi-)definite matrices.The following lemma shows that this structure is preserved in the stochastic Galerkin method; see [20, Lem.1] and its proof.These properties of the matrix S and thus A will be essential for our analysis of shifted Laplace preconditioners in Section 4.
and the stochastic Galerkin projection We then obtain for i, j = 0, 1, . . ., m where the inner product is taken component-wise.Additionally, A ij = A ji .Moreover, if A(ξ) is symmetric, then A is symmetric, and if A(ξ) is symmetric positive (semi-)definite for almost all ξ ∈ Ξ then A is symmetric positive (semi-)definite.
Corollary 3.3.In the notation of Lemma 3.2, if A(ξ) = A 0 is independent of ξ, then A ij = δ ij A 0 and A = I m+1 ⊗ A 0 , with the identity matrix I m+1 ∈ K m+1,m+1 and the Kronecker product.
Finally, we obtain the following result on the structure of the matrix A in (3.8).Proof.The FD discretizations leading to S(ξ)U (ξ) = F 0 in (3.4) are given in Section A.
The statement of the theorem follows in each case by applying the stochastic Galerkin approximation as described above and using Lemma 3.2 as well as Corollary 3.3 separately for each term composing S(ξ).
The matrix L results essentially from the discretization of the Laplacian, B from the (absorbing) boundary conditions, and K is the discretization of the term including the wavenumber; see Section A for details.

Stochastic Galerkin method and FDM
Alternatively, we can begin with the stochastic Galerkin method.This leads to a system of deterministic PDEs, which are subsequently discretized by a FDM.The PC expansion (2.7) suggests a stochastic Galerkin approximation of u(x, ξ) of the form (3.14) The coefficient functions v i,m in the stochastic Galerkin method are in general distinct from the coefficients v i in (2.8).Nevertheless, we will usually write v i instead of v i,m in the sequel for notational convenience.The coefficients in the Galerkin approach are determined by the orthogonality of the residual to the subspace span{φ 0 , φ 1 , . . ., φ m }, i.e., by R m (x, ξ), φ j (ξ) = 0 for j = 0, 1, . . ., m and each x ∈ Q.The latter is equivalent to for j = 0, 1, . . ., m in Q.Thus we obtain a system of PDEs for the unknown coefficient (3.16) Since by assumption k(x, ξ) > 0 for all x and ξ, the matrix C(x) is symmetric positive definite (as Gramian of an inner product with weight function k(x, ξ) 2 ρ(ξ)).Setting we write the system of PDEs (3.15) as which is a larger deterministic system of linear PDEs.Still we require boundary conditions for the system (3.18).The homogeneous Dirichlet boundary condition (3.2) implies v j (x) = 0 for x ∈ ∂Q and j = 0, 1, . . ., m, hence v(x) = 0 on ∂Q.
Inserting the Galerkin approximation (3.14) into the absorbing boundary conditions (3.3) yields the residual By the orthogonality R m (x, ξ), φ j (ξ) = 0 in the Galerkin approach, we obtain The matrix is symmetric and positive definite (since k(x, ξ) > 0 by assumption).The boundary condition (3.21) can be written with B(x) as The boundary value problem (3.18) with (3.19) or (3.23) is discretized in Section A.4 (in dimension d = 1).The resulting linear algebraic system is the same as the one obtained in Section 3.1.

Complex shifted Laplace preconditioner
Following the investigation in [9], we consider the Helmholtz equation (3.1) with a complex shift in the wavenumber with β ∈ R, together with either homogeneous Dirichlet boundary conditions (3.2) or absorbing boundary conditions (3.3).We discretize this boundary value problem as described in Section 3.1.For β = 0, we have the matrix (3.13) in Theorem 3.4, and for β ∈ R we obtain since only the constant term is multiplied by 1 + iβ.Motivated by [27, p. 1945], we call M a complex shifted Laplace preconditioner (CSL preconditioner).
For the deterministic Helmholtz equation, preconditioning with the CSL preconditioner is a widely studied and successful technique for solving the discretized Helmholtz equation; see, e.g., [6,1,23,3,7] and [9], as well as references therein.See also [5] for a survey and [17] for recent developments.In the deterministic case, the spectrum of the preconditioned matrix AM −1 lies in the disk (4.3), and the improved localization of the spectrum typically leads to a faster convergence of Krylov solvers.The CSL preconditioner M can be inverted efficiently, for example, by multigrid techniques.
Here, we focus on locating the spectrum of the preconditioned matrix in the stochastic case, in analogy to [27,8,9] for the deterministic Helmholtz equation.1.In the case of absorbing boundary conditions (3.3), the spectrum of the preconditioned matrix AM −1 is contained in the closed disk 2. In the case of homogeneous Dirichlet boundary conditions (3.2), the spectrum of the preconditioned matrix AM −1 lies on the circle Proof.We begin with the case of absorbing boundary conditions.The proof closely follows [27,Sect. 3] with minor modifications.We have with z 1 = 1 and z 2 = 1 + iβ and where L, B, K are symmetric, K is positive definite and L, B are positive semidefinite; see Theorem 3.4.Then A and M are of the form in [27,Sect. 3], except for the opposite sign of B. The opposite sign affects the positive semidefiniteness, but not the overall strategy of the proof.Nevertheless, we give a full proof here.
Step 1: Observe first that AM −1 and M −1 A have the same spectrum, and that M −1 Ax = σx is equivalent to the generalized eigenproblem Ax = σM x.
Step 2: x is an eigenvector of Ax = σM x if and only if (L − iB)x = λKx, which can be seen as follows: , implies that M is singular and thus not eligible as preconditioner.) Step 3: Location of λ in the generalized eigenvalue problem (L − iB)x = λKx.Since K is symmetric positive definite, it has a Cholesky factorization K = U U and the generalized eigenvalue problem is equivalent to where y = U x. Multiplication of (4.7) by y and division by y y yields This shows Re(λ) ≥ 0 and Im(λ) ≤ 0 since L and B are symmetric positive semidefinite.
Step 4: Estimate of the eigenvalues is a Möbius transformation.By step 2, σ = µ(λ) where λ is an eigenvalue of the generalized eigenvalue problem (L−iB)x = λKx which satisfies Im(λ) ≤ 0. To determine µ(R), we compute and Hence µ maps the real line onto the circle C in (4.4) (for any β = 0).For β > 0, the lower half-plane is mapped by µ onto the interior of C (for β < 0 onto the exterior); see Figure 1.This completes the proof in case of absorbing boundary conditions.The proof in the case of Dirichlet boundary conditions is very similar.The only difference is in the location of the eigenvalues λ in step 3. Since L is symmetric positive definite and B = 0, (4.7) implies λ > 0, hence σ = µ(λ) lies on the circle (4.4).Remark 4.2.In the proof of Theorem 4.1, we additionally have Re(λ) ≥ 0. Hence σ is located in the image of the (closed) fourth quadrant under µ in (4.9).To determine this image, note that µ maps the imaginary axis onto the circle which intersects µ(R) = C orthogonally in µ(0) and µ(∞) = 1.Considering the orientations shows that µ maps the right half-plane onto the exterior of C β ; see Figure 1.Thus the spectrum satisfies In case of Dirichlet boundary conditions, the eigenvalues of AM −1 lie on the arc of the circle C from µ(0) to µ(∞) = 1 that contains the origin.This observation further tightens the inclusion set of σ(AM −1 ), also in the case of a deterministic wavenumber.This tighter inclusion set is already visible in [8, Fig. 1, Fig. 2] and [9, Fig. 2.1] but we are not aware of a proof in the literature.

Mean value preconditioner
We consider the discretization from Section 3.1.Let S(ξ) ∈ K n,n be the coefficient matrix of a linear system resulting from a spatial discretization of the Helmholtz equation (2.1) including boundary conditions and wavenumber k(x, ξ).We assume that S(ξ) is nonsingular for almost all realizations ξ ∈ Ξ.Let ξ ∈ Ξ be the expected value of the multidimensional random variable ξ.It holds that The stochastic Galerkin method applied to S(ξ) yields a matrix A ∈ K (m+1)n,(m+1)n as shown in Section 3.1.Furthermore, we define the constant matrix Ā = I m+1 ⊗ S( ξ). (5.1) This matrix allows for the construction We employ the Frobenius matrix norm • F in the following.
Theorem 5.1.Using the Frobenius norm, it holds that with the constants provided that the L 2 -norm of the matrix norm is finite.
Proof.The definition (5.2) directly yields We obtain Ā−1 ∆A F ≤ Ā−1 F ∆A F .The properties of the Kronecker product and (5.1) imply Ā−1 2 F = (m + 1) S( ξ) −1 2 F .We estimate ∆A F using the Cauchy-Schwarz inequality with respect to the inner product (2.5) In the last step, we used that the square of an L 2 -norm is an integral and thus summation (with respect to µ, ν) and integration can be interchanged.Applying the square root to the above estimate yields the statement (5.3).Theorem 5.1 together with Remark 5.2 demonstrate that the matrix Ā is a good preconditioner for solving linear systems with coefficient matrix A. In this context, Ā is called the mean value preconditioner, as in [29] for the multi-element method.When Ā is used as a preconditioner (left-hand or right-hand), linear systems with coefficient matrix Ā have to be solved.The matrix Ā from (5.1) is block-diagonal with m + 1 identical blocks in this application.Thus just a single LU -decomposition of the matrix S( ξ) is required.Many linear systems with different right-hand sides are solved using this LU -decomposition in an iterative method like GMRES, for example.
Theorem 5.5.Let S(ξ) = S 0 + θT (ξ) with a non-singular constant matrix S 0 , a matrix T = t µ,ν µ,ν depending on a random variable ξ with components t µ,ν ∈ L 2 (Ξ, ρ) and a real parameter θ > 0. Using A 0 = I m+1 ⊗ S 0 , the Frobenius norm exhibits the asymptotic behavior A −1 0 A − I (m+1)n F = O(θ). (5.4) Proof.Since the entries of T (ξ) are assumed to be square-integrable, also the expected values are finite.Let T be the constant matrix containing the expected values of T (ξ).We apply the decomposition The matrix S 0 + θ T is non-singular for sufficiently small θ.Moreover, we obtain the relation which confirms (5.4).
An important case of Theorem 5.5 is T = 0, i.e., these expected values are zero.Then A 0 = Ā is the mean value preconditioner.for all sufficiently small ∆S.
Likewise, the Frobenius norm using A 0 instead of Ā is smaller than one if the parameter θ is sufficiently small in the context of Theorem 5.5.
A stationary iterative scheme for solving a linear system Ax = b reads as with a non-singular matrix B which should approximate A, see [26, p. 621].In each iteration step, we have to solve a linear system with coefficient matrix B. The property (5.5) is sufficient for the global convergence of the iteration (5.6) using B = Ā.The computational cost of an iteration step is much less than the steps in GMRES using Ā as preconditioner, because the construction of Krylov subspaces is avoided.In practice, we do not know if ∆S is sufficiently small such that the bound (5.5) is guaranteed.Nevertheless, it is worth to try this stationary iteration, as we will observe in Section 7.
In our numerical experiments in one and two spatial dimensions, we compute the mesh-size h = 1 q+1 in the FD discretization by lev = max(ceil(log2((15*maxk)/(2*pi))), 1); q = 2^lev -1; where maxk is the maximal value of the wavenumber.Then the relation 2π kh ≈ constant, advocated in [16, Sect.4.4.1], is satisfied.Indeed, the estimate x ≤ x ≤ x + 1 for x ∈ R implies 15k 2π ≤ q + 1 ≤ 2 15k 2π for large k.In particular, q grows linearly with k and thus the size of the matrices S(ξ) and A (see Section A) grows with k; see, e.g., Figure 3.Our choice for q can be adapted for a future use of a multigrid method (as in [9]).
Discretizing the model problem yields a linear algebraic system as given in Theorem A.2.This one-dimensional problem can be solved by a direct method, since the computational work is not too large.Nevertheless we also consider its solution with the GMRES method [24] and investigate the application of CSL and mean value preconditioners introduced in Sections 4 and 5, respectively.By Theorem A.2, the matrix A has the form If needed, we write A θ to indicate the dependence of A on θ, and in particular A 0 for θ = 0, which corresponds to the mean value preconditioner.Since the wavenumber in (6.1) is constant in space, the matrices [B ij ] and [C ij ] simplify to see Lemma A.4, and, by Lemma A.3, In other words, the matrices k(ξ)φ j (ξ), φ i (ξ) ij and k(ξ) 2 φ j (ξ), φ i (ξ) ij are tridiagonal and pentadiagonal, respectively, due to the orthogonality properties of the polynomials.

diagonal, and
with S(ξ) from Theorem A.2.This shows that A 0 is block-diagonal with m + 1 identical diagonal blocks.The latter are the FD-discretization of the deterministic Helmholtz equation with wavenumber k (associated to ξ = 0).
If not specified otherwise, we use m = 3 in the stochastic Galerkin method and θ = 0.1 in (6.1).Finally, we also consider the shifted Helmholtz equation (4.1) with shift β = 1 2 and denote the CSL preconditioner by M = M ( 1 2 ), see (4.2).As for A, we write M θ if we wish to emphasize the dependence on θ.
The numerical experiments have been performed in the software package MATLAB R2020b on an i7-7500U @ 2.70GHz CPU with 16 GB RAM.

Spectra
By Theorem 4.1, the eigenvalues of the CSL preconditioned matrix AM −1 lie in the closed disk (4.3).This is illustrated in the left panel of Figure 2, which displays the spectra of AM −1 (with θ = 0.1) and A 0 M −1 0 (i.e., with θ = 0).Each eigenvalue of  as functions of k (with θ = 0.1).Clearly, the condition numbers of M and AM −1 are much smaller than the condition number of A, which is beneficial when solving the preconditioned linear system AM −1 y = b, M x = y with the CSL preconditioner.In this example, κ 2 (M ) ≤ 205 for all k, which is very moderate, and κ 2 (AM −1 ) grows linearly in k from 2.6485 when k = 10 to only 36.5190 when k = 200.In contrast, κ 2 (A) is roughly 50 to 160 times larger than κ 2 (AM −1 ).The observed spikes of κ 2 (A) occur when more discretization points are used which leads to a larger size of A, compare the curve of size(A).The condition number of the mean value preconditioned matrix AA −1 0 is also moderate, growing from 2 to 141, which is beneficial for solving the preconditioned linear system, while κ 2 (A 0 ) is of the order of κ 2 (A).

GMRES
We solve the unpreconditioned system (6.2) and the right and left preconditioned systems with full GMRES (no restarts) and tolerance tol=1e-12, using MATLAB's built-in gmres command.The residual in the ith step is r (i) = b − Ax (i) for unpreconditioned and right preconditioned GMRES, and M −1 r (i) for left preconditioned GMRES.In particular, the stopping criterion for left and right preconditioning is in general different.We will consider the following three preconditioners: 1. the CSL preconditioner M , 2. the mean value preconditioner A 0 , 3. the mean value CSL preconditioner M 0 .
In preconditioned GMRES, we need to solve linear systems with the preconditioner, for which we use an LU -decomposition.In one spatial dimension, this is not competitive with the direct solution (see the end of Section 6.3), but in two spatial dimension the block structure of the preconditioners A 0 and M 0 leads to a competitive method.In MATLAB, the LU -decomposition of the sparse matrix M calls the associated routine from UMFPACK; see [4].The decomposition has the form with a lower triangular matrix L, upper triangular matrix U , and two permutation matrices P, Q.In our implementation, we use [L, U, p, q] = lu(M, 'vector'); qt = []; qt(q) = 1:numel(q); where, instead of the matrices P, Q, only vectors p, q representing the permutations are stored, and where the vector qt describes the inverse mapping of the permutation defined by q.Then, we implement M −1 x by x = U\(L\x(p,:)); x = x(qt,:); Table 1: Time in seconds (s) for solving (6.2) and (6.8) with GMRES.By Remark 6.1, A 0 = I m+1 ⊗ S(0) is block-diagonal with equal diagonal blocks so that, for fixed k, only a single LU -decomposition of S(0) ∈ K n,n is necessary to compute A −1 0 x for any vector x ∈ K (m+1)n .In our implementation, we partition and reshape x so that only one linear system with S(0) is solved: x = reshape(x, [n, m+1]); x = U\(L\x(p,:)); x = x(qt,:); x = reshape(x, [], 1); The preconditioner M 0 is implemented in the same way.
In a first experiment, we fix k = 50, θ = 0.1 and m = 3. Solving the unpreconditioned system (6.2) with GMRES suffers from a long delay of convergence; see Figure 4.In contrast, all three preconditioners M , M 0 , and A 0 lead to a significant decrease in the number of iteration steps from about 250 to 50 for M and M 0 (factor 5), and to about 25 for A 0 (factor 10); see Figure 4 (left panel).The computation times with the preconditioners M and M 0 reduce to about 6% of the computation time of unpreconditioned GMRES, while for A 0 it reduces to about 3%; see Table 1.The computation times for the preconditioned systems include the computation of the LU -decomposition (of M or of a diagonal block for A 0 or M 0 ).The differences between computed solutions are very small: x − x ∞ ≤ 1.7 • 10 −14 (and typically of order 10 −15 ), where x=A\b is the direct solution and x is a solution computed with GMRES (unpreconditioned or with one of the preconditioners).Left and right preconditioning lead to very similar relative residual norms and timings for each preconditioner.A heuristic explanation why A 0 performs better than M and M 0 , is that A is closer to A 0 than to M or M 0 .Indeed, we have A−A 0 ∞ < A−M ∞ < A−M 0 ∞ in this example.Repeating this experiment with θ = 0.2 leads to very similar results, see Figure 4 and Table 1, so we focus on θ = 0.1.
In a second experiment, we let k vary while θ = 0.1 and m = 3 are fixed.Figure 5 displays the number of GMRES iteration steps (top) and the computation time (bottom) as functions of k.For small k ∈ [10,50], the difference between unpreconditioned and preconditioned GMRES is not so pronounced, since the linear systems are rather small.For 60 ≤ k ≤ 200, the three preconditioners significantly reduce the number of iteration steps and the computation time compared to unpreconditioned GMRES.The number of iteration steps is reduced to 8-15% of the number of iteration steps in unpreconditioned GMRES when using M , to 9-16% when using M 0 and to only 3-6% when using A 0 as preconditioner.GMRES preconditioned with M or M 0 needs only 1-4% of the computa-Figure 6: Plots of the coefficients v 0 , v 1 , v 2 , v 3 for k = 50 and θ = 0.1.tion of unpreconditioned GMRES, and the computation time of GMRES preconditioned with A 0 is reduced to 0.5-1.1% of the computation time of unpreconditoned GMRES.The mean value preconditioner A 0 leads to the smallest number of GMRES iteration steps and computation time, which is likely due to the fact that A is closer to A 0 than to M or M 0 .Note, however, that the condition number of A 0 (and A) is much larger than that of M and M 0 .For k = 150, we have (rounded to the nearest integer) κ 2 (A) = 2428, κ 2 (A 0 ) = 2220, κ 2 (M ) = 109, κ 2 (M 0 ) = 91; see also Figure 3. Thus, if accuracy is an issue, it is preferable to work with the CSL preconditioners M or M 0 .
Finally, we note that the direct solution A\b with a sparse matrix in MATLAB calls an efficient algorithm from UMFPACK; see [4].In the above test example, solving the linear system (6.2) by GMRES (with or without preconditioner) is not competitive with this direct solution, as it is much faster; see the bottom right panel in Figure 5.

Solutions
Figure 6 displays the real and imaginary parts of the computed coefficients v 0 , v 1 , v 2 , v 3 in the Galerkin approximation for k = 50 and θ = 0.1 in (6.1).We recognize an effect of the point source at x = 1 2 in the real part of v 0 .We compute the solution for total polynomial degree m = 100.Figure 7 shows v i ∞  as a function of the polynomial degree i.In the left panel, θ = 0.1 is fixed and k varies, while in the right panel θ varies and k = 100 is fixed.We observe an exponential decay of the coefficients in all cases, which is related to the exponential convergence of the PC expansion (2.7).Larger wavenumbers and larger values of θ lead to a slower decay of the maximum-norm of the coefficients.The effect of larger θ on the convergence/decay is more pronounced, compare, for example, the curve for (k, θ) = (150, 0.1) in the left panel with the curve for (100, 0.5) in the right panel.
Next, we vary m (the maximal degree of the polynomials in the stochastic Galerkin method) and denote by x m the solution of (6.2), which consists of a discretization of the coefficients v 0,m , . . ., v m,m in a Galerkin approximation (3.14) of the solution u of the Helmholtz equation; see also Remark 3.1.The convergence of the stochastic Galerkin method is illustrated by the exponential decay of the norms x m − x m+1 2 in Figure 8
Let k 1 = 30, k 2 = 15, k 3 = 20, and θ = 0.1.Table 2 shows the size of A and the time (in seconds) for constructing the matrix A for polynomial degrees up to r = 0, 1, . . ., 8 in the stochastic Galerkin method.The computation time when solving Ax = b directly with x=A\b in MATLAB grows exponentially as a function of r, see Figure 9 (left panel).Thus we solve the linear algebraic system Ax = b with GMRES using the mean value preconditioner A 0 = I m+1 ⊗ S 0 from (5.1) as right preconditioner, that is, we solve Here S 0 denotes the FD discretization of the Helmholtz equation with absorbing boundary conditions and deterministic wavenumber (7.2); see Theorem A.6.The solution of linear systems with the preconditioner A 0 is implemented as described in Section 6.3.We solve (7.4) with full GMRES (no restarts), tol=1e-8 and maxit=200 for polynomial degrees up to r = 1, . . ., 8 in the stochastic Galerkin method.In contrast to the experiments in 1D in Section 6, preconditioned GMRES is significantly faster than the direct solution with MATLAB's 'backslash' command x=A\b; see Figure 9 (left panel).
The computation times for preconditioned GMRES include the computation of the LUdecomposition of a diagonal block of A 0 .Furthermore, the relative residual norms in GMRES for polynomial degree r = 8 are shown in Figure 9 (right panel).The mean value CSL preconditioner M 0 has a similar block-diagonal structure to A 0 and performs similarly well; see Figure 9.
Alternatively to GMRES or a direct solution of the linear system, we also investigate the stationary iteration (5.6) with B = A 0 , i.e., We take the starting vector x (0) = A −1 0 b.Linear systems with the matrix A 0 are solved as described above.For θ = 0.1, this iteration converges.Figure 10 displays the relative error norms in the maximum-norm for polynomial degrees r = 2, 4, 6, where we take the direct solution A\b as the 'exact' solution.The slower convergence for larger degree r in the stochastic Galerkin method is expected, since the matrix size also grows causing Figure 10: Relative error norms in the stationary iteration (7.5) for different polynomial degrees r in the stochastic Galerkin method.higher condition numbers.For θ = 0.2, the stationary iteration diverges.This behavior is in agreement to Theorem 5.5 and Corollary 5.6.
In Figure 11, the top row displays the expected value of the real and imaginary part of the computed stochastic Galerkin approximation u m (with polynomial degree r = 5).The variance is displayed in the bottom row of the figure.
Denote by x r the solution of Ax = b when using polynomials of degree up to r in the stochastic Galerkin method, where the number of basis polynomials is given in (7.3).The left panel of Figure 12 displays the differences x r−1 − x r 2 as a function of r (the vector x r−1 is padded with zeros at the end to match the size of x r ).Their exponential decay suggests convergence of the stochastic Galerkin method.
Next, we fix the degree r = 8 in the stochastic Galerkin method.Recall from (3.5) and (3.8) that the solution of Ax = b contains the coefficient vectors V 0 , . . ., V m of the polynomials φ 0 , . . ., φ m in the stochastic Galerkin method.We also examine the largest maximum norm of the coefficients associated to polynomials of total degree (exactly) j, i.e., the values The right panel of Figure 12 shows the magnitudes (7.6) for j = 0, 1, . . ., 8. The observed exponential decay stems from the exponential convergence of (2.7), since the wavenumber in (7.1) is analytic in ξ 1 , ξ 2 , ξ 3 .We repeat this experiment with θ = 0.2 instead of θ = 0.1.Overall, the behavior is similar as for θ = 0.1, but convergence is slower: the relative residual norms reach the prescribed tolerance in 60 instead of 20 iteration steps, and also x r−1 − x r 2 as well as the magnitudes (7.6) converge more slowly.

Conclusions
We investigated the Helmholtz equation including a random wavenumber.The combination of a stochastic Galerkin method and a finite difference method yielded a high-dimensional linear system of algebraic equations.We examined the iterative solution of these linear systems using three types of preconditioners: a complex shifted Laplace preconditioner, a mean value preconditioner, and a combined variant.Theoretical properties of the preconditioned linear systems were shown.Numerical computations demonstrate that the straightforward mean value preconditioner provides the most efficient iterative solution within these types.

1 ) 2 ∂x 2 j
with an (open) spatial domain Q ⊆ R d and given source term f : Q → R. The wavenumber k is either a positive constant or a function k : Q → R + .The unknown solution is u : Q → K with either K = R or K = C.Here ∆ = d j=1 ∂ denotes the Laplace operator with respect to x = [x 1 , . . ., x d ] ∈ R d .

Theorem 3 . 4 .
Let the spatial dimension be d ∈ {1, 2}.A finite difference and stochastic Galerkin approximation of the Helmholtz equation (3.1) with either homogeneous Dirichlet or absorbing boundary conditions leads to a linear system (3.8) with coefficient matrix A = L − iB − K (3.13) and real-valued matrices L, B, K.The matrix K is symmetric positive definite, B, L are symmetric positive semidefinite.In case of homogeneous Dirichlet boundary conditions, L is symmetric positive definite and B = 0.

Theorem 4 . 1 .
Let the notation be as in Theorem 3.4, let β > 0, let A be the discretization (3.13) of the stochastic Helmholtz equation (3.1) and M be the discretization (4.2) of the shifted Helmholtz equation (4.1).

Figure 1 :
Figure 1: Images under µ of the real and imaginary axis (solid and dashed circles, respectively) and of the four quadrants; see Remark 4.2 and the proof of Theorem 4.1.

Figure 3
Figure 3 displays the 2-norm condition numbers of A, M , AM −1 , A 0 and AA −1 0

Figure 5 :
Figure 5: Number of GMRES iteration steps (top) and computation time in seconds (bottom) as functions of k for different left and right preconditioners with fixed θ = 0.1 and m = 3.The right panels are zoom-ins.

Figure 7 :
Figure 7: Maximum-norms v i ∞ as a function of i = deg(φ i ).Left: For fixed θ = 0.1 and different values of k.Right: For fixed k = 100 and different values of θ.

Figure 8 :
Figure 8: Norms x m − x m+1 2 as a function of the maximal degree m in the stochastic Galerkin method.Left: For fixed θ = 0.1 and different values of k.Right: For fixed k = 100 and different values of θ.

Figure 9 :
Figure 9: Solving the linear system directly and with GMRES preconditioned by A 0 and M 0 ; see Section 7. Left: Computation time (in seconds) as a function of the polynomial degree r in the stochastic Galerkin method.Right: Relative residual norms in preconditioned GMRES for polynomial degree r = 8.

Figure 11 :
Figure 11: Expected value (top) and variance (bottom) of Re( u m ) and Im( u m ).

Figure 12 :
Figure 12: Left: Euclidean norms x r−1 − x r 2 as a function of the polynomial degree r = 2, . . ., 8, where x r is the solution of Ax = b using total degree r in the stochastic Galerkin method.Right: Magnitudes (7.6) as a function of the degree j.

Table 2 :
. For different total degrees r, number of basis polynomials (n.basis), size of the matrix A, number of non-zero elements in A (nnz), and time to generate A.